program OWmig3d_areal

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

implicit none

character*80                                                :: listagem
real						:: start, finish, start1, finish1, o3,a
integer						:: ixs,ino,ni,i,st,idip,ih,iz
real, allocatable					:: velint(:),taper(:,:)
real, allocatable					:: B_image(:,:,:,:,:)
real, allocatable					:: vel(:),vref(:,:,:),vel0(:)
complex, allocatable				:: recv(:,:,:),source(:,:,:)
real,allocatable					:: deltaI(:,:,:,:,:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()

 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)
 open(unit=13,file="refvel.txt")

 start=0.;finish=0.
 call owexec_time(start)
 call owparam()

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!
allocate(deltaI(kxm%n,nsecdimx,kym%n,nsecdimy,img%z%n));deltaI=0.
allocate(B_image(img%xm%n,nsecdimx,img%ym%n,nsecdimy,img%z%n));B_image=0.
allocate(taper(kxm%n,kym%n))
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))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(taper)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!

 call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
 call owgrab_velim(velint,vel0,ni)
 deallocate(velint)
 if (extrap(1:4).eq."pspi") then
    if (lloyd) then 
      allocate(vref(nref,2,img%z%n))
      call lloyd_vel1(vel0,vref)
    end if
 end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing bkg.image
!
       start1=0.;finish1=0.;
       if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n,nsecdimx,img%ym%n,nsecdimy,img%z%n));B_image=0.
       B_image=0.

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       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;deltaI=0

          call sreed("rec",recv,8*rec%hx%n*rec%w%n*rec%hy%n)
          
          if (.not.lloyd) then
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI)
          else
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,vref)
          end if
          B_image=B_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 owexec_time(finish1)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min"

       call srite("bkgimg",B_image,4*img%z%n*img%xm%n*nsecdimx*img%ym%n*nsecdimy)

 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a4)
 close(3)
 close(13)

 if (allocated(source)) deallocate(source)
 if (allocated(recv)) deallocate(recv)
 if (allocated(deltaI)) deallocate(deltaI)
 if (allocated(B_image)) deallocate(B_image)
 if (allocated(taper)) deallocate(taper)
 if (allocated(vel)) deallocate(vel)
 if (allocated(vref)) deallocate(vref)
 if (allocated(vel0)) deallocate(vel0)

 call sep_end_prog()

end program OWmig3d_areal
