module wemva_offdso_mod
  use sep
  implicit none

  integer                              , private :: nx, ny, nz
  integer                              , private :: nhx, nhy
  real                                 , private :: ohx, ohy
  real                                 , private :: dhx, dhy 
  real                                 , private :: hmax, hcut, hord
  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,                           &
                                     ohx_in, ohy_in, dhx_in, dhy_in, hcut_in, hord_in) result(stat)
    integer :: nx_in, ny_in, nz_in
    integer :: nhx_in, nhy_in
    real    :: ohx_in, ohy_in
    real    :: dhx_in, dhy_in
    real    :: hcut_in, hord_in
    real    :: hxmax, hymax, hzmax

    nx  = nx_in ; ny  = ny_in ; nz  = nz_in
    nhx = nhx_in; nhy = nhy_in
    ohx = ohx_in; ohy = ohy_in
    dhx = dhx_in; dhy = dhy_in

    hcut = hcut_in
    hord = hord_in

    allocate(modl(nhx,nhy,nx,ny))
    allocate(data(nhx,nhy,nx,ny))

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

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

hmax = 1.

write(0,*)"nx,ny,nz", nx,ny,nz
write(0,*)"nhx,nhy", nhx, nhy
write(0,*)"ohx,ohy,dhx,dhy", ohx,ohy,dhx,dhy


    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            :: iz, ihx, ihy
    real               :: hx, hy
    real               :: hh

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

    if (adj) then
       do iz=1, nz
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, iz, modl_, "reed")
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, iz, data_, "reed")
          do ihy=1, nhy; hy = ohy + (ihy-1)*dhy
          do ihx=1, nhx; hx = ohx + (ihx-1)*dhx
             hh = max(sqrt(hx*hx+hy*hy)-hcut, 0.)**hord
             modl(ihx,ihy,:,:) = modl(ihx,ihy,:,:) + hh*data(ihx,ihy,:,:)/hmax
          end do  
          end do
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, iz, modl_, "rite")
       end do

    else
       do iz=1, nz
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, iz, data_, "reed")
          call wemva_offdso_flio(modl, nx, ny, nz, nhx, nhy, iz, modl_, "reed")
          do ihy=1, nhy; hy = ohy + (ihy-1)*dhy
          do ihx=1, nhx; hx = ohx + (ihx-1)*dhx
             hh = max(sqrt(hx*hx+hy*hy)-hcut, 0.)**hord
             data(ihx,ihy,:,:) = data(ihx,ihy,:,:) + hh*modl(ihx,ihy,:,:)/hmax
          end do
          end do
          call wemva_offdso_flio(data, nx, ny, nz, nhx, nhy, iz, data_, "rite")
       end do 

    end if  
    stat = 0
  end function


  subroutine wemva_offdso_flio(cube, nx, ny, nz, nhx, nhy, indz, filetag, flag)
    real, dimension(nhx,nhy,nx,ny) :: cube
    integer                        :: nx, ny, nz, nhx, nhy, indz
    character(len=*)               :: filetag
    character(len=4)               :: flag

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

    n_dim_cube = 5
    n_cube = (/nhx,nhy,nx,ny,nz/)
    j_wind = (/1,1,1,1,1/)
    n_wind = (/nhx,nhy,nx,ny,1/)
    f_wind = (/0,0,0,0,indz-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
