 program OW3d_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_shaper
 use lloyd_vel
use OW_rndseed1

 implicit none

 logical                                                    :: shape
 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,igapx,igapy,ni
 integer						:: jgap,jxs,nnn,hij,k,ib,ic,ireal,igap1
 real, allocatable					:: velint(:),taper(:,:),vel(:)
 real,allocatable					:: vref(:,:,:),randz(:)
 complex,allocatable				:: source(:),GG(:,:,:)
 complex						:: phase

!-----------------------------------------
! Starting ...
!
write(0,*) "before sep_init"
 call sep_init()
write(0,*) "before sep_begin_prog"
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
!---------------------------------------
! Parameters
!
write(0,*) "before owparam"
 call owparam()
write(0,*) "before auxclose"
 call auxclose("rec_out")
 call auxclose("sou_out")
write(0,*) "before auxinout"
 call auxinout("rec_out")
 call auxinout("sou_out")
 shape=.true.
! call from_param("shape",shape,.true.)

!---------------------------------------------------
!
! Read velocities
!
write(0,*) "before allocating velocities"
 allocate(velint(v%z%n*v%xm%n*v%ym%n),vel(nz*kxm%n*kym%n))
 call sreed("vel",velint,v%z%n*v%xm%n*v%ym%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
!
write(0,*) "before reading source"
 allocate(source(sou%w%n))
 call sreed("sou",source,8*sou%w%n)
!--------------------------------------------------
! Computing taper
!
write(0,*) "before computing taper"
 allocate(taper(kxm%n,kym%n))
 call owtaper(taper)

!--------------------------------------------------
! Wavefield extrapolation
!
 allocate(randz(ncigx*ncigy*wblock%ns))
     do i=1,ngapx*ngapy
!        ihor=mod(i,nhorz);if (ihor==0) ihor=nhorz
!        igapy=int((i-ihor)/1./nhorz)-int((i-1)/1./ngapy/nhorz)*ngapy+1
!        igapx=int((i-igapy)/1./nhorz/ngapy)-int((i-1)/1./ngapx/ngapy/nhorz)*ngapx+1
!        ireal=int((i-igapx)/1./nhorz/ngapy/ngapx)-int((i-1)/1./ngapx/ngapy/nhorz/nreal)*nreal+1
!        write(0,*) i,ireal,igapx,igapy,ihor 
         igapy=mod(i,ngapy);if(igapy==0) igapy=ngapy
         igapx=int((i-igapy)/1./ngapy)-int((i-1)/1./ngapx/ngapy)*ngapx+1
         do ireal=1,nreal
           do ihor=1,nhorz
             call init_random_seed3(randz)
             if (lloyd) then
                call owexec_time(start)
                call owctrlwei_init(ihor,ireal,igapx,igapy)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,sign=+1,randz=randz,vref1=vref)
                call owexec_time(finish)
                write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' gapx ',igapx,' FROM ',ngapx,' gapy ',igapy,' FROM ',ngapy,  &
                        ' horz ',ihor,' FROM ',nhorz,' in ',int(finish-start),' s'
                call owexec_time(start)
                call owctrlwei_init(ihor,ireal,igapx,igapy)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,sign=-1,randz=randz,vref1=vref)
                call owexec_time(finish)
                write(0,4) "Source   wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' gapx ',igapx,' FROM ',ngapx,' gapy ',igapy,' FROM ',ngapy,  &
                        ' horz ',ihor,' FROM ',nhorz,' in ',int(finish-start),' s'
             else
                call owexec_time(start)
                call owctrlwei_init(ihor,ireal,igapx,igapy)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,sign=+1,randz=randz)
                call owexec_time(finish)
                write(0,4) "Receiver wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' gapx ',igapx,' FROM ',ngapx,' gapy ',igapy,' FROM ',ngapy,  &
                        ' horz ',ihor,' FROM ',nhorz,' in ',int(finish-start),' s'
                call owexec_time(start)
                call owctrlwei_init(ihor,ireal,igapx,igapy)
                call owctrlwei(taper=taper,vsal=vsal,vel=vel,sign=-1,randz=randz)
                call owexec_time(finish)
                write(0,4) "Source   wavefield   - realization # ",ireal,' FROM ',nreal,&
                        ' gapx ',igapx,' FROM ',ngapx,' gapy ',igapy,' FROM ',ngapy,  &
                        ' horz ',ihor,' FROM ',nhorz,' in ',int(finish-start),' s'
             end if
!          end do
          if (ihor==nhorz) then
             if (shape) then
                allocate(GG(kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),sou%w%n))
                call owctrlwei_init(ihor,ireal,igapx,igapy)
                call gathers_io(GG,"rec_out","reed")
                call owshaper(source,GG,+1)
                call gathers_io(GG,"rec_out","rite")
                call gathers_io(GG,"sou_out","reed")
                call owshaper(source,GG,-1)
                call gathers_io(GG,"sou_out","rite")
                deallocate(GG)
             end if
          end if
!          end do
!       end do
!    end do
    end do
4  format(A35,i4,7(A6,i4),A4,i5,a2)
!--------------------------------------------------
! Finishing ...
!
 if (allocated(taper))  deallocate(taper)
 if (allocated(vel))    deallocate(vel)
 if (allocated(vref))   deallocate(vref)
 if (allocated(source)) deallocate(source)
 if (allocated(GG))     deallocate(GG)
 if (allocated(randz))     deallocate(randz)

 call sep_end_prog()
 end program OW3d_arealsht
