diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx
index 7bd0c175e2..f70f349b2c 100644
--- a/src/mesh/difops.cxx
+++ b/src/mesh/difops.cxx
@@ -44,6 +44,7 @@
 #include <bout/solver.hxx>
 #include <bout/unused.hxx>
 #include <bout/utils.hxx>
+#include <bout/yboundary_regions.hxx>
 
 #include <cmath>
 #include <limits>
@@ -386,23 +387,32 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin,
     Field3D result{zeroFrom(fin)};
     flow_ylow = zeroFrom(fin);
 
+    const auto yboundary = coord->getYBoundary();
+
     BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) {
       const auto iyp = i.yp();
       const auto iym = i.ym();
 
       // Upper cell edge
-      const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary
-      const BoutReal gradient_up =
-          (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(coord->g_22_yhigh()[i]));
+      BoutReal flux_up = 0;
+      if (bndry_flux or not yboundary.contains<+1>(i)) {
+        const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary
+
+        const BoutReal gradient_up =
+            (f_up[iyp] - fin[i]) / (coord->dy()[i] * sqrt(coord->g_22_yhigh()[i]));
 
-      const BoutReal flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i];
+        flux_up = c_up * gradient_up * coord->cell_area_yhigh()[i];
+      }
 
       // Lower cell edge
-      const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary
-      const BoutReal gradient_down =
-          (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(coord->g_22_ylow()[i]));
+      BoutReal flux_down = 0;
+      if (bndry_flux or not yboundary.contains<-1>(i)) {
+        const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary
+        const BoutReal gradient_down =
+            (fin[i] - f_down[iym]) / (coord->dy()[i] * sqrt(coord->g_22_ylow()[i]));
 
-      const BoutReal flux_down = c_down * gradient_down * coord->cell_area_ylow()[i];
+        flux_down = c_down * gradient_down * coord->cell_area_ylow()[i];
+      }
 
       // Add the fluxes
       result[i] = (flux_up - flux_down) / (coord->cell_volume()[i]);
diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx
index d00187fedc..afe8088523 100644
--- a/src/mesh/fv_ops.cxx
+++ b/src/mesh/fv_ops.cxx
@@ -12,6 +12,7 @@
 #include "bout/msg_stack.hxx"
 #include "bout/region.hxx"
 #include "bout/utils.hxx"
+#include <bout/yboundary_regions.hxx>
 
 namespace {
 template <class T>
@@ -211,8 +212,12 @@ Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_fl
     const auto iyp = i.yp();
     const auto iym = i.ym();
 
-    if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x())
-        || (i.y() != mesh->yend)) {
+    const auto yboundary = coord->getYBoundary();
+
+    if (bndry_flux
+        || (not K.isFci()
+            and (mesh->periodicY(i.x()) || !mesh->lastY(i.x()) || (i.y() != mesh->yend)))
+        or (K.isFci() and yboundary.contains<+1>(i))) {
 
       const BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary
       const BoutReal J = 0.5 * (coord->J()[i] + coord->J()[iyp]); // Jacobian at boundary
@@ -227,8 +232,11 @@ Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_fl
     }
 
     // Calculate flux at lower surface
-    if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x())
-        || (i.y() != mesh->ystart)) {
+    if (bndry_flux
+        || (not K.isFci()
+            and (mesh->periodicY(i.x()) || !mesh->firstY(i.x())
+                 || (i.y() != mesh->ystart)))
+        or (K.isFci() and yboundary.contains<-1>(i))) {
       const BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary
       const BoutReal J = 0.5 * (coord->J()[i] + coord->J()[iym]); // Jacobian at boundary