program OWmod2d_par

use sep 
use OW_pspi
use OW_ffd
use OW_type
use OW_exec_time
use OW_parms
use OW_taper
use OW_grabvel

implicit none

real						:: start, finish, start1, finish1, owrk1
real, allocatable					:: velint(:,:),vel(:,:),vref(:,:,:),taper(:),illumination(:,:)
real, allocatable					:: vreflect(:,:),reflectivity(:,:,:),refl(:,:,:),illum1(:,:)
complex, allocatable				:: source(:),source1(:,:)
complex, allocatable				:: g1(:,:,:),g2(:,:,:,:),shot(:,:)
integer						:: ixs,ixr,ixm,iz,minx,maxx,nv,iw,ni1,nv1,iblock,ixxp
integer, external 					:: omp_get_thread_num,omp_get_num_threads

 call sep_init()
 call sep_begin_prog()

!---------------------------------------
! IO
!
 call owparam()
 allocate(taper(kxm%n));taper=1.
 call owtaper(taper)

!----------------------------------------------------
! Memory Allocation
!
allocate(source(sou%w%n))
allocate(vel(img%z%n,kxm%n),velint(v%z%n,v%xm%n))
!--------------------------------------------------
! Reading source fctn and velocity
 call sreed("sou",source,sou%w%n*8)

 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

!--------------------------------------------------
!
!
write(0,*) "=========================="

allocate(source1(kxm%n,sou%w%n))
!--------------------------------------------------
! CONVOLVING GFs AND REFLECTIVITY
!
allocate(reflectivity(v%z%n,noff,v%xm%n))
allocate(refl(img%z%n,noff,kxm%n))
if (exist_file("reflectivity")) then
   call sreed("reflectivity",reflectivity,4*v%z%n*v%xm%n*noff)
else
   allocate(vreflect(v%z%n,v%xm%n))
   call sreed("vreflect",vreflect,v%z%n*v%xm%n*4)
   call owreflc(vreflect,reflectivity)
   if (exist_file("refout")) then
      call srite("refout",reflectivity,v%z%n*v%xm%n*4)
   end if
   deallocate(vreflect)
end if

!--------------------------------------------------
!
!
if (light) then
   allocate(illumination(img%xm%n,img%z%n));illumination=0.
   allocate(illum1(kxm%n,img%z%n));illum1=0.
end if
   do ixs=1,rec%x%n
       start=0.;finish=0.;vel=0;refl=0.
       call owexec_time(start)
       if (marine) then
          call owgrab_vel(ixs,velint,vel,ni1)
          call owgrab_refl(ixs,reflectivity,refl,ni1)
          source1=0;source1(isp,:)=source
       else
          call owgrab_vel(1,velint,vel,ni1)
          call owgrab_refl(1,reflectivity,refl,ni1)
          ni1=1
          !--------------------------------------------------
          ! shot position in the acq.geom.
          ixxp=(tapertr+padtraces)+(ixs-1)*(rec%x%d/rec%h%d)+abs(rec%h%o)/rec%h%d+1
          source1=0;source1(ixxp,:)=source
       end if
       if (extrap(1:3).ne."ffd") then
          call owpspi_init(ixs,ni1,wbott)
          if (exist_file("vref")) then
             if (light) then
                illum1=0.
                call owpspi(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl,vref_in=vref,illum_in=illum1)
                illumination(ni1:ni1+kxm%n-1,:)=illumination(ni1:ni1+kxm%n-1,:)+illum1
             else
                call owpspi(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl,vref_in=vref)
             end if
          else
             if (light) then
                illum1=0.
                call owpspi(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl,illum_in=illum1)
                illumination(ni1:ni1+kxm%n-1,:)=illumination(ni1:ni1+kxm%n-1,:)+illum1
             else
                call owpspi(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl)
             end if
          end if
       else
          if (light) then
             illum1=0.
             call owffd(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl,illum_in=illum1)
             illumination(ni1:ni1+kxm%n-1,:)=illumination(ni1:ni1+kxm%n-1,:)+illum1
          else
             call owffd(source_in=source1,taper_in=taper,vel_in=vel,refl_in=refl)
          end if
       end if
       call owexec_time(finish)
       write(0,*) " GF # ",ixs," OUT OF ",rec%x%n," in ", (finish - start)," s",ni1
   end do

if (light) then
   call srite("illum",illumination,4*img%xm%n*img%z%n)
   deallocate(illum1,illumination)
end if
deallocate(taper,vel,velint,source1,source)

call sep_end_prog()
end program OWmod2d_par
