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(:,:)

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

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(:,1).ne.-1)
          if (allocated(sref1)) deallocate(sref1); allocate(sref1(nrefvel))
          sref1=sref(1:nrefvel,1)
          if (iwa==1) write(2,*) nrefvel,1/sref1
!       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
!             if(minv==maxv) nrefvel=1
!          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)
!             if (minv==maxv) nrefvel=1
!             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) then
          gv0=wv
          do ikx = 1,kxm%n
             sqa = (rec%w%all(iwa) * sref1(1))**2 - kxm%all(ikx)**2

                 kz2 = sqa/((rec%w%all(iwa) * sref1(1))**2)
                 if (kz2 >= 0. .and. kz2 < kz2_max) then
                    gv0(ikx)=0.
                 else if (kz2 >= kz2_max) then
                    sqa1=sqrt(sqa)*(rec%w%all(iwa) * sref1(1))
                    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 

!             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

!          if (present(sref)) then
             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
!          else
!             do ikx1=1,kxm%n
!                indref=0
!                indref=maxloc(1/sref1, mask = int(1/sref1).lt.int(vel(ikx1))) 
!                if (indref(1).ne.0.and.indref(1).ne.size(1/sref1)-1) then
!                   f=(vel(ikx1) - 1/sref1(indref(1)))/vel(ikx1)
!                   g=1-f
!                   wv(ikx1)=g*gv(ikx1,indref(1))+f*gv(ikx1,indref(1)+1)
!                else
!                   wv(ikx1)=gv(ikx1,indref(1)+1)
!                end if
!             end do
!          end if
       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
