module OW_ctrillwei

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 owctrlwei_init(iwa_in,ispa_in)

integer						:: iwa_in,ispa_in

iwa=iwa_in
ispa=ispa_in

end subroutine
!---------------------------------------------------------
!subroutine owctrlwei(taper,vsal,vel,recset,souset,GG1,GG2,rand1,rand2,window,wwr,wws)
subroutine owctrlwei(taper,vsal,vel,recset,souset,GG1,GG2,rand1,rand2,wwr,wws)

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

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
!
call owfftw(tfld)
do iz=sog%z%n,nwbott,-1
!--------------------------------------------------
! Compute vref
!
   call owvref(nref,vsal,vel(iz,:),vref,minv,maxv)

   tfld=0.
   tfld=cmplx(recset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=0.
   tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
   wr=wr+refl

   tfld=0.
   tfld=cmplx(souset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=0.
   tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
   ws=ws+refl

   if (minv .eq. maxv) then
      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,v%z%d*sqrt(sqa)))
            gs0(ikx)=gs0(ikx)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
         else
            gr0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*v%z%d)
            gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
         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)

   else
      do iref=1,nref                 
         gr(:,iref)=wr
         gs(:,iref)=ws
         do ikx = 1,kxm%n
!      phase shift
            sqa = sou%w%all(iwa)**2 / (vref(iref))**2 - kxm%all(ikx)**2;
            if (sqa > 0) then
               gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
               gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
            else
               gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
               gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
            end if   
         end do

! Operations in x            

         tfld=gr(:,iref); call sfftw_execute(planI); gr(:,iref)=tfld/sqrt(1.*kxm%n)
         tfld=gs(:,iref); call sfftw_execute(planI); gs(:,iref)=tfld/sqrt(1.*kxm%n)

 !      split step
         gr(:,iref) = gr(:,iref) * exp(cmplx(0,v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
         gs(:,iref) = gs(:,iref) * exp(cmplx(0,-v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
      end do

      do iref=1,nref-1
         where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
            wr=gr(:,iref)
            ws=gs(:,iref)
         end where
      end do
      where (vel(iz,:) >= vref(nref))
         wr=gr(:,nref)
         ws=gs(:,nref)
      end where 
   end if
   wr=wr*taper
   ws=ws*taper
   if (encode.and.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
!      if (nhorz.gt.0) then
!         where(iz.ge.window(:,1).and.iz.le.window(:,2).or.iz.ge.window(:,3).and.iz.le.window(:,4)) 
!            wr=wr*exp(cmplx(0,eps*2*pi*rnd_phase))
!            ws=ws*exp(cmplx(0,eps*2*pi*rnd_phase))
!         end where
!      else
            wr=wr*exp(cmplx(0,eps*2*pi*rnd_phase))
            ws=ws*exp(cmplx(0,eps*2*pi*rnd_phase))
!      end if
   end if
   if (exist_file("wwr")) then
      wwr(iz,:)=wr
      wws(iz,:)=ws
   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

if (nwbott.ne.1) then
call owvref(1,minval(vel(nwbott-1,:)),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)
!      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
         wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
         ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
      else
         wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
         ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
      end if
   end do
   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
end if
GG1=wr
GG2=ws

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

!---------------------------------------------------------
subroutine owctrlweiHD(taper,vsal,vel,recset,souset,GG1,GG2,HD1,HD2)

real						:: HD1(:),HD2(:)
!complex						:: HD1(:),HD2(:)
integer						:: iz,ikx,iref,ixw,nout,seed_1,seed_2
!integer						:: window(:,:)
real						:: minv,maxv,vsal,sqa,f,g,rnd_phase
real						:: vel(:,:),recset(:,:),souset(:,:),taper(:)
real, allocatable					:: vref(:),sqa0(:)!,rnd_phase(:)
complex, allocatable				:: ws1(:),gs01(:),gs1(:,:)
complex, allocatable				:: wr1(:),gr01(:),gr1(:,:)
complex, allocatable				:: tfld(:),refl(:)
complex, allocatable				:: ws2(:),gs02(:),gs2(:,:)
complex, allocatable				:: wr2(:),gr02(:),gr2(:,:)
complex						:: GG1(:,:),GG2(:,:)


allocate(ws1(kxm%n));	ws1=0
allocate(gs1(kxm%n,nref));	gs1=0
allocate(gs01(kxm%n));	gs01=0
allocate(wr1(kxm%n));	wr1=0
allocate(gr1(kxm%n,nref));	gr1=0
allocate(gr01(kxm%n));	gr01=0
allocate(ws2(kxm%n));	ws2=0
allocate(gs2(kxm%n,nref));	gs2=0
allocate(gs02(kxm%n));	gs02=0
allocate(wr2(kxm%n));	wr2=0
allocate(gr2(kxm%n,nref));	gr2=0
allocate(gr02(kxm%n));	gr02=0

allocate(sqa0(kxm%n));sqa0=0
allocate(vref(nref));vref=0
allocate(tfld(kxm%n));tfld=0
allocate(refl(kxm%n));refl=0

!--------------------------------------------------
! Setting FFTW plans
!
call owfftw(tfld)


do iz=sog%z%n,nwbott,-1

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

   tfld=cmplx(recset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=wr1; call sfftw_execute(planF); wr1=tfld/sqrt(1.*kxm%n)
   wr1=wr1+refl
   tfld=wr2; call sfftw_execute(planF); wr2=tfld/sqrt(1.*kxm%n)
   wr2=wr2+refl

   tfld=cmplx(souset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=ws1; call sfftw_execute(planF); ws1=tfld/sqrt(1.*kxm%n)
   ws1=ws1+refl
   tfld=ws2; call sfftw_execute(planF); ws2=tfld/sqrt(1.*kxm%n)
   ws2=ws2+refl

   if (minv .eq. maxv) then
      gr01=wr1;
      gs01=ws1;
      gr02=wr2;
      gs02=ws2;

!      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
            gr01(ikx)=gr01(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
            gs01(ikx)=gs01(ikx)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
            gr02(ikx)=gr02(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
            gs02(ikx)=gs02(ikx)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
         else
            gr01(ikx)=gr01(ikx)*exp(-sqrt(-sqa)*v%z%d)
            gs01(ikx)=gs01(ikx)*exp(-sqrt(-sqa)*v%z%d)
            gr02(ikx)=gr02(ikx)*exp(-sqrt(-sqa)*v%z%d)
            gs02(ikx)=gs02(ikx)*exp(-sqrt(-sqa)*v%z%d)
         end if
      end do

      tfld=gr01; call sfftw_execute(planI); wr1=tfld/sqrt(1.*kxm%n)
      tfld=gs01; call sfftw_execute(planI); ws1=tfld/sqrt(1.*kxm%n)
      tfld=gr02; call sfftw_execute(planI); wr2=tfld/sqrt(1.*kxm%n)
      tfld=gs02; call sfftw_execute(planI); ws2=tfld/sqrt(1.*kxm%n)

   else

      do iref=1,nref                 
         gr1(:,iref)=wr1
         gs1(:,iref)=ws1
         gr2(:,iref)=wr2
         gs2(:,iref)=ws2

         do ikx = 1,kxm%n
!      phase shift
            sqa = sou%w%all(iwa)**2 / (vref(iref))**2 - kxm%all(ikx)**2

            if (sqa > 0) then
               gr1(ikx,iref)=gr1(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
               gs1(ikx,iref)=gs1(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
               gr2(ikx,iref)=gr2(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
               gs2(ikx,iref)=gs2(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
            else
               gr1(ikx,iref)=gr1(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
               gs1(ikx,iref)=gs1(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
               gr2(ikx,iref)=gr2(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
               gs2(ikx,iref)=gs2(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
            end if   
         end do

! Operations in x            

         tfld=gr1(:,iref); call sfftw_execute(planI); gr1(:,iref)=tfld/sqrt(1.*kxm%n)
         tfld=gs1(:,iref); call sfftw_execute(planI); gs1(:,iref)=tfld/sqrt(1.*kxm%n)
         tfld=gr2(:,iref); call sfftw_execute(planI); gr2(:,iref)=tfld/sqrt(1.*kxm%n)
         tfld=gs2(:,iref); call sfftw_execute(planI); gs2(:,iref)=tfld/sqrt(1.*kxm%n)

 !      split step
         gr1(:,iref) = gr1(:,iref) * exp(cmplx(0,v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
         gs1(:,iref) = gs1(:,iref) * exp(cmplx(0,-v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
         gr2(:,iref) = gr2(:,iref) * exp(cmplx(0,v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
         gs2(:,iref) = gs2(:,iref) * exp(cmplx(0,-v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
      end do

      do iref=1,nref-1
         where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
            wr1=gr1(:,iref)
            ws1=gs1(:,iref)
            wr2=gr2(:,iref)
            ws2=gs2(:,iref)
         end where
      end do
      where (vel(iz,:) >= vref(nref))
         wr1=gr1(:,nref)
         ws1=gs1(:,nref)
         wr2=gr2(:,nref)
         ws2=gs2(:,nref)
      end where 
   end if

   wr1=wr1*exp(cmplx(0,eps*pi*HD1(iz)))*taper
   ws1=ws1*exp(cmplx(0,eps*pi*HD1(iz)))*taper
   wr2=wr2*exp(cmplx(0,eps*pi*HD2(iz)))*taper
   ws2=ws2*exp(cmplx(0,eps*pi*HD2(iz)))*taper

end do

if (nwbott.ne.1) then
   call owvref(1,minval(vel(nwbott-1,:)),vel(iz,:),vref,minv,maxv)
   tfld=wr1; call sfftw_execute(planF); wr1=tfld/sqrt(1.*kxm%n)
   tfld=ws1; call sfftw_execute(planF); ws1=tfld/sqrt(1.*kxm%n)
   tfld=wr2; call sfftw_execute(planF); wr2=tfld/sqrt(1.*kxm%n)
   tfld=ws2; call sfftw_execute(planF); ws2=tfld/sqrt(1.*kxm%n)
!      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
         wr1(ikx)=wr1(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
         ws1(ikx)=ws1(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
         wr2(ikx)=wr2(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
         ws2(ikx)=ws2(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
      else
         wr1(ikx)=wr1(ikx)*exp(-sqrt(-sqa)*dwbott)
         ws1(ikx)=ws1(ikx)*exp(-sqrt(-sqa)*dwbott)
         wr2(ikx)=wr2(ikx)*exp(-sqrt(-sqa)*dwbott)
         ws2(ikx)=ws2(ikx)*exp(-sqrt(-sqa)*dwbott)
      end if
   end do
   tfld=wr1; call sfftw_execute(planI); wr1=tfld/sqrt(1.*kxm%n)
   tfld=ws1; call sfftw_execute(planI); ws1=tfld/sqrt(1.*kxm%n)
   wr1=wr1*taper
   ws1=ws1*taper
   tfld=wr2; call sfftw_execute(planI); wr2=tfld/sqrt(1.*kxm%n)
   tfld=ws2; call sfftw_execute(planI); ws2=tfld/sqrt(1.*kxm%n)
   wr2=wr2*taper
   ws2=ws2*taper
end if

GG1(:,1)=wr1
GG2(:,1)=ws1
GG1(:,2)=wr2
GG2(:,2)=ws2

deallocate(ws1,gs01,gs1,wr1,gr01,gr1,sqa0,vref,tfld)
deallocate(ws2,gs02,gs2,wr2,gr02,gr2)

end subroutine

end module
