module OW_pspi_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal
use OWvelo_areal

implicit none

contains

!---------------------------------------------------------
subroutine owpspi(iw,inode,sign,ws,vel,tfld,gat,vref1) 

real,optional					:: vref1(:,:)
character*3					:: gat
integer						:: iw,iz,ikx,iref,iref1,ixw,i,j
integer						:: nwb,inode,nrefvel,sign
integer						:: indref(1)
real						:: sqa,a,dwb,f,g
real						:: minv, maxv
real, allocatable					:: vref(:)
complex						:: ws(:),gs0(kxm%n)
complex, allocatable				:: gs(:,:)
complex						:: tfld(:)
real						:: vel(:),vel1(kxm%n)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

!--------------------------------------------------
!nwb=floor(wb/img%z%d)+1
!dwb=nwb*img%z%d
!if (nwb.ne.1) then
!   vw=minval(vel(nwb-1,:))
!   write(0,*) "DATUMING WITH VELOCITY = ",vw
!end if
!--------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Compute vref
!
       vel1=1/vel
       if (present(vref1)) then
          nrefvel=count(mask = vref1(:,1).ne.-1)
          if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
          vref=1/vref1(:,1); 
          if (nrefvel==1) minv=vref(1)
       else
          if (allocated(vref)) deallocate(vref); allocate(vref(nref))
          call owvref(vel1,vref,minv,maxv)
          if (minv.eq.maxv) then
             nrefvel=1
             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
             vref=maxv
          else
             nrefvel=nref
             vref=vref
          end if
       end if

!          nrefvel=nref
!          if (maxval(vel1).ge.vsal) then
!             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel+1))
!             call owvref1(nrefvel,vel1,vref(1:nrefvel),minv,maxv)
!             vref=0;vref(nrefvel+1)=maxval(vel1)
!             nrefvel=nrefvel+1
!          else
!             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
!             vref=0
!             call owvref1(nrefvel,vel1,vref,minv,maxv)
!             if (minv.eq.maxv) then
!                nrefvel=1
!                if (allocated(vref)) then
!                   deallocate(vref)
!                   allocate(vref(nrefvel))
!                   vref=maxv
!                end if
!             end if
!          end if
!       end if

       if (nrefvel.gt.1) then
          if (allocated(gs)) deallocate(gs); allocate(gs(kxm%n,nrefvel));
          gs(:,:)=0
       else
          gs0=0
       end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  
       if (gat.eq."sou") then
          tfld=ws
          call owfftexes(tfld,inode)
          ws=tfld/sqrt(1.*kxm%n)
       else
          tfld=ws
          call owfftexer(tfld,inode)
          ws=tfld/sqrt(1.*kxm%n)
       end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!      Propagate wavefields when velocity is constant 
!
          if (nrefvel.eq.1) then
               gs0=ws
!
!      Operations in the wavenumber domain - phase shift
                do ikx = 1,kxm%n
                   sqa = rec%w%all(iw)**2 /minv**2 - kxm%all(ikx)**2
                   if (sqa > 0) then
                       gs0(ikx)=gs0(ikx)*exp(cmplx(0,sign*v%z%d*sqrt(sqa)))
                   else
                       gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
                   end if
                end do
                if (gat.eq."sou") then
                   tfld=gs0
                   call owifftexes(tfld,inode)
                   ws=tfld/sqrt(1.*kxm%n)
                else
                   tfld=gs0
                   call owifftexer(tfld,inode)
                   ws=tfld/sqrt(1.*kxm%n)
                end if
          else
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Propagate wavefields with reference velocities
!
                do iref=1,nrefvel
                    gs(:,iref)=ws
                    do ikx = 1,kxm%n
!      phase shift
                        sqa = rec%w%all(iw)**2 / vref(iref)**2 - kxm%all(ikx)**2
                        if (sqa > 0) then
                            gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,sign*v%z%d*sqrt(sqa)))
                        else
                            gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                        end if
                    end do
! Operations in x            
                    if (gat.eq."sou") then
                       tfld=gs(:,iref)
                       call owifftexes(tfld,inode)
                       gs(:,iref)=tfld/sqrt(1.*kxm%n)
                    else
                       tfld=gs(:,iref)
                       call owifftexer(tfld,inode)
                       gs(:,iref)=tfld/sqrt(1.*kxm%n)
                    end if
!      split step
                    gs(:,iref) = gs(:,iref) * exp(cmplx(0,sign*v%z%d*rec%w%all(iw))*(1./vel1-1./vref(iref)));
                 end do
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!       Interpolate between wavefields
!
!
                if (present(vref1)) then
                   do iref1=2,nrefvel
                      where (vel1 >= 1/vref1(iref1,2) .and. vel1 < 1/vref1(iref1-1,2))
                         ws=gs(:,iref1)
                      end where
                   end do
                   where (vel1 > 1/vref1(1,2))
                      ws=gs(:,1)
                   end where
                else
                   do ikx=1,kxm%n
                       indref=0
                       indref=maxloc(vref, mask = int(vref).lt.int(vel1(ikx))) 
                       if (indref(1).ne.0.and.indref(1).ne.size(vref)-1) then
                          f=(vel1(ikx) - vref(indref(1)))/vel1(ikx)
                          g=1-f
                          ws(ikx)=g*gs(ikx,indref(1))+f*gs(ikx,indref(1)+1)
                       else
                          ws(ikx)=gs(ikx,indref(1)+1)
                       end if
                   end do
                end if  
            end if 

       if (allocated(vref)) deallocate(vref) 
       if (allocated(gs)) deallocate(gs)

end subroutine owpspi

!======================================================================================
subroutine owpspi0(iwa,inode,sign,wv,vel,tfld,gat,sref)

real,optional					:: sref(:,:)
real,allocatable					:: sref1(:)
character*3					:: gat
complex						:: wv(:)
complex 						:: tfld(:),gv0(kxm%n)
integer						:: iwa,iz,ikx,iref,indref(1),ikx1,iref1
integer						:: inode,ixw,nrefvel,sign
real						:: sqa,maxv,minv,f,g,kz2,sqa1,kz2_max
real						:: vel(:),sqav(kxm%n),vel1(kxm%n)
complex, allocatable				:: gv(:,:)

       kz2_max = cos(maxang*4.*atan(1.)/180.)**2
       gv0=0
       sqav=0

       vel1=1/vel
       if (present(sref)) then
          nrefvel=count(mask = sref(:,1).ne.-1)
          if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
          sref1=sref(1:nrefvel,1)
       else


!
!          nrefvel=nref
!          if (minval(vel).le.1/vsal) then
!             if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel+1))
!             call owvref(nrefvel,1/vel(iz,:),sref1(2:nrefvel+1),minv,maxv)
!             vref=0;vref(nrefvel+1)=maxval(vel)
!             nrefvel=nrefvel+1 ! added in 10/24/08
!          else
!             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
!             vref=0
!             call owvref(nrefvel,vel(iz,:),vref,minv,maxv)
!             if (minv.eq.maxv) then
!                nrefvel=1
!                if (allocated(vref)) then
!                   deallocate(vref)
!                   allocate(vref(nrefvel))
!                   vref=maxv
!                end if
!             end if
!          end if
!
!
!
!
!

          if (allocated(sref1)) deallocate(sref1); allocate(sref1(nref))
          call owvref(vel1,sref1,minv,maxv)
          if (minv.eq.maxv) then
             nrefvel=1
             if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
             sref1=1/maxv
          else
             nrefvel=nref
             sref1=1/sref1
          end if
       end if

       if (nrefvel.gt.1) then
          if (allocated(gv)) deallocate(gv); allocate(gv(kxm%n,nrefvel));
          gv=0
       else
          gv0=0
       end if

       if (gat.eq."sou") then
          tfld=wv
          call owfftexes(tfld,inode)
          wv=tfld/sqrt(1.*kxm%n)
       else
          tfld=wv
          call owfftexer(tfld,inode)
          wv=tfld/sqrt(1.*kxm%n)
       end if

       if (nrefvel.eq.1) then
          gv0=wv
          sqav=(rec%w%all(iwa) * sref1(1))**2 - kxm%all**2
          where(sqav > 0) 
              gv0=gv0*exp(cmplx(0,sign*v%z%d*sqrt(sqav)))
          elsewhere
              gv0=gv0*exp(-sqrt(-sqav)*v%z%d)
          end where

          if (gat.eq."sou") then
             tfld=gv0
             call owifftexes(tfld,inode)
             wv=tfld/sqrt(1.*kxm%n)
          else
             tfld=gv0
             call owifftexer(tfld,inode)
             wv=tfld/sqrt(1.*kxm%n)
          end if
       else
          do iref=1,nrefvel                 
             gv(:,iref)=wv
             sqav=(rec%w%all(iwa) * sref1(iref))**2 - kxm%all**2
             where(sqav > 0) 
                 gv(:,iref)=gv(:,iref)*exp(cmplx(0,sign*v%z%d*sqrt(sqav)))
             elsewhere
                 gv(:,iref)=gv(:,iref)*exp(-sqrt(-sqav)*v%z%d)
             end where

             if (gat.eq."sou") then
                tfld=gv(:,iref)
                call owifftexes(tfld,inode)
                gv(:,iref)=tfld/sqrt(1.*kxm%n)
             else
                tfld=gv(:,iref)
                call owifftexer(tfld,inode)
                gv(:,iref)=tfld/sqrt(1.*kxm%n)
             end if

             gv(:,iref) = gv(:,iref) * exp(cmplx(0,sign*v%z%d*rec%w%all(iwa))*(vel-sref1(iref)))

          end do

          if (present(sref)) then
             do iref1=2,nrefvel
                where (vel < sref(iref1,2) .and. vel >= sref(iref1-1,2))
                   wv=gv(:,iref1)
                end where
             end do
             where (vel <= sref(1,2))
                wv=gv(:,1)
             end where
          else
             do ikx=1,kxm%n
                indref=0
                indref=maxloc(1/sref1, mask = int(1/sref1).lt.int(1/vel(ikx))) 
                if (indref(1).gt.0.and.indref(1).le.nrefvel-1) then
                   f=(1/vel(ikx) - 1/sref1(indref(1)))/(1/vel(ikx))
                   g=1-f
                   wv(ikx)=g*gv(ikx,indref(1))+f*gv(ikx,indref(1)+1)
                else if (indref(1).eq.0) then 
                   wv(ikx)=gv(ikx,1)
                else if (indref(1).ge.nrefvel-1) then 
                   wv(ikx)=gv(ikx,nrefvel)
                end if
             end do
          end if
       end if

!                   do ikx=1,kxm%n
!                       indref=0
!                       indref=maxloc(vref, mask = int(vref).lt.int(vel(iz,ikx))) 
!                       if (indref(1).ne.0.and.indref(1).ne.size(vref)-1) then
!                          f(inode)=(vel(iz,ikx) - vref(indref(1)))/vel(iz,ikx)
!                          g(inode)=1-f(inode)
!                          ws(ikx,inode)=g(inode)*gs(ikx,indref(1))+f(inode)*gs(ikx,indref(1)+1)
!                       else
!                          ws(ikx,inode)=gs(ikx,indref(1)+1)
!                       end if
!                   end do

       if (allocated(sref1)) deallocate(sref1) 
       if (allocated(gv)) deallocate(gv)
end subroutine owpspi0
!========================================================================
end module
