module OW_pspi_areal

use OW_type
use OW_parms_areal
use OW_wei_areal
use OW_image_areal

implicit none

!public						:: owpspi_init, owpspi
private						:: ni,ixs,nwb,node,vsal

integer						:: ni,ixs,nwb,node
real						:: wb

contains

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

integer						:: ni_in,ixs_in,node_in
real						:: wb_in

ni=ni_in
ixs=ixs_in
node=node_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)

integer					:: iw,inode,ino,i1,i2,nout,iwind
real					:: vel(:,:),vsal,taper(:)
complex					:: recv(:,:),source(:,:,:),G1(:,:,:)
complex, allocatable			:: GG(:,:,:,:),GGr(:,:,:),GGs(:,:,:),GG1(:,:,:),GG2(:,:,:),G2(:,:,:)

integer, external			:: omp_get_thread_num

!allocate(GG(img%z%n,img%sh%n,img%xm%n,node))
!if (allocated(GG1)) deallocate(GG1); allocate(GG1(img%z%n,img%xm%n,rec%w%n)); GG1=0
!if (allocated(GG2)) deallocate(GG2); allocate(GG2(img%z%n,img%xm%n,rec%w%n)); GG2=0

!write(0,*) size(GG,4)

!if (exist_file("wr")) then
!   if (allocated(rwr)) deallocate(rwr);allocate(rwr(img%z%n,rec%h%n));		rwr=0
!   if (allocated(rws)) deallocate(rws);allocate(rws(rec%w%n,img%z%n,rec%h%n));		rws=0
!   write(0,*) "rwr and rws allocated"
!end if

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

       inode=omp_get_thread_num()+1;

! $OMP do private(iw)
do iw = 1,rec%w%n
        if (allocated(GGr)) deallocate(GGr); allocate(GGr(img%z%n,img%xm%n,rec%w%n)); GGr=0
        if (allocated(GGs)) deallocate(GGs); allocate(GGs(img%z%n,img%xm%n,rec%w%n)); GGs=0

!        if (allocated(GG1)) deallocate(GG1); allocate(GG1(img%z%n,img%xm%n)); GG1=0
!        if (allocated(GG2)) deallocate(GG2); allocate(GG2(img%z%n,img%xm%n)); GG2=0

!        if (allocated(G2)) deallocate(G2); allocate(G2(img%z%n,img%sh%n,img%xm%n)); G2=0

        call owwei_init(ni,iw,wb,nwb,node,inode)
!        call owwei(taper,recv,source,vsal,vel,GG(:,:,:,inode),rwr)
        call owwei(taper,recv,source,vsal,vel,GGr(:,:,iw),GGs(:,:,iw))

!        do iwind=1,2*n_wind-1
!           if (iwind <= n_wind) then
!              GG1=GGr*windf(iwind)
!!              GG2=GGs*windf(iwind)
!           else
!              GG1=GGr*windf(iwind-n_wind+1)
!!              GG2=GGs*windf(iwind-n_wind+1)
!           end if
!!           call owimage1(GG1,GGs,G2)
!           G1=G1+G2
!        end do

           
!         if (exist_file("wr")) then
!           rws(iw,:,:)=rws(iw,:,:)+rwr
!        end if

! $OMP critical
!        G1(:,:,:)=G1(:,:,:)+GG(:,:,:,inode)
! $OMP end critical   
        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 do
! $OMP end parallel

if (allocated(GG1)) deallocate(GG1); allocate(GG1(img%z%n,img%xm%n,(rec%w%n+n_wind-1))); GG1=0
if (allocated(GG2)) deallocate(GG2); allocate(GG2(img%z%n,img%xm%n,(rec%w%n+n_wind-1))); GG2=0

! GG1=conv_fft(GGr,windf)
! GG2=conv_fft(GGs,windf)

!deallocate (GGr,GGs)

if (allocated(G2)) deallocate(G2); allocate(G2(img%z%n,img%sh%n,img%xm%n)); G2=0

do iw = 1,rec%w%n
   do i1=1,n_wind
      i2=iw+i1-1
      GG1(:,:,i2)=GG1(:,:,i2)+(GGr(:,:,iw))*windf(i1)
      GG2(:,:,i2)=GG2(:,:,i2)+(GGs(:,:,iw))*windf(i1)
   end do
end do

GG1(:,:,1:n_wind)=cmplx(real(GG1(:,:,1:n_wind)),aimag(GGr))
GG2(:,:,1:n_wind)=cmplx(real(GG2(:,:,1:n_wind)),aimag(GGs))


do iw = 1,rec%w%n
   call owimage1(GG1(:,:,iw),GG2(:,:,iw),G2)
   G1=G1+G2
end do

if (exist_file("wr")) call srite("wr",GG1,4*img%z%n*img%xm%n)
if (exist_file("ws")) call srite("ws",GG2,4*img%z%n*img%xm%n)

!if (exist_file("wr")) then
!   nout=8*rec%h%n*img%z%n*rec%w%n
!   write(0,*) "writing rws" 
!    call srite("wr",rws,nout)
!end if

!   if (iw < n_wind .and. (rec%w%n - iw + 1) < n_wind) then
!      m_wind=n_wind/2+iw
!   else


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

end subroutine owpspi


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function conv_fft(G,windf) result(temp1)

 complex, allocatable		:: temp1(:,:,:),temp2(:),tfld(:)
 complex, intent(in)		:: G(:,:,:)
 real 				:: windf(:)
 integer				:: ix, iz, ih


 allocate(temp1(img%z%n,img%xm%n,rec%w%n),temp2(rec%w%n),tfld(rec%w%n))
 call owfftw2(tfld)
 temp1=cmplx(0,0)
 temp2=cmplx(windf,0)
 tfld=temp2; call sfftw_execute(plan2F); temp2=tfld/sqrt(1.*rec%w%n)
 do iz=1,img%z%n
    do ih=1,img%xm%n
       tfld=G(iz,ih,:)
       call sfftw_execute(plan2F)
       temp1(iz,ih,:)=tfld/sqrt(1.*rec%w%n)*temp2
   end do
 end do
 
 do iz=1,img%z%n
    do ih=1,img%xm%n
       tfld=temp1(iz,ih,:); call sfftw_execute(plan2I); temp1(iz,ih,:)=tfld/sqrt(1.*rec%w%n)
   end do
 end do

 deallocate(temp2,tfld)

 end function conv_fft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


end module
