 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_pickbp
 use OW_rndseed1
 use OW_encoding
! use OW_gold

 implicit none

 integer						:: nhorz
 real						:: start,finish,start1,finish1,a(2),b,c
 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
 real, allocatable					:: velint(:,:),taper(:),vel(:,:)
 real,allocatable					:: randz(:,:,:,:),vref(:,:,:),correl(:,:),fun1(:)
 complex,allocatable				:: source(:),GGR(:,:),GGS(:,:)
 complex, allocatable				:: GGG(:,:,:),tfld(:)!,wave_R(:,:,:),wave_S(:,:,:)
 complex						:: phase

!-----------------------------------------
! Starting ...
!
 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)
!---------------------------------------
! Parameters
!
 call owparam(.false.)
 start=0.;finish=0. 
 call owexec_time(start)
!--------------------------------------------------
! Picks
!
 nhorz=2
! if (exist_file("hor1rec")) then
!    call owpick()
! end if

!--------------------------------------------------
! Encoding functions
!
 if (encode) then
    if (not(rand)) allocate(randz(sou%w%n,sog%xm%n,nhorz,nreal))
! Conventional random codes
    if(rand) then
       allocate(randz(sou%w%n,sog%xm%n,sog%z%n,nreal))
       call init_random_seed4(randz)
    ! Output correlations of conventional random phases
       if (exist_file("correlations")) then
          if (allocated(tfld)) deallocate(tfld);allocate(tfld(sou%w%n))
          call owfftw(tfld)
          allocate(fun1(sou%w%n))
          allocate(correl(sog%xm%n,sou%w%n))
          tfld=cmplx(randz(:,1,1,1),0.)
          call sfftw_execute(planF)
          fun1=conjg(tfld/sqrt(1.*(sou%w%n)))
          do ixs=1,sog%xm%n
             tfld=cmplx(randz(:,ixs,1,1),0.)
             call sfftw_execute(planF)
             tfld=tfld/sqrt(1.*(sou%w%n))*fun1
             call sfftw_execute(planI)
             correl(ixs,:)=tfld
             correl(ixs,:)=cshift(correl(ixs,:),int((sou%w%n)/2))
          end do
          correl=correl/maxval(correl(1,:))*2*sou%w%n
          call to_history("n1",(sog%xm%n),"correlations"); call to_history("d1",1.,"correlations"); 
          call to_history("o1",1.,"correlations")
          call to_history("n2",sou%w%n,"correlations") ; call to_history("d2",2.,"correlations") ;
          call to_history("o2",-float(sou%w%n),"correlations")
          call to_history("n3",1,"correlations"); call to_history("d3",1.,"correlations");
          call to_history("o3",1.,"correlations")  
          call to_history("n4",1,"correlations"); call to_history("d4",1.,"correlations");
          call to_history("o4",0.,"correlations")
          call to_history("n5",1,"correlations"); call to_history("d5",1.,"correlations");
          call to_history("o5",0.,"correlations")  
          call to_history("esize",4,"correlations")
          call srite("correlations",correl,4*sog%xm%n*sou%w%n)
          call owfftwdestroy()
          deallocate(tfld,fun1)
       end if
! Gold codes
    else if (gold) then
       call owgoldpw (randz,nhorz)
    else if(sgnpi) then
       randz=0.
       do iw=1,sou%w%n
          randz(iw,:,2,:)=pi*(-1)**iw
       end do
    else if (exist_file("InCode")) then
       nhorz=1
       if (allocated(randz)) deallocate(randz); allocate(randz(sou%w%n,sog%xm%n,nhorz,nreal))
       call sreed("InCode",randz,4*sou%w%n*sog%xm%n*nreal*nhorz)
    end if
    if(exist_file("OutCode")) then
      call srite("OutCode",randz,4*sou%w%n*sog%xm%n*nreal*nhorz)
    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)
!--------------------------------------------------
! Read source
!
 allocate(source(sou%w%n))
 call sreed("sou",source,8*sou%w%n)
 allocate(fun1(sou%w%n+14))
 fun1=0;fun1(8:sou%w%n+7)=cabs(source)
 fun1(1:7)=cabs(source(1))
 fun1(sou%w%n+8:sou%w%n+14)=cabs(source(sou%w%n))
source=0;
 do i=1,sou%w%n
    source(i)=cmplx(sum(fun1(i:i+14))/15,0.)
 end do
 deallocate(fun1)
call to_history("n1",sou%w%n,"wavelet"); call to_history("d1",sou%w%d/(2*pi),"wavelet"); call to_history("o1",sou%w%o/(2*pi),"wavelet")
call to_history("esize",8,"wavelet")
call srite("wavelet",source,8*sou%w%n)
!--------------------------------------------------
! Computing taper
!
 allocate(taper(kxm%n))
 call owtaper(int(padtraces/4),taper)
!--------------------------------------------------
! Wavefield extrapolation
!
 allocate(GGS(sou%w%n,kxm%n),GGR(sou%w%n,kxm%n))
 call owctrlwei_init()
 if (encode.and.(rand.or.exist_file("InCode"))) then
!    allocate(wave_S(sou%w%n,nx,nreal),wave_R(sou%w%n,nx,nreal));wave_S=0;wave_R=0
    do ireal=1,nreal
       GGS=0;GGR=0
       if (exist_file("vref")) then
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGR,sign=+1,rand1=randz(:,:,:,ireal),wavelet=source,vref1=vref)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGS,sign=-1,rand1=randz(:,:,:,ireal),wavelet=source,vref1=vref)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
       else
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGR,sign=+1,rand1=randz(:,:,:,ireal),wavelet=source)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGS,sign=-1,rand1=randz(:,:,:,ireal),wavelet=source)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
       end if
!       wave_R(:,:,ireal)=GGR
!       wave_S(:,:,ireal)=GGS
!       do i=1,sou%w%n
!          GGR(i,:)=GGR(i,:)*cabs(source(i))
!          GGS(i,:)=GGS(i,:)*cabs(source(i))
!       end do
       call srite("rec_out",GGR,8*sou%w%n*nx)
       call srite("sou_out",GGS,8*sou%w%n*nx)
    end do
 else if (encode.and.gold.or.encode.and.sgnpi) then
!    allocate(wave_S(sou%w%n,nx,ngap),wave_R(sou%w%n,nx,ngap));wave_S=0;wave_R=0    
    do igap=1,ngap
       GGS=0;GGR=0
       if (exist_file("vref")) then
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGR,sign=+1,rand1=randz(:,:,:,1),wavelet=source,vref1=vref,igap=igap)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGS,sign=-1,rand1=randz(:,:,:,1),wavelet=source,vref1=vref,igap=igap)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
       else
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGR,sign=+1,rand1=randz(:,:,:,1),wavelet=source,igap=igap)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGS,sign=-1,rand1=randz(:,:,:,1),wavelet=source,igap=igap)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
       end if
!       wave_R(:,:,igap)=GGR
!       wave_S(:,:,igap)=GGS
       do i=1,sou%w%n
          GGR(i,:)=GGR(i,:)*cabs(source(i))
          GGS(i,:)=GGS(i,:)*cabs(source(i))
       end do
       call srite("rec_out",GGR,8*sou%w%n*nx)
       call srite("sou_out",GGS,8*sou%w%n*nx)
    end do
 else if (.not.encode) then
!    allocate(wave_S(sou%w%n,nx,ngap),wave_R(sou%w%n,nx,ngap));wave_S=0;wave_R=0
    do igap=1,ngap
       GGS=0;GGR=0
       if (exist_file("vref")) then
          call owexec_time(start)
          call owctrlwei(taper=taper,vel=vel,GG=GGR,sign=+1,wavelet=source,vref1=vref,igap=igap)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vel=vel,GG=GGS,sign=-1,wavelet=source,vref1=vref,igap=igap)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
       else
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGR,sign=+1,wavelet=source,igap=igap)
          call owexec_time(finish)
          write(0,3) "Receiver wavefield - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
          call owexec_time(start)
          call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GGS,sign=-1,wavelet=source,igap=igap)
          call owexec_time(finish)
          write(0,3) "Source wavefield   - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
       end if
!       do i=1,sou%w%n
!          wave_R(i,:,igap)=GGR(i,:)*abs(source(i))
!          wave_S(i,:,igap)=GGS(i,:)*abs(source(i))
!       end do
       do i=1,sou%w%n
          GGR(i,:)=GGR(i,:)*cabs(source(i))
          GGS(i,:)=GGS(i,:)*cabs(source(i))
       end do
       call srite("rec_out",GGR,8*sou%w%n*nx)
       call srite("sou_out",GGS,8*sou%w%n*nx)
    end do
 end if
3  format(A35,i4,A6,i4,A4,f6.2,a2)

!if (encode.and.rand) then
!   call srite("rec_out",wave_R,8*sou%w%n*nx*nreal)
!   call srite("sou_out",wave_S,8*sou%w%n*nx*nreal)
!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
!--------------------------------------------------
! Finishing ...
!
 close(2)
 close(3)
 close(4)
 close(5)

 deallocate(GGR,GGS)
 deallocate(vel)

 call sep_end_prog()
 end program OW_arealsht
