diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index a4e8f908bf..1b1379d5ad 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -1108,22 +1108,14 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
         // so calculate inside the k loop.
 
         // For right cell boundaries
-        BoutReal common_factor =
-            (coord->J(i, j, k) + coord->J(i, j + 1, k))
-            / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k)));
+        const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k);
 
-        const BoutReal flux_factor_rc =
-            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal area_rp =
-            common_factor * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k);
+        const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j, k);
 
         // For left cell boundaries
-        common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k))
-                        / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k)));
+        const BoutReal area_lc = coord->cell_area_ylow()(i, j, k);
 
-        const BoutReal flux_factor_lc =
-            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal area_lc = common_factor * coord->dx(i, j, k) * coord->dz(i, j, k);
+        const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k);
 
         ////////////////////////////////////////////
         // Reconstruct f at the cell faces
@@ -1185,7 +1177,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           result(i, j, k) += (actual_ke - expected_ke) * flux_factor_rc;
 
           // Final flow through boundary is the expected value
-          flow_ylow(i, j + 1, k) += expected_ke * area_rp; //expected_ke * area_rp;
+          flow_ylow(i, j + 1, k) += expected_ke * area_rp;
 
         } else {
           // Maximum wave speed in the two cells