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_taper_areal
use OW_picksh
use OW_hadam1

implicit none

real						:: start, finish,start1, finish1
integer						:: node,i,nv,ixs,j,ni,nsy,iw,ih,icomb,ix,stat,iz,ihor
integer						:: ind_trac,i_pos,i_ini,i_fin,igap,jxs,nnn,hij,k
integer,allocatable					:: window(:,:),HD(:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:),vel1(:,:),rectmp(:,:),soutmp(:,:)
real,allocatable					:: reccig(:,:),soucig(:,:),soucig_rot(:,:),source1(:)
real,allocatable					:: tempr(:,:),temps(:,:),CHD(:,:)
complex,allocatable					:: source(:),GGR(:,:),GGS(:,:),GGRHD(:,:,:),GGSHD(:,:,:)
complex, allocatable				:: GGG(:,:,:),wave_R(:,:,:),wave_S(:,:,:),tfld(:),CCHD(:,:)

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")
node=omp_get_num_threads()

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

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

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

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

allocate(source(sou%w%n),source1(sou%w%n))
call sreed("sou",source,sou%w%n*8)
source1=cabs(source)
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

if (allocated(vel1)) deallocate(vel1);allocate(vel1(sog%z%n,kxm%n));vel1=0

!--------------------------------------------------
!
! Optionally uses Golay codes
!
  if (hadamard) then
      if (allocated(wave_S)) deallocate(wave_S);allocate(wave_S(sou%w%n,nx,2*ngap));wave_S=0
      if (allocated(wave_R)) deallocate(wave_R);allocate(wave_R(sou%w%n,nx,2*ngap));wave_R=0
      if (allocated(GGSHD)) deallocate(GGSHD);allocate(GGSHD(sou%w%n,kxm%n,2));GGSHD=0
      if (allocated(GGRHD)) deallocate(GGRHD);allocate(GGRHD(sou%w%n,kxm%n,2));GGRHD=0
      allocate(CHD(2*sog%xm%n*sou%w%n,2*sog%z%n))
      start=0.;finish=0.
      call owexec_time(start)
      call owhadamard(CHD)
      call owexec_time(finish)
      write(0,*) "Finished generating Golay codes in ",finish-start,' s'
   else
      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
   end if

   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. 
      if (hadamard) then
         GGSHD=0.;GGRHD=0.
      else
         GGS=0.;GGR=0.
      end if
      vel1=0.
      igap=ixs-int((ixs-1)/ngap)*ngap

!      call owgrab_vel_areal(ixs,vel,vel1,ni)

      jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
      write(0,*) jxs
      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 owctrl_init(ixs,node)

      if (hadamard) then
         call owctrlHD(taper,source1,reccig,soucig,vsal,vel,GGRHD,GGSHD,CHD(2*(ixs-1)*sou%w%n+1:2*ixs*sou%w%n,1:sog%z%n))
      else
         call owctrl(taper,source1,reccig,soucig,vsal,vel,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)

      if (hadamard) then
         wave_R(:,:,2*igap-1)=wave_R(:,:,2*igap-1)+GGRHD(:,:,1)
         wave_S(:,:,2*igap-1)=wave_S(:,:,2*igap-1)+GGSHD(:,:,1)
         wave_R(:,:,2*igap)=wave_R(:,:,2*igap)+GGRHD(:,:,2)
         wave_S(:,:,2*igap)=wave_S(:,:,2*igap)+GGSHD(:,:,2)
      else
         wave_R(:,:,igap)=wave_R(:,:,igap)+GGR
         wave_S(:,:,igap)=wave_S(:,:,igap)+GGS
      end if
   end do

if (hadamard) then
   call srite("rec_out",wave_R,8*sou%w%n*nx*2*ngap)
   call srite("sou_out",wave_S,8*sou%w%n*nx*2*ngap)
else
   call srite("rec_out",wave_R,8*sou%w%n*nx*ngap)
   call srite("sou_out",wave_S,8*sou%w%n*nx*ngap)
end if
close(2)
close(3)
close(4)
close(5)

if (hadamard) then
   deallocate(GGRHD,GGSHD,wave_R,wave_S)
else
   deallocate(GGR,GGS,wave_R,wave_S)
end if
deallocate(source1,vel)

call sep_end_prog()
end program OW_arealsht
