module wei_fft
  !! FFTW
  !! used for data(mx,my,hx,hy,w)

  use wei_pro
  use wei_par
  use wei_flags
  use wei_attr
  !---------------
  implicit none

  integer(kind=dbl), private, save :: planF(:) !! forward FFT
  integer(kind=dbl), private, save :: planI(:) !! inverse FFT
  real,              private, save :: scale

  !! temp  array
  complex, dimension(:,:,:,:,:), allocatable, private :: mmhht

contains
  !----------------------------------------------------------------
  subroutine weifft_init()
    integer, dimension(5) :: ni
    logical :: fftmeasure
    integer :: er
    integer :: mpierr
    integer :: iomp
    character(len=128), parameter :: unit='weifft_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    if(wpar%impi==0) call from_param("fftmeasure",fftmeasure,.false.)
    if(wpar%nmpi >0) call MPI_BCAST(fftmeasure, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr)

    allocate(planF(wpar%nomp))
    allocate(planI(wpar%nomp))

    !! FFT
    ni = (/ wpro%amx%n, wpro%amy%n, wpro%ahx%n, wpro%ahy%n, 1 /)
    scale = 1./sqrt(1.*product(ni))

    call weimem_add(8,ni(1)*ni(2)*ni(3)*ni(4)*wpar%nomp,unit)
    allocate(mmhht(   ni(1),ni(2),ni(3),ni(4),wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate mmhht')

    if(wpar%impi==0 .and. wpro%verb) call weimessage("init FFT plan")
    do iomp=1,wpar%nomp
       if(fftmeasure) then
          call fftwnd_f77_create_plan(planF(iomp),5,ni,FFTW_FORWARD  &
          ,    FFTW_MEASURE + FFTW_THREADSAFE + FFTW_IN_PLACE)
          call fftwnd_f77_create_plan(planI(iomp),5,ni,FFTW_BACKWARD &
          ,    FFTW_MEASURE + FFTW_THREADSAFE + FFTW_IN_PLACE) 
       else
          call fftwnd_f77_create_plan(planF(iomp),5,ni,FFTW_FORWARD  &
          ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
          call fftwnd_f77_create_plan(planI(iomp),5,ni,FFTW_BACKWARD &
          ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
       end if
    end do
    if(wpar%impi==0 .and. wpro%verb) call weimessage("OK")

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weifft_init
  !----------------------------------------------------------------
  !! forward FFT
  subroutine weifft(iws,ith,data)
    integer, intent(in) :: iws,ith
    complex, dimension(:,:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weifft'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    mmhht(:,:,:,:,ith) = data(:,:,:,:,iws)
    call fftwnd_f77_one(planF(ith),mmhht(:,:,:,:,ith),0)
    data(:,:,:,:,iws) = mmhht(:,:,:,:,ith) * scale

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weifft
  !----------------------------------------------------------------
  !! inverse FFT
  subroutine weiift(iws,ith,data)
    integer, intent(in) :: iws, ith
    complex, dimension(:,:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weiift'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    mmhht(:,:,:,:,ith)=data(:,:,:,:,iws)
    call fftwnd_f77_one(planI(ith),mmhht(:,:,:,:,ith),0)
    data(:,:,:,:,iws) = mmhht(:,:,:,:,ith) * scale

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weiift
  !----------------------------------------------------------------
  !! forward FFT
  subroutine weifft4d(data)
    complex, dimension(:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weifft4d'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    mmhht(:,:,:,:,1) = data (:,:,:,:)
    call fftwnd_f77_one(planF(1),mmhht(:,:,:,:,1),0)
    data(:,:,:,:) = mmhht(:,:,:,:,1) * scale

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weifft4d
  !----------------------------------------------------------------
  !! inverse FFT
  subroutine weiift4d(data)
    complex, dimension(:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weiift4d'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    mmhht(:,:,:,:,1) = data(:,:,:,:)
    call fftwnd_f77_one(planI(1),mmhht(:,:,:,:,1),0)
    data(:,:,:,:) = mmhht(:,:,:,:,1) * scale

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weiift4d
  !----------------------------------------------------------------
end module wei_fft
