!Hessian module...modified from Yaxun Tang's original code
module wei_hessian_mod
  use wei_type_mod
  use sep
  implicit none

  type(dat_grid)                         ,private :: dgrid !data grid
  type(gre_grid)                         ,private :: ggrid !greens function grid
  type(hes_grid)                         ,private :: hgrid !hessian grid

  real   , dimension(:,:,:,:,:), pointer, private :: hessian 
  complex, dimension(:,:,:,:,:), pointer, private :: greens
  complex, dimension(:)        , pointer, private :: wavelet

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

  integer, external                     , private :: srite_window, sreed_window, getch

  character(len=256)                    , private :: greens_ ="greens"
  character(len=256)                    , private :: hessian_="hessian"
  character(len=256)                    , private :: wavelet_="wavelet"
  character(len=256)                    , private :: hdr_sou_="hdr_sou"
  character(len=256)                    , private :: hdr_rec_="hdr_rec"

  logical                               , private :: upward, offset, sou_sig

  logical, parameter                    , private :: T=.true., F=.false.

integer :: index_null_gxmin, index_null_gxmax
real    :: null_gxmin, null_gxmax
contains
!=======================================================================================
  integer function wei_hessian_init() result(stat)
    call wei_hessian_input_history()  !read greens history
    call wei_hessian_param_init()     !initialize hessian parameters 
    call wei_hessian_output_history() !output the hessian history
    call wei_hessian_array_alloc()    !allocate hessian memory space for the hessian
    call wei_hessian_io_init()        !initialize hessian input/output params

    if(dgrid%hdr_sou) call wei_hessian_io_head_init("hsou")
    if(dgrid%hdr_rec) call wei_hessian_io_head_init("hrec")

    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(0,*)"dgrid%hdr_sou=", dgrid%hdr_sou
    write(0,*)"dgrid%hdr_rec=", dgrid%hdr_rec
    write(0,*)"nsh, ntr", dgrid%nsh, dgrid%ntr !nsh num. of shots ......ntr num. of traces
    stat = 0
  end function
!=======================================================================================
  integer function wei_hessian_clos() result(stat)
    call wei_hessian_array_dealloc()
    stat = 0
  end function
!=======================================================================================
  integer function wei_hessian_oper() result(stat)
    integer   :: ix_modl, iy_modl, iz_modl
    integer   :: ihx_modl, ihy_modl, iws
    integer   :: jhx_modl, jhy_modl
    integer   :: jxd_modl, jxu_modl
    integer   :: jyd_modl, jyu_modl
    integer   :: index_x_modl, index_y_modl, index_z_modl
    real      :: x_modl, y_modl, z_modl, hx_modl, hy_modl

    write(0,*)"hgrid%ax_modl", hgrid%ax_modl
    write(0,*)"hgrid%ay_modl", hgrid%ay_modl
    write(0,*)"hgrid%az_modl", hgrid%az_modl
    write(0,*)"ggrid%aw", ggrid%aw

    wavelet = 1.

    do iws=1, ggrid%aw%n                                                            ! freq. count.
       call wei_hessian_io_greens_test(greens, iws)
       write(0,*)"iws----->, w", iws , ggrid%aw%o + (iws-1)*ggrid%aw%d
       if (sou_sig) call wei_hessian_io_wavelet(wavelet, iws, 1)                    ! read freq. bloack of source wavelet
    !  write(0,*)"hollaaaa"; STOP 
       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
write(0,*)"iz_modl", iz_modl, z_modl
          do iy_modl=1, hgrid%ay_modl%n                                             ! x-line 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  
             do ix_modl=1, hgrid%ax_modl%n                                          ! in-line 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     
!write(0,*)"z_modl, xmodl, ymodli: ", z_modl, x_modl,y_modl
                do ihy_modl=1, hgrid%ahy_modl%n                                     ! x-line sub-surface offset loop
                   hy_modl  = hgrid%ahy_modl%o + (ihy_modl-1)*hgrid%ahy_modl%d      ! x-line sub-offset
                   jhy_modl = nint(hy_modl / ggrid%ay_surf%d)
                   !jyd_modl = max(1, min(iy_modl - jhy_modl, ggrid%ay_surf%n))
                   !jyu_modl = max(1, min(iy_modl + jhy_modl, ggrid%ay_surf%n))
                   jyd_modl = index_y_modl - jhy_modl                               !-ve offset
                   jyu_modl = index_y_modl + jhy_modl                               !+ve offset
                   if (jyd_modl < 1 .or. jyd_modl > ggrid%ay_surf%n .or. &          ! dont go out of surface grid.
                       jyu_modl < 1 .or. jyu_modl > ggrid%ay_surf%n) cycle

                   do ihx_modl=1, hgrid%ahx_modl%n                                  ! in-line sub-surface offset loop
                      hx_modl  = hgrid%ahx_modl%o + (ihx_modl-1)*hgrid%ahx_modl%d
                      jhx_modl = nint(hx_modl / ggrid%ax_surf%d)
                      !jxd_modl = max(1, min(ix_modl - jhx_modl, ggrid%ax_surf%n))
                      !jxu_modl = max(1, min(ix_modl + jhx_modl, ggrid%ax_surf%n))
                      jxd_modl = index_x_modl - jhx_modl                            !-ve offset
                      jxu_modl = index_x_modl + jhx_modl                            !+ve offset
                      if (jxd_modl < 1 .or. jxd_modl > ggrid%ax_surf%n .or. &       !dont go out of the grid
                          jxu_modl < 1 .or. jxu_modl > ggrid%ax_surf%n) cycle
                                                                                    ! read the greens of model point and relative model point...etc 
                      call wei_hessian_io_hessian("reed", hessian, ihx_modl, ihy_modl, ix_modl, iy_modl, iz_modl)
                                                                                    ! compute the hessian relative to the specified model point
                      call wei_hessian_compute(hessian, greens, T, index_x_modl, index_y_modl, index_z_modl, &
                                               jxd_modl, jyd_modl, jxu_modl, jyu_modl)
                                                                                    ! write out the computed hessian param.
                      call wei_hessian_io_hessian("rite", hessian, ihx_modl, ihy_modl, ix_modl, iy_modl, iz_modl)
                   end do
                end do
             end do
          end do
       end do
    end do
    stat = 0
  end function

!=========================================================================================
  subroutine wei_hessian_compute(hessian, greens, add, index_x_modl, index_y_modl, index_z_modl, &
                                 jxd_modl, jyd_modl, jxu_modl, jyu_modl)
    real   , dimension(:,:,:,:,:), pointer :: hessian
    complex, dimension(:,:,:,:,:), pointer :: greens
    logical                                :: add
    integer                                :: index_x_modl, index_y_modl, index_z_modl
    integer                                :: jxd_modl, jyd_modl, jxu_modl, jyu_modl
    integer                                :: ish, index_sx, index_sy
    complex :: greens_sou, greens_rec_sum
    integer :: igx, igy, index_gx, index_gy

    do ish=1, dgrid%nsh

       if(dgrid%hdr_sou) call wei_hessian_io_head_oper("hsou", ish)
       if(dgrid%hdr_rec) call wei_hessian_io_head_oper("hrec", ish)
       call wei_hessian_geometry(ish)
       index_sx = hgrid%index_datgre_sx
       index_sy = hgrid%index_datgre_sy
       call wei_hessian_one_step(hessian, greens, add, index_sx, index_sy, &
                                 index_x_modl, index_y_modl, index_z_modl, &
                                 jxd_modl, jyd_modl, jxu_modl, jyu_modl)
    end do
  end subroutine
!==========================================================================================
! compute Hessian for one model point and one shot
  subroutine wei_hessian_one_step(hessian, greens, add, index_sx, index_sy, &
                                  index_x_modl, index_y_modl, index_z_modl, &
                                  jxd_modl, jyd_modl, jxu_modl, jyu_modl)

    real   , dimension(:,:,:,:,:), pointer :: hessian
    complex, dimension(:,:,:,:,:), pointer :: greens
    logical                                :: add
    integer                                :: index_sx, index_sy                            !!.. relative index on Green's grid
    integer                                :: index_x_modl, index_y_modl, index_z_modl      !!.. relative index on ggrid
    integer                                :: jxd_modl, jyd_modl, jxu_modl, jyu_modl        !!.. relative index with rsp to index_x, y
    integer                                :: ix_offd, iy_offd, iz_offd, ihx_offd, ihy_offd
    integer                                :: jhx_offd, jhy_offd
    integer                                :: jxd_offd, jyd_offd, jxu_offd, jyu_offd
    integer                                :: jx_offd, jy_offd, jz_offd
    integer                                :: index_x_offd, index_y_offd, index_z_offd
    integer                                :: itr, igx, igy, index_gx, index_gy
    integer                                :: ndim1, ndim2, ndim3, ndim4, ndim5
    real                                   :: x_offd, y_offd, z_offd, hx_offd, hy_offd
    real                                   :: hessian_one_shot
    complex                                :: greens_rec_sum, greens_rec_mul, greens_sou

    ndim1 = size(greens, 1)
    ndim2 = size(greens, 2)
    ndim3 = size(greens, 3)
    ndim4 = size(greens, 4)
    ndim5 = size(greens, 5)
    if (.not. add) then
       hessian = 0.
    end if

    do iz_offd=1, hgrid%az_offd%n                                             !loop over offdiag. in depth
       z_offd = hgrid%az_offd%o + (iz_offd-1)*hgrid%az_offd%d                 !compute actual offdiag depth 
 !       write(0,*)z_offd
       jz_offd = nint(z_offd / ggrid%az_modl%d)                               !compute index of depth on grid 
       index_z_offd = index_z_modl + jz_offd                                  !relative index to our image point....
       if (index_z_offd < 1 .or. index_z_offd > ndim5) cycle                  !must be within our hessian depth window

       do iy_offd=1, hgrid%ay_offd%n                                          !xline 
          y_offd = hgrid%ay_offd%o + (iy_offd-1)*hgrid%ay_offd%d
          jy_offd = nint(y_offd / ggrid%ay_modl%d)

          do ix_offd=1, hgrid%ax_offd%n                                       !inline 
             x_offd = hgrid%ax_offd%o + (ix_offd-1)*hgrid%ax_offd%d
             jx_offd = nint(x_offd / ggrid%ax_modl%d)

             do ihy_offd=1, hgrid%ahy_offd%n                                  !xline sub-offset 
                hy_offd  = hgrid%ahy_offd%o + (ihy_offd-1)*hgrid%ahy_offd%d
                jhy_offd = nint(hy_offd / ggrid%ay_surf%d)

                jyd_offd = jyd_modl + jy_offd - jhy_offd     !!.. jyd_offd: the index for y'- hy' = y + yt - (hy + hyt) = y - hy + yt - hyt
                jyu_offd = jyu_modl + jy_offd + jhy_offd     !!.. jyu_offd: the index for y'+ hy' = y + yt + (hy + hyt) = y + hy + yt + hyt

                if (jyd_offd < 1 .or. jyd_offd > ndim4 .or. &
                    jyu_offd < 1 .or. jyu_offd > ndim4) cycle

                do ihx_offd=1, hgrid%ahx_offd%n                                !inline sub-offset
                   hx_offd = hgrid%ahx_offd%o + (ihx_offd-1)*hgrid%ahx_offd%d
                   jhx_offd = nint(hx_offd / ggrid%ax_surf%d)

                   jxd_offd = jxd_modl + jx_offd - jhx_offd   !!.. jxd_offd: the index for x'- hx' = x + xt - (hx + hxt) = x - hx + xt - hxt
                   jxu_offd = jxu_modl + jx_offd + jhx_offd   !!.. jxu_offd: the index for x'+ hx' = x + xt + (hx + hxt) = x + hx + xt + hxt

                   if (jxd_offd < 1 .or. jxd_offd > ndim3 .or. &
                       jxu_offd < 1 .or. jxu_offd > ndim3) cycle

                   greens_rec_sum = 0.

                   do itr=1, dgrid%ntr                           !loop over receivers/traces
                      index_gx = hgrid%index_datgre_rec(1, itr)
                      index_gy = hgrid%index_datgre_rec(2, itr)
                      if (index_gx < 1 .or. index_gx > ndim1) cycle
                      if (index_gy < 1 .or. index_gy > ndim2) cycle
                      greens_rec_mul = greens(index_gx, index_gy, jxu_modl, jyu_modl, index_z_modl) &
                                     * conjg(greens(index_gx, index_gy, jxu_offd, jyu_offd, index_z_offd))
                      greens_rec_sum = greens_rec_sum + greens_rec_mul
                   end do

                   greens_sou = greens(index_sx, index_sy, jxd_modl, jyd_modl, index_z_modl) &
                              *conjg(greens(index_sx, index_sy, jxd_offd, jyd_offd, index_z_offd))

                   hessian_one_shot = real(greens_sou * greens_rec_sum) * cabs(wavelet(1))

                   hessian(ihx_offd, ihy_offd, ix_offd, iy_offd, iz_offd) = &
                   hessian(ihx_offd, ihy_offd, ix_offd, iy_offd, iz_offd) + hessian_one_shot
!write(0,*)"-----------------------------------------"
!write(0,*)"jxd_offd, jyd_offd, jxu_offd, jyu_offd", jxd_offd, jyd_offd, jxu_offd, jyu_offd
!write(0,*)"jxd_modl, jyd_modl, jxu_modl, jyu_modl", jxd_modl, jyd_modl, jxu_modl, jyu_modl
                end do
             end do
          end do
       end do
    end do
  end subroutine

!==========================================================================================
! compute the recording geometry, if model points are out of recording geometry, cycle
!==========================================================================================
  subroutine wei_hessian_geometry(ish)
    integer                                 :: ish, itr, igx, igy
    integer                                 :: ix_modl, iy_modl, iz_modl
    integer                                 :: ishx, ishy, index_gx, index_gy
    real                                    :: gx, gy 

    if (dgrid%hdr_sou) then !!.. if source header presented
       dgrid%sx = dgrid%header_sou(1)
       dgrid%sy = dgrid%header_sou(2)
    else
       ishx = mod(ish-1, dgrid%asx%n)+1
       ishy = int(1. + real(ish-1)/real(dgrid%asx%n))
       dgrid%sx = dgrid%asx%o + (ishx-1)*dgrid%asx%d
       dgrid%sy = dgrid%asy%o + (ishy-1)*dgrid%asy%d
    end if


    !!.. if the input receiver coordinate is in offset space
    if (offset) then
       if (dgrid%hdr_rec) then !!.. not working properly
          dgrid%header_rec(1,:) = dgrid_header_rec_save(1,:) + dgrid%sx
          dgrid%header_rec(2,:) = dgrid_header_rec_save(2,:) + dgrid%sy
       else
          dgrid%agx%o = dgrid_gxmin_save + dgrid%sx
          dgrid%agy%o = dgrid_gymin_save + dgrid%sy
       end if
    end if

    !!.. determine the index of source location on the Green's grid
    hgrid%index_datgre_sx = 1.5 + (dgrid%sx - ggrid%ax_surf%o) / ggrid%ax_surf%d
    hgrid%index_datgre_sy = 1.5 + (dgrid%sy - ggrid%ay_surf%o) / ggrid%ay_surf%d

    !!.. determine the index of receivers on the Green's grid
    !!.. for the case of unsorted data, all indices should be calculated to
    !!.. map the receivers onto the propagation grid
    if (dgrid%hdr_rec) then
       do itr=1, dgrid%ntr
          hgrid%index_datgre_rec(1, itr) = 1.5 + (dgrid%header_rec(1, itr) - ggrid%ax_surf%o) / ggrid%ax_surf%d
          hgrid%index_datgre_rec(2, itr) = 1.5 + (dgrid%header_rec(2, itr) - ggrid%ay_surf%o) / ggrid%ay_surf%d
       end do
    !!.. for the case of regular data, only the first receiver location on the
    !!.. propagation grid need to be calculated
    else
       !hgrid%index_datgre_gxmin = 1.5 + (dgrid%agx%o - ggrid%ax_surf%o) / ggrid%ax_surf%d
       !hgrid%index_datgre_gymin = 1.5 + (dgrid%agy%o - ggrid%ay_surf%o) / ggrid%ay_surf%d
       itr = 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
             itr = itr + 1
             hgrid%index_datgre_rec(1, itr) = index_gx
             hgrid%index_datgre_rec(2, itr) = index_gy
          end do
       end do
    end if
!----------------------------------------
!! should include the propagation grid !!
!-----------------------------------------
!    if (hgrid%ax_surf%n > 1) then
!       if (x_modl < pgrid%pxmin .or. x_modl > pgrid%pxmax) hgrid%skip = .true.
!    end if
!    if (hgrid%ay_surf%n > 1) then
!       if (y_modl < pgrid%pymin .or. y_modl > pgrid%pymax) hgrid%skip = .true.
!    end if
    
  end subroutine
!=================================================================================================
  subroutine wei_hessian_array_alloc()
    integer :: memory
    allocate(hessian(hgrid%ahx_offd%n, hgrid%ahy_offd%n, hgrid%ax_offd%n, hgrid%ay_offd%n, hgrid%az_offd%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))
    allocate(hgrid%index_datgre_rec(2, dgrid%ntr))

    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 wei_hessian_array_dealloc()
    deallocate(hessian)
    deallocate(greens)
    deallocate(wavelet)
  end subroutine
!==========================================================================================
  subroutine wei_hessian_param_init()
    real                :: xmin_modl_default, xmax_modl_default
    real                :: ymin_modl_default, ymax_modl_default
    real                :: zmin_modl_default, zmax_modl_default
    integer             :: num_match
    character(len=1024) :: string

    dgrid%hdr_sou=.false. !set source header to zero...no source header
    dgrid%hdr_rec=.false. !set receiver header to zero...no receiver header 

    num_match = getch(hdr_rec_, "s", string) + getch(hdr_sou_, "s", string)
    if (num_match == 2) dgrid%hdr_rec = .true. !both source & receiver header given

    num_match = getch(hdr_sou_, "s", string)
    if (num_match == 1) dgrid%hdr_sou = .true. !source header given?

    num_match = getch(wavelet_, "s", string)   !source wavelet given?
    if (num_match == 1) sou_sig = .true.

    if (dgrid%hdr_sou) then
       call auxin(hdr_sou_)                       !read source ppties
       call from_aux(hdr_sou_, "n2", dgrid%nsh)   !num. of sources   
    else
       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..
       write(0,*)"=================dgrid%nsh: ", dgrid%nsh
    end if

    if (dgrid%hdr_rec) then
       call auxin(hdr_rec_)                     !same as for sources
       call from_aux(hdr_rec_, "n2", dgrid%ntr) !read num. of traces/receivers
    else
       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
    end if
    
    !============================================================
    !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.)
    !============================================================
    !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)
    !============================================================
    ! Freq. params of the greens and data grid must be same.
    if (dgrid%aw%n /= ggrid%aw%n .or. dgrid%aw%d /= ggrid%aw%d .or. dgrid%aw%o /= ggrid%aw%o) &
       call seperr("dgrid%aw .ne. ggrid%aw, error!!!")

    !============================================================
    !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
    !============================================================
    !check whether to calc. sub. offset
    call from_param("offset", offset, .false.)
    if (offset .and. dgrid%hdr_rec) allocate(dgrid_header_rec_save(2, dgrid%ntr)) !receiver header provided
    if (offset .and. .not. dgrid%hdr_rec) then 
       dgrid_gxmin_save = dgrid%agx%o    !define data grid for geophones
       dgrid_gymin_save = dgrid%agy%o
    end if
  end subroutine
!========================================================================================== 
  subroutine wei_hessian_input_history()
    !!!!!Reads the Greens functions properties from the history file...

    call auxin(greens_)

    call from_aux(greens_, "n1", ggrid%ax_surf%n) !surface location num. 
    call from_aux(greens_, "d1", ggrid%ax_surf%d) !surface location samp.
    call from_aux(greens_, "o1", ggrid%ax_surf%o) !surface location orig.

    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) !model location num.
    call from_aux(greens_, "d3", ggrid%ax_modl%d) !model location samp.
    call from_aux(greens_, "o3", ggrid%ax_modl%o) !model location orig.

    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) !greens freq. num.
    call from_aux(greens_, "d6",      ggrid%aw%d) !greens freq. samp.
    call from_aux(greens_, "o6",      ggrid%aw%o) !greens freq. orig.

    if (sou_sig) then
       call from_aux(wavelet_, "n1", dgrid%aw%n)  !source wave. freq. num. 
       call from_aux(wavelet_, "d1", dgrid%aw%d)  !source wave. freq. samp.
       call from_aux(wavelet_, "o1", dgrid%aw%o)  !source wave. freq. orig.
    else
       dgrid%aw = ggrid%aw                        !else use same as green's grid  
    end if
  end subroutine
!==========================================================================================
  subroutine wei_hessian_output_history()
    !!!!!Outputs the Hessain details into the history file...
    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 wei_hessian_io_init()
    real, allocatable, dimension(:) :: ctmp
    integer                         :: i, nloop
    allocate(ctmp(hgrid%ahx_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*hgrid%ahy_modl%n
    do i=1, nloop
       call sep_write(ctmp, hessian_, esize=4)
    end do
    deallocate(ctmp)
  end subroutine
!==============================================================================
!initialize the source/receiver header
subroutine wei_hessian_io_head_init(flag)
    character(len=4), intent(in) :: flag
    if(flag == "hrec") then
       allocate(dgrid%header_rec(2, dgrid%ntr))
       dgrid%header_rec = 0.
    else if(flag == "hsou") then
       allocate(dgrid%header_sou(2))
       dgrid%header_sou = 0.
    else
       call seperr("wei_io_head_init: you have to specify the type of header files")
    end if
  end subroutine
!==============================================================================
!read the reciever/source header
  subroutine wei_hessian_io_head_oper(flag, ish)
    character(len=4), intent(in) :: flag
    integer                      :: ish
    integer, dimension(2)        :: n_cube_sou, j_wind_sou, n_wind_sou, f_wind_sou
    integer, dimension(3)        :: n_cube_rec, j_wind_rec, n_wind_rec, f_wind_rec
    integer                      :: n_dim_cube_sou, n_dim_cube_rec, err
    if(flag == "hrec") then
       n_dim_cube_rec = 3
       n_cube_rec = (/2, dgrid%ntr, dgrid%nsh/)
       j_wind_rec = (/1, 1, 1/)
       n_wind_rec = (/2, dgrid%ntr, 1/)
       f_wind_rec = (/0, 0, ish-1/)
       err = sreed_window(hdr_rec_,n_dim_cube_rec,n_cube_rec,n_wind_rec,f_wind_rec,j_wind_rec,4,dgrid%header_rec)

    else if (flag == "hsou") then
       n_dim_cube_sou = 2
       n_cube_sou = (/2, dgrid%nsh/)
       j_wind_sou = (/1, 1/)
       n_wind_sou = (/2, 1/)
       f_wind_sou = (/0, ish-1/)
       err = sreed_window(hdr_sou_,n_dim_cube_sou,n_cube_sou,n_wind_sou,f_wind_sou,j_wind_sou,4,dgrid%header_sou)
    else
       call seperr("wei_io_head_lop: you have to specify the type of header files")

    end if
  end subroutine

!====================================================================================================================
!write out the offdiagonal terms for a given image point.
subroutine wei_hessian_io_hessian(flag, hessian, index_hx_modl, index_hy_modl, index_x_modl, index_y_modl, index_z_modl)
    character(len=4)                    :: flag
    integer                             :: index_hx_modl, index_hy_modl
    integer                             :: index_x_modl, index_y_modl, index_z_modl
    real, dimension(:,:,:,:,:), pointer :: hessian
    integer, dimension(10)              :: n_cube, j_wind, n_wind, f_wind
    integer                             :: n_dim_cube, err

    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, 1, 1, 1, 1, 1/)
    f_wind = (/0, 0, 0, 0, 0, index_hx_modl-1, index_hy_modl-1, index_x_modl-1, index_y_modl-1, index_z_modl-1/)
    if (flag == "rite") then
       err = srite_window(hessian_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessian)
    else if (flag == "reed") then
       hessian = 0.
       err = sreed_window(hessian_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessian)
    else
       call seperr("you have to specify the flag for wei_hessian_io_hessian()")

    end if
  end subroutine
!====================================================================================================================
! read a block of Green's function, nxs_modl = max(hgrid%ahx_modl%n, hgrid%ahx_offd%n)
  subroutine wei_hessian_io_greens(greens, nxs_modl, nys_modl, nzs_modl, nws, &
                                   index_xs_modl, index_ys_modl, index_zs_modl, index_ws)
    integer                                   :: nxs_modl, nys_modl, nzs_modl, nws
    integer                                   :: index_xs_modl, index_ys_modl, index_zs_modl, index_ws
    integer, dimension(6)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    integer                                   :: indx_xmin, indx_xmax
    integer                                   :: indx_ymin, indx_ymax
    complex, dimension(:,:,:,:,:), pointer    :: 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, nxs_modl, nys_modl, nzs_modl, nws/)
    f_wind = (/0, 0, index_xs_modl-1, index_ys_modl-1, index_zs_modl-1, index_ws-1/)

    greens = 0.
    err = sreed_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,greens)
    if (err /= 0) call seperr("wei_hessian_io_greens: error reading Green's function")

  end subroutine
!-------------------------------------------------------------------------------------------------------------
  subroutine wei_hessian_io_greens_test(greens, index_ws)
    integer                                   :: nxs_modl, nys_modl, nzs_modl, nws
    integer                                   :: index_xs_modl, index_ys_modl, index_zs_modl, index_ws
    integer, dimension(6)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    integer                                   :: indx_xmin, indx_xmax
    integer                                   :: indx_ymin, indx_ymax
    complex, dimension(:,:,:,:,:), pointer    :: 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.
    err = sreed_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,greens)
    if (err /= 0) call seperr("wei_hessian_io_greens: error reading Green's function")

    if (.not. upward) greens = conjg(greens)
  end subroutine

!========================================================================================================================
! For a given model point, read Green's function for one receiver location
  subroutine wei_hessian_io_greens_otr(greens_opt, index_xs_surf, index_ys_surf, index_zs_surf, &
                                       index_xs_modl, index_ys_modl, index_zs_modl, index_ws)
    complex                                   :: greens_opt
    
    integer, dimension(6)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    integer                                   :: index_xs_surf, index_ys_surf, index_zs_surf
    integer                                   :: index_xs_modl, index_ys_modl, index_zs_modl, index_ws
    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 = (/1, 1, 1, 1, 1, 1/)
    f_wind = (/index_xs_surf-1, index_ys_surf-1, index_xs_modl-1, index_ys_modl-1, index_zs_modl-1, index_ws-1/)

    greens = 0.
    err = sreed_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,greens_opt)
  end subroutine
!======================================================================================
!read the source wavelet
subroutine wei_hessian_io_wavelet(wavelet, indw_in, nws_in)
    complex, dimension(:), pointer :: wavelet
    integer                        :: indw_in, nws_in
    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 = nws_in
    f_wind = indw_in - 1
    wavelet = 0.
!    write(0,*)"n_dim_cube, n_cube, n_wind, j_wind,f_wind",n_dim_cube, n_cube, n_wind, j_wind, f_wind, wavelet_ ; 
    err = sreed_window(wavelet_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,wavelet)
!   write(0,*)'hmmmmmmmmmmm2';
    if (err.lt.0) call seperr("wei_io_csou: trouble reading source wavefield")
!    write(0,*)'hmmmmmmmmmmm3'; 
  end subroutine


end module




