program OW_arealsht

use sep 
use OW_ctrillwei
use OW_type
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_fftw_areal
use OW_picksh
use OW_rndseed1

use OW_encoding

implicit none

real						:: start, finish,start1, finish1,a(2),b,c
real						:: rw(2)=0.,rz(2,2)=0.
integer						:: i,nv,ixs,j,ni,nsy,iw,ih,icomb,ix,stat,iz,ihor,n_enc_sog,iq,nq
integer						:: ind_trac,i_pos,i_ini,i_fin,igap
integer						:: jgap,jxs,nnn,hij,k,ib,ic,ireal,igap1
integer,allocatable					:: window(:,:),HD(:,:),q(:),prefer(:,:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:),rectmp(:,:),soutmp(:,:)
real,allocatable					:: reccig(:,:),soucig(:,:),soucig_rot(:,:),source1(:),temp(:)
real,allocatable					:: tempr(:,:),temps(:,:),randz(:,:,:)
real, allocatable					:: randw(:,:),rand4(:,:,:,:),rand1(:),vref(:,:,:),randx(:)
real,allocatable					:: hold(:,:),gol(:,:,:),randq(:,:),coscode(:,:),sincode(:,:)
complex,allocatable					:: source(:),GGR(:,:),GGS(:,:)
complex, allocatable				:: GGG(:,:,:),wave_R(:,:,:),wave_S(:,:,:),tfld(:)

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

open(unit=2,file="logfile.txt")
open(unit=3,file="logfile1.txt")
open(unit=4,file="logfile2.txt")
open(unit=5,file="logfile3.txt")

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

!---------------------------------------
! IO
!
 call owparam(.false.)

!--------------------------------------------------
! Encoding functions
!
start=0.;finish=0. 
call owexec_time(start)

if (encode) then
   allocate(randz(sou%w%n,sog%xm%n,nreal))
   if(rand) then
      call init_random_seed3(randz)
   else if (gold.and.not(perfect)) then
      call owgoldw (randz)
   else if (gold.and.perfect) then
      call owgoldpw (randz)
   else if (exist_file("InCode")) then
      call sreed("InCode",randz,4*sou%w%n*sog%xm%n*nreal)
   end if
   if(exist_file("OutCode")) then
     call srite("OutCode",randz,4*sou%w%n*sog%xm%n*nreal)
   end if
   call owexec_time(finish)
   write(0,*) "Finished generating encode functions in ",finish-start,"(s)"
end if

!---------------------------------------------------
!
! Read velocities
!
allocate(velint(v%z%n,v%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

nv=int((xmin-v%xm%o)/v%xm%d)+1.5
allocate(vel(sog%z%n,nx))
vel(1:sog%z%n,:)=velint(1:sog%z%n,nv:nv+nx-1)
deallocate(velint)

!--------------------------------------------------
! Picks
!
if (exist_file("hor1")) then
   call owpick()
end if

!--------------------------------------------------
! Computing taper
!
allocate(taper(kxm%n))
call owtaper(int(padtraces/4),taper)

!--------------------------------------------------
!
write(0,*) "GENERATING AREAL DATA"
write(0,*) "====================="


!--------------------------------------------------
!
! Read rotated CIGs
!

allocate(reccig(sog%z%n,nx),soucig(sog%z%n,nx));reccig=0;soucig=0
allocate(wave_S(sou%w%n,nx,ngap*nreal),wave_R(sou%w%n,nx,ngap*nreal));wave_S=0;wave_R=0
allocate(GGS(sou%w%n,kxm%n),GGR(sou%w%n,kxm%n));GGS=0;GGR=0

do ixs=1,sog%xm%n

   allocate(soucig_rot(sog%z%n,sog%sh%n));soucig_rot=0.
   allocate(rectmp(sog%z%n,sog%sh%n));rectmp=0.
   allocate(soutmp(sog%z%n,sog%sh%n));soutmp=0.
   allocate(tempr(sog%z%n,sog%sh%n));tempr=0.
   allocate(temps(sog%z%n,sog%sh%n));temps=0.
   
   reccig=0.;soucig=0.
   start=0.;finish=0. 
   GGS=0.;GGR=0.

!--------------------------------------------------
! Selection scheme: cigs separated by ngap
   igap=ixs-int((ixs-1)/ngap)*ngap
   jgap=floor((ixs-1)/ngap*1.)+1

   jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
   call sreed("rec_cig",rectmp,4*sog%z%n*sog%sh%n)
   call sreed("sou_cig",soutmp,4*sog%z%n*sog%sh%n)

   do ih=1,sog%sh%n
      soucig_rot(:,ih)=soutmp(:,sog%sh%n-ih+1)
   end do

   deallocate(soutmp)

!--------------------------------------------------
!
! Upward continuing CIGs
!
!--------------------------------------------------     

   if (nhorz .ne. 0) then
!      allocate(window(sog%xm%n,2*nhorz))
      do ihor=1,nhorz
!         window(ixs,2*ihor-1)=minval(picks(ihor,ixs,:))
!         window(ixs,2*ihor)=maxval(picks(ihor,ixs,:))
!         np=minval(pick(ihor,:), mask=pick(ihor,:).ne.0)
!         if (np
         do iz=1,sog%z%n
            do ih=1,sog%sh%n
!               if (iz .ge. picks(ihor,ixs,ih)-wavz .and. iz .le. picks(ihor,ixs,ih)+wavz) then
               if (iz .ge. picks(ihor,ixs)-wavz .and. iz .le. picks(ihor,ixs)+wavz) then
                  tempr(iz,ih)=rectmp(iz,ih)
                  temps(iz,ih)=soucig_rot(iz,ih)
               end if
            end do
         end do
      end do
      reccig(:,jxs:jxs+sog%sh%n-1)=tempr
      soucig(:,jxs:jxs+sog%sh%n-1)=temps
   else
      reccig(:,jxs:jxs+sog%sh%n-1)=rectmp
      soucig(:,jxs:jxs+sog%sh%n-1)=soucig_rot
   end if


   deallocate(rectmp,soucig_rot,tempr,temps)

   if(exist_file("soucig")) then
      call srite("soucig",soucig,4*sog%z%n*nx)
   end if

   call owctrlwei_init(ixs)

   do ireal=1,nreal

      if (exist_file("vref")) then
         call owexec_time(start)
         call owctrlwei(taper,vsal,vel,reccig,GGR,+1,randz(:,ixs,ireal),vref)
         call owexec_time(finish)
      write(0,3) "COMPUTED REC.WAVEFIELD # ",ixs," OUT OF ",sog%xm%n,' / IGAP ',igap,' realization # ',ireal,' FROM ',nreal,' in ',finish-start,' s'
         call owexec_time(start)
         call owctrlwei(taper,vsal,vel,soucig,GGS,-1,randz(:,ixs,ireal),vref)
         call owexec_time(finish)
      write(0,3) "COMPUTED SOU.WAVEFIELD # ",ixs," OUT OF ",sog%xm%n,' / IGAP ',igap,' realization # ',ireal,' FROM ',nreal,' in ',finish-start,' s'
      else
         call owexec_time(start)
         call owctrlwei(taper,vsal,vel,reccig,GGR,+1,randz(:,ixs,ireal))
         call owexec_time(finish)
      write(0,3) "COMPUTED REC.WAVEFIELD # ",ixs," OUT OF ",sog%xm%n,' / IGAP ',igap,' realization # ',ireal,' FROM ',nreal,' in ',finish-start,' s'
         call owexec_time(start)
         call owctrlwei(taper,vsal,vel,soucig,GGS,-1,randz(:,ixs,ireal))
         call owexec_time(finish)
      write(0,3) "COMPUTED SOU.WAVEFIELD # ",ixs," OUT OF ",sog%xm%n,' / IGAP ',igap,' realization # ',ireal,' FROM ',nreal,' in ',finish-start,' s'
      end if

      igap1=(ireal-1)*ngap+igap

      wave_R(:,:,igap1)=wave_R(:,:,igap1)+GGR
      wave_S(:,:,igap1)=wave_S(:,:,igap1)+GGS
   end do
end do

3  format(A25,i4,A8,i4,A8,i4,a15,i4,a6,i4,a4,f6.2,a2)

   call srite("rec_out",wave_R,8*sou%w%n*nx*ngap*nreal)
   call srite("sou_out",wave_S,8*sou%w%n*nx*ngap*nreal)

close(2)
close(3)
close(4)
close(5)

deallocate(GGR,GGS,wave_R,wave_S)
deallocate(vel)

call sep_end_prog()
end program OW_arealsht
