!!!!!!!!!!!!!!!!!!!!!!!
!

program OWwemva_opers

use sep 
use OW_type
use OW_parms_areal
use OW_grabvel2_areal
use OW_exec_time
use OW_taper_areal
use OW_wei_areal

implicit none

character*80                            		:: hor1
logical						:: dotimg,pspi,dotdso
real						:: start, finish,start1,finish1, o3,a
integer						:: ixs,ino,ni,i,st,idip,n3,t3,iz,k,j,stat
real, allocatable					:: velint(:),taper(:),deltaSint(:)
real, allocatable					:: R_image(:,:,:),scatS(:)
real, allocatable					:: vel(:),vref(:,:,:),vel0(:)
complex, allocatable				:: recv(:,:),source(:,:),source1(:)
real,allocatable					:: deltaI(:,:,:),deltaS(:),deltaS0(:)

 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()
 if (adj) then
    allocate(deltaS(img%z%n*kxm%n));deltaS=0.
    allocate(scatS(img%z%n*img%xm%n));scatS=0.
    allocate(deltaI(img%xm%n,nsecdim,img%z%n));deltaI=0.
    call sreed("deltaI",deltaI,4*img%z%n*nsecdim*img%xm%n)
 else
    allocate(deltaSint(v%z%n*v%xm%n),deltaS(img%z%n*kxm%n));deltaSint=0;deltaS=0
    allocate(deltaS0(img%z%n*img%xm%n));deltaS0=0
    allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=0
    allocate(R_image(img%xm%n,nsecdim,img%z%n));R_image=0.
    call sreed("deltaS",deltaSint,v%z%n*v%xm%n*4)
    call owgrab_velim(deltaSint,deltaS0)
    deallocate(deltaSint)
 end if

 allocate(taper(img%xm%n))
 call owtaper(tapertr,taper)

 allocate(velint(v%z%n*v%xm%n),vel(img%z%n*kxm%n),vel0(img%z%n*img%xm%n))
 allocate(vref(nref,2,v%z%n))
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_velim(velint,vel0)
 deallocate(velint)
 if(exist_file("vref")) then
   call sreed("vref",vref,v%z%n*2*nref*4)
 end if
 call owexec_time(finish)
 write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish - start),' (sec)'

!----------------------------------------------

    start1=0.;finish1=0.;
    call sseek("rec",0,0)
    call sseek("sou",0,0)
    if (not(areal)) then
       allocate(source1(sou%w%n));source1=0
       allocate(source(rec%h%n,sou%w%n));source=0
       call sreed("sou",source1,8*sou%w%n)
       source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       deallocate(source1)
    end if
    call owexec_time(start1)
    if (adj) then
       scatS=0.
       do idip=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n));source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=vel0
             ni=1
          else
             call owgrab_velsht(vel0,vel,ni,idip)
          end if
          if (allocated(recv)) deallocate(recv); allocate(recv(rec%h%n,rec%w%n));recv=0;deltaS=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)
          call owwemva_init(idip)
          if (exist_file("vref")) then
             call owwemva(.true.,taper,recv,source,vel,deltaS,deltaI(ni:ni+kxm%n-1,:,:),vref)
          else
             call owwemva(.true.,taper,recv,source,vel,deltaS,deltaI(ni:ni+kxm%n-1,:,:))
          end if
          do j=ni,ni+kxm%n-1;do k=1,img%z%n
             scatS((k-1)*img%xm%n+j)=scatS((k-1)*img%xm%n+j)+deltaS((k-1)*kxm%n+j-ni+1)
          end do;end do
          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
!       scatS=scatS/maxval(abs(scatS))
       call srite("out",scatS,4*img%z%n*img%xm%n)
       deallocate(deltaI,scatS,recv,source,vel,vref,deltaS)
   else
       R_image=0
       do idip=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n));source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=vel0;deltaS=deltaS0
             ni=1
          else
             call owgrab_velsht(vel0,vel,ni,idip)
             call owgrab_velsht(deltaS0,deltaS,ni,idip)
          end if
          if (allocated(recv)) deallocate(recv); allocate(recv(rec%h%n,rec%w%n));recv=0;deltaI=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)
          call owwemva_init(idip)
          if (exist_file("vref")) then
             call owwemva(.false.,taper,recv,source,vel,deltaS,deltaI,vref)
          else
             call owwemva(.false.,taper,recv,source,vel,deltaS,deltaI)
          end if
          R_image(ni:ni+kxm%n-1,:,:)=R_image(ni:ni+kxm%n-1,:,:)+deltaI
          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       call srite("out",R_image,4*img%z%n*nsecdim*img%xm%n)
       deallocate(deltaI,R_image,recv,source,vel,vref,deltaS)
   end if
   111 format(a26,i4,a8,i4,a4,f6.2,a4)
    close(2)
 call sep_end_prog()
end program OWwemva_opers
       
