module OW_datum

use OW_type
use OW_fftw
use OW_grabvel
use OW_parms

implicit none

contains

!---------------------------------------------------------
subroutine owdatum(isp,wbott,vel,recv,source)

integer						:: iw,iz,ikx,nwbott,i1,isp
real						:: sign, sqa
real						:: wbott,dwbott,vwbott,wbot
real						:: vel(:,:)
real, allocatable				:: sqa0(:),vwaux(:)
complex, allocatable				:: wr(:),gr0(:)
complex, allocatable				:: ws(:),gs0(:)
complex, pointer				:: tfld1(:)
complex	, intent(inout)				:: recv(:,:),source(:,:,:)

sign=-1;

!--------------------------------------------------
! Setting FFTW plans
!
allocate(tfld1(kxm%n))

call owfftw(tfld1)

wbot=wbott
wbott=floor(wbott/img%z%d)*img%z%d
nwbott=floor(wbott/img%z%d)+1
dwbott=wbott!/2
allocate(vwaux(nwbott));vwaux=0.
do i1=1,nwbott
   vwaux(i1)=sum(vel(i1,:))
end do
vwbott=sum(vwaux)/(nwbott*kxm%n)
write(0,*) wbot,wbott,nwbott,dwbott,vwbott

omega:do iw = 1,rec%w%n

        if (allocated(wr)) deallocate(wr);	allocate(wr(rec%h%n));		wr=0
        if (allocated(gr0)) deallocate(gr0);	allocate(gr0(rec%h%n));		gr0=0
        if (allocated(ws)) deallocate(ws);	allocate(ws(rec%h%n));		ws=0
        if (allocated(gs0)) deallocate(gs0);	allocate(gs0(rec%h%n));		gs0=0
        if (allocated(sqa0)) deallocate(sqa0);	allocate(sqa0(rec%h%n));	sqa0=0

        ws(isp)=source(1,iw,1)
        wr=recv(:,iw) 

 !       depth:do iz=1,2
            tfld1=0
            tfld1=wr; call sfftw_execute(planF); wr=tfld1/sqrt(1.*kxm%n)
            tfld1=ws; call sfftw_execute(planF); ws=tfld1/sqrt(1.*kxm%n)
            gr0=wr; gs0=ws;

!      Operations in the wavenumber domain
!      phase shift
                wavenumber1:do ikx = 1,kxm%n
                   sqa = rec%w%all(iw)**2 / vwbott**2 - kxm%all(ikx)**2

                   evanescent1:if (sqa > 0) then
                       gs0(ikx)=gr0(ikx)*exp(cmplx(0,-sign*dwbott*sqrt(sqa)))
                       gr0(ikx)=gs0(ikx)*exp(cmplx(0,sign*dwbott*sqrt(sqa)))
                   else
                       gs0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*dwbott)
                       gr0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*dwbott)
                   end if evanescent1
                end do wavenumber1

                   tfld1=gr0; call sfftw_execute(planI); wr=tfld1/sqrt(1.*kxm%n)
                   tfld1=gs0; call sfftw_execute(planI); ws=tfld1/sqrt(1.*kxm%n)

 !       end do depth
        
        recv(:,iw)=wr
        source(:,iw,1)=ws

end do omega

write(0,*)  "DONE DATUMING"
write(0,*)  "============="

end subroutine owdatum

end module