diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 056294bbe0..6c9dea8b11 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -609,13 +609,19 @@ class Mesh {
   /// Derivative functions of a velocity field, and field stencil v, f
   typedef BoutReal (*flux_func)(stencil&, stencil &);
 
+  /// Calculate yup/ydown fields in case they can be computed for an
+  /// intermediate variable without communicating
+  void calcYUpDown(Field3D &f) {
+    getParallelTransform().calcYUpDown(f);
+  }
+
   /// Transform a field into field-aligned coordinates
-  const Field3D toFieldAligned(const Field3D &f) {
-    return getParallelTransform().toFieldAligned(f);
+  const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
+    return getParallelTransform().toFieldAligned(f, region);
   }
   /// Convert back into standard form
-  const Field3D fromFieldAligned(const Field3D &f) {
-    return getParallelTransform().fromFieldAligned(f);
+  const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
+    return getParallelTransform().fromFieldAligned(f, region);
   }
 
   bool canToFromFieldAligned() {
diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index a25540354a..1546bf57ab 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -38,11 +38,11 @@ public:
   
   /// Convert a 3D field into field-aligned coordinates
   /// so that the y index is along the magnetic field
-  virtual const Field3D toFieldAligned(const Field3D &f) = 0;
+  virtual const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
   
   /// Convert back from field-aligned coordinates
   /// into standard form
-  virtual const Field3D fromFieldAligned(const Field3D &f) = 0;
+  virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
   virtual bool canToFromFieldAligned() = 0;
 };
@@ -59,13 +59,16 @@ public:
    * Merges the yup and ydown() fields of f, so that
    * f.yup() = f.ydown() = f
    */ 
-  void calcYUpDown(Field3D &f) override {f.mergeYupYdown();}
+  void calcYUpDown(Field3D &f) override {
+    f.mergeYupYdown();
+    f.setHasValidYUpDown(true);
+  }
   
   /*!
    * The field is already aligned in Y, so this
    * does nothing
    */ 
-  const Field3D toFieldAligned(const Field3D &f) override {
+  const Field3D toFieldAligned(const Field3D &f, const REGION UNUSED(region)) override {
     return f;
   }
   
@@ -73,7 +76,7 @@ public:
    * The field is already aligned in Y, so this
    * does nothing
    */
-  const Field3D fromFieldAligned(const Field3D &f) override {
+  const Field3D fromFieldAligned(const Field3D &f, const REGION UNUSED(region)) override {
     return f;
   }
 
@@ -109,13 +112,13 @@ public:
    * in X-Z, and the metric tensor will need to be changed 
    * if X derivatives are used.
    */
-  const Field3D toFieldAligned(const Field3D &f) override;
+  const Field3D toFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
   /*!
    * Converts a field back to X-Z orthogonal coordinates
    * from field aligned coordinates.
    */
-  const Field3D fromFieldAligned(const Field3D &f) override;
+  const Field3D fromFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
   bool canToFromFieldAligned() override{
     return true;
@@ -142,7 +145,7 @@ private:
    * Shift a 2D field in Z. 
    * Since 2D fields are constant in Z, this has no effect
    */
-  const Field2D shiftZ(const Field2D &f, const Field2D &UNUSED(zangle)){return f;};
+  const Field2D shiftZ(const Field2D &f, const Field2D &UNUSED(zangle), const REGION UNUSED(region)=RGN_NOX){return f;};
 
   /*!
    * Shift a 3D field \p f in Z by the given \p zangle
@@ -151,7 +154,7 @@ private:
    * @param[in] zangle   Toroidal angle (z)
    *
    */ 
-  const Field3D shiftZ(const Field3D &f, const Field2D &zangle);
+  const Field3D shiftZ(const Field3D &f, const Field2D &zangle, const REGION region=RGN_NOX);
 
   /*!
    * Shift a 3D field \p f by the given phase \p phs in Z
@@ -162,7 +165,7 @@ private:
    * @param[in] f  The field to shift
    * @param[in] phs  The phase to shift by
    */
-  const Field3D shiftZ(const Field3D &f, const arr3Dvec &phs);
+  const Field3D shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region=RGN_NOX);
 
   /*!
    * Shift a given 1D array, assumed to be in Z, by the given \p zangle
diff --git a/include/field3d.hxx b/include/field3d.hxx
index e9f090bea6..a8a2e8ee74 100644
--- a/include/field3d.hxx
+++ b/include/field3d.hxx
@@ -223,10 +223,23 @@ class Field3D : public Field, public FieldData {
    * Ensure that yup and ydown refer to this field
    */
   void mergeYupYdown();
+
+  /*!
+   * Delete all yup/ydown fields and field_fa
+   */
+  void deleteYupYdown();
   
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
-    return (yup_field != nullptr) && (ydown_field != nullptr);
+    return hasValidYUpDown && (yup_field != nullptr) && (ydown_field != nullptr);
+  }
+
+  /// Set validity state of yup/ydown fields
+  void setHasValidYUpDown(bool set) {
+    if (set) {
+      ASSERT1(yup_field != nullptr && ydown_field != nullptr);
+    }
+    hasValidYUpDown = set;
   }
 
   /// Return reference to yup field
@@ -457,6 +470,9 @@ private:
   
   Field3D *deriv; ///< Time derivative (may be NULL)
 
+  /// flag to record whether yup_field and ydown_field have been set correctly
+  bool hasValidYUpDown;
+
   /// Pointers to fields containing values along Y
   Field3D *yup_field, *ydown_field;
 };
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 7d9dda7cea..624878ba46 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -146,12 +146,8 @@ Field3D::~Field3D() {
     // Now delete them as part of the deriv vector
     delete deriv;
   }
-  
-  if((yup_field != this) && (yup_field != nullptr))
-    delete yup_field;
-  
-  if((ydown_field != this) && (ydown_field != nullptr))
-    delete ydown_field;
+
+  deleteYupYdown();
 }
 
 void Field3D::allocate() {
@@ -195,18 +191,25 @@ void Field3D::mergeYupYdown() {
   if(yup_field == this && ydown_field == this)
     return;
 
-  if(yup_field != nullptr){
-    delete yup_field;
-  }
-
-  if(ydown_field != nullptr) {
-    delete ydown_field;
-  }
+  deleteYupYdown();
 
   yup_field = this;
   ydown_field = this;
 }
 
+void Field3D::deleteYupYdown() {
+  // Delete auxiliary fields if they have been set
+  if (yup_field == this && ydown_field == this) {
+    return;
+  }
+
+  delete yup_field;
+  yup_field = nullptr;
+
+  delete ydown_field;
+  ydown_field = nullptr;
+}
+
 Field3D& Field3D::ynext(int dir) {
   switch(dir) {
   case +1:
@@ -307,6 +310,7 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
 
   setLocation(rhs.location);
 
+  setHasValidYUpDown(false);
   return *this;
 }
 
@@ -328,7 +332,9 @@ Field3D & Field3D::operator=(const Field2D &rhs) {
   
   /// Only 3D fields have locations for now
   //location = CELL_CENTRE;
-  
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -347,6 +353,8 @@ void Field3D::operator=(const FieldPerp &rhs) {
   BOUT_FOR(i, region_all) {
     (*this)(i, rhs.getIndex()) = rhs[i];
   }
+
+  setHasValidYUpDown(false);
 }
 
 Field3D & Field3D::operator=(const BoutReal val) {
@@ -366,6 +374,8 @@ Field3D & Field3D::operator=(const BoutReal val) {
   //location = CELL_CENTRE;
   // DON'T RE-SET LOCATION
 
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja
index af3c3afb90..adf9544cb1 100644
--- a/src/field/gen_fieldops.jinja
+++ b/src/field/gen_fieldops.jinja
@@ -95,6 +95,11 @@
   } else {
     (*this) = (*this) {{operator}} {{rhs.name}};
   }
+
+  {% if (out == "Field3D") %}
+    setHasValidYUpDown(false);
+  {% endif %}
+
   return *this;
 }
 {% endif %}
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 0456dbef1f..e030cb70b8 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -62,6 +62,9 @@ Field3D &Field3D::operator*=(const Field3D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -121,6 +124,9 @@ Field3D &Field3D::operator/=(const Field3D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -180,6 +186,9 @@ Field3D &Field3D::operator+=(const Field3D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -239,6 +248,9 @@ Field3D &Field3D::operator-=(const Field3D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -306,6 +318,9 @@ Field3D &Field3D::operator*=(const Field2D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -375,6 +390,9 @@ Field3D &Field3D::operator/=(const Field2D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -442,6 +460,9 @@ Field3D &Field3D::operator+=(const Field2D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -509,6 +530,9 @@ Field3D &Field3D::operator-=(const Field2D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -546,6 +570,9 @@ Field3D &Field3D::operator*=(const BoutReal rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -583,6 +610,9 @@ Field3D &Field3D::operator/=(const BoutReal rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -620,6 +650,9 @@ Field3D &Field3D::operator+=(const BoutReal rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -657,6 +690,9 @@ Field3D &Field3D::operator-=(const BoutReal rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -844,6 +880,7 @@ Field2D &Field2D::operator*=(const Field2D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
   return *this;
 }
 
@@ -903,6 +940,7 @@ Field2D &Field2D::operator/=(const Field2D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
   return *this;
 }
 
@@ -962,6 +1000,7 @@ Field2D &Field2D::operator+=(const Field2D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
   return *this;
 }
 
@@ -1021,6 +1060,7 @@ Field2D &Field2D::operator-=(const Field2D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
   return *this;
 }
 
@@ -1058,6 +1098,7 @@ Field2D &Field2D::operator*=(const BoutReal rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
   return *this;
 }
 
@@ -1095,6 +1136,7 @@ Field2D &Field2D::operator/=(const BoutReal rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
   return *this;
 }
 
@@ -1132,6 +1174,7 @@ Field2D &Field2D::operator+=(const BoutReal rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
   return *this;
 }
 
@@ -1169,6 +1212,7 @@ Field2D &Field2D::operator-=(const BoutReal rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
   return *this;
 }
 
diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index f684cce292..1b6ba045f9 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2201,7 +2201,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         }
       }
 
-      result = this->fromFieldAligned(result);
+      result = this->fromFieldAligned(result, RGN_NOBNDRY);
     }
   } else {
     // Non-staggered case
diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx
index b3b9285811..91ce7dcbf9 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -151,11 +151,9 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
           // coordinates
 
           Field3D var_fa = fieldmesh->toFieldAligned(var);
-          Field3D result_fa(fieldmesh);
           if (region != RGN_NOBNDRY) {
-            result_fa = fieldmesh->toFieldAligned(result);
+            result = fieldmesh->toFieldAligned(result);
           }
-          result_fa.allocate();
           if (fieldmesh->ystart > 1) {
 
             // More than one guard cell, so set pp and mm values
@@ -180,7 +178,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
                   s.m = s.c;
                 }
 
-                result_fa[i] = interp(s);
+                result[i] = interp(s);
               }
             }
           } else {
@@ -206,12 +204,12 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
                   s.m = s.c;
                 }
 
-                result_fa[i] = interp(s);
+                result[i] = interp(s);
               }
             }
           }
           
-          result = fieldmesh->fromFieldAligned(result_fa);
+          result = fieldmesh->fromFieldAligned(result, RGN_NOBNDRY);
         }
         break;
       }
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index 25002ffed3..018b5cdb26 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -302,6 +302,8 @@ void FCITransform::calcYUpDown(Field3D &f) {
   // Interpolate f onto yup and ydown fields
   f.ynext(forward_map.dir) = forward_map.interpolate(f);
   f.ynext(backward_map.dir) = backward_map.interpolate(f);
+
+  f.setHasValidYUpDown(true);
 }
 
 void FCITransform::integrateYUpDown(Field3D &f) {
@@ -313,4 +315,6 @@ void FCITransform::integrateYUpDown(Field3D &f) {
   // Integrate f onto yup and ydown fields
   f.ynext(forward_map.dir) = forward_map.integrate(f);
   f.ynext(backward_map.dir) = backward_map.integrate(f);
+
+  f.setHasValidYUpDown(true);
 }
diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx
index 293113b91d..747f311c03 100644
--- a/src/mesh/parallel/fci.hxx
+++ b/src/mesh/parallel/fci.hxx
@@ -73,11 +73,11 @@ public:
   
   void integrateYUpDown(Field3D &f) override;
   
-  const Field3D toFieldAligned(const Field3D &UNUSED(f)) override {
+  const Field3D toFieldAligned(const Field3D &UNUSED(f), const REGION UNUSED(region)) override {
     throw BoutException("FCI method cannot transform into field aligned grid");
   }
 
-  const Field3D fromFieldAligned(const Field3D &UNUSED(f)) override {
+  const Field3D fromFieldAligned(const Field3D &UNUSED(f), const REGION UNUSED(region)) override {
     throw BoutException("FCI method cannot transform into field aligned grid");
   }
 
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 80fdc37274..7dbf6382ff 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -123,36 +123,39 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
       shiftZ(&(f(jx,jy-1,0)), ydownPhs[jx][jy], &(ydown(jx,jy-1,0)));
     }
   }
+
+  f.setHasValidYUpDown(true);
 }
   
 /*!
  * Shift the field so that X-Z is not orthogonal,
  * and Y is then field aligned.
  */
-const Field3D ShiftedMetric::toFieldAligned(const Field3D &f) {
-  return shiftZ(f, toAlignedPhs);
+const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) {
+  return shiftZ(f, toAlignedPhs, region);
 }
 
 /*!
  * Shift back, so that X-Z is orthogonal,
  * but Y is not field aligned.
  */
-const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f) {
-  return shiftZ(f, fromAlignedPhs);
+const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) {
+  return shiftZ(f, fromAlignedPhs, region);
 }
 
-const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs) {
+const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region) {
   ASSERT1(&mesh == f.getMesh());
+  ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here
+  ASSERT1(f.getLocation() == CELL_CENTRE); // only have zShift for CELL_CENTRE, so can only deal with CELL_CENTRE inputs
   if(mesh.LocalNz == 1)
     return f; // Shifting makes no difference
 
   Field3D result(&mesh);
   result.allocate();
-  
-  for(int jx=0;jx<mesh.LocalNx;jx++) {
-    for(int jy=0;jy<mesh.LocalNy;jy++) {
-      shiftZ(f(jx,jy), phs[jx][jy], result(jx,jy));
-    }
+  result.setLocation(f.getLocation());
+
+  for (const auto &i : mesh.getRegion2D(REGION_STRING(region))) {
+    shiftZ(f(i.x(),i.y()), phs[i.x()][i.y()], result(i.x(),i.y()));
   }
   
   return result;
@@ -177,18 +180,25 @@ void ShiftedMetric::shiftZ(const BoutReal *in, const std::vector<dcomplex> &phs,
 }
 
 //Old approach retained so we can still specify a general zShift
-const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle) {
+const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, const REGION region) {
   ASSERT1(&mesh == f.getMesh());
+  ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here
+  ASSERT1(f.getLocation() == zangle.getLocation());
   if(mesh.LocalNz == 1)
     return f; // Shifting makes no difference
 
   Field3D result(&mesh);
   result.allocate();
-
-  for(int jx=0;jx<mesh.LocalNx;jx++) {
-    for(int jy=0;jy<mesh.LocalNy;jy++) {
-      shiftZ(f(jx,jy), mesh.LocalNz, zangle(jx,jy), result(jx,jy));
-    }
+  result.setLocation(f.getLocation());
+
+  // We only use methods in ShiftedMetric to get fields for parallel operations
+  // like interp_to or DDY.
+  // Therefore we don't need x-guard cells, so do not set them.
+  // (Note valgrind complains about corner guard cells if we try to loop over
+  // the whole grid, because zShift is not initialized in the corner guard
+  // cells.)
+  for(const auto &i : mesh.getRegion2D(REGION_STRING(region))) {
+    shiftZ(f(i.x(), i.y()), mesh.LocalNz, zangle(i.x(),i.y()), result(i.x(), i.y()));
   }
   
   return result;
diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx
index 2267c2f06f..c51c7ca52c 100644
--- a/src/solver/solver.cxx
+++ b/src/solver/solver.cxx
@@ -982,6 +982,7 @@ void Solver::load_vars(BoutReal *udata) {
   for(const auto& f : f3d) {
     f.var->allocate();
     f.var->setLocation(f.location);
+    f.var->setHasValidYUpDown(false); // flag yup/ydown fields as invalid f.var is updated
   }
 
   loop_vars(udata, LOAD_VARS);
diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx
index 649a3abb1b..b6f5e6cd35 100644
--- a/tests/unit/field/test_field3d.cxx
+++ b/tests/unit/field/test_field3d.cxx
@@ -215,6 +215,7 @@ TEST_F(Field3DTest, SplitYupYDown) {
   EXPECT_FALSE(field.hasYupYdown());
 
   field.splitYupYdown();
+  field.setHasValidYUpDown(true);
 
   EXPECT_TRUE(field.hasYupYdown());
 
@@ -242,6 +243,7 @@ TEST_F(Field3DTest, MergeYupYDown) {
   EXPECT_FALSE(field.hasYupYdown());
 
   field.mergeYupYdown();
+  field.setHasValidYUpDown(true);
 
   EXPECT_TRUE(field.hasYupYdown());