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

    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')

    if(wpar%impi==0) call from_param("nsc",nsc,1)
    if(wpar%nmpi >0) call MPI_BCAST(nsc,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)
    if(nsc>4) nsc=4

    allocate(csc(4),stat=er)
    if(er/=0) call seperr('cant allocate csc')
    csc(1)= 1./  2.
    csc(2)= 3./  8.
    csc(3)= 5./ 16.
    csc(4)=35./128.

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

    real    :: w
    complex :: iwdz
    integer :: isc
#ifdef DEBUG
    character(len=128), parameter :: unit='weiboex_s2w'
    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

    dwfld(:,:,:,:,ith) = iwdz * &
    dslow(:,:,:,:,ith)        * &
    bwfld(:,:,:,:,iws)

    do isc=1,nsc
       dwtmp(:,:,:,:,ith) = iwdz * &
       dslow(:,:,:,:,ith)        * &
       bwfld(:,:,:,:,iws)

       call weifft(ith,ith,dwtmp)
       where(kk/ko <.999) 
          dwtmp(:,:,:,:,ith) = &
          dwtmp(:,:,:,:,ith) * kk**(2*isc)
       elsewhere
          dwtmp(:,:,:,:,ith) = 0.
       end where
       call weiift(ith,ith,dwtmp)

       dwfld(:,:,:,:,ith) = &
       dwfld(:,:,:,:,ith) + &
       dwtmp(:,:,:,:,ith) / (w*bslow(:,:,:,:,ith))**(2*isc) * csc(isc)
    end do

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

    real    :: w
    complex :: iwdz
    integer :: isc
#ifdef DEBUG
    character(len=128), parameter :: unit='weiboex_w2s'
    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

    dslow(:,:,:,:,ith) = iwdz * &
    dwfld(:,:,:,:,ith) * conjg( &
    bwfld(:,:,:,:,iws))

    do isc=1,nsc
       dwtmp(:,:,:,:,ith) = &
       dwfld(:,:,:,:,ith) / (w*bslow(:,:,:,:,ith))**(2*isc) * csc(isc)

       call weiift(ith,ith,dwtmp)
       where(kk/ko <.999) 
          dwtmp(:,:,:,:,ith) = &
          dwtmp(:,:,:,:,ith) * kk**(2*isc)
       elsewhere
          dwtmp(:,:,:,:,ith) = 0.
       end where
       call weifft(ith,ith,dwtmp)

       dslow(:,:,:,:,ith) = &
       dslow(:,:,:,:,ith) + iwdz * &
       dwtmp(:,:,:,:,ith)        * &
       conjg(bwfld(:,:,:,:,iws))
    end do

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
















