diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 1524923515..4067835e30 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -609,6 +609,12 @@ 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, REGION region = RGN_NOX) {
+    getParallelTransform().calcYUpDown(f, region);
+  }
+
   /// Transform a field into field-aligned coordinates
   const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
     return getParallelTransform().toFieldAligned(f, region);
diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index bab43dd61e..85776aaee7 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -10,6 +10,7 @@
 #include <boutexception.hxx>
 #include <dcomplex.hxx>
 #include <unused.hxx>
+#include <utils.hxx>
 
 class Mesh;
 
@@ -28,7 +29,7 @@ public:
   virtual ~ParallelTransform() {}
 
   /// Given a 3D field, calculate and set the Y up down fields
-  virtual void calcYUpDown(Field3D &f) = 0;
+  virtual void calcYUpDown(Field3D &f, REGION region = RGN_NOX) = 0;
 
   /// Calculate Yup and Ydown fields by integrating over mapped points
   /// This should be used for parallel divergence operators
@@ -59,7 +60,7 @@ 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, REGION UNUSED(region) = RGN_NOX) override {f.mergeYupYdown();}
   
   /*!
    * The field is already aligned in Y, so this
@@ -99,7 +100,7 @@ public:
    * Calculates the yup() and ydown() fields of f
    * by taking FFTs in Z and applying a phase shift.
    */ 
-  void calcYUpDown(Field3D &f) override;
+  void calcYUpDown(Field3D &f, REGION UNUSED(region) = RGN_NOX) override;
   
   /*!
    * Uses FFTs and a phase shift to align the grid points
@@ -184,5 +185,115 @@ private:
   void shiftZ(const BoutReal *in, const std::vector<dcomplex> &phs, BoutReal *out);
 };
 
+/*!
+ * Alternative shifted metric method
+ * Fields are stored on a grid is orthogonal in X-Z, interpolated onto a
+ * field-aligned grid to calculate parallel derivatives or interpolations
+ *
+ * In this implementation the interpolation is done using FFTs in Z
+ */
+class ShiftToFieldAligned : public ParallelTransform {
+public:
+  ShiftToFieldAligned(Mesh &mesh);
+
+  /*!
+   * Do not calculate yup() and ydown() fields of f
+   * for this method.
+   * Instead use this method to calculate f.field_fa.
+   * This is a bit hacky, but we can rename ParallelTransform::calcYUpDown to
+   * something more generic later.
+   */
+  void calcYUpDown(Field3D &f, REGION region = RGN_NOX) override;
+
+  /*!
+   * Uses FFTs and a phase shift to align the grid points
+   * with the y coordinate (along magnetic field usually).
+   *
+   * Note that the returned field will no longer be orthogonal
+   * in X-Z, and the metric tensor will need to be changed
+   * if X derivatives are used.
+   */
+  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, const REGION region=RGN_NOX) override;
+
+  bool canToFromFieldAligned() override{
+    return true;
+  }
+
+private:
+  ShiftToFieldAligned();
+
+  Mesh &mesh; ///< The mesh this paralleltransform is part of
+
+  /// This is the shift in toroidal angle (z) which takes a point from
+  /// X-Z orthogonal to field-aligned along Y.
+  Field2D zShift_CENTRE, zShift_XLOW, zShift_YLOW; ///< The angles to shift from orthogonal to field-aligned coordinates
+  Array<dcomplex> cmplx; ///< A temporary array, used for input/output to fft routines
+  Array<dcomplex> cmplxLoc; ///< A temporary array, used for input/output to fft routines
+
+  Matrix< Array<dcomplex> > getToAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  Matrix< Array<dcomplex> > getFromAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+
+  bool has_toAligned_CENTRE, has_toAligned_XLOW, has_toAligned_YLOW; ///< Flags saying whether phases for shifts to field-aligned coordinates have been calculated.
+  bool has_fromAligned_CENTRE, has_fromAligned_XLOW, has_fromAligned_YLOW; ///< Flags saying whether phases for shifts from field-aligned coordinates have been calculated.
+
+  Matrix< Array<dcomplex> > toAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Cell centre version.
+  Matrix< Array<dcomplex> > fromAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Cell centre version.
+  Matrix< Array<dcomplex> > toAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > fromAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > toAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_YLOW.
+  Matrix< Array<dcomplex> > fromAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_YLOW.
+
+  /*!
+   * 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), const REGION UNUSED(region)=RGN_NOX){return f;};
+
+  /*!
+   * Shift a 3D field \p f in Z by the given \p zangle
+   *
+   * @param[in] f  The field to shift
+   * @param[in] zangle   Toroidal angle (z)
+   *
+   */
+  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
+   *
+   * Calculates FFT in Z, multiplies by the complex phase
+   * and inverse FFTS.
+   *
+   * @param[in] f  The field to shift
+   * @param[in] phs  The phase to shift by
+   */
+  const Field3D shiftZ(const Field3D &f, const Matrix< Array<dcomplex> > &phs, const REGION region=RGN_NOX);
+
+  /*!
+   * Shift a given 1D array, assumed to be in Z, by the given \p zangle
+   *
+   * @param[in] in  A 1D array of length \p len
+   * @param[in] len  Length of the in and out arrays
+   * @param[in] zangle  The angle (z coordinate) to shift by
+   * @param[out] out  A 1D array of length \p len, already allocated
+   */
+  void shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutReal *out);
+
+  /*!
+   * Shift a given 1D array, assumed to be in Z, by the given \p zangle
+   *
+   * @param[in] in  A 1D array of length mesh.LocalNz
+   * @param[in] phs Phase shift, assumed to have length (mesh.LocalNz/2 + 1) i.e. the number of modes
+   * @param[out] out  A 1D array of length mesh.LocalNz, already allocated
+   */
+  void shiftZ(const BoutReal *in, const Array<dcomplex> &phs, BoutReal *out);
+};
+
 
 #endif // __PARALLELTRANSFORM_H__
diff --git a/include/field3d.hxx b/include/field3d.hxx
index 8d4ad7f9ec..06f18337a5 100644
--- a/include/field3d.hxx
+++ b/include/field3d.hxx
@@ -154,6 +154,11 @@ class Mesh;  // #include "bout/mesh.hxx"
 
       f.yup()(0,1,0) // ok
 
+  It is also useful to be able to access a version of the field interpolated
+  (in the z-direction) onto a field-aligned grid. This can be
+  calculated/accessed with the Mesh::toFieldAligned() and the main field can be
+  interpolated from a field-aligned one with Mesh::fromFieldAligned().
+
  */
 class Field3D : public Field, public FieldData {
  public:
@@ -234,6 +239,19 @@ class Field3D : public Field, public FieldData {
     return (yup_field != nullptr) && (ydown_field != nullptr);
   }
 
+  /// Check if this field has a field-aligned version
+  bool hasFieldAligned() const {
+    return has_field_aligned;
+  }
+
+  /// Set has_field_aligned
+  void setHasFieldAligned(bool val) {
+    if (val) {
+      ASSERT1(field_fa != nullptr && field_fa->isAllocated());
+    }
+    has_field_aligned = val;
+  }
+
   /// Return reference to yup field
   Field3D& yup() { 
     ASSERT2(yup_field != nullptr); // Check for communicate
@@ -261,6 +279,12 @@ class Field3D : public Field, public FieldData {
   Field3D& ynext(int dir);
   const Field3D& ynext(int dir) const;
 
+  /// Return reference to field-aligned field
+  Field3D& fieldAligned();
+
+  /// Return const reference to field-aligned field
+  const Field3D& fieldAligned() const;
+
   /// Set variable location for staggered grids to @param new_location
   ///
   /// Throws BoutException if new_location is not `CELL_CENTRE` and
@@ -464,6 +488,10 @@ private:
 
   /// Pointers to fields containing values along Y
   Field3D *yup_field, *ydown_field;
+
+  /// Pointer to field containing the field-aligned version of the field
+  Field3D *field_fa;
+  bool has_field_aligned;
 };
 
 // Non-member overloaded operators
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index eb117cf734..697b098fd2 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -46,7 +46,8 @@
 /// Constructor
 Field3D::Field3D(Mesh *localmesh)
     : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr),
-      ydown_field(nullptr) {
+      ydown_field(nullptr),
+      field_fa(nullptr), has_field_aligned(false) {
 #ifdef TRACK
   name = "<F3D>";
 #endif
@@ -72,7 +73,8 @@ Field3D::Field3D(Mesh *localmesh)
 Field3D::Field3D(const Field3D &f)
     : Field(f.fieldmesh),                // The mesh containing array sizes
       background(nullptr), data(f.data), // This handles references to the data array
-      deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) {
+      deriv(nullptr), yup_field(nullptr), ydown_field(nullptr),
+      field_fa(nullptr), has_field_aligned(false) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -101,7 +103,8 @@ Field3D::Field3D(const Field3D &f)
 
 Field3D::Field3D(const Field2D &f)
     : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup_field(nullptr),
-      ydown_field(nullptr) {
+      ydown_field(nullptr),
+      field_fa(nullptr), has_field_aligned(false) {
 
   TRACE("Field3D: Copy constructor from Field2D");
 
@@ -119,7 +122,8 @@ Field3D::Field3D(const Field2D &f)
 
 Field3D::Field3D(const BoutReal val, Mesh *localmesh)
     : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr),
-      ydown_field(nullptr) {
+      ydown_field(nullptr),
+      field_fa(nullptr), has_field_aligned(false) {
 
   TRACE("Field3D: Copy constructor from value");
 
@@ -199,15 +203,20 @@ void Field3D::mergeYupYdown() {
 
 void Field3D::deleteYupYdown() {
   // Delete auxiliary fields if they have been set
-  if (yup_field == this && ydown_field == this) {
-    return;
-  }
+  if (!(yup_field == this && ydown_field == this)) {
+
+    delete yup_field;
+    yup_field = nullptr;
 
-  delete yup_field;
-  yup_field = nullptr;
+    delete ydown_field;
+    ydown_field = nullptr;
+  }
 
-  delete ydown_field;
-  ydown_field = nullptr;
+  if (field_fa != this) {
+    delete field_fa;
+  }
+  field_fa = nullptr;
+  has_field_aligned = false;
 }
 
 Field3D& Field3D::ynext(int dir) {
@@ -232,6 +241,25 @@ const Field3D& Field3D::ynext(int dir) const {
   }
 }
 
+
+Field3D& Field3D::fieldAligned() {
+  if (field_fa == nullptr) {
+    // has_field_aligned should be false now, until field_fa is set to a value.
+    // This should have been set either in the constructor or when field_fa was
+    // deleted.
+    ASSERT1(!has_field_aligned);
+
+    field_fa = new Field3D(getMesh());
+    field_fa->setLocation(location);
+  }
+  return *field_fa;
+}
+
+const Field3D& Field3D::fieldAligned() const {
+  ASSERT1(field_fa != nullptr);
+  return *field_fa;
+}
+
 void Field3D::setLocation(CELL_LOC new_location) {
   if (getMesh()->StaggerGrids) {
     if (new_location == CELL_VSHIFT) {
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 81c0e921e4..d120d46728 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -702,26 +702,30 @@ const Field3D Coordinates::Div_par(const Field3D &f, CELL_LOC outloc,
   
   // Need Bxy at location of f, which might be different from location of this
   // Coordinates object
-  Field2D Bxy_floc = f.getCoordinates()->Bxy;
-
-  if (!f.hasYupYdown()) {
-    // No yup/ydown fields. The Grad_par operator will
-    // shift to field aligned coordinates
-    return Bxy * Grad_par(f / Bxy_floc, outloc, method);
-  }
-
-  // Need to modify yup and ydown fields
+  Field2D& Bxy_floc = f.getCoordinates()->Bxy;
   Field3D f_B = f / Bxy_floc;
-  if (&f.yup() == &f) {
-    // Identity, yup and ydown point to same field
-    f_B.mergeYupYdown();
+
+  if (f.hasYupYdown()) {
+    // Need to modify yup and ydown fields
+    if (&f.yup() == &f) {
+      // Identity, yup and ydown point to same field
+      f_B.mergeYupYdown();
+    } else {
+      // Distinct fields
+      f_B.splitYupYdown();
+      f_B.yup() = f.yup() / Bxy_floc;
+      f_B.ydown() = f.ydown() / Bxy_floc;
+    }
+    return Bxy * Grad_par(f_B, outloc, method);
+  } else if (f.hasFieldAligned()) {
+    f_B.fieldAligned() = f.fieldAligned() / Bxy;
+    f_B.setHasFieldAligned(true);
+    return Bxy * Grad_par(f_B, outloc, method);
   } else {
-    // Distinct fields
-    f_B.splitYupYdown();
-    f_B.yup() = f.yup() / Bxy_floc;
-    f_B.ydown() = f.ydown() / Bxy_floc;
+    // No yup/ydown or field-aligned fields. The Grad_par operator will shift
+    // to field aligned coordinates
+    return Bxy * Grad_par(f_B, outloc, method);
   }
-  return Bxy * Grad_par(f_B, outloc, method);
 }
 
 /////////////////////////////////////////////////////////
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index 70112cfe91..97f4a4816f 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -292,15 +292,19 @@ void Mesh::setParallelTransform() {
 
   // Convert to lower case for comparison
   ptstr = lowercase(ptstr);
-    
+
   if(ptstr == "identity") {
     // Identity method i.e. no transform needed
     transform = std::unique_ptr<ParallelTransform>(new ParallelTransformIdentity());
-      
+
   }else if(ptstr == "shifted") {
     // Shifted metric method
     transform = std::unique_ptr<ParallelTransform>(new ShiftedMetric(*this));
-      
+
+  }else if(ptstr == "shifttofieldaligned") {
+    // Alternative shifted metric method
+    transform = std::unique_ptr<ParallelTransform>(new ShiftToFieldAligned(*this));
+
   }else if(ptstr == "fci") {
 
     Options *fci_options = Options::getRoot()->getSection("fci");
@@ -308,7 +312,7 @@ void Mesh::setParallelTransform() {
     bool fci_zperiodic;
     fci_options->get("z_periodic", fci_zperiodic, true);
     transform = std::unique_ptr<ParallelTransform>(new FCITransform(*this, fci_zperiodic));
-      
+
   }else {
     throw BoutException("Unrecognised paralleltransform option.\n"
                         "Valid choices are 'identity', 'shifted', 'fci'");
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index 25002ffed3..00c7b64b20 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -293,7 +293,7 @@ const Field3D FCIMap::integrate(Field3D &f) const {
   return result;
 }
 
-void FCITransform::calcYUpDown(Field3D &f) {
+void FCITransform::calcYUpDown(Field3D &f, REGION region) {
   TRACE("FCITransform::calcYUpDown");
 
   // Ensure that yup and ydown are different fields
diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx
index 747f311c03..a8ec08da21 100644
--- a/src/mesh/parallel/fci.hxx
+++ b/src/mesh/parallel/fci.hxx
@@ -69,7 +69,7 @@ public:
       : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic),
         zperiodic(zperiodic) {}
 
-  void calcYUpDown(Field3D &f) override;
+  void calcYUpDown(Field3D &f, REGION UNUSED(region) = RGN_NOX) override;
   
   void integrateYUpDown(Field3D &f) override;
   
diff --git a/src/mesh/parallel/makefile b/src/mesh/parallel/makefile
index 9c4817b008..97728dd3c7 100644
--- a/src/mesh/parallel/makefile
+++ b/src/mesh/parallel/makefile
@@ -2,7 +2,7 @@
 BOUT_TOP = ../../..
 
 DIRS            = 
-SOURCEC		= shiftedmetric.cxx fci.cxx
+SOURCEC		= shiftedmetric.cxx fci.cxx shifttofieldaligned.cxx
 TARGET		= lib
 
 include $(BOUT_TOP)/make.config
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 53a5fdea96..04767f7aed 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -103,7 +103,7 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
 /*!
  * Calculate the Y up and down fields
  */
-void ShiftedMetric::calcYUpDown(Field3D &f) {
+void ShiftedMetric::calcYUpDown(Field3D &f, REGION region) {
   f.splitYupYdown();
   
   Field3D& yup = f.yup();
diff --git a/src/mesh/parallel/shifttofieldaligned.cxx b/src/mesh/parallel/shifttofieldaligned.cxx
new file mode 100644
index 0000000000..7c057a116e
--- /dev/null
+++ b/src/mesh/parallel/shifttofieldaligned.cxx
@@ -0,0 +1,380 @@
+/*
+ * Implements the shifted metric method for parallel derivatives
+ *
+ * By default fields are stored so that X-Z are orthogonal,
+ * and so not aligned in Y.
+ *
+ */
+
+#include <bout/paralleltransform.hxx>
+#include <bout/mesh.hxx>
+#include <interpolation.hxx>
+#include <fft.hxx>
+#include <bout/constants.hxx>
+
+#include <cmath>
+#include <msg_stack.hxx>
+
+#include <output.hxx>
+
+ShiftToFieldAligned::ShiftToFieldAligned(Mesh &m) :
+  mesh(m), zShift_CENTRE(&m), zShift_XLOW(&m), zShift_YLOW(&m),
+  has_toAligned_CENTRE(false), has_toAligned_XLOW(false), has_toAligned_YLOW(false),
+  has_fromAligned_CENTRE(false), has_fromAligned_XLOW(false), has_fromAligned_YLOW(false) {
+
+  // Read the zShift angle from the mesh
+  if(mesh.get(zShift_CENTRE, "zShift")) {
+    // No zShift variable. Try qinty in BOUT grid files
+    mesh.get(zShift_CENTRE, "qinty");
+  }
+
+  if (mesh.xstart >=2) {
+    // Can interpolate in x-direction
+    // Calculate staggered field for zShift and apply boundary conditions
+    zShift_XLOW = interp_to(zShift_CENTRE, CELL_XLOW, RGN_ALL);
+    zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value
+  }
+  if (mesh.ystart >=2) {
+    // Can interpolate in y-direction
+    // Calculate staggered field for zShift and apply boundary conditions
+    zShift_YLOW = interp_to(zShift_CENTRE, CELL_YLOW, RGN_ALL);
+    zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value
+  }
+
+  int nmodes = mesh.LocalNz/2 + 1;
+  //Allocate storage for complex intermediate
+  cmplx = Array<dcomplex>(nmodes);
+  std::fill(cmplx.begin(), cmplx.end(), 0.0);
+}
+
+//As we're attached to a mesh we can expect the z direction to not change
+//once we've been created so cache the complex phases used in transformations
+//the first time they are needed
+Matrix< Array<dcomplex> > ShiftToFieldAligned::getFromAlignedPhs(CELL_LOC location) {
+  switch (location) {
+  case CELL_CENTRE: {
+    if (!has_toAligned_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_toAligned_CENTRE = true;
+      fromAlignedPhs_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : fromAlignedPhs_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            fromAlignedPhs_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*zShift_CENTRE(jx,jy)) , -sin(kwave*zShift_CENTRE(jx,jy)));
+          }
+        }
+      }
+    }
+    return fromAlignedPhs_CENTRE;
+    break;
+  }
+  case CELL_XLOW: {
+    if (!has_toAligned_XLOW) {
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_toAligned_XLOW = true;
+      fromAlignedPhs_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : fromAlignedPhs_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            fromAlignedPhs_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_XLOW(jx,jy)), -sin(kwave*zShift_XLOW(jx,jy)));
+          }
+        }
+      }
+    }
+    return fromAlignedPhs_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (!has_toAligned_YLOW) {
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_toAligned_YLOW = true;
+      fromAlignedPhs_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : fromAlignedPhs_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            fromAlignedPhs_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_YLOW(jx,jy)), -sin(kwave*zShift_YLOW(jx,jy)));
+          }
+        }
+      }
+    }
+    return fromAlignedPhs_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getFromAlignedPhs(CELL_CENTRE);
+    break;
+  }
+  default: {
+    // This should never happen
+    throw BoutException("Unsupported stagger of phase shifts\n"
+                        " - don't know how to interpolate to %s",strLocation(location));
+    break;
+  }
+  };
+}
+
+Matrix< Array<dcomplex> > ShiftToFieldAligned::getToAlignedPhs(CELL_LOC location) {
+  switch (location) {
+  case CELL_CENTRE: {
+    if (!has_fromAligned_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_fromAligned_CENTRE = true;
+      toAlignedPhs_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : toAlignedPhs_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            toAlignedPhs_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*zShift_CENTRE(jx,jy)), sin(kwave*zShift_CENTRE(jx,jy)));
+          }
+        }
+      }
+    }
+    return toAlignedPhs_CENTRE;
+    break;
+  }
+  case CELL_XLOW: {
+    if (!has_fromAligned_XLOW) {
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_fromAligned_XLOW = true;
+      toAlignedPhs_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : toAlignedPhs_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            toAlignedPhs_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_XLOW(jx,jy)), sin(kwave*zShift_XLOW(jx,jy)));
+          }
+        }
+      }
+    }
+    return toAlignedPhs_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (!has_fromAligned_YLOW) {
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.getCoordinates()->zlength();
+
+      has_fromAligned_YLOW = true;
+      toAlignedPhs_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : toAlignedPhs_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
+      }
+
+      //To/From field aligned phases
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
+        for(int jy=0;jy<mesh.LocalNy;jy++) {
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            toAlignedPhs_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_YLOW(jx,jy)), sin(kwave*zShift_YLOW(jx,jy)));
+          }
+        }
+      }
+    }
+    return toAlignedPhs_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getToAlignedPhs(CELL_CENTRE);
+    break;
+  }
+  default: {
+    // This should never happen
+    throw BoutException("Unsupported stagger of phase shifts\n"
+                        " - don't know how to interpolate to %s",strLocation(location));
+    break;
+  }
+  };
+}
+
+
+/*!
+ * Calculate the field-aligned field
+ */
+void ShiftToFieldAligned::calcYUpDown(Field3D &f, REGION region) {
+  ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY);
+  ASSERT1(&mesh == f.getMesh());
+
+  // We only use methods in ShiftToFieldAligned 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.)
+
+  Field3D& f_fa = f.fieldAligned();
+  f_fa = shiftZ(f, getToAlignedPhs(f.getLocation()), region);
+  f.setHasFieldAligned(true);
+}
+
+/*!
+ * Get the shifted field so that X-Z is not orthogonal,
+ * and Y is then field aligned.
+ */
+const Field3D ShiftToFieldAligned::toFieldAligned(const Field3D &f, const REGION region) {
+  ASSERT1(region==RGN_NOX || region==RGN_NOBNDRY);
+  if (f.hasFieldAligned()) {
+    return f.fieldAligned();
+  } else {
+#if CHECK > 1
+    static int count = 0;
+
+    if (count<100) {
+      // Get a trace of where we were called from
+      std::string message = msg_stack.getDump();
+
+      output<<"Warning:"<<endl;
+      output<<"Called toFieldAligned on a field without a field_fa member (i.e. f.hasFieldAligned()==false)."<<endl;
+      output<<"This may be suboptimal."<<endl;
+      output<<"Was called from:"<<endl;
+      output<<message<<endl;
+
+      count++;
+    } else if (count == 100) {
+      output<<"More than 100 warnings from ShiftToFieldAligned::toFieldAligned."<<endl
+        <<"Suppressing further output from here."<<endl;
+      count++;
+    }
+#endif
+    return shiftZ(f, getToAlignedPhs(f.getLocation()), region);
+  }
+}
+
+/*!
+ * Shift back, so that X-Z is orthogonal,
+ * but Y is not field aligned.
+ */
+const Field3D ShiftToFieldAligned::fromFieldAligned(const Field3D &f, const REGION region) {
+  return shiftZ(f, getFromAlignedPhs(f.getLocation()), region);
+}
+
+const Field3D ShiftToFieldAligned::shiftZ(const Field3D &f, const Matrix< Array<dcomplex> > &phs, const REGION region) {
+  ASSERT1(&mesh == f.getMesh());
+  ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here
+  if(mesh.LocalNz == 1)
+    return f; // Shifting makes no difference
+
+  Field3D result(&mesh);
+  result.allocate();
+  result.setLocation(f.getLocation());
+  //result = 0.; // Set to value to avoid uninitialized value errors from Valgrind
+
+  invalidateGuards(result); // Won't set x-guard cells, so allow checking to throw exception if they are used.
+
+  // We only use methods in ShiftToFieldAligned 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()), phs(i.x(), i.y()), result(i.x(), i.y()));
+  }
+
+  return result;
+
+}
+
+void ShiftToFieldAligned::shiftZ(const BoutReal *in, const Array<dcomplex> &phs, BoutReal *out) {
+  // Take forward FFT
+  rfft(in, mesh.LocalNz, cmplx.begin());
+
+  //Following is an algorithm approach to write a = a*b where a and b are
+  //vectors of dcomplex.
+  //  std::transform(cmplxOneOff.begin(),cmplxOneOff.end(), ptr.begin(),
+  //                 cmplxOneOff.begin(), std::multiplies<dcomplex>());
+
+  const int nmodes = cmplx.size();
+  for(int jz=1;jz<nmodes;jz++) {
+    cmplx[jz] *= phs[jz];
+  }
+
+  irfft(cmplx.begin(), mesh.LocalNz, out); // Reverse FFT
+}
+
+//Old approach retained so we can still specify a general zShift
+const Field3D ShiftToFieldAligned::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();
+  invalidateGuards(result); // Won't set x-guard cells, so allow checking to throw exception if they are used.
+
+  // We only use methods in ShiftToFieldAligned 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;
+}
+
+void ShiftToFieldAligned::shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutReal *out) {
+  int nmodes = len/2 + 1;
+
+  // Complex array used for FFTs
+  cmplxLoc = Array<dcomplex>(nmodes);
+
+  // Take forward FFT
+  rfft(in, len, cmplxLoc.begin());
+
+  // Apply phase shift
+  BoutReal zlength = mesh.getCoordinates()->zlength();
+  for(int jz=1;jz<nmodes;jz++) {
+    BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+    cmplxLoc[jz] *= dcomplex(cos(kwave*zangle) , -sin(kwave*zangle));
+  }
+
+  irfft(cmplxLoc.begin(), len, out); // Reverse FFT
+}
diff --git a/tests/integrated/test-yupdown/runtest b/tests/integrated/test-yupdown/runtest
index c1da2a7501..99f9fe1bef 100755
--- a/tests/integrated/test-yupdown/runtest
+++ b/tests/integrated/test-yupdown/runtest
@@ -1,7 +1,5 @@
 #!/usr/bin/env python3
 
-#requires: all_tests
-
 from boututils.run_wrapper import shell, shell_safe, launch_safe, getmpirun
 from boutdata.collect import collect
 from sys import stdout, exit
@@ -18,18 +16,18 @@ s, out = launch_safe("./test_yupdown", runcmd=MPIRUN, nproc=1, pipe=True, verbos
 with open("run.log", "w") as f:
   f.write(out)
 
-vars = [ ("ddy", "ddy2") ]
+vars = [ ("ddy", "ddy_check"), ("ddy_tofieldaligned", "ddy_check") ]
 for v1, v2 in vars:
   stdout.write("Testing %s and %s ... " % (v1, v2) )
   ddy = collect(v1, path="data", xguards=False, yguards=False, info=False)
-  ddy2 = collect(v2, path="data", xguards=False, yguards=False, info=False)
+  ddy_check = collect(v2, path="data", xguards=False, yguards=False, info=False)
 
-  diff = max(abs(ddy - ddy2))
+  diff = max(abs(ddy - ddy_check))
 
   if diff < 1e-8:
-    print("Passed (Max difference %e)" % (diff))
+    print("Passed %s=%s (Max difference %e)" % (v1, v2, diff))
   else:
-    print("Failed (Max difference %e)" % (diff))
+    print("Failed %s=%s (Max difference %e)" % (v1, v2, diff))
     exit(1)
   
 exit(0)
diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx
index e414c0ba88..4e58ff81e0 100644
--- a/tests/integrated/test-yupdown/test_yupdown.cxx
+++ b/tests/integrated/test-yupdown/test_yupdown.cxx
@@ -42,6 +42,8 @@ int main(int argc, char** argv) {
   // Read variable from mesh
   Field3D var;
   mesh->get(var, "var");
+
+  Field3D var_tofieldaligned = copy(var);
   
   // Var starts in orthogonal X-Z coordinates
 
@@ -51,16 +53,19 @@ int main(int argc, char** argv) {
   // Calculate d/dy ysing yup() and ydown() fields
   Field3D ddy = DDY_yud(var);
 
+  // Calculate d/dy using ShiftToFieldAligned
+  Field3D ddy_tofieldaligned = DDY(var_tofieldaligned);
+
   // Change into field-aligned coordinates
   Field3D var_aligned = mesh->toFieldAligned(var);
   
   // var now field aligned
-  Field3D ddy2 = DDY_aligned(var_aligned);
+  Field3D ddy_check = DDY_aligned(var_aligned);
   
   // Shift back to orthogonal X-Z coordinates
-  ddy2 = mesh->fromFieldAligned(ddy2);
+  ddy_check = mesh->fromFieldAligned(ddy_check);
   
-  SAVE_ONCE2(ddy, ddy2);
+  SAVE_ONCE3(ddy, ddy_tofieldaligned, ddy_check);
   dump.write();
 
   BoutFinalise();