module OW_wei_areal2

use sep
use OWvelo_areal
use OW_type
use OW_image_areal
use OW_fftw_areal
use OW_parms_areal

implicit none

contains

subroutine owwei_f(taper,recv,source,sref,vel,deltaS,deltaI,GGr,GGs)

complex,optional					:: GGr(:,:,:),GGs(:,:,:)
complex						:: recv(:,:),source(:,:)
integer						:: iwa,iz,ikx,iref,ixw,j,inode,i,ish,nrefvel,nodes
real						:: sqa,taper(:),vel(:,:),sref(:,:,:),a
real						:: c(5),deltaS(:,:),deltaI(:,:,:)
integer,allocatable					:: indw(:,:)
real, allocatable					:: GGG(:,:,:,:),sqa0(:,:)
complex, allocatable				:: ws(:,:),wr(:,:),gs0(:,:),gr0(:,:)
complex, allocatable				:: ddws(:,:),ddwr(:,:),gs(:,:),gr(:,:)
complex, allocatable				:: dws(:,:),dwr(:,:),dgs0(:,:),dgr0(:,:)
complex, allocatable				:: dwss(:,:),dwrr(:,:),dgs(:,:),dgr(:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

if (node.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
else
  nodes=node
end if
write(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"

allocate(wr(kxm%n,nodes));                   wr=0
allocate(ws(kxm%n,nodes));                   ws=0
allocate(ddwr(kxm%n,nodes));                 ddwr=0
allocate(ddws(kxm%n,nodes));                 ddws=0
allocate(gs0(kxm%n,nodes));                  gs0=0
allocate(gr0(kxm%n,nodes));                  gr0=0
allocate(dwr(kxm%n,nodes));                  dwr=0
allocate(dws(kxm%n,nodes));                  dws=0
allocate(dwrr(kxm%n,nodes));                 dwrr=0
allocate(dwss(kxm%n,nodes));                 dwss=0
allocate(dgs0(kxm%n,nodes));                 dgs0=0
allocate(dgr0(kxm%n,nodes));                 dgr0=0
allocate(sqa0(kxm%n,nodes));                 sqa0=0
allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes)); GGG=0

 c=(/1.,0.5,0.375,0.3125,0.2734375/)

allocate(indw(rec%w%n,2));indw=0

!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

allocate(tfldr(kxm%n,nodes));tfldr=0
call owfftwr(tfldr,nodes)

!$OMP parallel do private(inode,&
                            iwa,&
                             iz,&
                            ikx,&
                           iref,&
                            ixw,&
                              j,&
                        nrefvel,&
                             gr,&
                             gs,&
                            dgr,&
                            dgs,&
                             sqa) ordered schedule(dynamic,1) num_threads(nodes)

do iwa = 1,rec%w%n

    inode=omp_get_thread_num()+1; 
    ws(padtraces+1:sou%h%n+padtraces,inode)=source(:,iwa)
    wr(padtraces+1:rec%h%n+padtraces,inode)=recv(:,iwa)
    dws(:,inode)=0
    dwr(:,inode)=0

    do iz=nwbott,img%z%n

!--------------------------------------------------
! Compute vref
!
       nrefvel=count(mask = sref(:,1,iz).ne.-1)
      
       if (nrefvel.gt.1) then
          if (allocated(gs)) deallocate(gs); allocate(gs(kxm%n,nrefvel));
          if (allocated(gr)) deallocate(gr); allocate(gr(kxm%n,nrefvel));
          if (allocated(dgs)) deallocate(dgs); allocate(dgs(kxm%n,nrefvel));
          if (allocated(dgr)) deallocate(dgr); allocate(dgr(kxm%n,nrefvel));
          gs=0
          gr=0
          dgs=0
          dgr=0
       else
          gs0(:,inode)=0
          gr0(:,inode)=0
          dgs0(:,inode)=0
          dgr0(:,inode)=0
       end if

       tflds(:,inode)=ws(:,inode)
       call owfftexes(tflds(:,inode),inode)
       ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

       tfldr(:,inode)=wr(:,inode)
       call owfftexer(tfldr(:,inode),inode)
       wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

       ddwr(:,inode)=0;ddws(:,inode)=0

       do j=0,4
          dwss(:,inode)=kxm%all**(2*j)*ws(:,inode)
          dwrr(:,inode)=-kxm%all**(2*j)*wr(:,inode)

          tflds(:,inode)=dwss(:,inode)
          call owifftexes(tflds(:,inode),inode)
!          dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
          dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1)

          tfldr(:,inode)=dwrr(:,inode)
          call owifftexer(tfldr(:,inode),inode)
!          dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
          dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1)

          ddws(:,inode)=ddws(:,inode)+dwss(:,inode)
          ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode)
       end do
    
!       dws(:,inode)=dws(:,inode)+ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz-1,:)
!       dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz-1,:)
       dws(:,inode)=dws(:,inode)+ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz,:)
       dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz,:)

       tflds(:,inode)=dws(:,inode)
       call owfftexes(tflds(:,inode),inode)
       dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

       tfldr(:,inode)=dwr(:,inode)
       call owifftexer(tfldr(:,inode),inode)
       dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

       if (nrefvel .eq. 1) then
          gs0(:,inode)=ws(:,inode);    gr0(:,inode)=wr(:,inode);
          dgs0(:,inode)=dws(:,inode);  dgr0(:,inode)=dwr(:,inode);

          do ikx = 1,kxm%n

             sqa = (rec%w%all(iwa) * sref(1,1,iz))**2 - kxm%all(ikx)**2

             if (sqa > 0) then
                gs0(ikx,inode)=gs0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                gr0(ikx,inode)=gr0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
                dgs0(ikx,inode)=dgs0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                dgr0(ikx,inode)=dgr0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
             else
                gs0(ikx,inode)=gs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                gr0(ikx,inode)=gr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                dgs0(ikx,inode)=dgs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                dgr0(ikx,inode)=dgr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
             end if
          end do

          tflds(:,inode)=gs0(:,inode)
          call owifftexes(tflds(:,inode),inode)
          ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

          tfldr(:,inode)=gr0(:,inode)
          call owifftexer(tfldr(:,inode),inode)
          wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

          tflds(:,inode)=dgs0(:,inode)
          call owifftexes(tflds(:,inode),inode)
          dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

          tfldr(:,inode)=dgr0(:,inode)
          call owifftexer(tfldr(:,inode),inode)
          dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
 
       else

          do iref=1,nrefvel                 

             gs(:,iref)=ws(:,inode)
             gr(:,iref)=wr(:,inode) 
             dgs(:,iref)=dws(:,inode)
             dgr(:,iref)=dwr(:,inode) 

             do ikx = 1,kxm%n

                sqa = (rec%w%all(iwa) * sref(iref,1,iz))**2 - kxm%all(ikx)**2

                if (sqa > 0) then
                   gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                   gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                   dgr(ikx,iref)=dgr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                else
                   gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                   gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                   dgr(ikx,iref)=dgr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                end if
             end do

             tflds(:,inode)=gs(:,iref)
             call owifftexes(tflds(:,inode),inode)
             gs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n)

             tfldr(:,inode)=gr(:,iref)
             call owifftexer(tfldr(:,inode),inode)
             gr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n)

             tflds(:,inode)=dgs(:,iref)
             call owifftexes(tflds(:,inode),inode)
             dgs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n)

             tfldr(:,inode)=dgr(:,iref)
             call owifftexer(tfldr(:,inode),inode)
             dgr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n)

             gs(:,iref) = gs(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)));
             gr(:,iref) = gr(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)));
             dgs(:,iref) = dgs(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)));
             dgr(:,iref) = dgr(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)));

          end do

          do iref=1,nrefvel-1
             where (vel(iz,:) >= 1/sref(iref,1,iz) .and. vel(iz,:) < 1/sref(iref+1,1,iz))
                ws(:,inode)=gs(:,iref)
                wr(:,inode)=gr(:,iref)
                dws(:,inode)=dgs(:,iref)
                dwr(:,inode)=dgr(:,iref)
             end where
          end do
          where (vel(iz,:) >= 1/sref(nrefvel,1,iz))
             ws(:,inode)=gs(:,nrefvel)
             wr(:,inode)=gr(:,nrefvel)
             dws(:,inode)=dgs(:,nrefvel)
             dwr(:,inode)=dgr(:,nrefvel)
          end where 
 
       end if

       ws(:,inode)=ws(:,inode)*taper
       wr(:,inode)=wr(:,inode)*taper
       dws(:,inode)=dws(:,inode)*taper
       dwr(:,inode)=dwr(:,inode)*taper
 
       if (exist_file("ws")) then
!$OMP ORDERED
          indw(iwa,2)=inode
          GGR(iz,:,iwa)=wr(:,indw(iwa,2))
          GGS(iz,:,iwa)=ws(:,indw(iwa,2))
!$OMP END ORDERED
       end if

       call owimage_f(iwa,ws(:,inode),wr(:,inode),dws(:,inode),dwr(:,inode),GGG(iz,:,:,inode))

!========================================================================||
!========================================================================||
!    call owimage0(iwa,ws,wr,GG(iz,:,:))
! if (nwbott.ne.1) then
!    call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),vref,minv,maxv)
!    tfld=0.; tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
!    tfld=0.; tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
! !      Operations in the wavenumber domain
! !      phase shift
!    do ikx = 1,kxm%n
!       sqa = rec%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2
!       if (sqa > 0) then
!          wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
!          ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
!       else
!          wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
!          ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
!       end if
!    end do
!    tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
!    tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
!    wr=wr*taper
!  ws=ws*taper
! end if
!========================================================================
!========================================================================
	
    end do
      if (mod(iwa,20)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, &
					" FROM NODE: ",inode
end do

!$OMP END PARALLEL DO

call owfftdestroy()
deallocate(tfldr,tflds)

a=0
do i=1,nodes
   !$OMP parallel do private(ixw,iz,ish)
   do iz=1,img%z%n
      do ish=1,img%sh%n
         do ixw=1,kxm%n
            a=a+abs(GGG(iz,ish,ixw,i))
         end do
      end do
   end do
   !$OMP END PARALLEL DO
end do
write(0,*) "******   ",a,"   ******"



do i=1,nodes
   !$OMP parallel do private(ixw,iz,ish)
   do iz=1,img%z%n
      do ish=1,img%sh%n
         do ixw=1,kxm%n
            deltaI(iz,ish,ixw)=deltaI(iz,ish,ixw)+GGG(iz,ish,ixw,i)
         end do
      end do
   end do
   !$OMP END PARALLEL DO
end do
6 format(A13,i3,A8,i3,A14,i3,A12,i3)
deallocate(wr,ws,ddwr,ddws,gs0,gr0)
deallocate(dwr,dws,dwrr,dwss)
deallocate(dgs0,dgr0,sqa0,GGG)

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

complex						:: recv(:,:,:),source(:,:,:)
real						:: deltaI(:,:,:)
integer						:: iwa,iz,ikx,iref,ixw,j,inode,i,ish,ihso,nrefvel,nodes
real						:: sqa,a
real						:: taper(:),vel(:,:),sref(:,:,:),c(5),deltaS(:,:)
integer,allocatable					:: indw(:,:)
real, allocatable					:: GGG(:,:,:),sqa0(:,:)
complex, allocatable				:: ws(:,:),wr(:,:),gs0(:,:),gr0(:,:)
complex, allocatable				:: ddws(:,:),ddwr(:,:),gs(:,:),gr(:,:)
complex, allocatable				:: dws(:,:,:),dwr(:,:,:),dgs0(:,:),dgr0(:,:)
complex, allocatable				:: dwss(:,:),dwrr(:,:),dgs(:,:),dgr(:,:)
complex,allocatable					:: awr(:,:),aws(:,:),adwr(:,:),adws(:,:)
complex,allocatable					:: tflds(:,:),tfldr(:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

if (node.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
else
  nodes=node
end if
write(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"

allocate(wr(kxm%n,nodes));                   wr=0
allocate(ws(kxm%n,nodes));                   ws=0
allocate(ddwr(kxm%n,nodes));                 ddwr=0
allocate(ddws(kxm%n,nodes));                 ddws=0
allocate(gs0(kxm%n,nodes));                  gs0=0
allocate(gr0(kxm%n,nodes));                  gr0=0
allocate(dwr(img%z%n,kxm%n,nodes));          dwr=0
allocate(dws(img%z%n,kxm%n,nodes));          dws=0
allocate(adwr(kxm%n,nodes));                 adwr=0
allocate(adws(kxm%n,nodes));                 adws=0
allocate(awr(kxm%n,nodes));                  awr=0
allocate(aws(kxm%n,nodes));                  aws=0
allocate(dwrr(kxm%n,nodes));                 dwrr=0
allocate(dwss(kxm%n,nodes));                 dwss=0
allocate(dgs0(kxm%n,nodes));                 dgs0=0
allocate(dgr0(kxm%n,nodes));                 dgr0=0
allocate(sqa0(kxm%n,nodes));                 sqa0=0
allocate(GGG(img%z%n,kxm%n,nodes));          GGG=0

 c=(/1.,0.5,0.375,0.3125,0.2734375/)

!--------------------------------------------------
! Setting FFTW plans
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

allocate(tfldr(kxm%n,nodes));tfldr=0
call owfftwr(tfldr,nodes)

!$OMP parallel do private(inode,&
                            iwa,&
                             iz,&
                            ikx,&
                           iref,&
                            ixw,&
                              j,&
                        nrefvel,&
                             gr,&
                             gs,&
                            dgr,&
                            dgs,&
                             sqa) ordered schedule(dynamic,1) num_threads(nodes)
!
! Freq
!
do iwa = 1,rec%w%n

    source=conjg(source)
    inode=omp_get_thread_num()+1
    dws(:,:,inode)=0.
    dwr(:,:,inode)=0.
    dws_z(:,inode)=0.
    dws_r(:,inode)=0.
    aws(:,inode)=0.
    awr(:,inode)=0.
    ddwr(:,inode)=0
    ddws(:,inode)=0
    dwrr(:,inode)=0
    dwss(:,inode)=0

    do iz=1,img%z%n
       do ixw=1,kxm%n
          do ish=1,img%sh%n
             ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
!             if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then
!                dwr(iz,ixw,inode)=dwr(iz,ixw,inode)+source(iz,ixw-ihso,iwa)*deltaI(iz,ish,ixw)
!                dws(iz,ixw,inode)=dws(iz,ixw,inode)+recv(iz,ixw+ihso,iwa)*deltaI(iz,ish,ixw)
             if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then
                dwr(iz,ixw,inode)=dwr(iz,ixw,inode)+source(iz,ixw-ish,iwa)*deltaI(iz,ish,ixw)
                dws(iz,ixw,inode)=dws(iz,ixw,inode)+(recv(iz,ixw+ish,iwa))*deltaI(iz,ish,ixw)
             end if
          end do
       end do
    end do

!
! Depth
!
    do iz=img%z%n,nwbott+1,-1
       dws_z(:,inode)=dws_z(:,inode)+dws(iz,:,inode)
       dwr_z(:,inode)=dwr_z(:,inode)+dwr(iz,:,inode)

!--------------------------------------------------
! Compute vref
!
       nrefvel=count(mask = sref(:,1,iz).ne.-1)

       if (nrefvel.gt.1) then
          if (allocated(dgs)) deallocate(dgs); allocate(dgs(kxm%n,nrefvel));
          if (allocated(dgr)) deallocate(dgr); allocate(dgr(kxm%n,nrefvel));
          dgs=0
          dgr=0
       else
          dgs0(:,inode)=0
          dgr0(:,inode)=0
       end if

       tflds(:,inode)=dws_z(:,inode)
       call owfftexes(tflds(:,inode),inode)
       dws_z(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

       tfldr(:,inode)=dwr_z(:,inode)
       call owfftexer(tfldr(:,inode),inode)
       dwr_z(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

       if (nrefvel .eq. 1) then
          dgs0(:,inode)=dws_z(:,inode)
          dgr0(:,inode)=dwr_z(:,inode)

          do ikx = 1,kxm%n

             sqa = (rec%w%all(iwa) * sref(1,1,iz))**2 - kxm%all(ikx)**2

             if (sqa > 0) then
                dgs0(ikx,inode)=dgs0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
                dgr0(ikx,inode)=dgr0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa)))
             else
                dgs0(ikx,inode)=dgs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
                dgr0(ikx,inode)=dgr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d)
             end if
          end do

          tflds(:,inode)=dgs0(:,inode)
          call owifftexes(tflds(:,inode),inode)
          dws_z(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

          tfldr(:,inode)=dgr0(:,inode)
          call owifftexer(tfldr(:,inode),inode)
          dwr_z(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

       else

          do iref=1,nrefvel                 

             dgs(:,iref)=dws(:,inode)
             dgr(:,iref)=dwr(:,inode) 

             do ikx = 1,kxm%n

                sqa = (rec%w%all(iwa) * sref(iref,1,iz))**2 - kxm%all(ikx)**2;

                if (sqa > 0) then
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
                   dgr(ikx,iref)=dgr(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                else
                   dgs(ikx,iref)=dgs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                   dgr(ikx,iref)=dgr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                end if
             end do

             tflds(:,inode)=dgs(:,iref)
             call owifftexes(tflds(:,inode),inode)
             dgs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n)

             tfldr(:,inode)=dgr(:,iref)
             call owifftexer(tfldr(:,inode),inode)
             dgr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n)

             dgs(:,iref) = dgs(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)))
             dgr(:,iref) = dgr(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz)))

          end do

          do iref=1,nrefvel-1
             where (vel(iz,:) >= 1/sref(iref,1,iz) .and. vel(iz,:) < 1/sref(iref+1,1,iz))
                dws(:,inode)=dgs(:,iref)
                dwr(:,inode)=dgr(:,iref)
             end where
          end do
          where (vel(iz,:) >= 1/sref(nrefvel,1,iz))
             dws(:,inode)=dgs(:,nrefvel)
             dwr(:,inode)=dgr(:,nrefvel)
          end where 
 
       end if

       dws(:,inode)=dws(:,inode)*taper
       dwr(:,inode)=dwr(:,inode)*taper

       ddwr(:,inode)=0
       ddws(:,inode)=0
       aws(:,inode)=0.
       awr(:,inode)=0.

       tflds(:,inode)=source(iz-1,:,iwa)
       call owfftexes(tflds(:,inode),inode)
       aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)

       tfldr(:,inode)=recv(iz-1,:,iwa)
       call owfftexer(tfldr(:,inode),inode)
       awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)

       do j=0,4
          dwss(:,inode)=kxm%all**(2*j)*aws(:,inode)
          dwrr(:,inode)=kxm%all**(2*j)*awr(:,inode)

          tflds(:,inode)=dwss(:,inode)
          call owifftexes(tflds(:,inode),inode)
          dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)

          tfldr(:,inode)=dwrr(:,inode)
          call owifftexer(tfldr(:,inode),inode)
          dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)

          ddws(:,inode)=ddws(:,inode)+dwss(:,inode)
          ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode)
       end do

       ddws(:,inode)=ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
       ddwr(:,inode)=ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))

       aws(:,inode)=0.
       awr(:,inode)=0.

       do ixw=1,kxm%n
          do ish=1,img%sh%n
             ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
!             if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then
!                awr(ixw,inode)=aws(ixw,inode)+source(iz,ixw-ihso,iwa)*deltaI(iz,ish,ixw)
!                aws(ixw,inode)=awr(ixw,inode)+recv(iz,ixw+ihso,iwa)*deltaI(iz,ish,ixw)
             if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then
                awr(ixw,inode)=aws(ixw,inode)+source(iz,ixw-ish,iwa)*deltaI(iz,ish,ixw)
                aws(ixw,inode)=awr(ixw,inode)+(recv(iz,ixw+ish,iwa))*deltaI(iz,ish,ixw)
             end if
          end do
          dws(:,inode)=dws(:,inode)+aws(:,inode)
          dwr(:,inode)=dwr(:,inode)+awr(:,inode)
       end do
!       GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real(conjg(ddws(:,inode))*dws(:,inode)+conjg(ddwr(:,inode))*dwr(:,inode))
       GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real((ddws(:,inode))*dws(:,inode)+(ddwr(:,inode))*dwr(:,inode))
!    call owimage_b(iwa,dws(:,inode),dwr(:,inode),ddws(:,inode),ddwr(:,inode),source(img%z%n,:,iwa),recv(img%z%n,:,iwa), &
!                   deltaI(img%z%n,ish,ixw),GGG(img%z%n-1,:,inode)) 

!========================================================================||
!========================================================================||
!    call owimage0(iwa,ws,wr,GG(iz,:,:))
! if (nwbott.ne.1) then
!    call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),vref,minv,maxv)
!    tfld=0.; tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
!    tfld=0.; tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
! !      Operations in the wavenumber domain
! !      phase shift
!    do ikx = 1,kxm%n
!       sqa = rec%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2
!       if (sqa > 0) then
!          wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
!          ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
!       else
!          wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
!          ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
!       end if
!    end do
!    tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
!    tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
!    wr=wr*taper
!  ws=ws*taper
! end if
!========================================================================
!========================================================================
	
    end do
      if (mod(iwa,20)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, &
					" FROM NODE: ",inode
6 format(A13,i3,A8,i3,A14,i3,A12,i3)

end do

!$OMP END PARALLEL DO

 call owfftdestroy()
 deallocate(tfldr,tflds)

a=0
do i=1,nodes
   !$OMP parallel do private(ixw,iz,ish)
   do iz=1,img%z%n
      do ixw=1,kxm%n
         a=a+abs(GGG(iz,ixw,i))
      end do
   end do
   !$OMP END PARALLEL DO
end do
write(0,*) "******   ",a,"   ******"

do i=1,nodes
   !$OMP parallel do private(ixw,iz,ish)
   do iz=1,img%z%n
      do ixw=1,kxm%n
         deltaS(iz,ixw)=deltaS(iz,ixw)+GGG(iz,ixw,i)
      end do
   end do
   !$OMP END PARALLEL DO
end do

deallocate(wr,ws,ddwr,ddws,gs0,gr0)
deallocate(dwr,dws,adwr,adws,dwrr,dwss)
deallocate(dgs0,dgr0,sqa0,GGG)

end subroutine

!
!    tflds(:,inode)=source(img%z%n-1,:,iwa)
!    call owfftexes(tflds(:,inode),inode)
!    aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
!
!    tfldr(:,inode)=recv(img%z%n-1,:,iwa)
!    call owfftexer(tfldr(:,inode),inode)
!    awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)
!
!    do j=0,4
!       dwss(:,inode)=kxm%all**(2*j)*aws(:,inode)
!       dwrr(:,inode)=kxm%all**(2*j)*awr(:,inode)
!
!       tflds(:,inode)=dwss(:,inode)
!       call owifftexes(tflds(:,inode),inode)
!       dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
!
!       tfldr(:,inode)=dwrr(:,inode)
!       call owifftexer(tfldr(:,inode),inode)
!       dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1)
!
!       ddws(:,inode)=ddws(:,inode)+dwss(:,inode)
!       ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode)
!    end do
!
!    ddws(:,inode)=ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
! !   ddwr(:,inode)=ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
!
!
!    do ixw=1,kxm%n
!       do ish=1,img%sh%n
!          ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
!!          if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then
!!             dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ihso,iwa)*deltaI(img%z%n,ish,ixw)
!!             dws(ixw,inode)=dws(ixw,inode)+recv(img%z%n,ixw+ihso,iwa)*deltaI(img%z%n,ish,ixw)
!          if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then
!             dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ish,iwa)*deltaI(img%z%n,ish,ixw)
!             dws(ixw,inode)=dws(ixw,inode)+(recv(img%z%n,ixw+ish,iwa))*deltaI(img%z%n,ish,ixw)
!          end if
!       end do
!    end do
!    GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real(conjg(ddws(:,inode))*dws(:,inode)+conjg(ddwr(:,inode))*dwr(:,inode))
!    GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real((ddws(:,inode))*dws(:,inode)+(ddwr(:,inode))*dwr(:,inode))


end module
