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
integer						:: inode,ixw,nrefvel,sign
real						:: sqa,maxv,minv
real						:: vel(:)
complex, allocatable				:: gv0(:),gv(:,:)

allocate(gv0(kxm%n));       gv0=0
allocate(gv(kxm%n,2))
!allocate(tflds(kxm%n),tfldr(kxm%n)); tflds=0;tfldr=0.
!--------------------------------------------------
! Compute vref
!
       if (present(sref)) then
          nrefvel=count(mask = sref.ne.-1)
          if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
          sref1=sref
       else
          nrefvel=nref
          if (maxval(vel).ge.vsal) then
             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel+1))
             if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel+1))
             vref=0;sref1=0
             call owvref(nrefvel,vsal,vel,vref(1:nrefvel),minv,maxv)
             vref(nrefvel+1)=maxval(vel)
             sref1=1/vref
          else
             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
             if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
             vref=0
             call owvref(nrefvel,vsal,vel,vref,minv,maxv)
             sref1=1/vref
          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.or.maxv.eq.minv) 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

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

             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))*(1./vel-sref1(iref)))

          end do
          do iref=1,nrefvel-1
             where (vel >= 1/sref1(iref) .and. vel < 1/sref1(iref+1))
                wv=gv(:,iref)
             end where
          end do
          where (vel >= 1/sref1(nrefvel))
             wv=gv(:,nrefvel)
          end where 
       end if

deallocate(gv0,gv)!,tfldr,tflds)
!========================================================================||
!  Datuming
!========================================================================||
!    call owimage0(iwa,ws,wr,GG(iz,:,:))
! 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
!========================================================================
!========================================================================

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