program OWmig2d2

use sep 
use OW_type
use OW_parms
use OW_pspi3

integer						:: time0(8), time1(8)
real						:: start_time, finish_time

integer						:: ixs,iz,ikx,iref,ixw,ino
integer						:: ni,nv,isp,inode,nnode
real						:: sign, sqa
real						:: minv, maxv, dvel, owrk
real, allocatable				:: vref(:,:),vel(:,:),velint(:,:)
complex, allocatable				:: recv(:,:),source(:,:,:),image(:,:,:),GG(:,:,:),G1(:,:,:,:)
complex, pointer				:: tfield(:)

!type (axis) 					:: kxm, h, w, x, xm, z, sh
!type (dat_grid)				:: rec, sou
!type (mod_grid)				:: img
!type (velo)					:: v

integer, external 				:: omp_get_thread_num,omp_get_num_procs,omp_get_num_threads

integer, parameter   				:: FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
integer, parameter   				:: FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
integer, parameter   				:: FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
integer, parameter   				:: FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
integer, parameter   				:: FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
integer, parameter   				:: FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
integer, parameter   				:: FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
integer, parameter   				:: FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536
 
!integer*8, pointer,private   			:: planF,planI
!integer*8, private	   			:: planF,planI
integer*8		   			:: planF,planI  

 call sep_init()
 call sep_begin_prog()

start_time=0.;finish_time=0.
call date_and_time(values=time0)
start_time = time0 (5) * 3600 + time0 (6) * 60 + time0 (7) + 0.001 * time0 (8)
!---------------------------------------
! IO
!
 call owparam()


nnode=omp_get_num_procs()
write(0,*) "# of nodes = ",nnode

 allocate(recv(rec%h%n,rec%w%n))
 allocate(source(sou%h%n,sou%w%n,sou%x%n))
 allocate(velint(v%z%n,v%xm%n))
 allocate(image(img%z%n,img%sh%n,img%xm%n))
 allocate(vref(nref,v%z%n))
 allocate(vel(v%z%n,rec%h%n))
 allocate(GG(v%z%n,img%sh%n,rec%h%n))
 allocate(G1(v%z%n,img%sh%n,rec%h%n,rec%w%n))

!--------------------------------------------------
! Setting FFTW plans
!

allocate(tfield(kxm%n))
call sfftw_plan_dft(planF,1,kxm%n,tfield,tfield,fftW_FORWARD,fftW_PATIENT)
call sfftw_plan_dft(planI,1,kxm%n,tfield,tfield,fftw_BACKWARD,fftW_PATIENT)
deallocate(tfield)

!--------------------------------------------------
! Reading source fctn and velocity

 call sreed("sou",source,sou%h%n*sou%w%n*sou%x%n*8)
 call sreed("vel",velint,v%z%n*v%xm%n*4)

!--------------------------------------------------
! Pre-compute vref
!

call date_and_time(values=time1)
finish_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)
write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish_time - start_time)/60./60.,' (hrs)'

!--------------------------------------------------
!
!
write(0,*) "DOWNWARD CONTINUES SOURCES AND RECEIVERS"
write(0,*) "========================================"
!
!--------------------------------------------------
! shot position in the acq.geom.
!
isp = abs(floor(rec%h%o/rec%h%d)) + 1

!--------------------------------------------------
! Migrating
!

shot:do ixs=1,rec%x%n

start_time=0.;finish_time=0.
call date_and_time(values=time0)
start_time = time0 (5) * 3600 + time0 (6) * 60 + time0 (7) + 0.001 * time0 (8)

owrk=(ixs-1)*rec%x%d+rec%x%o+rec%h%o
if (owrk < img%xm%o .or. owrk > img%xm%o+img%xm%d*(img%xm%n-1)) call seperr('WORKING SPACE GREATER THAN IMAGE SPACE FOR SHOT # ',ixs)
! grab velocity in the workspace       
nv=floor((owrk-v%xm%o)/v%xm%d)+1;
vref=0;vel=0
vel=velint(:,nv:nv+rec%h%n-1);
!
do iz=1,v%z%n
    if (nref .eq. 1) then
       vref(nref,iz)=sum(vel(iz,:))/rec%h%n
    else
       minv=minval(vel(iz,:))
       maxv=maxval(vel(iz,:))
       dvel=(maxv-minv)/(nref-1) 
       do iref=1,nref
          vref(iref,iz)=minv+(iref-1)*dvel
       end do 
    end if
end do 

! call sreed_window("in",5,(/rec%h%n,1,rec%w%n,rec%x%n,1/),(/rec%h%n,1,rec%w%n,1,1/),(/0,0,ixs-1,0,0/),(/1,1,1,1,1/),8,recv)
  call sreed("in",recv,8*rec%h%n*rec%w%n)
write(0,*) "Read shot # ",ixs

GG=0.

write(0,*) "Get velocity "
!---------------------------------------------------------      
! $OMP parallel do default(shared) private(inode,iw) !firstprivate(planF,planI)

!        write(0,*) inode,isp,planF,planI
! $OMP do private(iw)

    omega:do iw = 1,rec%w%n
        inode=omp_get_thread_num()

!        write(0,5)  "ENTER FREQUENCY # ",iw," OUT OF ",rec%w%n

        call owpspi3(iw,isp,planF,planI,recv,source,vref,vel,GG)
        G1(:,:,:,iw)=GG

!        if (mod(iw,10)==1) write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",rec%w%n," / SHOT # ",ixs," ",inode,maxval(cabs(G1(:,1,isp,iw)))
         write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",rec%w%n," / SHOT # ",ixs," ",inode,maxval(cabs(G1(:,1,isp,iw)))
5       format(a17,i3,a8,i3,a10,i3,a8,i3,f12.5)

 
   end do omega

! $OMP end do
! $OMP end parallel do

write(0,*) "Wavefields ropagated for shot # ",ixs
ni = nv-floor((img%xm%o-v%xm%o)/img%xm%d)+1
do ino=1,rec%w%n
image(:,:,ni:ni+kxm%n-1)=image(:,:,ni:ni+kxm%n-1)+G1(:,:,:,ino)
end do

write(0,*) "Final image for shot # ",ixs

    call date_and_time(values=time1)
    finish_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)
    write(0,*) " DONE SHOT # ",ixs," OUT OF ",rec%x%n," in ", (finish_time - start_time)/60./60.," (hrs)"

end do shot
call srite("out",real(image),4*img%z%n*img%sh%n*img%xm%n)

 call sep_end_prog()
 end program OWmig2d2
