# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!!------------------------------------------------------------------------------
!!.. Author: Yaxun Tang, tang@sep.stanford.edu
!!------------------------------------------------------------------------------
module wem_wfdexp_mod
  use wem_wnbmap_mod
  use wem_extrpl_mod
  use wem_ssrepi_mod
  implicit none

  real , dimension(:), pointer, private :: kx2, ky2
  integer , private :: nx, ny
  integer , private :: pdip, nvrf
  real , private :: dx, dy, dz
  real , private :: trick
  character(len=5) , private :: oper
contains
  !-------------------------------------------------------------------------------------------------
  !.. oper = 1. "phsft", phase shift
  ! 2. "ssfou", split-step Fourier
  ! 3. "foufd", Fourier finite-difference
  ! 4. "expfd", explicit finite-difference (not ready yet)
  ! 5. "impfd", implicit finite-difference
  ! 6. "phspi", phase-shift plus interpolation
  ! 7. "ssfpi", split-step Fourier plus interpolation
  ! 8. "ffdpi", Fourier finite-difference plus interpolation (not ready yet)
  !-------------------------------------------------------------------------------------------------
  integer function wem_wfdexp_init(kx2_in, ky2_in, nx_in, ny_in, dx_in, dy_in, dz_in, pdip_in, nvrf_in, oper_in, trick_in) result(stat)
    real, dimension(:), target :: kx2_in, ky2_in
    integer :: nx_in, ny_in
    real :: dx_in, dy_in, dz_in
    integer :: nvrf_in, pdip_in
    real :: trick_in
    character(len=5) :: oper_in

    stat = 1

    kx2 => kx2_in
    ky2 => ky2_in
    nx = nx_in
    ny = ny_in
    dx = dx_in
    dy = dy_in
    dz = dz_in
    pdip = pdip_in
    nvrf = nvrf_in
    oper = oper_in
    trick = trick_in

    stat = 0
  end function

  integer function wem_wfdexp_clos() result(stat)
    stat = 1
    nx = 0; ny=0
    dx = 0; dy=0; dz=0
    kx2 => null()
    ky2 => null()
    stat = 0
  end function

  integer function wem_wfdexp_oper(wfld_slice, vmod_slice, vref_slice, ws, opsign, iwrk2d, rwrk1d, cwrk1d, cwrk3d, plan_2dxF, plan_2dxI) result(stat)
    complex, dimension(:,:) :: wfld_slice
    real , dimension(:,:) :: vmod_slice
    real , dimension(:) :: vref_slice
    real :: ws
    integer :: opsign
    integer, dimension(:,:) :: iwrk2d
    real , dimension(:) :: rwrk1d
    complex, dimension(:) :: cwrk1d
    complex, dimension(:,:,:) :: cwrk3d
    integer*8 :: plan_2dxF
    integer*8 :: plan_2dxI

    real :: vmin, vmax, vavr

    stat = 1

    vmin = minval(vmod_slice)
    vmax = maxval(vmod_slice)
    vavr = sum(vmod_slice)/nx/ny

    select case(oper)

      case ("phsft")
           stat = ssr_phs(wfld_slice, vmod_slice, kx2, ky2, vavr, ws, dz, nx, ny, opsign, plan_2dxF, plan_2dxI)

      case ("ssfou")
           stat = ssr_ssf(wfld_slice, vmod_slice, kx2, ky2, vavr, ws, dz, nx, ny, opsign, plan_2dxF, plan_2dxI)

      case ("impfd")
           stat = ssr_ifd(wfld_slice, vmod_slice, kx2, ky2, vmin, ws, dx, dy, dz, trick, nx, ny, pdip, opsign, rwrk1d, cwrk1d, plan_2dxF, plan_2dxI)

      case ("foufd")
           stat = ssr_ffd(wfld_slice, vmod_slice, kx2, ky2, vmin, ws, dx, dy, dz, trick, nx, ny, pdip, opsign, rwrk1d, cwrk1d, plan_2dxF, plan_2dxI)

      case ("phspi")
           stat = ssr_epi(ssr_phs, wfld_slice, vmod_slice, vref_slice, kx2, ky2, ws, dz, nx, ny, nvrf, opsign, iwrk2d, cwrk3d, plan_2dxF, plan_2dxI)

      case ("ssfpi")
           stat = ssr_epi(ssr_ssf, wfld_slice, vmod_slice, vref_slice, kx2, ky2, ws, dz, nx, ny, nvrf, opsign, iwrk2d, cwrk3d, plan_2dxF, plan_2dxI)

    end select

    stat = 0

  end function



end module
