diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx
index 0bbfad8cfc..971d866c16 100644
--- a/include/bout/index_derivs_interface.hxx
+++ b/include/bout/index_derivs_interface.hxx
@@ -96,6 +96,9 @@ T flowDerivative(const T& vel, const T& f, CELL_LOC outloc, const std::string& m
   T result(localmesh);
   result.allocate(); // Make sure data allocated
   result.setLocation(outloc);
+  if (std::is_base_of<Field3D, T>::value) {
+    result.setCoordinateSystem(f.getCoordinateSystem());
+  }
 
   // Apply method
   derivativeMethod(vel, f, result, region);
@@ -162,6 +165,9 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method,
   T result(localmesh);
   result.allocate(); // Make sure data allocated
   result.setLocation(outloc);
+  if (std::is_base_of<Field3D, T>::value) {
+    result.setCoordinateSystem(f.getCoordinateSystem());
+  }
 
   // Apply method
   derivativeMethod(f, result, region);
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 0508de4cc3..3fe31a840c 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -689,7 +689,7 @@ class Mesh {
 
   /// Transform a field into field-aligned coordinates
   const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
-    return getParallelTransform().toFieldAligned(f, region);
+    return transform->toFieldAligned(f, region);
   }
   const Field2D toFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) {
     return f;
@@ -697,14 +697,24 @@ class Mesh {
   
   /// Convert back into standard form
   const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
-    return getParallelTransform().fromFieldAligned(f, region);
+    return transform->fromFieldAligned(f, region);
   }
   const Field2D fromFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) {
     return f;
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const {
+    if (transform) {
+      return transform->getCoordinateSystem();
+    } else {
+      // ParallelTransform not initialized yet. Any Field3D created better not
+      // need to know its coordinate system.
+      return COORDINATE_SYSTEM::None;
+    }
+  }
+
   bool canToFromFieldAligned() {
-    return getParallelTransform().canToFromFieldAligned();
+    return transform->canToFromFieldAligned();
   }
 
   /*!
@@ -725,7 +735,7 @@ class Mesh {
   void setParallelTransform();
 
   /*!
-   * Return the parallel transform, setting it if need be
+   * Return the parallel transform
    */
   ParallelTransform& getParallelTransform();
   
diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index f313c978c3..fd7ca135ef 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -43,6 +43,11 @@ public:
   /// into standard form
   virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
+  /*!
+   * Return the default coordinate system for Field3Ds when using this ParallelTransform
+   */
+  virtual COORDINATE_SYSTEM getCoordinateSystem() const = 0;
+
   virtual bool canToFromFieldAligned() = 0;
 };
 
@@ -76,6 +81,10 @@ public:
     return f;
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const override {
+    return COORDINATE_SYSTEM::FieldAligned;
+  }
+
   bool canToFromFieldAligned() override{
     return true;
   }
@@ -116,6 +125,10 @@ public:
    */
   const Field3D fromFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
+  COORDINATE_SYSTEM getCoordinateSystem() const override {
+    return COORDINATE_SYSTEM::Orthogonal;
+  }
+
   bool canToFromFieldAligned() override{
     return true;
   }
diff --git a/include/field.hxx b/include/field.hxx
index 13c288f635..9c17b055ac 100644
--- a/include/field.hxx
+++ b/include/field.hxx
@@ -46,6 +46,28 @@ extern Mesh * mesh; ///< Global mesh
 #include <string>
 #endif
 
+/*!
+ * Coordinate systems that Field3Ds can be stored in
+ *
+ * None - special value to be used for Field3Ds initialized before the Mesh and
+ * ParallelTransform have finished being initialized, because then we cannot
+ * get the COORDINATE_SYSTEM to use from the ParallelTransform, but Field3Ds
+ * created then must not need to know their COORDINATE_SYSTEM, or must set it
+ * explicitly later.
+ */
+enum class COORDINATE_SYSTEM { None, FieldAligned, Orthogonal, FCI };
+
+#define COORDENUMSTR(val) {COORDINATE_SYSTEM::val, #val}
+
+const std::map<COORDINATE_SYSTEM, std::string> COORDINATE_SYSTEMtoString = {
+  COORDENUMSTR(None),
+  COORDENUMSTR(FieldAligned),
+  COORDENUMSTR(Orthogonal),
+  COORDENUMSTR(FCI)
+};
+
+#undef COORDENUMSTR
+
 /*!
  * \brief Base class for fields
  *
diff --git a/include/field2d.hxx b/include/field2d.hxx
index adbf8b0800..d5174750c4 100644
--- a/include/field2d.hxx
+++ b/include/field2d.hxx
@@ -111,6 +111,24 @@ class Field2D : public Field, public FieldData {
    */
   int getNz() const override {return 1;};
 
+  /*!
+   * Return the coordinate system of this field
+   */
+  COORDINATE_SYSTEM getCoordinateSystem() const {
+#if CHECK > 0
+    throw BoutException("getCoordinateSystem() should not be called for Field2D");
+#endif
+  }
+
+  /*!
+   * Reset the coordinate system of this field
+   */
+  void setCoordinateSystem(COORDINATE_SYSTEM UNUSED(new_coords)) {
+#if CHECK > 0
+    throw BoutException("setCoordinateSystem() should not be called for Field2D");
+#endif
+  }
+
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
     return true;
diff --git a/include/field3d.hxx b/include/field3d.hxx
index f72d21f502..56269fece9 100644
--- a/include/field3d.hxx
+++ b/include/field3d.hxx
@@ -215,6 +215,16 @@ class Field3D : public Field, public FieldData {
    */
   int getNz() const override {return nz;};
 
+  /*!
+   * Return the coordinate system of this field
+   */
+  COORDINATE_SYSTEM getCoordinateSystem() const { return coordinate_system; }
+
+  /*!
+   * Reset the coordinate system of this field
+   */
+  void setCoordinateSystem(COORDINATE_SYSTEM new_coords) { coordinate_system = new_coords; }
+
   /*!
    * Ensure that this field has separate fields
    * for yup and ydown.
@@ -460,6 +470,7 @@ class Field3D : public Field, public FieldData {
     swap(first.nx, second.nx);
     swap(first.ny, second.ny);
     swap(first.nz, second.nz);
+    swap(first.coordinate_system, second.coordinate_system);
     swap(first.location, second.location);
     swap(first.deriv, second.deriv);
     swap(first.yup_field, second.yup_field);
@@ -478,6 +489,8 @@ private:
   /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null
   int nx{-1}, ny{-1}, nz{-1};
 
+  COORDINATE_SYSTEM coordinate_system; ///< Coordinate system used by this field, e.g. fieldaligned, orthogonal
+
   /// Internal data array. Handles allocation/freeing of memory
   Array<BoutReal> data;
 
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 6ab7a4ef79..428db8e274 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -54,11 +54,19 @@ Field3D::Field3D(Mesh* localmesh) : Field(localmesh) {
     ny = fieldmesh->LocalNy;
     nz = fieldmesh->LocalNz;
   }
+  if (fieldmesh) {
+    coordinate_system = fieldmesh->getCoordinateSystem();
+  } else {
+    coordinate_system = COORDINATE_SYSTEM::None;
+  }
+
 }
 
 /// Doesn't copy any data, just create a new reference to the same data (copy on change
 /// later)
-Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) {
+Field3D::Field3D(const Field3D& f)
+  : Field(f.fieldmesh), coordinate_system(f.coordinate_system), data(f.data),
+    location(f.location) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -72,7 +80,6 @@ Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) {
     nz = fieldmesh->LocalNz;
   }
 
-  location = f.location;
   fieldCoordinates = f.fieldCoordinates;
 }
 
@@ -86,6 +93,7 @@ Field3D::Field3D(const Field2D& f) : Field(f.getMesh()) {
 
   location = f.getLocation();
   fieldCoordinates = nullptr;
+  coordinate_system = fieldmesh->getCoordinateSystem();
 
   *this = f;
 }
@@ -98,6 +106,8 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field(localmesh) {
   ny = fieldmesh->LocalNy;
   nz = fieldmesh->LocalNz;
 
+  coordinate_system = fieldmesh->getCoordinateSystem();
+
   *this = val;
 }
 
@@ -275,6 +285,8 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
   fieldmesh = rhs.fieldmesh;
   nx = rhs.nx; ny = rhs.ny; nz = rhs.nz; 
   
+  coordinate_system = rhs.coordinate_system;
+
   data = rhs.data;
 
   setLocation(rhs.location);
@@ -614,6 +626,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) {
   TRACE("pow(Field3D, Field3D)");
 
   ASSERT1(lhs.getLocation() == rhs.getLocation());
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
 
   ASSERT1(lhs.getMesh() == rhs.getMesh());
   Field3D result(lhs.getMesh());
@@ -622,6 +635,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
   
   checkData(result);
   return result;
@@ -641,6 +655,7 @@ Field3D pow(const Field3D &lhs, const Field2D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
   
   checkData(result);
   return result;
@@ -653,6 +668,11 @@ FieldPerp pow(const Field3D &lhs, const FieldPerp &rhs, REGION rgn) {
   checkData(rhs);
   ASSERT1(lhs.getMesh() == rhs.getMesh());
 
+  // Assume FieldPerp is always on the default coordinate system of its Mesh,
+  // since at the moment FieldPerp cannot switch between orthogonal and field
+  // aligned coordinates
+  ASSERT2(lhs.getCoordinateSystem() == rhs.getMesh()->getCoordinateSystem());
+
   FieldPerp result{rhs.getMesh()};
   result.allocate();
   result.setIndex(rhs.getIndex());
@@ -679,6 +699,7 @@ Field3D pow(const Field3D &lhs, BoutReal rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -697,6 +718,7 @@ Field3D pow(BoutReal lhs, const Field3D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); }
 
   result.setLocation( rhs.getLocation() );
+  result.setCoordinateSystem( rhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -802,6 +824,7 @@ BoutReal mean(const Field3D &f, bool allpe, REGION rgn) {
     result.allocate();                                                                   \
     BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); }                       \
     result.setLocation(f.getLocation());                                                 \
+    result.setCoordinateSystem(f.getCoordinateSystem());                                 \
     checkData(result);                                                                   \
     return result;                                                                       \
   }
@@ -865,6 +888,7 @@ const Field3D filter(const Field3D &var, int N0, REGION rgn) {
 #endif
 
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;
@@ -912,6 +936,7 @@ const Field3D lowPass(const Field3D &var, int zmax, REGION rgn) {
     }
   }
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;
@@ -962,6 +987,7 @@ const Field3D lowPass(const Field3D &var, int zmax, int zmin, REGION rgn) {
   }
 
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;
diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx
index bc679e8a30..3daf96054d 100644
--- a/src/field/field_factory.cxx
+++ b/src/field/field_factory.cxx
@@ -229,6 +229,7 @@ const Field3D FieldFactory::create3D(const std::string &value, const Options *op
   if (localmesh->canToFromFieldAligned()){ // Ask wheter it is possible
     // Transform from field aligned coordinates, to be compatible with
     // older BOUT++ inputs. This is not a particularly "nice" solution.
+    result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
     result = localmesh->fromFieldAligned(result, RGN_ALL);
   }
 
diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja
index 58f07d0796..82dc3cc95e 100644
--- a/src/field/gen_fieldops.jinja
+++ b/src/field/gen_fieldops.jinja
@@ -16,6 +16,10 @@
     ASSERT1(localmesh == {{rhs.name}}.getMesh());
   {% endif %}
 
+  {% if lhs == "Field3D" and rhs == "Field3D" %}
+    ASSERT1({{lhs.name}}.getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem());
+  {% endif %}
+
   {{out.field_type}} {{out.name}}(localmesh);
   {{out.name}}.allocate();
   checkData({{lhs.name}});
@@ -46,8 +50,12 @@
 
   {% if out in ['Field3D','Field2D'] %}
     {{out.name}}.setLocation({{rhs.name if rhs in ["Field3D","Field2D"] else lhs.name}}.getLocation());
+
   {% endif %}
+  {% if out in ['Field3D'] %}
+    {{out.name}}.setCoordinateSystem({{rhs.name if rhs == "Field3D" else lhs.name}}.getCoordinateSystem());
 
+  {% endif %}
   checkData({{out.name}});
   return {{out.name}};
 }
@@ -58,7 +66,11 @@
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-    {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %}    
+    {% if ((lhs in ["Field3D"]) and (rhs in ["Field3D"])) %}
+    ASSERT1(this->getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem())
+
+    {% endif %}
+    {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %}
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException(
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 0aa552060f..f11e207659 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -7,7 +7,7 @@
 #include <interpolation.hxx>
 
 // Provide the C++ wrapper for multiplication of Field3D and Field3D
-Field3D operator*(const Field3D &lhs, const Field3D &rhs) {
+Field3D operator*(const Field3D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator*(Field3D, Field3D): fields at different "
@@ -16,10 +16,12 @@ Field3D operator*(const Field3D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
+
   Field3D result(localmesh);
   result.allocate();
   checkData(lhs);
@@ -31,15 +33,18 @@ Field3D operator*(const Field3D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by multiplication with Field3D
-Field3D &Field3D::operator*=(const Field3D &rhs) {
+Field3D& Field3D::operator*=(const Field3D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
+    ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem())
 
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
@@ -66,7 +71,7 @@ Field3D &Field3D::operator*=(const Field3D &rhs) {
 }
 
 // Provide the C++ wrapper for division of Field3D and Field3D
-Field3D operator/(const Field3D &lhs, const Field3D &rhs) {
+Field3D operator/(const Field3D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator/(Field3D, Field3D): fields at different "
@@ -75,10 +80,12 @@ Field3D operator/(const Field3D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
+
   Field3D result(localmesh);
   result.allocate();
   checkData(lhs);
@@ -90,15 +97,18 @@ Field3D operator/(const Field3D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by division with Field3D
-Field3D &Field3D::operator/=(const Field3D &rhs) {
+Field3D& Field3D::operator/=(const Field3D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
+    ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem())
 
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
@@ -125,7 +135,7 @@ Field3D &Field3D::operator/=(const Field3D &rhs) {
 }
 
 // Provide the C++ wrapper for addition of Field3D and Field3D
-Field3D operator+(const Field3D &lhs, const Field3D &rhs) {
+Field3D operator+(const Field3D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator+(Field3D, Field3D): fields at different "
@@ -134,10 +144,12 @@ Field3D operator+(const Field3D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
+
   Field3D result(localmesh);
   result.allocate();
   checkData(lhs);
@@ -149,15 +161,18 @@ Field3D operator+(const Field3D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by addition with Field3D
-Field3D &Field3D::operator+=(const Field3D &rhs) {
+Field3D& Field3D::operator+=(const Field3D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
+    ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem())
 
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
@@ -184,7 +199,7 @@ Field3D &Field3D::operator+=(const Field3D &rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of Field3D and Field3D
-Field3D operator-(const Field3D &lhs, const Field3D &rhs) {
+Field3D operator-(const Field3D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator-(Field3D, Field3D): fields at different "
@@ -193,10 +208,12 @@ Field3D operator-(const Field3D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
+
   Field3D result(localmesh);
   result.allocate();
   checkData(lhs);
@@ -208,15 +225,18 @@ Field3D operator-(const Field3D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by subtraction with Field3D
-Field3D &Field3D::operator-=(const Field3D &rhs) {
+Field3D& Field3D::operator-=(const Field3D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
+    ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem())
 
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
@@ -243,7 +263,7 @@ Field3D &Field3D::operator-=(const Field3D &rhs) {
 }
 
 // Provide the C++ wrapper for multiplication of Field3D and Field2D
-Field3D operator*(const Field3D &lhs, const Field2D &rhs) {
+Field3D operator*(const Field3D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator*(Field3D, Field2D): fields at different "
@@ -252,7 +272,7 @@ Field3D operator*(const Field3D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -270,16 +290,17 @@ Field3D operator*(const Field3D &lhs, const Field2D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by multiplication with Field2D
-Field3D &Field3D::operator*=(const Field2D &rhs) {
+Field3D& Field3D::operator*=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field3D::operator*=(Field2D): fields at different "
@@ -310,7 +331,7 @@ Field3D &Field3D::operator*=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for division of Field3D and Field2D
-Field3D operator/(const Field3D &lhs, const Field2D &rhs) {
+Field3D operator/(const Field3D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator/(Field3D, Field2D): fields at different "
@@ -319,7 +340,7 @@ Field3D operator/(const Field3D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -338,16 +359,17 @@ Field3D operator/(const Field3D &lhs, const Field2D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by division with Field2D
-Field3D &Field3D::operator/=(const Field2D &rhs) {
+Field3D& Field3D::operator/=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field3D::operator/=(Field2D): fields at different "
@@ -379,7 +401,7 @@ Field3D &Field3D::operator/=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for addition of Field3D and Field2D
-Field3D operator+(const Field3D &lhs, const Field2D &rhs) {
+Field3D operator+(const Field3D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator+(Field3D, Field2D): fields at different "
@@ -388,7 +410,7 @@ Field3D operator+(const Field3D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -406,16 +428,17 @@ Field3D operator+(const Field3D &lhs, const Field2D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by addition with Field2D
-Field3D &Field3D::operator+=(const Field2D &rhs) {
+Field3D& Field3D::operator+=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field3D::operator+=(Field2D): fields at different "
@@ -446,7 +469,7 @@ Field3D &Field3D::operator+=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of Field3D and Field2D
-Field3D operator-(const Field3D &lhs, const Field2D &rhs) {
+Field3D operator-(const Field3D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator-(Field3D, Field2D): fields at different "
@@ -455,7 +478,7 @@ Field3D operator-(const Field3D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -473,16 +496,17 @@ Field3D operator-(const Field3D &lhs, const Field2D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by subtraction with Field2D
-Field3D &Field3D::operator-=(const Field2D &rhs) {
+Field3D& Field3D::operator-=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field3D::operator-=(Field2D): fields at different "
@@ -513,9 +537,9 @@ Field3D &Field3D::operator-=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for multiplication of Field3D and BoutReal
-Field3D operator*(const Field3D &lhs, const BoutReal rhs) {
+Field3D operator*(const Field3D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -526,12 +550,14 @@ Field3D operator*(const Field3D &lhs, const BoutReal rhs) {
 
   result.setLocation(lhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by multiplication with BoutReal
-Field3D &Field3D::operator*=(const BoutReal rhs) {
+Field3D& Field3D::operator*=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -550,9 +576,9 @@ Field3D &Field3D::operator*=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for division of Field3D and BoutReal
-Field3D operator/(const Field3D &lhs, const BoutReal rhs) {
+Field3D operator/(const Field3D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -563,12 +589,14 @@ Field3D operator/(const Field3D &lhs, const BoutReal rhs) {
 
   result.setLocation(lhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by division with BoutReal
-Field3D &Field3D::operator/=(const BoutReal rhs) {
+Field3D& Field3D::operator/=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -587,9 +615,9 @@ Field3D &Field3D::operator/=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for addition of Field3D and BoutReal
-Field3D operator+(const Field3D &lhs, const BoutReal rhs) {
+Field3D operator+(const Field3D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -600,12 +628,14 @@ Field3D operator+(const Field3D &lhs, const BoutReal rhs) {
 
   result.setLocation(lhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by addition with BoutReal
-Field3D &Field3D::operator+=(const BoutReal rhs) {
+Field3D& Field3D::operator+=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -624,9 +654,9 @@ Field3D &Field3D::operator+=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of Field3D and BoutReal
-Field3D operator-(const Field3D &lhs, const BoutReal rhs) {
+Field3D operator-(const Field3D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -637,12 +667,14 @@ Field3D operator-(const Field3D &lhs, const BoutReal rhs) {
 
   result.setLocation(lhs.getLocation());
 
+  result.setCoordinateSystem(lhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ operator to update Field3D by subtraction with BoutReal
-Field3D &Field3D::operator-=(const BoutReal rhs) {
+Field3D& Field3D::operator-=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -661,7 +693,7 @@ Field3D &Field3D::operator-=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for multiplication of Field2D and Field3D
-Field3D operator*(const Field2D &lhs, const Field3D &rhs) {
+Field3D operator*(const Field2D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator*(Field2D, Field3D): fields at different "
@@ -670,7 +702,7 @@ Field3D operator*(const Field2D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -688,12 +720,14 @@ Field3D operator*(const Field2D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for division of Field2D and Field3D
-Field3D operator/(const Field2D &lhs, const Field3D &rhs) {
+Field3D operator/(const Field2D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator/(Field2D, Field3D): fields at different "
@@ -702,7 +736,7 @@ Field3D operator/(const Field2D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -720,12 +754,14 @@ Field3D operator/(const Field2D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for addition of Field2D and Field3D
-Field3D operator+(const Field2D &lhs, const Field3D &rhs) {
+Field3D operator+(const Field2D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator+(Field2D, Field3D): fields at different "
@@ -734,7 +770,7 @@ Field3D operator+(const Field2D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -752,12 +788,14 @@ Field3D operator+(const Field2D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for subtraction of Field2D and Field3D
-Field3D operator-(const Field2D &lhs, const Field3D &rhs) {
+Field3D operator-(const Field2D& lhs, const Field3D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator-(Field2D, Field3D): fields at different "
@@ -766,7 +804,7 @@ Field3D operator-(const Field2D &lhs, const Field3D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -784,12 +822,14 @@ Field3D operator-(const Field2D &lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for multiplication of Field2D and Field2D
-Field2D operator*(const Field2D &lhs, const Field2D &rhs) {
+Field2D operator*(const Field2D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator*(Field2D, Field2D): fields at different "
@@ -798,7 +838,7 @@ Field2D operator*(const Field2D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -818,11 +858,10 @@ Field2D operator*(const Field2D &lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ operator to update Field2D by multiplication with Field2D
-Field2D &Field2D::operator*=(const Field2D &rhs) {
+Field2D& Field2D::operator*=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field2D::operator*=(Field2D): fields at different "
@@ -848,7 +887,7 @@ Field2D &Field2D::operator*=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for division of Field2D and Field2D
-Field2D operator/(const Field2D &lhs, const Field2D &rhs) {
+Field2D operator/(const Field2D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator/(Field2D, Field2D): fields at different "
@@ -857,7 +896,7 @@ Field2D operator/(const Field2D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -877,11 +916,10 @@ Field2D operator/(const Field2D &lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ operator to update Field2D by division with Field2D
-Field2D &Field2D::operator/=(const Field2D &rhs) {
+Field2D& Field2D::operator/=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field2D::operator/=(Field2D): fields at different "
@@ -907,7 +945,7 @@ Field2D &Field2D::operator/=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for addition of Field2D and Field2D
-Field2D operator+(const Field2D &lhs, const Field2D &rhs) {
+Field2D operator+(const Field2D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator+(Field2D, Field2D): fields at different "
@@ -916,7 +954,7 @@ Field2D operator+(const Field2D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -936,11 +974,10 @@ Field2D operator+(const Field2D &lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ operator to update Field2D by addition with Field2D
-Field2D &Field2D::operator+=(const Field2D &rhs) {
+Field2D& Field2D::operator+=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field2D::operator+=(Field2D): fields at different "
@@ -966,7 +1003,7 @@ Field2D &Field2D::operator+=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of Field2D and Field2D
-Field2D operator-(const Field2D &lhs, const Field2D &rhs) {
+Field2D operator-(const Field2D& lhs, const Field2D& rhs) {
 #if CHECK > 0
   if (lhs.getLocation() != rhs.getLocation()) {
     throw BoutException("Error in operator-(Field2D, Field2D): fields at different "
@@ -975,7 +1012,7 @@ Field2D operator-(const Field2D &lhs, const Field2D &rhs) {
   }
 #endif
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   ASSERT1(localmesh == rhs.getMesh());
 
@@ -995,11 +1032,10 @@ Field2D operator-(const Field2D &lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ operator to update Field2D by subtraction with Field2D
-Field2D &Field2D::operator-=(const Field2D &rhs) {
+Field2D& Field2D::operator-=(const Field2D& rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException("Error in Field2D::operator-=(Field2D): fields at different "
@@ -1025,9 +1061,9 @@ Field2D &Field2D::operator-=(const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for multiplication of Field2D and BoutReal
-Field2D operator*(const Field2D &lhs, const BoutReal rhs) {
+Field2D operator*(const Field2D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1043,7 +1079,7 @@ Field2D operator*(const Field2D &lhs, const BoutReal rhs) {
 }
 
 // Provide the C++ operator to update Field2D by multiplication with BoutReal
-Field2D &Field2D::operator*=(const BoutReal rhs) {
+Field2D& Field2D::operator*=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -1062,9 +1098,9 @@ Field2D &Field2D::operator*=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for division of Field2D and BoutReal
-Field2D operator/(const Field2D &lhs, const BoutReal rhs) {
+Field2D operator/(const Field2D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1080,7 +1116,7 @@ Field2D operator/(const Field2D &lhs, const BoutReal rhs) {
 }
 
 // Provide the C++ operator to update Field2D by division with BoutReal
-Field2D &Field2D::operator/=(const BoutReal rhs) {
+Field2D& Field2D::operator/=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -1099,9 +1135,9 @@ Field2D &Field2D::operator/=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for addition of Field2D and BoutReal
-Field2D operator+(const Field2D &lhs, const BoutReal rhs) {
+Field2D operator+(const Field2D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1117,7 +1153,7 @@ Field2D operator+(const Field2D &lhs, const BoutReal rhs) {
 }
 
 // Provide the C++ operator to update Field2D by addition with BoutReal
-Field2D &Field2D::operator+=(const BoutReal rhs) {
+Field2D& Field2D::operator+=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -1136,9 +1172,9 @@ Field2D &Field2D::operator+=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of Field2D and BoutReal
-Field2D operator-(const Field2D &lhs, const BoutReal rhs) {
+Field2D operator-(const Field2D& lhs, const BoutReal rhs) {
 
-  Mesh *localmesh = lhs.getMesh();
+  Mesh* localmesh = lhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1154,7 +1190,7 @@ Field2D operator-(const Field2D &lhs, const BoutReal rhs) {
 }
 
 // Provide the C++ operator to update Field2D by subtraction with BoutReal
-Field2D &Field2D::operator-=(const BoutReal rhs) {
+Field2D& Field2D::operator-=(const BoutReal rhs) {
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
@@ -1173,9 +1209,9 @@ Field2D &Field2D::operator-=(const BoutReal rhs) {
 }
 
 // Provide the C++ wrapper for multiplication of BoutReal and Field3D
-Field3D operator*(const BoutReal lhs, const Field3D &rhs) {
+Field3D operator*(const BoutReal lhs, const Field3D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -1186,14 +1222,16 @@ Field3D operator*(const BoutReal lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for division of BoutReal and Field3D
-Field3D operator/(const BoutReal lhs, const Field3D &rhs) {
+Field3D operator/(const BoutReal lhs, const Field3D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -1204,14 +1242,16 @@ Field3D operator/(const BoutReal lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for addition of BoutReal and Field3D
-Field3D operator+(const BoutReal lhs, const Field3D &rhs) {
+Field3D operator+(const BoutReal lhs, const Field3D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -1222,14 +1262,16 @@ Field3D operator+(const BoutReal lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for subtraction of BoutReal and Field3D
-Field3D operator-(const BoutReal lhs, const Field3D &rhs) {
+Field3D operator-(const BoutReal lhs, const Field3D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field3D result(localmesh);
   result.allocate();
@@ -1240,14 +1282,16 @@ Field3D operator-(const BoutReal lhs, const Field3D &rhs) {
 
   result.setLocation(rhs.getLocation());
 
+  result.setCoordinateSystem(rhs.getCoordinateSystem());
+
   checkData(result);
   return result;
 }
 
 // Provide the C++ wrapper for multiplication of BoutReal and Field2D
-Field2D operator*(const BoutReal lhs, const Field2D &rhs) {
+Field2D operator*(const BoutReal lhs, const Field2D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1263,9 +1307,9 @@ Field2D operator*(const BoutReal lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for division of BoutReal and Field2D
-Field2D operator/(const BoutReal lhs, const Field2D &rhs) {
+Field2D operator/(const BoutReal lhs, const Field2D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1281,9 +1325,9 @@ Field2D operator/(const BoutReal lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for addition of BoutReal and Field2D
-Field2D operator+(const BoutReal lhs, const Field2D &rhs) {
+Field2D operator+(const BoutReal lhs, const Field2D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
@@ -1299,9 +1343,9 @@ Field2D operator+(const BoutReal lhs, const Field2D &rhs) {
 }
 
 // Provide the C++ wrapper for subtraction of BoutReal and Field2D
-Field2D operator-(const BoutReal lhs, const Field2D &rhs) {
+Field2D operator-(const BoutReal lhs, const Field2D& rhs) {
 
-  Mesh *localmesh = rhs.getMesh();
+  Mesh* localmesh = rhs.getMesh();
 
   Field2D result(localmesh);
   result.allocate();
diff --git a/src/fileio/datafile.cxx b/src/fileio/datafile.cxx
index b6ef8a1298..41dcb1332e 100644
--- a/src/fileio/datafile.cxx
+++ b/src/fileio/datafile.cxx
@@ -930,6 +930,12 @@ bool Datafile::write() {
   // Write 3D fields
   for (const auto& var : f3d_arr) {
     write_f3d(var.name, var.ptr, var.save_repeat);
+    // Add coordinate system
+    if (shiftOutput) {
+      file->setAttribute(var.name, "coordinate_system", COORDINATE_SYSTEMtoString.at(COORDINATE_SYSTEM::FieldAligned));
+    } else {
+      file->setAttribute(var.name, "coordinate_system", COORDINATE_SYSTEMtoString.at(var.ptr->getCoordinateSystem()));
+    }
   }
   
   // 2D vectors
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 0bf59de09a..4548a9368f 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -839,7 +839,10 @@ int BoutMesh::load() {
   // Add boundary regions
   addBoundaryRegions();
 
-  output_info.write(_("\tdone\n"));
+  // Set the ParallelTransform
+  setParallelTransform();
+
+  output_info.write("\tdone\n");
 
   return 0;
 }
diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index ea06fa66e7..71090eca61 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -45,7 +45,6 @@ STAGGER Mesh::getStagger(const CELL_LOC inloc, const CELL_LOC outloc,
   TRACE("Mesh::getStagger -- three arguments");
   ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == allowedStaggerLoc)
           || (outloc == allowedStaggerLoc && inloc == CELL_CENTRE));
-
   if ((!StaggerGrids) || outloc == inloc)
     return STAGGER::None;
   if (outloc == allowedStaggerLoc) {
diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx
index 45b7d47d9b..b687030617 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -161,6 +161,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
             result = var_fa; // NOTE: This is just for boundaries. FIX!
             result.allocate();
           }
+          result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
           if (fieldmesh->ystart > 1) {
 
             // More than one guard cell, so set pp and mm values
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index ea4147fa3a..521b94a911 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -321,10 +321,7 @@ void Mesh::setParallelTransform() {
 }
 
 ParallelTransform& Mesh::getParallelTransform() {
-  if(!transform) {
-    // No ParallelTransform object yet. Set from options
-    setParallelTransform();
-  }
+  ASSERT1(transform);
   
   // Return a reference to the ParallelTransform object
   return *transform;
diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx
index 747f311c03..6e3e7c90a3 100644
--- a/src/mesh/parallel/fci.hxx
+++ b/src/mesh/parallel/fci.hxx
@@ -81,6 +81,10 @@ public:
     throw BoutException("FCI method cannot transform into field aligned grid");
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const override {
+    return COORDINATE_SYSTEM::FCI;
+  }
+
   bool canToFromFieldAligned() override{
     return false;
   }
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 7a7c439c14..231d63188e 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -106,7 +106,10 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
  * and Y is then field aligned.
  */
 const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) {
-  return shiftZ(f, toAlignedPhs, region);
+  ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal);
+  Field3D result = shiftZ(f, toAlignedPhs, region);
+  result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
+  return result;
 }
 
 /*!
@@ -114,7 +117,10 @@ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION regio
  * but Y is not field aligned.
  */
 const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) {
-  return shiftZ(f, fromAlignedPhs, region);
+  ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned);
+  Field3D result = shiftZ(f, fromAlignedPhs, region);
+  result.setCoordinateSystem(f.getMesh()->getCoordinateSystem());
+  return result;
 }
 
 const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor<dcomplex>& phs,
diff --git a/src/mesh/parallel/shifttofieldaligned.cxx b/src/mesh/parallel/shifttofieldaligned.cxx
index 9b760fd929..8bb06f02c8 100644
--- a/src/mesh/parallel/shifttofieldaligned.cxx
+++ b/src/mesh/parallel/shifttofieldaligned.cxx
@@ -120,7 +120,14 @@ void ShiftToFieldAligned::Implementation::cachePhases() {
  */
 const Field3D ShiftToFieldAligned::Implementation::toFieldAligned(const Field3D& f,
                                                                   const REGION region) {
-  return shiftZ(f, toAlignedPhs, region);
+  if (f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned) {
+    return f;
+  } else {
+    ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal);
+    Field3D result = shiftZ(f, toAlignedPhs, region);
+    result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
+    return result;
+  }
 }
 
 /*!
@@ -129,7 +136,10 @@ const Field3D ShiftToFieldAligned::Implementation::toFieldAligned(const Field3D&
  */
 const Field3D ShiftToFieldAligned::Implementation::fromFieldAligned(const Field3D& f,
                                                                     const REGION region) {
-  return shiftZ(f, fromAlignedPhs, region);
+  ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned);
+  Field3D result = shiftZ(f, fromAlignedPhs, region);
+  result.setCoordinateSystem(COORDINATE_SYSTEM::Orthogonal);
+  return result;
 }
 
 const Field3D ShiftToFieldAligned::Implementation::shiftZ(
diff --git a/src/mesh/parallel/shifttofieldaligned.hxx b/src/mesh/parallel/shifttofieldaligned.hxx
index 57b8fcfc6c..31b404819e 100644
--- a/src/mesh/parallel/shifttofieldaligned.hxx
+++ b/src/mesh/parallel/shifttofieldaligned.hxx
@@ -61,6 +61,8 @@ public:
     return implementations.at(f.getLocation()).fromFieldAligned(f, region);
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const { return COORDINATE_SYSTEM::Orthogonal; }
+
   bool canToFromFieldAligned() override { return true; }
 
 private:
diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx
index f7cef52403..9fa03bb49a 100644
--- a/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx
+++ b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx
@@ -9,6 +9,7 @@ const Field3D DDY_aligned(const Field3D &f) {
   Field3D result;
   result.allocate();
   result = 0.0;
+  result.setCoordinateSystem(f.getCoordinateSystem());
 
   for(int i=0;i<mesh->LocalNx;i++)
     for(int j=mesh->ystart;j<=mesh->yend;j++)
diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx
index e05f5db587..1c54d68dea 100644
--- a/tests/integrated/test-yupdown/test_yupdown.cxx
+++ b/tests/integrated/test-yupdown/test_yupdown.cxx
@@ -24,6 +24,7 @@ const Field3D DDY_aligned(const Field3D& f) {
   Field3D result;
   result.allocate();
   result = 0.0;
+  result.setCoordinateSystem(f.getCoordinateSystem());
 
   for (int i = 0; i < mesh->LocalNx; i++)
     for (int j = mesh->ystart; j <= mesh->yend; j++)