program OWmig2d_areal

use sep 
!use OW_pspi_areal
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal

implicit none

real						:: start, finish, o3
integer						:: ixs,ino,ni,i,st,idip,n3,t3,iz
real, allocatable					:: velint(:,:),taper(:),deltaSint(:,:)
real, allocatable					:: vel(:,:),R_image(:,:,:),deltaS(:,:),vref(:,:,:),scatS(:,:)
complex, allocatable				:: recv(:,:),source(:,:),recv1(:,:,:),source1(:,:,:)
complex, allocatable				:: GGr(:,:,:),GGs(:,:,:)
real,allocatable					:: deltaI(:,:,:)
real,allocatable					:: G1(:,:,:)

 call sep_init()
 call sep_begin_prog()

 start=0.;finish=0.
 call owexec_time(start)
 call owparam(.true.)
 if (adj) then
    allocate(source1(img%z%n,img%xm%n,rec%w%n))
    allocate(recv1(img%z%n,img%xm%n,rec%w%n))
    allocate(deltaS(img%z%n,img%xm%n));deltaS=0.
    allocate(scatS(img%z%n,img%xm%n));scatS=0.
    allocate(deltaI(img%z%n,img%sh%n,img%xm%n));deltaI=0.
    call sreed("deltaI",deltaI(:,:,tapertr+padtraces+1:img%xm%n-(tapertr+padtraces)), &
                                 4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
 else
    allocate(source(sou%h%n,sou%w%n))
    allocate(recv(rec%h%n,rec%w%n))
    allocate(deltaSint(v%z%n,v%xm%n),deltaS(img%z%n,img%xm%n))
    allocate(deltaI(img%z%n,img%sh%n,img%xm%n))
    allocate(R_image(img%z%n,img%sh%n,img%xm%n));R_image=0.
    call sreed("deltaS",deltaSint,v%z%n*v%xm%n*4)
    call owgrab_vel(deltaSint,deltaS,ni)
    deallocate(deltaSint)
    if (exist_file("ws")) then
        allocate(GGr(img%z%n,img%xm%n,rec%w%n))
        allocate(GGs(img%z%n,img%xm%n,rec%w%n))
    end if     
 end if
 allocate(taper(img%xm%n))
 call owtaper(tapertr,taper)
 allocate(velint(v%z%n,v%xm%n),vel(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 sreed("vref",vref,v%z%n*2*nref*4)
 call owgrab_vel(velint,vel,ni)
 deallocate(velint)
 call owexec_time(finish)
 write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish - start),' (sec)'

 open(unit=2,file='lixo')

 if (adj) then
    do idip=1,rec%x%n
       start=0.;finish=0. 
       call owexec_time(start)
       recv1=0;source1=0;deltaS=0.
       call sreed("sou",source1(:,tapertr+padtraces+1:kxm%n-tapertr-padtraces,:), &
                              8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
       call sreed("in",recv1(:,tapertr+padtraces+1:kxm%n-tapertr-padtraces,:), &
                              8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
       call owwei_b(taper,recv1,source1,vref,vel,deltaS,deltaI)
       scatS=scatS+deltaS
       call owexec_time(finish)
       write(0,*) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%x%n," IN ",(finish-start)/60,"s"
!       if (exist_file("ws")) then
!          call srite("ws",source1,8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
!       end if
!       if (exist_file("wr")) then
!          call srite("wr",recv1,8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
!       end if
    end do
else
    do idip=1,rec%x%n
       start=0.;finish=0. 
       call owexec_time(start)
       recv=0;source=0;deltaI=0
       call sreed("sou",source,8*sou%h%n*sou%w%n)
       call sreed("in",recv,8*rec%h%n*rec%w%n)
       if (exist_file("ws")) then
          GGr=0.;GGs=0.
          call owwei_f(taper,recv,source,vref,vel,deltaS,deltaI,GGr,GGs)
          call srite("ws",GGs(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,:), &
                     8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
          call srite("wr",GGr(:,tapertr+padtraces+1:img%xm%n-tapertr-padtraces,:), &
                     8*img%z%n*rec%w%n*(img%xm%n-2*(tapertr+padtraces)))
       else
          call owwei_f(taper,recv,source,vref,vel,deltaS,deltaI)
       end if
       R_image=R_image+deltaI
       call owexec_time(finish)
       write(0,*) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%x%n," IN ",(finish-start)/60,"s"
    end do
end if

 if (adj) then
    call srite("out",scatS(:,tapertr+padtraces+1:img%xm%n-tapertr+padtraces), &
               4*img%z%n*(img%xm%n-2*(tapertr+padtraces)))
    deallocate(deltaI,scatS,recv1,source1,vel,vref,deltaS)
 else
    call srite("out",R_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr+padtraces), &
               4*img%z%n*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
    deallocate(deltaI,R_image,recv,source,vel,vref,deltaS)
 end if
 close(2)
 call sep_end_prog()
end program OWmig2d_areal
