module wei_rysc
  use wei_util
  use wei_pro
  use wei_fft
  use wei_cps
  !---------------
  implicit none
  complex, dimension(:,:,:,:,:), pointer,     private :: dwtmp
  real, private :: ryeps
contains
  !----------------------------------------------------------------
  integer function weirysc_init() result(st)
    integer :: er
    character(len=128), parameter :: unit='weirysc_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    call weicps_init()

    call weimem_add(8,wpro%amx%n*wpro%amy%n*wpro%ahx%n*wpro%ahy%n*wpar%nomp,unit)
    allocate(dwtmp(   wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n,wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate dwtmp')

    ryeps = 0.1

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weirysc_init
  !----------------------------------------------------------------
  !! ===============
  !! forward: dS->dF
  !! ===============
  integer function weirysc_s2w(bwfld,iws,ith,dslow,dphas,kk,ko) result(st)
    integer, intent(in)                    :: iws,ith
    complex, dimension(:,:,:,:,:), pointer :: dslow,bwfld,dphas
    real,    dimension(:,:,:,:),   pointer :: kk
    real,    intent(in)                    :: ko

    real    :: w, eps
    complex :: iwdz
    character(len=128), parameter :: unit='weirysc_s2w'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif
    w    = wpro%bw%o + (iws-1)*wpro%bw%d
    iwdz = (0,1)*wpro%opsign*w*wpro%az%d
    eps = ryeps*maxval(cabs(bwfld(:,:,:,:,iws)))

    dwtmp(:,:,:,:,ith) = iwdz * &
    dslow(:,:,:,:,ith) *        &
    bwfld(:,:,:,:,iws)              !!                 W dS i w dz

    call weifft(ith,ith,dwtmp)      !!             FFT[W dS i w dz]
    call weicps(dwtmp,ith,kk,ko)    !!      dkz/ds FFT[W dS i w dz]
    call weiift(ith,ith,dwtmp)      !! IFT[ dkz/ds FFT[W dS i w dz]]

    dphas(:,:,:,:,ith) = 0.
    where(cabs(bwfld(:,:,:,:,iws))>eps)
       dphas(:,:,:,:,ith) =        &
       dwtmp(:,:,:,:,ith) /        &
       bwfld(:,:,:,:,iws)           !! IFT[ dkz/ds FFT[W dS i w dz]]/W
    end where 

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
    st=weiOK
  end function weirysc_s2w
  !----------------------------------------------------------------
  !! ================
  !! adjoint:  dS<-dF
  !! ================
  integer function weirysc_w2s(bwfld,iws,ith,dslow,dphas,kk,ko) result(st)
    integer, intent(in)                    :: iws,ith
    complex, dimension(:,:,:,:,:), pointer :: dslow,bwfld,dphas
    real,    dimension(:,:,:,:),   pointer :: kk
    real,    intent(in)                    :: ko

    real    :: w, eps
    complex :: iwdz
    character(len=128), parameter :: unit='weirysc_w2s'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif
    w    = wpro%bw%o + (iws-1)*wpro%bw%d
    iwdz = (0,1)*wpro%opsign*w*wpro%az%d
    eps = ryeps*maxval(cabs(bwfld(:,:,:,:,iws)))

    dwtmp(:,:,:,:,ith) = 0.
    where(cabs(bwfld(:,:,:,:,iws))>eps)
       dwtmp(:,:,:,:,ith) =        &
       dphas(:,:,:,:,ith) / conjg( &
       bwfld(:,:,:,:,iws))         !!                         dF / W
    end where

    call weiift(ith,ith,dwtmp)     !!                     IFT[dF / W] 
    call weicps(dwtmp,ith,kk,ko)   !!            (dkz/ds) IFT[dF / W]
    call weifft(ith,ith,dwtmp)     !!       FFT[ (dkz/ds) IFT[dF / W]]

    dslow(:,:,:,:,ith) = iwdz * &
    dwtmp(:,:,:,:,ith) * conjg( &
    bwfld(:,:,:,:,iws))            !!i w dz FFT[ (dkz/ds) IFT[dF / W]] W* 

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
    st=weiOK
  end function weirysc_w2s
  !----------------------------------------------------------------
end module wei_rysc








