diff --git a/include/bout/field.hxx b/include/bout/field.hxx
index 27835ecbd7..77a24ea96a 100644
--- a/include/bout/field.hxx
+++ b/include/bout/field.hxx
@@ -133,6 +133,7 @@ public:
   virtual void setRegion(const std::string& UNUSED(region_name)) {}
   virtual void resetRegion() {}
   virtual std::optional<size_t> getRegionID() const { return {}; }
+  virtual int getYoffset() const { return 0; }
 
 private:
   /// Labels for the type of coordinate system this field is defined over
@@ -186,7 +187,7 @@ template <typename T>
 inline T emptyFrom(const T& f) {
   static_assert(bout::utils::is_Field_v<T>, "emptyFrom only works on Fields");
   return T(f.getMesh(), f.getLocation(), {f.getDirectionY(), f.getDirectionZ()},
-           f.getRegionID())
+           f.getRegionID(), f.getYoffset())
       .allocate();
 }
 
diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx
index 5eab330e8e..f51424f5bd 100644
--- a/include/bout/field2d.hxx
+++ b/include/bout/field2d.hxx
@@ -67,7 +67,7 @@ public:
   Field2D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE,
           DirectionTypes directions_in = {YDirectionType::Standard,
                                           ZDirectionType::Average},
-          std::optional<size_t> region = {});
+          std::optional<size_t> region = {}, int yoffset = 0);
 
   /*!
    * Copy constructor. After this both fields
diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 0643efbe0a..1a70fd32f1 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -29,6 +29,7 @@ class Field3D;
 #include "bout/array.hxx"
 #include "bout/assert.hxx"
 #include "bout/bout_types.hxx"
+#include "bout/boutexception.hxx"
 #include "bout/field.hxx"
 #include "bout/field2d.hxx"
 #include "bout/fieldperp.hxx"
@@ -167,7 +168,7 @@ public:
   Field3D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE,
           DirectionTypes directions_in = {YDirectionType::Standard,
                                           ZDirectionType::Standard},
-          std::optional<size_t> regionID = {});
+          std::optional<size_t> regionID = {}, int yoffset = 0);
 
   /*!
    * Copy constructor
@@ -271,23 +272,27 @@ public:
   /// Return reference to yup field
   Field3D& yup(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(yup_fields[index].yoffset == static_cast<int>(index) + 1);
     return yup_fields[index];
   }
   /// Return const reference to yup field
   const Field3D& yup(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(yup_fields[index].yoffset == static_cast<int>(index) + 1);
     return yup_fields[index];
   }
 
   /// Return reference to ydown field
   Field3D& ydown(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(ydown_fields[index].yoffset == -1 - static_cast<int>(index));
     return ydown_fields[index];
   }
 
   /// Return const reference to ydown field
   const Field3D& ydown(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(ydown_fields[index].yoffset == -1 - static_cast<int>(index));
     return ydown_fields[index];
   }
 
@@ -339,6 +344,8 @@ public:
   void setRegion(std::optional<size_t> id) override;
   std::optional<size_t> getRegionID() const override { return regionID; };
 
+  int getYoffset() const override { return yoffset; };
+
   /// Return a Region<Ind2D> reference to use to iterate over the x- and
   /// y-indices of this field
   const Region<Ind2D>& getRegion2D(REGION region) const;
@@ -351,8 +358,27 @@ public:
     return std::end(getRegion("RGN_ALL"));
   };
 
-  BoutReal& operator[](const Ind3D& d) { return data[d.ind]; }
-  const BoutReal& operator[](const Ind3D& d) const { return data[d.ind]; }
+  BoutReal& operator[](const Ind3D& d) {
+    if (d.yoffset != 0) {
+      if (yoffset == 0) {
+        if (hasParallelSlices()) {
+          return ynext(d.yoffset)[d];
+        }
+#if CHECK >= 2
+        if (isFci()) { // We probably should assert here that this is field aligned
+          throw BoutException(
+              "Tried to access parallel slices, but they are not calculated!");
+        }
+#endif
+      } else {
+        ASSERT2(d.yoffset == yoffset);
+      }
+    }
+    return data[d.ind];
+  }
+  const BoutReal& operator[](const Ind3D& d) const {
+    return (*const_cast<Field3D*>(this))[d];
+  }
 
   BoutReal& operator()(const IndPerp& d, int jy);
   const BoutReal& operator()(const IndPerp& d, int jy) const;
@@ -547,6 +573,10 @@ private:
   template <class T>
   Options* track(const T& change, std::string operation);
   Options* track(const BoutReal& change, std::string operation);
+
+public:
+  int yoffset{0};
+  void setParallelRegions();
 };
 
 // Non-member overloaded operators
diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx
index b50eef1991..8cb649a29d 100644
--- a/include/bout/fieldperp.hxx
+++ b/include/bout/fieldperp.hxx
@@ -59,7 +59,7 @@ public:
             int yindex_in = -1,
             DirectionTypes directions_in = {YDirectionType::Standard,
                                             ZDirectionType::Standard},
-            std::optional<size_t> regionID = {});
+            std::optional<size_t> regionID = {}, int yoffset = 0);
 
   /*!
    * Copy constructor. After this the data
diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index f441b3edd7..e7750747ed 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -168,9 +168,12 @@ template <IND_TYPE N>
 struct SpecificInd {
   int ind = -1;         ///< 1D index into Field
   int ny = -1, nz = -1; ///< Sizes of y and z dimensions
+  int yoffset = 0;
 
   SpecificInd() = default;
   SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){};
+  SpecificInd(int i, int ny, int nz, int yoffset)
+      : ind(i), ny(ny), nz(nz), yoffset(yoffset){};
   explicit SpecificInd(int i) : ind(i){};
 
   /// Allow explicit conversion to an int
@@ -275,7 +278,9 @@ struct SpecificInd {
     }
   }
 
-  inline SpecificInd xp(int dx = 1) const { return {ind + (dx * ny * nz), ny, nz}; }
+  inline SpecificInd xp(int dx = 1) const {
+    return {ind + (dx * ny * nz), ny, nz, yoffset};
+  }
   /// The index one point -1 in x
   inline SpecificInd xm(int dx = 1) const { return xp(-dx); }
   /// The index one point +1 in y
@@ -286,8 +291,18 @@ struct SpecificInd {
     }
 #endif
     ASSERT3(std::abs(dy) < ny);
-    return {ind + (dy * nz), ny, nz};
+    return {ind + (dy * nz), ny, nz, yoffset + dy};
   }
+  inline SpecificInd yp_no_parallel_shift(int dy = 1) const {
+#if CHECK >= 4
+    if (y() + dy < 0 or y() + dy >= ny) {
+      throw BoutException("Offset in y ({:d}) would go out of bounds at {:d}", dy, ind);
+    }
+#endif
+    ASSERT3(std::abs(dy) < ny);
+    return {ind + (dy * nz), ny, nz, yoffset};
+  }
+
   /// The index one point -1 in y
   inline SpecificInd ym(int dy = 1) const { return yp(-dy); }
   /// The index one point +1 in z. Wraps around zend to zstart
@@ -297,7 +312,7 @@ struct SpecificInd {
   inline SpecificInd zp(int dz = 1) const {
     ASSERT3(dz >= 0);
     dz = dz <= nz ? dz : dz % nz; //Fix in case dz > nz, if not force it to be in range
-    return {(ind + dz) % nz < dz ? ind - nz + dz : ind + dz, ny, nz};
+    return {(ind + dz) % nz < dz ? ind - nz + dz : ind + dz, ny, nz, yoffset};
   }
   /// The index one point -1 in z. Wraps around zstart to zend
   /// An alternative, non-branching calculation is :
@@ -306,7 +321,7 @@ struct SpecificInd {
   inline SpecificInd zm(int dz = 1) const {
     dz = dz <= nz ? dz : dz % nz; //Fix in case dz > nz, if not force it to be in range
     ASSERT3(dz >= 0);
-    return {(ind) % nz < dz ? ind + nz - dz : ind - dz, ny, nz};
+    return {(ind) % nz < dz ? ind + nz - dz : ind - dz, ny, nz, yoffset};
   }
   /// Automatically select zm or zp depending on sign
   inline SpecificInd zpm(int dz) const { return dz > 0 ? zp(dz) : zm(-dz); }
diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx
index 4591d228f7..e6674eb7e9 100644
--- a/src/field/field2d.cxx
+++ b/src/field/field2d.cxx
@@ -49,7 +49,7 @@
 #include <bout/assert.hxx>
 
 Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
-                 std::optional<size_t> UNUSED(regionID))
+                 std::optional<size_t> UNUSED(regionID), int UNUSED(yoffset))
     : Field(localmesh, location_in, directions_in) {
 
   if (fieldmesh) {
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 74a6f0853c..b66aa9c58d 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -49,8 +49,8 @@
 
 /// Constructor
 Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
-                 std::optional<size_t> regionID)
-    : Field(localmesh, location_in, directions_in), regionID{regionID} {
+                 std::optional<size_t> regionID, const int yoffset)
+    : Field(localmesh, location_in, directions_in), regionID{regionID}, yoffset{yoffset} {
 #if BOUT_USE_TRACK
   name = "<F3D>";
 #endif
@@ -66,7 +66,7 @@ Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes direction
 /// later)
 Field3D::Field3D(const Field3D& f)
     : Field(f), data(f.data), yup_fields(f.yup_fields), ydown_fields(f.ydown_fields),
-      regionID(f.regionID) {
+      regionID(f.regionID), yoffset(f.yoffset) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -157,7 +157,9 @@ void Field3D::splitParallelSlices() {
     // Note the fields constructed here will be fully overwritten by the
     // ParallelTransform, so we don't need a full constructor
     yup_fields.emplace_back(fieldmesh);
+    yup_fields[i].yoffset = i + 1;
     ydown_fields.emplace_back(fieldmesh);
+    ydown_fields[i].yoffset = -1 - i;
     if (isFci()) {
       yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1));
       ydown_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", -i - 1));
@@ -224,6 +226,15 @@ bool Field3D::requiresTwistShift(bool twist_shift_enabled) {
                                                                      getDirectionY());
 }
 
+void Field3D::setParallelRegions() {
+  for (int i = 0; i < fieldmesh->ystart; ++i) {
+    yup_fields[i].yoffset = i + 1;
+    yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1));
+    ydown_fields[i].yoffset = -1 - i;
+    ydown_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", -1 - i));
+  }
+}
+
 // Not in header because we need to access fieldmesh
 BoutReal& Field3D::operator()(const IndPerp& d, int jy) {
   return operator[](fieldmesh->indPerpto3D(d, jy));
@@ -277,6 +288,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) {
   yup_fields = rhs.yup_fields;
   ydown_fields = rhs.ydown_fields;
   regionID = rhs.regionID;
+  yoffset = rhs.yoffset;
 
   // Copy the data and data sizes
   nx = rhs.nx;
@@ -295,12 +307,13 @@ Field3D& Field3D::operator=(Field3D&& rhs) {
   // Move parallel slices or delete existing ones.
   yup_fields = std::move(rhs.yup_fields);
   ydown_fields = std::move(rhs.ydown_fields);
+  regionID = rhs.regionID;
+  yoffset = rhs.yoffset;
 
   // Move the data and data sizes
   nx = rhs.nx;
   ny = rhs.ny;
   nz = rhs.nz;
-  regionID = rhs.regionID;
 
   data = std::move(rhs.data);
 
@@ -896,6 +909,7 @@ void swap(Field3D& first, Field3D& second) noexcept {
   swap(first.deriv, second.deriv);
   swap(first.yup_fields, second.yup_fields);
   swap(first.ydown_fields, second.ydown_fields);
+  swap(first.yoffset, second.yoffset);
 }
 
 const Region<Ind3D>&
diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx
index c7a196a57a..a3458543c9 100644
--- a/src/field/fieldperp.cxx
+++ b/src/field/fieldperp.cxx
@@ -35,7 +35,8 @@
 #include <bout/utils.hxx>
 
 FieldPerp::FieldPerp(Mesh* localmesh, CELL_LOC location_in, int yindex_in,
-                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID))
+                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID),
+                     int UNUSED(yoffset))
     : Field(localmesh, location_in, directions), yindex(yindex_in) {
   if (fieldmesh) {
     nx = fieldmesh->LocalNx;
diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja
index 88e877c197..4cea18f4c8 100644
--- a/src/field/gen_fieldops.jinja
+++ b/src/field/gen_fieldops.jinja
@@ -12,6 +12,8 @@
     {% if lhs == rhs == "Field3D" %}
       {{out.name}}.setRegion({{lhs.name}}.getMesh()->getCommonRegion({{lhs.name}}.getRegionID(),
                                                                      {{rhs.name}}.getRegionID()));
+      ASSERT2({{lhs.name}}.yoffset == {{rhs.name}}.yoffset);
+      {{out.name}}.yoffset = {{lhs.name}}.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices() and {{rhs.name}}.hasParallelSlices()) {
         {{out.name}}.splitParallelSlices();
@@ -23,6 +25,7 @@
 #endif
     {% elif lhs == "Field3D" %}
       {{out.name}}.setRegion({{lhs.name}}.getRegionID());
+      {{out.name}}.yoffset = {{lhs.name}}.yoffset;
       {% if rhs == "BoutReal" %}
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices()) {
@@ -36,6 +39,7 @@
       {% endif %}
     {% elif rhs == "Field3D" %}
       {{out.name}}.setRegion({{rhs.name}}.getRegionID());
+      {{out.name}}.yoffset = {{rhs.name}}.yoffset;
       {% if lhs == "BoutReal" %}
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{rhs.name}}.isFci() and {{rhs.name}}.hasParallelSlices()) {
diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py
index 29631ff7aa..bf71a37671 100755
--- a/src/field/gen_fieldops.py
+++ b/src/field/gen_fieldops.py
@@ -65,6 +65,7 @@ def smart_open(filename, mode="r"):
 )
 
 header = """// This file is autogenerated - see gen_fieldops.py
+#include <bout/assert.hxx>
 #include <bout/field2d.hxx>
 #include <bout/field3d.hxx>
 #include <bout/globals.hxx>
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 75d2ede82d..eeb32225b4 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -1,4 +1,5 @@
 // This file is autogenerated - see gen_fieldops.py
+#include <bout/assert.hxx>
 #include <bout/field2d.hxx>
 #include <bout/field3d.hxx>
 #include <bout/globals.hxx>
@@ -15,6 +16,8 @@ Field3D operator*(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -87,6 +90,8 @@ Field3D operator/(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -159,6 +164,8 @@ Field3D operator+(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -231,6 +238,8 @@ Field3D operator-(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -303,6 +312,7 @@ Field3D operator*(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -364,6 +374,7 @@ Field3D operator/(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -427,6 +438,7 @@ Field3D operator+(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -488,6 +500,7 @@ Field3D operator-(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -640,6 +653,7 @@ Field3D operator*(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -708,6 +722,7 @@ Field3D operator/(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -778,6 +793,7 @@ Field3D operator+(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -846,6 +862,7 @@ Field3D operator-(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -915,6 +932,7 @@ Field3D operator*(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -941,6 +959,7 @@ Field3D operator/(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -967,6 +986,7 @@ Field3D operator+(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -993,6 +1013,7 @@ Field3D operator-(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -2209,6 +2230,7 @@ Field3D operator*(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2238,6 +2260,7 @@ Field3D operator/(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2267,6 +2290,7 @@ Field3D operator+(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2296,6 +2320,7 @@ Field3D operator-(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx
index aff546c2b0..aca0c319e8 100644
--- a/src/mesh/coordinates_accessor.cxx
+++ b/src/mesh/coordinates_accessor.cxx
@@ -55,9 +55,15 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) {
     COPY_STRIPE(J);
 
     data[stripe_size * ind.ind + static_cast<int>(Offset::B)] = coords->Bxy[ind];
-    data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] = coords->Bxy.yup()[ind];
-    data[stripe_size * ind.ind + static_cast<int>(Offset::Bydown)] =
-        coords->Bxy.ydown()[ind];
+    {
+      auto indc = ind;
+      indc.yoffset = 1;
+      data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] =
+          coords->Bxy.yup()[indc];
+      indc.yoffset = -1;
+      data[stripe_size * ind.ind + static_cast<int>(Offset::Bydown)] =
+          coords->Bxy.ydown()[indc];
+    }
 
     COPY_STRIPE(G1, G3);
     COPY_STRIPE(g11, g12, g13, g22, g23, g33);
diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx
index 03a062df42..394002424e 100644
--- a/src/mesh/interpolation/hermite_spline_xz.cxx
+++ b/src/mesh/interpolation/hermite_spline_xz.cxx
@@ -378,6 +378,7 @@ Field3D XZHermiteSplineBase<monotonic>::interpolate(const Field3D& f,
 
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
+  f_interp.yoffset = y_offset;
 
   const auto region2 = y_offset != 0 ? fmt::format("RGN_YPAR_{:+d}", y_offset) : region;
 
@@ -395,6 +396,7 @@ Field3D XZHermiteSplineBase<monotonic>::interpolate(const Field3D& f,
   MatMult(petscWeights, rhs, result);
   VecGetArrayRead(result, &cptr);
   BOUT_FOR(i, f.getRegion(region2)) {
+    i.yoffset = y_offset;
     f_interp[i] = cptr[int(i)];
     if constexpr (monotonic) {
       const auto iyp = i;
@@ -508,6 +510,7 @@ Field3D XZMonotonicHermiteSplineLegacy::interpolate(const Field3D& f,
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp(f.getMesh());
   f_interp.allocate();
+  f_interp.yoffset = y_offset;
 
   // Derivatives are used for tension and need to be on dimensionless
   // coordinates
diff --git a/src/mesh/interpolation/hermite_spline_z.cxx b/src/mesh/interpolation/hermite_spline_z.cxx
index c4c44cb7b4..02d8d45a61 100644
--- a/src/mesh/interpolation/hermite_spline_z.cxx
+++ b/src/mesh/interpolation/hermite_spline_z.cxx
@@ -151,6 +151,7 @@ Field3D ZHermiteSpline::interpolate(const Field3D& f,
 
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
+  f_interp.yoffset = y_offset;
 
   std::string local_fz_region;
   if (region_str == "DEFAULT") {
@@ -168,7 +169,7 @@ Field3D ZHermiteSpline::interpolate(const Field3D& f,
   Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", local_fz_region);
 
   BOUT_FOR(i, local_region) {
-    const auto corner = k_corner[i.ind].yp(y_offset);
+    const auto corner = k_corner[i.ind].yp_no_parallel_shift(y_offset);
     const auto corner_zp1 = corner.zp();
 
     // Interpolate in Z
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index cb6cb8410c..8cc204ff41 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -685,6 +685,15 @@ void Mesh::createDefaultRegions() {
                                 + getRegion3D("RGN_YGUARDS") + getRegion3D("RGN_ZGUARDS"))
                                    .unique());
 
+  // Add region for parallel slices
+  for (int i = 1; i <= ystart; ++i) {
+    for (const int offset : {-i, i}) {
+      const auto region = fmt::format("RGN_YPAR_{:+d}", offset);
+      addRegion3D(region, Region<Ind3D>(xstart, xend, ystart + offset, yend + offset, 0,
+                                        LocalNz - 1, LocalNy, LocalNz));
+    }
+  }
+
   //2D regions
   addRegion2D("RGN_ALL", Region<Ind2D>(0, LocalNx - 1, 0, LocalNy - 1, 0, 0, LocalNy, 1,
                                        maxregionblocksize));
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index 580897f47a..9d6997f660 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -99,6 +99,7 @@ bool load_parallel_metric_component(std::string name, Field3D& component, int of
   pcom.setRegion(fmt::format("RGN_YPAR_{:+d}", offset));
   pcom.name = name;
   BOUT_FOR(i, component.getRegion("RGN_NOBNDRY")) { pcom[i.yp(offset)] = tmp[i]; }
+  ASSERT2(pcom.yoffset == offset);
   return isValid;
 }
 #endif
@@ -339,15 +340,6 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
   region_no_boundary = region_no_boundary.mask(to_remove);
 
   interp->setRegion(region_no_boundary);
-
-  const auto region = fmt::format("RGN_YPAR_{:+d}", offset);
-  if (not map_mesh.hasRegion3D(region)) {
-    // The valid region for this slice
-    map_mesh.addRegion3D(
-        region, Region<Ind3D>(map_mesh.xstart, map_mesh.xend, map_mesh.ystart + offset,
-                              map_mesh.yend + offset, 0, map_mesh.LocalNz - 1,
-                              map_mesh.LocalNy, map_mesh.LocalNz));
-  }
 }
 
 Field3D FCIMap::integrate(Field3D& f) const {
@@ -430,8 +422,8 @@ void FCITransform::calcParallelSlices(Field3D& f) {
   // Interpolate f onto yup and ydown fields
   for (const auto& map : field_line_maps) {
     f.ynext(map.offset) = map.interpolate(f);
-    f.ynext(map.offset).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset));
   }
+  f.setParallelRegions();
 }
 
 void FCITransform::integrateParallelSlices(Field3D& f) {
@@ -449,6 +441,7 @@ void FCITransform::integrateParallelSlices(Field3D& f) {
   for (const auto& map : field_line_maps) {
     f.ynext(map.offset) = map.integrate(f);
   }
+  f.setParallelRegions();
 }
 
 void FCITransform::loadParallelMetrics(Coordinates* coords) {
diff --git a/src/mesh/parallel/identity.cxx b/src/mesh/parallel/identity.cxx
index 90fcaef45e..b9869a22eb 100644
--- a/src/mesh/parallel/identity.cxx
+++ b/src/mesh/parallel/identity.cxx
@@ -23,7 +23,9 @@ void ParallelTransformIdentity::calcParallelSlices(Field3D& f) {
   f.splitParallelSlices();
 
   for (int i = 0; i < f.getMesh()->ystart; ++i) {
+    f_copy.yoffset = i + 1;
     f.yup(i) = f_copy;
+    f_copy.yoffset = -1 - i;
     f.ydown(i) = f_copy;
   }
 }
diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx
index dfb397c626..bfe53b38e6 100644
--- a/src/mesh/parallel/shiftedmetricinterp.cxx
+++ b/src/mesh/parallel/shiftedmetricinterp.cxx
@@ -220,6 +220,7 @@ void ShiftedMetricInterp::calcParallelSlices(Field3D& f) {
   for (const auto& interp : parallel_slice_interpolators) {
     f.ynext(interp->y_offset) = interp->interpolate(f);
   }
+  f.setParallelRegions();
 }
 
 /*!
diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx
index 0a1391907f..036c475cd7 100644
--- a/src/physics/smoothing.cxx
+++ b/src/physics/smoothing.cxx
@@ -350,7 +350,7 @@ BoutReal Average_XY(const Field2D& var) {
   return Vol_Glb;
 }
 
-BoutReal Vol_Integral(const Field2D& var) {
+BoutReal Vol_Integral([[maybe_unused]] const Field2D& var) {
 #if BOUT_USE_METRIC_3D
   AUTO_TRACE();
   throw BoutException("Vol_Intregral currently incompatible with 3D metrics");
diff --git a/tests/integrated/test-laplace-hypre3d/CMakeLists.txt b/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
index 2645c18c67..b09b416feb 100644
--- a/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
+++ b/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
@@ -7,4 +7,5 @@ bout_add_integrated_test(test-laplace-hypre3d
     data_slab_sol/BOUT.inp
   USE_RUNTEST
   REQUIRES BOUT_HAS_HYPRE
+  REQUIRES BOUT_RUN_ALL_TESTS
   )
diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx
index 2758dbe416..06f865b66c 100644
--- a/tests/unit/fake_mesh_fixture.hxx
+++ b/tests/unit/fake_mesh_fixture.hxx
@@ -58,8 +58,13 @@ public:
     test_coords->d1_dx = test_coords->d1_dy = 0.2;
     test_coords->d1_dz = 0.0;
 #if BOUT_USE_METRIC_3D
+    auto Bup = test_coords->Bxy;
+    Bup.yoffset = 1;
+    auto Bdown = test_coords->Bxy;
+    Bdown.yoffset = -1;
     test_coords->Bxy.splitParallelSlices();
-    test_coords->Bxy.yup() = test_coords->Bxy.ydown() = test_coords->Bxy;
+    test_coords->Bxy.yup() = Bup;
+    test_coords->Bxy.ydown() = Bdown;
 #endif
 
     // No call to Coordinates::geometry() needed here
diff --git a/tests/unit/include/bout/test_petsc_indexer.cxx b/tests/unit/include/bout/test_petsc_indexer.cxx
index 0e9f597eae..70745e7203 100644
--- a/tests/unit/include/bout/test_petsc_indexer.cxx
+++ b/tests/unit/include/bout/test_petsc_indexer.cxx
@@ -474,9 +474,9 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
         const int global = index->getGlobal(i.xm()),
                   otherXind = (this->pe_xind == 0) ? this->pe_xind : this->pe_xind - 1,
                   otherYind = this->pe_yind - 1;
-        const Ind3D otherInd =
-            i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
-                     this->yend - this->ystart + 1, 0);
+        Ind3D otherInd = i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
+                                  this->yend - this->ystart + 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -486,7 +486,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field3D>(indexers, this->pe_xind, this->pe_yind - 1)
-              ->getGlobal(i.yp(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(this->yend - this->ystart + 1));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -494,9 +494,10 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
                   otherXind = (this->pe_xind == this->nxpe - 1) ? this->pe_xind
                                                                 : this->pe_xind + 1,
                   otherYind = this->pe_yind - 1;
-        const Ind3D otherInd =
+        Ind3D otherInd =
             i.offset((this->pe_xind == this->nxpe - 1) ? 1 : this->xstart - this->xend,
                      this->yend - this->ystart + 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -509,9 +510,9 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
         const int global = index->getGlobal(i.xm()),
                   otherXind = (this->pe_xind == 0) ? this->pe_xind : this->pe_xind - 1,
                   otherYind = this->pe_yind + 1;
-        const Ind3D otherInd =
-            i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
-                     this->ystart - this->yend - 1, 0);
+        Ind3D otherInd = i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
+                                  this->ystart - this->yend - 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -521,7 +522,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field3D>(indexers, this->pe_xind, this->pe_yind + 1)
-              ->getGlobal(i.ym(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(-(this->yend - this->ystart + 1)));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -529,9 +530,10 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
                   otherXind = (this->pe_xind == this->nxpe - 1) ? this->pe_xind
                                                                 : this->pe_xind + 1,
                   otherYind = this->pe_yind + 1;
-        const Ind3D otherInd =
+        Ind3D otherInd =
             i.offset((this->pe_xind == this->nxpe - 1) ? 1 : this->xstart - this->xend,
                      this->ystart - this->yend - 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -577,7 +579,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex2D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field2D>(indexers, this->pe_xind, this->pe_yind - 1)
-              ->getGlobal(i.yp(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(this->yend - this->ystart + 1));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -612,7 +614,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex2D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field2D>(indexers, this->pe_xind, this->pe_yind + 1)
-              ->getGlobal(i.ym(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(-(this->yend - this->ystart + 1)));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
diff --git a/tests/unit/include/bout/test_petsc_matrix.cxx b/tests/unit/include/bout/test_petsc_matrix.cxx
index 4f8a4291d6..7200504aa8 100644
--- a/tests/unit/include/bout/test_petsc_matrix.cxx
+++ b/tests/unit/include/bout/test_petsc_matrix.cxx
@@ -55,7 +55,7 @@ class PetscMatrixTest : public FakeMeshFixture {
     if constexpr (std::is_same_v<F, FieldPerp>) {
       indexB = indexA.zp();
     } else {
-      indexB = indexA.yp();
+      indexB = indexA.yp_no_parallel_shift();
     }
     iWU0 = indexB.xm();
     iWU1 = indexB;
@@ -65,9 +65,9 @@ class PetscMatrixTest : public FakeMeshFixture {
       iWD1 = indexB;
       iWD2 = indexB.zp();
     } else {
-      iWD0 = indexB.ym();
+      iWD0 = indexB.yp_no_parallel_shift(-1);
       iWD1 = indexB;
-      iWD2 = indexB.yp();
+      iWD2 = indexB.yp_no_parallel_shift();
     }
     std::unique_ptr<MockTransform> transform =
         bout::utils::make_unique<MockTransform>(*bout::globals::mesh);
diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
index 7e111e9c4e..3d7a52607f 100644
--- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx
+++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
@@ -284,25 +284,6 @@ TEST_F(ShiftedMetricTest, ToFromFieldAlignedFieldPerp) {
 }
 
 TEST_F(ShiftedMetricTest, CalcParallelSlices) {
-  // We don't shift in the guard cells, and the parallel slices are
-  // stored offset in y, therefore we need to make new regions that we
-  // can compare the expected and actual outputs over
-  output_info.disable();
-  mesh->addRegion3D("RGN_YUP",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart + 1, mesh->yend + 1,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  mesh->addRegion3D("RGN_YUP2",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart + 2, mesh->yend + 2,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-
-  mesh->addRegion3D("RGN_YDOWN",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart - 1, mesh->yend - 1,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  mesh->addRegion3D("RGN_YDOWN2",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart - 2, mesh->yend - 2,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  output_info.enable();
-
   // Actual interesting bit here!
   input.getCoordinates()->getParallelTransform().calcParallelSlices(input);
   // Expected output values
@@ -412,9 +393,10 @@ TEST_F(ShiftedMetricTest, CalcParallelSlices) {
                               {0., 0., 0., 0., 0.},
                               {0., 0., 0., 0., 0.}}});
 
-  EXPECT_TRUE(IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YUP", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YUP2", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YDOWN", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YDOWN2", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YPAR_+1", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YPAR_+2", FFTTolerance));
+  EXPECT_TRUE(
+      IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YPAR_-1", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YPAR_-2", FFTTolerance));
 }
 #endif
diff --git a/tests/unit/mesh/test_coordinates_accessor.cxx b/tests/unit/mesh/test_coordinates_accessor.cxx
index 16985981bb..1b96bb731e 100644
--- a/tests/unit/mesh/test_coordinates_accessor.cxx
+++ b/tests/unit/mesh/test_coordinates_accessor.cxx
@@ -81,8 +81,13 @@ TEST_F(CoordinatesAccessorTest, ClearBoth) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());
@@ -121,8 +126,13 @@ TEST_F(CoordinatesAccessorTest, ClearOneTwo) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());
@@ -163,8 +173,13 @@ TEST_F(CoordinatesAccessorTest, ClearTwoOneNone) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());