module OW_wei_areal

use sep
use OW_type
use OW_image_areal
use OW_fftw_areal
use OW_parms_areal
use OW_pspi_areal
use OW_ffd_areal

implicit none

integer, private                             :: idip,nodes

contains
!======================================================================
subroutine owwei_init()
integer, external			:: omp_get_thread_num,omp_get_num_threads
if (node.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
else
  nodes=node
end if
write(3,*) "****************"
write(3,*) "* NUM.THREADS = ",nodes
write(3,*) "****************"

end subroutine
!======================================================================
subroutine owwei(taper,recv,source,vel,image,sref)

real, optional					:: sref(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,inode,i,ixw,ish,iz
real						:: taper(:),vel(kxm%n,img%z%n)
real						:: image(:,:,:)
real						:: GGG(kxm%n,nsecdim,nz,nodes)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes),dwr(kxm%n,nz,nodes),bws(rec%h%n,img%z%n,rec%w%n)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes),GG(kxm%n,rec%w%n)
integer, external					:: omp_get_thread_num,omp_get_num_threads

wr=0
ws=0
GGG=0
GG=0
!    write(0,*) image
!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,nodes)
tfldr=0
call owfftwr(tfldr,nodes)
if (adj) then
   !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
   do iwa = 1,rec%w%n
      inode=omp_get_thread_num()+1
      ws=0;wr=0
      ws(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=source(:,iwa)
      wr(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=recv(:,iwa)
      call owimage(.false., iwa, ws(:,inode),  wr(:,inode), GGG(:,:,nzmin,inode))
      do iz=max(2,nzmin),nzmax
         if (extrap(1:3).eq."ffd") then
            if (rec%w%all(iwa).ne.0) then
               call owffd( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou")
               call owffd( iwa, inode,-1,  wr(:,inode), vel(:,iz-1), tfldr(:,inode), "rec")
            end if
         else if (extrap(1:4).eq."pspi") then
            if (present(sref)) then
               call owpspi_init(iz)
               call owpspi( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
               call owpspi( iwa, inode,-1,  wr(:,inode), vel(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
            else
               call owpspi_init(iz)
               call owpspi( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou")
               call owpspi( iwa, inode,-1,  wr(:,inode), vel(:,iz-1), tfldr(:,inode), "rec")
            end if
         end if
         ws(:,inode)=ws(:,inode)*taper
         wr(:,inode)=wr(:,inode)*taper
         call owimage(.false., iwa, ws(:,inode),  wr(:,inode), GGG(:,:,iz-nzmin+1,inode))
      end do
      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
				" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
   do i=1,nodes
      !$OMP parallel do private(ixw,iz,ish)
      do iz=1,nz
         do ish=1,nsecdim
            do ixw=1,kxm%n
               image(ixw,ish,iz)=image(ixw,ish,iz)+GGG(ixw,ish,iz,i)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end do
   6 format(A13,i3,A8,i3,A14,i3,A12,i3)
   call owfftdestroy()
else
   !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
   do iwa = 1,rec%w%n
      inode=omp_get_thread_num()+1
      ws=0
      ws(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=source(:,iwa)
      do iz=max(2,nzmin),nzmax
         if (extrap(1:3).eq."ffd") then
            if (rec%w%all(iwa).ne.0) then
               call owffd( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou")
            end if
         else if (extrap(1:4).eq."pspi") then
            if (present(sref)) then
               call owpspi_init(iz)
               call owpspi( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
            else
               call owpspi_init(iz)
               call owpspi( iwa, inode, 1,  ws(:,inode), vel(:,iz-1), tflds(:,inode), "sou")
            end if
         end if
         ws(:,inode)=ws(:,inode)*taper
         bws(:,iz,iwa)=ws(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)
         call owimage(.true.,  iwa, ws(:,inode),  dwr(:,iz,inode), image(:,:,iz), .false.)
      end do
      ws=0.
      do iz=nzmax,max(2,nzmin),-1
         ws(:,inode)=ws(:,inode)*taper
         ws(:,inode)=ws(:,inode)+dwr(:,iz,inode)

         if (extrap(1:3).eq."ffd") then
            if (rec%w%all(iwa).ne.0) then
               call owffd( iwa, inode,1,  ws(:,inode), vel(:,iz-1), tfldr(:,inode), "rec")
            end if
         else if (extrap(1:4).eq."pspi") then
            if (present(sref)) then
               call owpspi_init(iz)
               call owpspi( iwa, inode,1,  ws(:,inode), vel(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
            else
               call owpspi_init(iz)
               call owpspi( iwa, inode,1,  ws(:,inode), vel(:,iz-1), tfldr(:,inode), "rec")
            end if
         end if
      end do
      recv(:,iwa)=ws(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)
      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
				" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
   call owfftdestroy()
end if
   if (exist_file("bws")) call srite("bws",bws,8*rec%h%n*img%z%n*rec%w%n)
end subroutine
!=============================================================
end module
