 program kxkh_flt

 use sep
 use kxkhflt_mod

 implicit none

 logical					:: apply
 complex, allocatable				:: transf_cig(:,:,:)
 real, allocatable				:: sodcig(:,:,:), Rot_s_cig(:,:,:)
 integer					:: ih,ix,iz
 integer					:: ikh
 real						:: k1h
 integer					:: npadx, npadz, npadh, stat
 real						:: khflt1,kxflt1,khflt2,kxflt2,fatkx
 real						:: khflt3,kxflt3,khflt4,kxflt4
 type(axis)					:: x,z,h


 call sep_init()
 call sep_begin_prog()

 call from_history("n1",z%n); call from_history("d1",z%d); call from_history("o1",z%o)
 call from_history("n2",h%n); call from_history("d2",h%d); call from_history("o2",h%o)
 call from_history("n3",x%n); call from_history("d3",x%d); call from_history("o3",x%o)

 call from_param("npadz",npadz,int(z%n/4))
 call from_param("npadx",npadx,int(x%n/4))
 call from_param("npadh",npadh,int(h%n/4))
 call from_param("apply",apply,.true.)
 call from_param("khflt1",khflt1)
 call from_param("kxflt1",kxflt1)
 call from_param("khflt2",khflt2)
 call from_param("kxflt2",kxflt2)
 call from_param("khflt3",khflt3)
 call from_param("kxflt3",kxflt3)
 call from_param("khflt4",khflt4)
 call from_param("kxflt4",kxflt4)
 call from_param("fatkx",fatkx,1.)

 if (mod(x%n+npadx,2) == 1) then
    npadx=npadx+1
 end if
 if (mod(h%n+npadh,2) == 1) then
    npadh=npadh+1
 end if
 if (mod(z%n+npadz,2) == 1) then
    npadz=npadz+1
 end if

 call kxkh_init(npadz,npadx,npadh)
 stat = fftvar_init(x, z, h)

 allocate(sodcig(z%n,h%n,x%n),transf_cig(z%n+npadz,h%n+npadh,x%n+npadx))
 sodcig=0.;transf_cig=0.

 call sreed("in", sodcig, 4*z%n*h%n*x%n)

 if(apply) then
   allocate(Rot_s_cig(z%n,h%n,x%n))
   call to_history("label1","z(ft)")
   call to_history("label2","h(ft)")
   call to_history("label3","x(ft)")
   call to_history("n1",z%n,"Rot_s_cig"); call to_history("d1",z%d,"Rot_s_cig"); call to_history("o1",z%o,"Rot_s_cig")
   call to_history("n2",h%n,"Rot_s_cig"); call to_history("d2",h%d,"Rot_s_cig"); call to_history("o2",h%o,"Rot_s_cig")
   call to_history("n3",x%n,"Rot_s_cig"); call to_history("d3",x%d,"Rot_s_cig"); call to_history("o3",x%o,"Rot_s_cig")
   transf_cig = kxkh_fft(sodcig)
   do ix=1,kx%n
      do ih=1,kh%n
!         if (kx%all(ix).ge.kxflt1.and.kx%all(ix).lt.kxflt2.and.kh%all(ih).ge.khflt1.and.kh%all(ih).lt.khflt2) then
!	    k1h=(khflt1-khflt2)/(kxflt1-kxflt2)*(kx%all(ix)-kxflt2)+khflt2
!	    ikh=int((k1h-kh%all(1))/kh%d)+1
!  	    if (ih.ge.ikh) transf_cig(:,ih,ix)=0.
!         else if (kx%all(ix).lt.kxflt2.and.kh%all(ih).ge.khflt2.and.kh%all(ih).lt.khflt3) then
!            transf_cig(:,ih,ix)=0.
!         else if (kx%all(ix).lt.kxflt3.and.kx%all(ix).ge.kxflt4.and.kh%all(ih).ge.khflt3.and.kh%all(ih).le.khflt4) then
!	    k1h=(khflt3-khflt4)/(kxflt3-kxflt4)*(kx%all(ix)-kxflt4)+khflt4
!	    ikh=int((k1h-kh%all(1))/kh%d)+1
!   	    if (ih.le.ikh) transf_cig(:,ih,ix)=0.
         if (kx%all(ix).lt.-kxflt1.and.kx%all(ix).ge.-kxflt2.and.kh%all(ih).lt.-khflt1.and.kh%all(ih).ge.-khflt2) then
	    k1h=(khflt1-khflt2)/(kxflt1-kxflt2)*(kx%all(ix)+kxflt2)-khflt2
	    ikh=int((k1h-kh%all(1))/kh%d)+1
   	    if (ih.le.ikh) transf_cig(:,ih,ix)=0.
         else if (kx%all(ix).ge.-kxflt2.and.kh%all(ih).lt.-khflt2.and.kh%all(ih).ge.-khflt3) then
            transf_cig(:,ih,ix)=0.
         else if (kx%all(ix).ge.-kxflt3.and.kx%all(ix).lt.-kxflt4.and.kh%all(ih).lt.-khflt3.and.kh%all(ih).gt.-khflt4) then
	    k1h=(khflt3-khflt4)/(kxflt3-kxflt4)*(kx%all(ix)+kxflt4)-khflt4
	    ikh=int((k1h-kh%all(1))/kh%d)+1
   	    if (ih.ge.ikh) transf_cig(:,ih,ix)=0.
         end if
      end do
   end do
!   do ix=1,kx%n
!      if (kx%all(ix).lt.fatkx*kxflt) then
!	 k1h=khflt/(kx%all(1)-kxflt)*(kx%all(ix)-kxflt);ikh=int((k1h-kh%all(1))/kh%d)+1
!         do ih=1,kh%n
!	    if (ih.ge.ikh.and.ih.le.kh%n-ikh) then
!               transf_cig(:,ih,ix)=0.
!            end if
!	 end do
!      else if (kx%all(ix).gt.-fatkx*kxflt) then
!   	 k1h=khflt/(-kx%all(1)+kxflt)*(kx%all(ix)+kxflt);ikh=int((k1h-kh%all(1))/kh%d)+1
!         do ih=1,kh%n
!   	    if (ih.ge.ikh.and.ih.le.kh%n-ikh) then
!               transf_cig(:,ih,ix)=0.
!            end if
!         end do
!      end if
!   end do
   Rot_s_cig = -kxkh_ifft(transf_cig)
   call srite("Rot_s_cig",Rot_s_cig, 4*z%n*h%n*x%n)
 else   
   allocate(Rot_s_cig(z%n,kh%n,kx%n))
   call to_history("label1","z(ft)")
   call to_history("label2","kh(1/ft)")
   call to_history("label3","kx(1/ft)")
   call to_history("n1",z%n,"Rot_s_cig"); call to_history("d1",z%d,"Rot_s_cig"); call to_history("o1",z%o,"Rot_s_cig")
   call to_history("n2",kh%n,"Rot_s_cig"); call to_history("d2",kh%d,"Rot_s_cig"); call to_history("o2",kh%o,"Rot_s_cig")
   call to_history("n3",kx%n,"Rot_s_cig"); call to_history("d3",kx%d,"Rot_s_cig"); call to_history("o3",kx%o,"Rot_s_cig")
   transf_cig = kxkh_fft(sodcig)
   do ix=1,kx%n
      do ih=1,kh%n
!         if (kx%all(ix).ge.kxflt1.and.kx%all(ix).lt.kxflt2.and.kh%all(ih).ge.khflt1.and.kh%all(ih).lt.khflt2) then
!	    k1h=(khflt1-khflt2)/(kxflt1-kxflt2)*(kx%all(ix)-kxflt2)+khflt2
!	    ikh=int((k1h-kh%all(1))/kh%d)+1
!  	    if (ih.ge.ikh) transf_cig(:,ih,ix)=0.
!         else if (kx%all(ix).lt.kxflt2.and.kh%all(ih).ge.khflt2.and.kh%all(ih).lt.khflt3) then
!            transf_cig(:,ih,ix)=0.
!         else if (kx%all(ix).lt.kxflt3.and.kx%all(ix).ge.kxflt4.and.kh%all(ih).ge.khflt3.and.kh%all(ih).le.khflt4) then
!	    k1h=(khflt3-khflt4)/(kxflt3-kxflt4)*(kx%all(ix)-kxflt4)+khflt4
!	    ikh=int((k1h-kh%all(1))/kh%d)+1
!   	    if (ih.le.ikh) transf_cig(:,ih,ix)=0.
         if (kx%all(ix).lt.-kxflt1.and.kx%all(ix).ge.-kxflt2.and.kh%all(ih).lt.-khflt1.and.kh%all(ih).ge.-khflt2) then
	    k1h=(khflt1-khflt2)/(kxflt1-kxflt2)*(kx%all(ix)+kxflt2)-khflt2
	    ikh=int((k1h-kh%all(1))/kh%d)+1
   	    if (ih.le.ikh) transf_cig(:,ih,ix)=0.
         else if (kx%all(ix).ge.-kxflt2.and.kh%all(ih).lt.-khflt2.and.kh%all(ih).ge.-khflt3) then
            transf_cig(:,ih,ix)=0.
         else if (kx%all(ix).ge.-kxflt3.and.kx%all(ix).lt.-kxflt4.and.kh%all(ih).lt.-khflt3.and.kh%all(ih).gt.-khflt4) then
	    k1h=(khflt3-khflt4)/(kxflt3-kxflt4)*(kx%all(ix)+kxflt4)-khflt4
	    ikh=int((k1h-kh%all(1))/kh%d)+1
   	    if (ih.ge.ikh) transf_cig(:,ih,ix)=0.
         end if
      end do
   end do
   Rot_s_cig = cabs(transf_cig)
   call srite("Rot_s_cig",Rot_s_cig, 4*z%n*kh%n*kx%n)
 end if
5  format(3(A7,I4))
6  format(2(A7,I4),2(A7,F10.6))

 call sep_end_prog()
 end program kxkh_flt
