module OW_wei_areal

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

implicit none

integer,private					:: ni,iwa

contains

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

integer						:: ni_in,iwa_in

ni=ni_in
iwa=iwa_in

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

complex,optional					:: GGr(:,:),GGs(:,:)
complex						:: recv(:),source(:)
real						:: GG(:,:,:)
integer						:: iz,ikx,iref,ixw,j
real						:: sqa,f,g,minv, maxv,vsal
real						:: taper(:),vel(:,:),c(5),deltaS(:,:)
real, allocatable					:: vref(:),sqa0(:)
complex, allocatable				:: ws(:),wr(:),gs0(:),gr0(:)
complex, allocatable				:: ddws(:),ddwr(:),gs(:,:),gr(:,:),tfld(:)
complex, allocatable				:: dws(:),dwr(:),dgs0(:),dgr0(:)
complex, allocatable				:: dwss(:),dwrr(:),dgs(:,:),dgr(:,:)
complex,allocatable					:: awr(:),aws(:),adwr(:),adws(:)


 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(ddwr)) deallocate(ddwr);allocate(ddwr(img%xm%n));	ddwr=0
 if (allocated(ddws)) deallocate(ddws);allocate(ddws(img%xm%n));	ddws=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(dwr)) deallocate(dwr);allocate(dwr(img%xm%n));	dwr=0
 if (allocated(dws)) deallocate(dws);allocate(dws(img%xm%n));	dws=0
 if (allocated(adwr)) deallocate(adwr);allocate(adwr(img%xm%n));	adwr=0
 if (allocated(adws)) deallocate(adws);allocate(adws(img%xm%n));	adws=0
 if (allocated(dwrr)) deallocate(dwrr);allocate(dwrr(img%xm%n));	dwrr=0
 if (allocated(dwss)) deallocate(dwss);allocate(dwss(img%xm%n));	dwss=0
 if (allocated(dgs)) deallocate(dgs);allocate(dgs(img%xm%n,nref));	dgs=0
 if (allocated(dgr)) deallocate(dgr);allocate(dgr(img%xm%n,nref));	dgr=0
 if (allocated(dgs0)) deallocate(dgs0);allocate(dgs0(img%xm%n));	dgs0=0
 if (allocated(dgr0)) deallocate(dgr0);allocate(dgr0(img%xm%n));	dgr0=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

 c=(/1.,0.5,0.375,0.3125,0.2734375/)

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

 call owfftw(tfld)

    ws=source 
    wr=recv
    dws=0
    dwr=0

    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)

       ddwr=0;ddws=0

       do j=0,4
          dwss=kxm%all**(2*j)*ws
          dwrr=-kxm%all**(2*j)*wr
          tfld=dwss; call sfftw_execute(planI);
          dwss=tfld/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1)
          tfld=dwrr; call sfftw_execute(planI)
          dwrr=tfld/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1)
          ddws=ddws+dwss
          ddwr=ddwr+dwrr
       end do
    
!    dws=dws+ddws*((iz-1)*v%z%d+v%z%o)*rec%w%all(iwa)*exp(cmplx(0,pi/2))
!    dwr=dwr+ddwr*((iz-1)*v%z%d+v%z%o)*rec%w%all(iwa)*exp(cmplx(0,pi/2))
       dws=dws+ddws*img%z%all(iz)*rec%w%all(iwa)*exp(cmplx(0,pi/2))
       dwr=dwr+ddwr*img%z%all(iz)*rec%w%all(iwa)*exp(cmplx(0,pi/2))

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

       if (minv .eq. maxv) then
          gs0=ws;      gr0=wr;
          dgs0=dws;  dgr0=dwr;

!      Operations in the wavenumber domain
!      phase shift

          do ikx = 1,kxm%n

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

             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)))
                dgs0(ikx)=dgs0(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                dgr0(ikx)=dgr0(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)
                dgs0(ikx)=dgs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
                dgr0(ikx)=dgr0(ikx)*exp(-sqrt(-sqa)*v%z%d)
             end if
          end do
          tfld=gs0; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
          tfld=gr0; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
          tfld=dgs0; call sfftw_execute(planI); dws=tfld/sqrt(1.*kxm%n)
          tfld=dgr0; call sfftw_execute(planI); dwr=tfld/sqrt(1.*kxm%n)
 
       else
          do iref=1,nref                 
             gs(:,iref)=ws;  gr(:,iref)=wr; 
             dgs(:,iref)=dws;  dgr(:,iref)=dwr; 
             do ikx = 1,kxm%n
!      phase shift
                sqa = rec%w%all(iwa)**2 / vref(iref)**2 - kxm%all(ikx)**2;
                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)));
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                   dgr(ikx,iref)=dgr(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)
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                   dgr(ikx,iref)=dgr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                end if
             end do

! 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)
             tfld=dgs(:,iref); call sfftw_execute(planI); dgs(:,iref)=tfld/sqrt(1.*kxm%n)
             tfld=dgr(:,iref); call sfftw_execute(planI); dgr(:,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)));
             dgs(:,iref) = dgs(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-1./vref(iref)));
             dgr(:,iref) = dgr(:,iref) * exp(cmplx(0,-v%z%d*rec%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))
                ws=gs(:,iref)
                wr=gr(:,iref)
                dws=dgs(:,iref)
                dwr=dgr(:,iref)
             end where
          end do
          where (vel(iz,:) >= vref(nref))
             ws=gs(:,nref)
             wr=gr(:,nref)
             dws=dgs(:,nref)
             dwr=dgr(:,nref)
          end where 
 
       end if

       ws=ws*taper
       wr=wr*taper
       dws=dws*taper
       dwr=dwr*taper
 
       if (exist_file("ws")) then
          GGr(iz,:)=wr
          GGs(iz,:)=ws
       end if

       adwr=dwr*deltaS(iz,:)
       adws=dws*deltaS(iz,:)

       call owimage(iwa,ws,wr,adws,adwr,GG(iz,:,:))
!    call owimage0(iwa,ws,wr,GG(iz,:,:))
	
    end do

! if (nwbott.ne.1) then
!    call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),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 = rec%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

    deallocate(ws,wr,gs0,gr0,gs,gr,tfld,vref,sqa0)
    deallocate(dws,dwr,ddws,ddwr,dgs0,dgr0,dgs,dgr)
end subroutine
end module
