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(:,:,:)
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,shtidx,shtblm,i_r
integer                                 :: ierr,status
real				:: taper(:,:),vel(kxm%n,kym%n,img%z%n),ww(nodes)
real,allocatable                        :: image(:,:,:,:,:),image0(:,:,:,:,:),GGG(:,:,:,:,:,:)
complex				:: ws(kxm%n,kym%n,wblock%ns),wr(kxm%n,kym%n,wblock%ns)
complex				:: tflds(kxm%n,kym%n,nodes)
integer, external			:: omp_get_thread_num,omp_get_num_threads

allocate(image(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns))
!--------------------------------------------------
! Setting FFTW plans
if (wpar%i_mpi > 0) then
   tflds=0;call owfftws(tflds,nodes)
end if 
if (shtblk - idip +1 <= shtres) then
   shtidx = wpar%n_mpi
else
   shtidx = wpar%n_mpi-1
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((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
         call gathers_io((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)
      end do
      shtblm=shtidx
   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)
      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(nodes)
         do iwa = 1,wblock%blen(iwb)
            inode=omp_get_thread_num()+1; 
            ww(inode) = wblock%os(iwb) + (iwa-1)*wblock%ds
            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+finz), tflds(:,:,inode))
                     call owffd( ww(inode), 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(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,inode)&
                                , sref(:,:,iz-1+finz))
                     call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,inode)&
                                , sref  (:,:,iz-1+finz))
                  else
                     call owpspi( ww(inode), inode, 1,  ws(:,:,iwa), vel(:,:,iz-1+finz), tflds(:,:,inode))
                     call owpspi( ww(inode), inode,-1,  wr(:,:,iwa), vel(:,:,iz-1+finz), 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(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(nodes)
         do i=1,nodes
            do iz=initz,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
      if (allocated(GGG)) deallocate(GGG)
   end if
end do
if (allocated(image0)) deallocate(image0); allocate(image0(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns))
do izb=1,zblock%nb
   if (izb==1) then
      finz=0
   else
      finz=(izb-1)*zblock%blen(izb-1)
   end if
   if(wpar%i_mpi > 0) call image_io("bkgimg","reed",izb,finz,image)
   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) then
      call image_io("bkgimg","rite",izb,finz,image0)
   end if
!write(0,*) " 11 passei" , wpar%i_mpi
end do
deallocate(image,image0)
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
if (wpar%i_mpi > 0) call owfftdestroy()
!write(0,*) " 12 passei" , wpar%i_mpi
end subroutine
!========================================
subroutine gathers_io(iwbb,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwbb
character*3          :: tag
character*4          :: op

if(op.eq."reed") then 
!write(0,*) "before read gathers io", iwbb, wblock%blen(iwbb), wpar%i_mpi
  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(iwbb),       1,      1/),&
                          (/       0,       0, wblock%indx(iwbb)-1,  idip-1,ireal-1/),&
                          (/       1,       1,                   1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwbb))) 
!write(0,*) "after read gathers io", iwbb, wblock%blen(iwbb), wpar%i_mpi
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(iwbb),       1,      1/),&
                          (/       0,       0, wblock%indx(iwbb)-1,  idip-1,ireal-1/),&
                          (/       1,       1,                   1,       1,      1/),&
                          8,GG(:,:,1:wblock%blen(iwbb))) 
end if
end subroutine
!========================================
subroutine wavefield_io(iwb,izb,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwb,izb
character*3          :: 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,izb,GG,tag,op)
complex              :: GG(:,:)
integer              :: iwb,izb
character*3          :: tag
character*4          :: op

if(op.eq."reed") then 
  call sreed_window(tag,4,(/rec%hx%n,rec%hy%n, img%z%n,           sou%w%n/),&
                          (/rec%hx%n,rec%hy%n,       1,                 1/),&
                          (/       0,       0,   izb-1, (iwb-1)*wblock%ns/),&
                          (/       1,       1,       1,                 1/),&
                          8,GG) 
else
  call srite_window(tag,4,(/rec%hx%n,rec%hy%n, img%z%n,           sou%w%n/),&
                          (/rec%hx%n,rec%hy%n,       1,                 1/),&
                          (/       0,       0,   izb-1, (iwb-1)*wblock%ns/),&
                          (/       1,       1,       1,                 1/),&
                          8,GG) 
end if
end subroutine
!========================================
subroutine image_io(tag,op,izb,finz,image)
real,optional                 :: image(:,:,:,:,:)
real,allocatable              :: temp(:,:,:,:,:)
integer,optional              :: izb,finz
real                          :: tmp,a,b
integer                       :: i,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)
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,              finz/),&
                          (/  1,  1,  1,  1,                 1/),&
                          4,temp) 
  image(:,:,:,:,1:zblock%blen(izb))=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,              finz/),&
                          (/  1,  1,  1,  1,                 1/),&
                          4,image(:,:,:,:,1:zblock%blen(izb)))
else
  tmp=0.
  do i=1,nn1*nn2*nn3*nn4*img%z%n
    call srite(tag,tmp,4)
  end do
end if
end subroutine
!========================================
subroutine image_io1(tag,op,izb,image)
real,optional                 :: image(:,:,:,:)
real,allocatable              :: temp(:,:,:,:)
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,      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,      1/),&
                          (/  1,  1,  1,  1,      1/),&
                          4,image) 
end if
end subroutine
!!======================================================================
subroutine owdso(fun,riteD,pick1)
real,optional                                               :: pick1(:)
real						:: wgcig(kxm%n,nsecdimx,kym%n,nsecdimy,zblock%ns)
integer						:: ix,ihx,iy,ihy,izb,finz,iz
real						:: hx, hy, hhmax,hh
double precision                                            :: fun
logical                                                     :: riteD

   hx=max(abs(osecdimx),(nsecdimx-1)*dsecdimx+osecdimx)
   hy=max(abs(osecdimy),(nsecdimy-1)*dsecdimy+osecdimy)
   hhmax=sqrt(hx*hx+hy*hy)
   call sseek("bkgimg",0,0)
   call sseek("Dimage",0,0)
   do izb=1,zblock%nb
      wgcig=0
      if (izb==1) then
         finz=0
      else
         finz=(izb-1)*zblock%blen(izb-1)
      end if
      call image_io("bkgimg","reed",izb,finz,wgcig)
      do iz=1,zblock%blen(izb)
         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 ((iz-1+finz)*img%z%d+zblock%os(izb) < int(pick1((iy-1)*img%xm%n+ix))) &
                            wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     wgcig(ix,ihx,iy,ihy,iz)=wgcig(ix,ihx,iy,ihy,iz)*hh
                     fun = fun + wgcig(ix,ihx,iy,ihy,iz)*wgcig(ix,ihx,iy,ihy,iz)
                  end do
               end do
            end do
         end do
      end do
      if (riteD) call image_io("Dimage","rite",izb,finz,wgcig)
   end do

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(deltaS,taper,vel,sref)

real, optional					:: sref(:,:,:)
logical						:: adj
integer						:: iwa,iz,inode,i,ixw,ish,n0,iyw,ihx,ihy
real						:: deltaS(kxm%n,kym%n,img%z%n)
real						:: deltaI(kxm%n,nsecdimx,kym%n,nsecdimy)
real                                                        :: taper(kxm%n,kym%n),vel(kxm%n,kym%n,img%z%n)
real						:: GGG(kxm%n,kym%n,img%z%n,nodes),ww(nodes)
complex						:: ws(kxm%n,kym%n,nodes),wr(kxm%n,kym%n,nodes)
complex						:: ddws(kxm%n,kym%n,nodes),ddwr(kxm%n,kym%n,nodes)
complex						:: dws(kxm%n,kym%n,nodes),dwr(kxm%n,kym%n,nodes)
complex						:: tflds(kxm%n,kym%n,nodes)
complex						:: dws_z(kxm%n,kym%n,nodes),dwr_z(kxm%n,kym%n,nodes)
integer, external					:: omp_get_thread_num,omp_get_num_threads
integer                                                     :: ninx,niny,nfnx,nfny
call sseek("bkws",0,0)
call sseek("bkwr",0,0)
call sseek("Dimage",0,0)

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

call sseek("bkws",0,0)
call sseek("bkwr",0,0)
call sseek("Dimage",0,0)
!if (adj) then
 do iwb=1, wblock%nb
   wr=0;call gathers_io(iwb,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
   ws=0;call gathers_io(iwb,ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"sou","reed")
   GGG=0.
  !$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; 
         dws(:,:,inode)=0.
         dwr(:,:,inode)=0.
         dws_z(:,:,inode)=0.
         dwr_z(:,:,inode)=0.
         ddwr(:,:,node)=0
         ddws(:,:,inode)=0
         ww(inode) = wblock%os(iwb) + (iwa-1)*wblock%ds
!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Background wavefields
!! bws and bwr
!!
        !$OMP ORDERED
         call wavefield_io1(iwb,1,ws(:,:,iwa),"bkws","rite")
         call wavefield_io1(iwb,1,wr(:,:,iwa),"bkwr","rite")
        !$OMP END ORDERED
         do iz=2,nz
            if (extrap(1:3).eq."ffd") then
               if (ww(iwa).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
           !$OMP ORDERED
            call wavefield_io1(iwb,iz,ws(:,:,iwa),"bkws","rite")
            call wavefield_io1(iwb,iz,wr(:,:,iwa),"bkwr","rite")
           !$OMP END ORDERED
         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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Back project image perturbation
       ws=0;wr=0
       do iz=img%z%n,max(2,nzmin),-1
          !$OMP ORDERED
           call wavefield_io1(iwb,iz,wr(:,:,iwa),"bkwr","reed")
           call wavefield_io1(iwb,iz,ws(:,:,iwa),"bkws","reed")
           call image_io1("Dimage","reed",iz,deltaI)
          !$OMP END ORDERED
          call owimage(.true., ww(inode), dws(:,:,inode),  wr(:,:,iwa), deltaI, .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(iwa).ne.0) then
                call owffd( ww(inode), inode,-1,dws_z(:,:,inode), 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,dws_z(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), sref(:,:,iz-1))
             else
                call owpspi( ww(inode), inode,-1,dws_z(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode))
             end if
          end if
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(ww(iwa), 1,ddws(:,:,inode),ws(:,:,iwa),vel(:,:,iz-1))
          end if
          if (renorm) then
             ws(:,:,iwa)=ws(:,:,iwa)+ddws(:,:,inode)
          end if
          GGG(:,:,iz-1,inode)=GGG(:,:,iz-1,inode)+conjg(ddws(:,:,inode))*dws_z(:,:,inode)

          call owimage(.true., ww(inode), ws(:,:,inode),  dwr(:,:,inode), deltaI, .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(iwa).ne.0) then
                call owffd( ww(inode), inode,+1,dwr_z(:,:,inode), 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,dwr_z(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), sref(:,:,iz-1))
             else
                call owpspi( ww(inode), inode,+1,dwr_z(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode))
             end if
          end if
          if (rec%w%all(iwa).ne.0.) then
             call owscctrx(ww(iwa),-1,ddwr(:,:,inode),wr(:,:,iwa),vel(:,:,iz-1))
          end if
          if (renorm) then
             wr(:,:,iwa)=wr(:,:,iwa)+ddwr(:,:,inode)
          end if
          GGG(:,:,iz-1,inode)=GGG(:,:,iz-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
   call owfftdestroy()
   do i=1,nodes
      !$OMP parallel do private(ixw,iz,iyw)
      do iz=max(2,nzmin),img%z%n !max(1,nzmin),nzmax
         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
      !$OMP END PARALLEL DO
   end do

!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
!
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!=============================================================
end module




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