module OW_pspi_areal2

use OW_type
use OW_parms_areal
use OWvelo_areal
use OW_image_areal

implicit none

private						:: ni,ixs,nwb,vsal
integer						:: ni,ixs,nwb
real						:: wb

contains

!---------------------------------------------------------
subroutine owpspi_init(ni_in,ixs_in,wb_in)

integer						:: ni_in,ixs_in
real						:: wb_in

ni=ni_in
ixs=ixs_in
wb=wb_in
wb=floor(wb/img%z%d)*img%z%d
nwb=floor(wb/img%z%d)+1

end subroutine
!---------------------------------------------------------
subroutine owpspi(taper,recv,source,vsal,vel,G1)

!optional					:: rwr
integer						:: iw,inode,ino,i1,nout,node
integer						:: iz,ikx,iref,ixw,ish
integer*8,allocatable	   			:: planF(:),planI(:)
real						:: sqa,f,g
real						:: minv, maxv
real						:: taper(:),vel(:,:),vsal
complex						:: G1(:,:,:),recv(:,:),source(:,:,:)
complex, allocatable				:: GG(:,:,:,:),rwr(:,:),rws(:,:,:)
real, allocatable				:: vref(:,:),sqa0(:,:)
complex, allocatable				:: ws(:,:),wr(:,:),gs0(:,:),gr0(:,:)
complex, allocatable				:: gs(:,:,:),gr(:,:,:),tfld(:,:)

integer, external 				:: omp_get_thread_num,omp_get_num_threads,omp_get_num_procs,omp_set_num_threads

integer, parameter   				:: FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
integer, parameter   				:: FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
integer, parameter   				:: FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
integer, parameter   				:: FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
integer, parameter   				:: FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
integer, parameter   				:: FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
integer, parameter   				:: FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
integer, parameter   				:: FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536

node=omp_get_num_procs();node=1
allocate(planF(node));	planF=0
allocate(planI(node));	planI=0

if (allocated(GG))	deallocate(GG);		allocate(GG(img%z%n,img%sh%n,img%xm%n,node));	GG=0
if (allocated(wr))	deallocate(wr);		allocate(wr(img%xm%n,node));			wr=0
if (allocated(ws))	deallocate(ws);		allocate(ws(img%xm%n,node));			ws=0
if (allocated(gs))	deallocate(gs);		allocate(gs(img%xm%n,nref,node));		gs=0
if (allocated(gr))	deallocate(gr);		allocate(gr(img%xm%n,nref,node));		gr=0
if (allocated(gs0))	deallocate(gs0);	allocate(gs0(img%xm%n,node));			gs0=0
if (allocated(gr0))	deallocate(gr0);	allocate(gr0(img%xm%n,node));			gr0=0
if (allocated(sqa0))	deallocate(sqa0);	allocate(sqa0(img%xm%n,node));			sqa0=0
if (allocated(vref))	deallocate(vref);	allocate(vref(nref,node));			vref=0
if (allocated(tfld))	deallocate(tfld);	allocate(tfld(img%xm%n,node));			tfld=0

do i1=1,node
	call sfftw_plan_dft(planF(i1),1,kxm%n,tfld(:,i1),tfld(:,i1),fftW_FORWARD,fftW_PATIENT)
	call sfftw_plan_dft(planI(i1),1,kxm%n,tfld(:,i1),tfld(:,i1),fftw_BACKWARD,fftW_PATIENT)
end do

! $OMP parallel default(shared) private(inode)

        inode=omp_get_thread_num()+1;

! $OMP do private(iw,iz,ikx,minv,maxv,sqa,iref)

do iw = 1,rec%w%n
!kxm%n=700

!write(0,*) iw,inode

	do ikx=1,kxm%n
		ws(ikx,inode)=source(ikx,iw,1) 
		wr(ikx,inode)=recv(ikx,iw) 
	end do

!write(0,*) "Passei 01"

        depth:do iz=1,img%z%n

	    ! Compute vref
            call owvref(nref,vsal,vel(iz,:),vref(:,inode),minv,maxv)

	do ikx=1,kxm%n
            tfld(ikx,inode)=ws(ikx,inode)
	    call sfftw_execute(planF(inode));
	    ws(ikx,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
            tfld(ikx,inode)=wr(ikx,inode)
	    call sfftw_execute(planF(inode))
	    wr(ikx,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
	end do

!write(0,*) "Passei 02"

           vxz:if (minv .eq. maxv) then
		! Operations in the wavenumber domain
		! phase shift

                wavenumber1:do ikx = 1,kxm%n
                   gs0(ikx,inode)=ws(ikx,inode);  gr0(ikx,inode)=wr(ikx,inode);

                   sqa = rec%w%all(iw)**2 / minv**2 - kxm%all(ikx)**2

                   evanescent1:if (sqa > 0) then
                       gs0(ikx,inode)=gs0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                       gr0(ikx,inode)=gr0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
                   else
                       gs0(ikx,inode)=gs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                       gr0(ikx,inode)=gr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                   end if evanescent1
                end do wavenumber1

!write(0,*) "Passei 03"

		do ikx=1,kxm%n
	                tfld(ikx,inode)=gs0(ikx,inode)
			call sfftw_execute(planI(inode))
			ws(ikx,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
	                tfld(ikx,inode)=gr0(ikx,inode)
			call sfftw_execute(planI(inode))
			wr(ikx,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
		end do
 
!write(0,*) "Passei 04"

           else
                velref: do iref=1,nref                 
			do ikx=1,kxm%n
	                    gs(ikx,iref,inode)=ws(ikx,inode)
			    gr(ikx,iref,inode)=wr(ikx,inode)
			end do

!write(0,*) "Passei 05"

                    wavenumber2: do ikx = 1,kxm%n
			! phase shift
                        sqa = rec%w%all(iw)**2 / vref(iref,inode)**2 - kxm%all(ikx)**2;
                        evanescent2: if (sqa > 0) then
                            gs(ikx,iref,inode)=gs(ikx,iref,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                            gr(ikx,iref,inode)=gr(ikx,iref,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                        else
                            gs(ikx,iref,inode)=gs(ikx,iref,inode)*exp(-sqrt(-sqa)*v%z%d)
                            gr(ikx,iref,inode)=gr(ikx,iref,inode)*exp(-sqrt(-sqa)*v%z%d)
                        end if evanescent2
                    end do wavenumber2

!write(0,*) "Passei 06"

		      ! Operations in x            
			do ikx=1,kxm%n
	                      tfld(ikx,inode)=gs(ikx,iref,inode)
			      call sfftw_execute(planI(inode))
			      gs(ikx,iref,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
	                      tfld(ikx,inode)=gr(ikx,iref,inode)
			      call sfftw_execute(planI(inode))
			      gr(ikx,iref,inode)=tfld(ikx,inode)/sqrt(1.*kxm%n)
			end do

!write(0,*) "Passei 07"

		      ! split step
			do ikx=1,kxm%n
	                      gs(ikx,iref,inode) = gs(ikx,iref,inode) * exp(cmplx(0,v%z%d*rec%w%all(iw))*(1./vel(iz,ikx)-1./vref(iref,inode)))
	                      gr(ikx,iref,inode) = gr(ikx,iref,inode) * exp(cmplx(0,-v%z%d*rec%w%all(iw))*(1./vel(iz,ikx)-1./vref(iref,inode)))
			end do

!write(0,*) "Passei 08"

                end do velref

                do iref=1,nref-1
		   do ikx=1,kxm%n
                     if (vel(iz,ikx) >= vref(iref,inode) .and. vel(iz,ikx) < vref(iref+1,inode)) then
                	     ws(ikx,inode)=gs(ikx,iref,inode)
        	             wr(ikx,inode)=gr(ikx,iref,inode)
                     elseif (vel(iz,ikx) >= vref(nref,inode)) then
                             ws(ikx,inode)=gs(ikx,nref,inode)
                             wr(ikx,inode)=gr(ikx,nref,inode)
		     end if
		   end do
                end do

           end if vxz
	   do ikx=1,kxm%n
	           ws(ikx,inode)=ws(ikx,inode)*taper(ikx)
	           wr(ikx,inode)=wr(ikx,inode)*taper(ikx)
	   end do
           call owimage(iw,ws(:,inode),wr(:,inode),GG(iz,:,:,inode))

! $OMP critical
	   do ikx=1,kxm%n
		do ish=1,img%sh%n 
	             G1(iz,ish,ikx)=G1(iz,ish,ikx)+GG(iz,ish,ikx,inode)
		end do
	   end do
! $OMP end critical   

        end do depth

        if (mod(iw,10)==1) write(0,5) " FREQUENCY # ",iw," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, &
					" FROM NODE: ",inode

end do

!! $OMP end parallel do

! $OMP End do
! $OMP end parallel

5 format(A13,i3,A8,i3,A14,i3,A12,i3)

end subroutine owpspi
end module
