module weiop_box
  !! Born WEMVA
  !! calls the "eXpansion" scattering operator

  use wei_util
  use wei_pro
  use wei_par
  use wei_S
  use wei_kmap
  use wei_boex
  !---------------
  implicit none

  complex, dimension(:,:,:,:,:), pointer, private :: dw_s,dw_r
  complex, dimension(:,:,:,:,:), pointer, private :: ds_s,ds_r
  complex, dimension(:,:,:,:,:), pointer, private :: bs_s,bs_r

  real,    dimension(:,:,:,:),   pointer, private :: kss,krr
contains
  !----------------------------------------------------------------
  integer function weiop_box1_init() result(st)
    integer :: imx,ihx, imy,ihy
    integer :: er
    character(len=128), parameter :: unit='weiop_box1_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

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

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

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

    !! ks(mx,my,hx,hy)
    call weimem_add(4,wpro%amx%n*wpro%amy%n*wpro%ahx%n*wpro%ahy%n,unit)
    allocate(kss(     wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate kss')
    !! kr(mx,my,hx,hy)
    call weimem_add(4,wpro%amx%n*wpro%amy%n*wpro%ahx%n*wpro%ahy%n,unit)
    allocate(krr(     wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate krr')

    do ihy=1,wpro%ahy%n
       do ihx=1,wpro%ahx%n
          do imy=1,wpro%amy%n
             do imx=1,wpro%amx%n
                kss(imx,imy,ihx,ihy) = &
                sqrt(kmap%ksx(imx,ihx)**2   &
                +    kmap%ksy(imy,ihy)**2)

                krr(imx,imy,ihx,ihy) = &
                sqrt(kmap%krx(imx,ihx)**2   &
                +    kmap%kry(imy,ihy)**2)
             end do
          end do
       end do
    end do

    !! init scattering operator
    st=weiboex_init()

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_box1_init
  !----------------------------------------------------------------
  integer function weiop_box1(bwfld,dwfld,dssou,dsrec,iws,izs,ith) result(st)
    complex, dimension(:,:,:,:,:), pointer :: bwfld,dwfld
    complex, dimension(:,:,:,:,:), pointer :: dssou,dsrec
    integer, intent(in)       :: iws,izs,ith
    real                      :: w,kos,kor
    real                      :: sos,sor
    integer :: i1,i2,i3,i4
    integer :: jfk
    character(len=128), parameter :: unit='weiop_box1'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif
    w  = wpro%bw%o + (iws-1)*wpro%bw%d

    sos=   maxval(slo%rs_s(izs,:))
    do jfk=1,size(slo%rs_s(izs,:))
       if(        slo%rs_s(izs,jfk) /=0.) then
          if(     slo%rs_s(izs,jfk)< sos) sos=slo%rs_s(izs,jfk)
       end if
    end do

    sor=   maxval(slo%rs_r(izs,:))
    do jfk=1,size(slo%rs_r(izs,:))
       if(        slo%rs_r(izs,jfk) /=0.) then
          if(     slo%rs_r(izs,jfk)< sor) sor=slo%rs_r(izs,jfk)
       end if
    end do
    
    kos = 2*w * sos
    kor = 2*w * sor

    bs_s(:,:,:,:,ith) = slo%ss_s(:,:,:,:,izs)
    bs_r(:,:,:,:,ith) = slo%ss_r(:,:,:,:,izs)

    if(wpro%adj) then !! velocity analysis
#ifdef DEBUG
       if(wpro%dbg>9) call  in('dS<-dW')
#endif
       dw_s(:,:,:,:,ith) = dwfld(:,:,:,:,iws)
       dw_r(:,:,:,:,ith) = dwfld(:,:,:,:,iws)

       st= weiboex_w2s(bwfld,bs_s,iws,ith,ds_s,dw_s,kss,kos)
       st= weiboex_w2s(bwfld,bs_r,iws,ith,ds_r,dw_r,krr,kor)

       !$OMP ATOMIC
       dssou(:,:,:,:,izs) = &
       dssou(:,:,:,:,izs) + ds_s(:,:,:,:,ith)

       !$OMP ATOMIC
       dsrec(:,:,:,:,izs) = &
       dsrec(:,:,:,:,izs) + ds_r(:,:,:,:,ith)

#ifdef DEBUG
       if(wpro%dbg>9) call out('dS<-dW')
#endif
    else !! scaterring
#ifdef DEBUG
       if(wpro%dbg>9) call  in('dS->dW')
#endif

       ds_s(:,:,:,:,ith) = dssou(:,:,:,:,izs)
       ds_r(:,:,:,:,ith) = dsrec(:,:,:,:,izs)

       st= weiboex_s2w(bwfld,bs_s,iws,ith,ds_s,dw_s,kss,kos)
       st= weiboex_s2w(bwfld,bs_r,iws,ith,ds_r,dw_r,krr,kor)

       !! dW += (dWs + dWr)
       dwfld(:,:,:,:,iws) = &
       dwfld(:,:,:,:,iws) + &
       dw_s (:,:,:,:,ith) + &
       dw_r (:,:,:,:,ith)

#ifdef DEBUG
       if(wpro%dbg>9) call out('dS->dW')
#endif
    end if

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