module OW_weicg_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
use OW_scttr_areal
use OW_grabvel2_areal
use OW_exec_time
use lloyd_vel
use adj_mod
!use OW_cstret

implicit none

integer, private                             :: idip,nodes
real, allocatable, private                   :: taper2(:), vel2(:,:), velin(:)

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,iz,inode,i,ixw,ish
real						:: taper(:),vel(kxm%n,img%z%n),a
complex						:: image(kxm%n,nsecdim,nz,rec%w%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads

wr=0
ws=0

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,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)
      call owimage_cmplx(.false., iwa, ws(:,inode),  wr(:,inode), image(:,:,1,iwa), vel(:,iz-1))
      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_cmplx(.false., iwa, ws(:,inode),  wr(:,inode), image(:,:,iz-nzmin+1,iwa), vel(:,iz-1))
      end do
      if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
				" FROM NODE: ",inode
   end do
!$OMP END PARALLEL DO
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
call owfftdestroy()
end subroutine
!======================================================================
subroutine owdso(adj,wgcig)
logical						:: adj
real						:: wgcig(img%xm%n,nsecdim,img%z%n)
integer						:: iz,ix,ish
real						:: h(nsecdim)

   do ish=1,nsecdim
      h(ish)=(ish-1)*dsecdim+osecdim
   end do
   h = (abs(h / maxval(h)))!**2

!if (adj) then
!   wgcig=0.
   do ish=1,nsecdim;do ix=1,img%xm%n
      wgcig(ix,ish,:)=wgcig(ix,ish,:)*h(ish)!*ofcig(ix,ish,:)
   end do;end do
!else
!   ofcig=0.
!   do ish=1,nsecdim;do ix=1,img%xm%n
!      ofcig(ix,ish,:)=ofcig(ix,ish,:)+h(ish)*wgcig(ix,ish,:)
!   end do;end do
!end if
end subroutine
!======================================================================
subroutine owwemva_init()
!integer      :: idip_in
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,*) "****************"
!idip=idip_in
end subroutine
!======================================================================
   subroutine owwemva_lop_init ( taper_in,vel_in)
   real                       :: taper_in(:),vel_in(:)
   integer                    :: i,j
      if (allocated(taper2)) deallocate(taper2)
      if (allocated(vel2)) deallocate(vel2)
      if (allocated(velin)) deallocate(velin)
      allocate (taper2(kxm%n), velin(img%xm%n*img%z%n))
      taper2=taper_in
      if (not(areal)) then
         velin=vel_in
         allocate(vel2(kxm%n,img%z%n))
      else
         allocate(vel2(img%xm%n,img%z%n))
         do i=1,img%z%n
            do j=1,img%xm%n
               vel2(j,i)=vel_in((i-1)*img%xm%n+j)
            end do
         end do
      end if

   end subroutine
!======================================================================
   function owwemva_lop ( adj, add, modl, data) result(stat)
    integer               :: stat 
    logical,intent(in)    :: adj,add 
    complex,dimension(:)     :: data 
    complex,dimension(:)  :: modl 
    stat=1
    if (adj) modl=0.
    if (not(adj)) data=0.
    call owwemva_lop2(adj,add,modl,data )
    stat=0
  end function 
!======================================================================
   function owwemva_lopw ( adj, add, modl, data) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    complex,dimension(:)  :: modl,data 
    stat=1
    if (adj) modl=0.
    if (not(adj)) data=0.
    call owwemva_lop2w(adj,add,modl,data )
    stat=0
  end function 
!======================================================================
subroutine owwemva_lop2(adj,add,deltaS,deltaI)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: vel(kxm%n,img%z%n)
complex						:: deltaS(img%xm%n,img%z%n,sou%w%n),deltaSS(kxm%n,img%z%n,sou%w%n)
complex						:: deltaI(img%xm%n,nsecdim,img%z%n)
complex                                                     :: deltaII(kxm%n,nsecdim,img%z%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish

wr=0
ws=0
ddwr=0
ddws=0
bwr=0
bws=0

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,nodes)
tfldr=0
call owfftwr(tfldr,nodes)


call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      if (adj) then
         if (allocated(dwr)) deallocate(dwr); allocate(dwr(kxm%n,img%z%n,nodes))
         if (allocated(dws)) deallocate(dws); allocate(dws(kxm%n,img%z%n,nodes))
         dwr=0;dws=0
         deltaSS=0.;deltaII=0
         deltaII=deltaI(ni:ni+kxm%n-1,:,:)
         !$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)
!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Back project slowness perturbation
! 
! 
            do iz=nzmax,max(2,nzmin),-1
               call owimage_cmplx(.true.,  iwa, dws(:,iz,inode),  bwr(:,iz,inode), deltaII(:,:,iz), vel2(:,iz), .true. )
               dws_z(:,inode)=dws_z(:,inode)+dws(:,iz,inode)
               dws_z(:,inode)=dws_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               ddws(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               else
                  ddws(:,inode)=0
               end if
               if (renorm) then
                  dws_z(:,inode)=dws_z(:,inode)+bws(:,iz,inode)
               end if
               deltaSS(:,iz-1,iwa)=deltaSS(:,iz-1,iwa)+conjg(ddws(:,inode))*dws_z(:,inode)

               call owimage_cmplx(.true.,  iwa, bws(:,iz,inode),  dwr(:,iz,inode), deltaII(:,:,iz), vel2(:,iz), .false.)
               dwr_z(:,inode)=dwr_z(:,inode)+dwr(:,iz,inode)
               dwr_z(:,inode)=dwr_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, +1,  dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ddwr(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               else
                  ddwr(:,inode)=0
               end if
               if (renorm) then
                  dwr_z(:,inode)=dwr_z(:,inode)+bwr(:,iz,inode)
               end if
               deltaSS(:,iz-1,iwa)=deltaSS(:,iz-1,iwa)+conjg(ddwr(:,inode))*dwr_z(:,inode)
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
   !$OMP END PARALLEL DO
         deallocate(dwr,dws)
         deltaS(ni:ni+kxm%n-1,:,:)=deltaS(ni:ni+kxm%n-1,:,:)+deltaSS
     else
         deltaII=0;deltaSS=0.
         allocate(dwr(1,kxm%n,nodes));               dwr=0
         allocate(dws(1,kxm%n,nodes));               dws=0
         deltaSS=deltaS(ni:ni+kxm%n-1,:,:)
         !$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
            bws(:,:,inode)=0.
            bwr(:,:,inode)=0.
            bws(:,1,inode)=ws(:,inode)
            bwr(:,1,inode)=wr(:,inode)
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
            do iz=max(2,nzmin),nzmax
               ddws(:,inode)=0
               if (renorm) then
                  bws(:,iz,inode)=bws(:,iz,inode)+dws(1,:,inode)
               end if
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               end if
               dws(1,:,inode)=dws(1,:,inode)+ddws(:,inode)*deltaSS(:,iz-1,iwa)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, 1,  dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  else
                     dws(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               dws(1,:,inode)=dws(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, dws(1,:,inode), bwr(:,iz,inode), deltaII(:,:,iz),vel2(:,iz-1))

               ddwr(:,inode)=0
               if (renorm) then
                  bwr(:,iz,inode)=bwr(:,iz,inode)+dwr(1,:,inode)
               end if
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               end if
               dwr(1,:,inode)=dwr(1,:,inode)+ddwr(:,inode)*deltaSS(:,iz-1,iwa)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  else
                     dwr(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               dwr(1,:,inode)=dwr(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, bws(:,iz,inode), dwr(1,:,inode), deltaII(:,:,iz),vel2(:,iz-1))
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
         !$OMP END PARALLEL DO
         deallocate(dwr,dws)
         deltaI(ni:ni+kxm%n-1,:,:)=deltaI(ni:ni+kxm%n-1,:,:)+deltaII
      end if
end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,2(i4))
end subroutine
!======================================================================
subroutine owwemva_lop2w(adj,add,deltaS,deltaI)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: vel(kxm%n,img%z%n)
complex						:: deltaS(img%xm%n,img%z%n,sou%w%n),deltaSS(kxm%n,img%z%n,sou%w%n)
complex						:: deltaI(img%xm%n,nsecdim,img%z%n,sou%w%n)
complex                                                     :: deltaII(kxm%n,nsecdim,img%z%n,sou%w%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish

wr=0
ws=0
ddwr=0
ddws=0
bwr=0
bws=0

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,nodes)
tfldr=0
call owfftwr(tfldr,nodes)


call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      if (adj) then
         if (allocated(dwr)) deallocate(dwr); allocate(dwr(kxm%n,img%z%n,nodes))
         if (allocated(dws)) deallocate(dws); allocate(dws(kxm%n,img%z%n,nodes))
         dwr=0;dws=0
         deltaSS=0.;deltaII=0
         deltaII=deltaI(ni:ni+kxm%n-1,:,:,:)
         !$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)
!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Back project slowness perturbation
! 
! 
            do iz=img%z%n,max(2,nzmin),-1
               call owimage_cmplx(.true.,  iwa, dws(:,iz,inode),  bwr(:,iz,inode), deltaII(:,:,iz,iwa), vel2(:,iz), .true. )
               dws_z(:,inode)=dws_z(:,inode)+dws(:,iz,inode)
               dws_z(:,inode)=dws_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               ddws(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               else
                  ddws(:,inode)=0
               end if
               deltaSS(:,iz-1,iwa)=deltaSS(:,iz-1,iwa)+conjg(ddws(:,inode))*dws_z(:,inode)

               call owimage_cmplx(.true.,  iwa, bws(:,iz,inode),  dwr(:,iz,inode), deltaII(:,:,iz,iwa), vel2(:,iz), .false.)
               dwr_z(:,inode)=dwr_z(:,inode)+dwr(:,iz,inode)
               dwr_z(:,inode)=dwr_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, +1,  dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ddwr(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               else
                  ddwr(:,inode)=0
               end if
               deltaSS(:,iz-1,iwa)=deltaSS(:,iz-1,iwa)+conjg(ddwr(:,inode))*dwr_z(:,inode)
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
   !$OMP END PARALLEL DO
   deallocate(dwr,dws)
         deltaS(ni:ni+kxm%n-1,:,:)=deltaS(ni:ni+kxm%n-1,:,:)+deltaSS
     else
         deltaII=0;deltaSS=0.
         allocate(dwr(1,kxm%n,nodes));               dwr=0
         allocate(dws(1,kxm%n,nodes));               dws=0
         deltaSS=deltaS(ni:ni+kxm%n-1,:,:)
         !$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
            bws(:,:,inode)=0.
            bwr(:,:,inode)=0.
            bws(:,1,inode)=ws(:,inode)
            bwr(:,1,inode)=wr(:,inode)
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
            do iz=max(2,nzmin),nzmax
               ddws(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               end if
               dws(1,:,inode)=dws(1,:,inode)+ddws(:,inode)*deltaSS(:,iz-1,iwa)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, 1,  dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  else
                     dws(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               dws(1,:,inode)=dws(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, dws(1,:,inode), bwr(:,iz,inode), deltaII(:,:,iz,iwa),vel2(:,iz-1))

               ddwr(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               end if
               dwr(1,:,inode)=dwr(1,:,inode)+ddwr(:,inode)*deltaSS(:,iz-1,iwa)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  else
                     dwr(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               dwr(1,:,inode)=dwr(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, bws(:,iz,inode), dwr(1,:,inode), deltaII(:,:,iz,iwa),vel2(:,iz-1))
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
         !$OMP END PARALLEL DO
         deallocate(dwr,dws)
         deltaI(ni:ni+kxm%n-1,:,:,:)=deltaI(ni:ni+kxm%n-1,:,:,:)+deltaII
      end if
end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,2(i4))
end subroutine
!======================================================================
subroutine owwemva_forward(deltaS,deltaI)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: deltaS(img%xm%n,img%z%n),deltaSS(kxm%n,img%z%n)
real						:: deltaI(img%xm%n,nsecdim,img%z%n)
real                                                        :: deltaII(kxm%n,nsecdim,img%z%n,nodes)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex						:: dws(kxm%n,nodes),dwr(kxm%n,nodes)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish, a, b

!--------------------------------------------------
! Setting FFTW plans
tflds=0;call owfftws(tflds,nodes)
tfldr=0;call owfftwr(tfldr,nodes)

call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      deltaII=0;deltaSS=0.
      deltaSS=deltaS(ni:ni+kxm%n-1,:)
      a=-999 ; b=999
      do i=1,kxm%n
         if (maxval(deltaSS(i,:)).gt.a) a=maxval(deltaSS(i,:))
         if (minval(deltaSS(i,:)).lt.b) b=minval(deltaSS(i,:))
      end do
      write(0,*) a,b
   !$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.
       bws(:,1,inode)=ws(:,inode)
       bwr(:,1,inode)=wr(:,inode)
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
       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), vel2(:,iz-1), tflds(:,inode), "sou")
                call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ws(:,inode)=ws(:,inode)*taper2
          wr(:,inode)=wr(:,inode)*taper2
          bws(:,iz,inode)=ws(:,inode)
          bwr(:,iz,inode)=wr(:,inode)
       end do
      do iz=max(2,nzmin),nzmax
          ddws(:,inode)=0
          bws(:,iz,inode)=bws(:,iz,inode)+dws(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
          end if
          dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaS(:,iz-1)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          dws(:,inode)=dws(:,inode)*taper2
          call owimage(.false., iwa, dws(:,inode), bwr(:,iz,inode), deltaII(:,:,iz,inode),vel2(:,iz))

          ddwr(:,inode)=0
          bwr(:,iz,inode)=bwr(:,iz,inode)+dwr(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
          end if
          dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaS(:,iz-1)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, -1,  dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          dwr(:,inode)=dwr(:,inode)*taper2
          call owimage(.false., iwa, bws(:,iz,inode), dwr(:,inode), deltaII(:,:,iz,inode),vel2(:,iz-1))

       end do
       if (mod(iwa,50)==1) write(3,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=max(1,nzmin),nzmax
         do ish=1,nsecdim
            do ixw=ni,ni+kxm%n-1
               deltaI(ixw,ish,iz)=deltaI(ixw,ish,iz)+deltaII(ixw-ni+1,ish,iz,i)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end do
end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!======================================================================
subroutine owwemva_forward1(deltaS,deltaI)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: deltaS(img%xm%n,img%z%n),deltaSS(kxm%n,img%z%n)
complex						:: deltaI(img%xm%n,nsecdim,img%z%n,sou%w%n)
complex                                                     :: deltaII(kxm%n,nsecdim,img%z%n,sou%w%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex						:: dws(kxm%n,nodes),dwr(kxm%n,nodes)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish, a, b

!--------------------------------------------------
! Setting FFTW plans
tflds=0;call owfftws(tflds,nodes)
tfldr=0;call owfftwr(tfldr,nodes)

call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      deltaII=0;deltaSS=0.
      deltaSS=deltaS(ni:ni+kxm%n-1,:)
   !$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.
       bws(:,1,inode)=ws(:,inode)
       bwr(:,1,inode)=wr(:,inode)
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
       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), vel2(:,iz-1), tflds(:,inode), "sou")
                call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ws(:,inode)=ws(:,inode)*taper2
          wr(:,inode)=wr(:,inode)*taper2
          bws(:,iz,inode)=ws(:,inode)
          bwr(:,iz,inode)=wr(:,inode)
       end do
      do iz=max(2,nzmin),nzmax
          ddws(:,inode)=0
          bws(:,iz,inode)=bws(:,iz,inode)+dws(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
          end if
          dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaS(:,iz-1)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          dws(:,inode)=dws(:,inode)*taper2
          call owimage_cmplx(.false., iwa, dws(:,inode), bwr(:,iz,inode), deltaII(:,:,iz,iwa),vel2(:,iz))

          ddwr(:,inode)=0
          bwr(:,iz,inode)=bwr(:,iz,inode)+dwr(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
          end if
          dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaS(:,iz-1)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, -1,  dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          dwr(:,inode)=dwr(:,inode)*taper2
          call owimage_cmplx(.false., iwa, bws(:,iz,inode), dwr(:,inode), deltaII(:,:,iz,iwa),vel2(:,iz-1))

       end do
       if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
					" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
   deltaI(ni:ni+kxm%n-1,:,:,:)=deltaI(ni:ni+kxm%n-1,:,:,:)+deltaII
end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!======================================================================
subroutine owwolfe(deltaS,celtaS,n)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k,n
real,optional					:: deltaS(img%xm%n,img%z%n,n)
real 						:: deltaSS(kxm%n,img%z%n,n)
complex,optional					:: celtaS(img%xm%n,img%z%n,n)
complex						:: celtaSS(kxm%n,img%z%n,n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex						:: dws(kxm%n,nodes),dwr(kxm%n,nodes)
complex						:: bws(kxm%n,img%z%n,sou%w%n),bwr(kxm%n,img%z%n,sou%w%n)
complex						:: zws(kxm%n,img%z%n,sou%w%n),zwr(kxm%n,img%z%n,sou%w%n)
complex						:: tws(kxm%n,img%z%n,nodes),twr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish, a, b

!--------------------------------------------------
! Setting FFTW plans
tflds=0;call owfftws(tflds,nodes)
tfldr=0;call owfftwr(tfldr,nodes)

call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      deltaSS=0.
      if (n.eq.1) then
         deltaSS(:,:,1)=deltaS(ni:ni+kxm%n-1,:,1)
      else
         celtaSS(:,:,:)=celtaS(ni:ni+kxm%n-1,:,:)
      end if
   !$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(:,:,iwa)=0.
       bwr(:,:,iwa)=0.
       bws(:,1,iwa)=ws(:,inode)
       bwr(:,1,iwa)=wr(:,inode)
       zws(:,:,iwa)=0.
       zwr(:,:,iwa)=0.
       tws(:,:,inode)=0.
       twr(:,:,inode)=0.
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
       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), vel2(:,iz-1), tflds(:,inode), "sou")
                call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ws(:,inode)=ws(:,inode)*taper2
          wr(:,inode)=wr(:,inode)*taper2
          bws(:,iz,iwa)=ws(:,inode)
          bwr(:,iz,iwa)=wr(:,inode)
          tws(:,iz,inode)=ws(:,inode)
          twr(:,iz,inode)=wr(:,inode)
       end do
      do iz=max(2,nzmin),nzmax
          ddws(:,inode)=0
          tws(:,iz,inode)=tws(:,iz,inode)+dws(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa, 1,ddws(:,inode),tws(:,iz,inode),vel2(:,iz-1))
          end if
          if (n.eq.1) then
             dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaSS(:,iz-1,1)
          else
             dws(:,inode)=dws(:,inode)+ddws(:,inode)*celtaSS(:,iz-1,iwa)
          end if
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          dws(:,inode)=dws(:,inode)*taper2
          zws(:,iz,iwa)=dws(:,inode)

          ddwr(:,inode)=0
          twr(:,iz,inode)=twr(:,iz,inode)+dwr(:,inode)
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa,-1,ddwr(:,inode),twr(:,iz,inode),vel2(:,iz-1))
          end if
          if (n.eq.1) then
             dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*celtaSS(:,iz-1,1)
          else
             dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaSS(:,iz-1,iwa)
          end if
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, -1,  dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          dwr(:,inode)=dwr(:,inode)*taper2
          zwr(:,iz,iwa)=dwr(:,inode)
       end do
       if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
					" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(bwr(:,:,iwa),"bkwr","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(bws(:,:,iwa),"bkws","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(zwr(:,:,iwa),"dkwr","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(zws(:,:,iwa),"dkws","rite",ishot)
    end do
 end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!======================================================================
 subroutine owwemva_close()
 if (allocated(vel2)) deallocate(vel2)
 if (allocated(taper2)) deallocate(taper2)
 if (allocated(velin)) deallocate(velin)
 end subroutine
!!======================================================================
subroutine ow_wvfld_io_init(tag)
character*4    :: tag

    call to_history("n1",img%xm%n,tag); call to_history("d1",img%xm%d,tag); call to_history("o1",img%xm%o,tag)  
    call to_history("n2",nz,tag) ; call to_history("d2",img%z%d,tag) ; call to_history("o2",zmin,tag)
    call to_history("n3",1,tag); call to_history("d3",1,tag); call to_history("o3",0,tag)  
    call to_history("n4",1,tag); call to_history("d4",1,tag); call to_history("o4",1,tag)
    call to_history("n5",rec%w%n,tag) ; call to_history("d5",rec%w%d,tag) ; call to_history("o5",rec%w%o,tag)
    call to_history("n6",rec%xd%n,tag); call to_history("d6",1,tag); call to_history("o6",1,tag)
    call to_history("esize",8,tag)
end subroutine
!!======================================================================
subroutine ow_solver_io_init(tag)
character*4    :: tag

    call to_history("n1",img%xm%n,tag); call to_history("d1",img%xm%d,tag); call to_history("o1",img%xm%o,tag)  
    call to_history("n2",nz,tag) ; call to_history("d2",img%z%d,tag) ; call to_history("o2",zmin,tag)
    call to_history("n3",rec%w%n,tag); call to_history("d3",rec%w%d/(2*pi),tag); call to_history("o3",rec%w%o/(2*pi),tag)
    call to_history("n4",nitercg1,tag); call to_history("d4",1,tag); call to_history("o4",1,tag)  
    call to_history("n5",1,tag); call to_history("d5",1,tag); call to_history("o5",0,tag)
    call to_history("esize",8,tag)

end subroutine
!!======================================================================
subroutine ow_err_io_init()

    call to_history("n1",niter,"err"); call to_history("d1",1,"err"); call to_history("o1",1,"err")  
    call to_history("n2",1,"err"); call to_history("d2",1,"err"); call to_history("o2",0,"err")
    call to_history("esize",4,"err")

end subroutine
!!======================================================================
subroutine ow_wvfld_io(CCC,BBB,tag,op,iii)
character*4       :: tag,op
complex,optional  :: CCC(:,:),BBB(:,:,:)
complex           :: DDD(img%xm%n,nz)
integer           :: n0,nnn
integer, optional :: iii

nnn=idip
if(present(iii)) nnn=iii
if (op.eq."rite") then
   DDD=0
   n0 = (((nnn-1)*rec%xd%d+rec%xd%o+rec%h%o-(padtraces+tapertr)*rec%h%d)-img%xm%o)/img%xm%d+1.5
   DDD(n0:kxm%n+n0-1,:)=CCC
   call srite(tag,DDD,8*img%xm%n*nz)
else if (op.eq."reed") then
   call sreed(tag,BBB,8*img%xm%n*nz*rec%w%n)
end if

end subroutine
!!======================================================================
subroutine ow_wvfld_io2(CCC,BBB,tag,op,iii)
character*4       :: tag,op
complex,optional  :: CCC(:,:),BBB(:)
complex           :: DDD(img%xm%n,nz)
integer           :: n0,nnn
integer, optional :: iii

nnn=idip
if(present(iii)) nnn=iii
if (op.eq."rite") then
   DDD=0
   n0 = (((nnn-1)*rec%xd%d+rec%xd%o+rec%h%o-(padtraces+tapertr)*rec%h%d)-img%xm%o)/img%xm%d+1.5
   DDD(n0:kxm%n+n0-1,:)=CCC
   call srite(tag,DDD,8*img%xm%n*nz)
else if (op.eq."reed") then
   call sreed(tag,BBB,8*img%xm%n*nz*rec%w%n)
end if

end subroutine
!!======================================================================
subroutine ow_wvfld_io1(CCC,tag,op,ishot)
character*4    :: tag,op
complex        :: CCC(:,:),DDD(img%xm%n,nz)
integer        :: n0,ishot


if (op.eq."rite") then
   DDD=0
   n0 = (((ishot-1)*rec%xd%d+rec%xd%o+rec%h%o-(padtraces+tapertr)*rec%h%d)-img%xm%o)/img%xm%d+1.5
   DDD(n0:kxm%n+n0-1,:)=CCC
   call srite(tag,DDD,8*img%xm%n*nz)
else if (op.eq."reed") then
   call sreed(tag,CCC,8*img%xm%n*nz*rec%w%n)
end if

end subroutine

!======================================================================
subroutine owdso_renorm(adj,wgcig)!,ofcig)
logical						:: adj
real						:: wgcig(img%xm%n,nsecdim,img%z%n,sou%w%n)
integer						:: iz,ix,ish
real						:: h(nsecdim)

   do ish=1,nsecdim
      h(ish)=(ish-1)*dsecdim+osecdim
   end do
   h = (abs(h / maxval(h)))**expdso

   do iz=1,sou%w%n;do ish=1,nsecdim;do ix=1,img%xm%n
      wgcig(ix,ish,:,iz)=wgcig(ix,ish,:,iz)*h(ish)
   end do;end do;end do
end subroutine
!======================================================================
   function owwemva_renorm_lop ( adj, add, modl, data) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    complex,dimension(:)  :: modl 
    real,dimension(:)  :: data 
    stat=1
    if (adj) modl=0.
    if(not(adj)) data=0.
    call owwemva_renorm(adj,add,modl,data )
    stat=0
  end function 
!======================================================================
subroutine owwei_renorm(taper,recv,source,vel,image,sref)

real, optional					:: sref(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,inode,i,ixw,ish
real						:: taper(:),vel(kxm%n,img%z%n),a
real						:: image(:,:,:)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,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 (extrap(1:3).eq."ffd") then
            if (rec%w%all(iwa).ne.0) then
               call owffd( iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
               call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,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), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
               call owpspi( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
            else
               call owpspi_init(iz)
               call owpspi( iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
               call owpspi( iwa, inode,-1,  wr(:,inode), vel2(:,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), image(:,:,iz-nzmin+1), vel2(:,iz-1))
      end do
      if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
				" FROM NODE: ",inode
   end do
!$OMP END PARALLEL DO
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
call owfftdestroy()
end subroutine
!======================================================================
subroutine owwemva_renorm(adj,add,deltaSS,deltaII)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k,mm
real						:: vel(kxm%n,img%z%n)
complex						:: deltaS(kxm%n,img%z%n,sou%w%n),deltaSS(img%xm%n,img%z%n,sou%w%n)
real						:: deltaI(kxm%n,nsecdim,img%z%n),deltaII(img%xm%n,nsecdim,img%z%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex						:: dws(kxm%n,nodes),dwr(kxm%n,nodes)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish,a

!--------------------------------------------------
! Setting FFTW plans
tflds=0;call owfftws(tflds,nodes)
tfldr=0;call owfftwr(tfldr,nodes)
call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
do ishot=1,rec%xd%n
 start=0.;finish=0.
 call owexec_time(start)
 if (areal) then
   source=0;
   call sreed("sou",source,8*sou%h%n*sou%w%n)
   ni=1
 else
   call owgrab_velsht(velin,vel2,ni,ishot)
   if (lloyd) then
      sref=0.
      call lloyd_vel1(1/vel2,sref)
   end if
 end if
 recv=0
 call sreed("rec",recv,8*rec%h%n*rec%w%n)
 if (adj) then
   deltaS=0;deltaI=0
   do j=ni,ni+kxm%n-1;do ish=1,img%sh%n;do iz=1,nz
      deltaI(j-ni+1,ish,iz)=deltaI(j-ni+1,ish,iz)+deltaII(j,ish,iz)
   end do;end do;end do
   !$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(:,inode)=0;wr(:,inode)=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)
       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)*taper2
       bwr(:,1,inode)=wr(:,inode)*taper2

!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
      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), vel2(:,iz-1), tflds(:,inode), "sou")
               call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
            end if
         else if (extrap(1:4).eq."pspi") then
            if (lloyd) then
               call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
               call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
            else
               call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
               call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
            end if
         end if
         ws(:,inode)=ws(:,inode)*taper2
         wr(:,inode)=wr(:,inode)*taper2
         bws(:,iz,inode)=ws(:,inode)
         bwr(:,iz,inode)=wr(:,inode)
       end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Back project image perturbation
! 
! 
       do iz=nzmax,max(2,nzmin),-1
          dws(:,inode)=0
          call owimage(.true.,  iwa, dws(:,inode),  bwr(:,iz,inode), deltaI(:,:,iz), vel2(:,iz-1), .true. )
          dws_z(:,inode)=dws_z(:,inode)+dws(:,inode)
          dws_z(:,inode)=dws_z(:,inode)*taper2
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ge.tiny(1.)) then
                call owffd( iwa, inode, -1,  dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if(lloyd) then
                call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          ddws(:,inode)=0
          if (rec%w%all(iwa).ge.tiny(1.)) then
             call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
          end if
          deltaS(:,iz-1,iwa)=deltaS(:,iz-1,iwa)+conjg(ddws(:,inode))*dws_z(:,inode)

          dwr(:,inode)=0.
          call owimage(.true.,  iwa, bws(:,iz,inode),  dwr(:,inode), deltaI(:,:,iz), vel2(:,iz-1), .false.)
          dwr_z(:,inode)=dwr_z(:,inode)+dwr(:,inode)
          dwr_z(:,inode)=dwr_z(:,inode)*taper2
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ge.tiny(1.)) then
                call owffd( iwa, inode, +1,  dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if(lloyd) then
                call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ddwr(:,inode)=0
          if (rec%w%all(iwa).ge.tiny(1.)) then
             call owscctrx(iwa, -1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
          end if
          deltaS(:,iz-1,iwa)=deltaS(:,iz-1,iwa)+conjg(ddwr(:,inode))*dwr_z(:,inode)
          
       end do
       if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
					" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
    do j=ni,ni+kxm%n-1;do k=1,nz;do iwa=1,rec%w%n
       deltaSS(j,k,iwa)=deltaSS(j,k,iwa)+deltaS(j-ni+1,k,iwa)
    end do;end do;end do
 else
    deltaI=0;deltaS=0
    do j=ni,ni+kxm%n-1;do iz=1,nz;do iwa=1,rec%w%n
       deltaS(j-ni+1,iz,iwa)=deltaS(j-ni+1,iz,iwa)+deltaSS(j,iz,iwa)
    end do;end do ;end do
    !$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(:,inode)=0;wr(:,inode)=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)
       dws(:,inode)=0
       dwr(:,inode)=0
       bws(:,:,inode)=0.
       bwr(:,:,inode)=0.
       bws(:,1,inode)=ws(:,inode)*taper2
       bwr(:,1,inode)=wr(:,inode)*taper2
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
       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), vel2(:,iz-1), tflds(:,inode), "sou")
                call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ws(:,inode)=ws(:,inode)*taper2
          wr(:,inode)=wr(:,inode)*taper2
          bws(:,iz,inode)=ws(:,inode)
          bwr(:,iz,inode)=wr(:,inode)
       end do
      do iz=max(2,nzmin),nzmax
          ddws(:,inode)=0
          if (rec%w%all(iwa).ne.0.) then
              call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
          end if
          dws(:,inode)=dws(:,inode)+ddws(:,inode)*deltaS(:,iz-1,iwa)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          dws(:,inode)=dws(:,inode)*taper2
          call owimage(.false., iwa, dws(:,inode), bwr(:,iz,inode), deltaI(:,:,iz),vel2(:,iz))

          ddwr(:,inode)=0
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa, -1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
          end if
          dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*deltaS(:,iz-1,iwa)
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, -1,  dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          dwr(:,inode)=dwr(:,inode)*taper2
          call owimage(.false., iwa, bws(:,iz,inode), dwr(:,inode), deltaI(:,:,iz),vel2(:,iz-1))
       end do
       if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
					" FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
   do j=ni,ni+kxm%n-1;do iwa=1,img%sh%n;do iz=1,nz
      deltaII(j,iwa,iz)=deltaII(j,iwa,iz)+deltaI(j-ni+1,iwa,iz)
   end do;end do;end do

 end if
end do
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!======================================================================
   function owdeltaS_renorm_lop ( adj, add, modl, data) result(stat)
    integer               :: stat 
    logical,intent(in)    :: adj,add 
    complex,dimension(:)  :: modl 
    complex,dimension(:)  :: data 
    stat=1
    if (adj) modl=0.
    if(not(adj)) data=0.
    call owrenorm_exe(adj,add,modl,data )
    stat=0
  end function 
!============================================================
subroutine owrenorm_exe(adj,add,mm,dd)
logical             :: adj,add
complex             :: bwr(img%xm%n*nz*sou%w%n),bws(img%xm%n*nz*sou%w%n),br,bs,dr,ds
complex             :: dwr(img%xm%n*nz*sou%w%n),dws(img%xm%n*nz*sou%w%n)
complex             :: mm(img%xm%n,nz,sou%w%n)
complex             :: dd(img%xm%n,nz,sou%w%n)
integer             :: i,j,k,l,inode
real                :: as,ar
integer, external	:: omp_get_thread_num,omp_get_num_threads

do i=1,rec%xd%n
   bwr=0;bws=0;dwr=0;dws=0
   call ow_wvfld_io2(BBB=bwr,tag="bkwr",op="reed",iii=i);ar=0.01*sum(abs(bwr)/sqrt(float(img%xm%n*nz*sou%w%n)))
   call ow_wvfld_io2(BBB=bws,tag="bkws",op="reed",iii=i);as=0.01*sum(abs(bws)/sqrt(float(img%xm%n*nz*sou%w%n)))
   call ow_wvfld_io2(BBB=dwr,tag="dkwr",op="reed",iii=i)
   call ow_wvfld_io2(BBB=dws,tag="dkws",op="reed",iii=i)
   if (adj) then
      !$OMP parallel do private(inode,j,k,l,br,bs,dr,ds) ordered schedule(dynamic,1) num_threads(nodes)
      do j=1,rec%w%n
         do k=1,nz
            do l=1,img%xm%n
               br=bwr((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               bs=bws((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               dr=dwr((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               ds=dws((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               mm(l,k,j)=mm(l,k,j) + (1+dr*conjg(br)/(conjg(br)*br+ar)+ds*conjg(bs)/(conjg(bs)*bs+as))*dd(l,k,j)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   else
      !$OMP parallel do private(inode,j,k,l,br,bs,dr,ds) ordered schedule(dynamic,1) num_threads(nodes)
      do j=1,rec%w%n
         do k=1,nz
            do l=1,img%xm%n
               br=bwr((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               bs=bws((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               dr=dwr((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               ds=dws((j-1)*nz*img%xm%n+(k-1)*img%xm%n+l)
               dd(l,k,j)=dd(l,k,j) + (1+dr*conjg(br)/(conjg(br)*br+ar)+ds*conjg(bs)/(conjg(bs)*bs+as))*mm(l,k,j)
            end do
         end do
      end do
      !$OMP END PARALLEL DO
   end if
end do   
call sseek("bkws",0,0)
call sseek("bkwr",0,0)
call sseek("dkws",0,0)
call sseek("dkwr",0,0)
end subroutine
!======================================================================
subroutine owwfld_renorm(dd,ddw)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: vel(kxm%n,img%z%n)
complex,optional					:: dd(img%xm%n,img%z%n),ddw(img%xm%n,img%z%n,sou%w%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex						:: dws(kxm%n,nodes),dwr(kxm%n,nodes)
complex						:: bws(kxm%n,img%z%n,sou%w%n),bwr(kxm%n,img%z%n,sou%w%n)
complex						:: zws(kxm%n,img%z%n,sou%w%n),zwr(kxm%n,img%z%n,sou%w%n)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,nodes)
tfldr=0
call owfftwr(tfldr,nodes)
call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
if (not(areal)) then
   source1=0;
   source=0;
   call sreed("sou",source1,8*sou%w%n)
   source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
end if
!---------------------------------
!
! Shot loop
!
do ishot=1,rec%xd%n
  start=0.;finish=0. 
  call owexec_time(start)
  if (areal) then
     source=0;
     call sreed("sou",source,8*sou%h%n*sou%w%n)
     ni=1
  else
     call owgrab_velsht(velin,vel2,ni,ishot)
     if (lloyd) then
        sref=0.
        call lloyd_vel1(1/vel2,sref)
     end if
  end if
  recv=0
  call sreed("rec",recv,8*rec%h%n*rec%w%n)

  !$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(:,:,iwa)=0.
     bwr(:,:,iwa)=0.
     bws(:,1,iwa)=ws(:,inode)
     bwr(:,1,iwa)=wr(:,inode)
     zws(:,:,iwa)=0.
     zwr(:,:,iwa)=0.
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
       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), vel2(:,iz-1), tflds(:,inode), "sou")
                call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          ws(:,inode)=ws(:,inode)*taper2
          wr(:,inode)=wr(:,inode)*taper2
          bws(:,iz,iwa)=ws(:,inode)
          bwr(:,iz,iwa)=wr(:,inode)
       end do
       do iz=max(2,nzmin),nzmax
          ddws(:,inode)=0
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,iwa),vel2(:,iz-1))
          end if
          if (present(ddw)) then
             dws(:,inode)=dws(:,inode)+ddws(:,inode)*ddw(:,iz-1,iwa)
          else
             dws(:,inode)=dws(:,inode)+ddws(:,inode)*dd(:,iz-1)
          end if
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
             end if
          end if
          dws(:,inode)=dws(:,inode)*taper2
          zws(:,iz,iwa)=dws(:,inode)

          ddwr(:,inode)=0
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,iwa),vel2(:,iz-1))
          end if
          if (present(ddw)) then
             dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*ddw(:,iz-1,iwa)
          else
             dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*dd(:,iz-1)
          end if
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, -1,  dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (lloyd) then
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode,-1, dwr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
             end if
          end if
          dwr(:,inode)=dwr(:,inode)*taper2
          zwr(:,iz,iwa)=dwr(:,inode)
       end do

   end do   
  !$OMP END PARALLEL DO
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(bwr(:,:,iwa),"bkwr","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(bws(:,:,iwa),"bkws","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(zwr(:,:,iwa),"dkwr","rite",ishot)
    end do
    do iwa = 1,rec%w%n
       call ow_wvfld_io1(zws(:,:,iwa),"dkws","rite",ishot)
    end do
 end do
 call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!!===============================================================
  subroutine owdotst_image()

    integer				:: i,j,k,iwa
    complex				:: data1a(img%xm%n,img%z%n,rec%w%n),data1b(img%xm%n,img%z%n,rec%w%n)
    complex				:: data2a(img%xm%n,img%z%n,rec%w%n),data2b(img%xm%n,img%z%n,rec%w%n)
    real            		          :: vel(img%xm%n,img%z%n)
    double precision                              :: tempm(img%xm%n,nsecdim,img%z%n,rec%w%n),tempd(img%xm%n,img%z%n,rec%w%n)
    complex				:: modl1(img%xm%n,nsecdim,img%z%n,rec%w%n),modl2(img%xm%n,nsecdim,img%z%n,rec%w%n)
    double precision, dimension (2)	 	:: dot

    tempd=0;data1a=0;data1b=0;data2a=0
    data2b=0;tempm=0;modl1=0;modl2=0;vel=0

    call random_number(tempm)
    modl1=cmplx(cosd(3.14159*tempm),sind(3.14159*tempm))
    call random_number(tempd)
    data2a=cmplx(tempd,sin(3.14159*tempd))
    call random_number(tempd)
    data2b=cmplx(tempd,cos(3.14159*tempd))
    vel=tempd(:,:,1)
    
    modl2=0.;data1a=0.;data1b=0.;dot=0.

    data2a=data2a/100;data2b=data2b/100
    modl1=modl1/10000;

!    !$OMP parallel do private(iwa,k) schedule(dynamic,1) num_threads(12)
    do iwa=1,rec%w%n
       do k=1,img%z%n
          call owimage_cmplx(.false., iwa, data2a(:,k,iwa), data2b(:,k,iwa), modl2(:,:,k,iwa), vel(:,k), .false.)
          call owimage_cmplx(.true.,  iwa, data1a(:,k,iwa), data2b(:,k,iwa), modl1(:,:,k,iwa), vel(:,k),  .true.)
          call owimage_cmplx(.true.,  iwa, data2a(:,k,iwa), data1b(:,k,iwa), modl1(:,:,k,iwa), vel(:,k), .false.)
       end do
    end do
!    !$OMP end parallel do
!    !$OMP parallel do private(iwa,i,j,k) schedule(dynamic,1) num_threads(12)
    do iwa=1,rec%w%n
       do k=1,img%z%n
          do i=1,img%xm%n
             dot(1)=dot(1)+(conjg(data1a(i,k,iwa))*data2a(i,k,iwa)+conjg(data1b(i,k,iwa))*data2b(i,k,iwa))/2
             do j=1,nsecdim
                dot(2)=dot(2)+conjg(modl1(i,j,k,iwa))*modl2(i,j,k,iwa)
             end do
          end do
       end do
    end do
!    !$OMP end parallel do
    write (0,*) dot
    write (0,*) abs(dot)

end subroutine
!!===============================================================
  subroutine owdotst_image1()

    integer				:: i,j,k,iwa,inode
    complex				:: data1a(img%xm%n,img%z%n,rec%w%n),data1b(img%xm%n,img%z%n,rec%w%n)
    complex				:: data2a(img%xm%n,img%z%n,rec%w%n),data2b(img%xm%n,img%z%n,rec%w%n)
    double precision        		          :: tempd(img%xm%n,img%z%n,rec%w%n)
    double precision                              :: tempm(img%xm%n,nsecdim,img%z%n)
    real					:: vel(img%xm%n,img%z%n)
    complex				:: modl1(img%xm%n,nsecdim,img%z%n),modl2(img%xm%n,nsecdim,img%z%n)
    double precision, dimension (2)	 	:: dot
    integer, external			:: omp_get_thread_num,omp_get_num_threads   

    tempd=0;data1a=0;data1b=0;data2a=0
    data2b=0;tempm=0;modl1=0;modl2=0;vel=0

    call random_number(tempm)
    modl1=cosd(3.14159*tempm)
    call random_number(tempd)
    data2a=cmplx(tempd,sin(3.14159*tempd))
    call random_number(tempd)
    data2b=cmplx(tempd,cos(3.14159*tempd))
    vel=tempd(:,:,1)
    
    modl2=0.;data1a=0.;data1b=0.;dot=0.

    data2a=data2a/100;data2b=data2b/100
    modl1=modl1/10000;

!    !$OMP parallel do private(iwa,k,inode) schedule(dynamic,1) num_threads(1)
    do iwa=1,rec%w%n
!       inode=omp_get_thread_num()+1
       do k=1,img%z%n
          call owimage_cmplx(.false., iwa, data2a(:,k,iwa), data2b(:,k,iwa), modl2(:,:,k), vel(:,k), .false.)
          call owimage_cmplx(.true.,  iwa, data1a(:,k,iwa), data2b(:,k,iwa), modl1(:,:,k), vel(:,k),  .true.)
          call owimage_cmplx(.true.,  iwa, data2a(:,k,iwa), data1b(:,k,iwa), modl1(:,:,k), vel(:,k), .false.)
       end do
    end do
!    !$OMP end parallel do
!    !$OMP parallel do private(iwa,i,j,k) schedule(dynamic,1) num_threads(1)
     do iwa=1,rec%w%n
        do k=1,img%z%n
           do i=1,img%xm%n
              dot(1)=dot(1)+(conjg(data1a(i,k,iwa))*data2a(i,k,iwa)+conjg(data1b(i,k,iwa))*data2b(i,k,iwa))/2
           end do
        end do
     end do
!    !$OMP end parallel do
!    !$OMP parallel do private(inode,i,j,k) schedule(dynamic,1) num_threads(1)
!     do inode=1,12
        do k=1,img%z%n
           do i=1,img%xm%n
              do j=1,nsecdim
                 dot(2)=dot(2)+modl1(i,j,k)*modl2(i,j,k)
              end do
           end do
        end do
!     end do
!    !$OMP end parallel do
    write (0,*) dot

end subroutine
!!
  function igrad1_lop ( adj, add, xx, yy) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    complex,dimension(:)  :: xx,yy 
!    if (adj) xx=0.
!    if(not(adj)) yy=0.
    call igrad1_lop2(adj,add,xx,yy )
    stat=0
  end function 
!!
  subroutine igrad1_lop2(adj,add,xx,yy)
    logical,intent(in)   :: adj,add 
    complex              :: xx(img%xm%n,img%z%n,rec%w%n),yy(img%xm%n,img%z%n,rec%w%n)
    integer              :: i,j,k
    do i = 1, rec%w%n-1
       do j = 1,img%z%n 
          do k = 1,img%xm%n 
             if ( adj) then
                xx(k,j,i+1) = xx(k,j,i+1) + yy(k,j,i)
                xx(k,j,i  ) = xx(k,j,i  ) - yy(k,j,i)
             else
                yy(k,j,i) = yy(k,j,i) + xx(k,j,i+1) - xx(k,j,i)
             end if
          end do
       end do
    end do
  end subroutine 
!!===============================================================
!======================================================================
   function owwemva_ren ( adj, add, modl, data) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    complex,dimension(:)  :: modl,data 
    stat=1
    if (adj) modl=0.
    if (not(adj)) data=0.
    call owwemva_exe(adj,add,modl,data )
    stat=0
  end function 
!======================================================================
subroutine owwemva_exe(adj,add,deltaS,deltaI)

real      					:: sref(nref,2,img%z%n)
logical						:: adj,add
complex						:: recv(rec%h%n,sou%w%n),source(rec%h%n,sou%w%n)
complex                                                     :: source1(sou%w%n)
integer						:: iwa,iz,inode,i,ixw,ish,ni,ishot,j,k
real						:: vel(kxm%n,img%z%n)
complex						:: deltaS(img%xm%n,img%z%n),deltaSS(kxm%n,img%z%n)
complex						:: deltaI(img%xm%n,nsecdim,img%z%n)
complex                                                     :: deltaII(kxm%n,nsecdim,img%z%n)
complex						:: ws(kxm%n,nodes),wr(kxm%n,nodes)
complex						:: ddws(kxm%n,nodes),ddwr(kxm%n,nodes)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:)
complex						:: bws(kxm%n,img%z%n,nodes),bwr(kxm%n,img%z%n,nodes)
complex						:: tflds(kxm%n,nodes),tfldr(kxm%n,nodes)
complex						:: dws_z(kxm%n,nodes),dwr_z(kxm%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
real                                                        :: start, finish

wr=0
ws=0
ddwr=0
ddws=0
bwr=0
bws=0

!--------------------------------------------------
! Setting FFTW plans
tflds=0
call owfftws(tflds,nodes)
tfldr=0
call owfftwr(tfldr,nodes)


call sseek("rec",0,0)
call sseek("sou",0,0)
!---------------------------------
   if (not(areal)) then
      source1=0;
      source=0;
      call sreed("sou",source1,8*sou%w%n)
      source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
   end if
!---------------------------------
!
! Shot loop
!
   do ishot=1,rec%xd%n
      start=0.;finish=0. 
      call owexec_time(start)
      if (areal) then
         source=0;
         call sreed("sou",source,8*sou%h%n*sou%w%n)
         ni=1
      else
         call owgrab_velsht(velin,vel2,ni,ishot)
         if (lloyd) then
            sref=0.
            call lloyd_vel1(1/vel2,sref)
         end if
      end if
      recv=0
      call sreed("rec",recv,8*rec%h%n*rec%w%n)
      if (adj) then
         if (allocated(dwr)) deallocate(dwr); allocate(dwr(kxm%n,img%z%n,nodes))
         if (allocated(dws)) deallocate(dws); allocate(dws(kxm%n,img%z%n,nodes))
         dwr=0;dws=0
         deltaSS=0.;deltaII=0
         deltaII=deltaI(ni:ni+kxm%n-1,:,:)
         !$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)
!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Back project slowness perturbation
! 
! 
            do iz=nzmax,max(2,nzmin),-1
               call owimage_cmplx(.true.,  iwa, dws(:,iz,inode),  bwr(:,iz,inode), deltaII(:,:,iz), vel2(:,iz), .true. )
               dws_z(:,inode)=dws_z(:,inode)+dws(:,iz,inode)
               dws_z(:,inode)=dws_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, -1, dws_z(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               ddws(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               else
                  ddws(:,inode)=0
               end if
               if (renorm) then
!                  dws_z(:,inode)=dws_z(:,inode)+bws(:,iz,inode)
                  bws(:,iz,inode)=bws(:,iz,inode)+dws_z(:,inode)
               end if
               deltaSS(:,iz-1)=deltaSS(:,iz-1)+conjg(ddws(:,inode))*dws_z(:,inode)

               call owimage_cmplx(.true.,  iwa, bws(:,iz,inode),  dwr(:,iz,inode), deltaII(:,:,iz), vel2(:,iz), .false.)
               dwr_z(:,inode)=dwr_z(:,inode)+dwr(:,iz,inode)
               dwr_z(:,inode)=dwr_z(:,inode)*taper2
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, +1,  dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if(lloyd) then
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz))
                  else
                     call owpspi(  iwa, inode, +1, dwr_z(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ddwr(:,inode)=0
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               else
                  ddwr(:,inode)=0
               end if
               if (renorm) then
!                  dwr_z(:,inode)=dwr_z(:,inode)+bwr(:,iz,inode)
                  bwr(:,iz,inode)=bwr(:,iz,inode)+dwr_z(:,inode)
               end if
               deltaSS(:,iz-1)=deltaSS(:,iz-1)+conjg(ddwr(:,inode))*dwr_z(:,inode)
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
         !$OMP END PARALLEL DO
         deltaS(ni:ni+kxm%n-1,:)=deltaS(ni:ni+kxm%n-1,:)+deltaSS
     else
         deltaII=0;deltaSS=0.
         allocate(dwr(1,kxm%n,nodes));               dwr=0
         allocate(dws(1,kxm%n,nodes));               dws=0
         deltaSS=deltaS(ni:ni+kxm%n-1,:)
         !$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
            bws(:,:,inode)=0.
            bwr(:,:,inode)=0.
            bws(:,1,inode)=ws(:,inode)
            bwr(:,1,inode)=wr(:,inode)
!!!!!!!!!!!!!!!!!!!!!!!
!
! Background wavefields
! bws and bwr
!
            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), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owffd( iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1,  ws(:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                     call owpspi(  iwa, inode,-1,  wr(:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               ws(:,inode)=ws(:,inode)*taper2
               wr(:,inode)=wr(:,inode)*taper2
               bws(:,iz,inode)=ws(:,inode)
               bwr(:,iz,inode)=wr(:,inode)
            end do
            do iz=max(2,nzmin),nzmax
               ddws(:,inode)=0
               if (renorm) then
                  bws(:,iz,inode)=bws(:,iz,inode)+dws(1,:,inode)
               end if
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa, 1,ddws(:,inode),bws(:,iz,inode),vel2(:,iz-1))
               end if
               dws(1,:,inode)=dws(1,:,inode)+ddws(:,inode)*deltaSS(:,iz-1)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, 1,  dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  else
                     dws(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode, 1, dws(1,:,inode), vel2(:,iz-1), tflds(:,inode), "sou")
                  end if
               end if
               dws(1,:,inode)=dws(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, dws(1,:,inode), bwr(:,iz,inode), deltaII(:,:,iz),vel2(:,iz-1))

               ddwr(:,inode)=0
               if (renorm) then
                  bwr(:,iz,inode)=bwr(:,iz,inode)+dwr(1,:,inode)
               end if
               if (rec%w%all(iwa).ne.0.) then
                  call owscctrx(iwa,-1,ddwr(:,inode),bwr(:,iz,inode),vel2(:,iz-1))
               end if
               dwr(1,:,inode)=dwr(1,:,inode)+ddwr(:,inode)*deltaSS(:,iz-1)
               if (extrap(1:3).eq."ffd") then
                  if (rec%w%all(iwa).ne.0) then
                     call owffd( iwa, inode, -1,  dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  else
                     dwr(1,:,inode)=0
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (lloyd) then
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec", sref(:,:,iz-1))
                  else
                     call owpspi(  iwa, inode,-1, dwr(1,:,inode), vel2(:,iz-1), tfldr(:,inode), "rec")
                  end if
               end if
               dwr(1,:,inode)=dwr(1,:,inode)*taper2
               call owimage_cmplx(.false., iwa, bws(:,iz,inode), dwr(1,:,inode), deltaII(:,:,iz),vel2(:,iz-1))
            end do
            if (mod(iwa,50)==1) write(3,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                           " FROM NODE: ",inode,ishot
         end do
         !$OMP END PARALLEL DO
         deltaI(ni:ni+kxm%n-1,:,:)=deltaI(ni:ni+kxm%n-1,:,:)+deltaII
      end if
end do
deallocate(dwr,dws)
call owfftdestroy()
6 format(A13,i3,A8,i3,A14,i3,A12,2(i4))
end subroutine
!======================================================================
end module
