module mva_phsscr_mod
  implicit none

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

  integer function mva_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 mva_phsscr_clos() result(stat)
    dx = 0.; dy = 0.; dz = 0.
    nx = 0 ; ny = 0
    stat = 0
  end function

  !!.. G = G(P,S)
  integer function mva_phsscr_oper(G, P, S, ws) result(stat)
    complex, dimension(:,:) :: G, P
    real   , dimension(:,:) :: S
    real                    :: ws

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

    iwdz = cmplx(0., -1.)*ws*dz
    dx2  = dx*dx
    w2   = ws*ws
    do ix=1, nx
       S2    = S(ix,1)*S(ix,1)
       if (ix == 1) then
          lapxo = (P(ix+1,1)-P(ix,1))/dx
       else if (ix == nx) then
          lapxo = (P(ix,1)-P(ix-1,1))/dx
       else
          lapxo = (P(ix-1,1)-2.*P(ix,1)+P(ix+1,1))/dx2
       end if
       G(ix,1) = iwdz*(P(ix,1)+0.5/w2/S2*lapxo) 
    end do

    stat = 0
  end function


end module
