module OW_arealsh_mod

use OW_ctrill
use OW_type
use OW_grabvel_areal
use OW_parms_areal
!use OW_datum
use OW_taper_areal

implicit none

contains 

subroutine owarealsh_init(nsurf_in,dang_in,angmin_in,nang_in

subroutine owareal_sou(source,velint,)

integer						:: node,i,nv,ixs,j,ni,nsy,iw,ih
real						:: owrk
real, allocatable				:: velint(:,:),taper(:),vel(:,:)
real,allocatable				:: reflect(:,:)!,surf(:,:),reflect(:,:,:)
complex, allocatable				:: recv(:,:),source(:,:,:),synshot(:,:)
complex,allocatable				:: G1(:,:),refl(:,:),tfld(:),G2(:,:)

integer, external 				:: omp_get_thread_num,omp_get_num_threads,omp_get_num_procs,omp_set_num_threads


!----------------------------------------------------
! Memory Allocation
!
!--------------------------------------------------
!
!
write(0,*) "   COMPUTES THE SYNTHEIS OPERATOR BY UPWARD   "
write(0,*) "CONTINUING A WAVEFIELD WITH A PREDEFINED SHAPE"
write(0,*) "=============================================="

!---------------------------------------------------------------------------
!
! Predefine the wavefield
!

! For more than 1 surface
!open(unit=3,file="surfs.txt")
!allocate(surf(nsurf,img%xm%n))
!read(3,*) ((surf(i1,i2), i2=1,img%xm%n), i1=1,nsurf)
!close(3)
!allocate(reflect(nsurf,img%z%n,img%xm%n));reflect=0
!surf=(surf-img%z%o)/img%z%d
!do i1=1,nsurf
!   do i2=1,img%xm%n
!      reflect(i1,floor(surf((i1,i2)))+1,i2)=1
!   end do
!end do

!open(unit=3,file="surfs.txt")
!allocate(surf(img%xm%n))
!read(3,*) (surf(i1), i1=1,img%xm%n)
!close(3)

allocate(reflect(nang,img%z%n,img%xm%n));reflect=0
dzdx=img%z%d/img%xm%d
do i1=1,nang
   tang=tand((i1-1)*dang+angmin)
   nzx=floor(tang/dzdx)
   do i2=1,img%xm%n
      indz=deprefl/img%z%d+1-(i2-floor(img%xm%n/2.)-1)*tang/dzdx+1.5
      reflect(i1,indz,i2)=1.
      write(2,*) i1,indz,i2
   end do
end do

call owtaper(tapertr,taper)

!--------------------------------------------------
!
! Generating the Synthesis Operator
!
!-----------------------------------------------------------------------------------------------      

allocate(G1(rec%w%n,img%xm%n)); G1=0
call owctrl_init(wbott,node)
call owctrl(taper,reflect,source,vsal,vel,G1)
allocate(G2(rec%w%n,img%xm%n));G2=G1
deallocate(G1);allocate(G1(rec%w%n,img%xm%n-2*padtraces)) 
G1=G2(:,padtraces+1:img%xm%n-padtraces)
deallocate(G2)

allocate(synshot(rec%w%n,img%xm%n))

call owarealexe(G1,synshot)

deallocate(source,vel,G1)
close(2)
call sep_end_prog()
end program OWareal
