program OW_arealsht

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


implicit none

real						:: start, finish
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

 call sep_init()
 call sep_begin_prog()


node=omp_get_num_threads()

write(0,*) 'node = ',node

start=0.;finish=0.
call owexec_time(start)
!---------------------------------------
! IO
!
 call owparam(.false.)

!----------------------------------------------------
! Memory Allocation
!
allocate(source(sou%h%n,sou%w%n,sou%x%n))
allocate(velint(v%z%n,v%xm%n))
allocate(vel(img%z%n,img%xm%n))
allocate(taper(img%xm%n)) 

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

call sreed("sou",source,sou%h%n*sou%w%n*sou%x%n*8)

!---------------------------------------------------
!
! Read and extend velocities
!
call sreed("vel",velint,v%z%n*v%xm%n*4)
nv = floor((img%xm%o-v%xm%o)/v%xm%d)+1
vel(:,:)=velint(1:img%z%n,nv+1:nv+img%xm%n)

if (exist_file("velout")) then
   call to_history("n1",img%z%n,"velout");call to_history("o1",v%z%o,"velout");call to_history("d1",v%z%d,"velout"); 
   call to_history("n2",img%xm%n,"velout");call to_history("o2",img%xm%o,"velout");call to_history("d2",img%xm%d,"velout");
   call srite("velout",vel,4*v%z%n*img%xm%n)
end if
deallocate(velint)

!--------------------------------------------------
!
!
write(0,*) "UPWARD CONTINUES THE PREDEFINED WAVEFIELD"
write(0,*) "========================================="

open(unit=2,file='lixo')

!---------------------------------------------------------------------------
!
! 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)

write(0,*) "2 ",dref
allocate(reflect(img%z%n,img%xm%n));reflect=0
reflect(dref,:)=10;

if(exist_file("refkx")) then
   allocate(tfld(img%xm%n),refl(img%z%n,img%xm%n))
   call owfftw(tfld)
   do i=1,img%z%n
      tfld=cmplx(reflect(i,:),0); call sfftw_execute(planF); refl(i,:)=tfld/sqrt(1.*img%xm%n)
   end do
   call to_history("n1",img%z%n,"refkx") ; call to_history("d1",img%z%d,"refkx") ; call to_history("o1",img%z%o,"refkx")
   call to_history("n2",1,"refkx"); call to_history("d2",1,"refkx"); call to_history("o2",1,"refkx")
   call to_history("n3",img%xm%n,"refkx"); call to_history("d3",kxm%d,"refkx"); call to_history("o3",kxm%all(kxm%n),"refkx")  
   call to_history("n4",1,"refkx"); call to_history("d4",1,"refkx"); call to_history("o4",1,"refkx")
   call to_history("n5",1,"refkx"); call to_history("d5",1,"refkx"); call to_history("o5",1,"refkx")  
   call to_history("esize",8,"refkx")
   call srite("refkx",refl,8*img%z%n*kxm%n)
   deallocate(tfld,refl)
end if

!surf=(surf-img%z%o)/img%z%d
!do i2=1,img%xm%n
!   reflect(floor(surf((i2)))+1,i2)=1.
!end do

!deallocate(surf)

call owtaper(tapertr,taper)

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

if (exist_file("synth_op")) then
    allocate(G1(rec%w%n,imgxn))
    call sreed("synth_op",G1,imgxn*rec%w%n*8)
else
    allocate(G1(rec%w%n,img%xm%n)); G1=0
    start=0.;finish=0. 
    call owexec_time(start)
    call owctrl_init(wbott,node)
    call owctrl(taper,reflect,source,vsal,vel,G1)
    call owexec_time(finish)
!    allocate(G2(rec%w%n,img%xm%n));G2=G1
!    write(0,*) img%xm%n,imgxn,imgxn+padtraces
!    G1=G2(:,padtraces+1:imgxn+padtraces)
!    call srite("out",G1,8*rec%w%n*imgxn)
    call srite("out",G1,8*rec%w%n*img%xm%n)
!    deallocate(G1);allocate(G1(rec%w%n,img%xm%n))
!    G1=G2;deallocate(G2)
end if

!---------------------------------------------------
!
! Synthesizing areal shots
!
!---------------------------------------------------
if(exist_file("shot_in")) then
write(0,*) "Synthesizing areal shots"
write(0,*) "========================"

allocate(recv(rec%h%n,rec%w%n),synshot(rec%w%n,nsh_img))

recv=0;synshot=0
write(0,*) "SHOT   SHOT-X     OWRK  FIndI LIinI  G1ind  S2SYN  img_x_n"
    do ixs=1,rec%x%n
       owrk=(ixs-1)*rec%x%d+rec%x%o
       ni = (ixs-1)*floor(rec%x%d/img%xm%d) + 1
       nsy=floor(((ixs-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)+1
       call sreed("shot_in",recv,rec%w%n*rec%h%n*8)
       do ih = 1,rec%h%n
          do iw=1,rec%w%n
             synshot(iw,ni+ih-1)=synshot(iw,ni+ih-1)+recv(ih,iw)*(G1(iw,nsy))
          end do
       end do
    end do
!    do ixs=1,nsh_img
!       synshot(:,ixs)=synshot(:,ixs)*cmplx(0,2*rec%w%all(iwa)*rec%w%n))
!    end do
    call srite("synsh",synshot,8*rec%w%n*nsh_img)
end if
!----------------------------------------------------
!
! Migrating the areal shots
!
!----------------------------------------------------

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