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_phsscr_areal
use OW_cstret

implicit none

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

real, optional					:: sref(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,inode,nodes,i,ixw,ish
real						:: taper(:),vel(:,:),a
real						:: image(:,:,:)
real, allocatable					:: GGG(:,:,:,:)
!complex						:: image(:,:,:)
!complex, allocatable				:: GGG(:,:,:,:)
complex, allocatable				:: ws(:,:),wr(:,:),wws(:,:),wwr(:,:)
real, allocatable					:: Gst(:,:,:),Grt(:,:,:)
complex, allocatable				:: Gsw(:,:,:),Grw(:,:,:)
complex, allocatable				:: GGs(:,:,:),GGr(:,:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
complex,allocatable					:: tfldst(:,:),tfldrt(:,:)
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,*) "****************"

if (exist_file("ws")) then
   allocate(GGr(img%z%n,kxm%n,rec%w%n)); GGr=0
   allocate(GGs(img%z%n,kxm%n,rec%w%n)); GGs=0
end if
allocate(wr(kxm%n,nodes));                  wr=0
allocate(ws(kxm%n,nodes));                  ws=0
allocate(wwr(kxm%n,nodes));                  wwr=0
allocate(wws(kxm%n,nodes));                  wws=0
if(not(icwind)) then
!   allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes))
   allocate(GGG(nz,img%sh%n,kxm%n,nodes))
   GGG=0
end if

if (icwind) then
   allocate(Grw(img%z%n,img%xm%n,rec%w%n)); Grw=0
   allocate(Gsw(img%z%n,img%xm%n,rec%w%n)); Gsw=0
end if

!--------------------------------------------------
! 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)

!$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)

       do iz=max(2,nzmin),nzmax

          if (present(sref)) then
!             call owpspi1( iz, iwa, inode, 1,  ws(:,inode), wws(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
!             call owpspi1( iz, iwa, inode,-1,  wr(:,inode), wwr(:,inode), vel(iz-1,:), tfldr(:,inode), "rec", sref(:,:,iz-1))
!write(0,*) "iz  ",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 owpspi1( iz, iwa, inode, 1,  ws(:,inode), wws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
!             call owpspi1( iz, iwa, inode,-1,  wr(:,inode), wwr(:,inode), vel(iz-1,:), tfldr(:,inode), "rec")
             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
!          wws(:,inode)=wws(:,inode)*taper
!          wwr(:,inode)=wwr(:,inode)*taper


          if (icwind) then
             Grw(iz,:,iwa)=wr(:,inode)
             Gsw(iz,:,iwa)=ws(:,inode)
          else
             call owimage(.false., iwa, ws(:,inode),  wr(:,inode), GGG(iz-nzmin+1,:,:,inode))
          end if

          if (exist_file("ws")) then
!             !$OMP ORDERED
             GGR(iz,:,iwa)=wr(:,inode)
             GGS(iz,:,iwa)=ws(:,inode)
!             !$OMP END ORDERED
          end if

       end do

          if (not(icwind)) then
             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

if (exist_file("ws")) then
   call srite("ws",GGs, 8*img%z%n*rec%w%n*img%xm%n)
   call srite("wr",GGr, 8*img%z%n*rec%w%n*img%xm%n)
!   call srite("ws",GGs(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,:), &
!             8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
!   call srite("wr",GGr(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,:), &
!             8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
end if
deallocate(wr,ws)

!if (icwind) then
!   allocate(Grt(nodes,kxm%n,2*n_wind-1)); Grt=0 
!   allocate(Gst(nodes,kxm%n,2*n_wind-1)); Gst=0 
!   allocate(tfldrt(n_sampt,nodes));            tfldrt=0
!   allocate(tfldst(n_sampt,nodes));            tfldst=0
!
!   call owfftts(tfldst,nodes)
!   call owffttr(tfldrt,nodes)
!
!!$OMP parallel do private(inode,&
!                            iwa,&
!                             iz,&
!                            ixw) ordered schedule(dynamic,1) num_threads(nodes)
!
!   do iz=nwbott+1,img%z%n
!      inode=omp_get_thread_num()+1
!      do ixw=1,kxm%n
!         tfldst(:,inode)=0.;tfldrt(:,inode)=0.
!         tfldst(1:rec%w%n,inode)=Gsw(iz,ixw,:)
!         tfldrt(1:rec%w%n,inode)=Grw(iz,ixw,:)
!
!         call owfftexest(tfldst(:,inode),inode)
!         call owfftexert(tfldrt(:,inode),inode)
!
!         Gst(inode,ixw,n_wind:2*n_wind-1)=real(tfldst(1:n_wind,inode))
!         Grt(inode,ixw,n_wind:2*n_wind-1)=real(tfldrt(1:n_wind,inode))
!
!         do iwa=1,n_wind-1
!            Grt(inode,ixw,iwa)=real(tfldrt(n_sampt-iwa+1,inode))
!            Gst(inode,ixw,iwa)=real(tfldst(n_sampt-iwa+1,inode))
!         end do
!      end do
!
!      do iwa=1,2*n_wind-1
!         call owimage2(Gst(inode,:,iwa),Grt(inode,:,iwa),image(iz,:,:))
!      end do
!   end do
!
!!$OMP END PARALLEL DO
!
!   call owfftdestroy()
!   deallocate(tfldr,tflds)
!   deallocate(Grt,Gst,tfldrt,tfldst)
!
!else
   do i=1,nodes
      !$OMP parallel do private(ixw,iz,ish)
      do iz=1,nz
         do ish=1,img%sh%n
            do ixw=1,kxm%n
               image(iz,ish,ixw)=image(iz,ish,ixw)+GGG(iz,ish,ixw,i)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end do
   deallocate(GGG)
!end if
if (exist_file("ws")) then
   deallocate(GGr,GGs)
end if
6 format(A13,i3,A8,i3,A14,i3,A12,i3)

end subroutine

!======================================================================
subroutine owdso(adj,wgcig,ofcig)

logical						:: adj
real						:: wgcig(:,:,:),ofcig(:,:,:),a,b
!real						:: a,b
!complex						:: wgcig(:,:,:),ofcig(:,:,:)
integer						:: iz,ix,ish
real, allocatable					:: h(:),w(:,:,:),temp(:)

allocate(h(img%sh%n),w(img%z%n,img%sh%n,img%xm%n),temp(img%z%n))

do ish=1,img%sh%n
   h(ish)=(ish-1)*img%sh%d+img%sh%o
end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! So para teste. Comentar
! para aplicacao do DSO
!h=1
!ish=int(img%sh%n/2)+1.5
!h(ish)=0.;h(ish-1)=0.5;h(ish+1)=0.5
h = abs(h / maxval(h))
b=0
if (adj) then

   do ix=1,img%xm%n
      do ish=1,img%sh%n
         temp(1)=ofcig(1,ish,ix)
         temp(img%z%n)=ofcig(img%z%n,ish,ix)
         do iz=2,img%z%n-1
            temp(iz)=sum(abs(ofcig(iz-1:iz+1,ish,ix)))/3
         end do
         w(:,ish,ix)=temp
      end do
   end do

   do ix=1,img%xm%n
      do ish=1,img%sh%n
         do iz=1,img%z%n
            wgcig(iz,ish,ix)=wgcig(iz,ish,ix)+h(ish)*ofcig(iz,ish,ix)
            a=abs(w(iz,ish,ix))
            if (a > b) then
               b=a
            end if
         end do
      end do
   end do
!   wgcig=wgcig*ofcig*ofcig/b/b
!   wgcig=wgcig*w*w/b/b
else

   do ix=1,img%xm%n
      do ish=1,img%sh%n
         temp(1)=wgcig(1,ish,ix)
         temp(img%z%n)=wgcig(img%z%n,ish,ix)
         do iz=2,img%z%n-1
            temp(iz)=sum(abs(wgcig(iz-1:iz+1,ish,ix)))/3
         end do
         w(:,ish,ix)=temp
      end do
   end do

   do ix=1,img%xm%n
      do ish=1,img%sh%n
         do iz=1,img%z%n
            ofcig(iz,ish,ix)=ofcig(iz,ish,ix)+h(ish)*wgcig(iz,ish,ix)
            a=abs(w(iz,ish,ix))
            if (a > b) then
               b=a
            end if
         end do
      end do
   end do
!   ofcig=ofcig*wgcig*wgcig/b/b
!   ofcig=ofcig*w*w/b/b
end if
deallocate(h,w,temp)
end subroutine
!======================================================================
subroutine ow_resmig(imgi, rho) 

 complex, dimension(:,:,:) :: imgi
 complex, dimension(:), allocatable :: tri, tro
 real   , dimension(:), allocatable :: str
 real                      :: rho
 integer                   :: ikx, ish, iz, stat
 real                      :: ks, kg
 real                      :: kz1, kz2, kzs, kzg

   allocate(tri(kz%n), tro(kz%n), str(kz%n))
   eps=1.
   stat = base_cstret_init(kz%n, kz%n, kz%o, kz%d, eps)

   do ikx=1,kxm%n
      do ish=1,ksh%n
         ks = sqrt((kxm%all(ikx)-ksh%all(ish))**2)
         kg = sqrt((kxm%all(ikx)+ksh%all(ish))**2)
         do iz=1,kz%n
            kz1 = (4.*kz%all(iz)**2+(kg-ks)**2)/4.
            kz2 = (4.*kz%all(iz)**2+(kg+ks)**2)/4.
            kzs = 1./rho**2 * kz1*kz2 - kz%all(iz)**2*ks**2
            kzg = 1./rho**2 * kz1*kz2 - kz%all(iz)**2*kg**2
            if (kzs>=0. .and. kzg>=0) then
               str(iz) = 0.5/kz%all(iz)*(sqrt(kzs)+sqrt(kzg))
            else
               str(iz) = -2*kz%d
            end if
         end do
         tri = imgi(:,ish,ikx)
         stat = base_cstret_oper(str, tri, tro)
         imgi(:,ish,ikx) = tro
      end do
   end do
deallocate(tri,tro,str)
end subroutine
!======================================================================
subroutine ow_resmig_lin(imgi, rho) 

 complex, dimension(:,:,:) :: imgi
 real                      :: rho,rho2
 integer                   :: ikx, ish, iz, stat
 real                      :: ks, kg,mu,ks2,kg2,mu2
 real                      :: kz1, kz2, kzs, kzg,dkzdrho

   eps=1.

   rho2=rho*rho

   do ikx=1,kxm%n
      do ish=1,ksh%n
         ks = sqrt((kxm%all(ikx)-ksh%all(ish))**2);ks2=ks*ks
         kg = sqrt((kxm%all(ikx)+ksh%all(ish))**2);kg2=kg*kg
         do iz=1,kz%n
            kz1 = (4.*kz%all(iz)**2+(kg-ks)**2)/4.
            kz2 = (4.*kz%all(iz)**2+(kg+ks)**2)/4.
            if (iz.eq.1) then
               mu=huge(1.);mu2=huge(1.)
               kzs = huge(1.)
               kzg = huge(1.)
               dkzdrho=0.
            else
               mu=kz1*kz2/(16*kz%all(iz));mu2=mu*mu
               kzs = 0.5*sqrt(rho2*mu2-ks2)
               kzg = 0.5*sqrt(rho2*mu2-kg2)
               if ((mu2-ks2.ne.0.).or.(mu2-kg2.ne.0.)) then
                  dkzdrho=0.5*mu2*(1/sqrt(mu2-ks2)+1/sqrt(mu2-kg2))
               end if
            end if
!            write(0,*) ikx,ish,iz,dkzdrho
            imgi(iz,ish,ikx)=imgi(iz,ish,ikx)*dkzdrho
         end do
      end do
   end do

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

!!!!!!!!!!!!!!!!!!!!!!
!
! Scattered wavefield
!
       do iz=1,img%z%n
          call owimage(.true.,  iwa, bws(iz,:,inode),  dwr(iz,:,inode), deltaI(iz,:,:), .false.)
          call owimage(.true.,  iwa, dws(iz,:,inode),  bwr(iz,:,inode), deltaI(iz,:,:), .true. )
       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.,iwa, dws(1,:,inode), bwr(iz,:,inode), GGG(iz,:,:,inode))
!          call owimage(.false.,iwa, bws(iz,:,inode), dwr(1,:,inode), GGG(iz,:,:,inode))
          call owimage(.false., iwa, dws(1,:,inode), bwr(iz,:,inode), GG(iz,:,:,inode))
          call owimage(.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
!========================================================================||
!  Datuming
!========================================================================||
!    call owimage0(iwa,ws,wr,GG(iz,:,:))
! if (nwbott.ne.1) then
!    call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),vref,minv,maxv)
!    tfld=0.; tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
!    tfld=0.; tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
! !      Operations in the wavenumber domain
! !      phase shift
!    do ikx = 1,kxm%n
!       sqa = rec%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2
!       if (sqa > 0) then
!          wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
!          ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
!       else
!          wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
!          ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
!       end if
!    end do
!    tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
!    tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
!    wr=wr*taper
!  ws=ws*taper
! end if
!========================================================================
!========================================================================
!
! Derivative in k-domain
!
! c=(/1.,0.5,0.375,0.3125,0.2734375/)
!
!    tflds(:,inode)=source(img%z%n-1,:,iwa)
!    call owfftexes(tflds(:,inode),inode)
!    aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!
!    tfldr(:,inode)=recv(img%z%n-1,:,iwa)
!    call owfftexer(tfldr(:,inode),inode)
!    awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
!
!    do j=0,4
!       dwss(:,inode)=kxm%all**(2*j)*aws(:,inode)
!       dwrr(:,inode)=kxm%all**(2*j)*awr(:,inode)
!
!       tflds(:,inode)=dwss(:,inode)
!       call owifftexes(tflds(:,inode),inode)
!       dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
!
!       tfldr(:,inode)=dwrr(:,inode)
!       call owifftexer(tfldr(:,inode),inode)
!       dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
!
!       ddws(:,inode)=ddws(:,inode)+dwss(:,inode)
!       ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode)
!    end do
!
!    ddws(:,inode)=ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
!    ddwr(:,inode)=ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
!
!
!    do ixw=1,kxm%n
!       do ish=1,img%sh%n
!          ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
!!          if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then
!!             dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ihso,iwa)*deltaI(img%z%n,ish,ixw)
!!             dws(ixw,inode)=dws(ixw,inode)+recv(img%z%n,ixw+ihso,iwa)*deltaI(img%z%n,ish,ixw)
!          if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then
!             dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ish,iwa)*deltaI(img%z%n,ish,ixw)
!             dws(ixw,inode)=dws(ixw,inode)+(recv(img%z%n,ixw+ish,iwa))*deltaI(img%z%n,ish,ixw)
!          end if
!       end do
!    end do
!    GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real(conjg(ddws(:,inode))*dws(:,inode)+conjg(ddwr(:,inode))*dwr(:,inode))
!    GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real((ddws(:,inode))*dws(:,inode)+(ddwr(:,inode))*dwr(:,inode))
!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   Compute derivative in k-domain
!
!       aws(:,inode)=0.
!       tflds(:,inode)=source(iz-1,:,iwa)
!       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")
!
!
!       awr(:,inode)=0.
!       tfldr(:,inode)=recv(iz-1,:,iwa)
!       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")
!
!!!!!!!!!!!!!!!!
!
!  derivative in k
!       call owphsscr(iwa, inode,  1,  ddws(:,inode), ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
!       call owphsscr(iwa, inode, -1,  ddwr(:,inode), wr(:,inode), vel(iz-1,:), tfldr(:,inode), "rec")
!       dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaS(iz,:)
!       dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaS(iz,:)
!       tflds(:,inode)=dws(:,inode)
!       call owfftexes(tflds(:,inode),inode)
!       dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!       tfldr(:,inode)=dwr(:,inode)
!       call owfftexer(tfldr(:,inode),inode)
!       dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Output some wavefields
!
!
!if (exist_file("ws")) then
!!$OMP ORDERED
!    indw(iwa,2)=inode
!    call srite("ws",dws(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,indw(iwa,2)), &
!                8*img%z%n*(img%xm%n-2*(tapertr+padtraces)))
!    call srite("wr",dwr(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,indw(iwa,2)), &
!                8*img%z%n*(img%xm%n-2*(tapertr+padtraces)))
!!$OMP END ORDERED
!end if
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Changes
!
!
!
!
!
!
end module
