!!-------------------------------------------------
!!.. FFTW3
!!-------------------------------------------------
!!.. Author: Yaxun Tang, tang@sep.stanford.edu
!!-------------------------------------------------
module base_fftwv3_mod
  use  base_fftwfl_mod
  implicit none

contains

  !!--------------------------------------------------------------------------------
  !!.. 1D FFT along time
  !!.. flag = +1: forward , exp(+iwt)
  !!..        -1: backward, exp(-iwt)
  !!--------------------------------------------------------------------------------
  subroutine sfftw_dft_1dt_init(plan_1dt, vec1d, n1, flag)
    integer*8                :: plan_1dt   
    complex  , dimension(:)  :: vec1d
    integer                  :: n1
    integer                  :: flag

    if (flag == +1) then
       call sfftw_plan_dft_1d(plan_1dt, n1, vec1d, vec1d, FFTW_BACKWARD, FFTW_MEASURE)
    elseif (flag == -1) then
       call sfftw_plan_dft_1d(plan_1dt, n1, vec1d, vec1d, FFTW_FORWARD , FFTW_MEASURE)
    else
       write(0,*)"ERROR: flag"
       stop
    end if
    return
  end subroutine

  !!--------------------------------------------------------------------------------
  !!.. 1D FFT along x
  !!.. flag = +1: forward , exp(-ikx)
  !!..        -1: backward, exp(+ikx)
  !!--------------------------------------------------------------------------------
  subroutine sfftw_dft_1dx_init(plan_1dx, vec1d, n1, flag)
    integer*8                :: plan_1dx        
    complex  , dimension(:)  :: vec1d
    integer                  :: n1
    integer                  :: flag

    if (flag == +1) then
       call sfftw_plan_dft_1d(plan_1dx, n1, vec1d, vec1d, FFTW_FORWARD , FFTW_MEASURE)
    elseif (flag == -1) then
       call sfftw_plan_dft_1d(plan_1dx, n1, vec1d, vec1d, FFTW_BACKWARD, FFTW_MEASURE)
    else
       write(0,*)"ERROR: flag"
       stop
    end if
    return
  end subroutine

  !!---------------------------------------------------------------------------------
  !!.. 2D FFT along x, y
  !!--------------------------------------------------------------------------------- 
  subroutine sfftw_dft_2dx_init(plan_2dx, vec2d, n1, n2, flag)
    integer*8                 :: plan_2dx
    complex  , dimension(:,:) :: vec2d
    integer                   :: n1, n2
    integer                   :: flag

    if (flag == +1) then
       call sfftw_plan_dft_2d(plan_2dx, n1, n2, vec2d, vec2d, FFTW_FORWARD , FFTW_MEASURE)
    elseif (flag == -1) then
       call sfftw_plan_dft_2d(plan_2dx, n1, n2, vec2d, vec2d, FFTW_BACKWARD, FFTW_MEASURE)
    else
       write(0,*)"ERROR: flag"
       stop
    end if
    return
  end subroutine

  !!----------------------------------------------------------------------------------
  !!.. 3D FFT along x, y, z
  !!----------------------------------------------------------------------------------
  subroutine sfftw_dft_3dx_init(plan_3dx, vec3d, n1, n2, n3, flag)
    integer*8                   :: plan_3dx
    complex  , dimension(:,:,:) :: vec3d
    integer                     :: n1, n2, n3
    integer                     :: flag

    if (flag == +1) then
       call sfftw_plan_dft_3d(plan_3dx, n1, n2, n3, vec3d, vec3d, FFTW_FORWARD , FFTW_MEASURE)
    elseif (flag == -1) then
       call sfftw_plan_dft_3d(plan_3dx, n1, n2, n3, vec3d, vec3d, FFTW_BACKWARD, FFTW_MEASURE)
    else
       write(0,*)"ERROR: flag"
    end if
    return
  end subroutine

  !!-----------------------------------------------------------------------------------
  subroutine sfftw_dft_1dt_oper(plan_1dt, vec1d, n1)
    integer*8               :: plan_1dt
    complex  , dimension(:) :: vec1d
    integer                 :: n1
    call sfftw_execute(plan_1dt)
    vec1d = vec1d / sqrt(float(n1)) 
    return
  end subroutine

  !!------------------------------------------------------------------------------------
  subroutine sfftw_dft_1dx_oper(plan_1dx, vec1d, n1)
    integer*8               :: plan_1dx
    complex  , dimension(:) :: vec1d
    integer                 :: n1
    call sfftw_execute(plan_1dx)
    vec1d = vec1d / sqrt(float(n1))
    return
  end subroutine

  !!------------------------------------------------------------------------------------
  subroutine sfftw_dft_2dx_oper(plan_2dx, vec2d, n1, n2)
    integer*8                 :: plan_2dx
    complex  , dimension(:,:) :: vec2d
    integer                   :: n1, n2
    call sfftw_execute(plan_2dx)
    vec2d = vec2d / sqrt(float(n1*n2))
    return
  end subroutine

  !!------------------------------------------------------------------------------------
  subroutine sfftw_dft_3dx_oper(plan_3dx, vec3d, n1, n2, n3)
    integer*8                   :: plan_3dx
    complex  , dimension(:,:,:) :: vec3d
    integer                     :: n1, n2, n3
    call sfftw_execute(plan_3dx)
    vec3d = vec3d / sqrt(float(n1*n2*n3))
    return
  end subroutine



end module
