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_mseq

implicit none

real						:: start, finish,start1, finish1,a(2),b,c,m
integer						:: i,nv,ixs,j,ni,nsy,iw,ih,icomb,ix,stat,iz,ihor,n_enc_sog
integer						:: ind_trac,i_pos,i_ini,i_fin,igap,jgap,jxs,nnn,hij,k,ib,ic
integer,allocatable					:: window(:,:),HD(:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:),rectmp(:,:),soutmp(:,:)
real,allocatable					:: reccig(:,:),soucig(:,:),soucig_rot(:,:),source1(:),temp(:)
real,allocatable					:: tempr(:,:),temps(:,:),randz(:,:,:),randw(:,:),rand1(:)
real,allocatable					:: hold(:,:),gol(:,:,:)
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.)

!--------------------------------------------------

allocate(source(sou%w%n))
call sreed("sou",source,sou%w%n*8)
deallocate(source)

!---------------------------------------------------
!
! Read velocities
!
allocate(velint(v%z%n,v%xm%n))
call sreed("vel",velint,v%z%n*v%xm%n*4)

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)!/2
deallocate(velint)

!--------------------------------------------------
! Picks
!
call owpick()

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

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

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

if (allocated(reccig)) deallocate(reccig);allocate(reccig(sog%z%n,nx));reccig=0
if (allocated(soucig)) deallocate(soucig);allocate(soucig(sog%z%n,nx));soucig=0

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

if (encode) then
if(rand.and.encz) then 
   allocate(randz(sog%z%n,sou%w%n,sog%xm%n))
   call init_random_seed3(randz)
else if (rand.and.encw) then
   allocate(randw(sou%w%n,sog%xm%n))
   call init_random_seed2(randw)
else if (gold.and.encz) then
   allocate(hold(50,2*sog%z%n-subf))
   allocate(gol(2*sog%z%n,sou%w%n,sog%xm%n))
   allocate(randz(sou%w%n,sog%xm%n,2))
   allocate(rand1(2*sog%z%n-subf))
   do i=1,50
      rand1=0
      call owmseq(2,9,rand1,2*sog%z%n,i-int((i-1)/50)*50)
      hold(i,:)=rand1
   end do
   deallocate(rand1)
   hold=hold*2-1
   write(0,*) size(randz,1),size(randz,2),size(randz,3)
   call init_random_seed3(randz)
   write(0,*) size(randz,1),size(randz,2),size(randz,3)
   do iw=1,sou%w%n
      do i=1,sog%xm%n
         b=(randz(iw,i,1)+1)/2;ib=int(b*50)
         c=(randz(iw,i,2)+1)/2;ic=int(c*50)
         if (ib==0) ib=1;if (ic==0) ic=1
         if (b==c) then
            gol(1:2*sog%z%n-subf,iw,i)=hold(ib,:)*cshift(hold(ic+1,:),i+shift)
         else
            gol(1:2*sog%z%n-subf,iw,i)=hold(ib,:)*cshift(hold(ic,:),i+shift+ib)
         end if
      end do
   end do
   deallocate(hold)

   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         if (gol(1,iw,ixs).eq.1) then
            gol(2*sog%z%n,iw,ixs)=-1
         else
            gol(1,iw,ixs)=1;gol(2*sog%z%n,iw,ixs)=1
         end if
      end do
   end do

   allocate(tfld(2*sog%z%n))

!--------------------------------------------------
! Setting FFTW plans
!
   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         tfld=0
         tfld=cmplx(gol(:,iw,ixs),0)
         call sfftw_execute(planF)
         randz(:,iw,ixs)=atan(imag(tfld(sog%z%n+1:2*sog%z%n))*real(tfld(sog%z%n+1:2*sog%z%n))/ &
                            (real(tfld(sog%z%n+1:2*sog%z%n))*real(tfld(sog%z%n+1:2*sog%z%n))+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   
else if(gold.and.encw) then
   allocate(hold(50,2*sou%w%n-subf))
   allocate(gol(sog%z%n,2*sou%w%n,sog%xm%n))
   allocate(randz(sog%z%n,sog%xm%n,2))
   allocate(rand1(2*sou%w%n-subf))
   do i=1,50
      rand1=0
      call owmseq(2,9,rand1,2*sou%w%n,i-int((i-1)/50)*50)
      hold(i,:)=rand1
   end do
   deallocate(rand1)
   hold=hold*2-1
   call init_random_seed3(randz)
   do ixs=1,sog%xm%n
      do iw=1,sog%z%n
         b=(randz(iw,ixs,1)+1)/2;ib=int(b*50)
         c=(randz(iw,ixs,2)+1)/2;ic=int(c*50)
         if (ib==0) ib=1
         if (ic==0) ic=1
         if (b==c) then
            gol(iw,1:2*sou%w%n-subf,ixs)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
         else
            gol(iw,1:2*sou%w%n-subf,ixs)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib)
         end if
      end do
   end do
   deallocate(hold)

   do ixs=1,sog%xm%n
      do iw=1,sog%z%n
         if (gol(iw,1,ixs).eq.1) then
            gol(iw,2*sou%w%n,ixs)=-1
         else
            gol(iw,1,ixs)=1;gol(iw,2*sou%w%n,ixs)=1
         end if
      end do
   end do

   allocate(tfld(2*sou%w%n))

!--------------------------------------------------
! Setting FFTW plans
!
   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,sog%z%n
         tfld=0
         tfld=cmplx(gol(iw,:,ixs),0)
         call sfftw_execute(planF)
         randz(iw,:,ixs)=atan(imag(tfld(sou%w%n+1:2*sou%w%n))*real(tfld(sou%w%n+1:2*sou%w%n))/ &
                            (real(tfld(sou%w%n+1:2*sou%w%n))*real(tfld(sou%w%n+1:2*sou%w%n))+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   

end if
call owexec_time(finish)
write(0,*) "Finished generating encode functions in ",finish-start,"(s)"
end if

   if (allocated(wave_S)) deallocate(wave_S);allocate(wave_S(sou%w%n,nx,ngap));wave_S=0
   if (allocated(wave_R)) deallocate(wave_R);allocate(wave_R(sou%w%n,nx,ngap));wave_R=0
   if (allocated(GGS)) deallocate(GGS);allocate(GGS(sou%w%n,kxm%n));GGS=0
   if (allocated(GGR)) deallocate(GGR);allocate(GGR(sou%w%n,kxm%n));GGR=0

   n_enc_sog=int(sog%xm%n/ngap)

   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: adjacent cigs, in total ngap groups
! The old one is better
      if (encode.and.xiao) then
         igap=int((ixs-1)/n_enc_sog)+1
         jgap=ixs-int(ixs/n_enc_sog)*n_enc_sog
      else
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Selection scheme: cigs separated by ngap
         igap=ixs-int((ixs-1)/ngap)*ngap
         jgap=floor((ixs-1)/ngap*1.)+1
      end if

      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,:))
            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
                   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 owexec_time(start)
      call owctrlwei_init(ixs)

      if ((rand.and.encz).or.(gold.and.encz)) then
         call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,(/0.,0./),randz(:,:,ixs))

         write(0,*) "*******************"

      else if((rand.and.encw).or.(gold.and.encw)) then
         call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,randw(:,ixs))
      else
         call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS)
      end if

      call owexec_time(finish)

      write(0,3) "CIGs UPWARD PROPAGATED AT ",(ixs-1)*sog%xm%d+sog%xm%o, &
                 " # ",ixs," OUT OF ",sog%xm%n,' FROM IGAP ',igap,' in ',finish-start,' s'
3     format(A26,f10.2,A3,i3,A8,i5,a11,i4,a4,f10.2,a2)

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

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

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

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

call sep_end_prog()
end program OW_arealsht
