module wemva_offdso_incore_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

contains


  integer function wemva_offdso_incore_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

    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_incore_clos() result(stat)
    stat = 0
  end function

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

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

    call wemva_offdso_incore(adj, modl, data)

    stat = 0
  end function 

  subroutine wemva_offdso_incore(adj, modl, data)
    logical                               :: adj, add
    real, dimension(nx,ny,nz,nhx,nhy,nhz) :: modl
    real, dimension(nx,ny,nz,nhx,nhy,nhz) :: data

    integer                               :: ihx, ihy, ihz
    real                                  :: hx, hy, hz, hh

    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
          hh = sqrt(hx*hx+hy*hy+hz*hz)
          modl(:,:,:,ihx,ihy,ihz) = modl(:,:,:,ihx,ihy,ihz) + hh*data(:,:,:,ihx,ihy,ihz)/hmax
       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
          hh = sqrt(hx*hx+hy*hy+hz*hz)
          data(:,:,:,ihx,ihy,ihz) = data(:,:,:,ihx,ihy,ihz) + hh*modl(:,:,:,ihx,ihy,ihz)/hmax
       end do
       end do
       end do
    end if  
  end subroutine
end module
