module wemva_offdso_mod
  use sep
  implicit none

  integer                            , private :: nx, ny, nz
  integer                            , private :: nhx, nhy, nhz
  real                               , private :: ohx, ohy, ohz
  real                               , private :: dhx, dhy, dhz 
  real                               , private :: hmax
  real, dimension(:,:,:), allocatable, private :: modl
  real, dimension(:,:,:), allocatable, private :: data

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

contains


  integer function wemva_offdso_init(nx_in, ny_in, nz_in, nhx_in, nhy_in, nhz_in, &
                                     ohx_in, ohy_in, ohz_in, dhx_in, dhy_in, dhz_in) result(stat)
    integer :: nx_in, ny_in, nz_in
    integer :: nhx_in, nhy_in, nhz_in
    real    :: ohx_in, ohy_in, ohz_in
    real    :: dhx_in, dhy_in, dhz_in
    real    :: hxmax, hymax, hzmax

    nx  = nx_in ; ny  = ny_in ; nz  = nz_in
    nhx = nhx_in; nhy = nhy_in; nhz = nhz_in
    ohx = ohx_in; ohy = ohy_in; ohz = ohz_in
    dhx = dhx_in; dhy = dhy_in; dhz = dhz_in
    allocate(modl(nx, ny, nz))
    allocate(data(nx, ny, nz))

    hxmax = max(abs(ohx+(nhx-1)*dhx), abs(ohx))
    hymax = max(abs(ohy+(nhy-1)*dhy), abs(ohy))
    hzmax = max(abs(ohz+(nhz-1)*dhz), abs(ohz))

    hmax = sqrt(hxmax*hxmax+hymax*hymax+hzmax*hzmax)

write(0,*)"hmax = ", hmax

    stat = 0
  end function

  integer function wemva_offdso_clos() result(stat)
    deallocate(modl)
    deallocate(data)
    stat = 0
  end function

  integer function wemva_offdso_oper(adj, add, modl_, data_) result(stat)
    logical            :: adj, add
    character(len=*)   :: modl_, data_
    integer            :: ihx, ihy, ihz
    real               :: hx, hy, hz

    if (.not.add) then
       if (adj) then
          modl = 0.
          do ihz=1, nhz
          do ihy=1, nhy
          do ihx=1, nhx
             call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, modl_, "rite")
          end do
          end do
          end do
       else
          data = 0.
          do ihz=1, nhz
          do ihy=1, nhy
          do ihx=1, nhx
             call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, data_, "rite")
          end do
          end do
          end do 
       end if
    end if

    if (adj) then
       do ihz=1, nhz; hz = ohz + (ihz-1)*dhz
       do ihy=1, nhy; hy = ohy + (ihy-1)*dhy
       do ihx=1, nhx; hx = ohx + (ihx-1)*dhx
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, data_, "reed") 
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, modl_, "reed")
          modl = modl + sqrt(hx*hx+hy*hy+hz*hz)*data/hmax
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, modl_, "rite")
       end do
       end do
       end do 
    else
       do ihz=1, nhz; hz = ohz + (ihz-1)*dhz
       do ihy=1, nhy; hy = ohy + (ihy-1)*dhy
       do ihx=1, nhx; hx = ohx + (ihx-1)*dhx
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, modl_, "reed")
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, data_, "reed")
          data = data + sqrt(hx*hx+hy*hy+hz*hz)*modl/hmax
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, nhz, ihx, ihy, ihz, data_, "rite")
       end do
       end do
       end do
    end if  
    stat = 0
  end function


  subroutine wemva_offdso_flio(cube, nx, ny, nz, nhx, nhy, nhz, index_hx, index_hy, index_hz, filetag, flag)
    real, dimension(nx,ny,nz) :: cube
    integer                   :: nx, ny, nz, nhx, nhy, nhz
    integer                   :: index_hx, index_hy, index_hz
    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 = (/nx, ny, nz, nhx, nhy, nhz/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/nx,ny,nz,1,1,1/)
    f_wind = (/0,0,0,index_hx-1,index_hy-1,index_hz-1/)

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

    return
  end subroutine
end module
