diff --git a/CarpetX/param.ccl b/CarpetX/param.ccl
index 2d805b4ed..2b1adc979 100644
--- a/CarpetX/param.ccl
+++ b/CarpetX/param.ccl
@@ -433,7 +433,8 @@ KEYWORD prolongation_type "Prolongation type"
   "interpolate" :: "interpolate between data points"
   "conservative" :: "interpolate cell averages, ensuring conservation"
   "ddf" :: "interpolate in vertex centred and conserve (with one order lower) in cell centred directions"
-  "eno" :: "interpolate in vertex centred and ENO-conserve in cell centred directions"
+  "eno" :: "interpolate in vertex centred and minmod-conserve in cell centred directions"
+  "minmod" :: "interpolate in vertex centred and ENO-conserve in cell centred directions"
   "hermite" :: "Hermite-interpolate in vertex centred and conserve in cell centred directions"
   "natural" :: "interpolate in vertex centred and conserve in cell centred directions, using the same order"
   "poly-cons3lfb" :: "interpolate polynomially in vertex centred directions and conserve with 3rd order accuracy and a linear fallback in cell centred directions"
diff --git a/CarpetX/src/driver.cxx b/CarpetX/src/driver.cxx
index b4f3bc4c4..c0e9a7e4c 100644
--- a/CarpetX/src/driver.cxx
+++ b/CarpetX/src/driver.cxx
@@ -541,6 +541,7 @@ amrex::Interpolater *get_interpolator(const std::array<int, dim> indextype) {
     conservative,
     ddf,
     eno,
+    minmod,
     hermite,
     natural,
     poly_cons3lfb,
@@ -555,6 +556,8 @@ amrex::Interpolater *get_interpolator(const std::array<int, dim> indextype) {
       return interp_t::ddf;
     else if (CCTK_EQUALS(prolongation_type, "eno"))
       return interp_t::eno;
+    else if (CCTK_EQUALS(prolongation_type, "minmod"))
+      return interp_t::minmod;
     else if (CCTK_EQUALS(prolongation_type, "hermite"))
       return interp_t::hermite;
     else if (CCTK_EQUALS(prolongation_type, "natural"))
@@ -894,6 +897,75 @@ amrex::Interpolater *get_interpolator(const std::array<int, dim> indextype) {
     }
     break;
 
+  case interp_t::minmod:
+
+    switch (prolongation_order) {
+
+    case 1:
+      switch ((indextype[0] << 2) | (indextype[1] << 1) | (indextype[2] << 0)) {
+      case 0b000:
+        return &prolongate_minmod_3d_rf2_c000_o1;
+      case 0b001:
+        return &prolongate_minmod_3d_rf2_c001_o1;
+      case 0b010:
+        return &prolongate_minmod_3d_rf2_c010_o1;
+      case 0b011:
+        return &prolongate_minmod_3d_rf2_c011_o1;
+      case 0b100:
+        return &prolongate_minmod_3d_rf2_c100_o1;
+      case 0b101:
+        return &prolongate_minmod_3d_rf2_c101_o1;
+      case 0b110:
+        return &prolongate_minmod_3d_rf2_c110_o1;
+      case 0b111:
+        return &prolongate_minmod_3d_rf2_c111_o1;
+      }
+      break;
+
+    case 3:
+      switch ((indextype[0] << 2) | (indextype[1] << 1) | (indextype[2] << 0)) {
+      case 0b000:
+        return &prolongate_minmod_3d_rf2_c000_o3;
+      case 0b001:
+        return &prolongate_minmod_3d_rf2_c001_o3;
+      case 0b010:
+        return &prolongate_minmod_3d_rf2_c010_o3;
+      case 0b011:
+        return &prolongate_minmod_3d_rf2_c011_o3;
+      case 0b100:
+        return &prolongate_minmod_3d_rf2_c100_o3;
+      case 0b101:
+        return &prolongate_minmod_3d_rf2_c101_o3;
+      case 0b110:
+        return &prolongate_minmod_3d_rf2_c110_o3;
+      case 0b111:
+        return &prolongate_minmod_3d_rf2_c111_o3;
+      }
+      break;
+
+    case 5:
+      switch ((indextype[0] << 2) | (indextype[1] << 1) | (indextype[2] << 0)) {
+      case 0b000:
+        return &prolongate_minmod_3d_rf2_c000_o5;
+      case 0b001:
+        return &prolongate_minmod_3d_rf2_c001_o5;
+      case 0b010:
+        return &prolongate_minmod_3d_rf2_c010_o5;
+      case 0b011:
+        return &prolongate_minmod_3d_rf2_c011_o5;
+      case 0b100:
+        return &prolongate_minmod_3d_rf2_c100_o5;
+      case 0b101:
+        return &prolongate_minmod_3d_rf2_c101_o5;
+      case 0b110:
+        return &prolongate_minmod_3d_rf2_c110_o5;
+      case 0b111:
+        return &prolongate_minmod_3d_rf2_c111_o5;
+      }
+      break;
+    }
+    break;
+
   case interp_t::hermite:
 
     switch (prolongation_order) {
diff --git a/CarpetX/src/make.code.defn b/CarpetX/src/make.code.defn
index 32f01f91e..f5e43189c 100644
--- a/CarpetX/src/make.code.defn
+++ b/CarpetX/src/make.code.defn
@@ -47,6 +47,7 @@ SRCS =							\
 	prolongate_3d_rf2_impl_ddf.cxx			\
 	prolongate_3d_rf2_impl_eno.cxx			\
 	prolongate_3d_rf2_impl_hermite.cxx		\
+	prolongate_3d_rf2_impl_minmod.cxx		\
 	prolongate_3d_rf2_impl_natural.cxx		\
 	prolongate_3d_rf2_impl_poly.cxx			\
 	prolongate_3d_rf2_impl_poly_cons3lfb.cxx	\
diff --git a/CarpetX/src/prolongate_3d_rf2.cxx b/CarpetX/src/prolongate_3d_rf2.cxx
index 68cd24b83..543cae67a 100644
--- a/CarpetX/src/prolongate_3d_rf2.cxx
+++ b/CarpetX/src/prolongate_3d_rf2.cxx
@@ -25,6 +25,8 @@ std::ostream &operator<<(std::ostream &os, const interpolation_t intp) {
     return os << "cons";
   case interpolation_t::eno:
     return os << "eno";
+  case interpolation_t::minmod:
+    return os << "minmod";
   }
   assert(false);
 }
diff --git a/CarpetX/src/prolongate_3d_rf2.hxx b/CarpetX/src/prolongate_3d_rf2.hxx
index 5eced083f..0eae1da04 100644
--- a/CarpetX/src/prolongate_3d_rf2.hxx
+++ b/CarpetX/src/prolongate_3d_rf2.hxx
@@ -15,12 +15,13 @@ std::ostream &operator<<(std::ostream &os, centering_t cent);
 constexpr auto VC = centering_t::vertex;
 constexpr auto CC = centering_t::cell;
 
-enum class interpolation_t { poly, hermite, cons, eno };
+enum class interpolation_t { poly, hermite, cons, eno, minmod };
 std::ostream &operator<<(std::ostream &os, interpolation_t cent);
 constexpr auto POLY = interpolation_t::poly;
 constexpr auto HERMITE = interpolation_t::hermite;
 constexpr auto CONS = interpolation_t::cons;
 constexpr auto ENO = interpolation_t::eno;
+constexpr auto MINMOD = interpolation_t::minmod;
 
 enum class fallback_t { none, linear };
 std::ostream &operator<<(std::ostream &os, fallback_t fb);
@@ -39,17 +40,22 @@ class prolongate_3d_rf2 final : public amrex::Interpolater {
 
   // Conservative must be one of the possible choices
   static_assert(INTPI == POLY || INTPI == HERMITE || INTPI == CONS ||
-                INTPI == ENO);
+                INTPI == ENO || INTPI == MINMOD);
   static_assert(INTPJ == POLY || INTPJ == HERMITE || INTPJ == CONS ||
-                INTPJ == ENO);
+                INTPJ == ENO || INTPJ == MINMOD);
   static_assert(INTPK == POLY || INTPK == HERMITE || INTPK == CONS ||
-                INTPK == ENO);
+                INTPK == ENO || INTPK == MINMOD);
 
   // Order must be nonnegative
   static_assert(ORDERI >= 0);
   static_assert(ORDERJ >= 0);
   static_assert(ORDERK >= 0);
 
+  // Minmod is always linear
+  static_assert(INTPI == MINMOD ? ORDERI == 1 : true);
+  static_assert(INTPJ == MINMOD ? ORDERJ == 1 : true);
+  static_assert(INTPK == MINMOD ? ORDERK == 1 : true);
+
   // Fallback must be one of the possible choices
   static_assert(FB == FB_NONE || FB == FB_LINEAR);
 
@@ -401,6 +407,59 @@ extern prolongate_3d_rf2<CC, CC, VC, ENO, ENO, POLY, 2, 2, 5, FB_NONE>
 extern prolongate_3d_rf2<CC, CC, CC, ENO, ENO, ENO, 2, 2, 2, FB_NONE>
     prolongate_eno_3d_rf2_c111_o5;
 
+// Minmod (tensor product) interpolation
+
+extern prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o1;
+extern prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o1;
+extern prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o1;
+extern prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o1;
+extern prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o1;
+extern prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o1;
+extern prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o1;
+extern prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o1;
+
+extern prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 3, 3, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o3;
+extern prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 3, 3, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o3;
+extern prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 3, 1, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o3;
+extern prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 3, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o3;
+extern prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 3, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o3;
+extern prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 3, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o3;
+extern prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o3;
+extern prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o3;
+
+extern prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 5, 5, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o5;
+extern prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 5, 5, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o5;
+extern prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 5, 1, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o5;
+extern prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 5, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o5;
+extern prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 5, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o5;
+extern prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 5, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o5;
+extern prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o5;
+extern prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o5;
+
 // Hermite interpolation
 
 extern prolongate_3d_rf2<VC, VC, VC, HERMITE, HERMITE, HERMITE, 1, 1, 1,
diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx
index 608b12293..4ffa1633d 100644
--- a/CarpetX/src/prolongate_3d_rf2_impl.hxx
+++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx
@@ -838,6 +838,42 @@ template <int ORDER> struct interp1d<CC, ENO, ORDER> {
   }
 };
 
+// off=0: left sub-cell
+// off=1: right sub-cell
+template <> struct interp1d<CC, MINMOD, 1> {
+  static constexpr int required_ghosts = 1;
+  CCTK_DEVICE
+  CCTK_HOST constexpr
+      __attribute__((__always_inline__, __flatten__)) std::array<int, 2>
+      stencil_radius(const int shift, const int off) const {
+#ifdef CCTK_DEBUG
+    assert(off == 0 || off == 1);
+#endif
+    return {-1, +1};
+  }
+  template <typename F, typename T = std::invoke_result_t<F, int> >
+  CCTK_DEVICE CCTK_HOST inline __attribute__((__always_inline__, __flatten__)) T
+  operator()(const F &crse, const int shift, const int off) const {
+#ifdef CCTK_DEBUG
+    assert(off == 0 || off == 1);
+#endif
+    // Inspired by Athena-K's prolongation operator; see
+    // <https://github.com/IAS-Astrophysics/athenak/blob/main/src/mesh/prolongation.hpp>
+    using std::copysign, std::fabs, std::fmin;
+    const T cval = crse(0);
+    const T cval_minus = crse(-1);
+    const T cval_plus = crse(+1);
+    const T delta_minus = cval - cval_minus;
+    const T delta_plus = cval_plus - cval;
+    // minmod / 4
+    const T delta =
+        (copysign(T(0.125), delta_minus) + copysign(T(0.125), delta_plus)) *
+        fmin(fabs(delta_minus), fabs(delta_plus));
+    const T sign = off == 0 ? -1 : +1;
+    return cval + sign * delta;
+  }
+};
+
 // Test 1d interpolators
 
 template <centering_t CENT, interpolation_t INTP, int ORDER, typename T>
@@ -948,53 +984,53 @@ template <int ORDER, typename T> struct test_interp1d<CC, CONS, ORDER, T> {
       // Function f, a polynomial
       // const auto f{[&](T x) __attribute__((__always_inline__, __flatten__)) {
       // return (order + 1) * pown(x, order); }}; Integral of f (antiderivative)
-      const auto fint{[&](T x) __attribute__((__always_inline__, __flatten__)){
-          return pown(x, order + 1);
-    }
-  };
-  std::array<T, 2> x1;
-  std::array<T, 2> y1;
-  for (int off = 0; off < 2; ++off) {
-    const T rmin = stencil1d.stencil_radius(0, off)[0];
-    const T rmax = stencil1d.stencil_radius(0, off)[1];
-    assert(rmin <= 0 && rmin >= -nghosts);
-    assert(rmax >= 0 && rmax <= +nghosts);
-    for (int i = -(nghosts + 1); i <= +(nghosts + 1); ++i) {
-      if (i < rmin || i > rmax) {
-        xs[i] = 0 / T(0);
-        ys[i] = 0 / T(0);
-      } else {
-        T x = i + int(CC) / T(2);
-        // T y = f(x);
-        const T dx = 1;
-        const T xlo = x - dx / 2;
-        const T xhi = x + dx / 2;
-        const T y = fint(xhi) - fint(xlo); // average of f over cell
-        xs[i] = x;
-        ys[i] = y;
-      }
-    }
+      const auto fint = [&](T x)
+          __attribute__((__always_inline__, __flatten__)) {
+        return pown(x, order + 1);
+      };
+      std::array<T, 2> x1;
+      std::array<T, 2> y1;
+      for (int off = 0; off < 2; ++off) {
+        const T rmin = stencil1d.stencil_radius(0, off)[0];
+        const T rmax = stencil1d.stencil_radius(0, off)[1];
+        assert(rmin <= 0 && rmin >= -nghosts);
+        assert(rmax >= 0 && rmax <= +nghosts);
+        for (int i = -(nghosts + 1); i <= +(nghosts + 1); ++i) {
+          if (i < rmin || i > rmax) {
+            xs[i] = 0 / T(0);
+            ys[i] = 0 / T(0);
+          } else {
+            T x = i + int(CC) / T(2);
+            // T y = f(x);
+            const T dx = 1;
+            const T xlo = x - dx / 2;
+            const T xhi = x + dx / 2;
+            const T y = fint(xhi) - fint(xlo); // average of f over cell
+            xs[i] = x;
+            ys[i] = y;
+          }
+        }
 
-    x1[off] = int(CC) / T(4) + off / T(2);
-    y1[off] = stencil1d(
-        [&](const int i) __attribute__((__always_inline__, __flatten__)) {
-          assert(i >= rmin);
-          assert(i <= rmax);
-          return ys[i];
-        },
-        0, off);
-    assert(isfinite(y1[off]));
-  } // for off
-  assert(y1[0] / 2 + y1[1] / 2 == ys[0]);
-  const T dx = x1[1] - x1[0];
-  const T xlo = x1[0] - dx / 2;
-  const T xhi = x1[1] + dx / 2;
-  const T yint = fint(xhi) - fint(xlo);
-  assert(y1[0] * dx + y1[1] * dx == yint);
-}
-} // namespace CarpetX
-}
-;
+        x1[off] = int(CC) / T(4) + off / T(2);
+        y1[off] = stencil1d(
+            [&](const int i) __attribute__((__always_inline__, __flatten__)) {
+              assert(i >= rmin);
+              assert(i <= rmax);
+              return ys[i];
+            },
+            0, off);
+        assert(isfinite(y1[off]));
+      } // for off
+      // Ensure conservation
+      assert(y1[0] / 2 + y1[1] / 2 == ys[0]);
+      const T dx = x1[1] - x1[0];
+      const T xlo = x1[0] - dx / 2;
+      const T xhi = x1[1] + dx / 2;
+      const T yint = fint(xhi) - fint(xlo);
+      assert(y1[0] * dx + y1[1] * dx == yint);
+    }
+  } // namespace CarpetX
+};
 
 template <int ORDER, typename T> struct test_interp1d<VC, CONS, ORDER, T> {
   test_interp1d() {
@@ -1022,52 +1058,112 @@ template <int ORDER, typename T> struct test_interp1d<CC, ENO, ORDER, T> {
         // const auto f{[&](T x) __attribute__((__always_inline__, __flatten__))
         // { return (order + 1) * pown(x, order); }}; Integral of f
         // (antiderivative)
-        const auto fint{[&](T x) __attribute__((
-            __always_inline__, __flatten__)){return pown(x, order + 1);
+        const auto fint = [&](T x)
+            __attribute__((__always_inline__, __flatten__)) {
+          return pown(x, order + 1);
+        };
+        std::array<T, 2> x1;
+        std::array<T, 2> y1;
+        for (int off = 0; off < 2; ++off) {
+          const T rmin = stencil1d.stencil_radius(0, off)[0];
+          const T rmax = stencil1d.stencil_radius(0, off)[1];
+          assert(rmin >= -nghosts && rmax <= +nghosts);
+          assert(rmin == -(ORDER / 2) && rmax == +(ORDER / 2));
+          for (int i = -1; i < n + 1; ++i) {
+            if (i - i0 < rmin || i - i0 > rmax) {
+              xs[i + 1] = 0 / T(0);
+              ys[i + 1] = 0 / T(0);
+            } else {
+              T x = (i - i0) + int(CC) / T(2);
+              // T y = f(x);
+              const T dx = 1;
+              const T xlo = x - dx / 2;
+              const T xhi = x + dx / 2;
+              const T y = fint(xhi) - fint(xlo); // average of f over cell
+              xs[i + 1] = x;
+              ys[i + 1] = y;
+            }
+          }
+
+          x1[off] = int(CC) / T(4) + off / T(2);
+          y1[off] = stencil1d(
+              [&](int i) __attribute__((__always_inline__, __flatten__)) {
+                return ys[i0 + 1 + i];
+              },
+              0, off);
+          assert(isfinite(y1[off]));
+        } // for off
+        // Ensure conservation
+        assert(y1[0] / 2 + y1[1] / 2 == ys[i0 + 1]);
+        const T dx = x1[1] - x1[0];
+        const T xlo = x1[0] - dx / 2;
+        const T xhi = x1[1] + dx / 2;
+        const T yint = fint(xhi) - fint(xlo);
+        assert(y1[0] * dx + y1[1] * dx == yint);
       }
-    };
-    std::array<T, 2> x1;
-    std::array<T, 2> y1;
-    for (int off = 0; off < 2; ++off) {
-      const T rmin = stencil1d.stencil_radius(0, off)[0];
-      const T rmax = stencil1d.stencil_radius(0, off)[1];
-      assert(rmin >= -nghosts && rmax <= +nghosts);
-      assert(rmin == -(ORDER / 2) && rmax == +(ORDER / 2));
-      for (int i = -1; i < n + 1; ++i) {
-        if (i - i0 < rmin || i - i0 > rmax) {
-          xs[i + 1] = 0 / T(0);
-          ys[i + 1] = 0 / T(0);
-        } else {
-          T x = (i - i0) + int(CC) / T(2);
-          // T y = f(x);
-          const T dx = 1;
-          const T xlo = x - dx / 2;
-          const T xhi = x + dx / 2;
-          const T y = fint(xhi) - fint(xlo); // average of f over cell
-          xs[i + 1] = x;
-          ys[i + 1] = y;
+    }
+  } // namespace CarpetX
+};
+
+template <typename T> struct test_interp1d<CC, MINMOD, 1, T> {
+  test_interp1d() {
+    constexpr interp1d<CC, MINMOD, 1> stencil1d;
+    constexpr int nghosts = stencil1d.required_ghosts;
+    static_assert(nghosts == 1);
+    constexpr int n = 1 + 2 * nghosts;
+    constexpr int i0 = n / 2;
+    std::array<T, n + 2> xs, ys;
+
+    for (int order = 0; order <= 1; ++order) {
+      // Function f, a polynomial
+      // const auto f{[&](T x) __attribute__((__always_inline__, __flatten__))
+      // { return (order + 1) * pown(x, order); }}; Integral of f
+      // (antiderivative)
+      const auto fint = [&](T x)
+          __attribute__((__always_inline__, __flatten__)) {
+        return pown(x, order + 1);
+      };
+      std::array<T, 2> x1;
+      std::array<T, 2> y1;
+      for (int off = 0; off < 2; ++off) {
+        const T rmin = stencil1d.stencil_radius(0, off)[0];
+        const T rmax = stencil1d.stencil_radius(0, off)[1];
+        assert(rmin >= -nghosts && rmax <= +nghosts);
+        assert(rmin == -1 && rmax == +1);
+        for (int i = -1; i < n + 1; ++i) {
+          if (i - i0 < rmin || i - i0 > rmax) {
+            xs[i + 1] = 0 / T(0);
+            ys[i + 1] = 0 / T(0);
+          } else {
+            T x = (i - i0) + int(CC) / T(2);
+            // T y = f(x);
+            const T dx = 1;
+            const T xlo = x - dx / 2;
+            const T xhi = x + dx / 2;
+            const T y = fint(xhi) - fint(xlo); // average of f over cell
+            xs[i + 1] = x;
+            ys[i + 1] = y;
+          }
         }
-      }
 
-      x1[off] = int(CC) / T(4) + off / T(2);
-      y1[off] = stencil1d(
-          [&](int i) __attribute__((__always_inline__, __flatten__)) {
-            return ys[i0 + 1 + i];
-          },
-          0, off);
-      assert(isfinite(y1[off]));
-    } // for off
-    assert(y1[0] / 2 + y1[1] / 2 == ys[i0 + 1]);
-    const T dx = x1[1] - x1[0];
-    const T xlo = x1[0] - dx / 2;
-    const T xhi = x1[1] + dx / 2;
-    const T yint = fint(xhi) - fint(xlo);
-    assert(y1[0] * dx + y1[1] * dx == yint);
+        x1[off] = int(CC) / T(4) + off / T(2);
+        y1[off] = stencil1d(
+            [&](int i) __attribute__((__always_inline__, __flatten__)) {
+              return ys[i0 + 1 + i];
+            },
+            0, off);
+        assert(isfinite(y1[off]));
+      } // for off
+      // Ensure conservation
+      assert(y1[0] / 2 + y1[1] / 2 == ys[i0 + 1]);
+      const T dx = x1[1] - x1[0];
+      const T xlo = x1[0] - dx / 2;
+      const T xhi = x1[1] + dx / 2;
+      const T yint = fint(xhi) - fint(xlo);
+      assert(y1[0] * dx + y1[1] * dx == yint);
+    }
   }
-}
-}
-}
-;
+};
 
 ////////////////////////////////////////////////////////////////////////////////
 
diff --git a/CarpetX/src/prolongate_3d_rf2_impl_minmod.cxx b/CarpetX/src/prolongate_3d_rf2_impl_minmod.cxx
new file mode 100644
index 000000000..eeac6311d
--- /dev/null
+++ b/CarpetX/src/prolongate_3d_rf2_impl_minmod.cxx
@@ -0,0 +1,58 @@
+#include "prolongate_3d_rf2_impl.hxx"
+
+namespace CarpetX {
+
+// Minmod (tensor product) interpolation
+
+prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o1;
+prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o1;
+prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o1;
+prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o1;
+prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o1;
+prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o1;
+prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o1;
+prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o1;
+
+prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 3, 3, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o3;
+prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 3, 3, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o3;
+prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 3, 1, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o3;
+prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 3, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o3;
+prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 3, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o3;
+prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 3, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o3;
+prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 3, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o3;
+prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o3;
+
+prolongate_3d_rf2<VC, VC, VC, POLY, POLY, POLY, 5, 5, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c000_o5;
+prolongate_3d_rf2<VC, VC, CC, POLY, POLY, MINMOD, 5, 5, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c001_o5;
+prolongate_3d_rf2<VC, CC, VC, POLY, MINMOD, POLY, 5, 1, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c010_o5;
+prolongate_3d_rf2<VC, CC, CC, POLY, MINMOD, MINMOD, 5, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c011_o5;
+prolongate_3d_rf2<CC, VC, VC, MINMOD, POLY, POLY, 1, 5, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c100_o5;
+prolongate_3d_rf2<CC, VC, CC, MINMOD, POLY, MINMOD, 1, 5, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c101_o5;
+prolongate_3d_rf2<CC, CC, VC, MINMOD, MINMOD, POLY, 1, 1, 5, FB_NONE>
+    prolongate_minmod_3d_rf2_c110_o5;
+prolongate_3d_rf2<CC, CC, CC, MINMOD, MINMOD, MINMOD, 1, 1, 1, FB_NONE>
+    prolongate_minmod_3d_rf2_c111_o5;
+
+} // namespace CarpetX