!!.. Space domain non-stationary deconvolution
module hes_sdecon_mod
  use  base_ioutil_mod
  use  base_struct_mod
  use  base_fftwv3_mod
  use  hes_locprd_mod
  use  csolver_tiny_mod
  use  ccgstep_mod
!  use  solver_tiny_mod
!  use  cgstep_mod
  implicit none


  type(hes_grid)                        , private :: hgrid
  character(len=256)                    , private :: filetag

  real   , dimension(:,:,:), allocatable, private :: ldat, lmod, lhes
  real   , dimension(:,:,:), allocatable, private :: hess

  integer, external                     , private :: sreed_window, srite_window
 

contains

  integer function hes_sdecon_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(ldat(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))
    allocate(lmod(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))
    allocate(lhes(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))
    allocate(hess(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))

    ldat = 0.; lmod = 0.; lhes = 0.

    stat = 0
  end function

  integer function hes_sdecon_clos() result(stat)
    deallocate(ldat)
    deallocate(lmod)
    deallocate(lhes)
    deallocate(hess)
    stat = 0
  end function

  integer function hes_sdecon_oper(modl, data) result(stat)
    real, dimension(:,:,:) :: modl  !!.. the totoal image size, output
    real, dimension(:,:,:) :: data  !!.. the totoal image size, input

    stat = 1

    call sdecon_kernel(modl, data)

    stat = 0
  end function




  subroutine sdecon_kernel(modl, data)
    real, dimension(:,:,:) :: modl  !!.. the totoal image size, output
    real, dimension(:,:,:) :: data  !!.. the totoal image size, input

    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                :: stat


    modl = 0.

    do iz_modl=1, hgrid%az_modl%n; write(0,*)"iz_modl=",iz_modl 
    do iy_modl=1, hgrid%ay_modl%n
    do ix_modl=1, hgrid%ax_modl%n

       ldat = 0. 


       call deconv_fileio(hess, hgrid, ix_modl, iy_modl, iz_modl, filetag, "reed")


       !!.. extract the locally windowd images
       do iz_offd=1, hgrid%az_offd%n; index_z_offd = iz_modl+iz_offd-hgrid%az_offd%n/2-1; if(index_z_offd<1.or.index_z_offd>hgrid%az_modl%n) cycle
       do iy_offd=1, hgrid%ay_offd%n; index_y_offd = iy_modl+iy_offd-hgrid%ay_offd%n/2-1; if(index_y_offd<1.or.index_y_offd>hgrid%ay_modl%n) cycle
       do ix_offd=1, hgrid%ax_offd%n; index_x_offd = ix_modl+ix_offd-hgrid%ax_offd%n/2-1; if(index_x_offd<1.or.index_x_offd>hgrid%ax_modl%n) cycle

          modl(ix_modl,iy_modl,iz_modl) = modl(ix_modl,iy_modl,iz_modl) + hess(ix_offd,iy_offd,iz_offd)*data(index_x_ofd,index_y_offd,index_z_offd)
          ldat(ix_offd,iy_offd,iz_offd) = data(index_x_offd,index_y_offd,index_z_offd)

       end do
       end do
       end do


       !!.. read the off-diagonals of the Hessian
       call deconv_fileio(hess, hgrid, ix_modl, iy_modl, iz_modl, filetag, "reed")
       lhes = hess






       
       index_x_modl = hgrid%ax_offd%n/2 + 1
       index_y_modl = hgrid%ay_offd%n/2 + 1
       index_z_modl = hgrid%az_offd%n/2 + 1

       !!.. assign back the inversion result
       modl(ix_modl,iy_modl,iz_modl) = real(lmod(index_x_modl,index_y_modl,index_z_modl))

!       do iz_offd=1, hgrid%az_offd%n; index_z_offd = iz_modl+iz_offd-hgrid%az_offd%n/2-1; if(index_z_offd<1.or.index_z_offd>hgrid%az_modl%n) cycle
!       do iy_offd=1, hgrid%ay_offd%n; index_y_offd = iy_modl+iy_offd-hgrid%ay_offd%n/2-1; if(index_y_offd<1.or.index_y_offd>hgrid%ay_modl%n) cycle
!       do ix_offd=1, hgrid%ax_offd%n; index_x_offd = ix_modl+ix_offd-hgrid%ax_offd%n/2-1; if(index_x_offd<1.or.index_x_offd>hgrid%ax_modl%n) cycle

!          modl(index_x_offd,index_y_offd,index_z_offd) = modl(index_x_offd,index_y_offd,index_z_offd)+real(lmod(ix_offd,iy_offd,iz_offd))

!       end do
!       end do
!       end do


    end do
    end do
    end do
    return
  end subroutine


  subroutine deconv_fileio(hess, hgrid, index_x_modl, index_y_modl, index_z_modl, filetag, flag)
    real, dimension(:,:,:) :: hess
    type(hes_grid)         :: hgrid
    integer                :: index_x_modl, index_y_modl, index_z_modl
    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_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,1,1,1/)
    f_wind = (/0,0,0,index_x_modl-1,index_y_modl-1,index_z_modl-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


  subroutine deconv_fileio_test(hess, hgrid, index_x_modl, index_y_modl, index_z_modl, filetag, flag)
    complex, dimension(:,:,:) :: hess
    type(hes_grid)            :: hgrid
    integer                   :: index_x_modl, index_y_modl, index_z_modl
    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_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,1,1,1/)
    f_wind = (/0,0,0,index_x_modl-1,index_y_modl-1,index_z_modl-1/)

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

    return
  end subroutine




end module
