box_ahelmh Subroutine

  subroutine box_ahelmh(mg, id, nc, i_out)
    type(mg_t), intent(inout) :: mg
    integer, intent(in)       :: id
    integer, intent(in)       :: nc
    integer, intent(in)       :: i_out
    !< Index of variable to store Helmholtz in
    integer                   :: i, j, k
    real(dp)                  :: idr2(2*3), a0(2*3), u0, u(2*3), a(2*3)

    ! Duplicate 1/dr^2 array to multiply neighbor values
    idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
    idr2(2:2*3:2) = idr2(1:2*3:2)

    associate (cc => mg%boxes(id)%cc, n => mg_iphi, i_eps1 => mg_iveps1,&
        i_eps2 => mg_iveps2, i_eps3 => mg_iveps3)
      do k = 1, nc
         do j = 1, nc
            do i = 1, nc
               u0 = cc(i, j, k, n)
               a0(1:2) = cc(i, j, k, i_eps1)
               a0(3:4) = cc(i, j, k, i_eps2)
               a0(5:6) = cc(i, j, k, i_eps3)
               u(1:2) = cc(i-1:i+1:2, j, k, n)
               u(3:4) = cc(i, j-1:j+1:2, k, n)
               u(5:6) = cc(i, j, k-1:k+1:2, n)
               a(1:2) = cc(i-1:i+1:2, j, k, i_eps1)
               a(3:4) = cc(i, j-1:j+1:2, k, i_eps2)
               a(5:6) = cc(i, j, k-1:k+1:2, i_eps3)

               cc(i, j, k, i_out) = sum(2 * idr2 * a0(:)*a(:)/(a0(:) + a(:)) * &
                  (u(:) - u0)) - ahelmholtz_lambda * u0
            end do
         end do
      end do
    end associate
  end subroutine box_ahelmh