module hes_convol_mod2
  use  base_struct_mod
  use  base_ioutil_mod
  use  base_codtrs_mod

  implicit none

  type(hes_grid)                        , private :: hgrid
  character(len=256)                    , private :: filetag
  real   , dimension(:,:,:), allocatable, private :: hess
  integer, external                     , private :: sreed_window, srite_window

contains

  integer function hes_convol_init(hgrid_in, filetag_in) result(stat)
    type(hes_grid)        :: hgrid_in
    character(len=*)      :: filetag_in

    stat = 1

    hgrid = hgrid_in

    filetag = filetag_in

    allocate(hess(hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n))
    hess = 0.

    stat = 0

  end function

  integer function hes_convol_clos() result(stat)
    deallocate(hess)
    stat = 0
  end function

  integer function hes_convol_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    real, dimension(:) :: modl, data

    stat = 1

    if (.not. add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if  
    end if  

    call convol_kernel(adj, modl, data)
    stat = 0
  end function

  subroutine convol_kernel(adj, modl, data)
    logical                :: adj
    real, dimension(hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n) :: modl
    real, dimension(hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n) :: data

    integer                :: ix_modl, iy_modl, iz_modl
    integer                :: ix_offd, iy_offd, iz_offd
    integer                :: jx_offd, jy_offd, jz_offd
    integer                :: index_x_modl, index_y_modl, index_z_modl
    real                   :: x_offd, y_offd, z_offd


    do iz_offd=1, hgrid%az_offd%n; z_offd = hgrid%az_offd%o+(iz_offd-1)*hgrid%az_offd%d; jz_offd = nint(z_offd/hgrid%az_modl%d)
    do iy_offd=1, hgrid%ay_offd%n; y_offd = hgrid%ay_offd%o+(iy_offd-1)*hgrid%ay_offd%d; jy_offd = nint(y_offd/hgrid%ay_modl%d)
    do ix_offd=1, hgrid%ax_offd%n; x_offd = hgrid%ax_offd%o+(ix_offd-1)*hgrid%ax_offd%d; jx_offd = nint(x_offd/hgrid%ax_modl%d)

       call convol_fileio(hess, hgrid, ix_offd, iy_offd, iz_offd, filetag, "reed")

       do iz_modl=1, hgrid%az_modl%n; index_z_modl = max(1, min(iz_modl-jz_offd,hgrid%az_modl%n))
       do iy_modl=1, hgrid%ay_modl%n; index_y_modl = max(1, min(iy_modl-jy_offd,hgrid%ay_modl%n))
       do ix_modl=1, hgrid%ax_modl%n; index_x_modl = max(1, min(ix_modl-jx_offd,hgrid%ax_modl%n))
          if (adj) then
             modl(index_x_modl,index_y_modl,index_z_modl) = modl(index_x_modl,index_y_modl,index_z_modl) &
                                                          + hess(ix_modl,iy_modl,iz_modl)*data(ix_modl    ,iy_modl      ,iz_modl     )
          else 
             data(ix_modl     ,iy_modl     ,iz_modl     ) = data(ix_modl     ,iy_modl     ,iz_modl     ) &
                                                          + hess(ix_modl,iy_modl,iz_modl)*modl(index_x_modl,index_y_modl,index_z_modl)
          end if
       end do
       end do
       end do

    end do
    end do
    end do

    return
  end subroutine

  subroutine convol_fileio(hess, hgrid, index_x_offd, index_y_offd, index_z_offd, filetag, flag)
    real, dimension(:,:,:) :: hess
    type(hes_grid)         :: hgrid
    integer                :: index_x_offd, index_y_offd, index_z_offd
    character(len=*)       :: filetag
    character(len=4)       :: flag

    integer, dimension(6)  :: n_cube, j_wind, n_wind, f_wind
    integer                :: n_dim_cube, err

    n_dim_cube = 6
    n_cube = (/hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n,hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n,1,1,1/)
    f_wind = (/0,0,0,index_x_offd-1,index_y_offd-1,index_z_offd-1/)

    if (flag == "rite") then
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hess)
    else if (flag == "reed") then
       hess = 0.
       err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hess)
    end if

    return
  end subroutine

end module
