diff --git a/examples/6field-simple/elm_6f.cxx b/examples/6field-simple/elm_6f.cxx
index de5d52c47e..7afdb6a31b 100644
--- a/examples/6field-simple/elm_6f.cxx
+++ b/examples/6field-simple/elm_6f.cxx
@@ -36,7 +36,8 @@ class Elm_6f : public PhysicsModel {
   // 2D inital profiles
   /// Current and pressure
   Field2D J0, P0;
-
+  /// Curvature term
+  Vector2D b0xcv;
   /// When diamagnetic terms used
   Field2D phi0;
 
@@ -388,6 +389,10 @@ class Elm_6f : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     //////////////////////////////////////////////////////////////
     // Read parameters from the options file
     //
@@ -663,10 +668,30 @@ class Elm_6f : public PhysicsModel {
     phi_curv = options["phi_curv"].doc("Compressional ExB terms").withDefault(true);
     g = options["gamma"].doc("Ratio of specific heats").withDefault(5.0 / 3.0);
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (not mesh->IncIntShear) {
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
+    }
+
     //////////////////////////////////////////////////////////////
     // NORMALISE QUANTITIES
 
@@ -812,6 +837,10 @@ class Elm_6f : public PhysicsModel {
 
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
     if ((!T0_fake_prof) && n0_fake_prof) {
       N0 = N0tanh(n0_height * Nbar, n0_ave * Nbar, n0_width, n0_center, n0_bottom_x);
 
@@ -999,7 +1028,7 @@ class Elm_6f : public PhysicsModel {
     }
 
     /**************** CALCULATE METRICS ******************/
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
     
     //////////////////////////////////////////////////////////////
     // SHIFTED RADIAL COORDINATES
@@ -1358,7 +1387,7 @@ class Elm_6f : public PhysicsModel {
 
       ddt(U) = -SQ(tokamak_coordinates_factory.get_Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.get_Bxy(); // Grad j term
 
-      ddt(U) += 2.0 * Upara1 * tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+      ddt(U) += 2.0 * Upara1 * b0xcv * Grad(P); // curvature term
 
       ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar); // b dot grad j
 
diff --git a/examples/conducting-wall-mode/cwm.cxx b/examples/conducting-wall-mode/cwm.cxx
index 98637059e0..a3b1720da7 100644
--- a/examples/conducting-wall-mode/cwm.cxx
+++ b/examples/conducting-wall-mode/cwm.cxx
@@ -130,7 +130,7 @@ class CWM : public PhysicsModel {
 
     // Normalise geometry
     tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     // Set nu
     nu = nu_hat * Ni0 / pow(Te0, 1.5);
diff --git a/examples/constraints/alfven-wave/alfven.cxx b/examples/constraints/alfven-wave/alfven.cxx
index 3ac5ecb067..2b9b2beb91 100644
--- a/examples/constraints/alfven-wave/alfven.cxx
+++ b/examples/constraints/alfven-wave/alfven.cxx
@@ -172,7 +172,7 @@ class Alfven : public PhysicsModel {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
     tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
   }
 };
 
diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx
index 2f800ad825..8bba66d656 100644
--- a/examples/dalf3/dalf3.cxx
+++ b/examples/dalf3/dalf3.cxx
@@ -53,6 +53,7 @@ class DALF3 : public PhysicsModel {
   Field3D phi, apar, jpar;
 
   Field2D B0, Pe0, Jpar0;
+  Vector2D b0xcv;
 
   Field2D eta; // Collisional damping (resistivity)
   BoutReal beta_hat, mu_hat;
@@ -99,6 +100,10 @@ class DALF3 : public PhysicsModel {
     Pe0 = 2. * Charge * Ni0 * Te0; // Electron pressure in Pascals
     SAVE_ONCE(Pe0);
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     //////////////////////////////////////////////////////////////
     // Options
 
@@ -156,6 +161,8 @@ class DALF3 : public PhysicsModel {
         Options::root()["mesh"]["paralleltransform"]["type"].withDefault("identity");
 
     if (lowercase(ptstr) == "shifted") {
+      // Dimits style, using local coordinate system
+      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
       noshear = true;
     }
 
@@ -210,9 +217,13 @@ class DALF3 : public PhysicsModel {
     hyper_viscosity /= wci * SQ(SQ(rho_s));
     viscosity_par /= wci * SQ(rho_s);
 
+    b0xcv.x /= Bnorm;
+    b0xcv.y *= rho_s * rho_s;
+    b0xcv.z *= rho_s * rho_s;
+
     // Metrics
     tokamak_coordinates_factory.normalise(rho_s, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     SOLVE_FOR3(Vort, Pe, Vpar);
     comms.add(Vort, Pe, Vpar);
@@ -264,7 +275,7 @@ class DALF3 : public PhysicsModel {
   Field3D Kappa(const Field3D& f) {
     if (curv_kappa) {
       // Use the b0xcv vector from grid file
-      return -2. * tokamak_coordinates_factory.get_b0xcv() * Grad(f) / B0;
+      return -2. * b0xcv * Grad(f) / B0;
     }
 
     return 2. * bracket(log(B0), f, bm);
diff --git a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
index 05fb9363e0..7edfbe8e38 100644
--- a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
+++ b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
@@ -96,6 +96,7 @@ class ELMpb : public PhysicsModel {
 private:
   // 2D inital profiles
   Field2D J0, P0; // Current and pressure
+  Vector2D b0xcv; // Curvature term
   Field2D beta;   // Used for Vpar terms
   Coordinates::FieldMetric gradparB;
   Field2D phi0; // When diamagnetic terms used
@@ -368,6 +369,10 @@ class ELMpb : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     mesh->get(Psixy, "psixy");        // get Psi
     mesh->get(Psiaxis, "psi_axis");   // axis flux
     mesh->get(Psibndry, "psi_bndry"); // edge flux
@@ -712,7 +717,7 @@ class ELMpb : public PhysicsModel {
       Lbar = 1.0;
     }
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     V0 = -tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Bpxy() * Dphi0 / tokamak_coordinates_factory.get_Bxy();
 
@@ -799,10 +804,30 @@ class ELMpb : public PhysicsModel {
       }
     }
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (not mesh->IncIntShear) {
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
+    }
+
     //////////////////////////////////////////////////////////////
     // NORMALISE QUANTITIES
 
@@ -914,6 +939,11 @@ class ELMpb : public PhysicsModel {
     V0 = V0 / Va;
     Dphi0 *= Tbar;
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
+
     if (constn0) {
       T0_fake_prof = false;
       n0_fake_prof = false;
@@ -1690,7 +1720,7 @@ class ELMpb : public PhysicsModel {
     //   ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     // }
 
-    ddt(U) += tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+    ddt(U) += b0xcv * Grad(P); // curvature term
 
     // if (!nogradparj) { // Parallel current term
     //   ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
@@ -1866,7 +1896,7 @@ class ELMpb : public PhysicsModel {
       ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE);
 
       if (phi_curv) {
-        ddt(P) -= 2. * beta * tokamak_coordinates_factory.get_b0xcv() * Grad(phi);
+        ddt(P) -= 2. * beta * b0xcv * Grad(phi);
       }
 
       // Vpar equation
@@ -1966,7 +1996,7 @@ class ELMpb : public PhysicsModel {
     Coordinates* metric = mesh->getCoordinates();
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * metric->Bxy() * metric->Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * tokamak_coordinates_factory.get_b0xcv()) * Grad(P);
+    U1 += (gamma * metric->Bxy() * metric->Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};
@@ -2025,7 +2055,7 @@ class ELMpb : public PhysicsModel {
     JPsi.setBoundary("Psi");
     JPsi.applyBoundary();
 
-    Field3D JU = tokamak_coordinates_factory.get_b0xcv() * Grad(ddt(P)) - SQ(metric->Bxy()) * Grad_par(Jpar, CELL_CENTRE)
+    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(metric->Bxy()) * Grad_par(Jpar, CELL_CENTRE)
                  + SQ(metric->Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
     JU.setBoundary("U");
     JU.applyBoundary();
diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx
index 4f19fb573c..2b4603e9c8 100644
--- a/examples/elm-pb/elm_pb.cxx
+++ b/examples/elm-pb/elm_pb.cxx
@@ -38,6 +38,7 @@ class ELMpb : public PhysicsModel {
 private:
   // 2D inital profiles
   Field2D J0, P0; // Current and pressure
+  Vector2D b0xcv; // Curvature term
   Field2D beta;   // Used for Vpar terms
   Coordinates::FieldMetric gradparB;
   Field2D phi0; // When diamagnetic terms used
@@ -310,6 +311,10 @@ class ELMpb : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     mesh->get(Psixy, "psixy");        // get Psi
     mesh->get(Psiaxis, "psi_axis");   // axis flux
     mesh->get(Psibndry, "psi_bndry"); // edge flux
@@ -725,13 +730,28 @@ class ELMpb : public PhysicsModel {
       }
     }
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
-    // Dimits style, using local coordinate system
-    if (include_curvature and !mesh->IncIntShear) {
-      noshear = true;
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (not mesh->IncIntShear) {
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      }
     }
 
     //////////////////////////////////////////////////////////////
@@ -852,8 +872,12 @@ class ELMpb : public PhysicsModel {
     V0 = V0 / Va;
     Dphi0 *= Tbar;
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
 
     //////////////////////////////////////////////////////////////
@@ -1499,7 +1523,7 @@ class ELMpb : public PhysicsModel {
       ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     }
 
-    ddt(U) += tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+    ddt(U) += b0xcv * Grad(P); // curvature term
 
     if (!nogradparj) {                                 // Parallel current term
       ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
@@ -1673,7 +1697,7 @@ class ELMpb : public PhysicsModel {
       ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE);
 
       if (phi_curv) {
-        ddt(P) -= 2. * beta * tokamak_coordinates_factory.get_b0xcv() * Grad(phi);
+        ddt(P) -= 2. * beta * b0xcv * Grad(phi);
       }
 
       // Vpar equation
@@ -1772,7 +1796,7 @@ class ELMpb : public PhysicsModel {
     mesh->communicate(Jrhs, ddt(P));
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * tokamak_coordinates_factory.get_b0xcv()) * Grad(P);
+    U1 += (gamma * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};
@@ -1829,7 +1853,7 @@ class ELMpb : public PhysicsModel {
     JPsi.setBoundary("Psi");
     JPsi.applyBoundary();
 
-    Field3D JU = tokamak_coordinates_factory.get_b0xcv() * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jpar, CELL_CENTRE)
+    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jpar, CELL_CENTRE)
                  + SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
     JU.setBoundary("U");
     JU.applyBoundary();
diff --git a/examples/gyro-gem/gem.cxx b/examples/gyro-gem/gem.cxx
index f577254e94..111898d4fa 100644
--- a/examples/gyro-gem/gem.cxx
+++ b/examples/gyro-gem/gem.cxx
@@ -253,7 +253,7 @@ class GEM : public PhysicsModel {
     SAVE_ONCE(Bbar);
 
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
 
     beta_e = 4.e-7 * PI * max(p_e, true) / (Bbar * Bbar);
     SAVE_ONCE(beta_e);
diff --git a/examples/laplacexy/alfven-wave/alfven.cxx b/examples/laplacexy/alfven-wave/alfven.cxx
index 617e7bee31..8b1b4a836e 100644
--- a/examples/laplacexy/alfven-wave/alfven.cxx
+++ b/examples/laplacexy/alfven-wave/alfven.cxx
@@ -174,7 +174,7 @@ class Alfven : public PhysicsModel {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
     tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
   }
 };
 
diff --git a/examples/laplacexy/laplace_perp/test.cxx b/examples/laplacexy/laplace_perp/test.cxx
index 5cd024650e..7d3c8972c4 100644
--- a/examples/laplacexy/laplace_perp/test.cxx
+++ b/examples/laplacexy/laplace_perp/test.cxx
@@ -15,7 +15,7 @@ int main(int argc, char** argv) {
   calc_metric = Options::root()["calc_metric"].withDefault(false);
   if (calc_metric) {
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
   }
 
   ///////////////////////////////////////
diff --git a/examples/wave-slab/wave_slab.cxx b/examples/wave-slab/wave_slab.cxx
index e1f1cf2c5f..6ac11568e3 100644
--- a/examples/wave-slab/wave_slab.cxx
+++ b/examples/wave-slab/wave_slab.cxx
@@ -18,7 +18,7 @@ class WaveTest : public PhysicsModel {
   int init(bool UNUSED(restarting)) {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
-    const auto& coords = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coords = tokamak_coordinates_factory.make_tokamak_coordinates(true);
 
     solver->add(f, "f");
     solver->add(g, "g");
diff --git a/include/bout/tokamak_coordinates_factory.hxx b/include/bout/tokamak_coordinates_factory.hxx
index cfa2706c62..5cc9765153 100644
--- a/include/bout/tokamak_coordinates_factory.hxx
+++ b/include/bout/tokamak_coordinates_factory.hxx
@@ -22,7 +22,6 @@ private:
   Field2D hthe_m;
   FieldMetric ShearFactor_m;
   FieldMetric dx_m;
-  Vector2D b0xcv_m; // Curvature term
 
 
 public:
@@ -38,10 +37,6 @@ public:
     mesh.get(hthe_m, "hthe");
     mesh.get(ShearFactor_m, "sinty");
     mesh.get(dx_m, "dpsi");
-
-    // Load magnetic curvature term
-    b0xcv_m.covariant = false;  // Read contravariant components
-    mesh.get(b0xcv_m, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
   }
 
   BoutReal get_sign_of_bp() {
@@ -51,11 +46,13 @@ public:
     return 1.0;
   }
 
-  Coordinates* make_tokamak_coordinates(const bool noshear, const bool include_curvature) {
+  Coordinates* make_tokamak_coordinates(const bool noshear) {
 
     const BoutReal sign_of_bp = get_sign_of_bp();
 
-    set_shearfactor_and_curvature_term(noshear, include_curvature);
+    if (noshear) {
+      ShearFactor_m = 0.0;
+    }
 
     auto* coord = mesh_m.getCoordinates();
 
@@ -83,28 +80,6 @@ public:
     return coord;
   }
 
-  void set_shearfactor_and_curvature_term(const bool noshear, const bool include_curvature) {
-
-    if (!include_curvature) {
-      b0xcv_m = 0.0;
-    }
-
-    if (noshear) {
-      if (include_curvature) {
-        b0xcv_m.z += ShearFactor_m * b0xcv_m.x;
-      }
-    }
-
-    //      if (ShiftXderivs and not!mesh->IncIntShear) {
-    //        // Dimits style, using local coordinate system
-    //        if (include_curvature) {
-    //          b0xcv.z += I * b0xcv.x;
-    //        }
-    //        I = 0.0; // I disappears from metric
-    //      }
-
-  }
-
   void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor=1.0) {
 
     Rxy_m /= Lbar;
@@ -114,18 +89,12 @@ public:
     hthe_m /= Lbar;
     ShearFactor_m *= Lbar * Lbar * Bbar * ShearFactor;
     dx_m /= Lbar * Lbar * Bbar;
-
-    // Normalise curvature term
-    b0xcv_m.x /= Bbar;
-    b0xcv_m.y *= Lbar * Lbar;
-    b0xcv_m.z *= Lbar * Lbar;
   }
 
   const Field2D& get_Rxy() const { return Rxy_m; }
   const Field2D& get_Bpxy() const { return Bpxy_m; }
   const Field2D& get_Btxy() const { return Btxy_m; }
   const Field2D& get_Bxy() const { return Bxy_m; }
-  const Vector2D& get_b0xcv() const { return b0xcv_m; }
   const Field2D& get_hthe() const { return hthe_m; }
   const FieldMetric& get_ShearFactor() const { return ShearFactor_m; }
 };
diff --git a/tests/MMS/GBS/gbs.cxx b/tests/MMS/GBS/gbs.cxx
index d6795fbe37..f0a74c88cf 100644
--- a/tests/MMS/GBS/gbs.cxx
+++ b/tests/MMS/GBS/gbs.cxx
@@ -317,7 +317,7 @@ void GBS::LoadMetric(BoutReal Lnorm, BoutReal Bnorm) {
 
   auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
   tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-  coords = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+  coords = tokamak_coordinates_factory.make_tokamak_coordinates(true);
 
 // just define a macro for V_E dot Grad
 #define vE_Grad(f, p) (bracket(p, f, bm_exb))
diff --git a/tests/MMS/GBS/mms-slab2d.py b/tests/MMS/GBS/mms-slab2d.py
index 443517ab07..215e711a53 100644
--- a/tests/MMS/GBS/mms-slab2d.py
+++ b/tests/MMS/GBS/mms-slab2d.py
@@ -78,7 +78,15 @@ def C(f):
 
 bxcvz = 100  # Curvature
 
-tokamak_coordinates_factory.normalise(rho_s0, Bnorm);
+# Normalise
+
+Rxy /= rho_s0
+Bpxy /= Bnorm
+Btxy /= Bnorm
+Bxy /= Bnorm
+hthe /= rho_s0
+
+dx /= rho_s0**2 * Bnorm
 
 bxcvz *= rho_s0**2
 
diff --git a/tests/MMS/GBS/mms-slab3d.py b/tests/MMS/GBS/mms-slab3d.py
index cafff0ba76..2f958cfa69 100644
--- a/tests/MMS/GBS/mms-slab3d.py
+++ b/tests/MMS/GBS/mms-slab3d.py
@@ -79,7 +79,15 @@ def C(f):
 
 estatic = True
 
-tokamak_coordinates_factory.normalise(rho_s0, Bnorm);
+# Normalise
+
+Rxy /= rho_s0
+Bpxy /= Bnorm
+Btxy /= Bnorm
+Bxy /= Bnorm
+hthe /= rho_s0
+
+dx /= rho_s0**2 * Bnorm
 
 bxcvz *= rho_s0**2
 
diff --git a/tests/MMS/elm-pb/elm_pb.cxx b/tests/MMS/elm-pb/elm_pb.cxx
index e8d946e3b6..cccd9ec2c1 100644
--- a/tests/MMS/elm-pb/elm_pb.cxx
+++ b/tests/MMS/elm-pb/elm_pb.cxx
@@ -29,6 +29,7 @@ class ELMpb : public PhysicsModel {
 private:
   // 2D inital profiles
   Field2D J0, P0;         // Current and pressure
+  Vector2D b0xcv;         // Curvature term
   Field2D beta, gradparB; // Used for Vpar terms
   Field2D phi0;           // When diamagnetic terms used
 
@@ -274,6 +275,10 @@ class ELMpb : public PhysicsModel {
 
     OPTION(globalOptions->getSection("solver"), mms, false);
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
@@ -283,6 +288,19 @@ class ELMpb : public PhysicsModel {
 
     phi_solver = Laplacian::create();
 
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    bool ShiftXderivs;
+    globalOptions->get("shiftXderivs", ShiftXderivs, false); // Read global flag
+    if (ShiftXderivs) {
+      if (not mesh->IncIntShear) {
+        // Dimits style, using local coordinate system
+        if (include_curvature) {
+          b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        }
+      }
+    }
 
     //////////////////////////////////////////////////////////////
     // NORMALISE QUANTITIES
@@ -339,6 +357,10 @@ class ELMpb : public PhysicsModel {
     J0 = -MU0 * Lbar * J0 / B0;
     P0 = 2.0 * MU0 * P0 / (Bbar * Bbar);
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
     BoutReal pnorm = max(P0, true); // Maximum over all processors
 
     vacuum_pressure *= pnorm; // Get pressure from fraction
@@ -369,7 +391,7 @@ class ELMpb : public PhysicsModel {
     }
 
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     //////////////////////////////////////////////////////////////
     // SHIFTED RADIAL COORDINATES
@@ -562,7 +584,7 @@ class ELMpb : public PhysicsModel {
 
     ddt(U) = SQ(B0) * b0xGrad_dot_Grad(Psi, J0, CELL_CENTRE); // Grad j term
 
-    ddt(U) += tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+    ddt(U) += b0xcv * Grad(P); // curvature term
 
     // Parallel current term
     ddt(U) -= SQ(B0) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
@@ -612,7 +634,7 @@ class ELMpb : public PhysicsModel {
       ddt(P) -= beta * Div_par_CtoL(Vpar);
 
       if (phi_curv) {
-        ddt(P) -= 2. * beta * tokamak_coordinates_factory.get_b0xcv() * Grad(phi);
+        ddt(P) -= 2. * beta * b0xcv * Grad(phi);
       }
 
       // Vpar equation
diff --git a/tests/MMS/tokamak/tokamak.cxx b/tests/MMS/tokamak/tokamak.cxx
index da3b75acb6..bf066179fb 100644
--- a/tests/MMS/tokamak/tokamak.cxx
+++ b/tests/MMS/tokamak/tokamak.cxx
@@ -47,7 +47,7 @@ class TokamakMMS : public PhysicsModel {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
     tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-    auto coords = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    auto coords = tokamak_coordinates_factory.make_tokamak_coordinates(true);
   }
 
 private:
diff --git a/tests/integrated/test-drift-instability/2fluid.cxx b/tests/integrated/test-drift-instability/2fluid.cxx
index dbe0dfca35..f2fa3ebb91 100644
--- a/tests/integrated/test-drift-instability/2fluid.cxx
+++ b/tests/integrated/test-drift-instability/2fluid.cxx
@@ -23,6 +23,7 @@ class TwoFluid : public PhysicsModel {
   Field2D Ni0, Ti0, Te0, Vi0, phi0, Ve0, rho0, Ajpar0;
   // Staggered versions of initial profiles
   Field2D Ni0_maybe_ylow, Te0_maybe_ylow;
+  Vector2D b0xcv; // for curvature terms
 
   // 3D evolving fields
   Field3D rho, Te, Ni, Ajpar, Vi, Ti;
@@ -81,6 +82,10 @@ class TwoFluid : public PhysicsModel {
     GRID_LOAD(rho0);
     GRID_LOAD(Ajpar0);
 
+    // Load magnetic curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // b0xkappa terms
+
     // Load normalisation values
     GRID_LOAD(Te_x);
     GRID_LOAD(Ti_x);
@@ -129,6 +134,7 @@ class TwoFluid : public PhysicsModel {
     const bool ShiftXderivs = (*globalOptions)["ShiftXderivs"].withDefault(false);
     if (ShiftXderivs) {
       ShearFactor = 0.0; // I disappears from metric
+      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
       noshear = true;
     }
 
@@ -178,12 +184,17 @@ class TwoFluid : public PhysicsModel {
     phi0 /= Te_x;
     Vi0 /= Vi_x;
 
+    // Normalise curvature term
+    b0xcv.x /= (bmag / 1e4);
+    b0xcv.y *= rho_s * rho_s;
+    b0xcv.z *= rho_s * rho_s;
+
     // calculate pressures
     pei0 = (Ti0 + Te0) * Ni0;
     pe0 = Te0 * Ni0;
 
     tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     /**************** SET EVOLVING VARIABLES *************/
 
@@ -338,8 +349,8 @@ class TwoFluid : public PhysicsModel {
     if (evolve_te) {
       ddt(Te) -= vE_Grad(Te0, phi) + vE_Grad(Te, phi0) + vE_Grad(Te, phi);
       ddt(Te) -= Vpar_Grad_par(Ve, Te0) + Vpar_Grad_par(Ve0, Te) + Vpar_Grad_par(Ve, Te);
-      ddt(Te) += 1.333 * Te0 * (V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), pe) / Ni0 - V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), phi));
-      ddt(Te) += 3.333 * Te0 * V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), Te);
+      ddt(Te) += 1.333 * Te0 * (V_dot_Grad(b0xcv, pe) / Ni0 - V_dot_Grad(b0xcv, phi));
+      ddt(Te) += 3.333 * Te0 * V_dot_Grad(b0xcv, Te);
       ddt(Te) += (0.6666667 / Ni0) * Div_par_K_Grad_par(kapa_Te, Te);
     }
 
@@ -350,8 +361,8 @@ class TwoFluid : public PhysicsModel {
       ddt(Ti) -= vE_Grad(Ti0, phi) + vE_Grad(Ti, phi0) + vE_Grad(Ti, phi);
       ddt(Ti) -= Vpar_Grad_par(Vi, Ti0) + Vpar_Grad_par(Vi0, Ti) + Vpar_Grad_par(Vi, Ti);
       ddt(Ti) +=
-          1.333 * (Ti0 * V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), pe) / Ni0 - Ti * V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), phi));
-      ddt(Ti) -= 3.333 * Ti0 * V_dot_Grad(tokamak_coordinates_factory.get_b0xcv(), Ti);
+          1.333 * (Ti0 * V_dot_Grad(b0xcv, pe) / Ni0 - Ti * V_dot_Grad(b0xcv, phi));
+      ddt(Ti) -= 3.333 * Ti0 * V_dot_Grad(b0xcv, Ti);
       ddt(Ti) += (0.6666667 / Ni0) * Div_par_K_Grad_par(kapa_Ti, Ti);
     }
 
diff --git a/tests/integrated/test-interchange-instability/2fluid.cxx b/tests/integrated/test-interchange-instability/2fluid.cxx
index 815bd85101..8df63b3301 100644
--- a/tests/integrated/test-interchange-instability/2fluid.cxx
+++ b/tests/integrated/test-interchange-instability/2fluid.cxx
@@ -18,6 +18,7 @@ class Interchange : public PhysicsModel {
 
   // 2D initial profiles
   Field2D Ni0, Ti0, Te0;
+  Vector2D b0xcv; // for curvature terms
 
   // 3D evolving fields
   Field3D rho, Ni;
@@ -47,6 +48,12 @@ class Interchange : public PhysicsModel {
     GRID_LOAD(Ti0);
     GRID_LOAD(Te0);
 
+    // Load magnetic curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // b0xkappa terms
+
+    b0xcv *= -1.0; // NOTE: THIS IS FOR 'OLD' GRID FILES ONLY
+
     // Load normalisation values
     GRID_LOAD(Te_x);
     GRID_LOAD(Ti_x);
@@ -93,6 +100,13 @@ class Interchange : public PhysicsModel {
           hthe0 / rho_s);
     }
 
+    tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
+
+    if (ShiftXderivs) {
+      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+    }
+
     /************** NORMALISE QUANTITIES *****************/
 
     output.write("\tNormalising to rho_s = {:e}\n", rho_s);
@@ -102,9 +116,10 @@ class Interchange : public PhysicsModel {
     Ti0 /= Te_x;
     Te0 /= Te_x;
 
-//    b0xcv *= -1.0; // NOTE: THIS IS FOR 'OLD' GRID FILES ONLY  // TODO: Check if needed
-    tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    // Normalise curvature term
+    b0xcv.x /= (bmag / 1e4);
+    b0xcv.y *= rho_s * rho_s;
+    b0xcv.z *= rho_s * rho_s;
 
     // Tell BOUT++ which variables to evolve
     SOLVE_FOR2(rho, Ni);
@@ -132,7 +147,7 @@ class Interchange : public PhysicsModel {
     ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy();
 
     // VORTICITY
-    ddt(rho) = 2.0 * coord->Bxy() * tokamak_coordinates_factory.get_b0xcv() * Grad(pei);
+    ddt(rho) = 2.0 * coord->Bxy() * b0xcv * Grad(pei);
 
     return (0);
   }
diff --git a/tests/integrated/test-laplacexy/loadmetric.cxx b/tests/integrated/test-laplacexy/loadmetric.cxx
index c3151abbec..ef76de8018 100644
--- a/tests/integrated/test-laplacexy/loadmetric.cxx
+++ b/tests/integrated/test-laplacexy/loadmetric.cxx
@@ -19,5 +19,5 @@ void LoadMetric(BoutReal Lnorm, BoutReal Bnorm) {
 
   auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
   tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-  auto coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+  auto coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 }