module OW_wei_areal

use sep
use OWvelo_areal
use OW_type
use OW_image_areal
use OW_fftw_areal
use OW_parms_areal
use OW_pspi_areal
use OW_phsscr_areal
use OW_cstret

implicit none

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

real, optional					:: sref(:,:,:)
complex						:: green(:,:,:),source(:,:)
integer						:: iwa,iz,inode,nodes,i,ixw,ish,idip,sign
real						:: taper(:),vel(:,:),a
complex, allocatable				:: ws(:,:), tflds(:,:)
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

!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

!$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(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=source(:,iwa)

       if (sign==1) then
          green(:,numzdat,iwa)=ws(:,inode)
          do iz=numzdat+1,nzmax
   
             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
             green(:,iz,iwa)=ws(:,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
       else
          do iz=numzdat,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
             green(:,iz-1,iwa)=ws(:,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 if
   end do

!$OMP END PARALLEL DO

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

end subroutine

!======================================================================
subroutine owwemva(adj,taper,recv,source,vel,deltaS,deltaI,sref)

real, optional					:: sref(:,:,:)
logical						:: adj
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,inode,nodes,i,ixw,ish
real						:: taper(:),vel(:,:)
real						:: deltaS(:,:)
real						:: deltaI(:,:,:)
real, allocatable					:: GGG(:,:,:,:)
!complex						:: deltaI(:,:,:)
real, allocatable					:: GG(:,:,:,:)
!complex, allocatable				:: GG(:,:,:,:)
complex, allocatable				:: ws(:,:),wr(:,:)
complex, allocatable				:: aws(:,:),awr(:,:)
complex, allocatable				:: ddws(:,:),ddwr(:,:)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:)
complex, allocatable				:: bws(:,:,:),bwr(:,:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
complex, allocatable				:: dws_z(:,:),dwr_z(:,:)
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(wr(kxm%n,nodes));                  wr=0
allocate(ws(kxm%n,nodes));                  ws=0
allocate(ddwr(kxm%n,nodes));                ddwr=0
allocate(ddws(kxm%n,nodes));                ddws=0
allocate(bwr(img%z%n,kxm%n,nodes));         bwr=0
allocate(bws(img%z%n,kxm%n,nodes));         bws=0

!allocate(awr(kxm%n,nodes));                 awr=0
!allocate(aws(kxm%n,nodes));                 aws=0


!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

allocate(tfldr(kxm%n,nodes));tfldr=0
call owfftwr(tfldr,nodes)

if (adj) then

   allocate(GGG(img%z%n,1,kxm%n,nodes)); GGG=0
   allocate(dwr_z(kxm%n,nodes));         dwr_z=0
   allocate(dws_z(kxm%n,nodes));         dws_z=0
   allocate(dwr(img%z%n,kxm%n,nodes));   dwr=0
   allocate(dws(img%z%n,kxm%n,nodes));   dws=0

   !$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(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=source(:,iwa)
       wr(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=recv(:,iwa)
       dws(:,:,inode)=0.
       dwr(:,:,inode)=0.
       bws(:,:,inode)=0.
       bwr(:,:,inode)=0.
       dws_z(:,inode)=0.
       dwr_z(:,inode)=0.
       ddwr(:,inode)=0
       ddws(:,inode)=0

       bws(1,:,inode)=ws(:,inode)
       bwr(1,:,inode)=wr(:,inode)
       
       do iz=nwbott+1,img%z%n

          if (present(sref)) then
             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(  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
          ws(:,inode)=ws(:,inode)*taper
          wr(:,inode)=wr(:,inode)*taper
          bws(iz,:,inode)=ws(:,inode)
          bwr(iz,:,inode)=wr(:,inode)

       end do

       do iz=1,img%z%n
          call owimage(.true., .false., iwa, bws(iz,:,inode),  dwr(iz,:,inode), deltaI(iz,:,:))
          call owimage(.true., .true.,  iwa, dws(iz,:,inode),  bwr(iz,:,inode), deltaI(iz,:,:))
       end do

       do iz=img%z%n,nwbott+1,-1
          dws_z(:,inode)=dws_z(:,inode)+dws(iz,:,inode)
          dws_z(:,inode)=dws_z(:,inode)*taper
          if(present(sref)) then
             call owpspi(  iwa, inode, -1, dws_z(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
          else
             call owpspi(  iwa, inode, -1, dws_z(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
          end if
!          ddws(:,inode)=0
          call owphsscr(iwa, 1,ddws(:,inode),bws(iz-1,:,inode),vel(iz-1,:))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          tflds(:,inode)=bws(iz-1,:,inode)
!          call owfftexes(tflds(:,inode),inode)
!          aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          call owphsscr(iwa,inode,+1,ddws(:,inode),aws(:,inode),vel(iz-1,:),tflds(:,inode),"sou")
          GGG(iz-1,1,:,inode)=GGG(iz-1,1,:,inode)+conjg(ddws(:,inode))*dws_z(:,inode)

          dwr_z(:,inode)=dwr_z(:,inode)+dwr(iz,:,inode)
          dwr_z(:,inode)=dwr_z(:,inode)*taper
          if(present(sref)) then
             call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel(iz-1,:), tfldr(:,inode), "rec", sref(:,:,iz-1))
          else
             call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel(iz-1,:), tfldr(:,inode), "rec")
          end if
!          ddwr(:,inode)=0
          call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz-1,:,inode),vel(iz-1,:))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          tfldr(:,inode)=bwr(iz-1,:,inode)
!          call owfftexer(tfldr(:,inode),inode)
!          awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          call owphsscr(iwa,inode,-1,ddwr(:,inode),awr(:,inode),vel(iz-1,:),tfldr(:,inode),"rec")
          GGG(iz-1,1,:,inode)=GGG(iz-1,1,:,inode)+conjg(ddwr(:,inode))*dwr_z(:,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

   call owfftdestroy()
   deallocate(tfldr,tflds)

   do i=1,nodes
      !$OMP parallel do private(ixw,iz)
      do iz=1,img%z%n
         do ixw=1,kxm%n
            deltaS(iz,ixw)=deltaS(iz,ixw)+GGG(iz,1,ixw,i)
         end do
      end do
      !$OMP END PARALLEL DO
   end do

   deallocate(ddwr,ddws,dwr_z,dws_z)
   deallocate(wr,ws,dwr,dws,GGG,bwr,bws)

else

!   allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes));          GGG=0
   allocate(GG(img%z%n,img%sh%n,kxm%n,nodes));           GG=0
   allocate(dwr(1,kxm%n,nodes));                         dwr=0
   allocate(dws(1,kxm%n,nodes));                         dws=0

   !$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(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=source(:,iwa)
       wr(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=recv(:,iwa)
       dws(1,:,inode)=0
       dwr(1,:,inode)=0

       do iz=nwbott+1,img%z%n

          if (present(sref)) then
             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(  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
          ws(:,inode)=ws(:,inode)*taper
          wr(:,inode)=wr(:,inode)*taper
          bws(iz,:,inode)=ws(:,inode)
          bwr(iz,:,inode)=wr(:,inode)

       end do

       do iz=nwbott+1,img%z%n

          ddws(:,inode)=0
          call owphsscr(iwa, 1,ddws(:,inode),bws(iz,:,inode),vel(iz-1,:))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          tflds(:,inode)=bws(iz,:,inode)
!          call owfftexes(tflds(:,inode),inode)
!          aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          call owphsscr(iwa,inode,+1,ddws(:,inode),aws(:,inode),vel(iz-1,:),tflds(:,inode),"sou")
          dws(1,:,inode)=dws(1,:,inode)+ddws(:,inode)*deltaS(iz-1,:)

          ddwr(:,inode)=0
          call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz,:,inode),vel(iz-1,:))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          tfldr(:,inode)=bwr(iz,:,inode)
!          call owfftexer(tfldr(:,inode),inode)
!          awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          call owphsscr(iwa,inode,-1,ddwr(:,inode),awr(:,inode),vel(iz-1,:),tfldr(:,inode),"rec")
          dwr(1,:,inode)=dwr(1,:,inode)+ddwr(:,inode)*deltaS(iz-1,:)

          if (present(sref)) then
             call owpspi(  iwa, inode, 1, dws(1,:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
             call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel(iz-1,:), tfldr(:,inode), "rec", sref(:,:,iz-1))
          else
             call owpspi(  iwa, inode, 1, dws(1,:,inode), vel(iz-1,:), tflds(:,inode), "sou")
             call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel(iz-1,:), tfldr(:,inode), "rec")
          end if

          dws(1,:,inode)=dws(1,:,inode)*taper
          dwr(1,:,inode)=dwr(1,:,inode)*taper
 
!          call owimage(.false., .false., iwa, dws(1,:,inode), bwr(iz,:,inode), GGG(iz,:,:,inode))
!          call owimage(.false., .false., iwa, bws(iz,:,inode), dwr(1,:,inode), GGG(iz,:,:,inode))
          call owimage(.false., .false., iwa, dws(1,:,inode), bwr(iz,:,inode), GG(iz,:,:,inode))
          call owimage(.false., .false., iwa, bws(iz,:,inode), dwr(1,:,inode), GG(iz,:,:,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

   call owfftdestroy()
   deallocate(tfldr,tflds)

   do i=1,nodes
      !$OMP parallel do private(ixw,iz,ish)
      do iz=1,img%z%n
         do ish=1,img%sh%n
            do ixw=1,kxm%n
!               deltaI(iz,ish,ixw)=deltaI(iz,ish,ixw)+GGG(iz,ish,ixw,i)
               deltaI(iz,ish,ixw)=deltaI(iz,ish,ixw)+GG(iz,ish,ixw,i)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end do

   deallocate(wr,ws,ddwr,ddws)
   deallocate(dwr,dws,bwr,bws,GG)
end if

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

end subroutine
end module
