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(:),vref(:)
character*3					:: gat
complex						:: wv(:)
complex 						:: tfld(:)
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(:)
complex, allocatable				:: gv0(:),gv(:,:),tflds(:)

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

allocate(gv0(kxm%n));       gv0=0
allocate(gv(kxm%n,2))
!allocate(tflds(kxm%n));tflds=0.
!--------------------------------------------------
! Compute vref
!
     nrefvel=count(mask = sref(:,1).ne.-1)
     if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
     sref1=sref(1:nrefvel,1)

       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
          do ikx = 1,kxm%n
             sqa = (rec%w%all(iwa) * sref1(1))**2 - kxm%all(ikx)**2
             if (sqa > 0.) then
                gv0(ikx)=gv0(ikx)*exp(cmplx(0,sign*v%z%d*sqrt(sqa)))
             else
                gv0(ikx)=gv0(ikx)*exp(-sqrt(-sqa)*v%z%d)
             end if 
          end do

          tfld=gv0
          call owifftexes(tfld,inode)
          wv=tfld/sqrt(1.*kxm%n)

       else
          do iref=1,nrefvel                 
             gv(:,iref)=wv
             do ikx = 1,kxm%n
                sqa = (rec%w%all(iwa) * sref1(iref))**2 - kxm%all(ikx)**2;
                if (sqa > 0) then
                   gv(ikx,iref)=gv(ikx,iref)*exp(cmplx(0,sign*v%z%d*sqrt(sqa)))
                else
                   gv(ikx,iref)=gv(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                end if
             end do

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

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

          end do

          do iref1=2,nrefvel
             where (vel >= 1/sref(iref1,2) .and. vel < 1/sref(iref1-1,2))
                wv=gv(:,iref1)
             end where
          end do
          where (vel > 1/sref(1,2))
             wv=gv(:,1)
          end where
       end if

deallocate(gv0,gv,sref1)!,tfldr,tflds)

end subroutine owpspi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module
