program main
  use base_codtrs_mod
  use base_struct_mod
  use base_ioutil_mod
  use sep
  implicit none


  type(axis)                         :: ax_modl, ay_modl, az_modl
  type(axis)                         :: ax_offd, ay_offd, az_offd
  real   , dimension(:,:,:), pointer :: hessin
  real   , dimension(:)    , pointer :: hessou
  integer                            :: cart(3), ndim(3)
  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
  integer                            :: index_x_offd, index_y_offd, index_z_offd
  integer                            :: nmodl, index_offd
  real                               :: x_offd, y_offd, z_offd

  integer, dimension(6)              :: n_cube, j_wind, n_wind, f_wind
  integer                            :: n_dim_cube, err
  integer, external                  :: sreed_window, srite_window
  character(len=256)                 :: hessin_, hessou_

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()


  hessin_ = "hessin"
  hessou_ = "hessou"

  call auxin(hessin_)
  call auxout(hessou_)

  call hinaxis(hessin_, ax_offd, 1)
  call hinaxis(hessin_, ay_offd, 2)
  call hinaxis(hessin_, az_offd, 3)
  call hinaxis(hessin_, ax_modl, 4)
  call hinaxis(hessin_, ay_modl, 5)
  call hinaxis(hessin_, az_modl, 6)

  nmodl = ax_modl%n*ay_modl%n*az_modl%n

  call auxputch("n1", "i", nmodl, hessou_)
  call auxputch("n2", "i", nmodl, hessou_)

  n_dim_cube = 6
  n_cube = (/ax_offd%n,ay_offd%n,az_offd%n,ax_modl%n,ay_modl%n,az_modl%n/)
  j_wind = (/1, 1, 1, 1, 1, 1/)
  n_wind = (/ax_offd%n,ay_offd%n,az_offd%n, 1, 1, 1/)


  allocate(hessin(ax_offd%n,ay_offd%n,az_offd%n))
  allocate(hessou(nmodl))

  ndim(1) = ax_modl%n; ndim(2) = ay_modl%n; ndim(3) = az_modl%n

  do iz_modl=1, az_modl%n
  do iy_modl=1, ay_modl%n
  do ix_modl=1, ax_modl%n

     f_wind = (/0, 0, 0, ix_modl-1, iy_modl-1, iz_modl-1/)

     err = sreed_window(hessin_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hessin)

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

do iz_offd=1, 1
do iy_offd=1, 1
do ix_offd=1, 1
  jx_offd=0; jy_offd=0; jz_offd=0
        index_x_offd = max(1, min(ix_modl+jx_offd,ax_modl%n))
        index_y_offd = max(1, min(iy_modl+jy_offd,ay_modl%n))
        index_z_offd = max(1, min(iz_modl+jz_offd,az_modl%n))

        cart(1) = index_x_offd; cart(2) = index_y_offd; cart(3) = index_z_offd 

        call cart2line(ndim, cart, index_offd)

        hessou(index_offd) = hessin(index_x_offd,index_y_offd,index_z_offd)

     end do
     end do
     end do 

     call sep_write(hessou, hessou_)

  end do
  end do
  end do


end program
