module OW_datuming

use sep
use OWvelo_areal
use OW_type
use OW_fftw_areal
use OW_parms_areal
use OW_exec_time
use OW_rndseed

implicit none

integer,private					:: iwa, ispa
contains

!---------------------------------------------------------
subroutine owdatum_init(iwa_in,ispa_in)

integer						:: iwa_in,ispa_in

iwa=iwa_in
ispa=ispa_in

end subroutine
!---------------------------------------------------------
subroutine owdatum(taper,vel,recset,souset,GG1,GG2)!,rand1,rand2)

!real,optional					:: rand1,rand2(:)
integer						:: iz,ikx,iref,ixw,nout
real						:: minv,maxv,vsal,sqa,f,g,rnd_phase
real						:: vel(:,:),taper(:)
real, allocatable					:: vref(:),sqa0(:)!,rnd_phase(:)
complex, allocatable				:: ws(:),gs0(:),gs(:,:)
complex, allocatable				:: wr(:),gr0(:),gr(:,:)
complex, allocatable				:: tfld(:),refl(:)
complex						:: GG1(:),GG2(:),recset(:),souset(:)

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

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

ws=souset
wr=recset

call owfftw(tfld)

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

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

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

!      gr0=wr;
!      gs0=ws;

!      Operations in the wavenumber domain
!      phase shift

      do ikx = 1,kxm%n

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

         if (sqa > 0) then
!            gr0(ikx)=gr0(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
!            gs0(ikx)=gs0(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
            wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
            ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
         else
!            gr0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*dwbott)
!            gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*dwbott)
            wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
            ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
         end if
      end do

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

   wr=wr*taper
   ws=ws*taper
!   if (encode) then
!      if(encz) then
!          if (funct.eq.1) then
!             rnd_phase=sin(2*pi*rand2(iz))
!          else if (funct.eq.2) then
!             rnd_phase=cos(2*pi*rand2(iz))
!          else if (funct.eq.3) then
!             rnd_phase=tan(2*pi*rand2(iz))
!          else
!             rnd_phase=rand2(iz)    
!          end if
!          wr=wr*exp(cmplx(0,eps*2*pi*rnd_phase))
!          ws=ws*exp(cmplx(0,eps*2*pi*rnd_phase))
!      end if
!   end if
end do

!if (encode) then
!   if (encw) then
!      wr=wr*exp(cmplx(0,(2*pi*rand1)))
!      ws=ws*exp(cmplx(0,(2*pi*rand1)))
!   else if(enclin) then
!      wr=wr*exp(cmplx(0,(t0*sou%w%all(iwa))))
!      ws=ws*exp(cmplx(0,(t0*sou%w%all(iwa))))
!   end if
!end if

GG1=wr
GG2=ws

deallocate(ws,gs0,gs,wr,gr0,gr,sqa0,vref,tfld)
end subroutine

end module
