From 090f306b6cf92e318b853a1d5715f57f4b2ac137 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 14 May 2018 21:40:25 +0100
Subject: [PATCH 01/21] Remove mixed field-aligned/non-field-aligned
 derivatives

Advective and flux derivatives, which take v and f as inputs, had cases
using yup/ydown fields for one of v and f, but not the other. These do
not make sense as multiplication of a field in field-aligned coordinates
with another in non-field-aligned coordinates is incorrect.

Delete these cases and fall back to converting both v and f to
field-aligned if either does not have yup/ydown fields.
---
 src/mesh/difops.cxx       |  58 +------------
 src/mesh/index_derivs.cxx | 177 ++++----------------------------------
 2 files changed, 19 insertions(+), 216 deletions(-)

diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx
index 7cd73dbaed..8abd65fbf7 100644
--- a/src/mesh/difops.cxx
+++ b/src/mesh/difops.cxx
@@ -338,63 +338,11 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg
       result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p;
     }
   }
-  else if (vUseUpDown) {
-    // Only v has up/down fields
-    // f must shift to field aligned coordinates
-    Field3D f_fa = mesh->toFieldAligned(f);
-
-    vval.mm = nan("");
-    vval.pp = nan("");
-
-    for (const auto &i : result.region(region)) {
-
-      fval.mm = f_fa[i.offset(0, -2, 0)];
-      fval.m = f_fa[i.ym()];
-      fval.c = f_fa[i];
-      fval.p = f_fa[i.yp()];
-      fval.pp = f_fa[i.offset(0, 2, 0)];
-
-      vval.m = v.ydown()[i.ym()];
-      vval.c = v[i];
-      vval.p = v.yup()[i.yp()];
-
-      // Left side
-      result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c;
-      // Right side
-      result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p;
-    }
-  }
-  else if (fUseUpDown) {
-    // Only f has up/down fields
-    // v must shift to field aligned coordinates
-    Field3D v_fa = mesh->toFieldAligned(v);
-
-    stencil vval;
-
-    stencil fval;
-    fval.mm = nan("");
-    fval.pp = nan("");
-
-    for (const auto &i : result.region(region)) {
-
-      fval.m = f.ydown()[i.ym()];
-      fval.c = f[i];
-      fval.p = f.yup()[i.yp()];
-
-      vval.mm = v_fa[i.offset(0,-2,0)];
-      vval.m = v_fa[i.ym()];
-      vval.c = v_fa[i];
-      vval.p = v_fa[i.yp()];
-      vval.pp = v_fa[i.offset(0,2,0)];
-
-      // Left side
-      result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c;
-      // Right side
-      result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p;
-    }
-  }
   else {
     // Both must shift to field aligned
+    // (even if one of v and f has yup/ydown fields, it doesn't make sense to
+    // multiply them with one in field-aligned and one in non-field-aligned
+    // coordinates)
     Field3D v_fa = mesh->toFieldAligned(v);
     Field3D f_fa = mesh->toFieldAligned(f);
 
diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index ace0fe0af8..d869e5fde4 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2252,12 +2252,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       func = lookupFluxFunc(table, method);
     }
 
-    // There are four cases, corresponding to whether or not f and v
-    // have yup, ydown fields.
-
-    // If vUseUpDown is true, field "v" has distinct yup and ydown fields which
-    // will be used to calculate a derivative along
-    // the magnetic field
+    // If *UseUpDown is true, field "*" has distinct yup and ydown fields which
+    // will be used to calculate a derivative along the magnetic field
     bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)));
     bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f)));
 
@@ -2277,74 +2273,6 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         fval.p = f.yup()[i.yp()];
         fval.m = f.ydown()[i.ym()];
 
-        if (diffloc != CELL_DEFAULT) {
-          // Non-centred stencil
-          if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-            // Producing a stencil centred around a lower Y value
-            vval.pp = vval.p;
-            vval.p = vval.c;
-          } else if (vloc == CELL_YLOW) {
-            // Stencil centred around a cell centre
-            vval.mm = vval.m;
-            vval.m = vval.c;
-          }
-          // Shifted in one direction -> shift in another
-          // Could produce warning
-        }
-        result[i] = func(vval, fval);
-      }
-    } else if (vUseUpDown) {
-      // Only v has up/down fields
-      // f must shift to field aligned coordinates
-      Field3D f_fa = this->toFieldAligned(f);
-
-      stencil vval, fval;
-      vval.pp = nan("");
-      vval.mm = nan("");
-      for (const auto &i : result.region(region)) {
-        vval.c = v[i];
-        vval.p = v.yup()[i.yp()];
-        vval.m = v.ydown()[i.ym()];
-        fval.c = f_fa[i];
-        fval.p = f_fa[i.yp()];
-        fval.m = f_fa[i.ym()];
-        fval.pp = f_fa[i.offset(0, 2, 0)];
-        fval.mm = f_fa[i.offset(0, -2, 0)];
-
-        if (diffloc != CELL_DEFAULT) {
-          // Non-centred stencil
-          if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-            // Producing a stencil centred around a lower Y value
-            vval.pp = vval.p;
-            vval.p = vval.c;
-          } else if (vloc == CELL_YLOW) {
-            // Stencil centred around a cell centre
-            vval.mm = vval.m;
-            vval.m = vval.c;
-          }
-          // Shifted in one direction -> shift in another
-          // Could produce warning
-        }
-        result[i] = func(vval, fval);
-      }
-    } else if (fUseUpDown) {
-      // Only f has up/down fields
-      // v must shift to field aligned coordinates
-      Field3D v_fa = this->toFieldAligned(v);
-
-      stencil vval, fval;
-      fval.pp = nan("");
-      fval.mm = nan("");
-      for (const auto &i : result.region(region)) {
-        vval.c = v_fa[i];
-        vval.p = v_fa[i.yp()];
-        vval.m = v_fa[i.ym()];
-        vval.pp = v_fa[i.offset(0, 2, 0)];
-        vval.mm = v_fa[i.offset(0, -2, 0)];
-        fval.c = f[i];
-        fval.p = f.yup()[i.yp()];
-        fval.m = f.ydown()[i.ym()];
-
         if (diffloc != CELL_DEFAULT) {
           // Non-centred stencil
           if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
@@ -2363,6 +2291,9 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       }
     } else {
       // Both must shift to field aligned
+      // (even if one of v and f has yup/ydown fields, it doesn't make sense to
+      // multiply them with one in field-aligned and one in non-field-aligned
+      // coordinates)
       Field3D v_fa = this->toFieldAligned(v);
       Field3D f_fa = this->toFieldAligned(f);
 
@@ -2407,8 +2338,9 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       func = lookupUpwindFunc(table, method);
     }
 
-    if (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))) {
-      // f has yup and ydown fields which are distinct
+    if ((f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))) 
+        && (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)))) {
+      // f and v have yup and ydown fields which are distinct
 
       stencil fs;
       fs.pp = nan("");
@@ -2428,6 +2360,9 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
     } else {
       // Not using yup/ydown fields, so first transform to field-aligned coordinates
+      // (even if one of v and f has yup/ydown fields, it doesn't make sense to
+      // multiply them with one in field-aligned and one in non-field-aligned
+      // coordinates)
 
       Field3D f_fa = this->toFieldAligned(f);
       Field3D v_fa = this->toFieldAligned(v);
@@ -3015,12 +2950,8 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
   Field3D result(this);
   result.allocate(); // Make sure data allocated
 
-  // There are four cases, corresponding to whether or not f and v
-  // have yup, ydown fields.
-
-  // If vUseUpDown is true, field "v" has distinct yup and ydown fields which
-  // will be used to calculate a derivative along
-  // the magnetic field
+  // If *UseUpDown is true, field "*" has distinct yup and ydown fields which
+  // will be used to calculate a derivative along the magnetic field
   bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)));
   bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f)));
 
@@ -3058,87 +2989,11 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       result[i] = func(vval, fval);
     }
   }
-  else if (vUseUpDown) {
-    // Only v has up/down fields
-    // f must shift to field aligned coordinates
-    Field3D f_fa = this->toFieldAligned(f);
-
-    stencil vval;
-    vval.mm = nan("");
-    vval.pp = nan("");
-
-    stencil fval;
-    for (const auto &i : result.region(region)) {
-
-      fval.mm = f_fa[i.offset(0, -2, 0)];
-      fval.m = f_fa[i.ym()];
-      fval.c = f_fa[i];
-      fval.p = f_fa[i.yp()];
-      fval.pp = f_fa[i.offset(0, 2, 0)];
-
-      vval.m = v.ydown()[i.ym()];
-      vval.c = v[i];
-      vval.p = v.yup()[i.yp()];
-
-      if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) {
-        // Non-centred stencil
-        if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-          // Producing a stencil centred around a lower Y value
-          vval.pp = vval.p;
-          vval.p  = vval.c;
-        }else if(vloc == CELL_YLOW) {
-          // Stencil centred around a cell centre
-          vval.mm = vval.m;
-          vval.m  = vval.c;
-        }
-        // Shifted in one direction -> shift in another
-        // Could produce warning
-      }
-      result[i] = func(vval, fval);
-    }
-  }
-  else if (fUseUpDown) {
-    // Only f has up/down fields
-    // v must shift to field aligned coordinates
-    Field3D v_fa = this->toFieldAligned(v);
-
-    stencil vval;
-
-    stencil fval;
-    fval.pp = nan("");
-    fval.mm = nan("");
-
-    for (const auto &i : result.region(region)) {
-
-      fval.m = f.ydown()[i.ym()];
-      fval.c = f[i];
-      fval.p = f.yup()[i.yp()];
-
-      vval.mm = v_fa[i.offset(0,-2,0)];
-      vval.m = v_fa[i.ym()];
-      vval.c = v_fa[i];
-      vval.p = v_fa[i.yp()];
-      vval.pp = v_fa[i.offset(0,2,0)];
-
-      if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) {
-        // Non-centred stencil
-        if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-          // Producing a stencil centred around a lower Y value
-          vval.pp = vval.p;
-          vval.p  = vval.c;
-        }else if(vloc == CELL_YLOW) {
-          // Stencil centred around a cell centre
-          vval.mm = vval.m;
-          vval.m  = vval.c;
-        }
-        // Shifted in one direction -> shift in another
-        // Could produce warning
-      }
-      result[i] = func(vval, fval);
-    }
-  }
   else {
     // Both must shift to field aligned
+    // (even if one of v and f has yup/ydown fields, it doesn't make sense to
+    // multiply them with one in field-aligned and one in non-field-aligned
+    // coordinates)
     Field3D v_fa = this->toFieldAligned(v);
     Field3D f_fa = this->toFieldAligned(f);
 

From 04d9984ea2b8a6cf8cef7ed0e8b762c9ea42918c Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 17 Jul 2018 12:08:49 +0100
Subject: [PATCH 02/21] Remove result_fa from interp_to

Remove the intermediate variable result_fa in branch of interp_to that
transforms to field-aligned variables. Just store the intermediate
result in 'result' instead, to save memory.
---
 src/mesh/interpolation.cxx | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx
index 131e1c4c21..2d47b668b8 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -146,8 +146,6 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
           // coordinates
 
           Field3D var_fa = fieldmesh->toFieldAligned(var);
-          Field3D result_fa;
-          result_fa.allocate();
           if (fieldmesh->ystart > 1) {
 
             // More than one guard cell, so set pp and mm values
@@ -170,7 +168,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 {
             // Only one guard cell, so no pp or mm values
@@ -195,11 +193,11 @@ 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);
         }
         break;
       }

From b8e6336ac9bc8217aaa7c02dea9f10df92bd5288 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 14 May 2018 21:31:39 +0100
Subject: [PATCH 03/21] Add missing fromFieldAligned calls

In several y-derivatives where we have to convert the input to
field-aligned coordinates, the result was not transformed back to
the original coordinates.
---
 src/mesh/difops.cxx       | 12 ++++++++----
 src/mesh/index_derivs.cxx |  4 ++++
 2 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx
index 8abd65fbf7..b63a57ef07 100644
--- a/src/mesh/difops.cxx
+++ b/src/mesh/difops.cxx
@@ -365,6 +365,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg
       // Right side
       result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p;
     }
+
+    result = mesh->fromFieldAligned(result);
   }
 
   result.setLocation(CELL_CENTRE);
@@ -376,10 +378,11 @@ const Field3D Grad_par_LtoC(const Field3D &var) {
     ASSERT1(var.getLocation() == CELL_YLOW);
   } 
 
-  Field3D result(var.getMesh());
+  Mesh* fieldmesh = var.getMesh();
+  Field3D result(fieldmesh);
   result.allocate();
 
-  Coordinates *metric = var.getMesh()->coordinates();
+  Coordinates *metric = fieldmesh->coordinates();
 
   if (var.hasYupYdown()) {
     for (auto &i : result.region(RGN_NOBNDRY)) {
@@ -388,12 +391,13 @@ const Field3D Grad_par_LtoC(const Field3D &var) {
   } else {
     // No yup/ydown field, so transform to field aligned
 
-    Field3D var_fa = var.getMesh()->toFieldAligned(var);
+    Field3D var_fa = fieldmesh->toFieldAligned(var);
 
     for(auto &i : result.region(RGN_NOBNDRY)) {
       result[i] = (var_fa[i.yp()] - var_fa[i]) / (metric->dy[i]*sqrt(metric->g_22[i]));
     }
-    result = var.getMesh()->fromFieldAligned(result);
+
+    result = fieldmesh->fromFieldAligned(result);
   }
 
   result.setLocation(CELL_CENTRE);
diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index d869e5fde4..5f07e97b81 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2326,6 +2326,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         }
         result[i] = func(vval, fval);
       }
+
+      result = this->fromFieldAligned(result);
     }
   } else {
     // Non-staggered case
@@ -3029,6 +3031,8 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       }
       result[i] = func(vval, fval);
     }
+
+    result = this->fromFieldAligned(result);
   }
 
   result.setLocation(diffloc);

From 9bbed30493dd282e2d3d024418b3199be453396e Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 2 Jul 2018 20:40:55 +0100
Subject: [PATCH 04/21] Remove unnecessary check for v.hasYupYdown()

In the unstaggered case, only the value of v is used, so yup/ydown
fields for v are not needed.
---
 src/mesh/index_derivs.cxx | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index 5f07e97b81..15b5ce3db2 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2340,9 +2340,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       func = lookupUpwindFunc(table, method);
     }
 
-    if ((f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))) 
-        && (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)))) {
-      // f and v have yup and ydown fields which are distinct
+    if (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))) {
+      // f has yup and ydown fields which are distinct
 
       stencil fs;
       fs.pp = nan("");

From ee2731b89f833d4696f5316214d581153a16db67 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 3 Jul 2018 09:58:48 +0100
Subject: [PATCH 05/21] Fix f->f_fa

Was previously a typo that resulted in using non-field-aligned f in
VDDY, which is incorrect.
---
 src/mesh/index_derivs.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index 15b5ce3db2..7eacfc99f9 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2304,7 +2304,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         vval.m = v_fa[i.ym()];
         vval.pp = v_fa[i.offset(0, 2, 0)];
         vval.mm = v_fa[i.offset(0, -2, 0)];
-        fval.c = f[i];
+        fval.c = f_fa[i];
         fval.p = f_fa[i.yp()];
         fval.m = f_fa[i.ym()];
         fval.pp = f_fa[i.offset(0, 2, 0)];

From ee9231ac1f2a50021222217b9a9471cc2a88e166 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 3 Jul 2018 11:40:57 +0100
Subject: [PATCH 06/21] Correct location setting in
 applyXdiff/applyYdiff/applyZdiff

Previously these were incorrectly setting the location of their result
to the location of the input variable. This commit changes this to use
the specified output location. This probably did not cause bugs as the
location was overridden in DDX/DDY/DDZ, etc.
---
 src/mesh/index_derivs.cxx | 17 ++++-------------
 1 file changed, 4 insertions(+), 13 deletions(-)

diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index 7eacfc99f9..4f5dca8d2f 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -677,8 +677,6 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func,
     return Field2D(0., this);
   }
 
-  CELL_LOC diffloc = var.getLocation();
-
   Field2D result(this);
   result.allocate(); // Make sure data allocated
 
@@ -765,7 +763,7 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func,
     }
   }
 
-  result.setLocation(diffloc);
+  result.setLocation(loc);
 
 #if CHECK > 0
   // Mark boundaries as invalid
@@ -786,8 +784,6 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func,
     return Field3D(0., this);
   }
 
-  CELL_LOC diffloc = var.getLocation();
-
   Field3D result(this);
   result.allocate(); // Make sure data allocated
 
@@ -931,7 +927,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L
     }
   }
 
-  result.setLocation(diffloc);
+  result.setLocation(loc);
 
 #if CHECK > 0
   // Mark boundaries as invalid
@@ -951,8 +947,6 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
     return Field3D(0., this);
   }
 
-  CELL_LOC diffloc = var.getLocation();
-
   Field3D result(this);
   result.allocate(); // Make sure data allocated
 
@@ -1122,9 +1116,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
     return Field3D(0., this);
   }
 
-
-  CELL_LOC diffloc = var.getLocation();
-
   if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) {
     // Staggered differencing
     throw BoutException("No one used this before. And no one implemented it.");
@@ -1147,7 +1138,7 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
     result[i] = func(s);
   }
 
-  result.setLocation(diffloc);
+  result.setLocation(loc);
 
   return result;
 }
@@ -3034,7 +3025,7 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
     result = this->fromFieldAligned(result);
   }
 
-  result.setLocation(diffloc);
+  result.setLocation(loc);
 
 #if CHECK > 0
   // Mark boundaries as invalid

From 4784ef71c7377e996f15143e514eba67458774a7 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 3 Jul 2018 12:16:41 +0100
Subject: [PATCH 07/21] Set location of result as soon as it is declared in
 index_derivs.cxx

The location of the result needs to be set before shifting it from
field-aligned coordinates, if this is necessary. So it is safer to set
the location to diffloc (or outloc where diffloc is not yet set) as soon
as 'Field3D result' is declared.
There are a couple of  exceptions where the location is set after
calling apply*diff (which should not be necessary, but does not hurt,
except that it may hide errors setting location in apply*diff) where the
location of 'result' would not be used anyway earlier in the method.
---
 src/mesh/index_derivs.cxx | 60 +++++++++++++--------------------------
 1 file changed, 20 insertions(+), 40 deletions(-)

diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index 4f5dca8d2f..a1fd2d0adc 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -679,6 +679,7 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func,
 
   Field2D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(loc);
 
   if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) {
     // Staggered differencing
@@ -763,8 +764,6 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func,
     }
   }
 
-  result.setLocation(loc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -786,6 +785,7 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func,
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(loc);
 
   if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) {
     // Staggered differencing
@@ -870,8 +870,6 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func,
     }
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -896,6 +894,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L
 
   Field2D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(diffloc);
 
   if (this->ystart > 1) {
     // More than one guard cell, so set pp and mm values
@@ -927,8 +926,6 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L
     }
   }
 
-  result.setLocation(loc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_yup = result.bndry_ydown = false;
@@ -949,6 +946,7 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(loc);
 
   if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) {
     // Field "var" has distinct yup and ydown fields which
@@ -1094,8 +1092,6 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
     result = this->fromFieldAligned(result);
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -1123,6 +1119,7 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(loc);
 
   // Check that the input variable has data
   ASSERT1(var.isAllocated());
@@ -1138,8 +1135,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
     result[i] = func(s);
   }
 
-  result.setLocation(loc);
-
   return result;
 }
 
@@ -1157,8 +1152,6 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth
   CELL_LOC inloc = f.getLocation(); // Input location
   CELL_LOC diffloc = inloc;         // Location of differential result
 
-  Field3D result(this);
-
   if (this->StaggerGrids && (outloc == CELL_DEFAULT)) {
     // Take care of CELL_DEFAULT case
     outloc = diffloc; // No shift (i.e. same as no stagger case)
@@ -1190,7 +1183,7 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth
       throw BoutException("Cannot use FFT for X derivatives");
   }
 
-  result = applyXdiff(f, func, diffloc, region);
+  Field3D result = applyXdiff(f, func, diffloc, region);
 
   result.setLocation(diffloc); // Set the result location
 
@@ -1214,8 +1207,6 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc,
   CELL_LOC inloc = f.getLocation(); // Input location
   CELL_LOC diffloc = inloc;         // Location of differential result
 
-  Field3D result(this);
-
   if (this->StaggerGrids && (outloc == CELL_DEFAULT)) {
     // Take care of CELL_DEFAULT case
     outloc = diffloc; // No shift (i.e. same as no stagger case)
@@ -1246,7 +1237,7 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc,
       throw BoutException("Cannot use FFT for Y derivatives");
   }
 
-  result = applyYdiff(f, func, diffloc, region);
+  Field3D result = applyYdiff(f, func, diffloc, region);
 
   result.setLocation(diffloc); // Set the result location
 
@@ -1417,8 +1408,6 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc,
 
   ASSERT1(this == f.getMesh());
 
-  Field3D result(this);
-
   if (StaggerGrids && (outloc == CELL_DEFAULT)) {
     // Take care of CELL_DEFAULT case
     outloc = diffloc; // No shift (i.e. same as no stagger case)
@@ -1450,7 +1439,7 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc,
       throw BoutException("Cannot use FFT for X derivatives");
   }
 
-  result = applyXdiff(f, func, diffloc, region);
+  Field3D result = applyXdiff(f, func, diffloc, region);
 
   result.setLocation(diffloc);
 
@@ -1496,8 +1485,6 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc,
 
   ASSERT1(this == f.getMesh());
 
-  Field3D result(this);
-
   if (StaggerGrids && (outloc == CELL_DEFAULT)) {
     // Take care of CELL_DEFAULT case
     outloc = diffloc; // No shift (i.e. same as no stagger case)
@@ -1529,7 +1516,7 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc,
       throw BoutException("Cannot use FFT for Y derivatives");
   }
 
-  result = applyYdiff(f, func, diffloc, region);
+  Field3D result = applyYdiff(f, func, diffloc, region);
 
   result.setLocation(diffloc);
 
@@ -1752,6 +1739,7 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo
 
   Field2D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   if (this->xstart > 1) {
     // Two or more guard cells
@@ -1791,8 +1779,6 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo
   result.bndry_xin = result.bndry_xout = false;
 #endif
 
-  result.setLocation(diffloc);
-
   return result;
 }
 
@@ -1807,6 +1793,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   CELL_LOC vloc = v.getLocation();
   CELL_LOC inloc = f.getLocation(); // Input location
@@ -1980,8 +1967,6 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo
     }
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -2003,6 +1988,7 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo
 
   Field2D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   CELL_LOC vloc = v.getLocation();
   CELL_LOC inloc = f.getLocation(); // Input location
@@ -2015,7 +2001,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo
 
   if (this->LocalNy == 1){
     result=0;
-    result.setLocation(outloc);
     return result;
   }
 
@@ -2178,8 +2163,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo
     }
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -2198,6 +2181,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   CELL_LOC vloc = v.getLocation();
   CELL_LOC inloc = f.getLocation(); // Input location
@@ -2210,7 +2194,6 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   if (this->LocalNy == 1){
     result=0;
-    result.setLocation(outloc);
     return result;
   }
 
@@ -2411,6 +2394,7 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   CELL_LOC vloc = v.getLocation();
   CELL_LOC inloc = f.getLocation(); // Input location
@@ -2502,8 +2486,6 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo
     }
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -2520,6 +2502,8 @@ const Field2D Mesh::indexFDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo
                               DIFF_METHOD method, REGION region) {
   TRACE("Mesh::::indexFDDX(Field2D, Field2D)");
 
+  CELL_LOC diffloc = f.getLocation();
+
   if ((method == DIFF_SPLIT) || ((method == DIFF_DEFAULT) && (fFDDX == nullptr))) {
     // Split into an upwind and a central differencing part
     // d/dx(v*f) = v*d/dx(f) + f*d/dx(v)
@@ -2651,6 +2635,7 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(diffloc);
 
   if (this->xstart > 1) {
     // Two or more guard cells
@@ -2787,8 +2772,6 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo
     throw BoutException("Error: Derivatives in X requires at least one guard cell");
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -2822,6 +2805,7 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo
 
   Field2D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(outloc);
 
   if (StaggerGrids &&
       ((v.getLocation() != CELL_CENTRE) || (f.getLocation() != CELL_CENTRE))) {
@@ -2877,8 +2861,6 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo
     throw BoutException("Error: Derivatives in Y requires at least one guard cell");
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = false;
@@ -2941,6 +2923,7 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(diffloc);
 
   // If *UseUpDown is true, field "*" has distinct yup and ydown fields which
   // will be used to calculate a derivative along the magnetic field
@@ -3025,8 +3008,6 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
     result = this->fromFieldAligned(result);
   }
 
-  result.setLocation(loc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
@@ -3087,6 +3068,7 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo
 
   Field3D result(this);
   result.allocate(); // Make sure data allocated
+  result.setLocation(diffloc);
 
   stencil vval, fval;
   for (const auto &i : result.region(region)) {
@@ -3123,8 +3105,6 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo
     result[i] = func(vval, fval);
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;

From ac726d0a0c3313fcc1c9ef440d40eb77ef3f0697 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 23 May 2018 18:00:11 +0100
Subject: [PATCH 08/21] Make a method to clean up deleting yup/ydown and
 field_fa members

---
 include/field3d.hxx              |  5 ++++
 src/field/field3d.cxx            | 32 ++++++++++++++---------
 src/field/gen_fieldops.jinja     |  5 ++++
 src/field/generated_fieldops.cxx | 44 ++++++++++++++++++++++++++++++++
 4 files changed, 74 insertions(+), 12 deletions(-)

diff --git a/include/field3d.hxx b/include/field3d.hxx
index e24ea582f3..028eccc364 100644
--- a/include/field3d.hxx
+++ b/include/field3d.hxx
@@ -233,6 +233,11 @@ class Field3D : public Field, public FieldData {
    * Ensure that yup and ydown refer to this field
    */
   void mergeYupYdown();
+
+  /*!
+   * Clear all yup/ydown fields and field_fa
+   */
+  void clearYupYdown();
   
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index e6f72dc12b..4b7b828f25 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -149,12 +149,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;
+
+  clearYupYdown();
 }
 
 void Field3D::allocate() {
@@ -198,16 +194,23 @@ void Field3D::mergeYupYdown() {
   if(yup_field == this && ydown_field == this)
     return;
 
-  if(yup_field != nullptr){
+  clearYupYdown();
+
+  yup_field = this;
+  ydown_field = this;
+}
+
+void Field3D::clearYupYdown() {
+  // Delete auxiliary fields if they have been set
+  if (yup_field != nullptr && yup_field != this) {
     delete yup_field;
   }
+  yup_field = nullptr;
 
-  if(ydown_field != nullptr) {
+  if (ydown_field != nullptr && ydown_field != this) {
     delete ydown_field;
   }
-
-  yup_field = this;
-  ydown_field = this;
+  ydown_field = nullptr;
 }
 
 Field3D& Field3D::ynext(int dir) {
@@ -357,6 +360,7 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
   
   location = rhs.location;
   
+  clearYupYdown();
   return *this;
 }
 
@@ -377,7 +381,9 @@ Field3D & Field3D::operator=(const Field2D &rhs) {
   
   /// Only 3D fields have locations for now
   //location = CELL_CENTRE;
-  
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -408,6 +414,8 @@ Field3D & Field3D::operator=(const BoutReal val) {
   //location = CELL_CENTRE;
   // DON'T RE-SET LOCATION
 
+  clearYupYdown();
+
   return *this;
 }
 
diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja
index 31ae591235..93a7890d50 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") %}
+    clearYupYdown();
+  {% endif %}
+
   return *this;
 }
 {% endif %}
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 8a8b052b92..b75f25a8a6 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -61,6 +61,9 @@ Field3D &Field3D::operator*=(const Field3D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -119,6 +122,9 @@ Field3D &Field3D::operator/=(const Field3D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -177,6 +183,9 @@ Field3D &Field3D::operator+=(const Field3D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -235,6 +244,9 @@ Field3D &Field3D::operator-=(const Field3D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -283,6 +295,9 @@ Field3D &Field3D::operator*=(const Field2D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -332,6 +347,9 @@ Field3D &Field3D::operator/=(const Field2D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -380,6 +398,9 @@ Field3D &Field3D::operator+=(const Field2D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -428,6 +449,9 @@ Field3D &Field3D::operator-=(const Field2D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -466,6 +490,9 @@ Field3D &Field3D::operator*=(const BoutReal rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -504,6 +531,9 @@ Field3D &Field3D::operator/=(const BoutReal rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -542,6 +572,9 @@ Field3D &Field3D::operator+=(const BoutReal rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -580,6 +613,9 @@ Field3D &Field3D::operator-=(const BoutReal rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
+  clearYupYdown();
+
   return *this;
 }
 
@@ -717,6 +753,7 @@ Field2D &Field2D::operator*=(const Field2D &rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
   return *this;
 }
 
@@ -758,6 +795,7 @@ Field2D &Field2D::operator/=(const Field2D &rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
   return *this;
 }
 
@@ -799,6 +837,7 @@ Field2D &Field2D::operator+=(const Field2D &rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
   return *this;
 }
 
@@ -840,6 +879,7 @@ Field2D &Field2D::operator-=(const Field2D &rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
   return *this;
 }
 
@@ -876,6 +916,7 @@ Field2D &Field2D::operator*=(const BoutReal rhs) {
   } else {
     (*this) = (*this) * rhs;
   }
+
   return *this;
 }
 
@@ -912,6 +953,7 @@ Field2D &Field2D::operator/=(const BoutReal rhs) {
   } else {
     (*this) = (*this) / rhs;
   }
+
   return *this;
 }
 
@@ -948,6 +990,7 @@ Field2D &Field2D::operator+=(const BoutReal rhs) {
   } else {
     (*this) = (*this) + rhs;
   }
+
   return *this;
 }
 
@@ -984,6 +1027,7 @@ Field2D &Field2D::operator-=(const BoutReal rhs) {
   } else {
     (*this) = (*this) - rhs;
   }
+
   return *this;
 }
 

From 8cc925ed3ffae73ecbec41525711ee9d7088912e Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 18 Jul 2018 11:10:28 +0100
Subject: [PATCH 09/21] Add region arguments for
 toFieldAligned/fromFieldAligned

Add region arguments (where it is valid to give RGN_NOX or RGN_NOBNDRY)
so functions like toFieldAligned/fromFieldAligned can be used both when
y-guard cells have been set and when they have not. e.g. when using
fromFieldAligned on the results of interpolation or derivatives, y-guard
cells cannot be included, but when using toFieldAligned on the inputs to
interpolation or derivatives, y-guard cells must be included.
---
 include/bout/mesh.hxx               |  8 +++---
 include/bout/paralleltransform.hxx  | 18 +++++++-------
 src/mesh/index_derivs.cxx           |  2 +-
 src/mesh/interpolation.cxx          |  4 ++-
 src/mesh/parallel/fci.hxx           |  4 +--
 src/mesh/parallel/shiftedmetric.cxx | 38 +++++++++++++++++------------
 6 files changed, 41 insertions(+), 33 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 139175051c..de1c009852 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -603,12 +603,12 @@ class Mesh {
   typedef BoutReal (*flux_func)(stencil&, stencil &);
 
   /// 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 afd296844d..6679e8d8bb 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;
 };
@@ -65,7 +65,7 @@ public:
    * 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 +73,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;
   }
 
@@ -108,13 +108,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;
@@ -143,7 +143,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
@@ -152,7 +152,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
@@ -163,7 +163,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/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index a1fd2d0adc..07e4084e08 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -2301,7 +2301,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         result[i] = func(vval, fval);
       }
 
-      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 2d47b668b8..b45f42d41c 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -197,7 +197,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
             }
           }
           
-          result = fieldmesh->fromFieldAligned(result);
+          result = fieldmesh->fromFieldAligned(result, RGN_NOBNDRY);
         }
         break;
       }
@@ -230,6 +230,8 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
       }
       };
 
+      invalidateGuards(result); // Fill guard cells with NaN so we can check they are not used when unset.
+
       if ((dir != CELL_ZLOW) && (region != RGN_NOBNDRY)) {
         fieldmesh->communicate(result);
       }
diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx
index 05bd0093d6..017375b668 100644
--- a/src/mesh/parallel/fci.hxx
+++ b/src/mesh/parallel/fci.hxx
@@ -78,11 +78,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 13067fe07c..607dcd593f 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -120,30 +120,29 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
  * 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
   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));
-    }
+
+  for(auto i : f.region2D(region)) {
+    shiftZ(f(i.x,i.y), phs[i.x][i.y], result(i.x,i.y));
   }
   
   return result;
@@ -168,18 +167,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));
-    }
+  invalidateGuards(result); // Won't set x-guard cells, so allow checking to throw exception if they are used.
+
+  // 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(auto i : f.region2D(region)) {
+    shiftZ(f(i.x, i.y), mesh.LocalNz, zangle(i.x,i.y), result(i.x, i.y));
   }
   
   return result;

From ee09524ac506791eb1ab2cf095c997f5b9effad9 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 22 May 2018 13:18:08 +0100
Subject: [PATCH 10/21] Write zShift to output file when using ShiftedMetric

---
 include/bout/paralleltransform.hxx  | 10 ++++++++++
 src/mesh/impls/bout/boutmesh.cxx    |  1 +
 src/mesh/parallel/shiftedmetric.cxx |  4 ++++
 3 files changed, 15 insertions(+)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index 6679e8d8bb..db626c7d04 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -6,6 +6,7 @@
 #ifndef __PARALLELTRANSFORM_H__
 #define __PARALLELTRANSFORM_H__
 
+#include <datafile.hxx>
 #include <field3d.hxx>
 #include <boutexception.hxx>
 #include <dcomplex.hxx>
@@ -45,6 +46,9 @@ public:
   virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
   virtual bool canToFromFieldAligned() = 0;
+
+  /// Write out ParallelTransform variables to file
+  virtual void outputVars(Datafile &UNUSED(file)) {};
 };
 
 
@@ -80,6 +84,9 @@ public:
   bool canToFromFieldAligned() override{
     return true;
   }
+
+  /// Write out ParallelTransform variables to file
+  virtual void outputVars(Datafile &UNUSED(file)) {};
 };
 
 /*!
@@ -183,6 +190,9 @@ private:
    * @param[out] out  A 1D array of length mesh.LocalNz, already allocated
    */
   void shiftZ(const BoutReal *in, const std::vector<dcomplex> &phs, BoutReal *out);
+
+  /// Write out ParallelTransform variables to file
+  void outputVars(Datafile &file);
 };
 
 
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 11b855c937..19114cdd1d 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -2548,4 +2548,5 @@ void BoutMesh::outputVars(Datafile &file) {
   file.add(jyseps2_2, "jyseps2_2", 0);
 
   coordinates()->outputVars(file);
+  getParallelTransform().outputVars(file);
 }
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 607dcd593f..614e3376dd 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -209,3 +209,7 @@ void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutRe
 
   irfft(&cmplxLoc[0], len, out); // Reverse FFT
 }
+
+void ShiftedMetric::outputVars(Datafile &file) {
+  file.add(zShift, "zShift", 0);
+}

From de63f2a113c7cbcd85b960ae6070a3ac895de461 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Thu, 3 May 2018 09:52:08 +0100
Subject: [PATCH 11/21] Allow 2 yup/ydown fields in ShiftedMetric

Replace yup_field and ydown_field with yup1_field, ydown1_field,
yup2_field and ydown2_field. If MYG>1 (tested using fieldmesh->ystart>1)
then calculate yup2_field and ydown2_field as well.

This allows derivatives which require 5-point stencils when using
shifted metric.
---
 include/bout/paralleltransform.hxx  |   6 +-
 include/field3d.hxx                 |  47 ++---
 src/field/field3d.cxx               | 176 +++++++++++++----
 src/mesh/coordinates.cxx            |   5 +
 src/mesh/index_derivs.cxx           | 292 ++++++++++++++++++++--------
 src/mesh/interpolation.cxx          |  62 ++++--
 src/mesh/parallel/shiftedmetric.cxx |  74 +++++--
 7 files changed, 469 insertions(+), 193 deletions(-)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index db626c7d04..2010d327c0 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -143,8 +143,10 @@ private:
   arr3Dvec toAlignedPhs; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates
   arr3Dvec fromAlignedPhs; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates
 
-  arr3Dvec yupPhs; ///< Cache of phase shifts for calculating yup fields
-  arr3Dvec ydownPhs; ///< Cache of phase shifts for calculating ydown fields
+  arr3Dvec yupPhs1; ///< Cache of phase shifts for calculating yup1 fields
+  arr3Dvec ydownPhs1; ///< Cache of phase shifts for calculating ydown1 fields
+  arr3Dvec yupPhs2; ///< Cache of phase shifts for calculating yup2 fields
+  arr3Dvec ydownPhs2; ///< Cache of phase shifts for calculating ydown2 fields
 
   /*!
    * Shift a 2D field in Z. 
diff --git a/include/field3d.hxx b/include/field3d.hxx
index 87fa570695..7ea5d5e6b5 100644
--- a/include/field3d.hxx
+++ b/include/field3d.hxx
@@ -151,21 +151,21 @@ template <typename F> class Flexible;
 
       Field3D f(0.0); // f allocated, set to zero
 
-      f.yup() // error; f.yup not allocated
+      f.yup(1) // error; f.yup not allocated
 
-      f.mergeYupYdown(); // f.yup() and f.ydown() now point to f
-      f.yup()(0,1,0)  // ok, gives value of f at (0,1,0)
+      f.mergeYupYdown(); // f.yup(i) and f.ydown(i) now point to f
+      f.yup(1)(0,1,0)  // ok, gives value of f at (0,1,0)
 
       To have separate fields for yup and ydown, first call
 
-      f.splitYupYdown(); // f.yup() and f.ydown() separate
+      f.splitYupYdown(); // f.yup(i) and f.ydown(i) separate
 
-      f.yup(); // ok
-      f.yup()(0,1,0) // error; f.yup not allocated
+      f.yup(1); // ok
+      f.yup(1)(0,1,0) // error; f.yup not allocated
 
-      f.yup() = 1.0; // Set f.yup() field to 1.0
+      f.yup(1) = 1.0; // Set f.yup() field to 1.0
 
-      f.yup()(0,1,0) // ok
+      f.yup(1)(0,1,0) // ok
 
  */
 class Field3D : public Field, public FieldData {
@@ -242,33 +242,22 @@ class Field3D : public Field, public FieldData {
   
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
-    return (yup_field != nullptr) && (ydown_field != nullptr);
+    return (yup1_field != nullptr) && (ydown1_field != nullptr);
   }
 
   /// Return reference to yup field
-  Field3D& yup() { 
-    ASSERT2(yup_field != nullptr); // Check for communicate
-    return *yup_field; 
-  }
-  /// Return const reference to yup field
-  const Field3D& yup() const { 
-    ASSERT2(yup_field != nullptr);
-    return *yup_field; 
-  }
+  Field3D& yup(const int i = 1);
   
+  /// Return const reference to yup field
+  const Field3D& yup(const int i = 1) const;
+
   /// Return reference to ydown field
-  Field3D& ydown() { 
-    ASSERT2(ydown_field != nullptr);
-    return *ydown_field;
-  }
+  Field3D& ydown(const int i = 1);
   
   /// Return const reference to ydown field
-  const Field3D& ydown() const { 
-    ASSERT2(ydown_field != nullptr);
-    return *ydown_field; 
-  }
+  const Field3D& ydown(const int i = 1) const;
 
-  /// Return yup if dir=+1, and ydown if dir=-1
+  /// Return yup(dir) if dir>0, and ydown(-dir) if dir<0
   Field3D& ynext(int dir);
   const Field3D& ynext(int dir) const;
 
@@ -533,8 +522,8 @@ private:
   
   Field3D *deriv; ///< Time derivative (may be NULL)
 
-  /// Pointers to fields containing values along Y
-  Field3D *yup_field, *ydown_field;
+  /// Arrays of pointers to fields containing values along Y
+  Field3D *yup1_field, *ydown1_field, *yup2_field, *ydown2_field;
 };
 
 // Non-member overloaded operators
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index cbed1eaa79..6435316388 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -45,8 +45,8 @@
 
 /// Constructor
 Field3D::Field3D(Mesh *localmesh)
-    : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr),
-      ydown_field(nullptr) {
+    : Field(localmesh), background(nullptr), deriv(nullptr), yup1_field(nullptr),
+      ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) {
 #ifdef TRACK
   name = "<F3D>";
 #endif
@@ -74,7 +74,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), yup1_field(nullptr),
+      ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -101,8 +102,8 @@ Field3D::Field3D(const Field3D &f)
 }
 
 Field3D::Field3D(const Field2D &f)
-    : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup_field(nullptr),
-      ydown_field(nullptr) {
+    : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup1_field(nullptr),
+      ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) {
 
   TRACE("Field3D: Copy constructor from Field2D");
 
@@ -119,8 +120,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) {
+    : Field(localmesh), background(nullptr), deriv(nullptr), yup1_field(nullptr),
+      ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) {
 
   TRACE("Field3D: Copy constructor from value");
 
@@ -141,10 +142,14 @@ Field3D::~Field3D() {
     // The ddt of the yup/ydown_fields point to the same place as ddt.yup_field
     // only delete once
     // Also need to check that separate yup_field exists
-    if ((yup_field != this) && (yup_field != nullptr))
-      yup_field->deriv = nullptr;
-    if ((ydown_field != this) && (ydown_field != nullptr))
-      ydown_field->deriv = nullptr;
+    if ((yup1_field != this) && (yup1_field != nullptr))
+      yup1_field->deriv = nullptr;
+    if ((ydown1_field != this) && (ydown1_field != nullptr))
+      ydown1_field->deriv = nullptr;
+    if ((yup2_field != this) && (yup2_field != nullptr))
+      yup2_field->deriv = nullptr;
+    if ((ydown2_field != this) && (ydown2_field != nullptr))
+      ydown2_field->deriv = nullptr;
 
     // Now delete them as part of the deriv vector
     delete deriv;
@@ -180,58 +185,150 @@ Field3D* Field3D::timeDeriv() {
 void Field3D::splitYupYdown() {
   TRACE("Field3D::splitYupYdown");
   
-  if((yup_field != this) && (yup_field != nullptr))
+  if((yup1_field != this) && (yup1_field != nullptr))
     return;
 
-  // yup_field and ydown_field null
-  yup_field = new Field3D(fieldmesh);
-  ydown_field = new Field3D(fieldmesh);
+  // yup_array and ydown_array null
+  yup1_field = new Field3D(fieldmesh);
+  ydown1_field = new Field3D(fieldmesh);
+  if (fieldmesh->ystart>1) {
+    yup2_field = new Field3D(fieldmesh);
+    ydown2_field = new Field3D(fieldmesh);
+  }
 }
 
 void Field3D::mergeYupYdown() {
   TRACE("Field3D::mergeYupYdown");
   
-  if(yup_field == this && ydown_field == this)
+  if (yup1_field == this && ydown1_field == this)
     return;
 
   clearYupYdown();
 
-  yup_field = this;
-  ydown_field = this;
+  yup1_field = this;
+  ydown1_field = this;
+  if (fieldmesh->ystart>1) {
+    yup2_field = this;
+    ydown2_field = this;
+  }
 }
 
 void Field3D::clearYupYdown() {
   // Delete auxiliary fields if they have been set
-  if (yup_field != nullptr && yup_field != this) {
-    delete yup_field;
+  if (yup1_field != nullptr && yup1_field != this) {
+    delete yup1_field;
+  }
+  yup1_field = nullptr;
+
+  if (ydown1_field != nullptr && ydown1_field != this) {
+    delete ydown1_field;
+  }
+  ydown1_field = nullptr;
+
+  if (yup2_field != nullptr && yup2_field != this) {
+    delete yup2_field;
+  }
+  yup2_field = nullptr;
+
+  if (ydown2_field != nullptr && ydown2_field != this) {
+    delete ydown2_field;
+  }
+  ydown2_field = nullptr;
+}
+
+Field3D& Field3D::yup(const int i) {
+  switch (i) {
+  case 1: {
+    ASSERT2(yup1_field != nullptr); // Check for communicate
+    return *yup1_field;
+  }
+  case 2: {
+    ASSERT2(yup2_field != nullptr); // Check for communicate
+    return *yup2_field;
+  }
+#if CHECK > 1
+  default: {
+    throw BoutException("There is no yup field for i=%i", i);
+   }
+#endif
+  }
+}
+
+const Field3D& Field3D::yup(const int i) const {
+  switch (i) {
+  case 1: {
+    ASSERT2(yup1_field != nullptr); // Check for communicate
+    return *yup1_field;
+  }
+  case 2: {
+    ASSERT2(yup2_field != nullptr); // Check for communicate
+    return *yup2_field;
+  }
+#if CHECK > 1
+  default: {
+    throw BoutException("There is no yup field for i=%i", i);
+   }
+#endif
   }
-  yup_field = nullptr;
+}
 
-  if (ydown_field != nullptr && ydown_field != this) {
-    delete ydown_field;
+Field3D& Field3D::ydown(const int i) {
+  switch (i) {
+  case 1: {
+    ASSERT2(ydown1_field != nullptr); // Check for communicate
+    return *ydown1_field;
+  }
+  case 2: {
+    ASSERT2(ydown2_field != nullptr); // Check for communicate
+    return *ydown2_field;
+  }
+#if CHECK > 1
+  default: {
+    throw BoutException("There is no ydown field for i=%i", i);
+  }
+#endif
+  }
+}
+
+const Field3D& Field3D::ydown(const int i) const {
+  switch (i) {
+  case 1: {
+    ASSERT2(ydown1_field != nullptr); // Check for communicate
+    return *ydown1_field;
+  }
+  case 2: {
+    ASSERT2(ydown2_field != nullptr); // Check for communicate
+    return *ydown2_field;
+  }
+#if CHECK > 1
+  default: {
+    throw BoutException("There is no ydown field for i=%i", i);
+  }
+#endif
   }
-  ydown_field = nullptr;
 }
 
 Field3D& Field3D::ynext(int dir) {
-  switch(dir) {
-  case +1:
-    return yup();
-  case -1:
-    return ydown();
-  default:
-    throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir);
+  if (dir>0) {
+    return yup(dir);
+  } else if (dir<0) {
+    return ydown(-dir);
+  } else {
+#if CHECK > 1
+    throw BoutException("Field3D: Call to ynext with strange direction %d. Should not be zero.", dir);
+#endif
   }
 }
 
 const Field3D& Field3D::ynext(int dir) const {
-  switch(dir) {
-  case +1:
-    return yup();
-  case -1:
-    return ydown();
-  default:
-    throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir);
+  if (dir>0) {
+    return yup(dir);
+  } else if (dir<0) {
+    return ydown(-dir);
+  } else {
+#if CHECK > 1
+    throw BoutException("Field3D: Call to ynext with strange direction %d. Should not be zero.", dir);
+#endif
   }
 }
 
@@ -361,7 +458,8 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
   location = rhs.location;
   
   clearYupYdown();
-  return *this;
+
+    return *this;
 }
 
 Field3D & Field3D::operator=(const Field2D &rhs) {
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 7a8ba661a5..954230331e 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -796,6 +796,11 @@ const Field3D Coordinates::Div_par(const Field3D &f, CELL_LOC outloc,
       f_B.splitYupYdown();
       f_B.yup() = f.yup() / Bxy;
       f_B.ydown() = f.ydown() / Bxy;
+      if (mesh->ystart > 1) {
+        // Have a second yup/ydown field
+        f_B.yup(2) = f.yup(2) / Bxy;
+        f_B.ydown(2) = f.ydown(2) / Bxy;
+      }
     }
     return Bxy * Grad_par(f_B, outloc, method);
   }
diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx
index 07e4084e08..c50f43b180 100644
--- a/src/mesh/index_derivs.cxx
+++ b/src/mesh/index_derivs.cxx
@@ -960,38 +960,78 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
       CELL_LOC location = var.getLocation();
 
       stencil s;
-      s.pp = nan("");
-      s.mm = nan("");
-      for (const auto &i : result.region(region)) {
-        // Set stencils
-        s.c = var[i];
-        s.p = var.yup()[i.yp()];
-        s.m = var.ydown()[i.ym()];
+      if (mesh->ystart==1) {
+        // Only one guard cell, so can only use 3-point stencils
+        s.pp = nan("");
+        s.mm = nan("");
+        for (const auto &i : result.region(region)) {
+          // Set stencils
+          s.c = var[i];
+          s.p = var.yup()[i.yp()];
+          s.m = var.ydown()[i.ym()];
 
-        if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) {
-          // Producing a stencil centred around a lower Y value
-          s.pp = s.p;
-          s.p = s.c;
-        } else if (location == CELL_YLOW) {
-          // Stencil centred around a cell centre
-          s.mm = s.m;
-          s.m = s.c;
+          if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) {
+            // Producing a stencil centred around a lower Y value
+            s.pp = s.p;
+            s.p = s.c;
+          } else if (location == CELL_YLOW) {
+            // Stencil centred around a cell centre
+            s.mm = s.m;
+            s.m = s.c;
+          }
+
+          result[i] = func(s);
         }
+      } else {
+        // Can use 5-point stencils
+        for (const auto &i : result.region(region)) {
+          // Set stencils
+          s.c = var[i];
+          s.p = var.yup()[i.yp()];
+          s.m = var.ydown()[i.ym()];
+          s.pp = var.yup(2)[i.offset(0, 2, 0)];
+          s.mm = var.ydown(2)[i.offset(0, -2, 0)];
 
-        result[i] = func(s);
+          if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) {
+            // Producing a stencil centred around a lower Y value
+            s.pp = s.p;
+            s.p = s.c;
+          } else if (location == CELL_YLOW) {
+            // Stencil centred around a cell centre
+            s.mm = s.m;
+            s.m = s.c;
+          }
+
+          result[i] = func(s);
+        }
       }
     } else {
       // Non-staggered
       stencil s;
-      s.pp = nan("");
-      s.mm = nan("");
-      for (const auto &i : result.region(region)) {
-        // Set stencils
-        s.c = var[i];
-        s.p = var.yup()[i.yp()];
-        s.m = var.ydown()[i.ym()];
+      if (mesh->ystart == 1) {
+        // Only one guard cell, so can only use 3-point stencils
+        s.pp = nan("");
+        s.mm = nan("");
+        for (const auto &i : result.region(region)) {
+          // Set stencils
+          s.c = var[i];
+          s.p = var.yup()[i.yp()];
+          s.m = var.ydown()[i.ym()];
 
-        result[i] = func(s);
+          result[i] = func(s);
+        }
+      } else {
+        // Can use 5-point stencils
+        for (const auto &i : result.region(region)) {
+          // Set stencils
+          s.c = var[i];
+          s.p = var.yup()[i.yp()];
+          s.m = var.ydown()[i.ym()];
+          s.pp = var.yup(2)[i.offset(0, 2, 0)];
+          s.mm = var.ydown(2)[i.offset(0, -2, 0)];
+
+          result[i] = func(s);
+        }
       }
     }
   } else {
@@ -1095,6 +1135,7 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
+  invalidateGuards(result); // extra check: set guard cells to NaN if CHECK>2
 #endif
 
   return result;
@@ -2235,33 +2276,66 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       // Both v and f have up/down fields
 
       stencil vval, fval;
-      vval.pp = nan("");
-      vval.mm = nan("");
-      fval.pp = nan("");
-      fval.mm = nan("");
-      for (const auto &i : result.region(region)) {
-        vval.c = v[i];
-        vval.p = v.yup()[i.yp()];
-        vval.m = v.ydown()[i.ym()];
-        fval.c = f[i];
-        fval.p = f.yup()[i.yp()];
-        fval.m = f.ydown()[i.ym()];
-
-        if (diffloc != CELL_DEFAULT) {
-          // Non-centred stencil
-          if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-            // Producing a stencil centred around a lower Y value
-            vval.pp = vval.p;
-            vval.p = vval.c;
-          } else if (vloc == CELL_YLOW) {
-            // Stencil centred around a cell centre
-            vval.mm = vval.m;
-            vval.m = vval.c;
+      if (mesh->ystart==1) {
+        // Only one guard cell, so can only use 3-point stencils
+        vval.pp = nan("");
+        vval.mm = nan("");
+        fval.pp = nan("");
+        fval.mm = nan("");
+        for (const auto &i : result.region(region)) {
+          vval.c = v[i];
+          vval.p = v.yup()[i.yp()];
+          vval.m = v.ydown()[i.ym()];
+          fval.c = f[i];
+          fval.p = f.yup()[i.yp()];
+          fval.m = f.ydown()[i.ym()];
+
+          if (diffloc != CELL_DEFAULT) {
+            // Non-centred stencil
+            if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
+              // Producing a stencil centred around a lower Y value
+              vval.pp = vval.p;
+              vval.p = vval.c;
+            } else if (vloc == CELL_YLOW) {
+              // Stencil centred around a cell centre
+              vval.mm = vval.m;
+              vval.m = vval.c;
+            }
+            // Shifted in one direction -> shift in another
+            // Could produce warning
           }
-          // Shifted in one direction -> shift in another
-          // Could produce warning
+          result[i] = func(vval, fval);
+        }
+      } else {
+        // Can use 5-point stencils
+        for (const auto &i : result.region(region)) {
+          vval.c = v[i];
+          vval.p = v.yup()[i.yp()];
+          vval.m = v.ydown()[i.ym()];
+          vval.pp = v.yup(2)[i.offset(0, 2, 0)];
+          vval.mm = v.ydown(2)[i.offset(0, -2, 0)];
+          fval.c = f[i];
+          fval.p = f.yup()[i.yp()];
+          fval.m = f.ydown()[i.ym()];
+          fval.pp = f.yup(2)[i.offset(0, 2, 0)];
+          fval.mm = f.ydown(2)[i.offset(0, -2, 0)];
+
+          if (diffloc != CELL_DEFAULT) {
+            // Non-centred stencil
+            if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
+              // Producing a stencil centred around a lower Y value
+              vval.pp = vval.p;
+              vval.p = vval.c;
+            } else if (vloc == CELL_YLOW) {
+              // Stencil centred around a cell centre
+              vval.mm = vval.m;
+              vval.m = vval.c;
+            }
+            // Shifted in one direction -> shift in another
+            // Could produce warning
+          }
+          result[i] = func(vval, fval);
         }
-        result[i] = func(vval, fval);
       }
     } else {
       // Both must shift to field aligned
@@ -2318,21 +2392,32 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       // f has yup and ydown fields which are distinct
 
       stencil fs;
-      fs.pp = nan("");
-      fs.mm = nan("");
+      if (mesh->ystart==1) {
+        // Only one guard cell, so can only use 3-point stencils
+        fs.pp = nan("");
+        fs.mm = nan("");
 
-      Field3D f_yup = f.yup();
-      Field3D f_ydown = f.ydown();
+        for (const auto &i : result.region(region)) {
 
-      for (const auto &i : result.region(region)) {
+          fs.c = f[i];
+          fs.p = f.yup()[i.yp()];
+          fs.m = f.ydown()[i.ym()];
 
-        fs.c = f[i];
-        fs.p = f_yup[i.yp()];
-        fs.m = f_ydown[i.ym()];
+          result[i] = func(v[i], fs);
+        }
+      } else {
+        // Can use 5-point stencils
+        for (const auto &i : result.region(region)) {
 
-        result[i] = func(v[i], fs);
-      }
+          fs.c = f[i];
+          fs.p = f.yup()[i.yp()];
+          fs.m = f.ydown()[i.ym()];
+          fs.pp = f.yup(2)[i.offset(0, 2, 0)];
+          fs.mm = f.ydown(2)[i.offset(0, -2, 0)];
 
+          result[i] = func(v[i], fs);
+        }
+      }
     } else {
       // Not using yup/ydown fields, so first transform to field-aligned coordinates
       // (even if one of v and f has yup/ydown fields, it doesn't make sense to
@@ -2368,15 +2453,14 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
         }
       }
       // Shift result back
-      result = this->fromFieldAligned(result);
+      result = this->fromFieldAligned(result, RGN_NOBNDRY);
     }
   }
 
-  result.setLocation(diffloc);
-
 #if CHECK > 0
   // Mark boundaries as invalid
   result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false;
+  invalidateGuards(result); // extra check: set guard cells to NaN if CHECK>2
 #endif
 
   return result;
@@ -2933,35 +3017,71 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
   if (vUseUpDown && fUseUpDown) {
     // Both v and f have up/down fields
     stencil vval, fval;
-    vval.mm = nan("");
-    vval.pp = nan("");
-    fval.mm = nan("");
-    fval.pp = nan("");
-    for (const auto &i : result.region(region)) {
+    if (mesh->ystart==1) {
+      // Only one guard cell, so can only use 3-point stencils
+      vval.mm = nan("");
+      vval.pp = nan("");
+      fval.mm = nan("");
+      fval.pp = nan("");
+      for (const auto &i : result.region(region)) {
 
-      fval.m = f.ydown()[i.ym()];
-      fval.c = f[i];
-      fval.p = f.yup()[i.yp()];
+        fval.m = f.ydown()[i.ym()];
+        fval.c = f[i];
+        fval.p = f.yup()[i.yp()];
 
-      vval.m = v.ydown()[i.ym()];
-      vval.c = v[i];
-      vval.p = v.yup()[i.yp()];
+        vval.m = v.ydown()[i.ym()];
+        vval.c = v[i];
+        vval.p = v.yup()[i.yp()];
 
-      if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) {
-        // Non-centred stencil
-        if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
-          // Producing a stencil centred around a lower Y value
-          vval.pp = vval.p;
-          vval.p  = vval.c;
-        }else if(vloc == CELL_YLOW) {
-          // Stencil centred around a cell centre
-          vval.mm = vval.m;
-          vval.m  = vval.c;
+        if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) {
+          // Non-centred stencil
+          if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
+            // Producing a stencil centred around a lower Y value
+            vval.pp = vval.p;
+            vval.p  = vval.c;
+          }else if(vloc == CELL_YLOW) {
+            // Stencil centred around a cell centre
+            vval.mm = vval.m;
+            vval.m  = vval.c;
+          }
+          // Shifted in one direction -> shift in another
+          // Could produce warning
         }
-        // Shifted in one direction -> shift in another
-        // Could produce warning
+        result[i] = func(vval, fval);
       }
-      result[i] = func(vval, fval);
+    } else {
+      // Can use 5-point stencils
+      for (const auto &i : result.region(region)) {
+
+        fval.mm = f.ydown(2)[i.offset(0, -2, 0)];
+        fval.m = f.ydown()[i.ym()];
+        fval.c = f[i];
+        fval.p = f.yup()[i.yp()];
+        fval.pp = f.yup(2)[i.offset(0, 2, 0)];
+
+        vval.mm = v.ydown(2)[i.offset(0, -2, 0)];
+        vval.m = v.ydown()[i.ym()];
+        vval.c = v[i];
+        vval.p = v.yup()[i.yp()];
+        vval.pp = v.yup(2)[i.offset(0, 2, 0)];
+
+        if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) {
+          // Non-centred stencil
+          if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) {
+            // Producing a stencil centred around a lower Y value
+            vval.pp = vval.p;
+            vval.p  = vval.c;
+          }else if(vloc == CELL_YLOW) {
+            // Stencil centred around a cell centre
+            vval.mm = vval.m;
+            vval.m  = vval.c;
+          }
+          // Shifted in one direction -> shift in another
+          // Could produce warning
+        }
+        result[i] = func(vval, fval);
+      }
+
     }
   }
   else {
diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx
index d018041a74..93baef9ad5 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -120,29 +120,53 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
           // Field "var" has distinct yup and ydown fields which
           // will be used to calculate a derivative along
           // the magnetic field
-          throw BoutException("At the moment, fields with yup/ydown cannot use interp_to.\n"
-                              "If we implement a 3-point stencil for interpolate or double-up\n"
-                              "/double-down fields, then we can use this case.");
-          s.pp = nan("");
-          s.mm = nan("");
+          // More than one guard cell, so set pp and mm values
+          // This allows higher-order methods to be used
+          if (fieldmesh->ystart > 1) {
+            for(const auto &i : result.region(RGN_NOY)) {
+              // Set stencils
+              s.c = var[i];
+              s.p = var.yup()[i.yp()];
+              s.m = var.ydown()[i.ym()];
+              s.pp = var.yup(2)[i.offset(0,2,0)];
+              s.mm = var.ydown(2)[i.offset(0,-2,0)];
 
-          for (const auto &i : result.region(RGN_NOBNDRY)) {
-            // Set stencils
-            s.c = var[i];
-            s.p = var.yup()[i.yp()];
-            s.m = var.ydown()[i.ym()];
+              if (location == CELL_CENTRE) {
+                // Producing a stencil centred around a lower Y value
+                s.pp = s.p;
+                s.p  = s.c;
+              } else {
+                // Stencil centred around a cell centre
+                s.mm = s.m;
+                s.m  = s.c;
+              }
 
-            if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) {
-              // Producing a stencil centred around a lower Y value
-              s.pp = s.p;
-              s.p = s.c;
-            } else if (location == CELL_YLOW) {
-              // Stencil centred around a cell centre
-              s.mm = s.m;
-              s.m = s.c;
+              result[i] = interp(s);
             }
+          } else {
+            // Note: at the moment we cannot reach this case because of the
+            // 'ASSERT0(mesh->ystart >=2)' above, but if we implement a 3-point
+            // stencil for interp, then this will be useful
+            s.pp = nan("");
+            s.mm = nan("");
+            for(const auto &i : result.region(RGN_NOY)) {
+              // Set stencils
+              s.c = var[i];
+              s.p = var.yup()[i.yp()];
+              s.m = var.ydown()[i.ym()];
 
-            result[i] = interp(s);
+              if (location == CELL_CENTRE) {
+                // Producing a stencil centred around a lower Y value
+                s.pp = s.p;
+                s.p  = s.c;
+              } else {
+                // Stencil centred around a cell centre
+                s.mm = s.m;
+                s.m  = s.c;
+              }
+
+              result[i] = interp(s);
+            }
           }
         } else {
           // var has no yup/ydown fields, so we need to shift into field-aligned
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 614e3376dd..581ea21419 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -45,21 +45,27 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
   fromAlignedPhs.resize(mesh.LocalNx);
   toAlignedPhs.resize(mesh.LocalNx);
   
-  yupPhs.resize(mesh.LocalNx);
-  ydownPhs.resize(mesh.LocalNx);
+  yupPhs1.resize(mesh.LocalNx);
+  ydownPhs1.resize(mesh.LocalNx);
+  yupPhs2.resize(mesh.LocalNx);
+  ydownPhs2.resize(mesh.LocalNx);
 
   for(int jx=0;jx<mesh.LocalNx;jx++){
     fromAlignedPhs[jx].resize(mesh.LocalNy);
     toAlignedPhs[jx].resize(mesh.LocalNy);
 
-    yupPhs[jx].resize(mesh.LocalNy);
-    ydownPhs[jx].resize(mesh.LocalNy);
+    yupPhs1[jx].resize(mesh.LocalNy);
+    ydownPhs1[jx].resize(mesh.LocalNy);
+    yupPhs2[jx].resize(mesh.LocalNy);
+    ydownPhs2[jx].resize(mesh.LocalNy);
     for(int jy=0;jy<mesh.LocalNy;jy++){
       fromAlignedPhs[jx][jy].resize(nmodes);
       toAlignedPhs[jx][jy].resize(nmodes);
       
-      yupPhs[jx][jy].resize(nmodes);
-      ydownPhs[jx][jy].resize(nmodes);
+      yupPhs1[jx][jy].resize(nmodes);
+      ydownPhs1[jx][jy].resize(nmodes);
+      yupPhs2[jx][jy].resize(nmodes);
+      ydownPhs2[jx][jy].resize(nmodes);
     }
   }
 	
@@ -77,14 +83,30 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
   //Yup/Ydown phases -- note we don't shift in the boundaries/guards
   for(int jx=0;jx<mesh.LocalNx;jx++){
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++){
-      BoutReal yupShift = zShift(jx,jy) - zShift(jx,jy+1);
-      BoutReal ydownShift = zShift(jx,jy) - zShift(jx,jy-1);
+      BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
+      BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
       
       for(int jz=0;jz<nmodes;jz++) {
   	BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
 
-  	yupPhs[jx][jy][jz] = dcomplex(cos(kwave*yupShift) , -sin(kwave*yupShift));
-  	ydownPhs[jx][jy][jz] = dcomplex(cos(kwave*ydownShift) , -sin(kwave*ydownShift));
+        yupPhs1[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+        ydownPhs1[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+      }
+    }
+  }
+  if (mesh.ystart>1) {
+    // Need phases for the second yup/ydown fields too
+    for(int jx=0;jx<mesh.LocalNx;jx++){
+      for(int jy=mesh.ystart;jy<=mesh.yend;jy++){
+        BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
+        BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
+
+        for(int jz=0;jz<nmodes;jz++) {
+          BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+
+          yupPhs2[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+          ydownPhs2[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+        }
       }
     }
   }
@@ -97,21 +119,37 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
 void ShiftedMetric::calcYUpDown(Field3D &f) {
   f.splitYupYdown();
   
-  Field3D& yup = f.yup();
-  yup.allocate();
-
+  Field3D& yup1 = f.yup();
+  yup1.allocate();
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy+1,0)), yupPhs[jx][jy], &(yup(jx,jy+1,0)));
+      shiftZ(&(f(jx,jy+1,0)), yupPhs1[jx][jy], &(yup1(jx,jy+1,0)));
+    }
+  }
+  if (mesh.ystart>1) {
+    Field3D& yup2 = f.yup(2);
+    yup2.allocate();
+    for(int jx=0;jx<mesh.LocalNx;jx++) {
+      for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
+        shiftZ(&(f(jx,jy+2,0)), yupPhs2[jx][jy], &(yup2(jx,jy+2,0)));
+      }
     }
   }
 
-  Field3D& ydown = f.ydown();
-  ydown.allocate();
-
+  Field3D& ydown1 = f.ydown();
+  ydown1.allocate();
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy-1,0)), ydownPhs[jx][jy], &(ydown(jx,jy-1,0)));
+      shiftZ(&(f(jx,jy-1,0)), ydownPhs1[jx][jy], &(ydown1(jx,jy-1,0)));
+    }
+  }
+  if (mesh.ystart > 1) {
+    Field3D& ydown2 = f.ydown(2);
+    ydown2.allocate();
+    for(int jx=0;jx<mesh.LocalNx;jx++) {
+      for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
+        shiftZ(&(f(jx,jy-2,0)), ydownPhs2[jx][jy], &(ydown2(jx,jy-2,0)));
+      }
     }
   }
 }

From 4d603bcc55e0468371ae190beb080621a9ec5a5c Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Fri, 27 Apr 2018 15:48:46 +0100
Subject: [PATCH 12/21] Add check that MYG=1 to FCITransform

FCITransform does not support more than one yup/ydown field, so require
it to use only one y-guard cell. Check for this in the constructor.

Fix test-fci-slab to use myg=1.
---
 src/mesh/parallel/fci.hxx                    | 18 ++++++++++++++++--
 tests/integrated/test-fci-slab/data/BOUT.inp |  2 ++
 tests/integrated/test-fci-slab/mms/BOUT.inp  |  2 ++
 3 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx
index 017375b668..915356b50e 100644
--- a/src/mesh/parallel/fci.hxx
+++ b/src/mesh/parallel/fci.hxx
@@ -69,10 +69,24 @@ public:
 class FCITransform : public ParallelTransform {
 public:
   DEPRECATED(FCITransform(Mesh &mesh, bool UNUSED(yperiodic), bool zperiodic))
-      : FCITransform(mesh, zperiodic) {}
+      : FCITransform(mesh, zperiodic) {
+        if (mesh.ystart > 1)
+          // FCITransform can only use myg=1 because it only loads grid
+          // information for one point forward or back along the magnetic
+          // field, so it cannot set Field3D::yup2_field or
+          // Field3D::ydown2_field
+          throw BoutException("FCI method must use only one y-guard cell: set option myg=1");
+      }
   FCITransform(Mesh &mesh, bool zperiodic = true)
       : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic),
-        zperiodic(zperiodic) {}
+        zperiodic(zperiodic) {
+          if (mesh.ystart > 1)
+            // FCITransform can only use myg=1 because it only loads grid
+            // information for one point forward or back along the magnetic
+            // field, so it cannot set Field3D::yup2_field or
+            // Field3D::ydown2_field
+            throw BoutException("FCI method must use only one y-guard cell: set option myg=1");
+        }
 
   void calcYUpDown(Field3D &f) override;
   
diff --git a/tests/integrated/test-fci-slab/data/BOUT.inp b/tests/integrated/test-fci-slab/data/BOUT.inp
index 218946e676..dc0ff71b81 100644
--- a/tests/integrated/test-fci-slab/data/BOUT.inp
+++ b/tests/integrated/test-fci-slab/data/BOUT.inp
@@ -6,6 +6,8 @@ timestep = 0.2
 
 MZ = 64
 
+myg = 1
+
 [mesh]
 paralleltransform = fci
 
diff --git a/tests/integrated/test-fci-slab/mms/BOUT.inp b/tests/integrated/test-fci-slab/mms/BOUT.inp
index 3c1e5f9559..5280f748be 100644
--- a/tests/integrated/test-fci-slab/mms/BOUT.inp
+++ b/tests/integrated/test-fci-slab/mms/BOUT.inp
@@ -7,6 +7,8 @@ MZ = 64
 
 NXPE = 1
 
+myg = 1
+
 [mesh]
 paralleltransform = fci
 

From db7a2d9e0fda2cb7f3508131470810f04307c095 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Fri, 27 Apr 2018 17:53:52 +0100
Subject: [PATCH 13/21] Add unit tests for second yup/ydown fields

---
 tests/unit/field/test_field3d.cxx | 95 ++++++++++++++++++++++++++++++-
 1 file changed, 94 insertions(+), 1 deletion(-)

diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx
index 2fb7e612c5..bd9a1b4872 100644
--- a/tests/unit/field/test_field3d.cxx
+++ b/tests/unit/field/test_field3d.cxx
@@ -12,8 +12,9 @@
 #include <set>
 #include <vector>
 
-/// Global mesh
+/// Global meshes
 extern Mesh *mesh;
+Mesh *mesh2; // mesh2 has 2 guard cells
 
 /// Test fixture to make sure the global mesh is our fake one
 class Field3DTest : public ::testing::Test {
@@ -26,11 +27,24 @@ class Field3DTest : public ::testing::Test {
     }
     mesh = new FakeMesh(nx, ny, nz);
     mesh->createDefaultRegions();
+
+    if (mesh2 != nullptr) {
+      delete mesh2;
+      mesh2 = nullptr;
+    }
+    mesh2 = new FakeMesh(nx, ny, nz);
+    mesh2->createDefaultRegions();
+    mesh2->xstart += 1;
+    mesh2->xend -= 1;
+    mesh2->ystart += 1;
+    mesh2->yend -=1;
   }
 
   static void TearDownTestCase() {
     delete mesh;
     mesh = nullptr;
+    delete mesh2;
+    mesh2 = nullptr;
   }
 
 public:
@@ -275,6 +289,34 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) {
   EXPECT_EQ(&field, &ydown2);
 }
 
+TEST_F(Field3DTest, SplitThenMergeYupYDown2) {
+  Field3D field(mesh2);
+
+  field = 0.;
+  field.splitYupYdown();
+
+  auto& yup1 = field.yup();
+  EXPECT_NE(&field, &yup1);
+  auto& ydown1 = field.ydown();
+  EXPECT_NE(&field, &ydown1);
+  auto& yup2 = field.yup(2);
+  EXPECT_NE(&field, &yup2);
+  auto& ydown2 = field.ydown(2);
+  EXPECT_NE(&field, &ydown2);
+
+  field.mergeYupYdown();
+
+  auto& yup1_2 = field.yup();
+  EXPECT_EQ(&field, &yup1_2);
+  auto& ydown1_2 = field.ydown();
+  EXPECT_EQ(&field, &ydown1_2);
+
+  auto& yup2_2 = field.yup(2);
+  EXPECT_EQ(&field, &yup2_2);
+  auto& ydown2_2 = field.ydown(2);
+  EXPECT_EQ(&field, &ydown2_2);
+}
+
 TEST_F(Field3DTest, Ynext) {
   Field3D field;
 
@@ -287,7 +329,9 @@ TEST_F(Field3DTest, Ynext) {
   EXPECT_NE(&field, &ydown);
   EXPECT_NE(&yup, &ydown);
 
+#if CHECK > 1
   EXPECT_THROW(field.ynext(99), BoutException);
+#endif
 }
 
 TEST_F(Field3DTest, ConstYnext) {
@@ -303,7 +347,56 @@ TEST_F(Field3DTest, ConstYnext) {
   EXPECT_NE(&field2, &ydown);
   EXPECT_NE(&yup, &ydown);
 
+#if CHECK > 1
   EXPECT_THROW(field2.ynext(99), BoutException);
+#endif
+}
+
+TEST_F(Field3DTest, Ynext2) {
+  Field3D field(mesh2);
+
+  field = 0.;
+  field.splitYupYdown();
+
+  auto& yup = field.ynext(1);
+  EXPECT_NE(&field, &yup);
+  auto& ydown = field.ynext(-1);
+  EXPECT_NE(&field, &ydown);
+  EXPECT_NE(&yup, &ydown);
+
+  auto& yup2 = field.ynext(2);
+  EXPECT_NE(&field, &yup2);
+  auto& ydown2 = field.ynext(-2);
+  EXPECT_NE(&field, &ydown2);
+  EXPECT_NE(&yup2, &ydown2);
+
+#if CHECK > 1
+  EXPECT_THROW(field.ynext(99), BoutException);
+#endif
+}
+
+TEST_F(Field3DTest, ConstYnext2) {
+  Field3D field(0., mesh2);
+
+  field.splitYupYdown();
+
+  const Field3D& field2 = field;
+
+  auto& yup = field2.ynext(1);
+  EXPECT_NE(&field2, &yup);
+  auto& ydown = field2.ynext(-1);
+  EXPECT_NE(&field2, &ydown);
+  EXPECT_NE(&yup, &ydown);
+
+  auto& yup2 = field2.ynext(2);
+  EXPECT_NE(&field2, &yup2);
+  auto& ydown2 = field2.ynext(-2);
+  EXPECT_NE(&field2, &ydown2);
+  EXPECT_NE(&yup2, &ydown2);
+
+#if CHECK > 1
+  EXPECT_THROW(field2.ynext(99), BoutException);
+#endif
 }
 
 TEST_F(Field3DTest, GetGlobalMesh) {

From bd908a5c7a622d5cba7f7fefbbc3b8270f71973b Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Fri, 4 May 2018 16:12:28 +0100
Subject: [PATCH 14/21] Implemented staggered grids in ShiftedMetric

---
 include/bout/paralleltransform.hxx  |  33 +-
 src/mesh/parallel/shiftedmetric.cxx | 728 +++++++++++++++++++++++++---
 2 files changed, 682 insertions(+), 79 deletions(-)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index 2010d327c0..102468a1b3 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -140,13 +140,32 @@ private:
   std::vector<dcomplex> cmplx; ///< A temporary array, used for input/output to fft routines
   std::vector<dcomplex> cmplxLoc; ///< A temporary array, used for input/output to fft routines
 
-  arr3Dvec toAlignedPhs; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates
-  arr3Dvec fromAlignedPhs; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates
-
-  arr3Dvec yupPhs1; ///< Cache of phase shifts for calculating yup1 fields
-  arr3Dvec ydownPhs1; ///< Cache of phase shifts for calculating ydown1 fields
-  arr3Dvec yupPhs2; ///< Cache of phase shifts for calculating yup2 fields
-  arr3Dvec ydownPhs2; ///< Cache of phase shifts for calculating ydown2 fields
+  arr3Dvec getToAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  arr3Dvec getFromAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  arr3Dvec getYupPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  arr3Dvec getYdownPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  arr3Dvec getYupPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  arr3Dvec getYdownPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+
+  arr3Dvec toAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Cell centre version.
+  arr3Dvec fromAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Cell centre version.
+  arr3Dvec toAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_XLOW.
+  arr3Dvec fromAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_XLOW.
+  arr3Dvec toAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_YLOW.
+  arr3Dvec fromAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_YLOW.
+
+  arr3Dvec yupPhs1_CENTRE; ///< Cache of phase shifts for calculating yup1 fields. Cell centre version.
+  arr3Dvec ydownPhs1_CENTRE; ///< Cache of phase shifts for calculating ydown1 fields. Cell centre version.
+  arr3Dvec yupPhs2_CENTRE; ///< Cache of phase shifts for calculating yup2 fields. Cell centre version.
+  arr3Dvec ydownPhs2_CENTRE; ///< Cache of phase shifts for calculating ydown2 fields. Cell centre version.
+  arr3Dvec yupPhs1_XLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_XLOW.
+  arr3Dvec ydownPhs1_XLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_XLOW.
+  arr3Dvec yupPhs2_XLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_XLOW.
+  arr3Dvec ydownPhs2_XLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_XLOW.
+  arr3Dvec yupPhs1_YLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_YLOW.
+  arr3Dvec ydownPhs1_YLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_YLOW.
+  arr3Dvec yupPhs2_YLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_YLOW.
+  arr3Dvec ydownPhs2_YLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_YLOW.
 
   /*!
    * Shift a 2D field in Z. 
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 581ea21419..65ff0a2c8b 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -8,6 +8,7 @@
 
 #include <bout/paralleltransform.hxx>
 #include <bout/mesh.hxx>
+#include <interpolation.hxx>
 #include <fft.hxx>
 #include <bout/constants.hxx>
 
@@ -23,94 +24,673 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
     mesh.get(zShift, "qinty");
   }
 
-  //If we wanted to be efficient we could move the following cached phase setup
-  //into the relevant shifting routines (with static bool first protection)
-  //so that we only calculate the phase if we actually call a relevant shift 
-  //routine -- however as we're only going to do this initialisation once I 
-  //think it's cleaner to put it in the constructor here.
-
-  //As we're attached to a mesh we can expect the z direction to
-  //not change once we've been created so precalculate the complex
-  //phases used in transformations
   int nmodes = mesh.LocalNz/2 + 1;
-  BoutReal zlength = mesh.coordinates()->zlength();
-
   //Allocate storage for complex intermediate
   cmplx.resize(nmodes);
   std::fill(cmplx.begin(), cmplx.end(), 0.0);
+}
 
-  //Allocate storage for our 3d vector structures.
-  //This could be made more succinct but this approach is fairly
-  //verbose --> transparent
-  fromAlignedPhs.resize(mesh.LocalNx);
-  toAlignedPhs.resize(mesh.LocalNx);
-  
-  yupPhs1.resize(mesh.LocalNx);
-  ydownPhs1.resize(mesh.LocalNx);
-  yupPhs2.resize(mesh.LocalNx);
-  ydownPhs2.resize(mesh.LocalNx);
+//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
+ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
+
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      first_CENTRE = false;
+      fromAlignedPhs_CENTRE.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        fromAlignedPhs_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          fromAlignedPhs_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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(jx,jy)) , -sin(kwave*zShift(jx,jy)));
+          }
+        }
+      }
+    }
+    return fromAlignedPhs_CENTRE;
+    break;
+  }
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      fromAlignedPhs_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        fromAlignedPhs_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          fromAlignedPhs_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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 (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      fromAlignedPhs_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        fromAlignedPhs_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          fromAlignedPhs_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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;
+  }
+  };
+}
 
-  for(int jx=0;jx<mesh.LocalNx;jx++){
-    fromAlignedPhs[jx].resize(mesh.LocalNy);
-    toAlignedPhs[jx].resize(mesh.LocalNy);
+ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
-    yupPhs1[jx].resize(mesh.LocalNy);
-    ydownPhs1[jx].resize(mesh.LocalNy);
-    yupPhs2[jx].resize(mesh.LocalNy);
-    ydownPhs2[jx].resize(mesh.LocalNy);
-    for(int jy=0;jy<mesh.LocalNy;jy++){
-      fromAlignedPhs[jx][jy].resize(nmodes);
-      toAlignedPhs[jx][jy].resize(nmodes);
-      
-      yupPhs1[jx][jy].resize(nmodes);
-      ydownPhs1[jx][jy].resize(nmodes);
-      yupPhs2[jx][jy].resize(nmodes);
-      ydownPhs2[jx][jy].resize(nmodes);
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      first_CENTRE = false;
+      toAlignedPhs_CENTRE.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        toAlignedPhs_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          toAlignedPhs_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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(jx,jy)), sin(kwave*zShift(jx,jy)));
+          }
+        }
+      }
     }
+    return toAlignedPhs_CENTRE;
+    break;
   }
-	
-  //To/From field aligned phases
-  for(int jx=0;jx<mesh.LocalNx;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[jx][jy][jz] = dcomplex(cos(kwave*zShift(jx,jy)) , -sin(kwave*zShift(jx,jy)));
-  	toAlignedPhs[jx][jy][jz] =   dcomplex(cos(kwave*zShift(jx,jy)) ,  sin(kwave*zShift(jx,jy)));
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      toAlignedPhs_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        toAlignedPhs_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          toAlignedPhs_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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 (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      toAlignedPhs_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        toAlignedPhs_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          toAlignedPhs_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;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;
   }
+  };
+}
+
+ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
-  //Yup/Ydown phases -- note we don't shift in the boundaries/guards
-  for(int jx=0;jx<mesh.LocalNx;jx++){
-    for(int jy=mesh.ystart;jy<=mesh.yend;jy++){
-      BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
-      BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
-      
-      for(int jz=0;jz<nmodes;jz++) {
-  	BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
 
-        yupPhs1[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
-        ydownPhs1[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+      first_CENTRE = false;
+      yupPhs1_CENTRE.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs1_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs1_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs1_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+          }
+        }
       }
     }
+    return yupPhs1_CENTRE;
+    break;
   }
-  if (mesh.ystart>1) {
-    // Need phases for the second yup/ydown fields too
-    for(int jx=0;jx<mesh.LocalNx;jx++){
-      for(int jy=mesh.ystart;jy<=mesh.yend;jy++){
-        BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
-        BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      yupPhs1_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs1_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs1_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs1_XLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+          }
+        }
+      }
+    }
+    return yupPhs1_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      yupPhs1_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs1_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs1_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs1_YLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+          }
+        }
+      }
+    }
+    return yupPhs1_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getYupPhs1(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;
+  }
+  };
+}
+
+ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
+
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
 
-        for(int jz=0;jz<nmodes;jz++) {
-          BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+      first_CENTRE = false;
+      yupPhs2_CENTRE.resize(mesh.LocalNx);
 
-          yupPhs2[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
-          ydownPhs2[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs2_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs2_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs2_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+          }
         }
       }
     }
+    return yupPhs2_CENTRE;
+    break;
   }
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      yupPhs2_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs2_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs2_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs2_XLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+          }
+        }
+      }
+    }
+    return yupPhs2_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      yupPhs2_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        yupPhs2_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          yupPhs2_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal yupShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            yupPhs2_YLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+          }
+        }
+      }
+    }
+    return yupPhs2_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getYupPhs1(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;
+  }
+  };
+}
 
+ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
+
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      first_CENTRE = false;
+      ydownPhs1_CENTRE.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs1_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs1_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs1_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+          }
+        }
+      }
+    }
+    return ydownPhs1_CENTRE;
+    break;
+  }
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      ydownPhs1_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs1_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs1_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs1_XLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+          }
+        }
+      }
+    }
+    return ydownPhs1_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      ydownPhs1_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs1_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs1_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-1);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs1_YLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+          }
+        }
+      }
+    }
+    return ydownPhs1_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getYdownPhs1(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;
+  }
+  };
+}
+
+ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
+  // bools so we only calculate the cached values the first time for each location
+  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
+
+  switch (location) {
+  case CELL_CENTRE: {
+    if (first_CENTRE) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      first_CENTRE = false;
+      ydownPhs2_CENTRE.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs2_CENTRE[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs2_CENTRE[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs2_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+          }
+        }
+      }
+    }
+    return ydownPhs2_CENTRE;
+    break;
+  }
+  case CELL_XLOW: {
+    if (first_XLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to XLOW
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_XLOW = false;
+      ydownPhs2_XLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs2_XLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs2_XLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs2_XLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+          }
+        }
+      }
+    }
+    return ydownPhs2_XLOW;
+    break;
+  }
+  case CELL_YLOW: {
+    if (first_YLOW) {
+      int nmodes = mesh.LocalNz/2 + 1;
+      BoutReal zlength = mesh.coordinates()->zlength();
+
+      // interpolate zShift to YLOW
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+
+      first_YLOW = false;
+      ydownPhs2_YLOW.resize(mesh.LocalNx);
+
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        ydownPhs2_YLOW[jx].resize(mesh.LocalNy);
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          ydownPhs2_YLOW[jx][jy].resize(nmodes);
+        }
+      }
+      //To/From field aligned phases
+      for(int jx=0;jx<mesh.LocalNx;jx++){
+        for(int jy=0;jy<mesh.LocalNy;jy++){
+          BoutReal ydownShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-2);
+          for(int jz=0;jz<nmodes;jz++) {
+            BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+            ydownPhs2_YLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+          }
+        }
+      }
+    }
+    return ydownPhs2_YLOW;
+    break;
+  }
+  case CELL_ZLOW: {
+    // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
+    return getYdownPhs1(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;
+  }
+  };
 }
 
 /*!
@@ -118,37 +698,42 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
  */
 void ShiftedMetric::calcYUpDown(Field3D &f) {
   f.splitYupYdown();
+  CELL_LOC location = f.getLocation();
   
   Field3D& yup1 = f.yup();
   yup1.allocate();
+  arr3Dvec phases = getYupPhs1(location);
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy+1,0)), yupPhs1[jx][jy], &(yup1(jx,jy+1,0)));
+      shiftZ(&(f(jx,jy+1,0)), phases[jx][jy], &(yup1(jx,jy+1,0)));
     }
   }
   if (mesh.ystart>1) {
     Field3D& yup2 = f.yup(2);
     yup2.allocate();
+    phases = getYupPhs2(location);
     for(int jx=0;jx<mesh.LocalNx;jx++) {
       for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-        shiftZ(&(f(jx,jy+2,0)), yupPhs2[jx][jy], &(yup2(jx,jy+2,0)));
+        shiftZ(&(f(jx,jy+2,0)), phases[jx][jy], &(yup2(jx,jy+2,0)));
       }
     }
   }
 
   Field3D& ydown1 = f.ydown();
   ydown1.allocate();
+  phases = getYdownPhs1(location);
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy-1,0)), ydownPhs1[jx][jy], &(ydown1(jx,jy-1,0)));
+      shiftZ(&(f(jx,jy-1,0)), phases[jx][jy], &(ydown1(jx,jy-1,0)));
     }
   }
   if (mesh.ystart > 1) {
     Field3D& ydown2 = f.ydown(2);
     ydown2.allocate();
+    phases = getYdownPhs2(location);
     for(int jx=0;jx<mesh.LocalNx;jx++) {
       for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-        shiftZ(&(f(jx,jy-2,0)), ydownPhs2[jx][jy], &(ydown2(jx,jy-2,0)));
+        shiftZ(&(f(jx,jy-2,0)), phases[jx][jy], &(ydown2(jx,jy-2,0)));
       }
     }
   }
@@ -159,7 +744,7 @@ 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);
+  return shiftZ(f, getToAlignedPhs(f.getLocation()), region);
 }
 
 /*!
@@ -167,7 +752,7 @@ 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);
+  return shiftZ(f, getFromAlignedPhs(f.getLocation()), region);
 }
 
 const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region) {
@@ -176,8 +761,7 @@ const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const
   if(mesh.LocalNz == 1)
     return f; // Shifting makes no difference
 
-  Field3D result(&mesh);
-  result.allocate();
+  Field3D result(f); // Initialize from f, mostly so location get set correctly. (Does not copy data because of copy-on-change).
 
   for(auto i : f.region2D(region)) {
     shiftZ(f(i.x,i.y), phs[i.x][i.y], result(i.x,i.y));

From 8720e08a0ac030e1d144748b5bb138d2f158d2ae Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 17 Jul 2018 17:45:30 +0100
Subject: [PATCH 15/21] Fix test-yupdown

Due to using static bool to cache phases, can only have one
ShiftedMetric object. Therefore need to update test to use the
ShiftedMetric in mesh instead of creating a new one.
---
 tests/integrated/test-yupdown/test_yupdown.cxx | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx
index e414c0ba88..f7c9e02b79 100644
--- a/tests/integrated/test-yupdown/test_yupdown.cxx
+++ b/tests/integrated/test-yupdown/test_yupdown.cxx
@@ -37,8 +37,6 @@ int main(int argc, char** argv) {
 
   BoutInitialise(argc, argv);
 
-  ShiftedMetric s(*mesh);
-
   // Read variable from mesh
   Field3D var;
   mesh->get(var, "var");
@@ -46,7 +44,7 @@ int main(int argc, char** argv) {
   // Var starts in orthogonal X-Z coordinates
 
   // Calculate yup and ydown
-  s.calcYUpDown(var);
+  mesh->communicate(var);
   
   // Calculate d/dy ysing yup() and ydown() fields
   Field3D ddy = DDY_yud(var);

From 865d4293dc7c2daaf260c9950f8c5ea73a12b74b Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 14 May 2018 21:00:12 +0100
Subject: [PATCH 16/21] Use Flexible<Field2D> for ShiftedMetric::zShift

Also add checks for locations of Field2D variables.
---
 include/bout/paralleltransform.hxx  |  3 +-
 src/mesh/parallel/shiftedmetric.cxx | 89 +++++++++++++++++------------
 2 files changed, 55 insertions(+), 37 deletions(-)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index 102468a1b3..4be9ad8e70 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -8,6 +8,7 @@
 
 #include <datafile.hxx>
 #include <field3d.hxx>
+#include <flexible.hxx>
 #include <boutexception.hxx>
 #include <dcomplex.hxx>
 #include <unused.hxx>
@@ -136,7 +137,6 @@ private:
 
   /// This is the shift in toroidal angle (z) which takes a point from
   /// X-Z orthogonal to field-aligned along Y.
-  Field2D zShift;
   std::vector<dcomplex> cmplx; ///< A temporary array, used for input/output to fft routines
   std::vector<dcomplex> cmplxLoc; ///< A temporary array, used for input/output to fft routines
 
@@ -166,6 +166,7 @@ private:
   arr3Dvec ydownPhs1_YLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_YLOW.
   arr3Dvec yupPhs2_YLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_YLOW.
   arr3Dvec ydownPhs2_YLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_YLOW.
+  Flexible<Field2D> zShift;
 
   /*!
    * Shift a 2D field in Z. 
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 65ff0a2c8b..62817e2f84 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -24,6 +24,23 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
     mesh.get(zShift, "qinty");
   }
 
+  if (mesh.StaggerGrids) {
+    if (mesh.xstart >=2) {
+      // Can interpolate in x-direction
+      // Calculate staggered field for zShift and apply boundary conditions
+      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
+      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value
+      zShift.set(zShift_XLOW);
+    }
+    if (mesh.ystart >=2) {
+      // Can interpolate in y-direction
+      // Calculate staggered field for zShift and apply boundary conditions
+      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
+      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value
+      zShift.set(zShift_YLOW);
+    }
+  }
+
   int nmodes = mesh.LocalNz/2 + 1;
   //Allocate storage for complex intermediate
   cmplx.resize(nmodes);
@@ -67,12 +84,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       fromAlignedPhs_XLOW.resize(mesh.LocalNx);
@@ -98,12 +115,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       fromAlignedPhs_YLOW.resize(mesh.LocalNx);
@@ -175,12 +192,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       toAlignedPhs_XLOW.resize(mesh.LocalNx);
@@ -206,12 +223,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       toAlignedPhs_YLOW.resize(mesh.LocalNx);
@@ -284,12 +301,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       yupPhs1_XLOW.resize(mesh.LocalNx);
@@ -316,12 +333,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       yupPhs1_YLOW.resize(mesh.LocalNx);
@@ -395,12 +412,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       yupPhs2_XLOW.resize(mesh.LocalNx);
@@ -427,12 +444,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       yupPhs2_YLOW.resize(mesh.LocalNx);
@@ -506,12 +523,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       ydownPhs1_XLOW.resize(mesh.LocalNx);
@@ -538,12 +555,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       ydownPhs1_YLOW.resize(mesh.LocalNx);
@@ -617,12 +634,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
+      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to XLOW
-      Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL);
-      zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at XLOW
+      Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
       ydownPhs2_XLOW.resize(mesh.LocalNx);
@@ -649,12 +666,12 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
+      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      // interpolate zShift to YLOW
-      Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL);
-      zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells equal to nearest grid cell
+      // get zShift at YLOW
+      Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
       ydownPhs2_YLOW.resize(mesh.LocalNx);

From d486464474483cb9e96ee9483e69bb3ef8400add Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 14 May 2018 21:01:13 +0100
Subject: [PATCH 17/21] Use Matrix/Array instead of std::vector

In ShiftedMetric, replace std::vector<dcomplex> with Array<dcomplex> and
nested std::vector with Matrix< Array<dcomplex> >.
---
 include/bout/paralleltransform.hxx  |  65 +++---
 src/mesh/parallel/shiftedmetric.cxx | 302 ++++++++++++----------------
 2 files changed, 156 insertions(+), 211 deletions(-)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index 4be9ad8e70..2c71e4f0ff 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -12,6 +12,7 @@
 #include <boutexception.hxx>
 #include <dcomplex.hxx>
 #include <unused.hxx>
+#include <utils.hxx>
 
 class Mesh;
 
@@ -128,8 +129,6 @@ public:
     return true;
   }
 
-  /// A 3D array, implemented as nested vectors
-  typedef std::vector<std::vector<std::vector<dcomplex>>> arr3Dvec;
 private:
   ShiftedMetric();
 
@@ -137,36 +136,36 @@ private:
 
   /// This is the shift in toroidal angle (z) which takes a point from
   /// X-Z orthogonal to field-aligned along Y.
-  std::vector<dcomplex> cmplx; ///< A temporary array, used for input/output to fft routines
-  std::vector<dcomplex> cmplxLoc; ///< A temporary array, used for input/output to fft routines
-
-  arr3Dvec getToAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-  arr3Dvec getFromAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-  arr3Dvec getYupPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-  arr3Dvec getYdownPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-  arr3Dvec getYupPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-  arr3Dvec getYdownPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
-
-  arr3Dvec toAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Cell centre version.
-  arr3Dvec fromAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Cell centre version.
-  arr3Dvec toAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_XLOW.
-  arr3Dvec fromAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_XLOW.
-  arr3Dvec toAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_YLOW.
-  arr3Dvec fromAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_YLOW.
-
-  arr3Dvec yupPhs1_CENTRE; ///< Cache of phase shifts for calculating yup1 fields. Cell centre version.
-  arr3Dvec ydownPhs1_CENTRE; ///< Cache of phase shifts for calculating ydown1 fields. Cell centre version.
-  arr3Dvec yupPhs2_CENTRE; ///< Cache of phase shifts for calculating yup2 fields. Cell centre version.
-  arr3Dvec ydownPhs2_CENTRE; ///< Cache of phase shifts for calculating ydown2 fields. Cell centre version.
-  arr3Dvec yupPhs1_XLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_XLOW.
-  arr3Dvec ydownPhs1_XLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_XLOW.
-  arr3Dvec yupPhs2_XLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_XLOW.
-  arr3Dvec ydownPhs2_XLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_XLOW.
-  arr3Dvec yupPhs1_YLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_YLOW.
-  arr3Dvec ydownPhs1_YLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_YLOW.
-  arr3Dvec yupPhs2_YLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_YLOW.
-  arr3Dvec ydownPhs2_YLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_YLOW.
   Flexible<Field2D> zShift;
+  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;
+  Matrix< Array<dcomplex> > getYupPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  Matrix< Array<dcomplex> > getYdownPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  Matrix< Array<dcomplex> > getYupPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+  Matrix< Array<dcomplex> > getYdownPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
+
+  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.
+
+  Matrix< Array<dcomplex> > yupPhs1_CENTRE; ///< Cache of phase shifts for calculating yup1 fields. Cell centre version.
+  Matrix< Array<dcomplex> > ydownPhs1_CENTRE; ///< Cache of phase shifts for calculating ydown1 fields. Cell centre version.
+  Matrix< Array<dcomplex> > yupPhs2_CENTRE; ///< Cache of phase shifts for calculating yup2 fields. Cell centre version.
+  Matrix< Array<dcomplex> > ydownPhs2_CENTRE; ///< Cache of phase shifts for calculating ydown2 fields. Cell centre version.
+  Matrix< Array<dcomplex> > yupPhs1_XLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > ydownPhs1_XLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > yupPhs2_XLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > ydownPhs2_XLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_XLOW.
+  Matrix< Array<dcomplex> > yupPhs1_YLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_YLOW.
+  Matrix< Array<dcomplex> > ydownPhs1_YLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_YLOW.
+  Matrix< Array<dcomplex> > yupPhs2_YLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_YLOW.
+  Matrix< Array<dcomplex> > ydownPhs2_YLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_YLOW.
 
   /*!
    * Shift a 2D field in Z. 
@@ -192,7 +191,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 REGION region=RGN_NOX);
+  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
@@ -211,7 +210,7 @@ private:
    * @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 std::vector<dcomplex> &phs, BoutReal *out);
+  void shiftZ(const BoutReal *in, const Array<dcomplex> &phs, BoutReal *out);
 
   /// Write out ParallelTransform variables to file
   void outputVars(Datafile &file);
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 62817e2f84..fd522e64db 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -43,14 +43,14 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
 
   int nmodes = mesh.LocalNz/2 + 1;
   //Allocate storage for complex intermediate
-  cmplx.resize(nmodes);
+  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
-ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -61,20 +61,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      fromAlignedPhs_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        fromAlignedPhs_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          fromAlignedPhs_CENTRE[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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(jx,jy)) , -sin(kwave*zShift(jx,jy)));
+            fromAlignedPhs_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*zShift(jx,jy)) , -sin(kwave*zShift(jx,jy)));
           }
         }
       }
@@ -84,7 +81,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -92,20 +89,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      fromAlignedPhs_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        fromAlignedPhs_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          fromAlignedPhs_XLOW[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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)));
+            fromAlignedPhs_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_XLOW(jx,jy)), -sin(kwave*zShift_XLOW(jx,jy)));
           }
         }
       }
@@ -115,7 +109,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -123,20 +117,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      fromAlignedPhs_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        fromAlignedPhs_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          fromAlignedPhs_YLOW[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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)));
+            fromAlignedPhs_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_YLOW(jx,jy)), -sin(kwave*zShift_YLOW(jx,jy)));
           }
         }
       }
@@ -158,7 +149,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
   };
 }
 
-ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -169,20 +160,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      toAlignedPhs_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        toAlignedPhs_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          toAlignedPhs_CENTRE[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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(jx,jy)), sin(kwave*zShift(jx,jy)));
+            toAlignedPhs_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*zShift(jx,jy)), sin(kwave*zShift(jx,jy)));
           }
         }
       }
@@ -192,7 +180,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -200,20 +188,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      toAlignedPhs_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        toAlignedPhs_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          toAlignedPhs_XLOW[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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)));
+            toAlignedPhs_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_XLOW(jx,jy)), sin(kwave*zShift_XLOW(jx,jy)));
           }
         }
       }
@@ -223,7 +208,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -231,20 +216,17 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      toAlignedPhs_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        toAlignedPhs_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          toAlignedPhs_YLOW[jx][jy].resize(nmodes);
-        }
+      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=0;jx<mesh.LocalNx;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)));
+            toAlignedPhs_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*zShift_YLOW(jx,jy)), sin(kwave*zShift_YLOW(jx,jy)));
           }
         }
       }
@@ -266,7 +248,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
   };
 }
 
-ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -277,21 +259,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      yupPhs1_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs1_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs1_CENTRE[jx][jy].resize(nmodes);
-        }
+      yupPhs1_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs1_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs1_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+            yupPhs1_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
           }
         }
       }
@@ -301,7 +280,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -309,21 +288,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      yupPhs1_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs1_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs1_XLOW[jx][jy].resize(nmodes);
-        }
+      yupPhs1_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs1_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs1_XLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+            yupPhs1_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
           }
         }
       }
@@ -333,7 +309,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -341,21 +317,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      yupPhs1_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs1_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs1_YLOW[jx][jy].resize(nmodes);
-        }
+      yupPhs1_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs1_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs1_YLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
+            yupPhs1_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*yupShift1) , -sin(kwave*yupShift1));
           }
         }
       }
@@ -377,7 +350,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs1(CELL_LOC location) {
   };
 }
 
-ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -388,21 +361,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      yupPhs2_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs2_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs2_CENTRE[jx][jy].resize(nmodes);
-        }
+      yupPhs2_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs2_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs2_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+            yupPhs2_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
           }
         }
       }
@@ -412,7 +382,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -420,21 +390,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      yupPhs2_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs2_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs2_XLOW[jx][jy].resize(nmodes);
-        }
+      yupPhs2_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs2_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs2_XLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+            yupPhs2_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
           }
         }
       }
@@ -444,7 +411,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -452,21 +419,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      yupPhs2_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        yupPhs2_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          yupPhs2_YLOW[jx][jy].resize(nmodes);
-        }
+      yupPhs2_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : yupPhs2_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal yupShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            yupPhs2_YLOW[jx][jy][jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
+            yupPhs2_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*yupShift2) , -sin(kwave*yupShift2));
           }
         }
       }
@@ -476,7 +440,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   }
   case CELL_ZLOW: {
     // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
-    return getYupPhs1(CELL_CENTRE);
+    return getYupPhs2(CELL_CENTRE);
     break;
   }
   default: {
@@ -488,7 +452,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYupPhs2(CELL_LOC location) {
   };
 }
 
-ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -499,21 +463,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      ydownPhs1_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs1_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs1_CENTRE[jx][jy].resize(nmodes);
-        }
+      ydownPhs1_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs1_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs1_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+            ydownPhs1_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
           }
         }
       }
@@ -523,7 +484,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -531,21 +492,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      ydownPhs1_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs1_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs1_XLOW[jx][jy].resize(nmodes);
-        }
+      ydownPhs1_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs1_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs1_XLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+            ydownPhs1_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
           }
         }
       }
@@ -555,7 +513,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -563,21 +521,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      ydownPhs1_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs1_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs1_YLOW[jx][jy].resize(nmodes);
-        }
+      ydownPhs1_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs1_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs1_YLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
+            ydownPhs1_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift1) , -sin(kwave*ydownShift1));
           }
         }
       }
@@ -599,7 +554,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs1(CELL_LOC location) {
   };
 }
 
-ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
+Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   // bools so we only calculate the cached values the first time for each location
   static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
 
@@ -610,21 +565,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       BoutReal zlength = mesh.coordinates()->zlength();
 
       first_CENTRE = false;
-      ydownPhs2_CENTRE.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs2_CENTRE[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs2_CENTRE[jx][jy].resize(nmodes);
-        }
+      ydownPhs2_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs2_CENTRE) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs2_CENTRE[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+            ydownPhs2_CENTRE(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
           }
         }
       }
@@ -634,7 +586,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   }
   case CELL_XLOW: {
     if (first_XLOW) {
-      ASSERT1(mesh.xstart>2); //otherwise we cannot interpolate in the x-direction
+      ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -642,21 +594,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
       first_XLOW = false;
-      ydownPhs2_XLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs2_XLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs2_XLOW[jx][jy].resize(nmodes);
-        }
+      ydownPhs2_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs2_XLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs2_XLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+            ydownPhs2_XLOW(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
           }
         }
       }
@@ -666,7 +615,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   }
   case CELL_YLOW: {
     if (first_YLOW) {
-      ASSERT1(mesh.ystart>2); //otherwise we cannot interpolate in the y-direction
+      ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
@@ -674,21 +623,18 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
       first_YLOW = false;
-      ydownPhs2_YLOW.resize(mesh.LocalNx);
-
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        ydownPhs2_YLOW[jx].resize(mesh.LocalNy);
-        for(int jy=0;jy<mesh.LocalNy;jy++){
-          ydownPhs2_YLOW[jx][jy].resize(nmodes);
-        }
+      ydownPhs2_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
+      for (auto &element : ydownPhs2_YLOW) {
+        element = Array<dcomplex>(mesh.LocalNz);
       }
+
       //To/From field aligned phases
       for(int jx=0;jx<mesh.LocalNx;jx++){
         for(int jy=0;jy<mesh.LocalNy;jy++){
           BoutReal ydownShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
-            ydownPhs2_YLOW[jx][jy][jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
+            ydownPhs2_YLOW(jx, jy)[jz] = dcomplex(cos(kwave*ydownShift2) , -sin(kwave*ydownShift2));
           }
         }
       }
@@ -698,7 +644,7 @@ ShiftedMetric::arr3Dvec ShiftedMetric::getYdownPhs2(CELL_LOC location) {
   }
   case CELL_ZLOW: {
     // shifts don't depend on z, so are the same for CELL_ZLOW as for CELL_CENTRE
-    return getYdownPhs1(CELL_CENTRE);
+    return getYdownPhs2(CELL_CENTRE);
     break;
   }
   default: {
@@ -719,10 +665,10 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
   
   Field3D& yup1 = f.yup();
   yup1.allocate();
-  arr3Dvec phases = getYupPhs1(location);
+  Matrix< Array<dcomplex> > phases = getYupPhs1(location);
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy+1,0)), phases[jx][jy], &(yup1(jx,jy+1,0)));
+      shiftZ(&(f(jx,jy+1,0)), phases(jx, jy), &(yup1(jx,jy+1,0)));
     }
   }
   if (mesh.ystart>1) {
@@ -731,7 +677,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
     phases = getYupPhs2(location);
     for(int jx=0;jx<mesh.LocalNx;jx++) {
       for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-        shiftZ(&(f(jx,jy+2,0)), phases[jx][jy], &(yup2(jx,jy+2,0)));
+        shiftZ(&(f(jx,jy+2,0)), phases(jx, jy), &(yup2(jx,jy+2,0)));
       }
     }
   }
@@ -741,7 +687,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
   phases = getYdownPhs1(location);
   for(int jx=0;jx<mesh.LocalNx;jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy-1,0)), phases[jx][jy], &(ydown1(jx,jy-1,0)));
+      shiftZ(&(f(jx,jy-1,0)), phases(jx, jy), &(ydown1(jx,jy-1,0)));
     }
   }
   if (mesh.ystart > 1) {
@@ -750,7 +696,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
     phases = getYdownPhs2(location);
     for(int jx=0;jx<mesh.LocalNx;jx++) {
       for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-        shiftZ(&(f(jx,jy-2,0)), phases[jx][jy], &(ydown2(jx,jy-2,0)));
+        shiftZ(&(f(jx,jy-2,0)), phases(jx, jy), &(ydown2(jx,jy-2,0)));
       }
     }
   }
@@ -772,7 +718,7 @@ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION reg
   return shiftZ(f, getFromAlignedPhs(f.getLocation()), region);
 }
 
-const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region) {
+const Field3D ShiftedMetric::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)
@@ -781,16 +727,16 @@ const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const
   Field3D result(f); // Initialize from f, mostly so location get set correctly. (Does not copy data because of copy-on-change).
 
   for(auto i : f.region2D(region)) {
-    shiftZ(f(i.x,i.y), phs[i.x][i.y], result(i.x,i.y));
+    shiftZ(f(i.x,i.y), phs(i.x, i.y), result(i.x,i.y));
   }
   
   return result;
 
 }
 
-void ShiftedMetric::shiftZ(const BoutReal *in, const std::vector<dcomplex> &phs, BoutReal *out) {
+void ShiftedMetric::shiftZ(const BoutReal *in, const Array<dcomplex> &phs, BoutReal *out) {
   // Take forward FFT
-  rfft(in, mesh.LocalNz, &cmplx[0]);
+  rfft(in, mesh.LocalNz, cmplx.begin());
 
   //Following is an algorithm approach to write a = a*b where a and b are
   //vectors of dcomplex.
@@ -802,7 +748,7 @@ void ShiftedMetric::shiftZ(const BoutReal *in, const std::vector<dcomplex> &phs,
     cmplx[jz] *= phs[jz];
   }
 
-  irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT
+  irfft(cmplx.begin(), mesh.LocalNz, out); // Reverse FFT
 }
 
 //Old approach retained so we can still specify a general zShift
@@ -834,10 +780,10 @@ void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutRe
   int nmodes = len/2 + 1;
 
   // Complex array used for FFTs
-  cmplxLoc.resize(nmodes);
+  cmplxLoc = Array<dcomplex>(nmodes);
   
   // Take forward FFT
-  rfft(in, len, &cmplxLoc[0]);
+  rfft(in, len, cmplxLoc.begin());
   
   // Apply phase shift
   BoutReal zlength = mesh.coordinates()->zlength();
@@ -846,7 +792,7 @@ void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutRe
     cmplxLoc[jz] *= dcomplex(cos(kwave*zangle) , -sin(kwave*zangle));
   }
 
-  irfft(&cmplxLoc[0], len, out); // Reverse FFT
+  irfft(cmplxLoc.begin(), len, out); // Reverse FFT
 }
 
 void ShiftedMetric::outputVars(Datafile &file) {

From f30daa92cab7448973d432d4dc7f60f03c7f58d3 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 16 May 2018 19:13:21 +0100
Subject: [PATCH 18/21] Make loops go from ystart->yend, not 0->Ny-1

---
 src/mesh/parallel/shiftedmetric.cxx | 56 ++++++++++++++---------------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index fd522e64db..542081d0b6 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -265,8 +265,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -294,8 +294,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -323,8 +323,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -367,8 +367,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -396,8 +396,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -425,8 +425,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -469,8 +469,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -498,8 +498,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -527,8 +527,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -571,8 +571,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -600,8 +600,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -629,8 +629,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      for(int jx=0; jx<mesh.LocalNx; jx++) {
+        for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
             BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
@@ -666,7 +666,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
   Field3D& yup1 = f.yup();
   yup1.allocate();
   Matrix< Array<dcomplex> > phases = getYupPhs1(location);
-  for(int jx=0;jx<mesh.LocalNx;jx++) {
+  for(int jx=0; jx<mesh.LocalNx; jx++) {
     for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
       shiftZ(&(f(jx,jy+1,0)), phases(jx, jy), &(yup1(jx,jy+1,0)));
     }
@@ -675,7 +675,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
     Field3D& yup2 = f.yup(2);
     yup2.allocate();
     phases = getYupPhs2(location);
-    for(int jx=0;jx<mesh.LocalNx;jx++) {
+    for(int jx=0; jx<mesh.LocalNx; jx++) {
       for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
         shiftZ(&(f(jx,jy+2,0)), phases(jx, jy), &(yup2(jx,jy+2,0)));
       }
@@ -686,7 +686,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
   ydown1.allocate();
   phases = getYdownPhs1(location);
   for(int jx=0;jx<mesh.LocalNx;jx++) {
-    for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
+    for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
       shiftZ(&(f(jx,jy-1,0)), phases(jx, jy), &(ydown1(jx,jy-1,0)));
     }
   }
@@ -695,7 +695,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
     ydown2.allocate();
     phases = getYdownPhs2(location);
     for(int jx=0;jx<mesh.LocalNx;jx++) {
-      for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
+      for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
         shiftZ(&(f(jx,jy-2,0)), phases(jx, jy), &(ydown2(jx,jy-2,0)));
       }
     }

From 582d3375a5ce98c511118881a03b18fee87c1d2d Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 21 May 2018 21:08:53 +0100
Subject: [PATCH 19/21] Ensure all fields have correct location

In ParallelTransform classes and in interp_to. This is necessary to make
sure Flexible<> fields know the right location to give when they are
given as arguments to arithmetic operators.
---
 src/field/field3d.cxx      | 4 ++++
 src/mesh/interpolation.cxx | 1 +
 2 files changed, 5 insertions(+)

diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 6435316388..e76d28a353 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -190,10 +190,14 @@ void Field3D::splitYupYdown() {
 
   // yup_array and ydown_array null
   yup1_field = new Field3D(fieldmesh);
+  yup1_field->setLocation(location);
   ydown1_field = new Field3D(fieldmesh);
+  ydown1_field->setLocation(location);
   if (fieldmesh->ystart>1) {
     yup2_field = new Field3D(fieldmesh);
+    yup2_field->setLocation(location);
     ydown2_field = new Field3D(fieldmesh);
+    ydown2_field->setLocation(location);
   }
 }
 
diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx
index 93baef9ad5..684316c6e8 100644
--- a/src/mesh/interpolation.cxx
+++ b/src/mesh/interpolation.cxx
@@ -56,6 +56,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
 
   Mesh *fieldmesh = var.getMesh();
   Field3D result(fieldmesh);
+  result.setLocation(loc);
 
   if ((loc != CELL_CENTRE && loc != CELL_DEFAULT) && (fieldmesh->StaggerGrids == false)) {
     throw BoutException("Asked to interpolate, but StaggerGrids is disabled!");

From ad81d4a272e1baacff52289f0201821afddc91d0 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 18 Jul 2018 11:09:49 +0100
Subject: [PATCH 20/21] Fix use of unassigned values in ShiftedMetric

To ensure corner guard cells are not read, restrict all loops in
ShiftedMetric to not include x-guard cells.

Also add some extra calls to invalidateGuards to help catch future
errors.
---
 src/mesh/parallel/shiftedmetric.cxx | 88 ++++++++++++++++-------------
 1 file changed, 48 insertions(+), 40 deletions(-)

diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index 542081d0b6..b831bc4a6b 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -67,8 +67,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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(jx,jy)) , -sin(kwave*zShift(jx,jy)));
@@ -95,8 +95,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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)));
@@ -123,8 +123,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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)));
@@ -166,8 +166,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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(jx,jy)), sin(kwave*zShift(jx,jy)));
@@ -194,8 +194,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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)));
@@ -222,8 +222,8 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0;jx<mesh.LocalNx;jx++){
-        for(int jy=0;jy<mesh.LocalNy;jy++){
+      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)));
@@ -265,7 +265,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -294,7 +294,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -323,7 +323,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -367,7 +367,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -396,7 +396,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -425,7 +425,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal yupShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -469,7 +469,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -498,7 +498,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -527,7 +527,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-1);
           for(int jz=0;jz<nmodes;jz++) {
@@ -571,7 +571,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -600,7 +600,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -629,7 +629,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       }
 
       //To/From field aligned phases
-      for(int jx=0; jx<mesh.LocalNx; jx++) {
+      for(int jx=mesh.xstart; jx<=mesh.xend; jx++) {
         for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
           BoutReal ydownShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-2);
           for(int jz=0;jz<nmodes;jz++) {
@@ -660,44 +660,52 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
  * Calculate the Y up and down fields
  */
 void ShiftedMetric::calcYUpDown(Field3D &f) {
+  ASSERT1(&mesh == f.getMesh());
   f.splitYupYdown();
   CELL_LOC location = 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.)
+  // Also, only makes sense to calculate yup/ydown for the y-grid points (not
+  // guard cells) since, e.g. yup(yend) contains the (shifted) value for f in
+  // the 'guard cell'.
+  // Therefore, only loop over RGN_NOBNDRY here.
+
   Field3D& yup1 = f.yup();
   yup1.allocate();
+  invalidateGuards(yup1); // Won't set x-guard cells, so allow checking to throw exception if they are used.
   Matrix< Array<dcomplex> > phases = getYupPhs1(location);
-  for(int jx=0; jx<mesh.LocalNx; jx++) {
-    for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-      shiftZ(&(f(jx,jy+1,0)), phases(jx, jy), &(yup1(jx,jy+1,0)));
-    }
+  for(auto i : f.region2D(RGN_NOBNDRY)) {
+    shiftZ(f(i.x, i.y+1), phases(i.x, i.y), yup1(i.x, i.y+1));
   }
   if (mesh.ystart>1) {
     Field3D& yup2 = f.yup(2);
     yup2.allocate();
+    invalidateGuards(yup2); // Won't set x-guard cells, so allow checking to throw exception if they are used.
     phases = getYupPhs2(location);
-    for(int jx=0; jx<mesh.LocalNx; jx++) {
-      for(int jy=mesh.ystart;jy<=mesh.yend;jy++) {
-        shiftZ(&(f(jx,jy+2,0)), phases(jx, jy), &(yup2(jx,jy+2,0)));
-      }
+    for(auto i : f.region2D(RGN_NOBNDRY)) {
+      shiftZ(f(i.x, i.y+2), phases(i.x, i.y), yup2(i.x, i.y+2));
     }
   }
 
   Field3D& ydown1 = f.ydown();
   ydown1.allocate();
+  invalidateGuards(ydown1); // Won't set x-guard cells, so allow checking to throw exception if they are used.
   phases = getYdownPhs1(location);
-  for(int jx=0;jx<mesh.LocalNx;jx++) {
-    for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
-      shiftZ(&(f(jx,jy-1,0)), phases(jx, jy), &(ydown1(jx,jy-1,0)));
-    }
+  for(auto i : f.region2D(RGN_NOBNDRY)) {
+    shiftZ(f(i.x, i.y-1), phases(i.x, i.y), ydown1(i.x, i.y-1));
   }
   if (mesh.ystart > 1) {
     Field3D& ydown2 = f.ydown(2);
     ydown2.allocate();
+    invalidateGuards(ydown2); // Won't set x-guard cells, so allow checking to throw exception if they are used.
     phases = getYdownPhs2(location);
-    for(int jx=0;jx<mesh.LocalNx;jx++) {
-      for(int jy=mesh.ystart; jy<=mesh.yend; jy++) {
-        shiftZ(&(f(jx,jy-2,0)), phases(jx, jy), &(ydown2(jx,jy-2,0)));
-      }
+    for(auto i : f.region2D(RGN_NOBNDRY)) {
+      shiftZ(f(i.x, i.y-2), phases(i.x, i.y), ydown2(i.x, i.y-2));
     }
   }
 }

From 970d72df125e8d5f60d12e421f43fe634f7d1205 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 18 Jul 2018 17:05:14 +0100
Subject: [PATCH 21/21] Check if phases have been calculated with member bools
 not static bool

Using static bool to check for cached phases means the check is global
so if two ShiftedMetric objects are created the one which calls 'get...'
second will fail. This is unlikely in a normal simulation but might
happen if there are multiple meshes, for example in tests or
post-processing using boutcore.
---
 include/bout/paralleltransform.hxx  |  7 ++
 src/mesh/parallel/shiftedmetric.cxx | 99 +++++++++++++----------------
 2 files changed, 51 insertions(+), 55 deletions(-)

diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx
index 2c71e4f0ff..03d280454b 100644
--- a/include/bout/paralleltransform.hxx
+++ b/include/bout/paralleltransform.hxx
@@ -147,6 +147,13 @@ private:
   Matrix< Array<dcomplex> > getYupPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
   Matrix< Array<dcomplex> > getYdownPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary;
 
+  bool has_toAligned_CENTRE, has_toAligned_XLOW, has_toAligned_YLOW; ///< Have phase shifts for shift to field aligned coordinates been calculated
+  bool has_fromAligned_CENTRE, has_fromAligned_XLOW, has_fromAligned_YLOW; ///< Have phase shifts for shift from field aligned coordinates been calculated
+  bool has_yupPhs1_CENTRE, has_yupPhs1_XLOW, has_yupPhs1_YLOW; ///< Have phase shifts for yup1 been calculated
+  bool has_ydownPhs1_CENTRE, has_ydownPhs1_XLOW, has_ydownPhs1_YLOW; ///< Have phase shifts for ydown1 been calculated
+  bool has_yupPhs2_CENTRE, has_yupPhs2_XLOW, has_yupPhs2_YLOW; ///< Have phase shifts for yup2 been calculated
+  bool has_ydownPhs2_CENTRE, has_ydownPhs2_XLOW, has_ydownPhs2_YLOW; ///< Have phase shifts for ydown2 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.
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index b831bc4a6b..cb998d38c7 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -16,7 +16,14 @@
 
 #include <output.hxx>
 
-ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
+ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&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),
+  has_yupPhs1_CENTRE(false), has_yupPhs1_XLOW(false), has_yupPhs1_YLOW(false),
+  has_ydownPhs1_CENTRE(false), has_ydownPhs1_XLOW(false), has_ydownPhs1_YLOW(false),
+  has_yupPhs2_CENTRE(false), has_yupPhs2_XLOW(false), has_yupPhs2_YLOW(false),
+  has_ydownPhs2_CENTRE(false), has_ydownPhs2_XLOW(false), has_ydownPhs2_YLOW(false)
+{
   // Read the zShift angle from the mesh
   
   if(mesh.get(zShift, "zShift")) {
@@ -51,16 +58,13 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
 //once we've been created so cache the complex phases used in transformations
 //the first time they are needed
 Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (!has_fromAligned_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_fromAligned_CENTRE = true;
       fromAlignedPhs_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : fromAlignedPhs_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -80,7 +84,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_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.coordinates()->zlength();
@@ -88,7 +92,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_fromAligned_XLOW = true;
       fromAlignedPhs_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : fromAlignedPhs_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -108,7 +112,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_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.coordinates()->zlength();
@@ -116,7 +120,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_fromAligned_YLOW = true;
       fromAlignedPhs_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : fromAlignedPhs_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -150,16 +154,13 @@ Matrix< Array<dcomplex> > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) {
 }
 
 Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (!has_toAligned_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_toAligned_CENTRE = true;
       toAlignedPhs_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : toAlignedPhs_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -179,7 +180,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_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.coordinates()->zlength();
@@ -187,7 +188,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_toAligned_XLOW = true;
       toAlignedPhs_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : toAlignedPhs_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -207,7 +208,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_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.coordinates()->zlength();
@@ -215,7 +216,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_toAligned_YLOW = true;
       toAlignedPhs_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : toAlignedPhs_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -249,16 +250,13 @@ Matrix< Array<dcomplex> > ShiftedMetric::getToAlignedPhs(CELL_LOC location) {
 }
 
 Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (!has_yupPhs1_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_yupPhs1_CENTRE = true;
       yupPhs1_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs1_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -279,7 +277,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_XLOW) {
+    if (!has_yupPhs1_XLOW) {
       ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -287,7 +285,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_yupPhs1_XLOW = true;
       yupPhs1_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs1_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -308,7 +306,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_YLOW) {
+    if (!has_yupPhs1_YLOW) {
       ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -316,7 +314,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_yupPhs1_YLOW = true;
       yupPhs1_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs1_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -351,16 +349,13 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs1(CELL_LOC location) {
 }
 
 Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (!has_yupPhs2_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_yupPhs2_CENTRE = true;
       yupPhs2_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs2_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -381,7 +376,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_XLOW) {
+    if (!has_yupPhs2_XLOW) {
       ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -389,7 +384,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_yupPhs2_XLOW = true;
       yupPhs2_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs2_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -410,7 +405,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_YLOW) {
+    if (!has_yupPhs2_YLOW) {
       ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -418,7 +413,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_yupPhs2_YLOW = true;
       yupPhs2_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : yupPhs2_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -453,16 +448,13 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYupPhs2(CELL_LOC location) {
 }
 
 Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (has_ydownPhs1_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_ydownPhs1_CENTRE = true;
       ydownPhs1_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs1_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -483,7 +475,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_XLOW) {
+    if (!has_ydownPhs1_XLOW) {
       ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -491,7 +483,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_ydownPhs1_XLOW = true;
       ydownPhs1_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs1_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -512,7 +504,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_YLOW) {
+    if (!has_ydownPhs1_YLOW) {
       ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -520,7 +512,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_ydownPhs1_YLOW = true;
       ydownPhs1_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs1_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -555,16 +547,13 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs1(CELL_LOC location) {
 }
 
 Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
-  // bools so we only calculate the cached values the first time for each location
-  static bool first_CENTRE = true, first_XLOW=true, first_YLOW=true;
-
   switch (location) {
   case CELL_CENTRE: {
-    if (first_CENTRE) {
+    if (!has_ydownPhs2_CENTRE) {
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
 
-      first_CENTRE = false;
+      has_ydownPhs2_CENTRE = true;
       ydownPhs2_CENTRE = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs2_CENTRE) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -585,7 +574,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
     break;
   }
   case CELL_XLOW: {
-    if (first_XLOW) {
+    if (!has_ydownPhs2_XLOW) {
       ASSERT1(mesh.xstart>=2); //otherwise we cannot interpolate in the x-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -593,7 +582,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       // get zShift at XLOW
       Field2D zShift_XLOW = zShift.get(CELL_XLOW);
 
-      first_XLOW = false;
+      has_ydownPhs2_XLOW = true;
       ydownPhs2_XLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs2_XLOW) {
         element = Array<dcomplex>(mesh.LocalNz);
@@ -614,7 +603,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
     break;
   }
   case CELL_YLOW: {
-    if (first_YLOW) {
+    if (!has_ydownPhs2_YLOW) {
       ASSERT1(mesh.ystart>=2); //otherwise we cannot interpolate in the y-direction
       int nmodes = mesh.LocalNz/2 + 1;
       BoutReal zlength = mesh.coordinates()->zlength();
@@ -622,7 +611,7 @@ Matrix< Array<dcomplex> > ShiftedMetric::getYdownPhs2(CELL_LOC location) {
       // get zShift at YLOW
       Field2D zShift_YLOW = zShift.get(CELL_YLOW);
 
-      first_YLOW = false;
+      has_ydownPhs2_YLOW = true;
       ydownPhs2_YLOW = Matrix< Array<dcomplex> >(mesh.LocalNx, mesh.LocalNy);
       for (auto &element : ydownPhs2_YLOW) {
         element = Array<dcomplex>(mesh.LocalNz);