module kxkhflt_mod

!  use fft_data_mod
  use OW_fftw_areal
  implicit none

  type axis
     integer 					:: n
     real    					:: o, d, v ! v=(i-1)*d+o
     real, allocatable				:: all(:)
  end type axis

  public :: axis
  public :: kh,kx,kz
  private
  public :: fftvar_init
  public :: kfilt_init, kxkh_init, kzkh_init, kzkx_init, kxkh_fft, kxkh_ifft,kzkh_fft, kzkh_ifft!,kzkx_fft,kzkx_ifft
  public :: kxkh_exe,kzkh_exe,kzkx_exe,merge_init, merge !,kxkh_fft3,kxkh_ifft3, kxkh_flt
  
  real,    parameter 			:: pi=3.1415926535898
  complex, parameter 			:: i = cmplx(0.,1.)
  integer            			:: nx, nh, nz
  integer	     				:: npadx, npadz, npadh
  real					:: kzflt1,kzflt2,kzflt3,kzflt4
  real					:: kzflt5,kzflt6,kzflt7,kzflt8
  real					:: kxflt1,kxflt2,kxflt3,kxflt4
  real					:: kxflt5,kxflt6,kxflt7,kxflt8
  real					:: khflt1,khflt2,khflt3,khflt4
  real					:: khflt5,khflt6,khflt7,khflt8
  real					:: zmerge1,zmerge2,zramp
  real					:: xmerge1,xmerge2,xramp
  real               			:: dx, dh, dz
  real               			:: ox, oh, oz
!  type(fft_type)     			:: fft_1d_idx, fft_1d_idz, fft_1d_idh, fft_2d_idxh, fft_2d_idzh,fft_2d_idzx
  type(axis)         			:: kh, kx, kz

  ! Bug danger! The common variables should be modified ONLY by init & clean. Or at worst, by a single subroutine only.

contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fftvar_init(x, z, h) result(stat)

    type(axis), intent(in) :: x, z, h
    integer                :: ih, ikx, ikz, stat
    type(axis)             :: ks

    nx = x%n;    dx = x%d;    ox = x%o
    nz = z%n;    dz = z%d;    oz = z%o
    nh = h%n;    dh = h%d;    oh = h%o

    stat=1
!!!
    kh%n = nh+npadh
    kh%d = 2*pi/(kh%n*dh)
    kh%o = -kh%n/2*kh%d
    allocate(kh%all(kh%n)); kh%all=0.
    do ih = 1, kh%n
      kh%all(ih) = kh%o+(ih-1)*kh%d
    end do
!!!
    kx%n = nx+npadx
    kx%d = 2*pi/(kx%n*dx)
    kx%o = -kx%n/2*kx%d
    allocate(kx%all(kx%n)); kx%all=0.
    do ikx = 1, kx%n
      kx%all(ikx) = kx%o+(ikx-1)*kx%d
    end do
!!!
    kz%n = nz+npadz
    kz%d = 2*pi/(kz%n*dz)
    kz%o = -kz%n/2*kz%d
    allocate(kz%all(kz%n)); kz%all=0.
    do ikz = 1, kz%n
      kz%all(ikz) = kz%o+(ikz-1)*kz%d
    end do
!!!
!    fft_1d_idx = init_fft_1d(kx%n)
!    fft_1d_idz = init_fft_1d(kz%n)
!    fft_2d_idxh = init_fft_2d(kh%n,kx%n)
!    fft_2d_idzh = init_fft_2d(kz%n,kh%n)
!    fft_2d_idzx = init_fft_2d(kz%n,kx%n)

!    write(0,*) kz%n,kz%o,kz%d
!    write(0,*) kh%n,kh%o,kh%d
!    write(0,*) kx%n,kx%o,kx%d

!    write(0,*) kz%all
!    write(0,*) kh%all
!    write(0,*) kx%all

stat = 0

  end function fftvar_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kfilt_init(npadz_in,npadx_in,npadh_in)
 integer			:: npadz_in,npadx_in,npadh_in
 npadz=npadz_in; npadx=npadx_in; npadh=npadh_in
 end subroutine kfilt_init 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine merge_init(zmerge1_in,zmerge2_in,zramp_in,xmerge1_in,xmerge2_in,xramp_in)

 real				:: zmerge1_in,zmerge2_in,zramp_in
 real				:: xmerge1_in,xmerge2_in,xramp_in

 zmerge1=zmerge1_in;zmerge2=zmerge2_in;zramp=zramp_in
 xmerge1=xmerge1_in;xmerge2=xmerge2_in;xramp=xramp_in

 end subroutine merge_init 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kzkx_init(kzflt1_in,kzflt2_in,kzflt3_in,kzflt4_in, &
		      kzflt5_in,kzflt6_in,kzflt7_in,kzflt8_in, &
		      kxflt1_in,kxflt2_in,kxflt3_in,kxflt4_in, &
   		      kxflt5_in,kxflt6_in,kxflt7_in,kxflt8_in)

  real					:: kzflt1_in,kzflt2_in,kzflt3_in,kzflt4_in
  real					:: kzflt5_in,kzflt6_in,kzflt7_in,kzflt8_in
  real					:: kxflt1_in,kxflt2_in,kxflt3_in,kxflt4_in
  real					:: kxflt5_in,kxflt6_in,kxflt7_in,kxflt8_in

 kzflt1=kzflt1_in; kzflt2=kzflt2_in; kzflt3=kzflt3_in; kzflt4=kzflt4_in
 kzflt5=kzflt5_in; kzflt6=kzflt6_in; kzflt7=kzflt7_in; kzflt8=kzflt8_in
 kxflt1=kxflt1_in; kxflt2=kxflt2_in; kxflt3=kxflt3_in; kxflt4=kxflt4_in
 kxflt5=kxflt5_in; kxflt6=kxflt6_in; kxflt7=kxflt7_in; kxflt8=kxflt8_in

 end subroutine kzkx_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kzkh_init(kzflt1_in,kzflt2_in,kzflt3_in,kzflt4_in, &
		      kzflt5_in,kzflt6_in,kzflt7_in,kzflt8_in, &
		      khflt1_in,khflt2_in,khflt3_in,khflt4_in, &
   		      khflt5_in,khflt6_in,khflt7_in,khflt8_in)

  real					:: kzflt1_in,kzflt2_in,kzflt3_in,kzflt4_in
  real					:: kzflt5_in,kzflt6_in,kzflt7_in,kzflt8_in
  real					:: khflt1_in,khflt2_in,khflt3_in,khflt4_in
  real					:: khflt5_in,khflt6_in,khflt7_in,khflt8_in

 kzflt1=kzflt1_in; kzflt2=kzflt2_in; kzflt3=kzflt3_in; kzflt4=kzflt4_in
 kzflt5=kzflt5_in; kzflt6=kzflt6_in; kzflt7=kzflt7_in; kzflt8=kzflt8_in
 khflt1=khflt1_in; khflt2=khflt2_in; khflt3=khflt3_in; khflt4=khflt4_in
 khflt5=khflt5_in; khflt6=khflt6_in; khflt7=khflt7_in; khflt8=khflt8_in

 end subroutine kzkh_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kxkh_init(khflt1_in,khflt2_in,khflt3_in,khflt4_in, &
		      khflt5_in,khflt6_in,khflt7_in,khflt8_in, &
		      kxflt1_in,kxflt2_in,kxflt3_in,kxflt4_in, &
   		      kxflt5_in,kxflt6_in,kxflt7_in,kxflt8_in)

  real					:: khflt1_in,khflt2_in,khflt3_in,khflt4_in
  real					:: khflt5_in,khflt6_in,khflt7_in,khflt8_in
  real					:: kxflt1_in,kxflt2_in,kxflt3_in,kxflt4_in
  real					:: kxflt5_in,kxflt6_in,kxflt7_in,kxflt8_in

 khflt1=khflt1_in; khflt2=khflt2_in; khflt3=khflt3_in; khflt4=khflt4_in
 khflt5=khflt5_in; khflt6=khflt6_in; khflt7=khflt7_in; khflt8=khflt8_in
 kxflt1=kxflt1_in; kxflt2=kxflt2_in; kxflt3=kxflt3_in; kxflt4=kxflt4_in
 kxflt5=kxflt5_in; kxflt6=kxflt6_in; kxflt7=kxflt7_in; kxflt8=kxflt8_in

 end subroutine kxkh_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kxkh_exe(Ccig)

 complex				:: Ccig(:,:,:)
 integer				:: ix, ih, ikh
 real				:: k1h

   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) Ccig(:,ih,ix)=0.
  	 else if (kx%all(ix).le.kxflt2.and.kh%all(ih).ge.khflt2.and.kh%all(ih).lt.khflt3) then
  	    Ccig(:,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) Ccig(:,ih,ix)=0.
	 else if (kx%all(ix).le.kxflt5.and.kx%all(ix).gt.kxflt6.and.kh%all(ih).le.khflt5.and.kh%all(ih).gt.khflt6) then
	    k1h=(khflt5-khflt6)/(kxflt5-kxflt6)*(kx%all(ix)-kxflt6)+khflt6
	    ikh=int((k1h-kh%all(1))/kh%d)+1
	    if (ih.le.ikh) Ccig(:,ih,ix)=0.
	 else if (kx%all(ix).ge.kxflt6.and.kh%all(ih).le.khflt6.and.kh%all(ih).ge.khflt7) then
	    Ccig(:,ih,ix)=0.
	 else if (kx%all(ix).ge.kxflt7.and.kx%all(ix).le.kxflt8.and.kh%all(ih).lt.khflt7.and.kh%all(ih).ge.khflt8) then
	    k1h=(khflt7-khflt8)/(kxflt7-kxflt8)*(kx%all(ix)-kxflt8)+khflt8
	    ikh=int((k1h-kh%all(1))/kh%d)+1
   	    if (ih.ge.ikh) Ccig(:,ih,ix)=0.
   	 end if
      end do
   end do

 end subroutine kxkh_exe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kzkh_exe(Ccig)

 complex				:: Ccig(:,:,:)
 integer				:: iz, ih, ikz
 real					:: k1z

   do ih=1,kh%n
      do iz=1,kz%n
 	 if (kh%all(ih).ge.khflt1.and.kh%all(ih).lt.khflt2.and.kz%all(iz).le.kzflt1.and.kz%all(iz).gt.kzflt2) then
	    k1z=(kzflt1-kzflt2)/(khflt1-khflt2)*(kh%all(ih)-khflt2)+kzflt2
	    ikz=int((k1z-kz%all(1))/kz%d)+1
  	    if (iz.ge.ikz) Ccig(iz,ih,:)=0.
  	 else if (kh%all(ih).ge.khflt2.and.kh%all(ih).le.khflt3.and.kz%all(iz).ge.kzflt2.and.kz%all(iz).gt.kzflt3) then
  	    Ccig(iz,ih,:)=0.
  	 else if (kh%all(ih).ge.khflt3.and.kh%all(ih).le.khflt4.and.kz%all(iz).ge.kzflt3.and.kz%all(iz).le.kzflt4) then
	    k1z=(kzflt3-kzflt4)/(khflt3-khflt4)*(kh%all(ih)-khflt4)+kzflt4
	    ikz=int((k1z-kz%all(1))/kz%d)+1
	    if (iz.ge.ikz) Ccig(iz,ih,:)=0.
	 else if (kh%all(ih).ge.khflt5.and.kh%all(ih).lt.khflt6.and.kz%all(iz).ge.kzflt5.and.kz%all(iz).lt.kzflt6) then
	    k1z=(kzflt5-kzflt6)/(khflt5-khflt6)*(kh%all(ih)-khflt6)+kzflt6
	    ikz=int((k1z-kz%all(1))/kz%d)+1
	    if (iz.lt.ikz) Ccig(iz,ih,:)=0.
	 else if (kh%all(ih).ge.khflt6.and.kh%all(ih).le.khflt7.and.kz%all(iz).lt.kzflt6.and.kz%all(iz).le.kzflt7) then
	    Ccig(iz,ih,:)=0.
	 else if (kh%all(ih).ge.khflt7.and.kh%all(ih).le.khflt8.and.kz%all(iz).lt.kzflt7.and.kz%all(iz).ge.kzflt8) then
	    k1z=(kzflt7-kzflt8)/(khflt7-khflt8)*(kh%all(ih)-khflt8)+kzflt8
	    ikz=int((k1z-kz%all(1))/kz%d)+1
   	    if (iz.le.ikz) Ccig(iz,ih,:)=0.
   	 end if
      end do
   end do

 end subroutine kzkh_exe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine kzkx_exe(Ccig)

 complex				:: Ccig(:,:,:)
 integer				:: iz, ix, ikz
 real					:: k1z

   do ix=1,kx%n
      do iz=1,kz%n
 	 if (kx%all(ix).ge.kxflt1.and.kx%all(ix).lt.kxflt2.and.kz%all(iz).ge.kzflt1.and.kz%all(iz).lt.kzflt2) then
	    k1z=(kzflt1-kzflt2)/(kxflt1-kxflt2)*(kx%all(ix)-kxflt2)+kzflt2
	    ikz=int((k1z-kz%all(1))/kz%d)+1
  	    if (iz.le.ikz) Ccig(iz,:,ix)=0.
  	 else if (kx%all(ix).ge.kxflt2.and.kx%all(ix).lt.kxflt3.and.kz%all(iz).lt.kzflt2.and.kz%all(iz).lt.kzflt3) then
  	    Ccig(iz,:,ix)=0.
  	 else if (kx%all(ix).ge.kxflt3.and.kx%all(ix).lt.kxflt4.and.kz%all(iz).lt.kzflt3.and.kz%all(iz).ge.kzflt4) then
	    k1z=(kzflt3-kzflt4)/(kxflt3-kxflt4)*(kx%all(ix)-kxflt4)+kzflt4
	    ikz=int((k1z-kz%all(1))/kz%d)+1
	    if (iz.ge.ikz) Ccig(iz,:,ix)=0.
	 else if (kx%all(ix).ge.kxflt5.and.kx%all(ix).lt.kxflt6.and.kz%all(iz).lt.kzflt5.and.kz%all(iz).ge.kzflt6) then
	    k1z=(kzflt5-kzflt6)/(kxflt5-kxflt6)*(kx%all(ix)-kxflt6)+kzflt6
	    ikz=int((k1z-kz%all(1))/kz%d)+1
	    if (iz.ge.ikz) Ccig(iz,:,ix)=0.
	 else if (kx%all(ix).ge.kxflt6.and.kx%all(ix).lt.kxflt7.and.kz%all(iz).ge.kzflt6.and.kz%all(iz).ge.kzflt7) then
	    Ccig(iz,:,ix)=0.
	 else if (kx%all(ix).ge.kxflt7.and.kx%all(ix).le.kxflt8.and.kz%all(iz).ge.kzflt7.and.kz%all(iz).le.kzflt8) then
	    k1z=(kzflt7-kzflt8)/(kxflt7-kxflt8)*(kx%all(ix)-kxflt8)+kzflt8
	    ikz=int((k1z-kz%all(1))/kz%d)+1
   	    if (iz.ge.ikz) Ccig(iz,:,ix)=0.
   	 end if
      end do
   end do

 end subroutine kzkx_exe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function kxkh_fft(cig) result(Ccig)

 complex, allocatable			:: Ccig(:,:,:),tfld(:,:)
 real, allocatable				:: cig1(:,:,:)
 real 					:: cig(:,:,:)
 integer					:: ix, iz, ih

 allocate(Ccig(nz,kh%n,kx%n),cig1(nz,nh,nx))
 allocate(tfld(kh%n,kx%n))
 call owfftw2(tfld,1)
 cig1=cig;
 do ix=1,nx
    do ih=1,nh
          cig1(:,ih,ix)=cig(:,ih,ix)*(-1)**(ix+ih)
    end do
 end do
 Ccig=0.
 Ccig(1:nz,1:nh,1:nx)=cmplx(cig1,0)
write(0,*) "passei kxkh_fft 01"

 do iz=1,nz
    tfld=Ccig(iz,:,:)
    call owfftexe2(tfld)
    Ccig(iz,:,:)=tfld/sqrt(1.*kz%n*kh%n)
 end do
write(0,*) "passei kxkh_fft 02"
 deallocate(tfld)
 call owfftwdestroy2()
write(0,*) "passei kxkh_fft 03"
 end function kxkh_fft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function kzkh_fft(cig) result(Ccig)

 complex, allocatable			:: Ccig(:,:,:),tfld(:,:)
 real, allocatable				:: cig1(:,:,:)
 real 					:: cig(:,:,:)
 integer					:: ix, iz, ih

 allocate(Ccig(kz%n,kh%n,nx),cig1(nz,nh,nx))
 allocate(tfld(kz%n,kh%n))
 call owfftw2(tfld,1)
 cig1=cig;

    do ih=1,nh
       do iz=1,nz
          cig1(iz,ih,:)=cig(iz,ih,:)*(-1)**(iz+ih)
       end do
    end do
 Ccig=0.
 Ccig(1:nz,1:nh,1:nx)=cmplx(cig1,0)

 do ix=1,nx
    tfld=Ccig(:,:,ix)
    call owfftexe2(tfld)
    Ccig(:,:,ix)=tfld/sqrt(1.*kz%n*kh%n)
 end do
 deallocate(tfld)
 call owfftwdestroy2()
 end function kzkh_fft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! function kzkx_fft(cig) result(Ccig)
!
! complex, allocatable			:: Ccig(:,:,:)
! complex, allocatable			:: temp_h(:,:)
! real, allocatable			:: cig1(:,:,:)
! real 					:: cig(:,:,:)
! integer				:: ix, iz, ih
!
! allocate(Ccig(kz%n,nh,kx%n),cig1(nz,nh,nx))
! allocate(temp_h(kz%n,kx%n))
!
! cig1=cig;
!
!    do ix=1,nx
!       do iz=1,nz
!          cig1(iz,:,ix)=cig(iz,:,ix)*(-1)**(iz+ix)
!       end do
!    end do
! Ccig=0.;temp_h=0.
!! Ccig(1:nz,1:nh,1:nx)=cmplx(cig1,0)
!
! do ih=1,nh
!    do ix=1,nx
!       do iz=1,nz
!          temp_h(iz,ix)=cmplx(cig1(iz,ih,ix),0)
!       end do
!    end do
!!    call fft_2d_data(fft_2d_idzx,  1,  kz%n,kx%n,Ccig(:,ih,:))
!    call fft_2d_data(fft_2d_idzx,  1,  kz%n,kx%n,temp_h)
!    do ix=1,kx%n
!       do iz=1,kz%n
!          Ccig(iz,ih,ix)=temp_h(iz,ix)
!       end do
!    end do
! end do
!
! end function kzkx_fft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function kxkh_ifft(Ccig) result(cig) 
 complex				:: Ccig(:,:,:)
 real, allocatable			:: cig(:,:,:)
 complex, allocatable		:: temp(:,:,:),tfld(:,:)
 integer				:: ix, iz, ih

 allocate(cig(nz,nh,nx))
 allocate(temp(nz,kh%n,kx%n))
 allocate(tfld(kh%n,kx%n))
 call owfftw2(tfld,1)

 temp=Ccig
 cig=0.
 do ix=1,kx%n
    do ih=1,kh%n
          if (ih.lt.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
             temp(:,ih+int(kh%n/2),ix+int(kx%n/2))=Ccig(:,ih,ix)
          else if (ih.lt.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
             temp(:,ih+int(kh%n/2),ix-int(kx%n/2))=Ccig(:,ih,ix)
          else if (ih.ge.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
             temp(:,ih-int(kh%n/2),ix-int(kx%n/2))=Ccig(:,ih,ix)
          else if (ih.ge.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
             temp(:,ih-int(kh%n/2),ix+int(kx%n/2))=Ccig(:,ih,ix)
          end if
    end do
 end do
 Ccig=-temp

 do iz=1,nz
    tfld=Ccig(iz,:,:)
    call owfftexe2(tfld)
    Ccig(iz,:,:)=tfld/sqrt(1.*kz%n*kh%n)
 end do
 deallocate(tfld)
 call owfftwdestroy2()

 cig=real(Ccig(:,1:nh,1:nx))

 end function kxkh_ifft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function kzkh_ifft(Ccig) result(cig) 
 complex				:: Ccig(:,:,:)
 real, allocatable			:: cig(:,:,:)
 complex, allocatable		:: temp(:,:,:),tfld(:,:)
 integer				:: ix, iz, ih

 allocate(cig(nz,nh,nx))
 allocate(temp(kz%n,kh%n,nx))
 allocate(tfld(kz%n,kh%n))
 call owfftw2(tfld,1)
 temp=Ccig
 cig=0.
    do ih=1,kh%n
       do iz=1,kz%n
          if (iz.lt.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1) then
             temp(iz+int(kz%n/2),ih+int(kh%n/2),:)=Ccig(iz,ih,:)
          else if (iz.lt.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1) then
             temp(iz+int(kz%n/2),ih-int(kh%n/2),:)=Ccig(iz,ih,:)
          else if (iz.ge.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1) then
             temp(iz-int(kz%n/2),ih-int(kh%n/2),:)=Ccig(iz,ih,:)
          else if (iz.ge.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1) then
             temp(iz-int(kz%n/2),ih+int(kh%n/2),:)=Ccig(iz,ih,:)
          end if
       end do
    end do
 Ccig=-temp

 do ix=1,nx
    tfld=Ccig(:,:,ix)
    call owfftexe2(tfld)
    Ccig(:,:,ix)=tfld/sqrt(1.*kz%n*kh%n)
 end do
 deallocate(tfld)
 call owfftwdestroy2()

 cig=real(Ccig(1:nz,1:nh,:))

 end function kzkh_ifft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! function kzkx_ifft(Ccig) result(cig) 
! complex				:: Ccig(:,:,:)
! real, allocatable			:: cig(:,:,:)
! complex, allocatable			:: temp(:,:,:)
!! complex, allocatable			:: temp_h(:,:)
! integer				:: ix, iz, ih
!
! allocate(cig(nz,nh,nx))
! allocate(temp(kz%n,nh,kx%n))
!! allocate(temp_h(kz%n,kx%n))
!
! temp=Ccig
! cig=0.
!    do ix=1,kx%n
!       do iz=1,kz%n
!          if (iz.lt.int(kz%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),:,ix+int(kx%n/2))=Ccig(iz,:,ix)
!          else if (iz.lt.int(kz%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),:,ix-int(kx%n/2))=Ccig(iz,:,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),:,ix-int(kx%n/2))=Ccig(iz,:,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),:,ix+int(kx%n/2))=Ccig(iz,:,ix)
!          end if
!       end do
!    end do
! Ccig=-temp
!
! do ih=1,nh
!!    do ix=1,kx%n
!!       do iz=1,kz%n
!!          temp_h(iz,ix)=Ccig(iz,ih,ix)
!!       end do
!!    end do
!!    call fft_2d_data(fft_2d_idzx,  -1,  kz%n, kx%n,temp_h)
!    call fft_2d_data(fft_2d_idzx,  -1,  kz%n, kx%n,Ccig(:,ih,:))
!!    do ix=1,kx%n
!!       do iz=1,kz%n
!!          cig(iz,ih,ix)=real(temp_h(iz,ix))
!!       end do
!!    end do
! end do
!
! cig=real(Ccig(1:nz,:,1:nx))
!
! end function kzkx_ifft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! function kxkh_fft3(cig) result(Ccig)
!
! complex, allocatable			:: Ccig(:,:,:)
! real 					:: cig(:,:,:)
! integer				:: ix, iz, ih
!
! allocate(Ccig(kz%n,kh%n,kx%n))
!
! do ix=1,nx
!    do ih=1,nh
!       do iz=1,nz
!          cig(iz,ih,ix)=cig(iz,ih,ix)*(-1)**(ix+ih+iz)
!       end do
!    end do
! end do
! Ccig=0.
! Ccig(1:nz,1:nh,1:nx)=cmplx(cig,0)
!
! do ix=1,kx%n
!    do ih=1,kh%n
!       call fft_1d_data(fft_1d_idz,  1,  kz%n,Ccig(:,ih,ix))
!   end do
! end do
! 
! do iz=1,kz%n
!    call fft_2d_data(fft_2d_idxh,  1,  kh%n,kx%n,Ccig(iz,:,:))
! end do
!
! end function kxkh_fft3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! function kxkh_ifft3(Ccig) result(cig) !result(cig) 
!
! complex				:: Ccig(:,:,:)
! real, allocatable			:: cig(:,:,:)
! complex, allocatable			:: temp(:,:,:)
! integer				:: ix, iz, ih
!
! allocate(cig(nz,nh,nx))
! allocate(temp(kz%n,kh%n,kx%n))
!
! temp=Ccig
! cig=0.
! do ix=1,kx%n
!    do ih=1,kh%n
!       do iz=1,kz%n
!          if (iz.lt.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),ih+int(kh%n/2),ix+int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.lt.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),ih+int(kh%n/2),ix-int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.lt.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),ih-int(kh%n/2),ix-int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.lt.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz+int(kz%n/2),ih-int(kh%n/2),ix+int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),ih-int(kh%n/2),ix-int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ih.ge.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),ih-int(kh%n/2),ix+int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1.and.ix.lt.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),ih+int(kh%n/2),ix+int(kx%n/2))=Ccig(iz,ih,ix)
!          else if (iz.ge.int(kz%n/2)+1.and.ih.lt.int(kh%n/2)+1.and.ix.ge.int(kx%n/2)+1) then
!             temp(iz-int(kz%n/2),ih+int(kh%n/2),ix-int(kx%n/2))=Ccig(iz,ih,ix)
!          end if
!       end do
!    end do
! end do
! Ccig=-temp
!
! deallocate(temp)
!
! do ix=1,kx%n
!    do ih=1,kh%n
!       call fft_1d_data(fft_1d_idz,  -1,  kz%n,Ccig(:,ih,ix))
!   end do
! end do
!
! do iz=1,kz%n
!    call fft_2d_data(fft_2d_idxh,  -1,  kh%n, kx%n,Ccig(iz,:,:))
! end do
!
! cig=real(Ccig(1:nz,1:nh,1:nx))
!
! end function kxkh_ifft3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function merge(cig_flt,cig_or) result(sodmerge)

 real, intent(in)			:: cig_or(:,:,:), cig_flt(:,:,:)
 real, allocatable			:: sodmerge(:,:,:)
 integer				:: ih, ix, iz, izx
 integer				:: izmerge1,izmerge2,izramp
 integer				:: ixmerge1,ixmerge2,ixramp
 real					:: depth,coorx,wx,wz

 allocate(sodmerge(nz,nh,nx))
 
 izmerge1=int((zmerge1-oz)/dz)+1
 izmerge2=int((zmerge2-oz)/dz)+1
 ixmerge1=int((xmerge1-ox)/dx)+1
 ixmerge2=int((xmerge2-ox)/dx)+1
 ixramp=int(xramp/dx)
 izramp=int(zramp/dz)

 do ix=1,nx; do ih=1,nh; do iz=1,nz
    depth=(iz-1)*dz+oz;coorx=(ix-1)*dx+ox
    if (iz.lt.izmerge1.or.ix.lt.ixmerge1.or.iz.gt.izmerge2.or.ix.gt.ixmerge2) then
             sodmerge(iz,ih,ix) = cig_or(iz,ih,ix)
    else if (iz.ge.izmerge1+izramp.and.iz.le.izmerge2-izramp.and.ix.ge.ixmerge1+ixramp.and.ix.le.ixmerge2-ixramp) then
             sodmerge(iz,ih,ix) = cig_flt(iz,ih,ix)
    else if (iz.ge.izmerge1.and.iz.lt.izmerge1+izramp.and.ix.ge.ixmerge1.and.ix.lt.ixmerge1+ixramp) then
             wx=1/xramp*(coorx-xmerge1)
             wz=1/zramp*(depth-zmerge1)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.ge.izmerge1.and.iz.lt.izmerge1+izramp.and.ix.le.ixmerge2.and.ix.gt.ixmerge2-ixramp) then
             wx=1/xramp*(coorx-xmerge2+xramp)
             wz=1/zramp*(depth-zmerge1)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.gt.izmerge2-izramp.and.iz.le.izmerge2.and.ix.ge.ixmerge1.and.ix.lt.ixmerge1+ixramp) then
             wx=1/xramp*(coorx-xmerge1)
             wz=1/zramp*(depth-zmerge2+zramp)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.gt.izmerge2-izramp.and.iz.le.izmerge2.and.ix.le.ixmerge2.and.ix.gt.ixmerge2-ixramp ) then
             wx=1/xramp*(coorx-xmerge2+xramp)
             wz=1/zramp*(depth-zmerge2+zramp)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.ge.izmerge1.and.iz.lt.izmerge1+izramp.and.ix.ge.ixmerge1+ixramp.and.ix.le.ixmerge2-ixramp) then
             wx=1
             wz=1/zramp*(depth-zmerge1)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.gt.izmerge2-izramp.and.iz.le.izmerge2.and.ix.ge.ixmerge1+ixramp.and.ix.le.ixmerge2-ixramp) then
             wx=1
             wz=1/zramp*(depth-zmerge2+zramp)
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.ge.izmerge1+izramp.and.iz.le.izmerge2-izramp.and.ix.ge.ixmerge1.and.ix.lt.ixmerge1+ixramp) then
             wx=1/xramp*(coorx-xmerge1)
             wz=1
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    else if (iz.ge.izmerge1+izramp.and.iz.le.izmerge2-izramp.and.ix.gt.ixmerge2-izramp.and.ix.le.ixmerge2) then
             wx=1/xramp*(coorx-xmerge2+xramp)
             wz=1
             sodmerge(iz,ih,ix) = (wx*wz)*cig_flt(iz,ih,ix) + (1-wx*wz)*cig_or(iz,ih,ix)
    end if
 end do; end do ; end do

 end function merge
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end module kxkhflt_mod
