module OW_pspi_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal
use OWvelo_areal

implicit none

contains

!---------------------------------------------------------
subroutine owpspi(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)
complex, allocatable				:: gv(:,:)

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

       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
          if (allocated(sref1)) deallocate(sref1); allocate(sref1(nref))
          call owvref(nref,vsal,1/vel,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=1 - kxm%all**2/((rec%w%all(iwa) * sref1(1))**2)
          where(sqav >=0 .and. sqav < kz2_max) 
              gv0=0.
          elsewhere(sqav >= kz2_max) 
              gv0=gv0*exp(cmplx(0,sign*v%z%d*sqrt(sqav*(rec%w%all(iwa)*sref1(1))**2)))
          elsewhere
              gv0=gv0*exp(-sqrt(-sqav*(rec%w%all(iwa) * sref1(1))**2)*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=1 - kxm%all**2/((rec%w%all(iwa) * sref1(iref))**2)
             where(sqav >=0 .and. sqav < kz2_max) 
                 gv(:,iref)=0.
             elsewhere(sqav >= kz2_max) 
                 gv(:,iref)=gv(:,iref)*exp(cmplx(0,sign*v%z%d*sqrt(sqav*(rec%w%all(iwa)*sref1(iref))**2)))
             elsewhere
                 gv(:,iref)=gv(:,iref)*exp(-sqrt(-sqav*(rec%w%all(iwa) * sref1(iref))**2)*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).lt.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

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