diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx
index 6d803ab28a..a21f62cd4d 100644
--- a/include/bout/coordinates.hxx
+++ b/include/bout/coordinates.hxx
@@ -116,101 +116,115 @@ public:
   FieldMetric& g_22_yhigh();
   // Cell Areas
   const FieldMetric& cell_area_xlow() const {
-    if (!_cell_area_xlow.has_value()) {
-      _compute_cell_area_x();
+    if (_cell_area_xlow.has_value()) {
+      return *_cell_area_xlow;
     }
+    _compute_cell_area_x();
     ASSERT2(_cell_area_xlow.has_value());
     return *_cell_area_xlow;
   }
   const FieldMetric& cell_area_xhigh() const {
-    if (!_cell_area_xhigh.has_value()) {
-      _compute_cell_area_x();
+    if (_cell_area_xhigh.has_value()) {
+      return *_cell_area_xhigh;
     }
+    _compute_cell_area_x();
     ASSERT2(_cell_area_xhigh.has_value());
     return *_cell_area_xhigh;
   }
   const FieldMetric& cell_area_ylow() const {
-    if (!_cell_area_ylow.has_value()) {
-      _compute_cell_area_y();
+    if (_cell_area_ylow.has_value()) {
+      return *_cell_area_ylow;
     }
+    _compute_cell_area_y();
     ASSERT2(_cell_area_ylow.has_value());
     return *_cell_area_ylow;
   }
   const FieldMetric& cell_area_yhigh() const {
-    if (!_cell_area_yhigh.has_value()) {
-      _compute_cell_area_y();
+    if (_cell_area_yhigh.has_value()) {
+      return *_cell_area_yhigh;
     }
+    _compute_cell_area_y();
     ASSERT2(_cell_area_yhigh.has_value());
     return *_cell_area_yhigh;
   }
   const FieldMetric& cell_area_zlow() const {
-    if (!_cell_area_zlow.has_value()) {
-      _compute_cell_area_z();
+    if (_cell_area_zlow.has_value()) {
+      return *_cell_area_zlow;
     }
+    _compute_cell_area_z();
     ASSERT2(_cell_area_zlow.has_value());
     return *_cell_area_zlow;
   }
   const FieldMetric& cell_area_zhigh() const {
-    if (!_cell_area_zhigh.has_value()) {
-      _compute_cell_area_z();
+    if (_cell_area_zhigh.has_value()) {
+      return *_cell_area_zhigh;
     }
+    _compute_cell_area_z();
     ASSERT2(_cell_area_zhigh.has_value());
     return *_cell_area_zhigh;
   }
   FieldMetric& cell_area_xlow() {
-    if (!_cell_area_xlow.has_value()) {
-      _compute_cell_area_x();
+    if (_cell_area_xlow.has_value()) {
+      return *_cell_area_xlow;
     }
+    _compute_cell_area_x();
     ASSERT2(_cell_area_xlow.has_value());
     return *_cell_area_xlow;
   }
   FieldMetric& cell_area_xhigh() {
-    if (!_cell_area_xhigh.has_value()) {
-      _compute_cell_area_x();
+    if (_cell_area_xhigh.has_value()) {
+      return *_cell_area_xhigh;
     }
+    _compute_cell_area_x();
     ASSERT2(_cell_area_xhigh.has_value());
     return *_cell_area_xhigh;
   }
   FieldMetric& cell_area_ylow() {
-    if (!_cell_area_ylow.has_value()) {
-      _compute_cell_area_y();
+    if (_cell_area_ylow.has_value()) {
+      return *_cell_area_ylow;
     }
+    _compute_cell_area_y();
     ASSERT2(_cell_area_ylow.has_value());
     return *_cell_area_ylow;
   }
   FieldMetric& cell_area_yhigh() {
-    if (!_cell_area_yhigh.has_value()) {
-      _compute_cell_area_y();
+    if (_cell_area_yhigh.has_value()) {
+      return *_cell_area_yhigh;
     }
+    _compute_cell_area_y();
     ASSERT2(_cell_area_yhigh.has_value());
     return *_cell_area_yhigh;
   }
   FieldMetric& cell_area_zlow() {
-    if (!_cell_area_zlow.has_value()) {
-      _compute_cell_area_z();
+    if (_cell_area_zlow.has_value()) {
+      return *_cell_area_zlow;
     }
+    _compute_cell_area_z();
     ASSERT2(_cell_area_zlow.has_value());
     return *_cell_area_zlow;
   }
   FieldMetric& cell_area_zhigh() {
-    if (!_cell_area_zhigh.has_value()) {
-      _compute_cell_area_z();
+    if (_cell_area_zhigh.has_value()) {
+      return *_cell_area_zhigh;
     }
+    _compute_cell_area_z();
     ASSERT2(_cell_area_zhigh.has_value());
     return *_cell_area_zhigh;
   }
   // Cell Volume
   const FieldMetric& cell_volume() const {
-    if (!_cell_volume.has_value()) {
-      _compute_cell_volume();
+    if (_cell_volume.has_value()) {
+      return *_cell_volume;
     }
+    _compute_cell_volume();
     ASSERT2(_cell_volume.has_value());
     return *_cell_volume;
   }
   FieldMetric& cell_volume() {
-    if (!_cell_volume.has_value()) {
-      _compute_cell_volume();
+    if (_cell_volume.has_value()) {
+      return *_cell_volume;
     }
+    _compute_cell_volume();
     ASSERT2(_cell_volume.has_value());
     return *_cell_volume;
   }
diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx
index a8d9279378..14bc5c2824 100644
--- a/include/bout/derivs.hxx
+++ b/include/bout/derivs.hxx
@@ -35,6 +35,7 @@
 #include "bout/vector3d.hxx"
 
 #include "bout/bout_types.hxx"
+#include <type_traits>
 
 ////////// FIRST DERIVATIVES //////////
 
@@ -86,6 +87,13 @@ Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
             const std::string& method = "DEFAULT",
             const std::string& region = "RGN_NOBNDRY");
 
+template <typename E>
+std::enable_if_t<is_expr_field3d_v<E> && !bout::utils::is_Field3D_v<E>, const Field3D>
+DDY(const E& expr, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT",
+    const std::string& region = "RGN_NOBNDRY") {
+  return DDY(Field3D{expr}, outloc, method, region);
+}
+
 /// Calculate first partial derivative in Y
 ///
 ///   \f$\partial / \partial y\f$
diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 81d4c67c25..67efdd6476 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -334,41 +334,24 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp
       // Pre-calculate factors which multiply fluxes
 #if not(BOUT_USE_METRIC_3D)
       // For right cell boundaries
-      BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1))
-                               / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1)));
-
-      const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j));
-      const BoutReal flux_factor_rp =
-          common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1));
-
+      const BoutReal area_r = coord->cell_area_yhigh()(i, j);
+      const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j);
+      const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1);
       // For left cell boundaries
-      common_factor = (coord->J(i, j) + coord->J(i, j - 1))
-                      / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1)));
-
-      const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j));
-      const BoutReal flux_factor_lm =
-          common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1));
+      const BoutReal area_l = coord->cell_area_ylow()(i, j);
+      const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j);
+      const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1);
 #endif
       for (int k = mesh->zstart; k <= mesh->zend; k++) {
 #if BOUT_USE_METRIC_3D
         // 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)));
-
-        BoutReal flux_factor_rc =
-            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
-        BoutReal flux_factor_rp =
-            common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));
-
+        const BoutReal area_r = coord->cell_area_yhigh()(i, j, k);
+        const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k);
+        const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, 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)));
-
-        BoutReal flux_factor_lc =
-            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
-        BoutReal flux_factor_lm =
-            common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));
+        const BoutReal area_l = coord->cell_area_ylow()(i, j, k);
+        const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k);
+        const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k);
 #endif
 
         ////////////////////////////////////////////
@@ -656,12 +639,10 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
       const auto iym = i.ym();
 
       result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp])
-                       * (coord->J[i] + coord->J.yup()[iyp])
-                       / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp]))
+                       * coord->cell_area_yhigh()[i]
                    - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym])
-                         * (coord->J[i] + coord->J.ydown()[iym])
-                         / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym])))
-                  / (coord->dy[i] * coord->J[i]);
+                         * coord->cell_area_ylow()[i])
+                  / coord->cell_volume()[i];
     }
     return result;
   }
@@ -702,57 +683,30 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
       // Pre-calculate factors which multiply fluxes
 #if not(BOUT_USE_METRIC_3D)
       // For right cell boundaries
-      const BoutReal common_factor_r =
-          (coord->J(i, j) + coord->J(i, j + 1))
-          / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1)));
+      const BoutReal area_rp = coord->cell_area_yhigh()(i, j);
 
-      const BoutReal flux_factor_rc =
-          common_factor_r / (coord->dy(i, j) * coord->J(i, j));
-      const BoutReal flux_factor_rp =
-          common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1));
-
-      const BoutReal area_rp =
-          common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1);
+      const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j);
+      const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1);
 
       // For left cell boundaries
-      const BoutReal common_factor_l =
-          (coord->J(i, j) + coord->J(i, j - 1))
-          / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1)));
-
-      const BoutReal flux_factor_lc =
-          common_factor_l / (coord->dy(i, j) * coord->J(i, j));
-      const BoutReal flux_factor_lm =
-          common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1));
+      const BoutReal area_lc = coord->cell_area_ylow()(i, j);
 
-      const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j);
+      const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j);
+      const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1);
 #endif
       for (int k = 0; k < mesh->LocalNz; k++) {
 #if BOUT_USE_METRIC_3D
         // For right cell boundaries
-        const BoutReal common_factor_r =
-            (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 flux_factor_rc =
-            common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal flux_factor_rp =
-            common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));
+        const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k);
 
-        const BoutReal area_rp =
-            common_factor_r * 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);
+        const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1, k);
 
         // For left cell boundaries
-        const BoutReal common_factor_l =
-            (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_l / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal flux_factor_lm =
-            common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));
-
-        const BoutReal area_lc =
-            common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k);
+        const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k);
+        const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1, k);
 #endif
 
         ////////////////////////////////////////////
@@ -937,24 +891,16 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
 
       for (int k = 0; k < mesh->LocalNz; k++) {
         // For right cell boundaries
-        const BoutReal common_factor_r =
-            (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_r = coord->cell_area_yhigh()(i, j, k);
 
-        const BoutReal flux_factor_rc =
-            common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal flux_factor_rp =
-            common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k));
+        const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k);
+        const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k);
 
         // For left cell boundaries
-        const BoutReal common_factor_l =
-            (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 flux_factor_lc =
-            common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k));
-        const BoutReal flux_factor_lm =
-            common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k));
+        const BoutReal area_l = coord->cell_area_ylow()(i, j, k);
+
+        const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k);
+        const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k);
 
         ////////////////////////////////////////////
         // Reconstruct f at the cell faces
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 057ffa65b5..afb63f981d 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -2015,17 +2015,23 @@ const Coordinates::FieldMetric& Coordinates::g_22_ylow() const {
   if (_g_22_ylow.has_value()) {
     return *_g_22_ylow;
   }
-  _g_22_ylow.emplace(emptyFrom(g_22));
-  //_g_22_ylow->setLocation(CELL_YLOW);
-  auto* mesh = Bxy.getMesh();
-  if (Bxy.isFci()) {
-    if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) {
-      throw BoutException("The grid file does not contain `g_22_cell_ylow`.");
-    }
-  } else {
-    ASSERT0(mesh->ystart > 0);
-    BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
-      _g_22_ylow.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()])));
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_g_22_ylow.has_value()) {
+      _g_22_ylow.emplace(emptyFrom(g_22));
+      //_g_22_ylow->setLocation(CELL_YLOW);
+      auto* mesh = Bxy.getMesh();
+      if (Bxy.isFci()) {
+        if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) {
+          throw BoutException("The grid file does not contain `g_22_cell_ylow`.");
+        }
+      } else {
+        ASSERT0(mesh->ystart > 0);
+        BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
+          _g_22_ylow.value()[i] =
+              SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()])));
+        }
+      }
     }
   }
   return g_22_ylow();
@@ -2035,99 +2041,127 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const {
   if (_g_22_yhigh.has_value()) {
     return *_g_22_yhigh;
   }
-  _g_22_yhigh.emplace(emptyFrom(g_22));
-  auto* mesh = Bxy.getMesh();
-  if (Bxy.isFci()) {
-    if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) {
-      throw BoutException("The grid file does not contain `g_22_cell_yhigh`.");
-    }
-  } else {
-    ASSERT0(mesh->ystart > 0);
-    BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
-      _g_22_yhigh.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()])));
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_g_22_yhigh.has_value()) {
+      _g_22_yhigh.emplace(emptyFrom(g_22));
+      auto* mesh = Bxy.getMesh();
+      if (Bxy.isFci()) {
+        if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) {
+          throw BoutException("The grid file does not contain `g_22_cell_yhigh`.");
+        }
+      } else {
+        ASSERT0(mesh->ystart > 0);
+        BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
+          _g_22_yhigh.value()[i] =
+              SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()])));
+        }
+      }
     }
   }
   return g_22_yhigh();
 }
 
 void Coordinates::_compute_cell_area_x() const {
-  const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz;
-  _cell_area_xlow.emplace(emptyFrom(area_centre));
-  _cell_area_xhigh.emplace(emptyFrom(area_centre));
-  // We cannot setLocation, as that would trigger the computation of staggered
-  // metrics.
-  auto* mesh = Bxy.getMesh();
-  ASSERT0(mesh->xstart > 0);
-  BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) {
-    (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]);
-    (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]);
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_cell_area_xlow.has_value()) {
+      const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz;
+      _cell_area_xlow.emplace(emptyFrom(area_centre));
+      _cell_area_xhigh.emplace(emptyFrom(area_centre));
+      // We cannot setLocation, as that would trigger the computation of staggered
+      // metrics.
+      auto* mesh = Bxy.getMesh();
+      ASSERT0(mesh->xstart > 0);
+      BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) {
+        (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]);
+        (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]);
+      }
+    }
   }
 }
 
 void Coordinates::_compute_cell_area_y() const {
-  auto* mesh = Bxy.getMesh();
-  if (g_11.isFci()) {
-    const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13));
-    auto jxz_ylow = emptyFrom(jxz_centre);
-    auto jxz_yhigh = emptyFrom(jxz_centre);
-
-    auto By_c = emptyFrom(jxz_centre);
-    auto By_h = emptyFrom(jxz_yhigh);
-    auto By_l = emptyFrom(jxz_ylow);
-    if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) {
-      throw BoutException("The grid file does not contain `By`.");
-    }
-    if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) {
-      throw BoutException("The grid file does not contain `By_cell_ylow`.");
-    }
-    if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) {
-      throw BoutException("The grid file does not contain `By_cell_yhigh`.");
-    }
-    BOUT_FOR(i, By_c.getRegion("RGN_NOY")) {
-      jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i];
-      jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i];
-    }
-    ASSERT3(isUniform(dx, true, "RGN_ALL"));
-    ASSERT2(isUniform(dx, false, "RGN_ALL"));
-    ASSERT3(isUniform(dz, true, "RGN_ALL"));
-    ASSERT2(isUniform(dz, false, "RGN_ALL"));
-    _cell_area_ylow.emplace(jxz_ylow * dx * dz);
-    _cell_area_yhigh.emplace(jxz_yhigh * dx * dz);
-  } else {
-    // Field aligned
-    const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz;
-    _cell_area_ylow.emplace(emptyFrom(area_centre));
-    _cell_area_yhigh.emplace(emptyFrom(area_centre));
-    // We cannot setLocation, as that would trigger the computation of staggered
-    // metrics.
-    BOUT_FOR(i, mesh->getRegion("RGN_ALL")) {
-      if (i.y() > 0) {
-        (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]);
-      } else {
-        (*_cell_area_ylow)[i] = BoutNaN;
-      }
-      if (i.y() < mesh->LocalNy - 1) {
-        (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]);
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_cell_area_ylow.has_value()) {
+      auto* mesh = Bxy.getMesh();
+      if (g_11.isFci()) {
+        const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13));
+        auto jxz_ylow = emptyFrom(jxz_centre);
+        auto jxz_yhigh = emptyFrom(jxz_centre);
+
+        auto By_c = emptyFrom(jxz_centre);
+        auto By_h = emptyFrom(jxz_yhigh);
+        auto By_l = emptyFrom(jxz_ylow);
+        if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) {
+          throw BoutException("The grid file does not contain `By`.");
+        }
+        if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) {
+          throw BoutException("The grid file does not contain `By_cell_ylow`.");
+        }
+        if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) {
+          throw BoutException("The grid file does not contain `By_cell_yhigh`.");
+        }
+        BOUT_FOR(i, By_c.getRegion("RGN_NOY")) {
+          jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i];
+          jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i];
+        }
+        ASSERT3(isUniform(dx, true, "RGN_ALL"));
+        ASSERT2(isUniform(dx, false, "RGN_ALL"));
+        ASSERT3(isUniform(dz, true, "RGN_ALL"));
+        ASSERT2(isUniform(dz, false, "RGN_ALL"));
+        _cell_area_ylow.emplace(jxz_ylow * dx * dz);
+        _cell_area_yhigh.emplace(jxz_yhigh * dx * dz);
       } else {
-        (*_cell_area_yhigh)[i] = BoutNaN;
+        // Field aligned
+        const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz;
+        _cell_area_ylow.emplace(emptyFrom(area_centre));
+        _cell_area_yhigh.emplace(emptyFrom(area_centre));
+        // We cannot setLocation, as that would trigger the computation of staggered
+        // metrics.
+        BOUT_FOR(i, mesh->getRegion("RGN_ALL")) {
+          if (i.y() > 0) {
+            (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]);
+          } else {
+            (*_cell_area_ylow)[i] = BoutNaN;
+          }
+          if (i.y() < mesh->LocalNy - 1) {
+            (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]);
+          } else {
+            (*_cell_area_yhigh)[i] = BoutNaN;
+          }
+        }
       }
     }
   }
 }
 
 void Coordinates::_compute_cell_area_z() const {
-  const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy;
-  _cell_area_zlow.emplace(emptyFrom(area_centre));
-  _cell_area_zhigh.emplace(emptyFrom(area_centre));
-  // We cannot setLocation, as that would trigger the computation of staggered
-  // metrics.
-  BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) {
-    (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]);
-    (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]);
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_cell_volume.has_value()) {
+      const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy;
+      _cell_area_zlow.emplace(emptyFrom(area_centre));
+      _cell_area_zhigh.emplace(emptyFrom(area_centre));
+      // We cannot setLocation, as that would trigger the computation of staggered
+      // metrics.
+      BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) {
+        (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]);
+        (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]);
+      }
+    }
   }
 }
 
-void Coordinates::_compute_cell_volume() const { _cell_volume.emplace(J * dx * dy * dz); }
+void Coordinates::_compute_cell_volume() const {
+  BOUT_OMP_SAFE(critical)
+  {
+    if (!_cell_volume.has_value()) {
+      _cell_volume.emplace(J * dx * dy * dz);
+    }
+  }
+}
 
 std::shared_ptr<YBoundary> Coordinates::makeYBoundary(YBndryType type) const {
   return std::make_shared<YBoundary>(type, localoptions, *localmesh);