module OW_pspi_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal
use OWvelo_areal

implicit none

contains

!---------------------------------------------------------
subroutine owpspi(inode,sign,ww,ws,vel,tflds,vref1) 

real,optional					:: vref1(:,:)
integer						:: ikx,iref,iref1,ixw,i,j
integer						:: nwb,inode,nrefvel,sign
integer						:: indref(1)
real						:: sqa,a,dwb,f,g,sqav,sq,ww
real						:: minv, maxv, dvel
real, allocatable					:: vref(:)
complex						:: ws(:,:)
complex, allocatable				:: gs(:,:,:)
complex						:: tflds(:,:)
real						:: vel(:,:),vel1(kxm%n,kym%n)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 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(vel,vref,minv,maxv,dvel)
          if (int(minv).eq.int(maxv)) then
             nrefvel=1
             if (allocated(vref)) deallocate(vref); allocate(vref(nrefvel))
             vref=maxv
          else
             nrefvel=nref
          end if
       end if
       if (allocated(gs)) deallocate(gs); allocate(gs(kxm%n,kym%n,nrefvel));gs=0

!write(0,*) minval(vref),maxval(vref),nrefvel, minval(vel(:,1)),maxval(vel(:,1))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Propagate wavefields with reference velocities
!

       do iref=1,nrefvel
          do i=1,kym%n;do j=1,kxm%n
             sq = kxm%all(j)**2 + kym%all(i)**2
             gs(j,i,iref)=ws(j,i)
             sqav=(ww / vref(iref))**2 - sq
             if (sqav >= 0) then
                 gs(j,i,iref)=gs(j,i,iref)*exp(cmplx(0,sign*v%z%d*sqrt(sqav)))
             else
                 gs(j,i,iref)=0
             end if
          end do;end do
          tflds=gs(:,:,iref)
          call owifftexes(inode)
          gs(:,:,iref)=tflds/sqrt(1.*kxm%n*kym%n)
!      split step
          gs(:,:,iref) = gs(:,:,iref) * exp(cmplx(0,sign*v%z%d*ww)*(vel-1./vref(iref)));
       end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Interpolate wavefields
!
       if (nrefvel.eq.1) then
          ws=gs(:,:,1)
       else
          do iref1=2,nrefvel
             do i=1,kym%n;do j=1,kxm%n
!                if (vel1(j,i) >= 1/vref1(iref1,2) .and. vel1(j,i) < 1/vref1(iref1-1,2)) then
!                   ws(j,i)=gs(j,i,iref1)
!                else if (vel1(j,i) > 1/vref1(1,2)) then
!                   ws(j,i)=gs(j,i,1)
!                end if
                if (vel(j,i) >= 1/vref1(iref1,2) .and. vel(j,i) < 1/vref1(iref1-1,2)) then
                   ws(j,i)=gs(j,i,iref1)
                else if (vel(j,i) > 1/vref1(1,2)) then
                   ws(j,i)=gs(j,i,1)
                end if
             end do;end do
          end do
       end if

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

end subroutine owpspi

!======================================================================================
end module
