 program OW_arealsht

 use sep 
 use OW_type
 use OW_parms_areal
 use OW_exec_time
 use OW_taper_areal
 use OW_ctrillwei
 use OW_grabvel_areal
 use OW_rndseed1
 use OW_encoding
 use OW_shaper
 use lloyd_vel

 implicit none

 real						:: start,finish,start1,finish1,a(2),b,c
 integer						:: i,nv,ixs,j,nsy,iw,ih,icomb,ix,stat,iz,ihor,n_enc_sog,iq,nq
 integer						:: ind_trac,i_pos,i_ini,i_fin,igap,ni
 integer						:: jgap,jxs,nnn,hij,k,ib,ic,ireal,igap1
 real, allocatable					:: velint(:),taper(:),vel(:)
 real,allocatable					:: randz(:,:,:,:),vref(:,:,:),correl(:,:),fun1(:),temp(:)
 complex,allocatable				:: source(:),GG(:,:)
 complex, allocatable				:: GGG(:,:,:),tfld(:),wave_R(:,:),wave_S(:,:)
 complex, allocatable				:: 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()
 start=0.;finish=0. 
!--------------------------------------------------
! Encoding functions
!
 if (encode) then
  call owexec_time(start)
  allocate(randz(sou%w%n,sog%xm%n,nhorz,nreal*ngap))
    if(rand) then
       call init_random_seed4(randz)
    else if (gold) then
       call owgoldpw (randz,nhorz)
    else if (exist_file("InCode")) then
       call sreed("InCode",randz,4*sou%w%n*sog%xm%n*nreal*nhorz*ngap)
    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),vel(nz*kxm%n))
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_velim(velint,vel,ni)

 deallocate(velint)
 if (extrap(1:4).eq."pspi") then
    if (lloyd) then 
       allocate(vref(nref,2,nz))
       call lloyd_vel1(vel,vref)
    end if
 end if
!--------------------------------------------------
! Read source
!
 allocate(source(sou%w%n))
 call sreed("sou",source,8*sou%w%n)
 allocate(phase(sou%w%n))
 phase = cmplx(cos(pi/4.), -sin(pi/4.)) * sqrt(sou%w%all)
 source=source*phase
 deallocate(phase)
!--------------------------------------------------
! Computing taper
!
 allocate(taper(kxm%n))
 call owtaper_init()
 call owtaper(taper)
 if (exist_file("green_s")) call send_sou(source)
!--------------------------------------------------
! Wavefield extrapolation
!
 allocate(GG(kxm%n,sou%w%n))
 if (encode.and.(rand.or.exist_file("InCode"))) then
    allocate(wave_S(kxm%n,sou%w%n),wave_R(kxm%n,sou%w%n))
    do ireal=1,nreal
       do igap=1,ngap
          wave_S=0;wave_R=0
          do ihor=1,nhorz
             if (lloyd) then
                call owexec_time(start)
                GG=0
                call owctrlwei_init(ihor,ireal)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=+1,&
                            rand1=randz(:,:,ihor,(ireal-1)*ngap+igap),vref1=vref,igap=igap)
                wave_R=wave_R+GG
                call owexec_time(finish)
                write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
                call owexec_time(start)
                GG=0
                call owctrlwei_init(ihor,ireal)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=-1,&
                            rand1=randz(:,:,ihor,(ireal-1)*ngap+igap),vref1=vref,igap=igap)
                wave_S=wave_S+GG
                call owexec_time(finish)
                write(0,4) "Source wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
             else
                call owexec_time(start)
                GG=0
                call owctrlwei_init(ihor,ireal)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=+1,&
                            rand1=randz(:,:,ihor,(ireal-1)*ngap+igap),igap=igap)
                wave_R=wave_R+GG
                call owexec_time(finish)
                write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
                call owexec_time(start)
                GG=0
                call owctrlwei_init(ihor,ireal)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=-1,&
                            rand1=randz(:,:,ihor,(ireal-1)*ngap+igap),igap=igap)
                wave_S=wave_S+GG
                call owexec_time(finish)
                write(0,4) "Source wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
             end if
          end do
          if (shaper) call owshaper(wave_S,-1)
          if (shaper) call owshaper(wave_R,+1)
          call srite("rec_out",wave_R(padtraces+tapertr+1:kxm%n-(padtraces+tapertr),:),8*sou%w%n*(kxm%n-2*(padtraces+tapertr)))
          call srite("sou_out",wave_S(padtraces+tapertr+1:kxm%n-(padtraces+tapertr),:),8*sou%w%n*(kxm%n-2*(padtraces+tapertr)))
       end do
    end do
else if (encode.and.gold.or.encode.and.sgnpi) then
    allocate(wave_S(kxm%n,sou%w%n),wave_R(kxm%n,sou%w%n))
    do igap=1,ngap
      wave_S=0;wave_R=0
      do ihor=1,nhorz
         if (lloyd) then
            GG=0
            call owexec_time(start)
            call owctrlwei_init(ihor,1)
            call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=+1,rand1=randz(:,:,ihor,1),vref1=vref,igap=igap)
            call owexec_time(finish)
            wave_R=wave_R+GG
            write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
            GG=0
            call owexec_time(start)
            call owctrlwei_init(ihor,1)
            call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=-1,rand1=randz(:,:,ihor,1),vref1=vref,igap=igap)
            call owexec_time(finish)
            write(0,4) "Source wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
            wave_S=wave_S+GG
         else
            GG=0
            call owexec_time(start)
            call owctrlwei_init(ihor,1)
            call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=+1,rand1=randz(:,:,ihor,1),igap=igap)
            call owexec_time(finish)
            wave_R=wave_R+GG
            write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
            GG=0
            call owexec_time(start)
            call owctrlwei_init(ihor,1)
            call owctrlwei(taper=taper,vsal=vsal,vel=vel,GG=GG,sign=-1,rand1=randz(:,:,ihor,1),igap=igap)
            call owexec_time(finish)
            wave_S=wave_S+GG
            write(0,4) "Source wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' igap ',igap,' FROM ',ngap,' horz ',ihor,' FROM ',nhorz,' in ',finish-start,' s'
         end if
      end do
      call owshaper(wave_S,-1)
      call owshaper(wave_R,+1)
      call srite("rec_out",wave_R(padtraces+tapertr+1:kxm%n-(padtraces+tapertr),:),8*sou%w%n*(kxm%n-2*(padtraces+tapertr)))
      call srite("sou_out",wave_S(padtraces+tapertr+1:kxm%n-(padtraces+tapertr),:),8*sou%w%n*(kxm%n-2*(padtraces+tapertr)))
   end do
!else if (.not.encode) then
!   do igap=1,ngap
!      GGS=0;GGR=0
!      if (lloyd) then
!         call owexec_time(start)
!         call owctrlwei(taper=taper,vel=vel,GG=GGR,sign=+1,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,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,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,igap=igap)
!         call owexec_time(finish)
!         write(0,3) "Source wavefield   - super  data # ",igap,' FROM ',ngap,' in ',finish-start,' s'
!      end if
!
!      call owshaper(source,GGS)
!      call owshaper(source,GGR)
!    
!      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)
4  format(A35,i4,5(A6,i4),A4,f6.2,a2)
 if (exist_file("green_s")) call clos_sou()
!--------------------------------------------------
! Finishing ...
!
 close(2)
 close(3)
 close(4)
 close(5)

 deallocate(GG,wave_R,wave_S,taper)
 deallocate(vel)

 call sep_end_prog()
 end program OW_arealsht
