!!!!!!!!!!!!!!!!!!!!!!!
!
! need to be adapted to shotprofile
! and timeshift imaging condition
!
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
use OW_dotst
use base_spline_mod

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(:),pick1(:)
real, allocatable					:: vel(:),vref(:,:,:),vel0(:)
complex, allocatable				:: recv(:,:,:),source(:,:,:),source1(:)
real,allocatable					:: deltaI(:,:,:,:,:),deltaS(:),deltaS0(:)
double precision,allocatable				:: gspl1(:),gspl2(:)

 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()
 if (not(dot)) then
    if (adj) then
       allocate(deltaS(img%z%n*kxm%n*kym%n));deltaS=0.
       allocate(scatS(img%z%n*img%xm%n*img%ym%n));scatS=0.
       allocate(deltaI(img%xm%n,nsecdimx,img%ym%n,nsecdimy,img%z%n));deltaI=0.
       call sreed("deltaI",deltaI,4*img%z%n*nsecdimx*img%xm%n*nsecdimy*img%ym%n)
       call from_param("oxk", axko, 0.); call from_param("dxk", axkd, 1.); call from_param("nxk", axkn, 99999 )
       call from_param("oyk", ayko, 0.); call from_param("dyk", aykd, 1.); call from_param("nyk", aykn, 99999 )
       call from_param("ozk", azko, 0.); call from_param("dzk", azkd, 1.); call from_param("nzk", azkn, 99999 )
    else
       allocate(deltaSint(v%z%n*v%xm%n*v%ym%n),deltaS(img%z%n*kxm%n*kym%n));deltaSint=0;deltaS=0
       allocate(deltaS0(img%z%n*img%xm%n*img%ym%n));deltaS0=0
       allocate(deltaI(kxm%n,nsecdimx,kym%n,nsecdimy,img%z%n));deltaI=0
       allocate(R_image(img%xm%n,nsecdimx,img%ym%n,nsecdimy,img%z%n));R_image=0.
       call sreed("deltaS",deltaSint,v%z%n*v%xm%n*v%ym%n*4)
       call owgrab_velim(deltaSint,deltaS0,ni)
       deallocate(deltaSint)
    end if
 end if
 allocate(taper(kxm%n,kym%n))
 call owtaper(taper)

 allocate(velint(v%z%n*v%xm%n*v%ym%n),vel(img%z%n*kxm%n*kym%n),vel0(img%z%n*img%xm%n*img%ym%n))
 allocate(vref(nref,2,v%z%n))

 call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
 call owgrab_velim(velint,vel0,ni)
 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)'

!----------------------------------------------
! Horizon with no velocity update
!if (exist_file("horv0")) then
!   call from_param('horv0',hor1)
!   open(51,file=hor1,status='old')
!   allocate(pick1(img%xm%n));pick1=9999.
!   do i=1,img%xm%n
!      read(51,*,end=1) pick1(i)
!   end do
!1 close(51)
!end if
 open(unit=2,file='lixo')
 if (dot) then
    recv=0;source=0;
    if (allocated(source)) deallocate(source); allocate(source(sou%hx%n,sou%hy%n,sou%w%n))
    source=0;
    call sreed("sou",source,8*sou%hx%n*sou%w%n*sou%hy%n)
    vel=vel0
    ni=1
    if (allocated(recv))   deallocate(recv);   allocate(recv(rec%hx%n,rec%hy%n,rec%w%n))
    call sreed("rec",recv,8*rec%hx%n*rec%w%n*rec%hy%n)
    if(exist_file("vref")) then
       call owdotst_wemva(taper,recv,source,vel,vref)
    else
       call owdotst_wemva(taper,recv,source,vel)
    end if
    call owdotst_image()
 else
    start1=0.;finish1=0.;
    call sseek("rec",0,0)
    call sseek("sou",0,0)
    call owexec_time(start1)
    if (adj) then
       scatS=0.
       do idip=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (allocated(source)) deallocate(source); allocate(source(sou%hx%n,sou%hy%n,sou%w%n))
          source=0;
          call sreed("sou",source,8*sou%hx%n*sou%hy%n*sou%w%n)
          vel=vel0
          ni=1
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%hx%n,rec%hy%n,rec%w%n))
          recv=0;deltaS=0
          call sreed("rec",recv,8*rec%hx%n*rec%hy%n*rec%w%n)
          call owwemva_init(idip)
          if (exist_file("vref")) then
             call owwemva(.true.,taper,recv,source,vel,deltaS,deltaI,vref)
          else
             call owwemva(.true.,taper,recv,source,vel,deltaS,deltaI)
          end if
          scatS=scatS+deltaS
          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(scatS)
       if (axkn.ne.99999 .or. azkn.ne.99999) then
          axk%d=axkd
          ayk%d=aykd
          azk%d=azkd
          axk%n=axkn
          ayk%n=aykn
          azk%n=azkn
          axk%o=axko
          ayk%o=ayko
          azk%o=azko
          stat = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
          allocate(gspl1(axkn*aykn*azkn),gspl2(size(scatS)))
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          stat = base_spln3d_oper(.true. , .false., gspl1, dble(scatS))
          stat = base_spln3d_oper(.false. ,.false., gspl1, gspl2)
          scatS=sngl(gspl2)
          deallocate(gspl1,gspl2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if
!       if (exist_file("horv0")) then
!          allocate(gspl2(size(scatS)))    
!          gspl2=dble(scatS)
!          call oweditvel(pick1,gspl2)
!          scatS=sngl(gspl2)
!          deallocate(gspl2)    
!       end if
       call srite("out",scatS,4*img%z%n*img%xm%n*img%ym%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 (allocated(source)) deallocate(source); allocate(source(sou%hx%n,sou%hy%n,sou%w%n))
          source=0;
          call sreed("sou",source,8*sou%hx%n*sou%hy%n*sou%w%n)
          vel=vel0
          ni=1
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%hx%n,rec%hy%n,rec%w%n))
          recv=0;deltaS=0
          call sreed("rec",recv,8*rec%hx%n*rec%hy%n*rec%w%n)
          call owwemva_init(idip)
          deltaI=0
          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=R_image+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*nsecdimx*img%xm%n*nsecdimy*img%ym%n)
       deallocate(deltaI,R_image,recv,source,vel,vref,deltaS)
   end if
   111 format(a26,i4,a8,i4,a4,f6.2,a4)
    close(2)
 end if
 call sep_end_prog()
end program OWwemva_opers
