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


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,h
integer,allocatable					:: window(:,:),HD(:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:),vel1(:,:),rectmp(:,:),soutmp(:,:)
real,allocatable					:: reccig(:,:),soucig(:,:),soucig_rot(:,:),source1(:)
real,allocatable					:: tempr(:,:),temps(:,:)
complex,allocatable					:: source(:),GGR(:,:),GGS(:,:)
complex, allocatable				:: GGG(:,:,:),wave_R(:,:,:),wave_S(:,:,:)

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

   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


   if (hadamard) then
   allocate(HD(2*sog%xm%n,sog%z%n))
   do i=1,2*sog%xm%n
      do j=1,sog%z%n
         call owhadamard(i,j,k,n,h)
         HD(i,j)=h
      end do
   end do
   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. 
      GGS=0.;GGR=0.
      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
!            tempr(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:)=tempr(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:) &
!                       /maxval(tempr(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:))
!            temps(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:)=temps(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:) &
!                       /maxval(temps(picks(ihor,ixs)-wavz:picks(ihor,ixs)+wavz,:))
         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)

!      call owctrl(taper,source1,reccig,soucig,vsal,vel,GGR,GGS,window)
      call owctrl(taper,source1,reccig,soucig,vsal,vel,GGR,GGS)
 
      call owexec_time(finish)

!      write(0,*) sum(real(GGR(40,:))),sum(real(GGR(50,:))),sum(real(GGR(60,:))),sum(real(GGS(40,:))),sum(real(GGS(50,:))),sum(real(GGS(60,:)))

      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)

deallocate(source1,vel)

call sep_end_prog()
end program OW_arealsht
