program FKflt_
  use sep
  use fkvsp_util
  implicit none

  real    :: v1,v2,v3,v4
  integer :: pass

  integer :: nz,nt, pz,pt 
  real    :: oz,ot
  real    :: dz,dt
  integer :: iz,it

  real    :: ow,dw
  real    :: ok,dk

  real,    dimension(:,:), pointer :: rin, rou
  complex, dimension(:,:), pointer :: cin, cou

  real :: w,k,v, t1,t2, coef
  logical :: dofft,doift,fkfilt

  call sep_init()

  call from_history("n1",nz,1 )
  call from_history("o1",oz,0.)
  call from_history("d1",dz,1.)

  call from_history("n2",nt,1 )
  call from_history("o2",ot,0.)
  call from_history("d2",dt,1.)

  call from_param("v1",v1,1.)
  call from_param("v2",v2,1.)
  call from_param("v3",v3,1.)
  call from_param("v4",v4,1.)

  call from_param("pass",pass,0)

  call from_param("fkfilt",fkfilt,.false.)
  call from_param("dofft",dofft,.true.)
  call from_param("doift",doift,.true.)

  call sep_close()

  write(0,*) dofft,fkfilt,doift

  pz=pow2(nz)
  pt=pow2(nt)

  write(0,*) nz,nt
  write(0,*) pz,pt

  allocate(rin(nz,nt),rou(nz,nt))
  allocate(cin(pz,pt),cou(pz,pt))

  ow = -1./dt/2.
  dw =  1./(pt*dt)

  ok = -1./dz/2.
  dk =  1./(pz*dz)

  !----------------------------------------------------------------
  write(0,*) "read"
  call sep_read(rin,'in',nz,nt)
  call pad2(rin,cin)
  !----------------------------------------------------------------
  call ftinit(pz,pt)
  !----------------------------------------------------------------
  if(dofft) then
     write(0,*) "FFT"
     call fft(cin)
  end if

  if(fkfilt) then
     write(0,*) "FLT"
     do iz=1,pz;    k=ok+(iz-1)*dk
        do it=1,pt; w=ow+(it-1)*dw

           v=w/k
!           write(0,*) v

           if(     (v>=v4 .and. v<v3).or.(v<=v4 .and. v>v3)) then
              t2=(v-v4)/(v3-v4)
              coef=(1-t2)*(1-pass)+t2*pass
           else if((v>=v3 .and. v<v2).or.(v<=v3 .and. v>v2)) then
              coef=pass
           else if((v>=v2 .and. v<v1).or.(v<=v2 .and. v>v1)) then
              t1=(v-v1)/(v2-v1)
              coef=(1-t1)*(1-pass)+t1*pass
           else
              coef=(1-pass)
           end if

           cou(iz,it) = coef * cmplx(real(cin(iz,it)), aimag(cin(iz,it)))

        end do
     end do
  else
     cou=cin
  end if

  if(doift) then
     write(0,*) "IFT"
     call ift(cou)
  end if
  !----------------------------------------------------------------
  call unpad2(rou,cou)
  call sep_write(rou)
  !----------------------------------------------------------------
  call exit(0)
end program FKflt_





