bounded_to_local_neighbors Function

   real(rk) function bounded_to_local_neighbors(mesh, phi, candidate, cell, boundary_value, include_boundary) result(out)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(in) :: candidate
      integer, intent(in) :: cell
      real(rk), intent(in) :: boundary_value
      logical, intent(in) :: include_boundary

      integer :: lf, fid, nb
      real(rk) :: lo, hi

      lo = phi(cell)
      hi = phi(cell)
      do lf = 1, mesh%ncell_faces(cell)
         fid = mesh%cell_faces(lf, cell)
         nb = adjacent_cell(mesh, fid, cell)
         if (nb > 0) then
            lo = min(lo, phi(nb))
            hi = max(hi, phi(nb))
         end if
      end do
      if (include_boundary) then
         lo = min(lo, boundary_value)
         hi = max(hi, boundary_value)
      end if

      out = clamp_value(candidate, lo, hi)
   end function bounded_to_local_neighbors