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

implicit none

contains

subroutine owwei(adj,taper,recv,source,sref,vel,deltaS,deltaI)

logical						:: adj
complex,allocatable					:: bwr(:,:,:),bws(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,inode,nodes,i,ixw,ish
real						:: taper(:),vel(:,:),sref(:,:,:)
real						:: deltaS(:,:),deltaI(:,:,:)
integer,allocatable					:: indw(:,:)
real, allocatable					:: GGG(:,:,:,:),GG(:,:,:)
complex, allocatable				:: ws(:,:),wr(:,:)
complex, allocatable				:: ddws(:,:),ddwr(:,:)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
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(bwr(img%z%n,kxm%n,rec%w%n));        bwr=0
allocate(bws(img%z%n,kxm%n,rec%w%n));        bws=0
allocate(ddwr(kxm%n,nodes));                 ddwr=0
allocate(ddws(kxm%n,nodes));                 ddws=0
allocate(dwr(img%z%n,kxm%n,rec%w%n));          dwr=0
allocate(dws(img%z%n,kxm%n,rec%w%n));          dws=0
if (adj) then
   allocate(GG(img%z%n,kxm%n,nodes)); GG=0
else
   allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes)); GGG=0
end if
allocate(indw(rec%w%n,2));indw=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

!$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:sou%h%n+padtraces+tapertr,inode)=source(:,iwa)
      wr(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=recv(:,iwa)
      dws(:,:,iwa)=0.
      dwr(:,:,iwa)=0.
      ddwr(:,inode)=0
      ddws(:,inode)=0
   
      do iz=nwbott+1,img%z%n
   
         call owpspi(  iwa, inode, 1,  ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
         ws(:,inode)=ws(:,inode)*taper
   
         call owpspi(  iwa, inode,-1,  wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
         wr(:,inode)=wr(:,inode)*taper

         !$OMP ORDERED
         indw(iwa,2)=inode
         bws(iz,:,iwa)=ws(:,indw(iwa,2))
         bwr(iz,:,iwa)=wr(:,indw(iwa,2))
         !$OMP END ORDERED

         call owimage(.true., .false., iwa, bws(iz,:,iwa)  , dwr(iz,:,iwa), deltaI(iz,:,:))
         call owimage(.true., .true.,  iwa, dws(iz,:,iwa), bwr(iz,:,iwa)  , deltaI(iz,:,:))
   
      end do
   
      ws(:,inode)=0; wr(:,inode)=0

      do iz=img%z%n,nwbott+1,-1
         ws(:,inode)=ws(:,inode)+dws(iz,:,iwa)
         ws(:,inode)=ws(:,inode)*taper
         call owpspi(  iwa, inode, -1, ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
         ddws(:,inode)=0
         call owphsscr(iwa, 1,ddws(:,inode),bws(iz-1,:,iwa),vel(iz-1,:))
         GG(iz-1,:,inode)=GG(iz-1,:,inode)+real(conjg(ddws(:,inode))*ws(:,inode))
   
         wr(:,inode)=wr(:,inode)+dwr(iz,:,iwa)
         wr(:,inode)=wr(:,inode)*taper
         call owpspi(  iwa, inode, +1, wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
         ddwr(:,inode)=0
         call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz-1,:,iwa),vel(iz-1,:))
         GG(iz-1,:,inode)=GG(iz-1,:,inode)+real(conjg(ddwr(:,inode))*wr(:,inode))
   	
      end do

      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, &
   					" FROM NODE: ",inode
      6 format(A13,i3,A8,i3,A14,i3,A12,i3)

   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)+GG(iz,ixw,i)
         end do
      end do
      !$OMP END PARALLEL DO
   end do

   deallocate(ddwr,ddws,wr,ws)
   deallocate(dwr,dws,GGG,indw)

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(tapertr+padtraces+1:sou%h%n+padtraces+tapertr,inode)=source(:,iwa)
      wr(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=recv(:,iwa)
      dws(:,:,iwa)=0
      dwr(:,:,iwa)=0
      ddwr(:,inode)=0
      ddws(:,inode)=0
 
      do iz=nwbott+1,img%z%n

         call owpspi(  iwa, inode, 1,  ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
         ws(:,inode)=ws(:,inode)*taper
         call owpspi(  iwa, inode,-1,  wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
         wr(:,inode)=wr(:,inode)*taper
         !$OMP ORDERED
         indw(iwa,2)=inode
         bws(iz,:,iwa)=ws(:,indw(iwa,2))
         bwr(iz,:,iwa)=wr(:,indw(iwa,2))
         !$OMP END ORDERED

      end do

      ws(:,inode)=0; wr(:,inode)=0

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

         call owphsscr(iwa, 1,ddws(:,inode),bws(iz,:,iwa),vel(iz-1,:))
         ws(:,inode)=ws(:,inode)+ddws(:,inode)*deltaS(iz-1,:)
         call owpspi(  iwa, inode, 1, ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
         ws(:,inode)=ws(:,inode)*taper
         !$OMP ORDERED
         indw(iwa,2)=inode
         dws(iz,:,iwa)=ws(:,indw(iwa,2))
         !$OMP END ORDERED

         call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz,:,iwa),vel(iz-1,:))
         wr(:,inode)=wr(:,inode)+ddwr(:,inode)*deltaS(iz-1,:)
         call owpspi(  iwa, inode,-1, wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "sou")
         wr(:,inode)=wr(:,inode)*taper
         !$OMP ORDERED
         indw(iwa,2)=inode
         dwr(iz,:,iwa)=wr(:,indw(iwa,2))
         !$OMP END ORDERED

      end do

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

         call owimage(.false., .false., iwa, dws(iz,:,iwa), bwr(iz,:,iwa), GGG(iz,:,:,inode))
         call owimage(.false., .false., iwa, bws(iz,:,iwa), dwr(iz,:,iwa), GGG(iz,:,:,inode))

      end do

      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%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)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end do

   deallocate(wr,ws,ddwr,ddws)
   deallocate(dwr,dws,GGG,indw)

end if

end subroutine

!----------------------------------------------------------

subroutine owwei_f(taper,recv,source,sref,vel,deltaS,deltaI,bwr,bws)

complex,optional					:: bwr(:,:,:),bws(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,inode,nodes,i,ixw,ish
real						:: taper(:),vel(:,:),sref(:,:,:)
real						:: deltaS(:,:),deltaI(:,:,:)
integer,allocatable					:: indw(:,:)
real, allocatable					:: GGG(:,:,:,:)
complex, allocatable				:: ws(:,:),wr(:,:)
complex, allocatable				:: ddws(:,:),ddwr(:,:)
complex, allocatable				:: dws(:,:),dwr(:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
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(dwr(kxm%n,nodes));                  dwr=0
allocate(dws(kxm%n,nodes));                  dws=0
allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes)); GGG=0
allocate(indw(rec%w%n,2));indw=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)

!$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:sou%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

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

       call owpspi(  iwa, inode, 1,  ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
       ws(:,inode)=ws(:,inode)*taper

       call owpspi(  iwa, inode,-1,  wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
       wr(:,inode)=wr(:,inode)*taper
       !$OMP ORDERED
          indw(iwa,2)=inode
          bws(iz,:,iwa)=ws(:,indw(iwa,2))
          bwr(iz,:,iwa)=wr(:,indw(iwa,2))
       !$OMP END ORDERED

    end do

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

       ddws(:,inode)=0
       call owphsscr(iwa, 1,ddws(:,inode),bws(iz,:,iwa),vel(iz-1,:))
       dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaS(iz-1,:)

       ddwr(:,inode)=0
       call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz,:,iwa),vel(iz-1,:))
       dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaS(iz-1,:)

       call owpspi(  iwa, inode, 1, dws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
       call owpspi(  iwa, inode,-1, dwr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "sou")

       dws(:,inode)=dws(:,inode)*taper
       dwr(:,inode)=dwr(:,inode)*taper
 
       call owimage(.false., .false., iwa, dws(:,inode), bwr(iz,:,iwa), GGG(iz,:,:,inode))
       call owimage(.false., .false., iwa, bws(iz,:,iwa), dwr(:,inode), GGG(iz,:,:,inode))

    end do

      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%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)
         end do
      end do
   end do
   !$OMP END PARALLEL DO
end do
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
deallocate(wr,ws,ddwr,ddws)
deallocate(dwr,dws,GGG,indw)

end subroutine
!=========================================================================================
subroutine owwei_b(taper,recv,source,sref,vel,deltaS,deltaI)

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

!$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:sou%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.
    dws_z(:,inode)=0.
    dwr_z(:,inode)=0.
    ddwr(:,inode)=0
    ddws(:,inode)=0

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

       call owpspi(  iwa, inode, 1,  ws(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
       ws(:,inode)=ws(:,inode)*taper

       call owpspi(  iwa, inode,-1,  wr(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
       wr(:,inode)=wr(:,inode)*taper
       !$OMP ORDERED
       indw(iwa,2)=inode
       bws(iz,:,iwa)=ws(:,indw(iwa,2))
       bwr(iz,:,iwa)=wr(:,indw(iwa,2))
       !$OMP END ORDERED
       call owimage(.true., .false., iwa, bws(iz,:,iwa)  , dwr(iz,:,inode), deltaI(iz,:,:))
       call owimage(.true., .true.,  iwa, dws(iz,:,inode), bwr(iz,:,iwa)  , 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
       call owpspi(  iwa, inode, -1, dws_z(:,inode), sref(:,1,iz-1), vel(iz-1,:), tflds(:,inode), "sou")
       ddws(:,inode)=0
       call owphsscr(iwa, 1,ddws(:,inode),bws(iz-1,:,iwa),vel(iz-1,:))
       GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real(conjg(ddws(:,inode))*dws_z(:,inode))

       dwr_z(:,inode)=dwr_z(:,inode)+dwr(iz,:,inode)
       dwr_z(:,inode)=dwr_z(:,inode)*taper
       call owpspi(  iwa, inode, +1, dwr_z(:,inode), sref(:,1,iz-1), vel(iz-1,:), tfldr(:,inode), "rec")
       ddwr(:,inode)=0
       call owphsscr(iwa,-1,ddwr(:,inode),bwr(iz-1,:,iwa),vel(iz-1,:))
       GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real(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%x%n, &
					" FROM NODE: ",inode
6 format(A13,i3,A8,i3,A14,i3,A12,i3)

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,ixw,i)
      end do
   end do
   !$OMP END PARALLEL DO
end do

deallocate(ddwr,ddws,dwr_z,dws_z)
deallocate(dwr,dws,GGG,indw)

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
