program OWmig2d1_par

use sep 
use OW_pspi2
use OW_type
use OW_param
use OW_fftw

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

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

integer, external 				:: omp_get_thread_num, omp_get_num_procs

 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)

call ow_parms()

nnode=omp_get_num_procs()
write(0,*) nnode

 open(unit=2,file="lixo")

 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,rec%x%n))
 allocate(vel(v%z%n,rec%h%n))
 allocate(GG(v%z%n,rec%x%n,rec%h%n,rec%w%n))


write(0,*) "passei alocacao"

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

allocate(tfield(kxm%n))
call owfftw(tfield)
deallocate(tfield)

write(0,*) "passei fft"
!--------------------------------------------------
! 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)

write(0,*) "passei leitura"
!--------------------------------------------------
! 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;
    vel=0;
    vel=velint(:,nv:nv+rec%h%n-1);
    write(2,*) maxval(velint(:,ceiling(rec%h%n/2.)))
!
    do iz=1,v%z%n
        if (nref .eq. 1) then
           vref(nref,iz,ixs)=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,ixs)=minv+(iref-1)*dvel
           end do 
        end if
    end do 
end do 

write(0,*) "passei 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)

  call sreed("in",recv,8*rec%h%n*rec%w%n)

write(0,*) "passei leitura dado"

    GG=0.
! workspace origin 
    owrk=(ixs-1)*rec%x%d+rec%x%o+rec%h%o
! grab velocity in the workspace       
    nv = ceiling((owrk-v%xm%o)/v%xm%d)
    ni = nv-ceiling((img%xm%o-v%xm%o)/img%xm%d)

    vel=0;
    vel(:,:)=velint(:,nv:nv+rec%h%n-1);

write(0,*) "passei grabvel"
!      
!$OMP parallel do default(shared) private(iw,nref,ni,isp) schedule(dynamic,1) ordered !firstprivate(nref,ni,isp) 

        

! $OMP do private(iw)

    omega:do iw = 120,rec%w%n

!        inode=omp_get_thread_num();
        call owpspi_init(iw,nref,isp,ixs)
        call owpspi(vel,vref,recv,source,GG(:,:,:,iw))

        if (mod(iw,10)==1) write(0,*)  "DONE FREQUENCY # ",iw," OUT OF ",rec%w%n," / SHOT # ",ixs,maxval(cabs(GG(10:img%z%n,1,isp,iw)))
    end do omega

! $OMP end do
!$OMP end parallel do

write(0,*) "passei downcont"
    do ino=1,rec%w%n
       image(:,:,ni:ni+kxm%n-1)=image(:,:,ni:ni+kxm%n-1)+GG(:,:,:,ino);
    end do
write(0,*) "passei image"
    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)
close(2)
 call sep_end_prog()
 end program OWmig2d1_par
