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

implicit none

real						:: start, finish, start1, finish1, o3,a
integer						:: ixs,ino,ni,i,st,idip,ih
real, allocatable					:: velint(:),taper(:),wav(:),wav0(:)
real, allocatable					:: R_image(:,:,:)
real, allocatable					:: vel(:),vref(:,:,:),vel0(:)
complex, allocatable				:: recv(:,:),source(:,:),tfld(:),source1(:)
real,allocatable					:: deltaI(:,:,:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 open(unit=2,file='lixo')

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

 allocate(source(sou%h%n,sou%w%n))
 allocate(recv(rec%h%n,rec%w%n))
 allocate(R_image(nz,img%sh%n,img%xm%n));R_image=0.
 allocate(deltaI(nz,img%sh%n,img%xm%n));deltaI=0.
 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))
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 if (exist_file("vref")) then
   allocate(vref(nref,2,v%z%n))
   call sreed("vref",vref,v%z%n*2*nref*4)
 end if
 call owgrab_vel(velint,vel0,.true.)

 deallocate(velint)
 call owexec_time(finish)
 write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish - start),' (sec)'

! Computing bkg.image
       start1=0.;finish1=0.;
       call sseek("in",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
          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("in",recv,8*rec%h%n*rec%w%n)
          if (exist_file("vref")) then
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,vref)
          else
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI)
          end if
          if (areal) then
             R_image=R_image+deltaI
          else
             R_image(:,:,ni:ni+kxm%n-1)=R_image(:,:,ni:ni+kxm%n-1)+deltaI
          end if
          call owexec_time(finish)
          write(0,111) "EXTRAPOLATED WAVEFIELDS = ",idip," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       if (exist_file("wav")) then
          allocate(wav(3*nz),wav0(3*nz),tfld(3*nz))
          wav=0;wav0=0;tfld=0.
          call sreed("wav",wav(1:nz),4*nz)
          call owfftw(tfld)
          tfld(1:nz)=cmplx(wav(1:nz),0.)
          call sfftw_execute(planF)
          wav=abs(tfld/sqrt(3.*nz))
          wav=wav/maxval(wav)
          tfld=0.
          tfld(1:nz)=cmplx(R_image(:,int(abs(img%sh%o/img%sh%d)+1.),int(abs(img%xm%o/img%xm%d)+1.)),0.)
          call sfftw_execute(planF)
          wav0=abs(tfld/sqrt(3.*nz))
          wav0=wav0/maxval(wav0)
          wav=wav/(wav0+0.001)
          wav0=wav
          do i=4,3*nz-3
             wav(i)=sum(wav0(i-3:i+3))/7
          end do
          do ixs=1,img%xm%n
             do ih=1,img%sh%n
                tfld=0;
                tfld(1:nz)=cmplx(R_image(:,ih,ixs),0.)
                call sfftw_execute(planF)
                tfld=tfld/sqrt(3.*nz)
                tfld=tfld*wav
                call sfftw_execute(planI)
                R_image(:,ih,ixs)=real(tfld(1:nz)/sqrt(3.*nz))
             end do
          end do
          deallocate(wav,wav0,tfld)
       end if
!       call srite("out",R_image(:,:,tapertr+padtraces+1:img%xm%n-tapertr+padtraces), &
!                    4*nz*img%sh%n*(img%xm%n-2*(tapertr+padtraces)))
       call srite("out",R_image,4*nz*img%sh%n*img%xm%n)

 111 format(a26,i4,a8,i4,a4,f6.2,a4)

 close(2)

 deallocate(source)
 deallocate(recv)
 deallocate(R_image,deltaI)
 deallocate(taper)
 deallocate(vel)
 if (exist_file("vref")) then
    deallocate(vref)
 end if
 call sep_end_prog()

end program OWmig2d_areal
