module OW_wei_areal

use sep
use OW_type
use OW_image_areal
use OW_fftw_areal
use OW_parms_areal
use OW_pspi_areal
use OW_ffd_areal
use OW_scttr_areal
use OW_exec_time
use OW_mpi

implicit none

integer, private                             :: idip,nodes,ireal

contains
!======================================================================
subroutine owwei_init(idip_in,ireal_in)
integer                                 :: idip_in,ireal_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(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"
idip=idip_in
ireal=ireal_in
end subroutine
!======================================================================
subroutine owwei(taper,vel,sref)

real, optional			:: sref(:,:,:)
integer				:: iwa,iwb,izb,iz,inode,i,ixw,ihx,iyw,ihy,initz,finz,shtidx,i_r
integer                                 :: ierr,status
real				:: taper(:,:),vel(kxm%n,kym%n,img%z%n),ww(nodes),inic,fini
real,allocatable                        :: image(:,:,:,:,:),image0(:,:,:,:,:),GGG(:,:,:,:,:,:)
complex,allocatable			:: ws(:,:,:),wr(:,:,:),tflds(:,:,:)
integer, external			:: omp_get_thread_num,omp_get_num_threads

allocate(image(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns))
allocate(ws(kxm%n,kym%n,wblock%ns),wr(kxm%n,kym%n,wblock%ns),tflds(kxm%n,kym%n,nodes))
!--------------------------------------------------
! Setting FFTW plans
if (wpar%i_mpi > 0) then
   call owexec_time(inic)
   tflds=0;call owfftws(tflds,nodes)
   call owexec_time(fini)
end if 
if (idip == shtblk) then
   if (shtres.eq.0) then
      shtidx = shtblm
   else
      shtidx = shtres
   end if
else
   shtidx = shtblm
end if
999 format(a17,7i4)
do iwb=1, wblock%nb
   if (wpar%i_mpi==0) then
      do i_r=1,shtidx
write(0,999) "owwei: send data ",i_r,idip,shtblk,shtidx,shtblm,(idip-1)*shtblm+i_r
         ws=0;wr=0
write(0,*) "Sent shot ",(idip-1)*shtblm+i_r, " to node ", i_r," freq.block = ",iwb
         call gathers_io(iwb,(idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
         call gathers_io(iwb,(idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"sou","reed")
         call MPI_SEND_(wr,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,i_r,120,MPI_COMM_WORLD,ierr)
         call MPI_SEND_(ws,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,i_r,130,MPI_COMM_WORLD,ierr)
         call MPI_SEND_(i_r,1,MPI_INTEGER,i_r,131,MPI_COMM_WORLD,ierr)
      end do
   else
      allocate(GGG(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns,nodes))
      ws=0;wr=0
      call MPI_RECV_(wr,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,0,120,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV_(ws,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,0,130,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV_(i_r,1,MPI_INTEGER,0,131,MPI_COMM_WORLD,status,ierr)
write(0,*) "RECV shot ",(idip-1)*shtblm+i_r," in node ", wpar%i_mpi," freq.block = ",iwb

      do izb=1,zblock%nb
         GGG=0.;image=0;
         if (izb==1) then
            initz=2
         else
            initz=1
         end if
         !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
         do iwa = 1,wblock%blen(iwb)
            inode=omp_get_thread_num()+1; 
            ww(inode) = wblock%os(iwb) + (iwa-1)*wblock%ds
            if (initz==2) call owimage(.false., ww(inode), ws(:,:,iwa),  wr(:,:,iwa), GGG(:,:,:,:,1,inode))
            do iz=initz,zblock%blen(izb)
               if (extrap(1:3).eq."ffd") then
                  if (ww(inode).ne.0) then
                     call owffd( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode))
                     call owffd( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode))
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (present(sref)) then
                     call owpspi( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode)&
                                , sref(:,:,iz-1+zblock%indx(izb)))
                     call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode)&
                                , sref(:,:,iz-1+zblock%indx(izb)))
                  else
                     call owpspi( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode))
                     call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+zblock%indx(izb)), tflds(:,:,inode))
                  end if
               end if
               ws(:,:,iwa)=ws(:,:,iwa)*taper
               wr(:,:,iwa)=wr(:,:,iwa)*taper
               call owimage(.false., ww(inode), ws(:,:,iwa),  wr(:,:,iwa), GGG(:,:,:,:,iz,inode))
            end do
            if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",(idip-1)*shtblm+i_r, &
   				" FROM NODE: ",inode
         end do
         !$OMP END PARALLEL DO
         !$OMP parallel do private(i,ixw,iyw,iz,ihx,ihy) schedule(dynamic,1) num_threads(nodes)
         do i=1,nodes
            do iz=1,zblock%blen(izb)
               do ihy=1,nsecdimy
                  do iyw=1,kym%n
                     do ihx=1,nsecdimx
                        do ixw=1,kxm%n
                           image(ixw,ihx,iyw,ihy,iz)=image(ixw,ihx,iyw,ihy,iz)+GGG(ixw,ihx,iyw,ihy,iz,i)
                        end do
                     end do
                  end do
               end do
            end do
         end do
         !$OMP END PARALLEL DO
         call image_io("bkgimg","reed",izb,image,.true.)
         call image_io("bkgimg","rite",izb,image)
      end do
      if (allocated(GGG)) deallocate(GGG)
   end if
end do
if (allocated(image0)) deallocate(image0); allocate(image0(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns));image0=0

do izb=1,zblock%nb
   image0=0
   if(wpar%i_mpi > 0) call image_io("bkgimg","reed",izb,image,.false.)
   call MPI_REDUCE_(image(:,:,:,:,1:zblock%blen(izb)),image0(:,:,:,:,1:zblock%blen(izb)),&
            kxm%n*nsecdimx*kym%n*nsecdimy*zblock%blen(izb),MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
   if(wpar%i_mpi == 0) call image_io("bkgimg","rite",izb,image0)
end do


deallocate(image,image0)
deallocate(ws,wr,tflds)
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
if (wpar%i_mpi > 0) call owfftdestroy()
write(0,*) " 12 passei" , wpar%i_mpi,idip
end subroutine
!========================================
subroutine gathers_io(iwb,isht,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwb,isht
character*3          :: tag
character*4          :: op

if(op.eq."reed") then 

  call sreed_window(tag,5,(/rec%hx%n,rec%hy%n,            sou%w%n, rec%xd%n,  nreal/),&
                          (/rec%hx%n,rec%hy%n,   wblock%blen(iwb),        1,      1/),&
                          (/       0,       0, wblock%indx(iwb)-1,   isht-1,ireal-1/),&
                          (/       1,       1,                   1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwb))) 

else
  call srite_window(tag,5,(/rec%hx%n,rec%hy%n,            sou%w%n, rec%xd%n,  nreal/),&
                          (/rec%hx%n,rec%hy%n,   wblock%blen(iwb),        1,      1/),&
                          (/       0,       0, wblock%indx(iwb)-1,   isht-1,ireal-1/),&
                          (/       1,       1,                   1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwb))) 
end if
end subroutine
!========================================
subroutine wavefield_io(iwb,izb,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwb,izb
character*4          :: tag
character*4          :: op

if(op.eq."reed") then 
  call sreed_window(tag,4,(/rec%hx%n,rec%hy%n,          sou%w%n,rec%xd%n,  nreal/),&
                          (/rec%hx%n,rec%hy%n, wblock%blen(iwb),       1,      1/),&
                          (/       0,       0,(iwb-1)*wblock%ns,  idip-1,ireal-1/),&
                          (/       1,       1,                1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwb))) 
else
  call srite_window(tag,4,(/rec%hx%n,rec%hy%n,          sou%w%n,rec%xd%n,  nreal/),&
                          (/rec%hx%n,rec%hy%n, wblock%blen(iwb),       1,      1/),&
                          (/       0,       0,(iwb-1)*wblock%ns,  idip-1,ireal-1/),&
                          (/       1,       1,                1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwb))) 
end if
end subroutine
!========================================
subroutine wavefield_io1(iwb,iwa,izb,GG,tag,op)
complex              :: GG(:,:)
integer              :: iwb,izb,iwa,iw
character*4          :: tag
character*4          :: op

call sseek(tag,0,0)
iw=wblock%indx(iwb)-1+iwa
if(op.eq."reed") then 
!if (wpar%i_mpi == 1) write(0,*) "before reading wvfl ",tag, wblock%indx(iwb)+iwa-2, izb-1,iw
  call sreed_window(tag,4,(/kxm%n,kym%n, img%z%n, sou%w%n/),&
                          (/kxm%n,kym%n,       1,       1/),&
                          (/    0,    0,   izb-1,    iw-1/),&
                          (/    1,    1,       1,       1/),&
                          8,GG) 
!if (wpar%i_mpi == 1) write(0,*) "after reading wvfl ",tag
else
!if (wpar%i_mpi == 1) write(0,*) "before writing wvfl ",tag, wblock%indx(iwb)+iwa-2, izb-1,iw
  call srite_window(tag,4,(/kxm%n,kym%n, img%z%n, sou%w%n/),&
                          (/kxm%n,kym%n,       1,       1/),&
                          (/    0,    0,   izb-1,    iw-1/),&
                          (/    1,    1,       1,       1/),&
                          8,GG) 
!if (wpar%i_mpi == 1) write(0,*) "after writing wvfl ",tag
end if
end subroutine
!========================================
subroutine image_io(tag,op,izb,image,sum)
real                          :: image(:,:,:,:,:)
real,allocatable              :: temp(:,:,:,:,:)
integer                       :: izb
real                          :: tmp,a,b
integer                       :: i,nn1,nn2,nn3,nn4,nn5
character*4                   :: op
character*6                   :: tag
logical,optional              :: sum

nn1=size(image,1)
nn2=size(image,2)
nn3=size(image,3)
nn4=size(image,4)
nn5=size(image,5)

if(op.eq."reed") then
  allocate(temp(nn1,nn2,nn3,nn4,zblock%blen(izb)));temp=0.
  call sreed_window(tag,5,(/nn1,nn2,nn3,nn4,           img%z%n/),&
                          (/nn1,nn2,nn3,nn4,  zblock%blen(izb)/),&
                          (/  0,  0,  0,  0,zblock%indx(izb)-1/),&
                          (/  1,  1,  1,  1,                 1/),&
                          4,temp) 
  if(sum) image(:,:,:,:,1:zblock%blen(izb))=image(:,:,:,:,1:zblock%blen(izb))+temp
  if(not(sum)) image(:,:,:,:,1:zblock%blen(izb))=temp
  deallocate(temp)
else if(op.eq."rite") then
  call srite_window(tag,5,(/nn1,nn2,nn3,nn4,           img%z%n/),&
                          (/nn1,nn2,nn3,nn4,  zblock%blen(izb)/),&
                          (/  0,  0,  0,  0,zblock%indx(izb)-1/),&
                          (/  1,  1,  1,  1,                 1/),&
                          4,image(:,:,:,:,1:zblock%blen(izb)))
end if
end subroutine
!========================================
subroutine image_io1(tag,op,izb,image)
real,optional                 :: image(:,:,:,:)
integer,optional              :: izb
integer                       :: nn1,nn2,nn3,nn4,nn5
character*4                   :: op
character*6                   :: tag

nn1=size(image,1)
nn2=size(image,2)
nn3=size(image,3)
nn4=size(image,4)

if(op.eq."reed") then
  call sreed_window(trim(tag),5,(/nn1,nn2,nn3,nn4,img%z%n/),&
                          (/nn1,nn2,nn3,nn4,      1/),&
                          (/  0,  0,  0,  0,  izb-1/),&
                          (/  1,  1,  1,  1,      1/),&
                          4,image) 
else if(op.eq."rite") then
  call srite_window(trim(tag),5,(/nn1,nn2,nn3,nn4,img%z%n/),&
                          (/nn1,nn2,nn3,nn4,      1/),&
                          (/  0,  0,  0,  0,  izb-1/),&
                          (/  1,  1,  1,  1,      1/),&
                          4,image) 
end if
end subroutine
!========================================
subroutine grad_io(tag,op,image)
real                          :: image(:,:,:)
real, allocatable             :: tmp(:,:,:)
integer                       :: nn1,nn2,nn3
character*4                   :: op
character*5                   :: tag

nn1=size(image,1)
nn2=size(image,2)
nn3=size(image,3)

if(op.eq."reed") then
  allocate(tmp(nn1,nn2,nn3))
  call sreed_window(trim(tag),3,(/nn1,nn2,nn3/),&
                                (/nn1,nn2,nn3/),&
                                (/  0,  0,  0/),&
                                (/  1,  1,  1/),&
                                         4,tmp) 
  image=image+tmp
  deallocate(tmp)
else if(op.eq."rite") then
  call srite_window(trim(tag),3,(/nn1,nn2,nn3/),&
                                (/nn1,nn2,nn3/),&
                                (/  0,  0,  0/),&
                                (/  1,  1,  1/),&
                                       4,image) 
end if
end subroutine
!!======================================================================
subroutine owdso(fun,riteD,pick1)
real,optional                                               :: pick1(:)
real,allocatable					:: wgcig(:,:,:,:)
integer						:: ix,ihx,iy,ihy,izb,finz,iz,k,ierr,i_r,status
real						:: hx, hy, hhmax,hh
double precision                                            :: fun
logical                                                     :: riteD

allocate(wgcig(kxm%n,nsecdimx,kym%n,nsecdimy))
   hx=max(abs(osecdimx),(nsecdimx-1)*dsecdimx+osecdimx)
   hy=max(abs(osecdimy),(nsecdimy-1)*dsecdimy+osecdimy)
   hhmax=sqrt(hx*hx+hy*hy)
   k=(yline_-img%ym%o)/img%ym%d+1.5

   do izb=1,nz
      wgcig=0
      if(wpar%i_mpi==0) then
         call image_io1("bkgimg","reed",izb,wgcig)
         do ihy=1,nsecdimy
            hy=(ihy-1)*dsecdimy+osecdimy
            do iy=1,img%ym%n
               do ihx=1,nsecdimx
                  hx=(ihx-1)*dsecdimx+osecdimx
                  hh=(sqrt(hx*hx+hy*hy)/hhmax)**(2*expdso)
                  do ix=1,img%xm%n
                     if (present(pick1)) then
!                        if ((izb-1)*img%z%d+img%z%o < int(pick1((iy-1)*img%xm%n+ix))) wgcig(ix,ihx,iy,ihy)=0
                        if ((izb-1)*img%z%d+img%z%o < int(pick1((ix-1)*img%ym%n+iy))) wgcig(ix,ihx,iy,ihy)=0
                     end if
                     wgcig(ix,ihx,iy,ihy)=wgcig(ix,ihx,iy,ihy)*hh
                     fun = fun + wgcig(ix,ihx,iy,ihy)*wgcig(ix,ihx,iy,ihy)
                  end do
               end do
            end do
         end do
         if (riteD) call image_io1("Dimage","rite",izb,wgcig)
         do i_r=1,wpar%n_mpi-1
            call MPI_SEND_(wgcig,kxm%n*nsecdimx*kym%n*nsecdimy,MPI_REAL,i_r,7120,MPI_COMM_WORLD,ierr)
         end do
      else
         call MPI_RECV_(wgcig,kxm%n*nsecdimx*kym%n*nsecdimy,MPI_REAL,0,7120,MPI_COMM_WORLD,status,ierr)
         if (riteD) call image_io1("Dimage","rite",izb,wgcig)
      end if   
   end do
deallocate(wgcig)

end subroutine
!======================================================================
subroutine owwemva_init(idip_in)
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(deltaS0,taper,vel,sref)

real, optional					:: sref(:,:,:)
logical						:: adj
integer						:: iwa,iz,inode,i,ixw,ish,n0,iyw,ihx,ihy,i_r
integer                                                     :: initz,finz,izz,izb
real,allocatable					:: deltaS(:,:,:)
real                                                        :: deltaS0(kxm%n,kym%n,img%z%n)
real,allocatable					:: deltaI(:,:,:,:,:)
real                                                        :: taper(:,:),vel(kxm%n,kym%n,img%z%n)
real,allocatable					:: GGG(:,:,:,:)
real                                                        :: ww(nodes)
complex,allocatable					:: ws(:,:,:),wr(:,:,:)
complex,allocatable					:: bws(:,:,:,:),bwr(:,:,:,:)
complex,allocatable					:: ddws(:,:,:),ddwr(:,:,:)
complex,allocatable					:: dws(:,:,:),dwr(:,:,:)
complex,allocatable					:: tflds(:,:,:)
complex,allocatable					:: dws_z(:,:,:),dwr_z(:,:,:)
integer, external					:: omp_get_thread_num,omp_get_num_threads
integer                                                     :: status,ierr,SHTIDX

write(0,*) "passei aa"
allocate (deltaS(kxm%n,kym%n,img%z%n));deltaS=0
allocate (ws(kxm%n,kym%n,wblock%ns),wr(kxm%n,kym%n,wblock%ns));ws=0;wr=0
if (wpar%i_mpi > 0) then
   allocate (deltaI(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns));deltaI=0
!   allocate(GGG(kxm%n,kym%n,zblock%ns,nodes));GGG=0
   allocate(GGG(kxm%n,kym%n,nz,nodes));GGG=0
   allocate (bws(kxm%n,kym%n,nz,wblock%ns),bwr(kxm%n,kym%n,nz,wblock%ns));bws=0;bwr=0
   allocate (ddws(kxm%n,kym%n,nodes),ddwr(kxm%n,kym%n,nodes));ddws=0;ddwr=0
   allocate (dws(kxm%n,kym%n,nodes),dwr(kxm%n,kym%n,nodes),tflds(kxm%n,kym%n,nodes));dws=0;dwr=0;tflds=0
   allocate (dws_z(kxm%n,kym%n,nodes),dwr_z(kxm%n,kym%n,nodes));dws_z=0;dwr_z=0
end if
write(0,*) "passei bb"
!--------------------------------------------------
! Setting FFTW plans
if (wpar%i_mpi > 0) then
   tflds=0;call owfftws(tflds,nodes)
end if 
if (idip == shtblk) then
   if (shtres.eq.0) then
      shtidx = shtblm
   else
      shtidx = shtres
   end if
else
   shtidx = shtblm
end if

do iwb=1, wblock%nb
   if (wpar%i_mpi==0) then
      do i_r=1,shtidx
         ws=0;wr=0
         call gathers_io(iwb,(idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
         call gathers_io(iwb,(idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"sou","reed")
         call MPI_SEND_(wr,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,i_r,5120,MPI_COMM_WORLD,ierr)
         call MPI_SEND_(ws,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,i_r,5130,MPI_COMM_WORLD,ierr)
      end do
   else
      ws=0;wr=0
      call MPI_RECV_(wr,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,0,5120,MPI_COMM_WORLD,status,ierr)
      call MPI_RECV_(ws,kxm%n*kym%n*wblock%blen(iwb),MPI_COMPLEX,0,5130,MPI_COMM_WORLD,status,ierr)
      !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
      do iwa = 1,wblock%blen(iwb)
         inode=omp_get_thread_num()+1; 
         ww(inode) = wblock%os(iwb) + (iwa-1)*wblock%ds
         bwr(:,:,:,iwa)=0
         bws(:,:,:,iwa)=0
         !!!!!!!!!!!!!!!!!!!!!!!!
         !! Background wavefields
         !! bws and bwr
         do iz=2,nz
            if (extrap(1:3).eq."ffd") then
               if (ww(inode).ne.0) then
                  call owffd( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode))
                  call owffd( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode))
               end if
            else if (extrap(1:4).eq."pspi") then
               if (present(sref)) then
                  call owpspi( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode)&
                          , sref(:,:,iz-1))
                  call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode)&
                             , sref(:,:,iz-1))
               else
                  call owpspi( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode))
                  call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1), tflds(:,:,inode))
               end if
            end if
            ws(:,:,iwa)=ws(:,:,iwa)*taper
            wr(:,:,iwa)=wr(:,:,iwa)*taper
            bws(:,:,iz,iwa)=ws(:,:,iwa)
            bwr(:,:,iz,iwa)=wr(:,:,iwa)
         end do
      end do
      !$OMP END PARALLEL DO

write(0,*) "passei 00",wpar%i_mpi
      do iwa = 1,nodes
write(0,*) "$$$$$$$$$",iwa,wpar%i_mpi
         do ixw=1,kxm%n; do iyw=1,kym%n
            dws(ixw,iyw,iwa)=0.
            dwr(ixw,iyw,iwa)=0.
            dws_z(ixw,iyw,iwa)=0.
            dwr_z(ixw,iyw,iwa)=0.
            ddwr(ixw,iyw,iwa)=0.
            ddws(ixw,iyw,iwa)=0.
            GGG(ixw,iyw,:,iwa)=0.
          end do;end do
      end do
write(0,*) "passei 01",wpar%i_mpi
      do izb=zblock%nb,1,-1
write(0,*) "passei 01.1",wpar%i_mpi
         initz=1; if (izb.eq.1) initz=2
write(0,*) "passei 01.2",wpar%i_mpi
         finz=zblock%blen(izb)+zblock%indx(izb)
write(0,*) "passei 01.3",wpar%i_mpi
         deltaI=0;call image_io("Dimage","reed",izb,deltaI,.false.)
write(0,*) "passei 02",wpar%i_mpi
         !$OMP parallel do private(inode,iwa,iz,izz) ordered schedule(dynamic,1) num_threads(nodes)
         do iwa = 1,wblock%blen(iwb)
            inode=omp_get_thread_num()+1; 
            ww(inode) = wblock%os(iwb) + (iwa-1)*wblock%ds
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !! Back project image perturbation
            do iz=zblock%blen(izb),initz,-1
               izz=finz+iz-zblock%blen(izb)
               call owimage(.true., ww(inode), dws(:,:,inode),  bwr(:,:,izz,iwa), deltaI(:,:,:,:,iz), .true. )
               dws_z(:,:,inode)=dws_z(:,:,inode)+dws(:,:,inode)
               dws_z(:,:,inode)=dws_z(:,:,inode)*taper
               if (extrap(1:3).eq."ffd") then
                  if (ww(inode).ne.0) then
                     call owffd( ww(inode), inode,-1,dws_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode))
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (present(sref)) then
                     call owpspi( ww(inode), inode,-1,dws_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode), sref(:,:,izz-1))
                  else
                     call owpspi( ww(inode), inode,-1,dws_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode))
                  end if
               end if
               if (ww(inode).ne.0.) call owscctrx(ww(inode), 1,ddws(:,:,inode),bws(:,:,izz,iwa),vel(:,:,izz-1))

               GGG(:,:,izz-1,inode)=GGG(:,:,izz-1,inode)+conjg(ddws(:,:,inode))*dws_z(:,:,inode)

               call owimage(.true., ww(inode), bws(:,:,izz,iwa),  dwr(:,:,inode), deltaI(:,:,:,:,iz), .false.)
               dwr_z(:,:,inode)=dwr_z(:,:,inode)+dwr(:,:,inode)
               dwr_z(:,:,inode)=dwr_z(:,:,inode)*taper
               if (extrap(1:3).eq."ffd") then
                  if (ww(inode).ne.0) then
                     call owffd( ww(inode), inode,+1,dwr_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode))
                  end if
               else if (extrap(1:4).eq."pspi") then
                  if (present(sref)) then
                     call owpspi( ww(inode), inode,+1,dwr_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode), sref(:,:,izz-1))
                  else
                     call owpspi( ww(inode), inode,+1,dwr_z(:,:,inode), vel(:,:,izz-1), tflds(:,:,inode))
                  end if
               end if
               if (ww(inode).ne.0.) call owscctrx(ww(inode),-1,ddwr(:,:,inode),bwr(:,:,izz,iwa),vel(:,:,izz-1))

               GGG(:,:,izz-1,inode)=GGG(:,:,izz-1,inode)+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
      end do
      do i=1,nodes
         do iz=1,nz
            do ixw=1,kxm%n; do iyw=1,kym%n
               deltaS(ixw,iyw,iz)=deltaS(ixw,iyw,iz)+GGG(ixw,iyw,iz,i)
            end do; end do
         end do
      end do
   end if
end do
call MPI_REDUCE_(deltaS,deltaS0,kxm%n*kym%n*nz,MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
if(wpar%i_mpi == 0) then
   call grad_io("grad_","reed",deltaS0)
   call grad_io("grad_","rite",deltaS0)
end if
if (wpar%i_mpi > 0) call owfftdestroy()
write(0,*) " 13 passei" , wpar%i_mpi,idip

deallocate (deltaS)
if (wpar%i_mpi > 0) then
   deallocate (deltaI)
   deallocate(GGG)
   deallocate (ws,wr)
   deallocate (bws,bwr)
   deallocate (ddws,ddwr)
   deallocate (dws,dwr,tflds)
   deallocate (dws_z,dwr_z)
end if

6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!=============================================================
end module




!else
!   GG=0
!   allocate(dwr(1,kxm%n,kym%n,nodes));               dwr=0
!   allocate(dws(1,kxm%n,kym%n,nodes));               dws=0
!   !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
!    do iwa = 1,rec%w%n
!       inode=omp_get_thread_num()+1; 
!       ws(:,:,inode)=0.;wr(:,:,inode)=0.
!       ws(tapertrx+padtracesx+1:rec%hx%n+padtracesx+tapertrx,tapertry+padtracesy+1:rec%hy%n+padtracesy+tapertry,inode)=source(:,:,iwa)
!       wr(tapertrx+padtracesx+1:rec%hx%n+padtracesx+tapertrx,tapertry+padtracesy+1:rec%hy%n+padtracesy+tapertry,inode)=recev(:,:,iwa)
!       dws(:,:,:,inode)=0.
!       dwr(:,:,:,inode)=0.
!       bws(:,:,:,inode)=0.
!       bwr(:,:,:,inode)=0.
!       ddwr(:,:,node)=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), 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(  iwa, inode, 1,  ws(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou", sref(:,:,iz-1))
!               call owpspi(  iwa, inode,-1,  wr(:,:,inode), vel(:,:,iz-1), tfldr(:,:,inode), "rec", sref(:,:,iz-1))
!            else
!               call owpspi(  iwa, inode, 1,  ws(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou")
!               call owpspi(  iwa, inode,-1,  wr(:,:,inode), vel(:,:,iz-1), tfldr(:,:,inode), "rec")
!            end if
!         end if
!          ws(:,:,inode)=ws(:,:,inode)*taper
!          wr(:,:,inode)=wr(:,:,inode)*taper
!          bws(:,:,iz,inode)=ws(:,:,inode)
!          bwr(:,:,iz,inode)=wr(:,:,inode)
!       end do
!      do iz=max(2,nzmin),nzmax
!          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),vel(:,:,iz-1))
!          end if
!          dws(1,:,:,inode)=dws(1,:,:,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(1,:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou")
!             end if
!          else if (extrap(1:4).eq."pspi") then
!             if (present(sref)) then
!                call owpspi(  iwa, inode, 1, dws(1,:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou", sref(:,:,iz-1))
!             else
!                call owpspi(  iwa, inode, 1, dws(1,:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou")
!             end if
!          end if
!          dws(1,:,:,inode)=dws(1,:,:,inode)*taper
!
!          call owimage(.false., iwa, dws(1,:,:,inode), bwr(:,:,iz,inode), GG(:,:,:,:,iz,inode))
!
!          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),vel(:,:,iz-1))
!          end if
!          dwr(1,:,:,inode)=dwr(1,:,:,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(1,:,:,inode), vel(:,:,iz-1), tfldr(:,:,inode), "rec")
!             end if
!          else if (extrap(1:4).eq."pspi") then
!             if (present(sref)) then
!                call owpspi(  iwa, inode,-1, dwr(1,:,:,inode), vel(:,:,iz-1), tfldr(:,:,inode), "rec", sref(:,:,iz-1))
!             else
!                call owpspi(  iwa, inode,-1, dwr(1,:,:,inode), vel(:,:,iz-1), tfldr(:,:,inode), "rec")
!             end if
!          end if
!          dwr(1,:,:,inode)=dwr(1,:,:,inode)*taper
!
!          call owimage(.false., iwa, bws(:,:,iz,inode), dwr(1,:,:,inode), GG(:,:,:,:,iz,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
!   call owfftdestroy()
!   do i=1,nodes
!      !$OMP parallel do private(ixw,iyw,iz,ihx,ihy)
!      do iz=max(2,nzmin),nzmax
!         do ihy=1,nsecdimy
!            do iyw=1,kym%n
!               do ihx=1,nsecdimx
!                  do ixw=1,kxm%n
!                     deltaI(ixw,ihx,iyw,ihy,iz)=deltaI(ixw,ihx,iyw,ihy,iz)+GG(ixw,ihx,iyw,ihy,iz,i)
!                  end do
!               end do
!            end do
!         end do
!      end do
!      !$OMP END PARALLEL DO
!   end do
!   deallocate(dwr,dws)
!end if
!




!======================================================================
!subroutine owwei(taper,vel,sref)
!
!real, optional			:: sref(:,:,:)
!complex				:: recev(rec%hx%n,rec%hy%n,wblock%ns),source(sou%hx%n,sou%hy%n,wblock%ns)
!integer				:: iwa,iwb,izb,iz,inode,i,ixw,ihx,iyw,ihy,initz,finz
!real				:: taper(:,:),vel(kxm%n,kym%n,img%z%n),ww(wblock%ns)
!real				:: image(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns)
!real				:: GGG(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns,nodes)
!complex				:: ws(kxm%n,kym%n,nodes),wr(kxm%n,kym%n,nodes)
!complex				:: tflds(kxm%n,kym%n,nodes)
!integer, external			:: omp_get_thread_num,omp_get_num_threads
!integer                                 :: ninx,nfnx,niny,nfny
!
!!--------------------------------------------------
!! Setting FFTW plans
!tflds=0;call owfftws(tflds,nodes)
!ninx=tapertrx+padtracesx+1
!nfnx=rec%hx%n+padtracesx+tapertrx
!niny=tapertry+padtracesy+1
!nfny=rec%hy%n+padtracesy+tapertry
!do iwb=1, wblock%nb
!   RECV_=0;   call wavefield_io(iwb,  RECV_,"rec","reed")
!   source=0; call wavefield_io(iwb,source,"sou","reed")
!   ws=0;ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))=source(:,:,1:wblock%blen(iwb))
!   wr=0;wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))=recev(:,:,1:wblock%blen(iwb))
!   do izb=1,zblock%nb
!      GGG=0.;image=0;
!      if (izb==1) then
!         initz=2
!         finz=0
!      else
!         initz=1
!         finz=(izb-1)*zblock%blen(izb-1)
!      end if
!      !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(wblock%blen(iwb))
!      do iwa = 1,wblock%blen(iwb)
!         inode=omp_get_thread_num()+1; 
!         ww(iwa) = wblock%os(iwb) + (iwa-1)*wblock%ds
!         do iz=initz,zblock%blen(izb)
!            if (extrap(1:3).eq."ffd") then
!               if (ww(iwa).ne.0) then
!                  call owffd( ww(iwa), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,inode))
!                  call owffd( ww(iwa), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,inode))
!               end if
!            else if (extrap(1:4).eq."pspi") then
!               if (present(sref)) then
!                  call owpspi( ww(iwa), iwa, 1,  ws(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,iwa), sref(:,:,iz-1+finz))
!                  call owpspi( ww(iwa), iwa,-1,  wr(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,iwa), sref(:,:,iz-1+finz))
!               else
!                  call owpspi( ww(iwa), iwa, 1,  ws(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,iwa))
!                  call owpspi( ww(iwa), iwa,-1,  wr(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,iwa))
!               end if
!            end if
!            ws(:,:,iwa)=ws(:,:,iwa)*taper
!            wr(:,:,iwa)=wr(:,:,iwa)*taper
!            call owimage(.false., ww(iwa), ws(:,:,iwa),  wr(:,:,iwa), GGG(:,:,:,:,iz,iwa))
!         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
!      !$OMP parallel do private(i,ixw,iyw,iz,ihx,ihy) schedule(dynamic,1) num_threads(wblock%blen(iwb))
!      do i=1,wblock%blen(iwb)
!         do iz=1,zblock%blen(izb)
!            do ihy=1,nsecdimy
!               do iyw=1,kym%n
!                  do ihx=1,nsecdimx
!                     do ixw=1,kxm%n
!                        image(ixw,ihx,iyw,ihy,iz)=image(ixw,ihx,iyw,ihy,iz)+GGG(ixw,ihx,iyw,ihy,iz,i)
!                     end do
!                  end do
!               end do
!            end do
!         end do
!      end do
!      !$OMP END PARALLEL DO
!      call image_io("bkgimg","reed",izb,finz,image)
!      call image_io("bkgimg","rite",izb,finz,image)
!   end do
!end do
!6 format(A13,i3,A8,i3,A14,i3,A12,i3)
!call owfftdestroy()
!end subroutine
