module OW_wei_areal

use sep
use OWvelo_areal
use OW_type
use OW_image_areal
use OW_fftw_areal
use OW_parms_areal

integer,private					:: ni,iwa,node,inode

contains

!---------------------------------------------------------
subroutine owwei_init(ni_in,iwa_in,node_in,inode_in)

integer						:: ni_in,iwa_in,node_in,inode_in

ni=ni_in
iwa=iwa_in
node=node_in
inode=inode_in

end subroutine
!---------------------------------------------------------
subroutine owwei(taper,recv,source,vsal,vel,GG)

integer						:: iz,ikx,iref,ixw
real						:: sqa,f,g
real						:: minv, maxv
real						:: taper(:),vel(:,:),vsal
real, allocatable					:: vref(:),sqa0(:)
complex, allocatable				:: ws(:),wr(:),gs0(:),gr0(:)
complex, allocatable				:: gs(:,:),gr(:,:),tfld(:)
complex						:: recv(:),source(:),GG(:,:,:)

if (allocated(wr)) deallocate(wr);allocate(wr(img%xm%n));			wr=0
if (allocated(ws)) deallocate(ws);allocate(ws(img%xm%n));			ws=0
if (allocated(gs)) deallocate(gs);allocate(gs(img%xm%n,nref));		gs=0
if (allocated(gr)) deallocate(gr);allocate(gr(img%xm%n,nref));		gr=0
if (allocated(gs0)) deallocate(gs0);allocate(gs0(img%xm%n));		gs0=0
if (allocated(gr0)) deallocate(gr0);allocate(gr0(img%xm%n));		gr0=0
if (allocated(sqa0)) deallocate(sqa0);allocate(sqa0(img%xm%n));		sqa0=0
if (allocated(vref)) deallocate(vref);allocate(vref(nref));			vref=0
if (allocated(tfld)) deallocate(tfld);allocate(tfld(img%xm%n));		tfld=0

!!--------------------------------------------------
!! Setting FFTW plans
!!

call owfftw(tfld)

ws=source 
wr=recv

        depth:do iz=1,img%z%n

!--------------------------------------------------
! Compute vref
!
            call owvref(nref,vsal,vel(iz,:),vref,minv,maxv)

!!!!!!!
            tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
            tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)

            gs0=ws;  gr0=wr;

!      Operations in the wavenumber domain
!      phase shift

                wavenumber1:do ikx = 1,kxm%n

                   sqa = rec%w%all(iwa)**2 / minv**2 - kxm%all(ikx)**2

                   evanescent1:if (sqa > 0) then
                       gs0(ikx)=gs0(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                       gr0(ikx)=gr0(ikx)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
                   else
                       gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
                       gr0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*v%z%d)
                   end if evanescent1
                end do wavenumber1
                tfld=gs0; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
                tfld=gr0; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
 
           ws=ws*taper
           wr=wr*taper

           call owimage(ws,wr,GG(iz,:,:))

        end do depth

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine owwei1(taper,recv,source,vsal,vel,GGr,GGs)

integer						:: iz,ikx,iref,ixw
real						:: sqa,f,g
real						:: minv, maxv
real						:: taper(:),vel(:,:),vsal
real, allocatable					:: vref(:),sqa0(:)
complex, allocatable				:: ws(:),wr(:),gs0(:),gr0(:)
complex, allocatable				:: gs(:,:),gr(:,:),tfld(:)
complex						:: recv(:),source(:),GGr(:,:),GGs(:,:)

if (allocated(wr)) deallocate(wr);allocate(wr(img%xm%n));		wr=0
if (allocated(ws)) deallocate(ws);allocate(ws(img%xm%n));		ws=0
if (allocated(gs)) deallocate(gs);allocate(gs(img%xm%n,nref));		gs=0
if (allocated(gr)) deallocate(gr);allocate(gr(img%xm%n,nref));		gr=0
if (allocated(gs0)) deallocate(gs0);allocate(gs0(img%xm%n));		gs0=0
if (allocated(gr0)) deallocate(gr0);allocate(gr0(img%xm%n));		gr0=0
if (allocated(sqa0)) deallocate(sqa0);allocate(sqa0(img%xm%n));		sqa0=0
if (allocated(vref)) deallocate(vref);allocate(vref(nref));		vref=0
if (allocated(tfld)) deallocate(tfld);allocate(tfld(img%xm%n));		tfld=0

!!--------------------------------------------------
!! Setting FFTW plans
!!

call owfftw(tfld)

ws=source
wr=recv 

        depth:do iz=1,nwbott-1

!--------------------------------------------------
! Compute vref
!
            call owvref(1,minval(vel(iz,:)),vel(iz,:),vref,minv,maxv)

!!!!!!!
            tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
            tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)

            gs0=ws;  gr0=wr;

!      Operations in the wavenumber domain
!      phase shift

           wavenumber1:do ikx = 1,kxm%n

              sqa = rec%w%all(iwa)**2 / minv**2 - kxm%all(ikx)**2

              evanescent1:if (sqa > 0) then
                  gs0(ikx)=gs0(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
                  gr0(ikx)=gr0(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
              else
                  gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*dwbott)
                  gr0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*dwbott)
              end if evanescent1
           end do wavenumber1
           tfld=gs0; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
           tfld=gr0; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
 
           ws=ws*taper
           wr=wr*taper
           GGr(iz,:)=wr
           GGs(iz,:)=ws

        end do depth

end subroutine owwei1
end module
