program OWmod2d_par

use sep 
use OW_pspi
use OW_fftw
use OW_type
use OW_exec_time
use OW_parms
!use OW_imageIO
use OW_taper
use sep3d_struct_mod

implicit none

real						:: start, finish, start1, finish1, owrk1
real, allocatable					:: velint(:,:),vel(:,:),vref(:,:,:),taper(:)
real, allocatable					:: vreflect(:,:),reflectivity(:,:),refl(:,:)
complex, allocatable				:: source(:),source1(:,:)
complex, allocatable				:: g1(:,:,:),g2(:,:,:,:),shot(:,:)
integer						:: ixs,ixr,ixm,iz,minx,maxx,nv,iw,ni1,nv1,iblock
integer, external 					:: omp_get_thread_num,omp_get_num_threads

 call sep_init()
 call sep_begin_prog()

!---------------------------------------
! IO
!
 call owparam()
 allocate(taper(kxm%n));taper=1.
 call owtaper(taper)

!----------------------------------------------------
! Memory Allocation
!
allocate(source(sou%w%n))
allocate(vel(nzmax,kxm%n),velint(v%z%n,v%xm%n))
!--------------------------------------------------
! Reading source fctn and velocity
 call sreed("sou",source,sou%w%n*8)
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 if (exist_file("vref")) then
    allocate(vref(nref,2,v%z%n))
    call sreed("vref",vref,v%z%n*2*nref*4)
 end if

if (compute) then
!--------------------------------------------------
!
!
write(0,*) "DOWNWARD CONTINUES SOURCES"
write(0,*) "=========================="

allocate(source1(kxm%n,sou%w%n))
source1=0;source1(isp,:)=source
deallocate(source)
!--------------------------------------------------
! DOWNWARD CONTINUING SOURCES
!
   do ixs=1,rec%x%n+2*(ixp-1)
       start=0.;finish=0.
       call owexec_time(start)
       call owgrab_vel(ixs,velint,vel)
       call owpspi_init(ixs,wbott)
       if (exist_file("vref")) then
          call owpspi(source1,taper,vel,vref)
       else
          call owpspi(source1,taper,vel)
       end if
       call owexec_time(finish)
       write(0,*) " GF # ",ixs," OUT OF ",rec%x%n+2*(ixp-1)," in ", (finish - start)," s"
   end do

deallocate(taper,vel,velint,source1)

!--------------------------------------------------
!
!
write(0,*) "CONVOLVING GFs AND REFLECTIVITY"
write(0,*) "==============================="
if (node.eq.9999) then
  !$OMP PARALLEL
  node = omp_get_num_threads()
  !$OMP END PARALLEL
end if
write(0,*) " NUM.THREADS = ",node

!--------------------------------------------------
! CONVOLVING GFs AND REFLECTIVITY
!
allocate(reflectivity(v%z%n,v%xm%n))
allocate(refl(img%z%n,img%xm%n))
if (exist_file("reflectivity")) then
   call sreed("reflectivity",reflectivity,4*v%z%n*v%xm%n)
else
   allocate(vreflect(v%z%n,v%xm%n))
   call sreed("vreflect",vreflect,v%z%n*v%xm%n*4)
   call owreflc(vreflect,reflectivity)
   nv=(img%xm%o-v%xm%o)/img%xm%d+1.5
   refl=reflectivity(1:img%z%n,nv:nv+img%xm%n-1)
   deallocate(vreflect)!,reflectivity)
end if
 call sep3df_delete(axgr,"ws")
 call init_sep3d("ws",axgr,"SCRATCH")
!allocate(g1(img%xm%n,sou%w%n,img%z%n),g2(img%xm%n,sou%w%n,img%z%n,1),shot(rec%h%n,sou%w%n))
!do ixs=ixp,rec%x%n+ixp-1
!   start=0.;finish=0.;g1=0;shot=0.
!   call owexec_time(start)
!   call sseek("ws",0,0)
!   call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,1/), &
!        (/0,0,0,ixs-1/),(/1,1,1,1/),8,g1)
!   minx=ixs-ixp+1
!   maxx=minx+rec%h%n-1
!   do ixr=minx,maxx
!      g2=0
!      call sseek("ws",0,0)
!      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,1/),&
!           (/0,0,0,ixr-1/),(/1,1,1,1/),8,g2)
!      !$OMP parallel do private (iw,iz,ixm) schedule(dynamic,1) num_threads(node)
!      do iw=1,sou%w%n
!         do iz=1,img%z%n
!            do ixm=1,img%xm%n
!               shot(ixr-minx+1,iw)=shot(ixr-minx+1,iw)+refl(iz,ixm)*g1(ixm,iw,iz)*g2(ixm,iw,iz,1)
!            end do
!         end do
!      end do
!      !$OMP END PARALLEL DO
!      write(0,*) "Computed channel", ixr-minx+1, "from shot ",ixs-ixp+1
!   end do   
!   call owexec_time(finish)
!   call srite("shots",shot,8*rec%h%n*sou%w%n)
!   write(0,*) " GF # ",ixs-ixp+1," OUT OF ",rec%x%n," in ", (finish - start)," s"
!end do



!allocate(g1(img%xm%n,sou%w%n,img%z%n),g2(img%xm%n,sou%w%n,img%z%n,nblock),shot(rec%h%n,sou%w%n))
allocate(g1(img%xm%n,sou%w%n,img%z%n),g2(401,sou%w%n,img%z%n,nblock),shot(rec%h%n,sou%w%n))
do ixs=ixp,rec%x%n+ixp-1
   start=0.;finish=0.;g1=0;shot=0.
   call owexec_time(start)
   call sseek("ws",0,0)
   call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,1/), &
        (/0,0,0,ixs-1/),(/1,1,1,1/),8,g1)
   minx=ixs-ixp+1
   maxx=minx+rec%h%n-1+3
   do ixr=minx,maxx-nblock+1,nblock
      g2=0
      owrk1=(ixr-1-200)*rec%x%d+rec%x%o
      nv1 = (owrk-v%xm%o)/v%xm%d+1.5
      ni1 = (owrk-img%xm%o)/img%xm%d+1.5
      call sseek("ws",0,0)
!      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,nblock/),&
!           (/0,0,0,ixr-1/),(/1,1,1,1/),8,g2)
      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/401,sou%w%n,img%z%n,nblock/),&
           (/ni1-1,0,0,ixr-1/),(/1,1,1,1/),8,g2)
      !$OMP parallel do private (iw,iz,ixm) schedule(dynamic,1) num_threads(node)
      do iblock=1,nblock
      do iw=1,sou%w%n
         do iz=1,img%z%n
!            do ixm=1,img%xm%n
            do ixm=1,401
               shot(ixr-minx+1,iw)=shot(ixr-minx+1,iw)+reflectivity(iz,nv1+ixm-1)*g1(ni1+ixm-1,iw,iz)*g2(ixm,iw,iz,iblock)
            end do
         end do
      end do
      end do
      !$OMP END PARALLEL DO
      write(0,*) "Computed channels", ixr-minx+1,'-',ixr-minx+7, "from shot ",ixs-ixp+1
   end do   
   call owexec_time(finish)
   call srite("shots",shot,8*rec%h%n*sou%w%n)
   write(0,*) " GF # ",ixs-ixp+1," OUT OF ",rec%x%n," in ", (finish - start)," s"
end do

















else

!--------------------------------------------------
!
!
write(0,*) "CONVOLVING GFs AND REFLECTIVITY"
write(0,*) "==============================="
if (node.eq.9999) then
  !$OMP PARALLEL
  node = omp_get_num_threads()
  !$OMP END PARALLEL
end if
write(0,*) " NUM.THREADS = ",node

!--------------------------------------------------
! CONVOLVING GFs AND REFLECTIVITY
!
allocate(reflectivity(v%z%n,v%xm%n))
allocate(refl(img%z%n,img%xm%n))
if (exist_file("reflectivity")) then
   call sreed("reflectivity",reflectivity,4*v%z%n*v%xm%n)
else
   allocate(vreflect(v%z%n,v%xm%n))
   call sreed("vreflect",vreflect,v%z%n*v%xm%n*4)
   call owreflc(vreflect,reflectivity)
   nv=(img%xm%o-v%xm%o)/img%xm%d+1.5
   refl=reflectivity(1:img%z%n,nv:nv+img%xm%n-1)
   deallocate(vreflect)!,reflectivity)
end if
 call sep3df_delete(axgr,"ws")
 call init_sep3d("ws",axgr,"SCRATCH")
allocate(g1(img%xm%n,sou%w%n,img%z%n),g2(img%xm%n,sou%w%n,img%z%n,nblock),shot(rec%h%n,sou%w%n))
do ixs=ixp,rec%x%n+ixp-1
   start=0.;finish=0.;g1=0;shot=0.
   call owexec_time(start)
   call sseek("ws",0,0)
   call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,1/), &
        (/0,0,0,ixs-1/),(/1,1,1,1/),8,g1)
   minx=ixs-ixp+1
!   maxx=minx+rec%h%n-1+3
!   do ixr=minx,maxx-nblock,nblock
!      g2=0
!      call sseek("ws",0,0)
!      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,nblock/),&
!           (/0,0,0,ixr-1/),(/1,1,1,1/),8,g2)
!      !$OMP parallel do private (iw,iz,ixm) schedule(dynamic,1) num_threads(node)
!      do iblock=1,nblock
!      do iw=1,sou%w%n
!         do iz=1,img%z%n
!            do ixm=1,img%xm%n
!               shot(ixr-minx+1,iw)=shot(ixr-minx+1,iw)+refl(iz,ixm)*g1(ixm,iw,iz)*g2(ixm,iw,iz,iblock)
!            end do
!         end do
!      end do
!      end do
!      !$OMP END PARALLEL DO
!      write(0,*) "Computed channels", ixr-minx+1,'-',ixr-minx+7, "from shot ",ixs-ixp+1
!   end do   
!   call owexec_time(finish)
!   call srite("shots",shot,8*rec%h%n*sou%w%n)
!   write(0,*) " GF # ",ixs-ixp+1," OUT OF ",rec%x%n," in ", (finish - start)," s"
!end do


   maxx=minx+rec%h%n-1-7
   do ixr=minx,maxx-nblock+1,nblock
      g2=0
      owrk1=(ixr-1-200)*rec%x%d+rec%x%o
      nv1 = (owrk-v%xm%o)/v%xm%d+1.5
      ni1 = (owrk-img%xm%o)/img%xm%d+1.5
      call owexec_time(start1)      
      call sseek("ws",0,0)
!      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,nblock/),&
!           (/0,0,0,ixr-1/),(/1,1,1,1/),8,g2)
      call sreed_window("ws",4,(/img%xm%n,sou%w%n,img%z%n,rec%x%n+2*(ixp-1)/),(/img%xm%n,sou%w%n,img%z%n,nblock/),&
           (/0,0,0,ixr-1/),(/1,1,1,1/),8,g2)
      call owexec_time(finish1)      
      write(0,*) "Wasted ",(finish1-start1)," s to read in data"
      call owexec_time(start1)      
      !$OMP parallel do private (iw,iz,ixm) schedule(dynamic,1) num_threads(node)
      do iblock=1,nblock
      do iw=1,sou%w%n
         do iz=1,img%z%n
!            do ixm=1,img%xm%n
            do ixm=1,401
               shot(ixr-minx+iblock,iw)=shot(ixr-minx+iblock,iw)+reflectivity(iz,nv1+ixm-1)*&
                                               g1(ni1+ixm-1,iw,iz)*g2(ni1+ixm,iw,iz,iblock)
            end do
         end do
      end do
      end do
      !$OMP END PARALLEL DO
      write(0,*) "Computed channels", ixr-minx+1,'-',ixr-minx+nblock+1, "from shot ",ixs-ixp+1
      call owexec_time(finish1)      
      write(0,*) "Convolution in ",(finish1-start1)," s"
   end do   
   call owexec_time(finish)
   call srite("shots",shot,8*rec%h%n*sou%w%n)
   write(0,*) " GF # ",ixs-ixp+1," OUT OF ",rec%x%n," in ", (finish - start)," s"
end do



















end if

call sep_end_prog()
end program OWmod2d_par
