module OW_wei_areal

use sep
use OWvelo_areal
use OW_type
use OW_fftw_areal
use OW_parms_areal
use OW_pspi_areal


implicit none

contains
!======================================================================
subroutine owtgtwei(idip,taper,source,vel,sref)

real, optional					:: sref(:,:,:)
complex						:: source(:,:),phase
integer						:: iwa,iz,inode,nodes,i,ixw,ish,idip,sign,er
real						:: taper(:),vel(:,:),a
complex, allocatable				:: green(:,:,:,:),ws(:,:), tflds(:,:)
real, allocatable					:: sqav1(:,:)
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(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"

allocate(ws(kxm%n,nodes));                  ws=0
allocate(sqav1(kxm%n,nodes));                  sqav1=0.
allocate(green(img%xm%n,img%z%n,rec%w%n,2),stat=er);green=0
!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

!$OMP parallel do private(inode,&
                            iwa,&
                             iz,&
                           phase) ordered schedule(dynamic,1) num_threads(nodes)

   do iwa = 1,rec%w%n

       phase = cmplx(cos(pi/4.), -sin(pi/4.)) / sqrt(rec%w%all(iwa))
       inode=omp_get_thread_num()+1
       ws(:,inode)=source(:,iwa)
!!!!!!!!!!!!!!
!
! Using velocity at 
! datum to correct
! amplitudes
!
!       tflds(:,inode)=ws(:,inode)
!       call owfftexes(tflds(:,inode),inode)
!       sqav1(:,inode)=1 - kxm%all**2/((rec%w%all(iwa) * sref(1,1,numzdat))**2)
!       where(sqav1(:,inode) > 0.)
!          ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)/sqrt(sqav1(:,inode))
!       elsewhere
!          ws(:,inode)=0.
!       end where
!       tflds(:,inode)=ws(:,inode)
!       call owifftexes(tflds(:,inode),inode)
!       ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!
!
!
!!!!!!!!!!
          green(:,numzdat,iwa,1)=ws(:,inode)
          do iz=numzdat+1,nz
   
             if (present(sref)) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
             end if

             ws(:,inode)=ws(:,inode)*taper
!             !$OMP ordered
             green(:,iz,iwa,1)=ws(:,inode)!*phase
!             !$OMP end ordered

          end do

          ws(:,inode)=source(:,iwa)
          green(:,numzdat,iwa,2)=ws(:,inode)
          do iz=numzdat+1,2,-1
   
             if (present(sref)) then
                call owpspi(  iwa, inode, -1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, -1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
             end if

             ws(:,inode)=ws(:,inode)*taper
!             !$OMP ordered
             green(:,iz-1,iwa,2)=ws(:,inode)!*phase
!             !$OMP end ordered
          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

if (exist_file("greenup")) then
   call srite("greenup",green(:,:,:,2), 8*(img%xm%n*nz*rec%w%n))
end if
if (exist_file("greendown")) then
   call srite("greendown",green(:,:,:,1), 8*(img%xm%n*nz*rec%w%n))
end if

deallocate(green)
deallocate(ws,tflds)
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
222 format(a15,i4,a8,i4,a4,f6.2,a4)
end subroutine

end module
