module OW_pspi_areal

use OW_type
use OW_parms_areal
use OW_wei_areal
use OW_image_areal
use OW_exec_time

implicit none

private						:: ni,ixs,node,vsal
integer						:: ni,ixs,node

contains

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

integer						:: ni_in,ixs_in,node_in

ni=ni_in
ixs=ixs_in
node=node_in

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(:,:,:,:),G2(:,:,:)
real					:: start, finish

integer, external				:: omp_get_thread_num

!allocate(GG(img%z%n,img%sh%n,img%xm%n,node))

allocate(G2(img%z%n,img%sh%n,img%xm%n));G2=0.

! $OMP parallel default(shared) private(inode) ! num_threads(8)
       inode=omp_get_thread_num()+1;
! $OMP do private(iw)

   call owexec_time(start)
   do iw = 1,rec%w%n
       call owwei_init(ni,iw,node,inode)
       call owwei(taper,recv(:,iw),source(:,iw),vsal,vel,G2)
       G1=G1+G2

! $OMP critical
!      G1(:,:,:)=G1(:,:,:)+GG(:,:,:,inode)
! $OMP end critical   
       if (mod(iw,10)==1) write(0,6) " FREQUENCY # ",iw," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, &
					" FROM NODE: ",inode
   end do

   call owexec_time(finish)

! $OMP End do
! $OMP end parallel

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

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

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

integer, external				:: omp_get_thread_num

!allocate(GG(img%z%n,img%sh%n,img%xm%n,node))

! $OMP parallel default(shared) private(inode) ! num_threads(8)
       inode=omp_get_thread_num()+1;
! $OMP do private(iw)

   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

   call owexec_time(start)
   do iw = 1,rec%w%n
       call owwei_init(ni,iw,node,inode)
       call owwei1(taper,recv(:,iw),source(:,iw),vsal,vel,GGr(:,:,iw),GGs(:,:,iw))

! $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
   call owexec_time(finish)
   write(0,'(8x,a,f10.2,a)')' Propagating wavefields:', (finish - start),' (sec)'

! $OMP End do
! $OMP end parallel

   call owexec_time(start)
   if (allocated(GG1)) deallocate(GG1); allocate(GG1(img%z%n,img%xm%n,2*n_wind-1)); GG1=0 
   call tr_fft(GGr,GG1,.false.);deallocate (GGr)
   if (allocated(GG2)) deallocate(GG2); allocate(GG2(img%z%n,img%xm%n,2*n_wind-1)); GG2=0 
   call tr_fft(GGs,GG2,.false.);deallocate (GGs)
   call owexec_time(finish)
   write(0,'(8x,a,f10.2,a)')' Fourier Transform:', (finish - start),' (sec)'

   if (exist_file("wr") .and. ixs.eq.1) call srite("wr",GG1,4*img%z%n*img%xm%n*(2*n_wind-1))
   if (exist_file("ws") .and. ixs.eq.1) call srite("ws",GG2,4*img%z%n*img%xm%n*(2*n_wind-1))

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

   if (not(stereo)) then
      call owexec_time(start)
      do iw = 1,2*n_wind-1
         call owimage2(GG1(:,:,iw+itzero),GG2(:,:,iw+itzero),G2)
         G1=G1+G2
      end do
      call owexec_time(finish)
      write(0,'(8x,a,f10.2,a)')' Imaging condition:', (finish - start),' (sec)'
      deallocate (GG1,GG2)
   else
      call owexec_time(start)
      do iw = 1,2*n_wind-1
         call owimage_stereo(GG1(:,:,iw+itzero),GG2(:,:,iw+itzero),G2)
         G1=G1+G2
      end do
      call owexec_time(finish)
      write(0,'(8x,a,f10.2,a)')' Imaging condition:', (finish - start),' (sec)'
      deallocate (GG1,GG2)
end if

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

end subroutine owpspi1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine tr_fft(G,H,upward)

 logical				:: upward
 complex				:: G(:,:,:)
 real				:: H(:,:,:) 
 real, allocatable			:: temp(:,:,:)
 complex, allocatable		:: tfld(:)
 integer				:: iz, ix, iw

 allocate(tfld(n_sampt))
 call owfftw2(tfld)
 do ix=1,img%xm%n
    do iz=1,img%z%n
       tfld=0; 
       tfld(1:rec%w%n)=G(iz,ix,:)
       call sfftw_execute(plan2I)
       H(iz,ix,n_wind:2*n_wind-1)=real(tfld(1:n_wind))
       do iw=1,n_wind-1
          H(iz,ix,iw)=real(tfld(n_sampt-iw+1))
       end do
   end do
 end do
 if (upward) then
    allocate(temp(img%z%n,img%xm%n,2*n_wind-1))
    do iw=1,2*n_wind-1
       temp(:,:,iw)=H(:,:,2*n_wind-iw)
    end do
    H=temp
    deallocate(temp)
 end if 
 deallocate(tfld)

 end subroutine tr_fft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module
