module wemva_phsscr_mod
  implicit none

  real                                , private :: dx, dy, dz
  integer                             , private :: nx, ny
contains

  integer function wemva_phsscr_init(dx_in, dy_in, dz_in, nx_in, ny_in) result(stat)
    real    :: dx_in, dy_in, dz_in
    integer :: nx_in, ny_in
 
    dx = dx_in
    dy = dy_in
    dz = dz_in
    nx = nx_in
    ny = ny_in

    stat = 0
  end function 


  integer function wemva_phsscr_clos() result(stat)
    dx = 0.; dy = 0.; dz = 0.
    nx = 0 ; ny = 0
    stat = 0
  end function

  !!.. G = G(P,S)
  integer function wemva_phsscr_oper(dwfd, bwfd, bvel, ws, opsign) result(stat)
    complex, dimension(:,:) :: dwfd
    complex, dimension(:,:) :: bwfd
    real   , dimension(:,:) :: bvel
    real                    :: ws
    integer                 :: opsign 


    real                    :: dx2, w2, s2
    complex                 :: iwdz, lapxo
    integer                 :: ix

    iwdz = cmplx(0., 1.)*ws*dz*opsign 
    dx2  = dx*dx
    w2   = ws*ws
    do ix=1, nx
       s2    = 1./bvel(ix,1)/bvel(ix,1)
       if (ix == 1) then
          lapxo = (bwfd(ix+1,1)-bwfd(ix,1))/dx
       else if (ix == nx) then
          lapxo = (bwfd(ix,1)-bwfd(ix-1,1))/dx
       else
          lapxo = (bwfd(ix-1,1)-2.*bwfd(ix,1)+bwfd(ix+1,1))/dx2
       end if
       dwfd(ix,1) = iwdz*(bwfd(ix,1)+0.5/w2/s2*lapxo) 
    end do

    stat = 0
  end function


end module
