module OW_diffwei_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(adj,taper,recv,source,vref,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(:)


allocate(wr(kxm%n,nodes));		wr=0
allocate(ws(kxm%n,nodes));		ws=0
allocate(ddwr(kxm%n,nodes));		ddwr=0
allocate(ddws(kxm%n,nodes));		ddws=0
allocate(gs(kxm%n,nref,nodes));	gs=0
allocate(gr(kxm%n,nref,nodes));	gr=0
allocate(gs0(kxm%n,nodes));		gs0=0
allocate(gr0(kxm%n,nodes));		gr0=0
allocate(dwr(kxm%n,nodes));		dwr=0
allocate(dws(kxm%n,nodes));		dws=0
allocate(adwr(kxm%n,nodes));		adwr=0
allocate(adws(kxm%n,nodes));		adws=0
allocate(dwrr(kxm%n,nodes));		dwrr=0
allocate(dwss(kxm%n,nodes));		dwss=0
allocate(dgs(kxm%n,nref,nodes));	dgs=0
allocate(dgr(kxm%n,nref,nodes));	dgr=0
allocate(dgs0(kxm%n,nodes));		dgs0=0
allocate(dgr0(kxm%n,nodes));		dgr0=0
allocate(sqa0(kxm%n,nodes));		sqa0=0
allocate(vref(nref,nodes));		vref=0
allocate(refl(kxm%n,nodes));		refl=0

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

!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

allocate(tfldr(kxm%n,nodes));tfldr=0
call owfftwr(tfldr,nodes)

!$OMP parallel do private (inode,&
                             iw,&
                             iz,&
                            ikx,&
                           iref,&
                         indref,&
                            ixw,&
                            sqa,&
                           minv,&
                           maxv,&
                       rnd_phase) ordered schedule(dynamic,1) num_threads(nodes)


 call owfftw(tfld)

if (adj) then

    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)

else

real,optional					:: rand1(:),rand2(:,:),vref1(:,:,:)
integer						:: iw,indref(1),iz,ikx,iref,ixw,inode,nrefvel
real						:: minv,maxv,vsal,sqa,rnd_phase
real						:: vel(:,:),recset(:,:),souset(:,:),taper(:)
real, allocatable					:: vref(:,:),f(:),g(:)
integer,allocatable					:: indw(:,:)
complex, allocatable				:: ws(:,:),gs0(:,:),gs(:,:,:)
complex, allocatable				:: wr(:,:),gr0(:,:),gr(:,:,:),tfld(:,:)
complex, allocatable				:: tflds(:,:),tfldr(:,:),refl(:,:)
complex						:: GGR(:,:),GGS(:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

allocate(wr(kxm%n,nodes));            wr=0
allocate(ws(kxm%n,nodes));            ws=0
allocate(ddwr(kxm%n,nodes));          ddwr=0
allocate(ddws(kxm%n,nodes));          ddws=0
allocate(gs(kxm%n,nref,nodes));       gs=0
allocate(gr(kxm%n,nref,nodes));       gr=0
allocate(gs0(kxm%n,nodes));           gs0=0
allocate(gr0(kxm%n,nodes));           gr0=0
allocate(dwr(kxm%n,nodes));           dwr=0
allocate(dws(kxm%n,nodes));           dws=0
allocate(adwr(kxm%n,nodes));          adwr=0
allocate(adws(kxm%n,nodes));          adws=0
allocate(dwrr(kxm%n,nodes));          dwrr=0
allocate(dwss(kxm%n,nodes));          dwss=0
allocate(dgs(kxm%n,nref,nodes));      dgs=0
allocate(dgr(kxm%n,nref,nodes));      dgr=0
allocate(dgs0(kxm%n,nodes));          dgs0=0
allocate(dgr0(kxm%n,nodes));          dgr0=0
allocate(sqa0(kxm%n,nodes));          sqa0=0
allocate(vref(nref,nodes));           vref=0
allocate(refl(kxm%n,nodes));          refl=0

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

!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

allocate(tfldr(kxm%n,nodes));tfldr=0
call owfftwr(tfldr,nodes)

!$OMP parallel do private (inode,&
                             iw,&
                             iz,&
                            ikx,&
                           iref,&
                         indref,&
                            ixw,&
                            sqa,&
                           minv,&
                           maxv) ordered schedule(dynamic,1) num_threads(nodes)


do iw = 1,rec%w%n

        inode=omp_get_thread_num()+1; 
        ws(:,inode)=0;wr(:,inode)=0.

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

!!--------------------------------------------------
!! Compute vref
!!


            if (present(vref1)) then
               nrefvel=count(mask = vref1(:,1,iz).ne.-1)
               vref(:,inode)=0.
               vref(1:nrefvel,inode)=1/vref1(:,1,iz)
               gs(:,:,inode)=0; gr(:,:,inode)=0
               if (nrefvel==1) then
                  minv=vref(1,inode);maxv=minv
               else
                  minv=vref(nrefvel,inode);maxv=vref(1,inode)
               end if
            else
               minv=minval(vel(iz,:));maxv=minv;vref(:,inode)=minv
            end if

            tflds(:,inode)=0.
            tflds(:,inode)=cmplx(souset(iz,:),0)
            call owfftexes(tflds(:,inode),inode)
            refl(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
            tflds(:,inode)=0.
            tflds(:,inode)=ws(:,inode)
            call owfftexes(tflds(:,inode),inode)
            ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
            ws(:,inode)=ws(:,inode)+refl(:,inode)

            tfldr(:,inode)=0.
            tfldr(:,inode)=cmplx(recset(iz,:),0)
            call owfftexer(tfldr(:,inode),inode)
            refl(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
            tfldr(:,inode)=0.
            tfldr(:,inode)=wr(:,inode)
            call owfftexer(tfldr(:,inode),inode)
            wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
            wr(:,inode)=wr(:,inode)+refl(:,inode)

            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*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))

            if (minv .eq. maxv .and. maxv.eq.vref(1,inode)) then

               gs0(:,inode)=ws(:,inode)
               gr0(:,inode)=wr(:,inode)

               do ikx = 1,kxm%n
                  sqa = sou%w%all(iw)**2 / minv**2 - kxm%all(ikx)**2

                  if (sqa > 0) then
                     gs0(ikx,inode)=gs0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
                     gr0(ikx,inode)=gr0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                  else
                     gs0(ikx,inode)=gs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                     gr0(ikx,inode)=gr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                  end if
                end do

                tflds(:,inode)=gs0(:,inode)
                call owifftexes(tflds(:,inode),inode)
                ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

                tfldr(:,inode)=gr0(:,inode)
                call owifftexer(tfldr(:,inode),inode)
                wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

            else

                do iref=1,nrefvel                 

                   gs(:,iref,inode)=ws(:,inode)
                   gr(:,iref,inode)=wr(:,inode); 

                   do ikx = 1,kxm%n
                      sqa = sou%w%all(iw)**2 / vref(iref,inode)**2 - kxm%all(ikx)**2;
                      if (sqa > 0) then
                         gs(ikx,iref,inode)=gs(ikx,iref,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                         gr(ikx,iref,inode)=gr(ikx,iref,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                      else
                         gs(ikx,iref,inode)=gs(ikx,iref,inode)*exp(-sqrt(-sqa)*v%z%d)
                         gr(ikx,iref,inode)=gr(ikx,iref,inode)*exp(-sqrt(-sqa)*v%z%d)
                      end if
                   end do

                   tflds(:,inode)=gs(:,iref,inode)
                   call owifftexes(tflds(:,inode),inode)
                   gs(:,iref,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

                   tfldr(:,inode)=gr(:,iref,inode)
                   call owifftexer(tfldr(:,inode),inode)
                   gr(:,iref,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

                   gs(:,iref,inode) = gs(:,iref,inode) * exp(cmplx(0,-v%z%d*sou%w%all(iw)) & 
                                                       * (1./vel(iz,:)-1./vref(iref,inode)));
                   gr(:,iref,inode) = gr(:,iref,inode) * exp(cmplx(0,v%z%d*sou%w%all(iw))  &
                                                       * (1./vel(iz,:)-1./vref(iref,inode)));

                end do

                do iref=2,nrefvel
                   where (vel(iz,:) >= 1/vref1(iref,2,iz) .and. vel(iz,:) < 1/vref1(iref-1,2,iz))
                      ws(:,inode)=gs(:,iref,inode)
                      wr(:,inode)=gr(:,iref,inode)
                   end where
                end do
                where (vel(iz,:) > 1/vref1(1,2,iz))
                   ws(:,inode)=gs(:,1,inode)
                   wr(:,inode)=gr(:,1,inode)
                end where

           end if

           wr(:,inode)=wr(:,inode)*taper
           ws(:,inode)=ws(:,inode)*taper

        end do

        if (mod(iw,20)==1) write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",sou%w%n," / CIG # ",ixs
5       format(A17,i5,2(A8,i5))

        !$OMP ORDERED
        indw(iw,2)=inode
        GGR(iw,:)=wr(:,indw(iw,2))
        GGS(iw,:)=ws(:,indw(iw,2))
        !$OMP END ORDERED

end do

!$OMP END PARALLEL DO

call owfftdestroy()
deallocate(tfldr,tflds)
deallocate(ws,gs0,wr,gr0,refl,f,g,vref,gr,gs)




    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 if
end subroutine
end module
