program OWmig2d_areal

use sep 
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use OW_fftw_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(:),wav(:),wav0(:)
real, allocatable					:: B_image(:,:,:),illum(:,:),illumination(:,:)
real, allocatable					:: vel(:),vref(:,:,:),vel0(:)
complex, allocatable				:: recv(:,:),source(:,:),tfld(:),source1(:)
real,allocatable					:: deltaI(:,:,:)

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

 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!
allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=0.
allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
allocate(taper(kxm%n))
allocate(velint(v%z%n*v%xm%n),vel(img%z%n*kxm%n),vel0(img%z%n*img%xm%n))
allocate(illum(kxm%n,img%z%n));illum=0.
allocate(illumination(img%xm%n,img%z%n));illumination=0.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(tapertr,taper)

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

 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_vel(velint,vel0,.true.,ni)
 deallocate(velint)
 if (lloyd) then 
   allocate(vref(nref,2,img%z%n))
   call lloyd_vel1(vel0,vref)
 end if


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

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0;
          if (allocated(source)) deallocate(source); 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
       end if
       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_vel(vel0,vel,.false.,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)
          if (lloyd) then
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,illum,vref)
          else
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,illum)
          end if
          if (areal) then
             B_image=B_image+deltaI
             illumination=illumination+illum
          else
             B_image(ni:ni+kxm%n-1,:,:)=B_image(ni:ni+kxm%n-1,:,:)+deltaI
             illumination(ni:ni+kxm%n-1,:)=illumination(ni:ni+kxm%n-1,:)+illum
          end if
          if (light.and.not(timeshift)) then
             do ih=1,img%sh%n
                do iz=1,img%z%n
                   B_image(:,ih,iz)=B_image(:,ih,iz)*illumination(:,iz)/(illumination(:,iz)**2+0.5)
                end do
             end do
          end if

          call owexec_time(finish)

          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
!          call owimage_io(B_image,ixs)
       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*nsecdim*img%xm%n)
       if (light) call srite("illum",illumination,4*img%z%n*img%xm%n)
!

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

 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)
 if (allocated(illum)) deallocate(illum)
 if (allocated(illumination)) deallocate(illumination)

 call sep_end_prog()

end program OWmig2d_areal
