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

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

        depth:do iz=nwbott,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)

           vxz:if (minv .eq. maxv) then
                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)
 
           else
                velref: do iref=1,nref                 
                    gs(:,iref)=ws;  gr(:,iref)=wr; 
                    wavenumber2: do ikx = 1,kxm%n
!      phase shift
                        sqa = rec%w%all(iwa)**2 / vref(iref)**2 - kxm%all(ikx)**2;
                        evanescent2: if (sqa > 0) then
                            gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                            gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                        else
                            gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                            gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                        end if evanescent2
                    end do wavenumber2

! Operations in x            

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

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

                do iref=1,nref-1
                interpwvfld: where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
                     ws=gs(:,iref)
                     wr=gr(:,iref)
                     end where interpwvfld
                end do
                where (vel(iz,:) >= vref(nref))
                      ws=gs(:,nref)
                      wr=gr(:,nref)
                end where 

           end if vxz

           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 

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

        depth:do iz=nwbott,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)

           vxz:if (minv .eq. maxv) then
                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)
 
           else
                velref: do iref=1,nref                 
                    gs(:,iref)=ws;  gr(:,iref)=wr; 
                    wavenumber2: do ikx = 1,kxm%n
!      phase shift
                        sqa = rec%w%all(iwa)**2 / vref(iref)**2 - kxm%all(ikx)**2;
                        evanescent2: if (sqa > 0) then
                            gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                            gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                        else
                            gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                            gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                        end if evanescent2
                    end do wavenumber2

! Operations in x            

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

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

                do iref=1,nref-1
                interpwvfld: where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
                     ws=gs(:,iref)
                     wr=gr(:,iref)
                     end where interpwvfld
                end do
                where (vel(iz,:) >= vref(nref))
                      ws=gs(:,nref)
                      wr=gr(:,nref)
                end where 

           end if vxz

           ws=ws*taper
           wr=wr*taper
           GGr(iz,:)=wr
           GGs(iz,:)=ws

        end do depth

end subroutine owwei1
end module
