! ---------------------------------------------------------------------------
!           Lt_L_ng_3D.f90 - by Alejandro Valenciano 
!                        
!           January 2006.
! ---------------------------------------------------------------------------
!  Parameters: 
!              nax     x offset from the diagonal  (odd)
!              nay     y offset from the diagonal  (odd)
!              naz     z offset from the diagonal  (odd)
!         I/O:
!              in      green_in function (nox,noy,nz,nsux,nsuy) 
!              out     Hessian           (naz,nax,nay,nz,nx,ny) 
!                          
! Description: 
!              Computes the Hessian 
!              Parallelize in the omega dimension
! ---------------------------------------------------------------------------

program L_tL_ng_3D
use sep
use omp_lib
use omp_sector_real_bob
!use sep_timers_mod

implicit none

integer                                    :: n1,n2,n3,n4,n5,n6
real                                       :: o1,o2,o3,o4,o5,o6
real                                       :: d1,d2,d3,d4,d5,d6
integer                                    :: esize
integer                                    :: nsux,nsuy,nz,nx,ny,nax,nay,naz,nax2,nay2,naz2,nxd,nyd,ratiox,ratioy,nry,nrx,nsy,nsx
integer                                    :: ix,iy,iz,i1,i2,i3,i4,i5,ixdata,iydata
real                                       :: minix,miniy,maxix,maxiy,minidx,minidy,maxidx,maxidy,xdata,ydata
complex, dimension(:,:,:,:,:), allocatable :: green,green_input
real,    dimension(:)        , allocatable :: auxi
real   , dimension(:,:,:,:,:), allocatable :: LtL,auxi2 
integer, external                          :: srite,sreed
logical                                    :: add,marine,logic
integer                                    :: time_array_0(8), time_array_1(8),a1,a2,a3,a4,a5
real                                       :: finish_timet, start_timet

call sep_init()
call set_no_putch()

call from_history('n1',n1); call from_history('d1',d1); call from_history('o1',o1) !off_x
call from_history('n2',n2); call from_history('d2',d2); call from_history('o2',o2) !off_y
call from_history('n3',n3); call from_history('d3',d3); call from_history('o3',o3) 
call from_history('n4',n4); call from_history('d4',d4); call from_history('o4',o4) !z
call from_history('n5',n5); call from_history('d5',d5); call from_history('o5',o5) !sx
call from_history('n6',n6); call from_history('d6',d6); call from_history('o6',o6) !sy
call from_history("esize",esize)

call from_param('nax',nax,11)        ! allways odd 
call from_param('nay',nay,11)        ! allways odd 
call from_param('naz',naz,11)        ! allways odd
call from_param('minix',minix,-999999999.)  
call from_param('miniy',miniy,-999999999.)  
call from_param('nx',nx,-999999999)  
call from_param('ny',ny,-999999999)  
call from_param('ratiox',ratiox,1)       ! source to receiver ratio 
call from_param('ratioy',ratioy,1)       ! source to receiver ratio
call from_param('nrx',nrx,1)        
call from_param('nry',nry,1)
call from_param('nrx',nsx,1)        
call from_param('nry',nsy,1)
call from_param('marine',marine,.false.)  
!   logic= init_sep_timers()
!   logic= setup_next_timer("a1",a1)
!   logic= setup_next_timer("a2",a2)
!   logic= setup_next_timer("a3",a3)
!   logic= setup_next_timer("a4",a4)
!   logic= setup_next_timer("a5",a5)

     
! Parameters calculation out of the loop
if (marine) then 
  minidx=o4+o1
  minidy=o5+o2
  nxd=1.5+((n5-1)*d5+(n1-1)*d1)/d1                
  nyd=1.5+((n6-1)*d6+(n2-1)*d2)/d2                
else
  minidx=o1
  minidy=o2
  nxd=n1
  nyd=n2
endif

if (nx   ==-999999999 ) nx = nxd
if (ny   ==-999999999 ) ny = nyd
if (minix==-999999999.) minix = minidx
if (miniy==-999999999.) miniy = minidy
maxix = minix + (nx-1)*d1
maxiy = miniy + (ny-1)*d2
nsux=n5
nsuy=n6
nz=n4
naz2=(naz-1)/2
nax2=(nax-1)/2
nay2=(nay-1)/2

!check to see if this is the first time we are writing to the file----------------------------------------------------
call from_param("add",add,.false.)
if (.not. add) then
  call to_history('n1',naz,"outtag") ;  call to_history('d1',d4,"outtag")
!  call to_history('n1',naz,"outtag") ;  call to_history('d1',d3,"outtag")
  call to_history('o1',0.,"outtag")  ;  call to_history('label1',"az","outtag")
  call to_history('n2',nax,"outtag") ;  call to_history('d2',d1,"outtag"); 
  call to_history('o2',0.,"outtag" ) ;  call to_history('label2',"ax","outtag")
  call to_history('n3',nay,"outtag") ;  call to_history('d3',d2,"outtag"); 
  call to_history('o3',0.,"outtag")  ;  call to_history('label3',"ay","outtag")
  call to_history('n4',nz,"outtag")  ;  call to_history('d4',d4,"outtag"); 
!  call to_history('n4',nz,"outtag")  ;  call to_history('d4',d3,"outtag"); 
  call to_history('o4',o4,"outtag")  ;  call to_history('label4',"z","outtag")
!  call to_history('o4',o3,"outtag")  ;  call to_history('label4',"z","outtag")
  call to_history('n5',nx,"outtag")  ;  call to_history('d5',d1,"outtag"); 
  call to_history('o5',minix,"outtag"); call to_history('label5',"mx","outtag")
  call to_history('n6',ny,"outtag")  ;  call to_history('d6',d2,"outtag"); 
  call to_history('o6',miniy,"outtag"); call to_history('label6',"my","outtag")
  call to_history('n7',1,"outtag")   ;  call to_history('d7',1.,"outtag"); 
  call to_history('o7',0.,"outtag")  ;  call to_history('label7',"w","outtag")
  call to_history('esize',esize/2,"outtag")
  do iy=1,ny
    allocate(auxi(naz*nax*nay*nz*nx))
    auxi=0.
    if (srite("outtag",auxi,esize*naz*nax*nay*nz*nx/2) /= esize*naz*nax*nay*nz*nx/2) &
      call seperr("Error initializing Hessian")
    deallocate(auxi)
  enddo
  call auxclose("outtag")
end if
call auxinout("outtag")
!---------------------------------------------------------------------------------------------------------------------
call sep_begin_prog()


!call start_timer_num(a1)
!----------------------------------------------------------------------------------------------------------------------
call date_and_time(values=time_array_0)
start_timet = time_array_0 (3) * 86400.+ time_array_0 (5) * 3600  + time_array_0 (6) * 60 &
           + time_array_0 (7) + 0.001 * time_array_0 (8)
!----------------------------------------------------------------------------------------------------------------------
  allocate(green_input(n1,n2,n4,n5,n6),green(nsuy,nsux,nz+2,nx+2,ny+2))
!  allocate(green_input(n1,n2,n3,n4,n5),green(nsuy,nsux,nz+2,nx+2,ny+2))
  call sep_read(green_input)
  green=(0.,0.)
  do i1=1,n1
  do i2=1,n2
  do i3=1,n4
  do i4=1,n5
    xdata=o4+(i4-1)*d4+o1+(i1-1)*d1                                                        
    if ((xdata >= minix) .and. (xdata <= maxix)) then
    do i5=1,n6     
      ydata=o5+(i5-1)*d5+o2+(i2-1)*d2
      if ((ydata >= miniy) .and. (ydata <= maxiy)) then
        ixdata=1.5+(xdata-minix)/d1
        iydata=1.5+(ydata-miniy)/d2
        green(i5,i4,i3+1,ixdata+1,iydata+1)=green_input(i1,i2,i3,i4,i5) ! Very slow operation I need to look into this (Transpose outside?)
!        green(i5,i4,i3,ixdata,iydata)=green_input(i1,i2,i3,i4,i5) ! Very slow operation I need to look into this (Transpose outside?)
      endif
    end do
    end if
  end do
  end do
  end do
  end do

  deallocate(green_input)
!----------------------------------------------------------------------------------------------------------------------
  allocate(LtL(naz,nax,nay,nz,nx),auxi2(naz,nax,nay,nz,nx)) 
  write(0,*)"allocated memory green",nsux*nsuy*nz*nx*ny*esize 
  write(0,*)"allocated memory LtL"  ,naz*nax*nay*nz*nx*esize/2
  write(0,*)"allocated memory auxi2",naz*nax*nay*nz*nx*esize/2
  write(0,*)"n1,n2,n3,n4,n5,n6",n1,n2,n3,n4,n5,n6
  write(0,*)"nx,ny,nz,nsux,nsuy",nx,ny,nz,nsux,nsuy
  write(0,*)"d1,d2,d3,d4,d5,d6",d1,d2,d3,d4,d5,d6
  write(0,*)"o1,o2,o3,o4,o5,o6",o1,o2,o3,o4,o5,o6
  write(0,*)"marine,ratiox", marine,ratiox
  write(0,*)"marine,ratioy", marine,ratioy
!call stop_timer_num(a1)
!call start_timer_num(a2)

  do iy=1,ny
   write(0,*)"crossline #",iy
   LtL=0.
!call start_timer_num(a3)
   !$omp parallel do private(ix)
   do ix=1,nx
!     call omp_do (nx,ny,nz,nax,nay,naz,nax2,nay2,naz2,ix,iy,green,LtL,ratiox,ratioy,nsux,nsuy,a4,a5)
     call omp_do (nx,ny,nz,nax,nay,naz,nax2,nay2,naz2,ix,iy,green,LtL,ratiox,ratioy,nsux,nsuy,nry,nrx,nsy,nsx)
   enddo  
   !$omp end parallel do
!call stop_timer_num(a3)
   !---------------------------------------------------------------------------------------------------------------------
   call sep_open_lock("outtag")  !make sure only one process at a time writes
   call sseek("outtag",(iy-1)*esize*naz*nax*nay*nz*nx/2,0)
   if (sreed("outtag",auxi2,esize*naz*nax*nay*nz*nx/2) /= esize*naz*nax*nay*nz*nx/2) &
      call seperr("Error reading data")
   LtL= LtL+auxi2
   call sseek("outtag",0,0)
   call sseek("outtag",(iy-1)*esize*naz*nax*nay*nz*nx/2,0)
   if (srite("outtag",LtL,esize*naz*nax*nay*nz*nx/2) /= esize*naz*nax*nay*nz*nx/2) &
      call seperr("Error writing data")
   call sep_close_lock("outtag") !allow other processes to write
   !---------------------------------------------------------------------------------------------------------------------
  enddo

!call stop_timer_num(a2)
!----------------------------------------------------------------------------------------------------------------------
call date_and_time(values=time_array_1)
finish_timet = time_array_0 (3) * 86400.+ time_array_1 (5) * 3600  + time_array_1 (6) * 60 &
            + time_array_1 (7) + 0.001 * time_array_1 (8)
write(0,*)'total time', finish_timet - start_timet
!----------------------------------------------------------------------------------------------------------------------

deallocate(green,LtL,auxi2)
!call print_timers()

call sep_end_prog()

end program 
