program OWmig2d2p

use sep 
use OW_type
use OW_param
use OW_migrate

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

integer						:: ixs,iz,iref,nv
real						:: minv, maxv, dvel, owrk
real, allocatable				:: vref(:,:,:),vel(:,:),velint(:,:),R_image(:,:,:)
complex, allocatable				:: wav(:,:,:)
complex, allocatable				:: source(:,:,:),image(:,:,:)

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 ow_parms()

 allocate(velint(v%z%n,v%xm%n))
 allocate(image(img%z%n,img%sh%n,img%xm%n))
 allocate(R_image(img%z%n,img%sh%n,img%xm%n))
 allocate(vref(nref,v%z%n,rec%x%n))
 allocate(vel(v%z%n,rec%h%n))

!--------------------------------------------------
! Reading source fctn and velocity
  allocate(source(rec%h%n,sou%w%n,sou%x%n));source=0
  allocate(wav(sou%h%n,sou%w%n,sou%x%n))
  call sreed("sou",wav,sou%h%n*sou%w%n*sou%x%n*8)
  source(isp,:,:)=wav(1,:,:)
  deallocate(wav)

 call sreed("vel",velint,v%z%n*v%xm%n*4)

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

vref=0
do ixs=1,rec%x%n
    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.5
    vel=0;
    vel=velint(:,nv:nv+rec%h%n-1);
!
    do iz=1,v%z%n
!        k=count(vel(iz,:) < vmax))
!        where(vel(iz,:) < vmax)
        if (nref .eq. 1) then
           vref(nref,iz,ixs)=sum(vel(iz,:))/rec%h%n
        else
           minv=minval(vel(iz,:))
           if (maxval(vel(iz,:)) >= vsal) then
               maxv=maxval(vel(iz,:),1,vel(iz,:) < vsal)
               dvel=(maxv-minv)/(nref-2) 
               do iref=1,nref-1
                  vref(iref,iz,ixs)=minv+(iref-1)*dvel
               end do 
               vref(nref,iz,ixs)=vsal
           else
               maxv=maxval(vel(iz,:))
               dvel=(maxv-minv)/(nref-1) 
               do iref=1,nref
                  vref(iref,iz,ixs)=minv+(iref-1)*dvel
               end do 
           end if 
        end if
    end do 
end do 
deallocate(vel)

if(exist_file("vref")) then
 call to_history("n1",nref,"vref"); call to_history("d1",1.,"vref"); call to_history("o1",0.,"vref")
 call to_history("n2",img%z%n,"vref") ; call to_history("d2",img%z%d,"vref") ; call to_history("o2",img%z%o,"vref")
 call to_history("n3",rec%x%n,"vref"); call to_history("d3",rec%x%d,"vref"); call to_history("o3",rec%x%o,"vref")  
 call to_history("esize",4)
 call srite("vref",vref,4*nref*img%z%n*rec%x%n)
end if 

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,*) "========================================"
!

call owmigra(nref,wbott,vref,velint,source,image)
R_image=real(image)
call srite("out",R_image,4*img%z%n*img%sh%n*img%xm%n)

 call sep_end_prog()
 end program OWmig2d2p
