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

implicit none

logical,private                              :: doscale
integer, private                             :: idip,nodes,ireal,itera
real,private,allocatable                     :: asc(:,:),asc1(:,:,:)

contains
!======================================================================
subroutine owwei_init(idip_in,ireal_in,iter_in)
integer                                 :: idip_in,ireal_in,iter_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
itera=iter_in
end subroutine
!======================================================================
subroutine owwei_nozb(taper,vel,image,sref)

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

!--------------------------------------------------
! Setting FFTW plans
   allocate(tflds(kxm%n,kym%n,nodes))
   allocate(ws(kxm%n,kym%n,wblock%ns),wr(kxm%n,kym%n,wblock%ns))
   call owexec_time(inic)
   tflds=0;call owfftws(tflds,nodes)
   call owexec_time(fini)

if (multsc) then
   if (itera < 4) then
      nfreqb=3 !2*itera
   else if (itera < 6 .and. itera > 3) then
      nfreqb=6 !2*itera
   else
      nfreqb=wblock%nb
   end if
else
   nfreqb=wblock%nb
end if
do iwb=1, nfreqb
      ws=0;wr=0
      call gathers_io(iwb,idip,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
      call gathers_io(iwb,idip,ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"sou","reed")
      write(0,*) "Read shot ",idip," in node ",wpar%i_mpi," freq.block",iwb
      !$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
         call owimage(.false., ww(inode), ws(:,:,iwa),  wr(:,:,iwa), image(:,:,:,:,1))
         do iz=2,img%z%n
            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
            call owimage(.false., ww(inode), ws(:,:,iwa),  wr(:,:,iwa), image(:,:,:,:,iz))
         end do
         if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",idip, &
   	" FROM NODE: ",inode
      end do
      !$OMP END PARALLEL DO
end do

   deallocate(ws,wr,tflds)
   call owfftdestroy()

6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!!======================================================================
subroutine owdso_nozb0(pick1,wgcig,fun)
real                                                        :: pick1(:)
real    					:: wgcig(:,:,:,:,:)
integer						:: ix,ihx,iy,ihy,izb,finz,iz,k,ierr,i_r,status
real						:: hx, hy, hhmax,hh,zm
double precision                                :: fun
integer                                         :: nx_max,nx_min,ny_max,ny_min
   fun=0
   hx=max(abs(osecdimx),(nsecdimx-1)*dsecdimx+osecdimx)
   hy=max(abs(osecdimy),(nsecdimy-1)*dsecdimy+osecdimy)
   hhmax=sqrt(hx*hx+hy*hy)
   nx_max=(xv_max-img%xm%o)/img%xm%d+1.5
   nx_min=(xv_min-img%xm%o)/img%xm%d+1.5
   ny_max=(yv_max-img%ym%o)/img%ym%d+1.5
   ny_min=(yv_min-img%ym%o)/img%ym%d+1.5


      do iz=1,img%z%n
         zm=(iz-1)*img%z%d+zmin
         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)
                  do ix=1,img%xm%n
                     wgcig(ix,ihx,iy,ihy,iz)=wgcig(ix,ihx,iy,ihy,iz)*hh**2
                     if (zm < int(pick1((ix-1)*img%ym%n+iy))) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix>=nx_max.or.ix<=nx_min) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix<=nx_max.and.ix>=nx_min) then
                        fun=fun+wgcig(ix,ihx,iy,ihy,iz)**2
                     end if
                  end do
               end do
            end do
         end do
      end do
   
end subroutine
!!======================================================================
subroutine owdso_nozb1(pick1,wgcig,fun)
real                                                        :: pick1(:)
real        					:: wgcig(:,:,:,:,:)
real, allocatable                               :: ogcig(:,:,:,:)
integer						:: ix,ihx,iy,ihy,izb,finz,iz,k,ierr,i_r,status
real						:: hx, hy, hhmax,hh,zm
double precision                                            :: fun,a
integer                                         :: nx_max,nx_min,ny_max,ny_min

   allocate(ogcig(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)
   nx_max=(xv_max-img%xm%o)/img%xm%d+1.5
   nx_min=(xv_min-img%xm%o)/img%xm%d+1.5
   ny_max=(yv_max-img%ym%o)/img%ym%d+1.5
   ny_min=(yv_min-img%ym%o)/img%ym%d+1.5
   fun=0.
      do iz=1,img%z%n
         ogcig=0
         zm=(iz-1)*img%z%d+zmin
         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)
                  do ix=1,img%xm%n
                     if (hx < 0) then
                        ogcig(ix,ihx,iy,ihy)=(wgcig(ix,ihx+1,iy,ihy,iz)-wgcig(ix,ihx,iy,ihy,iz))/dsecdimx
                     else if (hx > 0) then
                        ogcig(ix,ihx,iy,ihy)=-(wgcig(ix,ihx,iy,ihy,iz)-wgcig(ix,ihx-1,iy,ihy,iz))/dsecdimx
                     else
                        ogcig(ix,ihx,iy,ihy)=(wgcig(ix,ihx+1,iy,ihy,iz)-wgcig(ix,ihx-1,iy,ihy,iz))/2/dsecdimx
                     end if
                     if (zm < int(pick1((ix-1)*img%ym%n+iy))) then
                        ogcig(ix,ihx,iy,ihy)=0
                     end if
                  end do
               end do
            end do
         end do
         wgcig(:,:,:,:,iz)=ogcig
      end do

      do iz=1,img%z%n
         do ihy=1,nsecdimy
            do iy=1,img%ym%n
               do ihx=1,nsecdimx
                  do ix=1,img%xm%n
                     if (ix>=nx_max.or.ix<=nx_min) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix<=nx_max.and.ix>=nx_min) then
                        fun=fun+wgcig(ix,ihx,iy,ihy,iz)**2
                     end if
                  end do
               end do
            end do
         end do
      end do

   deallocate(ogcig)
   
end subroutine
!!======================================================================
subroutine owdso_nozb2(pick1,wgcig,fun)
real                                                        :: pick1(:)
real    					:: wgcig(:,:,:,:,:)
integer						:: ix,ihx,iy,ihy,izb,finz,iz,k,ierr,i_r,status
real						:: hx, hy, hhmax,hh,zm
double precision                                :: fun
integer                                         :: nx_max,nx_min,ny_max,ny_min
   fun=0
   hx=max(abs(osecdimx),(nsecdimx-1)*dsecdimx+osecdimx)
   hy=max(abs(osecdimy),(nsecdimy-1)*dsecdimy+osecdimy)
   hhmax=sqrt(hx*hx+hy*hy)
   nx_max=(xv_max-img%xm%o)/img%xm%d+1.5
   nx_min=(xv_min-img%xm%o)/img%xm%d+1.5
   ny_max=(yv_max-img%ym%o)/img%ym%d+1.5
   ny_min=(yv_min-img%ym%o)/img%ym%d+1.5


      do iz=1,img%z%n
         zm=(iz-1)*img%z%d+zmin
         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)*sqrt(abs(cos(hx/hhmax*pi/2)))
                  do ix=1,img%xm%n
                     wgcig(ix,ihx,iy,ihy,iz)=wgcig(ix,ihx,iy,ihy,iz)*hh**2
                     if (zm < int(pick1((ix-1)*img%ym%n+iy))) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix>=nx_max.or.ix<=nx_min) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix<=nx_max.and.ix>=nx_min) then
                        fun=fun+wgcig(ix,ihx,iy,ihy,iz)**2
                     end if
                  end do
               end do
            end do
         end do
      end do
   
end subroutine
!!======================================================================
subroutine owdso_nozb3(pick1,wgcig,fun)
real                                                        :: pick1(:)
real        					:: wgcig(:,:,:,:,:)
real, allocatable                               :: ogcig(:,:,:,:)
integer						:: ix,ihx,iy,ihy,izb,finz,iz,k,ierr,i_r,status
real						:: hx, hy, hhmax,hh,zm
double precision                                            :: fun,a
integer                                         :: nx_max,nx_min,ny_max,ny_min

   allocate(ogcig(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)
   nx_max=(xv_max-img%xm%o)/img%xm%d+1.5
   nx_min=(xv_min-img%xm%o)/img%xm%d+1.5
   ny_max=(yv_max-img%ym%o)/img%ym%d+1.5
   ny_min=(yv_min-img%ym%o)/img%ym%d+1.5
   fun=0.
      do iz=1,img%z%n
         ogcig=0
         zm=(iz-1)*img%z%d+zmin
         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)*sqrt(abs(cos(hx/hhmax*pi/2)))
                  do ix=1,img%xm%n
                     if (hx < 0) then
                        ogcig(ix,ihx,iy,ihy)=(wgcig(ix,ihx+1,iy,ihy,iz)-wgcig(ix,ihx,iy,ihy,iz))/dsecdimx
                     else if (hx > 0) then
                        ogcig(ix,ihx,iy,ihy)=-(wgcig(ix,ihx,iy,ihy,iz)-wgcig(ix,ihx-1,iy,ihy,iz))/dsecdimx
                     else
                        ogcig(ix,ihx,iy,ihy)=(wgcig(ix,ihx+1,iy,ihy,iz)-wgcig(ix,ihx-1,iy,ihy,iz))/2/dsecdimx
                     end if
                     if (zm < int(pick1((ix-1)*img%ym%n+iy))) then
                        ogcig(ix,ihx,iy,ihy)=0
                     end if
                  end do

               end do
            end do
         end do
         wgcig(:,:,:,:,iz)=ogcig
      end do
      do iz=1,img%z%n
         do ihy=1,nsecdimy
            do iy=1,img%ym%n
               do ihx=1,nsecdimx
                  hx=(ihx-1)*dsecdimx+osecdimx
                  hh=(sqrt(hx*hx+hy*hy)/hhmax)*sqrt(abs(cos(hx/hhmax*pi/2)))
                  do ix=1,img%xm%n
                     wgcig(ix,ihx,iy,ihy,iz)=wgcig(ix,ihx,iy,ihy,iz)*hh**2
                     if (ix>=nx_max.or.ix<=nx_min) then
                        wgcig(ix,ihx,iy,ihy,iz)=0
                     end if
                     if (ix<=nx_max.and.ix>=nx_min) then
                        fun=fun+wgcig(ix,ihx,iy,ihy,iz)**2
                     end if
                  end do
               end do
            end do
         end do
      end do

   deallocate(ogcig)
   
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(taper,vel,iter,deltaI,deltaS,sref)

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

   allocate (ws(kxm%n,kym%n,wblock%ns),wr(kxm%n,kym%n,wblock%ns));ws=0;wr=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

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

if (multsc) then
   if (itera < 4) then
      nfreqb=3
   else if (itera < 6 .and. itera > 3) then
      nfreqb=6
   else
      nfreqb=wblock%nb
   end if
else
   nfreqb=wblock%nb
end if
do iwb=1, nfreqb
      ws=0;wr=0
      call gathers_io(iwb,idip,wr(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"rec","reed")
      call gathers_io(iwb,idip,ws(ninx:nfnx,niny:nfny,1:wblock%blen(iwb)),"sou","reed")

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

      dws=0.;dwr=0.
      dws_z=0.;dwr_z=0.
      ddwr=0.;ddws=0.
      GGG=0.
      !$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=nz,2,-1
            call owimage(.true., ww(inode), dws(:,:,inode),  bwr(:,:,iz,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(:,:,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 (ww(inode).ne.0.) call owscctrx(ww(inode), 1,ddws(:,:,inode),bws(:,:,iz,iwa),vel(:,:,iz-1))

            GGG(:,:,iz,inode)=GGG(:,:,iz,inode)+conjg(ddws(:,:,inode))*dws_z(:,:,inode)

            call owimage(.true., ww(inode), bws(:,:,iz,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(:,:,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 (ww(inode).ne.0.) call owscctrx(ww(inode),-1,ddwr(:,:,inode),bwr(:,:,iz,iwa),vel(:,:,iz-1))

            GGG(:,:,iz,inode)=GGG(:,:,iz,inode)+conjg(ddwr(:,:,inode))*dwr_z(:,:,inode)
         end do
         if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",idip, &
   	" FROM NODE: ",inode
      end do
      !$OMP END PARALLEL DO

      do i=1,nodes
         deltaS=deltaS+GGG(:,:,:,i)
      end do
end do

call owfftdestroy()

   deallocate(GGG)
   deallocate (bws,bwr)
   deallocate (ws,wr)
   deallocate (ddws,ddwr)
   deallocate (dws,dwr,tflds)
   deallocate (dws_z,dwr_z)

6 format(A13,i3,A8,i3,A14,i3,A12,i3)
end subroutine
!========================================
subroutine gathers_io(iwb,isht,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwb,isht
character*3          :: tag
character*4          :: op
!write(0,*) "#####  ",iwb,isht,op,ireal-1,"  #####  "
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 

  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) 

else

  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) 

end if
end subroutine
!========================================
subroutine wavefield_io2(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,(/kxm%n,kym%n,            img%z%n,            sou%w%n/),&
                          (/kxm%n,kym%n,   zblock%blen(izb),   wblock%blen(iwb)/),&
                          (/    0,    0, zblock%indx(izb)-1, wblock%indx(iwb)-1/),&
                          (/    1,    1,                  1,                  1/),&
                          8,GG(:,:,1:zblock%blen(izb),1:wblock%blen(iwb))) 
else
  call srite_window(tag,4,(/kxm%n,kym%n,            img%z%n,            sou%w%n/),&
                          (/kxm%n,kym%n,   zblock%blen(izb),   wblock%blen(iwb)/),&
                          (/    0,    0, zblock%indx(izb)-1, wblock%indx(iwb)-1/),&
                          (/    1,    1,                  1,                  1/),&
                          8,GG(:,:,1:zblock%blen(izb),1:wblock%blen(iwb))) 
end if
end subroutine
!========================================
subroutine image_io(tag,op,izb,image,soma)
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              :: soma

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(soma) image(:,:,:,:,1:zblock%blen(izb))=image(:,:,:,:,1:zblock%blen(izb))+temp
  if(not(soma)) 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,image,soma)
real                          :: image(:,:,:,:,:)
real,allocatable              :: temp(:,:,:,:,:)
integer                       :: nn1,nn2,nn3,nn4,nn5
character*4                   :: op
character*6                   :: tag
logical,optional              :: soma

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,nn5));temp=0.
  call sreed_window(trim(tag),5,(/nn1,nn2,nn3,nn4,nn5/),&
                                (/nn1,nn2,nn3,nn4,nn5/),&
                                (/  0,  0,  0,  0,  0/),&
                                (/  1,  1,  1,  1,  1/),&
                                4,temp) 
  if(soma) image=image+temp
  if(not(soma)) image=temp
  deallocate(temp)
else if(op.eq."rite") then
  call srite_window(trim(tag),5,(/nn1,nn2,nn3,nn4,nn5/),&
                                (/nn1,nn2,nn3,nn4,nn5/),&
                                (/  0,  0,  0,  0,  0/),&
                                (/  1,  1,  1,  1,  1/),&
                                4,image) 
end if
end subroutine
!========================================
subroutine image_io2(tag,op,iz,image)
real                          :: image(:,:,:,:)
integer                       :: nn1,nn2,nn3,nn4,nn5,iz
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,         iz-1/),&
                          (/  0,  0,  0,  0,            0/),&
                          (/  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,         iz-1/),&
                          (/  0,  0,  0,  0,            0/),&
                          (/  1,  1,  1,  1,            1/),&
                          4,image) 
end if
end subroutine
!========================================
subroutine grad_io(tag,op,nit,image,soma)
logical,optional              :: soma
real                          :: image(kxm%n,kym%n,nz)
real, allocatable             :: tmp(:,:,:)
integer                       :: nn1,nn2,nn3,nit
character*4                   :: op
character*5                   :: tag

nn1=kxm%n
nn2=kym%n
nn3=nz

if(op.eq."reed") then
  allocate(tmp(nn1,nn2,nn3))
  call sreed_window(trim(tag),4,(/nn1,nn2,nn3,max_iters/),&
                                (/nn1,nn2,nn3,1/),&
                                (/  0,  0,  0,nit-1/),&
                                (/  1,  1,  1,1/),&
                                         4,tmp) 
  if(soma) then
    image=image+tmp
  else
    image=tmp
  end if
  deallocate(tmp)
else if(op.eq."rite") then
  call srite_window(trim(tag),4,(/nn1,nn2,nn3,max_iters/),&
                                (/nn1,nn2,nn3,1/),&
                                (/  0,  0,  0,nit-1/),&
                                (/  1,  1,  1,1/),&
                                         4,image) 
end if
end subroutine
!========================================
subroutine grad_io1(tag,op,izb,image,soma)
real                          :: image(:,:,:)
real,allocatable              :: temp(:,:,:)
integer                       :: izb
real                          :: tmp,a,b
integer                       :: i,nn1,nn2,nn3,nn4,nn5
character*4                   :: op
character*4                   :: tag
logical,optional              :: soma

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

if(op.eq."reed") then
  allocate(temp(nn1,nn2,zblock%blen(izb)));temp=0.
  call sreed_window(tag,3,(/nn1,nn2,           img%z%n/),&
                          (/nn1,nn2,  zblock%blen(izb)/),&
                          (/  0,  0,zblock%indx(izb)-1/),&
                          (/  1,  1,                 1/),&
                          4,temp) 
  if(soma) image(:,:,1:zblock%blen(izb))=image(:,:,1:zblock%blen(izb))+temp
  if(not(soma)) image(:,:,1:zblock%blen(izb))=temp
  deallocate(temp)
else if(op.eq."rite") then
  call srite_window(tag,3,(/nn1,nn2,           img%z%n/),&
                          (/nn1,nn2,  zblock%blen(izb)/),&
                          (/  0,  0,zblock%indx(izb)-1/),&
                          (/  1,  1,                 1/),&
                          4,image(:,:,1:zblock%blen(izb)))
end if
end subroutine
!=============================================================
subroutine ow_scale(data)
real                                :: data(:,:,:,:,:),a(img%xm%n,img%ym%n),temp(img%xm%n,img%ym%n)
integer                             :: i0,i1,i2,i3,i4,i5

  a=0
  do i5=1,img%z%n;do i4=1,img%shy%n;do i3=1,img%ym%n;do i2=1,img%shx%n;do i1=1,img%xm%n
     a(i1,i3) = a(i1,i3) + data(i1,i2,i3,i4,i5)**2
  end do;end do;end do;end do;end do
  a = sqrt(a/(img%z%n*img%shy%n*img%ym%n*img%shx%n*img%xm%n))
  temp=0;
  do i3=1,img%ym%n;do i1=1,img%xm%n
     temp(i1,i3)=sum(a(max(1,i1-5):min(img%xm%n,i1+5),i3))/(min(img%xm%n,i1+5)-max(1,i1-5)+1)
  end do;end do
  do i3=1,img%ym%n;do i1=1,img%xm%n
     a(i1,i3)=sum(temp(i1,max(1,i3-5):min(img%ym%n,i3+5)))/(min(img%ym%n,i3+5)-max(1,i3-5)+1)
  end do;end do
  do i3=1,img%ym%n;do i1=1,img%xm%n
     if(a(i1,i3).ne.0)  data(i1,:,i3,:,:)=data(i1,:,i3,:,:)/a(i1,i3)
  end do;end do

  call to_history("n1",img%xm%n,"scale"); call to_history("d1",img%xm%d,"scale"); call to_history("o1",img%xm%o,"scale")  
  call to_history("n2",img%ym%n,"scale"); call to_history("d2",img%ym%d,"scale"); call to_history("o2",img%ym%o,"scale")  
  call to_history("n3",1,"scale") ; call to_history("d3",1.,"scale") ; call to_history("o5",1.,"scale")
  call to_history("esize",4,"scale")
  call srite("scale",a,4*img%xm%n*img%ym%n)

end subroutine
!=============================================================
subroutine scale_init()
  allocate(asc(img%xm%n,img%ym%n))
  allocate(asc1(img%xm%n,img%ym%n,img%z%n))
end subroutine
!=============================================================
subroutine scale_init1(doscale_in)
  logical    :: doscale_in
  doscale=doscale_in
end subroutine
!=============================================================
subroutine scale_close()
  deallocate(asc,asc1)
end subroutine
!=============================================================
subroutine ow_scale1(nnn,data,scalegr)
real, optional                      :: scalegr(img%xm%n,img%ym%n,img%z%n)
real                                :: data(:,:,:,:,:),temp(img%xm%n,img%ym%n,img%z%n)
integer                             :: i1,i2,i3,i4,i5,izb,ierr,nnn
real                                :: b


  if (multsc) then
     if (nnn==1.or.nnn==4.or.nnn==6.and.doscale) then
        asc=0;
        do i5=1,img%z%n
           do i4=1,img%shy%n;do i3=1,img%ym%n;do i2=1,img%shx%n;do i1=1,img%xm%n
              asc(i1,i3) = asc(i1,i3) + data(i1,i2,i3,i4,i5)**2
           end do;end do;end do;end do
        end do
        asc = sqrt(asc/(img%z%n*img%shy%n*img%ym%n*img%shx%n*img%xm%n))
        b=0
        do i3=1,img%ym%n;do i1=1,img%xm%n
          if (asc(i1,i3)>b) b=asc(i1,i3)
        end do;end do
        asc=asc/b
      end if
  else if (onlygrad) then
!      asc=0;
!      do i5=1,img%z%n
!         do i4=1,img%shy%n;do i3=1,img%ym%n;do i2=1,img%shx%n;do i1=1,img%xm%n
!            asc(i1,i3) = asc(i1,i3) + data(i1,i2,i3,i4,i5)**2
!         end do;end do;end do;end do
!      end do
!      asc = sqrt(asc/(img%z%n*img%shy%n*img%ym%n*img%shx%n*img%xm%n))
!      b=0
!      do i3=1,img%ym%n;do i1=1,img%xm%n
!        if (asc(i1,i3)>b) b=asc(i1,i3)
!      end do;end do
!      asc=asc/b
!      temp=0;
!      do i3=1,img%ym%n;do i1=1,img%xm%n
!         temp(i1,i3)=sum(asc(max(1,i1-25):min(img%xm%n,i1+25),i3))/(min(img%xm%n,i1+25)-max(1,i1-25)+1)
!      end do;end do
!      do i3=1,img%ym%n;do i1=1,img%xm%n
!         asc(i1,i3)=sum(temp(i1,max(1,i3-25):min(img%ym%n,i3+25)))/(min(img%ym%n,i3+25)-max(1,i3-25)+1)
!      end do;end do
!      scalegr=asc
      asc1=0;
      do i5=1,img%z%n
         do i4=1,img%shy%n;do i3=1,img%ym%n;do i2=1,img%shx%n;do i1=1,img%xm%n
            asc1(i1,i3,i5) = asc1(i1,i3,i5) + data(i1,i2,i3,i4,i5)**2
         end do;end do;end do;end do
      end do
      asc1 = sqrt(asc1/(img%shy%n*img%shx%n))
      temp=0;
      do i5=1,img%z%n;do i3=1,img%ym%n;do i1=1,img%xm%n
         temp(i1,i3,i5)=sum(asc1(max(1,i1-25):min(img%xm%n,i1+25),i3,i5))/(min(img%xm%n,i1+25)-max(1,i1-25)+1)
      end do;end do;end do
      do i5=1,img%z%n;do i3=1,img%ym%n;do i1=1,img%xm%n
         asc1(i1,i3,i5)=sum(temp(i1,max(1,i3-25):min(img%ym%n,i3+25),i5))/(min(img%ym%n,i3+25)-max(1,i3-25)+1)
      end do;end do;end do
      do i5=1,img%z%n;do i3=1,img%ym%n;do i1=1,img%xm%n
         temp(i1,i3,i5)=sum(asc1(i1,i3,max(1,i5-100):min(img%z%n ,i5+100)))/(min(img%z%n, i5+100)-max(1,i5-100)+1)
      end do;end do;end do
      b=0
      do i5=1,img%z%n;do i3=1,img%ym%n;do i1=1,img%xm%n
        if (asc1(i1,i3,i5)>b) b=asc1(i1,i3,i5)
      end do;end do;end do
      scalegr=asc1/b
  else 
     if (nnn==1.and.doscale) then
        asc=0;
        do i5=1,img%z%n
           do i4=1,img%shy%n;do i3=1,img%ym%n;do i2=1,img%shx%n;do i1=1,img%xm%n
              asc(i1,i3) = asc(i1,i3) + data(i1,i2,i3,i4,i5)**2
           end do;end do;end do;end do
        end do
        asc = sqrt(asc/(img%z%n*img%shy%n*img%ym%n*img%shx%n*img%xm%n))
        b=0
        do i3=1,img%ym%n;do i1=1,img%xm%n
          if (asc(i1,i3)>b) b=asc(i1,i3)
        end do;end do
        asc=asc/b
      end if
  end if

  do i3=1,img%ym%n;do i1=1,img%xm%n
     if(asc(i1,i3).ne.0)  data(i1,:,i3,:,:)=data(i1,:,i3,:,:)/asc(i1,i3)
  end do;end do

end subroutine
!=============================================================

end module

