module base_grad3d_mod
  use  base_adjnul_mod
  implicit none


  integer, private :: n1, n2, n3

contains


  integer function base_grad3d_init(n1_in, n2_in, n3_in) result(stat)
    integer :: n1_in, n2_in, n3_in
    n1 = n1_in
    n2 = n2_in
    n3 = n3_in
    stat = 0
  end function

  integer function base_grad3d_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    real, dimension(:) :: modl, data
    call base_adjnul_real(adj, add, modl, data) 
    call base_grad3d(adj, modl, data)
    stat = 0
  end function

  integer function base_grad3d_clos() result(stat)
    n1 = 0
    n2 = 0
    n3 = 0
    stat = 0
  end function


  subroutine base_grad3d(adj, modl, data)
    logical                     :: adj
    real, dimension(n1,n2,n3)   :: modl
    real, dimension(n1,n2,n3,3) :: data

    integer                     :: i, j, k

    if (n2 > 1) then

       do k=1, n3-1
       do j=1, n2-1
       do i=1, n1-1
          if (adj) then
             modl(i+1,j,k) = modl(i+1,j,k) + data(i,j,k,1)
             modl(i  ,j,k) = modl(i  ,j,k) - data(i,j,k,1)
             modl(i,j+1,k) = modl(i,j+1,k) + data(i,j,k,2)
             modl(i,j  ,k) = modl(i,j  ,k) - data(i,j,k,2)
             modl(i,j,k+1) = modl(i,j,k+1) + data(i,j,k,3)
             modl(i,j,k  ) = modl(i,j,k  ) - data(i,j,k,3)
          else
             data(i,j,k,1) = data(i,j,k,1) + modl(i+1,j,k) - modl(i,j,k)
             data(i,j,k,2) = data(i,j,k,2) + modl(i,j+1,k) - modl(i,j,k) 
             data(i,j,k,3) = data(i,j,k,3) + modl(i,j,k+1) - modl(i,j,k)
          end if
       end do
       end do
       end do
    else
       do k=1, n3-1
       do i=1, n1-1
          if (adj) then
             modl(i+1,1,k) = modl(i+1,1,k) + data(i,1,k,1)
             modl(i  ,1,k) = modl(i  ,1,k) - data(i,1,k,1)
             modl(i,1,k+1) = modl(i,1,k+1) + data(i,1,k,3)
             modl(i,1,k  ) = modl(i,1,k  ) - data(i,1,k,3)
          else
             data(i,1,k,1) = data(i,1,k,1) + modl(i+1,1,k) - modl(i,1,k)
             data(i,1,k,3) = data(i,1,k,3) + modl(i,1,k+1) - modl(i,1,k) 
          end if 
       end do
       end do
    end if
  end subroutine


end module








