!hessian compute
module hessian_mod
  use type_mod
  use global_defi_mod
  use ricker_mod
  use ctaper_mod
  use sep

  implicit none
  real   , dimension(:,:)      , pointer, private :: dgrid_header_rec_save
  real                                  , private :: dgrid_gxmin_save
  real                                  , private :: dgrid_gymin_save

  integer :: index_null_gxmin, index_null_gxmax
  real    :: null_gxmin, null_gxmax

    complex, dimension(:,:,:,:,:  ), allocatable :: greens1
    complex, dimension(:,:,:,:,:,:), allocatable :: hessian1
contains
!=======================================================================================
  integer function hess_init() result(stat)
    call greens_input_history()  !read greens history
    call param_init()             !initialize oneway parameters 
    call hess_output_history()    !output the oneway history
    call hess_array_alloc()       !allocate oneway memory space for the data
    !#    call hess_io_init()  
    call propagation_grid()
    
    call from_param("null_gxmin", null_gxmin, 0.) !min. geoph.
    call from_param("null_gxmax", null_gxmax, 0.) !max. geoph.

    index_null_gxmin = 1.5 + (null_gxmin - ggrid%ax_surf%o) / ggrid%ax_surf%d
    index_null_gxmax = 1.5 + (null_gxmax - ggrid%ax_surf%o) / ggrid%ax_surf%d

    write(*,*)'========================='
    write(0,*)"num shts.", dgrid%nsh!nsh num. of shots ......ntr num. of traces
    write(0,*)"num recv.", dgrid%ntr !nsh num. of shots ......ntr num. of traces
    write(*,*)'========================='
    
    stat = 0
  end function
!=======================================================================================
  integer function hess_clos() result(stat)
    call hess_array_dealloc()
    close(10)
    stat = 0
  end function
!=======================================================================================
  integer function hess_oper() result(stat)
    real                :: x_modl, y_modl, z_modl
    real                :: x_offd, y_offd, z_offd
    real                :: sx, sy, gx, gy
    integer             :: ix_modl, iy_modl, iz_modl
    integer             :: ix_offd, iy_offd, iz_offd
    integer             :: index_x_modl, index_y_modl, index_z_modl
    integer             :: index_x_offd, index_y_offd, index_z_offd
    integer             :: jxd_modl, jxu_modl
    integer             :: jyd_modl, jyu_modl
    integer             :: index_sx, index_sy                            !!.. relative index on Green's grid
    integer             :: jxd_offd, jyd_offd, jxu_offd, jyu_offd
    integer             :: jx_offd, jy_offd, jz_offd
    integer             :: igx, igy, index_gx, index_gy, isx,isy 
    complex             :: greens_rec_sum, greens_rec_mul, greens_so

    write(0,*)"num. freq.:", ggrid%aw
    write(*,*)'========================='

    open(unit=10, file=report_, status='unknown')
    
    hessian1=0.;
    nnode =omp_get_num_procs(); thread=omp_get_thread_num(); 
    call omp_set_num_threads(nodes); write(*,*) 'using', nodes,'of', nnode
!$OMP parallel do private(inode,iz_modl,z_modl,index_z_modl,ix_modl,x_modl,index_x_modl,z_offd,index_z_offd,ix_offd,x_offd,index_x_offd,iw,isx,sx,igx,gx,index_sx,index_gx,greens_so,greens_rec_sum, greens_rec_mul, greens1, hessian  ) schedule(dynamic,1)   
do iw=1, ggrid%aw%n, jw                                                            ! freq. count.
         inode=omp_get_thread_num();
         write(10,*)'iw: ', iw, 'of ',ggrid%aw%n ,' at inode: ', inode
         write(* ,*)'iw: ', iw, 'of ',ggrid%aw%n ,' at inode: ', inode
         !$OMP CRITICAL
         if (allocated(greens1)) then; write(*,*)size(greens1); deallocate(greens1); endif
         allocate(greens1(ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n)); greens1=0.
         call hess_io_greens(greens1, iw)
         !$OMP END CRITICAL
         do  iz_modl= 1,hgrid%az_modl%n                                                ! depth loop 
                z_modl = hgrid%az_modl%o + (iz_modl-1)*hgrid%az_modl%d
                index_z_modl = 1.5 + (z_modl - ggrid%az_modl%o)/ggrid%az_modl%d;
                if(inode==0) write(*,*)'iz: ', iz_modl, 'of ',hgrid%az_modl%n 
                !write(*,*)'iz: ', iz_modl, 'of ',hgrid%az_modl%n 
         do iy_modl= 1,hgrid%ay_modl%n                                                ! depth loop 
                y_modl = hgrid%ay_modl%o + (iy_modl-1)*hgrid%ay_modl%d
                index_y_modl = 1.5 + (y_modl - ggrid%ay_modl%o)/ggrid%ay_modl%d;
                ! write(*,*)'x_modl',x_modl,' at depth', z_modl, inode 
         do ix_modl= 1,hgrid%ax_modl%n                                               ! depth loop 
                x_modl = hgrid%ax_modl%o + (ix_modl-1)*hgrid%ax_modl%d
                index_x_modl = 1.5 + (x_modl - ggrid%ax_modl%o)/ggrid%ax_modl%d;
                !if (inode==1) write(*,*)iw,'x_modl:',x_modl,'depth', z_modl  
                !write(*,*)'x_modl',x_modl,' at depth', z_modl, inode 
                do  iz_offd= 2,hgrid%az_offd%n
                        z_offd = z_modl + hgrid%az_offd%o + (iz_offd-1)*hgrid%az_offd%d
                        if (z_offd < ggrid%az_modl%o .or. z_offd > ggrid%az_modl%o + (ggrid%az_modl%n-1)*ggrid%az_modl%d  ) cycle
                        index_z_offd = 1.5 + (z_offd - ggrid%az_modl%o)/ggrid%az_modl%d;
                do  iy_offd= 1,hgrid%ay_offd%n 
                        y_offd = y_modl + hgrid%ay_offd%o + (iy_offd-1)*hgrid%ay_offd%d
                        if (y_offd < ggrid%ay_modl%o .or. y_offd > ggrid%ay_modl%o + (ggrid%ay_modl%n-1)*ggrid%ay_modl%d  ) cycle
                        index_y_offd = 1.5 + (y_offd - ggrid%ay_modl%o)/ggrid%ay_modl%d;
                do  ix_offd= 2,hgrid%ax_offd%n-1 
                        x_offd = x_modl + hgrid%ax_offd%o + (ix_offd-1)*hgrid%ax_offd%d
                        if (x_offd < ggrid%ax_modl%o .or. x_offd > ggrid%ax_modl%o + (ggrid%ax_modl%n-1)*ggrid%ax_modl%d  ) cycle
                        index_x_offd = 1.5 + (x_offd - ggrid%ax_modl%o)/ggrid%ax_modl%d;
                        greens_so=0.
                        do isy= 1,dgrid%asy%n
                                sx = dgrid%asy%o + (isy-1)*dgrid%asy%d
                                index_sy = 1.5 + (sy - ggrid%ay_surf%o)/ggrid%ay_surf%d;
                        do isx= 1,dgrid%asx%n
                                sx = dgrid%asx%o + (isx-1)*dgrid%asx%d
                                index_sx = 1.5 + (sx - ggrid%ax_surf%o)/ggrid%ax_surf%d;
                                greens_so = greens1(index_sx, index_sy, index_x_modl, index_y_modl, index_z_modl)* & 
                                      conjg(greens1(index_sx, index_sy, index_x_offd, index_y_offd, index_z_offd)) 
                                greens_rec_sum = 0.
                        do igy= 1,dgrid%agy%n
                                gy = dgrid%agy%o + (igy-1)*dgrid%agy%d
                                index_gy = 1.5 + (gy - ggrid%ay_surf%o)/ggrid%ay_surf%d;
                        do igx= 1,dgrid%agx%n
                                gx = dgrid%agx%o + (igx-1)*dgrid%agx%d
                                index_gx = 1.5 + (gx - ggrid%ax_surf%o)/ggrid%ax_surf%d;
                                greens_rec_mul = greens1(index_gx, index_gy, index_x_modl, index_y_modl, index_z_modl)* & 
                                           conjg(greens1(index_gx, index_gy, index_x_offd, index_y_offd, index_z_offd)) 
                                greens_rec_sum = greens_rec_sum + real(greens_so*greens_rec_mul)
                        enddo
                        enddo !recv

                        !$OMP CRITICAL
                        hessian1(ix_offd, iy_offd, iz_offd, ix_modl, iy_modl, iz_modl )= & 
                        hessian1(ix_offd, iy_offd, iz_offd, ix_modl, iy_modl, iz_modl )+real(greens_rec_sum)
                        !$OMP END CRITICAL
                        enddo 
                        enddo !sou
                enddo
                enddo
                enddo         !offd
        enddo
        enddo
        enddo                 !modl  
        !deallocate(greens1)
        !hessian1(1,:,:,:,:,:)=0.0; hessian1(hgrid%ax_offd%n,:,:,:,:,:)=0.0; 
        !hessian1(:,:,1,:,:,:)=0.0; !hessian1(:,:,hgrid%az_offd%n,:,:,:)=0.0; 
        if (inode==0) then
               if (hgrid%ax_offd%n >1) then; hessian1(2,:,:,:,:,:)=hessian1(2,:,:,:,:,:)/1.5; hessian1(hgrid%ax_offd%n-1,:,:,:,:,:)= hessian1(hgrid%ax_offd%n-1,:,:,:,:,:)/1.5; endif
               if (hgrid%az_offd%n >1) then; hessian1(:,:,2,:,:,:)=hessian1(:,:,2,:,:,:)/1.5; endif !hessian1(:,:,hgrid%az_offd%n-1,:,:,:)= hessian1(:,:,hgrid%az_offd%n-1,:,:,:)/2;endif 
               call hess_io_hess(hessian1)
        endif
enddo
!$OMP end parallel do 
        
    call hess_io_hess(hessian1)
    stat = 0
  end function

!=========================================================================================
  subroutine hess_array_alloc()
    integer :: memory
    allocate(hessian1(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n))
    allocate(greens(ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n))
    allocate(wavelet(1))

    memory = 4 * size(hessian) + 8 * (size(greens) + size(wavelet)) ! 4*size(hessian)+ 4*size(source)+4*size(rec.)+size(wavelet)...men
    memory = memory / 1024 / 1024
    write(0,*)"Memory allocated: ", memory, "Mb"

  end subroutine
!=================================================================================================
subroutine hess_array_dealloc()
  if(associated(hessian)) deallocate(hessian)
  if(associated(wavelet)) deallocate(wavelet)
  if(associated(greens )) deallocate(greens)
  end subroutine
!==========================================================================================
  subroutine param_init()
    real    :: xmin_modl_default, xmax_modl_default
    real    :: ymin_modl_default, ymax_modl_default
    real    :: zmin_modl_default, zmax_modl_default
    real    :: surface_xmin_default, surface_xmax_default
    real    :: surface_ymin_default, surface_ymax_default

    call from_param("nsx", dgrid%asx%n, 1)   !num. sources inline
    call from_param("dsx", dgrid%asx%d, 1.)   !samp. sources inline 
    call from_param("osx", dgrid%asx%o, 0.)   !orig. sources inline
    call from_param("nsy", dgrid%asy%n, 1)   !              x-line
    call from_param("dsy", dgrid%asy%d, 1.)   !
    call from_param("osy", dgrid%asy%o, 0.)   !
    dgrid%nsh = dgrid%asx%n * dgrid%asy%n     !num  of sources..
    call from_param("ngx", dgrid%agx%n, 1)
    call from_param("dgx", dgrid%agx%d, 1.)
    call from_param("ogx", dgrid%agx%o, 0.)
    call from_param("ngy", dgrid%agy%n, 1)
    call from_param("dgy", dgrid%agy%d, 1.)
    call from_param("ogy", dgrid%agy%o, 0.)
    dgrid%ntr = dgrid%agx%n * dgrid%agy%n   !num. of traces/recievers
    write(*,*)'sx:n_d_o',dgrid%asx
    write(*,*)'sy:n_d_o',dgrid%asy
    write(*,*)'gx:n_d_o',dgrid%agx
    write(*,*)'gy:n_d_o',dgrid%agy
    write(*,*)'========================='
    !off-diags
    call from_param("nx_offd", hgrid%ax_offd%n, 1 )               !num. of off-diags.
    call from_param("dx_offd", hgrid%ax_offd%d, ggrid%ax_modl%d)  !samp 
    call from_param("ox_offd", hgrid%ax_offd%o, 0.)               !orig 
   
    call from_param("ny_offd", hgrid%ay_offd%n, 1 )
    call from_param("dy_offd", hgrid%ay_offd%d, ggrid%ay_modl%d)
    call from_param("oy_offd", hgrid%ay_offd%o, 0.)

    call from_param("nz_offd", hgrid%az_offd%n, 1 )
    call from_param("dz_offd", hgrid%az_offd%d, ggrid%az_modl%d)
    call from_param("oz_offd", hgrid%az_offd%o, 0.)
    write(*,*)'ax_offd:n_d_o',hgrid%ax_offd
    write(*,*)'ay_offd:n_d_o',hgrid%ay_offd
    write(*,*)'az_offd:n_d_o',hgrid%az_offd
    write(*,*)'========================='
    !============================================================
    !sub-offset
    call from_param("nhx_modl", hgrid%ahx_modl%n, 1 )
    call from_param("dhx_modl", hgrid%ahx_modl%d, ggrid%ax_surf%d)
    call from_param("ohx_modl", hgrid%ahx_modl%o, 0.)
    call from_param("nhy_modl", hgrid%ahy_modl%n, 1 )
    call from_param("dhy_modl", hgrid%ahy_modl%d, ggrid%ay_surf%d)
    call from_param("ohy_modl", hgrid%ahy_modl%o, 0.)
    !sub-offset off-diags
    call from_param("nhx_offd", hgrid%ahx_offd%n, 1 )
    call from_param("dhx_offd", hgrid%ahx_offd%d, ggrid%ax_surf%d)
    call from_param("ohx_offd", hgrid%ahx_offd%o, 0.)
    call from_param("nhy_offd", hgrid%ahy_offd%n, 1 )
    call from_param("dhy_offd", hgrid%ahy_offd%d, ggrid%ay_surf%d)
    call from_param("ohy_offd", hgrid%ahy_offd%o, 0.)
    !============================================================
    !set default hessian filter sizes
    xmin_modl_default = ggrid%ax_modl%o
    xmax_modl_default = ggrid%ax_modl%o + (ggrid%ax_modl%n-1)*ggrid%ax_modl%d
    ymin_modl_default = ggrid%ay_modl%o
    ymax_modl_default = ggrid%ay_modl%o + (ggrid%ay_modl%n-1)*ggrid%ay_modl%d
    zmin_modl_default = ggrid%az_modl%o
    zmax_modl_default = ggrid%az_modl%o + (ggrid%az_modl%n-1)*ggrid%az_modl%d

    call from_param("xmin_modl", hgrid%xmin_modl, xmin_modl_default)
    call from_param("xmax_modl", hgrid%xmax_modl, xmax_modl_default)
    call from_param("ymin_modl", hgrid%ymin_modl, ymin_modl_default)
    call from_param("ymax_modl", hgrid%ymax_modl, ymax_modl_default)
    call from_param("zmin_modl", hgrid%zmin_modl, zmin_modl_default)
    call from_param("zmax_modl", hgrid%zmax_modl, zmax_modl_default)
    !============================================================
    !define the local hessian grid...around image point
    hgrid%ax_modl%o = hgrid%xmin_modl
    hgrid%ax_modl%d = ggrid%ax_modl%d
    hgrid%ax_modl%n = 1.5 + (hgrid%xmax_modl-hgrid%xmin_modl)/hgrid%ax_modl%d

    hgrid%ay_modl%o = hgrid%ymin_modl
    hgrid%ay_modl%d = ggrid%ay_modl%d
    hgrid%ay_modl%n = 1.5 + (hgrid%ymax_modl-hgrid%ymin_modl)/hgrid%ay_modl%d

    hgrid%az_modl%o = hgrid%zmin_modl
    hgrid%az_modl%d = ggrid%az_modl%d
    hgrid%az_modl%n = 1.5 + (hgrid%zmax_modl-hgrid%zmin_modl)/hgrid%az_modl%d

    write(*,*)'ax_modl:n_d_o',hgrid%ax_modl
    write(*,*)'ay_modl:n_d_o',hgrid%ay_modl
    write(*,*)'az_modl:n_d_o',hgrid%az_modl
    write(*,*)'========================='

    !============================================================
    report_ = "report"    ; 
    call from_param("jw"          , jw           , 1      )
    call from_param("taper"       , taper        , .true. )
    call from_param("taper_pct"   , taper_pct    , 0.1    )
    call from_param("taper_type"  , taper_type   , 2      )
    call from_param("nodes     "  , nodes        , 1      )
  end subroutine
!========================================================================================== 
  subroutine greens_input_history()
    greens_= "greens"; 

    call auxin(greens_)
    call from_aux(greens_, "n1", ggrid%ax_surf%n)  
    call from_aux(greens_, "d1", ggrid%ax_surf%d) 
    call from_aux(greens_, "o1", ggrid%ax_surf%o) 

    call from_aux(greens_, "n2", ggrid%ay_surf%n)  
    call from_aux(greens_, "d2", ggrid%ay_surf%d) 
    call from_aux(greens_, "o2", ggrid%ay_surf%o) 

    call from_aux(greens_, "n3", ggrid%ax_modl%n)  
    call from_aux(greens_, "d3", ggrid%ax_modl%d) 
    call from_aux(greens_, "o3", ggrid%ax_modl%o) 

    call from_aux(greens_, "n4", ggrid%ay_modl%n)  
    call from_aux(greens_, "d4", ggrid%ay_modl%d) 
    call from_aux(greens_, "o4", ggrid%ay_modl%o) 

    call from_aux(greens_, "n5", ggrid%az_modl%n)  
    call from_aux(greens_, "d5", ggrid%az_modl%d) 
    call from_aux(greens_, "o5", ggrid%az_modl%o) 

    call from_aux(greens_, "n6", ggrid%aw%n)  
    call from_aux(greens_, "d6", ggrid%aw%d) 
    call from_aux(greens_, "o6", ggrid%aw%o) 

    write(*,*)'x_surf:n_d_o',ggrid%ax_surf
    write(*,*)'y_surf:n_d_o',ggrid%ay_surf
    write(*,*)'x_modl:n_d_o',ggrid%ax_modl
    write(*,*)'y_modl:n_d_o',ggrid%ay_modl
    write(*,*)'z_modl:n_d_o',ggrid%az_modl
    write(*,*)'========================='
  end subroutine
!==========================================================================================
  subroutine hess_output_history()
    !!!!!Outputs the Hessain details into the history file...
    hessian_= "hessian"; 

    call auxinout(hessian_)
    !sub-offset off-diagonals ==========================
    call auxputch("n1", "i", hgrid%ahx_offd%n, hessian_) !sub. offset off-diag num 
    call auxputch("d1", "r", hgrid%ahx_offd%d, hessian_) !                     samp  
    call auxputch("o1", "r", hgrid%ahx_offd%o, hessian_) !                     orig.
    call auxputch("label1", "s", "hx_offd", hessian_)

    call auxputch("n2", "i", hgrid%ahy_offd%n, hessian_)
    call auxputch("d2", "r", hgrid%ahy_offd%d, hessian_)
    call auxputch("o2", "r", hgrid%ahy_offd%o, hessian_)
    call auxputch("label2", "s", "hy_offd", hessian_)
    !===================================================
    !zero offset off-diagonals
    call auxputch("n3", "i", hgrid%ax_offd%n, hessian_) !off-diag num.
    call auxputch("d3", "r", hgrid%ax_offd%d, hessian_) !         samp.
    call auxputch("o3", "r", hgrid%ax_offd%o, hessian_) !         orig.
    call auxputch("label3", "s", "x_offd", hessian_)

    call auxputch("n4", "i", hgrid%ay_offd%n, hessian_)
    call auxputch("d4", "r", hgrid%ay_offd%d, hessian_)
    call auxputch("o4", "r", hgrid%ay_offd%o, hessian_)
    call auxputch("label4", "s", "y_offd", hessian_)

    call auxputch("n5", "i", hgrid%az_offd%n, hessian_)
    call auxputch("d5", "r", hgrid%az_offd%d, hessian_)
    call auxputch("o5", "r", hgrid%az_offd%o, hessian_)
    call auxputch("label5", "s", "z_offd", hessian_)
    !===================================================
    !sub-offset modl. points
    call auxputch("n6", "i", hgrid%ahx_modl%n, hessian_) !sub-offset model point num.
    call auxputch("d6", "r", hgrid%ahx_modl%d, hessian_) !                       samp.
    call auxputch("o6", "r", hgrid%ahx_modl%o, hessian_) !                       orig.
    call auxputch("label6", "s", "hx_modl", hessian_)

    call auxputch("n7", "i", hgrid%ahy_modl%n, hessian_)
    call auxputch("d7", "r", hgrid%ahy_modl%d, hessian_)
    call auxputch("o7", "r", hgrid%ahy_modl%o, hessian_)
    call auxputch("label7", "s", "hy_modl", hessian_)
    !===================================================
    !zero-offset modl. points...actual model points
    call auxputch("n8", "i", hgrid%ax_modl%n, hessian_) !model point num.
    call auxputch("d8", "r", hgrid%ax_modl%d, hessian_) !            samp.
    call auxputch("o8", "r", hgrid%ax_modl%o, hessian_) !            orig. 
    call auxputch("label8", "s", "x_modl", hessian_)

    call auxputch("n9", "i", hgrid%ay_modl%n, hessian_)
    call auxputch("d9", "r", hgrid%ay_modl%d, hessian_)
    call auxputch("o9", "r", hgrid%ay_modl%o, hessian_)
    call auxputch("label9", "s", "y_modl", hessian_)

    call auxputch("n10", "i", hgrid%az_modl%n, hessian_)
    call auxputch("d10", "r", hgrid%az_modl%d, hessian_)
    call auxputch("o10", "r", hgrid%az_modl%o, hessian_)
    call auxputch("label10", "s", "z_modl", hessian_)
    !===================================================
  end subroutine

!====================================================================================================================
!initialize the full hessian to zeros.
subroutine hess_io_init()
    real, allocatable, dimension(:) :: ctmp
    integer                         :: i, nloop
    allocate(ctmp(hgrid%ahx_modl%n*hgrid%ahy_modl%n*hgrid%ax_modl%n*hgrid%ay_modl%n*hgrid%az_modl%n))
    ctmp = 0.
          nloop = hgrid%ahx_offd%n*hgrid%ahy_offd%n*hgrid%ax_offd%n*hgrid%ay_offd%n*hgrid%az_offd%n
    do i=1, nloop
       call sep_write(ctmp, hessian_, esize=4)
    end do
    deallocate(ctmp)
  end subroutine
!====================================================================================================================
  subroutine hess_io_greens(greens, index_ws)
    integer                                   :: index_ws
    integer, dimension(6)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    complex, dimension(:,:,:,:,:), allocatable:: greens

    n_dim_cube = 6
    n_cube = (/ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n, ggrid%aw%n/)
    j_wind = (/1, 1, 1, 1, 1, 1/)
    n_wind = (/ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n, 1/)
    f_wind = (/0, 0, 0, 0, 0, index_ws-1/)

    greens = 0.
    !write(0,*)'reading greens frequency: ', index_ws, 'of ', ggrid%aw%n
    err = sreed_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,greens)
    if (err /= 0) call seperr("data_io_greens: error writing greens")

  end subroutine
!===================================================
  subroutine hess_io_wavelet(wavelet, index_ws)
    complex, dimension(:), pointer :: wavelet
    integer                        :: index_ws
    integer, dimension(1)          :: n_cube, j_wind, n_wind, f_wind
    integer                        :: n_dim_cube, err

    n_dim_cube = 1
    n_cube = dgrid%aw%n
    j_wind = 1 
    n_wind = 1
    f_wind = index_ws-1

    wavelet = 0.
    err = sreed_window(wavelet_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,wavelet)
    if (err .lt. 0) call seperr("io_wavelet: trouble reading source wavelet")
  end subroutine

!#####################################################
 subroutine hess_io_hess(hessian)
    integer                                     :: index_zs
    integer, dimension(10)                      :: n_cube, j_wind, n_wind, f_wind
    integer                                     :: n_dim_cube, err
    complex, dimension(:,:,:,:,:,:), allocatable:: hessian

    n_dim_cube = 10
    n_cube = (/hgrid%ahx_offd%n, hgrid%ahy_offd%n, hgrid%ax_offd%n, hgrid%ay_offd%n, hgrid%az_offd%n, &
               hgrid%ahx_modl%n, hgrid%ahy_modl%n, hgrid%ax_modl%n, hgrid%ay_modl%n, hgrid%az_modl%n/)
    j_wind = (/1, 1, 1, 1, 1, 1, 1, 1, 1, 1/)
    n_wind = (/hgrid%ahx_offd%n, hgrid%ahy_offd%n, hgrid%ax_offd%n, hgrid%ay_offd%n, hgrid%az_offd%n, & 
               hgrid%ahx_modl%n, hgrid%ahy_modl%n, hgrid%ax_modl%n, hgrid%ay_modl%n, hgrid%az_modl%n /)
    f_wind = (/0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)

!   write(0,*)'writing hessian depth: ', index_zs, 'of ', hgrid%az_modl%n
    err = srite_window(hessian_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,real(hessian))
    if (err /= 0) call seperr("hess_io_hess: error writing hessian")

  end subroutine

!#####################################################
subroutine propagation_grid()
    real    :: pgrid_xmin, pgrid_xmax
    real    :: pgrid_ymin, pgrid_ymax

    pgrid%ax%d  = vgrid%ax%d; pgrid%ay%d  = vgrid%ay%d; pgrid%az%d  = vgrid%az%d
 
    pgrid%gxmin = ggrid%surface_xmin; pgrid%gxmax = ggrid%surface_xmax
    pgrid%gymin = ggrid%surface_ymin; pgrid%gymax = ggrid%surface_ymax

    !################################################
    pgrid_xmin = pgrid%gxmin; pgrid_xmax = pgrid%gxmax
    pgrid_ymin = pgrid%gymin; pgrid_ymax = pgrid%gymax

    !!.. add the padding on both x and y direction
    pgrid_xmin = pgrid_xmin - pgrid%nxpad_beg*pgrid%ax%d
    pgrid_xmax = pgrid_xmax + pgrid%nxpad_end*pgrid%ax%d

    pgrid_ymin = pgrid_ymin - pgrid%nypad_beg*pgrid%ay%d
    pgrid_ymax = pgrid_ymax + pgrid%nypad_end*pgrid%ay%d

    !!.. make sure the propagation grid is within the velocity model 
    pgrid_xmin = max(pgrid_xmin, vgrid%ax%o)
    pgrid_xmax = min(pgrid_xmax, vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d)

    pgrid_ymin = max(pgrid_ymin, vgrid%ay%o)
    pgrid_ymax = min(pgrid_ymax, vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d)
    !################################################
    !!.. propagation grid 
    pgrid%ax%o = pgrid_xmin;     pgrid%ay%o = pgrid_ymin
    pgrid%ax%n = 1.5 + (pgrid_xmax - pgrid_xmin)/pgrid%ax%d
    pgrid%ay%n = 1.5 + (pgrid_ymax - pgrid_ymin)/pgrid%ay%d
    pgrid%az%n = vgrid%az%n

    !!.. the first and last index of surface x and y on the propagation grid (for output use)
    ggrid%index_surpro_xmin = 1.5 + (ggrid%surface_xmin - pgrid%ax%o) / pgrid%ax%d
    ggrid%index_surpro_xmax = 1.5 + (ggrid%surface_xmax - pgrid%ax%o) / pgrid%ax%d

    ggrid%index_surpro_ymin = 1.5 + (ggrid%surface_ymin - pgrid%ay%o) / pgrid%ay%d
    ggrid%index_surpro_ymax = 1.5 + (ggrid%surface_ymax - pgrid%ay%o) / pgrid%ay%d

    !!....first and last index of geoph. and source
    pgrid%index_datpro_gxmin = 1.5 + (pgrid%gxmin  - pgrid_xmin ) / pgrid%ax%d  
    pgrid%index_datpro_gymin = 1.5 + (pgrid%gymin  - pgrid_ymin ) / pgrid%ay%d  

    pgrid%index_datpro_gxmax = 1.5 + (pgrid%gxmax  - pgrid_xmin ) / pgrid%ax%d  
    pgrid%index_datpro_gymax = 1.5 + (pgrid%gymax  - pgrid_ymin ) / pgrid%ay%d  

    !!....first and last index of propagation grid in the velocity model
    pgrid%index_xmin = 1.5 + (pgrid_xmin  - vgrid%ax%o ) / pgrid%ax%d 
    pgrid%index_ymin = 1.5 + (pgrid_ymin  - vgrid%ay%o ) / pgrid%ay%d  

    pgrid%index_xmax = 1.5 + (pgrid_xmax  - vgrid%ax%o ) / pgrid%ax%d  
    pgrid%index_ymax = 1.5 + (pgrid_ymax  - vgrid%ay%o ) / pgrid%ay%d  

  end subroutine
!######################################################################
  end module

