diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs
index 43a5292e9d..945f769e70 100644
--- a/.git-blame-ignore-revs
+++ b/.git-blame-ignore-revs
@@ -16,6 +16,7 @@ d1cfb8abd6aa5c76e6c1a4d7ab20929c65f8afc2
 8d5cb39e03c2644715a50684f8cd0318b4e82676
 ec69e8838be2dde140a915e50506f8e1ce3cb534
 f2bc0488a298f136294c523bc5ab4086d090059b
+c59626a0498b0cfd1be46f84f0fa38b4855a43cc
 1b4707187a3a85126338303dc766280b8fb2dc56
 b2e2f3575e68f771be2a4341af5fd14e2870469e
 64e4285ec38569f66d31e589a4610cefd16d8076
diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx
new file mode 100644
index 0000000000..765f83e694
--- /dev/null
+++ b/include/bout/boundary_iterator.hxx
@@ -0,0 +1,218 @@
+#pragma once
+
+#include "bout/assert.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/build_defines.hxx"
+#include "bout/field2d.hxx"
+#include "bout/field3d.hxx"
+#include "bout/mesh.hxx"
+#include "bout/parallel_boundary_region.hxx"
+#include "bout/region.hxx"
+#include "bout/sys/parallel_stencils.hxx"
+#include "bout/sys/range.hxx"
+#include <algorithm>
+#include <functional>
+
+class BoundaryRegionIter {
+public:
+  virtual ~BoundaryRegionIter() = default;
+  BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh)
+      : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) {
+    ASSERT3(bx * by == 0);
+  }
+  bool operator!=(const BoundaryRegionIter& rhs) const { return ind() != rhs.ind(); }
+
+  Ind3D ind() const { return xyz2ind(x, y, z); }
+  BoundaryRegionIter& operator++() {
+    ASSERT3(z < nz());
+    z++;
+    if (z == nz()) {
+      z = 0;
+      _next();
+    }
+    return *this;
+  }
+  virtual void _next() = 0;
+  BoundaryRegionIter& operator*() { return *this; }
+
+  void dirichlet_o2(Field3D& f, BoutReal value) const {
+    ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value);
+  }
+
+  BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); }
+
+  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] * 3 - yprev(f)) * 0.5;
+  }
+
+  BoutReal extrapolate_next_o2(const Field3D& f) const {
+    return (2 * f[ind()]) - yprev(f);
+  }
+
+  BoutReal
+  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return (2 * f(0, ind())) - f(0, ind().yp(-by).xp(-bx));
+  }
+
+  BoutReal interpolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] + ynext(f)) * 0.5;
+  }
+
+  BoutReal
+  interpolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return (f(0, ind()) + f(0, ind().yp(-by).xp(-bx))) * 0.5;
+  }
+
+  BoutReal
+  extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)));
+  }
+
+  BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
+    const BoutReal val = ythis(f);
+    const BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac;
+    return 0.5 * (val + next);
+  }
+
+  void set_free(Field3D& f, SheathLimitMode mode) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
+    BoutReal val = ythis(f);
+    if (mode == SheathLimitMode::linear_free) {
+      for (int i = 1; i <= localmesh->ystart; ++i) {
+        val += fac;
+        f[ind().yp(by * i).xp(bx * i)] = val;
+      }
+    } else {
+      for (int i = 1; i <= localmesh->ystart; ++i) {
+        val *= fac;
+        f[ind().yp(by * i).xp(bx * i)] = val;
+      }
+    }
+  }
+
+  void limitFree(Field3D& f) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f));
+    BoutReal val = ythis(f);
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val *= fac;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  bool is_lower() const {
+    ASSERT2(bx == 0);
+    return by == -1;
+  }
+
+  void neumann_o1(Field3D& f, BoutReal grad) const {
+    BoutReal val = ythis(f);
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val += grad;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  void neumann_o2(Field3D& f, BoutReal grad) const {
+    BoutReal val = yprev(f) + grad;
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val += grad;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  void limit_at_least(Field3D& f, BoutReal value) const {
+    ynext(f) = std::max(ynext(f), value);
+  }
+
+  BoutReal& ynext(Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  BoutReal& ythis(Field3D& f) const { return f[ind()]; }
+  const BoutReal& ythis(const Field3D& f) const { return f[ind()]; }
+
+  void setYPrevIfValid(Field3D& f, BoutReal val) const { yprev(f) = val; }
+  void setAll(Field3D& f, const BoutReal val) const {
+    for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) {
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  static int abs_offset() { return 1; }
+
+#if BOUT_USE_METRIC_3D == 0
+  BoutReal& ynext(Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+#endif
+
+  const int dir;
+
+protected:
+  int z{0};
+  int x;
+  int y;
+  const int bx;
+  const int by;
+
+private:
+  Mesh* localmesh;
+  int nx() const { return localmesh->LocalNx; }
+  int ny() const { return localmesh->LocalNy; }
+  int nz() const { return localmesh->LocalNz; }
+
+  Ind3D xyz2ind(int x, int y, int z) const {
+    return Ind3D{((x * ny() + y) * nz()) + z, ny(), nz()};
+  }
+};
+
+class BoundaryRegionIterY : public BoundaryRegionIter {
+public:
+  ~BoundaryRegionIterY() override = default;
+  BoundaryRegionIterY(const RangeIterator& r, int y, int dir, bool is_end, Mesh* mesh)
+      : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {}
+
+  bool operator!=(const BoundaryRegionIterY& rhs) {
+    ASSERT2(y == rhs.y);
+    if (is_end) {
+      if (rhs.is_end) {
+        return false;
+      }
+      return !rhs.r.isDone();
+    }
+    if (rhs.is_end) {
+      return !r.isDone();
+    }
+    return x != rhs.x;
+  }
+
+  void _next() override {
+    ++r;
+    x = r.ind;
+  }
+
+private:
+  RangeIterator r;
+  bool is_end;
+};
+
+class NewBoundaryRegionY {
+public:
+  NewBoundaryRegionY(Mesh* mesh, bool lower, const RangeIterator& r)
+      : mesh(mesh), lower(lower), r(r) {}
+  BoundaryRegionIterY begin(bool begin = true) {
+    return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1,
+                               !begin, mesh);
+  }
+  BoundaryRegionIterY end() { return begin(false); }
+
+private:
+  Mesh* mesh;
+  bool lower;
+  RangeIterator r;
+};
diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 785b139fbb..81cba6d809 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -498,7 +498,8 @@ public:
   /// This uses 2nd order central differences to set the value
   /// on the boundary to the value on the boundary in field \p f3d.
   /// Note: does not just copy values in boundary region.
-  void setBoundaryTo(const Field3D& f3d);
+  void setBoundaryTo(const Field3D& f3d) { setBoundaryTo(f3d, true); }
+  void setBoundaryTo(const Field3D& f3d, bool copyParallelSlices);
 
   using FieldData::applyParallelBoundary;
   void applyParallelBoundary() override;
diff --git a/include/bout/parallel_boundary_op.hxx b/include/bout/parallel_boundary_op.hxx
index d8620e892b..9e551ebc17 100644
--- a/include/bout/parallel_boundary_op.hxx
+++ b/include/bout/parallel_boundary_op.hxx
@@ -49,7 +49,7 @@ protected:
   enum class ValueType { GEN, FIELD, REAL };
   const ValueType value_type{ValueType::REAL};
 
-  BoutReal getValue(const BoundaryRegionPar& bndry, BoutReal t);
+  BoutReal getValue(const BoundaryRegionParIter& bndry, BoutReal t);
 };
 
 template <class T, bool isNeumann = false>
@@ -95,12 +95,13 @@ public:
 
     auto dy = f.getCoordinates()->dy;
 
-    for (bndry->first(); !bndry->isDone(); bndry->next()) {
-      BoutReal value = getValue(*bndry, t);
+    for (auto pnt : *bndry) {
+      //for (bndry->first(); !bndry->isDone(); bndry->next()) {
+      BoutReal value = getValue(pnt, t);
       if (isNeumann) {
-        value *= dy[bndry->ind()];
+        value *= dy[pnt.ind()];
       }
-      static_cast<T*>(this)->apply_stencil(f, bndry, value);
+      static_cast<T*>(this)->apply_stencil(f, pnt, value);
     }
   }
 };
@@ -111,24 +112,27 @@ public:
 class BoundaryOpPar_dirichlet_o1 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o1> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o1(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o1(f, value);
   }
 };
 
 class BoundaryOpPar_dirichlet_o2 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o2> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o2(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o2(f, value);
   }
 };
 
 class BoundaryOpPar_dirichlet_o3 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o3> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o3(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o3(f, value);
   }
 };
 
@@ -136,8 +140,9 @@ class BoundaryOpPar_neumann_o1
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o1, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o1(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o1(f, value);
   }
 };
 
@@ -145,8 +150,9 @@ class BoundaryOpPar_neumann_o2
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o2, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o2(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o2(f, value);
   }
 };
 
@@ -154,8 +160,9 @@ class BoundaryOpPar_neumann_o3
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o3, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o3(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o3(f, value);
   }
 };
 
diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx
index 4d5278d00f..d42a2b2366 100644
--- a/include/bout/parallel_boundary_region.hxx
+++ b/include/bout/parallel_boundary_region.hxx
@@ -1,10 +1,22 @@
 #ifndef BOUT_PAR_BNDRY_H
 #define BOUT_PAR_BNDRY_H
 
+#include "bout/assert.hxx"
 #include "bout/boundary_region.hxx"
+#include "bout/bout_enum_class.hxx"
 #include "bout/bout_types.hxx"
+#include <algorithm>
+#include <cmath>
+#include <cstdlib>
+#include <functional>
+#include <string>
 #include <vector>
 
+#include "bout/build_defines.hxx"
+#include "bout/field2d.hxx"
+#include "bout/region.hxx"
+#include "bout/sys/parallel_stencils.hxx"
+#include "bout/utils.hxx"
 #include <bout/field3d.hxx>
 #include <bout/mesh.hxx>
 
@@ -14,99 +26,100 @@
  *
  */
 
-namespace parallel_stencil {
-// generated by src/mesh/parallel_boundary_stencil.cxx.py
-inline BoutReal pow(BoutReal val, int exp) {
-  // constexpr int expval = exp;
-  // static_assert(expval == 2 or expval == 3, "This pow is only for exponent 2 or 3");
-  if (exp == 2) {
-    return val * val;
-  }
-  ASSERT3(exp == 3);
-  return val * val * val;
-}
-inline BoutReal dirichlet_o1(BoutReal UNUSED(spacing0), BoutReal value0) {
-  return value0;
-}
-inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
-                             BoutReal value1) {
-  return (spacing0 * value1 - spacing1 * value0) / (spacing0 - spacing1);
-}
-inline BoutReal neumann_o2(BoutReal UNUSED(spacing0), BoutReal value0, BoutReal spacing1,
-                           BoutReal value1) {
-  return -spacing1 * value0 + value1;
-}
-inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
-                             BoutReal value1, BoutReal spacing2, BoutReal value2) {
-  return (pow(spacing0, 2) * spacing1 * value2 - pow(spacing0, 2) * spacing2 * value1
-          - spacing0 * pow(spacing1, 2) * value2 + spacing0 * pow(spacing2, 2) * value1
-          + pow(spacing1, 2) * spacing2 * value0 - spacing1 * pow(spacing2, 2) * value0)
-         / ((spacing0 - spacing1) * (spacing0 - spacing2) * (spacing1 - spacing2));
-}
-inline BoutReal neumann_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
-                           BoutReal value1, BoutReal spacing2, BoutReal value2) {
-  return (2 * spacing0 * spacing1 * value2 - 2 * spacing0 * spacing2 * value1
-          + pow(spacing1, 2) * spacing2 * value0 - pow(spacing1, 2) * value2
-          - spacing1 * pow(spacing2, 2) * value0 + pow(spacing2, 2) * value1)
-         / ((spacing1 - spacing2) * (2 * spacing0 - spacing1 - spacing2));
-}
-} // namespace parallel_stencil
+BOUT_ENUM_CLASS(SheathLimitMode, limit_free, exponential_free, linear_free);
 
-class BoundaryRegionPar : public BoundaryRegionBase {
+namespace bout {
+namespace parallel_boundary_region {
 
-  struct RealPoint {
-    BoutReal s_x;
-    BoutReal s_y;
-    BoutReal s_z;
-  };
-
-  struct Indices {
-    // Indices of the boundary point
-    Ind3D index;
-    // Intersection with boundary in index space
-    RealPoint intersection;
-    // Distance to intersection
-    BoutReal length;
-    // Angle between field line and boundary
-    // BoutReal angle;
-    // How many points we can go in the opposite direction
-    signed char valid;
-  };
-
-  using IndicesVec = std::vector<Indices>;
-  using IndicesIter = IndicesVec::iterator;
+struct RealPoint {
+  BoutReal s_x;
+  BoutReal s_y;
+  BoutReal s_z;
+};
 
-  /// Vector of points in the boundary
-  IndicesVec bndry_points;
-  /// Current position in the boundary points
-  IndicesIter bndry_position;
+struct Indices {
+  // Indices of the boundary point
+  Ind3D index;
+  // Intersection with boundary in index space
+  RealPoint intersection;
+  // Distance to intersection
+  BoutReal length;
+  // Angle between field line and boundary
+  // BoutReal angle;
+  // How many points we can go in the opposite direction
+  signed char valid;
+  signed char offset;
+  unsigned char abs_offset;
+  Indices(Ind3D index, RealPoint&& intersection, BoutReal length, signed char valid,
+          signed char offset, unsigned char abs_offset)
+      : index(index), intersection(intersection), length(length), valid(valid),
+        offset(offset), abs_offset(abs_offset) {};
+};
 
-public:
-  BoundaryRegionPar(const std::string& name, int dir, Mesh* passmesh)
-      : BoundaryRegionBase(name, passmesh), dir(dir) {
-    ASSERT0(std::abs(dir) == 1);
-    BoundaryRegionBase::isParallel = true;
+using IndicesVec = std::vector<Indices>;
+using IndicesIter = IndicesVec::iterator;
+using IndicesIterConst = IndicesVec::const_iterator;
+
+/// Limited free gradient of log of a quantity
+/// This ensures that the guard cell values remain positive
+/// while also ensuring that the quantity never increases
+///
+///  fm  fc | fp
+///         ^ boundary
+///
+/// exp( 2*log(fc) - log(fm) )
+inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc, SheathLimitMode mode) {
+  if ((fm < fc) && (mode == SheathLimitMode::limit_free)) {
+    return fc; // Neumann rather than increasing into boundary
   }
-  BoundaryRegionPar(const std::string& name, BndryLoc loc, int dir, Mesh* passmesh)
-      : BoundaryRegionBase(name, loc, passmesh), dir(dir) {
-    BoundaryRegionBase::isParallel = true;
-    ASSERT0(std::abs(dir) == 1);
+  if (fm < 1e-10) {
+    return fc; // Low / no density condition
   }
 
-  /// Add a point to the boundary
-  void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length,
-                 signed char valid) {
-    bndry_points.push_back({ind, {x, y, z}, length, valid});
+  BoutReal fp = 0;
+  switch (mode) {
+  case SheathLimitMode::limit_free:
+  case SheathLimitMode::exponential_free:
+    fp = SQ(fc) / fm; // Exponential
+    break;
+  case SheathLimitMode::linear_free:
+    fp = (2.0 * fc) - fm; // Linear
+    break;
   }
-  void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z,
-                 BoutReal length, signed char valid) {
-    bndry_points.push_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid});
+
+#if CHECKLEVEL >= 2
+  if (!std::isfinite(fp)) {
+    throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp);
   }
+#endif
 
-  // final, so they can be inlined
-  void first() final { bndry_position = begin(bndry_points); }
-  void next() final { ++bndry_position; }
-  bool isDone() final { return (bndry_position == end(bndry_points)); }
+  return fp;
+}
+
+inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) {
+  if (fm < fc) {
+    return 1; // Neumann rather than increasing into boundary
+  }
+  if (fm < 1e-10) {
+    return 1; // Low / no density condition
+  }
+  BoutReal fp = fc / fm;
+#if CHECKLEVEL >= 2
+  if (!std::isfinite(fp)) {
+    throw BoutException("SheathBoundaryParallel limitFree: {}, {} -> {}", fm, fc, fp);
+  }
+#endif
+  return fp;
+}
+
+template <class IndicesVec, class IndicesIter>
+class BoundaryRegionParIterBase {
+
+public:
+  BoundaryRegionParIterBase(IndicesVec& bndry_points, IndicesIter bndry_position, int dir,
+                            Mesh* localmesh)
+      : bndry_points(bndry_points), bndry_position(bndry_position), dir(dir),
+        localmesh(localmesh) {};
 
   // getter
   Ind3D ind() const { return bndry_position->index; }
@@ -115,37 +128,99 @@ public:
   BoutReal s_z() const { return bndry_position->intersection.s_z; }
   BoutReal length() const { return bndry_position->length; }
   signed char valid() const { return bndry_position->valid; }
+  signed char offset() const { return bndry_position->offset; }
+  unsigned char abs_offset() const { return bndry_position->abs_offset; }
 
   // setter
-  void setValid(signed char val) { bndry_position->valid = val; }
+  void setValid(signed char valid) { bndry_position->valid = valid; }
 
-  bool contains(const BoundaryRegionPar& bndry) const {
-    return std::binary_search(
-        begin(bndry_points), end(bndry_points), *bndry.bndry_position,
-        [](const Indices& i1, const Indices& i2) { return i1.index < i2.index; });
+  // extrapolate a given point to the boundary
+  BoutReal extrapolate_sheath_o1(const Field3D& f) const { return ythis(f); }
+  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
+    ASSERT3(valid() >= 0);
+    if (valid() < 1) {
+      return extrapolate_sheath_o1(f);
+    }
+    return ythis(f) * (1 + length()) - yprev(f) * length();
+  }
+  BoutReal
+  extrapolate_sheath_o1(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return ythis(f);
+  }
+  BoutReal
+  extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    ASSERT3(valid() >= 0);
+    if (valid() < 1) {
+      return extrapolate_sheath_o1(f);
+    }
+    return ythis(f) * (1 + length()) - yprev(f) * length();
   }
 
-  // extrapolate a given point to the boundary
-  BoutReal extrapolate_o1(const Field3D& f) const { return f[ind()]; }
-  BoutReal extrapolate_o2(const Field3D& f) const {
+  BoutReal interpolate_sheath_o2(const Field3D& f) const {
+    return ythis(f) * (1 - length()) + ynext(f) * length();
+  }
+  BoutReal
+  interpolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return ythis(f) * (1 - length()) + ynext(f) * length();
+  }
+
+  BoutReal extrapolate_next_o1(const Field3D& f) const { return ythis(f); }
+  BoutReal extrapolate_next_o2(const Field3D& f) const {
+    ASSERT3(valid() >= 0);
+    if (valid() < 1) {
+      return extrapolate_next_o1(f);
+    }
+    return ythis(f) * 2 - yprev(f);
+  }
+
+  BoutReal
+  extrapolate_next_o1(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return ythis(f);
+  }
+  BoutReal
+  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    ASSERT3(valid() >= 0);
+    if (valid() < 1) {
+      return extrapolate_next_o1(f);
+    }
+    return ythis(f) * 2 - yprev(f);
+  }
+
+  // extrapolate the gradient into the boundary
+  BoutReal extrapolate_grad_o1([[maybe_unused]] const Field3D& f) const { return 0; }
+  BoutReal extrapolate_grad_o2(const Field3D& f) const {
     ASSERT3(valid() >= 0);
     if (valid() < 1) {
-      return extrapolate_o1(f);
+      return extrapolate_grad_o1(f);
     }
-    return f[ind()] * (1 + length()) - f.ynext(-dir)[ind().yp(-dir)] * length();
+    return ythis(f) - ynext(f);
   }
 
+  BoundaryRegionParIterBase& operator*() { return *this; }
+
+  BoundaryRegionParIterBase& operator++() {
+    ++bndry_position;
+    return *this;
+  }
+
+  bool operator!=(const BoundaryRegionParIterBase& rhs) {
+    return bndry_position != rhs.bndry_position;
+  }
+
+#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i)
   // dirichlet boundary code
   void dirichlet_o1(Field3D& f, BoutReal value) const {
-    f.ynext(dir)[ind().yp(dir)] = value;
+    ITER() { getAt(f, i) = value; }
   }
 
   void dirichlet_o2(Field3D& f, BoutReal value) const {
     if (length() < small_value) {
       return dirichlet_o1(f, value);
     }
-    ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 1 - length(), value);
-    // ynext(f) = f[ind()] * (1 + 1/length()) + value / length();
+    ITER() {
+      getAt(f, i) =
+          parallel_stencil::dirichlet_o2(i + 1, ythis(f), i + 1 - length(), value);
+    }
   }
 
   void dirichlet_o3(Field3D& f, BoutReal value) const {
@@ -154,17 +229,34 @@ public:
       return dirichlet_o2(f, value);
     }
     if (length() < small_value) {
-      ynext(f) = parallel_stencil::dirichlet_o2(2, yprev(f), 1 - length(), value);
+      ITER() {
+        getAt(f, i) =
+            parallel_stencil::dirichlet_o2(i + 2, yprev(f), i + 1 - length(), value);
+      }
     } else {
-      ynext(f) =
-          parallel_stencil::dirichlet_o3(2, yprev(f), 1, f[ind()], 1 - length(), value);
+      ITER() {
+        getAt(f, i) = parallel_stencil::dirichlet_o3(i + 2, yprev(f), i + 1, ythis(f),
+                                                     i + 1 - length(), value);
+      }
+    }
+  }
+
+  void limit_at_least(Field3D& f, BoutReal value) const {
+    ITER() {
+      if (getAt(f, i) < value) {
+        getAt(f, i) = value;
+      }
     }
   }
 
+  bool is_lower() const { return dir == -1; }
+
   // NB: value needs to be scaled by dy
   // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil.
   // But in general we do not, and thus for normal C2 stencils, this is 1st order.
-  void neumann_o1(Field3D& f, BoutReal value) const { ynext(f) = f[ind()] + value; }
+  void neumann_o1(Field3D& f, BoutReal value) const {
+    ITER() { getAt(f, i) = ythis(f) + value * (i + 1); }
+  }
 
   // NB: value needs to be scaled by dy
   void neumann_o2(Field3D& f, BoutReal value) const {
@@ -172,34 +264,227 @@ public:
     if (valid() < 1) {
       return neumann_o1(f, value);
     }
-    ynext(f) = yprev(f) + 2 * value;
+    ITER() { getAt(f, i) = yprev(f) + (2 + i) * value; }
   }
 
   // NB: value needs to be scaled by dy
   void neumann_o3(Field3D& f, BoutReal value) const {
     ASSERT3(valid() >= 0);
     if (valid() < 1) {
-      return neumann_o1(f, value);
+      return neumann_o2(f, value);
+    }
+    ITER() {
+      getAt(f, i) = parallel_stencil::neumann_o3(i + 1 - length(), value, i + 1, ythis(f),
+                                                 2, yprev(f));
     }
-    ynext(f) =
-        parallel_stencil::neumann_o3(1 - length(), value, 1, f[ind()], 2, yprev(f));
   }
 
-  const int dir;
+  // extrapolate into the boundary using only monotonic decreasing values.
+  // f needs to be positive
+  void limitFree(Field3D& f) const {
+    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f)) : 1;
+    auto val = ythis(f);
+    ITER() {
+      val *= fac;
+      getAt(f, i) = val;
+    }
+  }
+
+  BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const {
+    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode)
+                                 : (mode == SheathLimitMode::linear_free ? 0 : 1);
+    auto val = ythis(f);
+    BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac;
+    return val * length() + next * (1 - length());
+  }
+
+  void set_free(Field3D& f, SheathLimitMode mode) const {
+    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode)
+                                 : (mode == SheathLimitMode::linear_free ? 0 : 1);
+    auto val = ythis(f);
+    if (mode == SheathLimitMode::linear_free) {
+      ITER() {
+        val += fac;
+        getAt(f, i) = val;
+      }
+    } else {
+      ITER() {
+        val *= fac;
+        getAt(f, i) = val;
+      }
+    }
+  }
+
+  void setAll(Field3D& f, const BoutReal val) const {
+    for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) {
+      f.ynext(i)[ind().yp(i)] = val;
+    }
+  }
+
+  template <bool check = true>
+  BoutReal& getAt(Field3D& f, int off) const {
+    ASSERT3(f.hasParallelSlices());
+    if constexpr (check) {
+      ASSERT3(valid() > -off - 2);
+    }
+    auto _off = offset() + off * dir;
+    return f.ynext(_off)[ind().yp(_off)];
+  }
+  template <bool check = true>
+  const BoutReal& getAt(const Field3D& f, int off) const {
+    ASSERT3(f.hasParallelSlices());
+    if constexpr (check) {
+      ASSERT3(valid() > -off - 2);
+    }
+    auto _off = offset() + off * dir;
+    return f.ynext(_off)[ind().yp(_off)];
+  }
+
+  const BoutReal& ynext(const Field3D& f) const { return getAt(f, 0); }
+  BoutReal& ynext(Field3D& f) const { return getAt(f, 0); }
+  const BoutReal& ythis(const Field3D& f) const { return getAt(f, -1); }
+  BoutReal& ythis(Field3D& f) const { return getAt(f, -1); }
+  const BoutReal& yprev(const Field3D& f) const { return getAt(f, -2); }
+  BoutReal& yprev(Field3D& f) const { return getAt(f, -2); }
+
+  template <bool check = true>
+  BoutReal getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
+                 int off) const {
+    if constexpr (check) {
+      ASSERT3(valid() > -off - 2);
+    }
+    auto _off = offset() + off * dir;
+    return f(_off, ind().yp(_off));
+  }
+  BoutReal ynext(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return getAt(f, 0);
+  }
+  BoutReal ythis(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return getAt(f, -1);
+  }
+  BoutReal yprev(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return getAt(f, -2);
+  }
+
+  void setYPrevIfValid(Field3D& f, BoutReal val) const {
+    if (valid() > 0) {
+      yprev(f) = val;
+    }
+  }
+
+#if BOUT_USE_METRIC_3D == 0
+  const BoutReal& ynext(const Field2D& f) const { return f.ynext(dir)[ind().yp(dir)]; }
+  BoutReal& ynext(Field2D& f) const { return f.ynext(dir)[ind().yp(dir)]; }
+
+  const BoutReal& yprev(const Field2D& f) const {
+    ASSERT3(valid() > 0);
+    return f.ynext(-dir)[ind().yp(-dir)];
+  }
+  BoutReal& yprev(Field2D& f) const {
+    ASSERT3(valid() > 0);
+    return f.ynext(-dir)[ind().yp(-dir)];
+  }
+#endif
 
 private:
+  const IndicesVec& bndry_points;
+  IndicesIter bndry_position;
+
   constexpr static BoutReal small_value = 1e-2;
 
-  // BoutReal get(const Field3D& f, int off)
-  const BoutReal& ynext(const Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; }
-  BoutReal& ynext(Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; }
-  const BoutReal& yprev(const Field3D& f) const { return f.ynext(-dir)[ind().yp(-dir)]; }
-  BoutReal& yprev(Field3D& f) const { return f.ynext(-dir)[ind().yp(-dir)]; }
+public:
+  const int dir;
+  Mesh* localmesh;
+};
+} // namespace parallel_boundary_region
+} // namespace bout
+using BoundaryRegionParIter = bout::parallel_boundary_region::BoundaryRegionParIterBase<
+    bout::parallel_boundary_region::IndicesVec,
+    bout::parallel_boundary_region::IndicesIter>;
+using BoundaryRegionParIterConst =
+    bout::parallel_boundary_region::BoundaryRegionParIterBase<
+        const bout::parallel_boundary_region::IndicesVec,
+        bout::parallel_boundary_region::IndicesIterConst>;
+
+class BoundaryRegionPar : public BoundaryRegionBase {
+public:
+  BoundaryRegionPar(const std::string& name, int dir, Mesh* passmesh)
+      : BoundaryRegionBase(name, passmesh), dir(dir) {
+    ASSERT0(std::abs(dir) == 1);
+    BoundaryRegionBase::isParallel = true;
+  }
+  BoundaryRegionPar(const std::string& name, BndryLoc loc, int dir, Mesh* passmesh)
+      : BoundaryRegionBase(name, loc, passmesh), dir(dir) {
+    BoundaryRegionBase::isParallel = true;
+    ASSERT0(std::abs(dir) == 1);
+  }
+
+  /// Add a point to the boundary
+  void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length,
+                 char valid, signed char offset) {
+    if (!bndry_points.empty() && bndry_points.back().index > ind) {
+      is_sorted = false;
+    }
+    bndry_points.emplace_back(ind, bout::parallel_boundary_region::RealPoint{x, y, z},
+                              length, valid, offset,
+                              static_cast<unsigned char>(std::abs(offset)));
+  }
+  void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z,
+                 BoutReal length, char valid, signed char offset) {
+    add_point(xyz2ind(ix, iy, iz, localmesh), x, y, z, length, valid, offset);
+  }
+
+  // final, so they can be inlined
+  void first() final { bndry_position = std::begin(bndry_points); }
+  void next() final { ++bndry_position; }
+  bool isDone() final { return (bndry_position == std::end(bndry_points)); }
+
+  bool contains(const BoundaryRegionPar& bndry) const {
+    ASSERT2(is_sorted);
+    return std::binary_search(std::begin(bndry_points), std::end(bndry_points),
+                              *bndry.bndry_position,
+                              [](const bout::parallel_boundary_region::Indices& i1,
+                                 const bout::parallel_boundary_region::Indices& i2) {
+                                return i1.index < i2.index;
+                              });
+  }
+
+  bool contains(const int ix, const int iy, const int iz) const {
+    const auto i2 = xyz2ind(ix, iy, iz, localmesh);
+    return std::ranges::any_of(bndry_points.begin(), bndry_points.end(),
+                               [&i2](auto i1) { return i1.index == i2; });
+  }
+
+  // setter
+  void setValid(char val) { bndry_position->valid = val; }
+
+  // BoundaryRegionParIterConst begin() const {
+  //   return BoundaryRegionParIterConst(bndry_points, bndry_points.begin(), dir);
+  // }
+  // BoundaryRegionParIterConst end() const {
+  //   return BoundaryRegionParIterConst(bndry_points, bndry_points.begin(), dir);
+  // }
+  BoundaryRegionParIter begin() {
+    return BoundaryRegionParIter(bndry_points, bndry_points.begin(), dir, localmesh);
+  }
+  BoundaryRegionParIter end() {
+    return BoundaryRegionParIter(bndry_points, bndry_points.end(), dir, localmesh);
+  }
+
+  const int dir;
+
+private:
+  /// Vector of points in the boundary
+  bout::parallel_boundary_region::IndicesVec bndry_points;
+  /// Current position in the boundary points
+  bout::parallel_boundary_region::IndicesIter bndry_position;
+
   static Ind3D xyz2ind(int x, int y, int z, Mesh* mesh) {
     const int ny = mesh->LocalNy;
     const int nz = mesh->LocalNz;
     return Ind3D{(x * ny + y) * nz + z, ny, nz};
   }
+  bool is_sorted{true};
 };
 
 #endif //  BOUT_PAR_BNDRY_H
diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx
new file mode 100644
index 0000000000..2cbec9c8b0
--- /dev/null
+++ b/include/bout/sys/parallel_stencils.hxx
@@ -0,0 +1,41 @@
+#pragma once
+#include <bout/assert.hxx>
+#include <bout/bout_types.hxx>
+
+namespace parallel_stencil {
+// generated by src/mesh/parallel_boundary_stencil.cxx.py
+inline BoutReal pow(BoutReal val, int exp) {
+  // constexpr int expval = exp;
+  // static_assert(expval == 2 or expval == 3, "This pow is only for exponent 2 or 3");
+  if (exp == 2) {
+    return val * val;
+  }
+  ASSERT3(exp == 3);
+  return val * val * val;
+}
+inline BoutReal dirichlet_o1([[maybe_unused]] BoutReal spacing0, BoutReal value0) {
+  return value0;
+}
+inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
+                             BoutReal value1) {
+  return (spacing0 * value1 - spacing1 * value0) / (spacing0 - spacing1);
+}
+inline BoutReal neumann_o2([[maybe_unused]] BoutReal spacing0, BoutReal value0,
+                           BoutReal spacing1, BoutReal value1) {
+  return (-spacing1 * value0) + value1;
+}
+inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
+                             BoutReal value1, BoutReal spacing2, BoutReal value2) {
+  return (pow(spacing0, 2) * spacing1 * value2 - pow(spacing0, 2) * spacing2 * value1
+          - spacing0 * pow(spacing1, 2) * value2 + spacing0 * pow(spacing2, 2) * value1
+          + pow(spacing1, 2) * spacing2 * value0 - spacing1 * pow(spacing2, 2) * value0)
+         / ((spacing0 - spacing1) * (spacing0 - spacing2) * (spacing1 - spacing2));
+}
+inline BoutReal neumann_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1,
+                           BoutReal value1, BoutReal spacing2, BoutReal value2) {
+  return (2 * spacing0 * spacing1 * value2 - 2 * spacing0 * spacing2 * value1
+          + pow(spacing1, 2) * spacing2 * value0 - pow(spacing1, 2) * value2
+          - spacing1 * pow(spacing2, 2) * value0 + pow(spacing2, 2) * value1)
+         / ((spacing1 - spacing2) * (2 * spacing0 - spacing1 - spacing2));
+}
+} // namespace parallel_stencil
diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx
new file mode 100644
index 0000000000..6254ef7862
--- /dev/null
+++ b/include/bout/yboundary_regions.hxx
@@ -0,0 +1,132 @@
+#pragma once
+
+#include "./boundary_iterator.hxx"
+#include "bout/assert.hxx"
+#include "bout/boutexception.hxx"
+#include "bout/globals.hxx"
+#include "bout/mask.hxx"
+#include "bout/options.hxx"
+#include "bout/parallel_boundary_region.hxx"
+#include "bout/region.hxx"
+#include <memory>
+#include <vector>
+
+/// This class allows to simplify iterating over y-boundaries.
+///
+/// It makes it easier to write code for FieldAligned boundaries, but if a bit
+/// care is taken the code also works with FluxCoordinateIndependent code.
+///
+/// An example how to replace old code is given here:
+/// ../../manual/sphinx/user_docs/boundary_options.rst
+
+class YBoundary {
+public:
+  template <class F>
+  void iter_regions(const F& f) {
+    ASSERT1(is_init);
+    for (auto& region : boundary_regions) {
+      f(*region);
+    }
+    for (auto& region : boundary_regions_par) {
+      f(*region);
+    }
+  }
+  template <class F>
+  void iter_pnts(const F& f) {
+    iter_regions([&](auto& region) {
+      for (auto& pnt : region) {
+        f(pnt);
+      }
+    });
+  }
+
+  template <class F>
+  void iter(const F& f) {
+    return iter_regions(f);
+  }
+
+  void init(Options& options, Mesh* mesh = nullptr) {
+    if (is_init) {
+      if (optptr == &options) {
+        return;
+      }
+      throw BoutException("YBoundary is already initialised!");
+    }
+    if (mesh == nullptr) {
+      mesh = bout::globals::mesh;
+    }
+
+    const bool lower_y =
+        options["lower_y"].doc("Boundary on lower y?").withDefault<bool>(true);
+    const bool upper_y =
+        options["upper_y"].doc("Boundary on upper y?").withDefault<bool>(true);
+    const bool outer_x =
+        options["outer_x"].doc("Boundary on outer x?").withDefault<bool>(true);
+    const bool inner_x =
+        options["inner_x"].doc("Boundary on inner x?").withDefault<bool>(false);
+
+    if (mesh->isFci()) {
+      if (outer_x) {
+        for (auto& bndry : mesh->getBoundariesPar(BoundaryParType::xout)) {
+          boundary_regions_par.push_back(bndry);
+        }
+      }
+      if (inner_x) {
+        for (auto& bndry : mesh->getBoundariesPar(BoundaryParType::xin)) {
+          boundary_regions_par.push_back(bndry);
+        }
+      }
+    } else {
+      if (lower_y) {
+        boundary_regions.push_back(
+            std::make_shared<NewBoundaryRegionY>(mesh, true, mesh->iterateBndryLowerY()));
+      }
+      if (upper_y) {
+        boundary_regions.push_back(std::make_shared<NewBoundaryRegionY>(
+            mesh, false, mesh->iterateBndryUpperY()));
+      }
+    }
+    is_init = true;
+    optptr = &options;
+    // Cache boundary regions
+    _contains.emplace_back(mesh, false);
+    _contains.emplace_back(mesh, false);
+    iter_pnts([&](const auto& pnt) {
+      if (pnt.dir == 1) {
+        _contains[1][pnt.ind()] = true;
+      } else if (pnt.dir == -1) {
+        _contains[0][pnt.ind()] = true;
+      }
+    });
+  }
+
+  bool contains_ylow(Ind3D ind) const { return _contains[0][ind]; }
+  bool contains_yhigh(Ind3D ind) const { return _contains[1][ind]; }
+  template <int dir>
+  bool contains(Ind3D ind) const {
+    static_assert(dir == 1 || dir == -1);
+    if constexpr (dir == 1) {
+      return _contains[1][ind];
+    }
+    if constexpr (dir == -1) {
+      return _contains[0][ind];
+    }
+  }
+  bool contains(int dir, Ind3D ind) const {
+    if (dir == 1) {
+      return contains<+1>(ind);
+    }
+    if (dir == -1) {
+      return contains<-1>(ind);
+    }
+    throw BoutException("only dir == 1 and dir == -1 are implemented, not {}", dir);
+  }
+
+private:
+  std::vector<std::shared_ptr<BoundaryRegionPar>> boundary_regions_par;
+  std::vector<std::shared_ptr<NewBoundaryRegionY>> boundary_regions;
+
+  std::vector<BoutMask> _contains;
+  bool is_init{false};
+  Options* optptr;
+};
diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst
index 8493049516..f04d09cabc 100644
--- a/manual/sphinx/user_docs/boundary_options.rst
+++ b/manual/sphinx/user_docs/boundary_options.rst
@@ -436,6 +436,130 @@ the upper Y boundary of a 2D variable ``var``::
 The `BoundaryRegion` class is defined in
 ``include/boundary_region.hxx``
 
+Y-Boundaries
+------------
+
+The sheath boundaries are often implemented in the physics model.
+Previously of they where implemented using a `RangeIterator`::
+
+    class yboundary_example_legacy {
+    public:
+      yboundary_example_legacy(Options* opt, const Field3D& N, const Field3D& V)
+          : N(N), V(V) {
+        Options& options = *opt;
+        lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault<bool>(lower_y);
+        upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault<bool>(upper_y);
+      }
+
+      void rhs() {
+        BoutReal totalFlux = 0;
+        if (lower_y) {
+          for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
+            for (int jz = 0; jz < mesh->LocalNz; jz++) {
+              // Calculate flux through surface [normalised m^-2 s^-1],
+              // should be positive since V < 0.0
+              BoutReal flux =
+                  -0.5 * (N(r.ind, mesh->ystart, jz) + N(r.ind, mesh->ystart - 1, jz)) * 0.5
+                  * (V(r.ind, mesh->ystart, jz) + V(r.ind, mesh->ystart - 1, jz));
+              totalFlux += flux;
+            }
+          }
+        }
+        if (upper_y) {
+          for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
+            for (int jz = 0; jz < mesh->LocalNz; jz++) {
+              // Calculate flux through surface [normalised m^-2 s^-1],
+              // should be positive since V < 0.0
+              BoutReal flux = -0.5 * (N(r.ind, mesh->yend, jz) + N(r.ind, mesh->yend + 1, jz))
+                              * 0.5
+                              * (V(r.ind, mesh->yend, jz) + V(r.ind, mesh->yend + 1, jz));
+              totalFlux += flux;
+            }
+          }
+        }
+      }
+
+    private:
+      bool lower_y{true};
+      bool upper_y{true};
+      const Field3D& N;
+      const Field3D& V;
+    }
+
+
+This can be replaced using the `YBoundary` class, which not only simplifies the
+code, but also allows to have the same code working with non-field-aligned
+geometries, as flux coordinate independent (FCI) method::
+
+    #include <bout/yboundary_regions.hxx>
+
+    class yboundary_example {
+    public:
+      yboundary_example(Options* opt, const Field3D& N, const Field3D& V) : N(N), V(V) {
+        // Set what kind of yboundaries you want to include
+        yboundary.init(opt);
+      }
+
+      void rhs() {
+        BoutReal totalFlux = 0;
+        yboundary.iter_pnts([&](auto& pnt) {
+          BoutReal flux = pnt.interpolate_sheath_o2(N) * pnt.interpolate_sheath_o2(V);
+        });
+      }
+
+    private:
+      YBoundary ybounday;
+      const Field3D& N;
+      const Field3D& V;
+    };
+
+
+
+There are several member functions of ``pnt``. ``pnt`` is of type
+`BoundaryRegionParIterBase` and `BoundaryRegionIter`, and both should provide
+the same interface. If they don't that is a bug, as the above code is a
+template, that gets instantiated for both types, and thus requires both
+classes to provide the same interface, one for FCI-like boundaries and one for
+field aligned boundaries.
+
+Here is a short summary of some members of ``pnt``, where ``f`` is a :
+
+.. list-table:: Members for boundary operation
+   :widths: 15 70
+   :header-rows: 1
+
+   * - Function
+     - Description
+   * - ``pnt.ythis(f)``
+     - Returns the value at the last point in the domain
+   * - ``pnt.ynext(f)``
+     - Returns the value at the first point in the domain
+   * - ``pnt.yprev(f)``
+     - Returns the value at the second to last point in the domain, if it is
+       valid. NB: this point may not be valid.
+   * - ``pnt.interpolate_sheath_o2(f)``
+     - Returns the value at the boundary, assuming the bounday value has been set
+   * - ``pnt.extrapolate_sheath_o1(f)``
+     - Returns the value at the boundary, extrapolating from the bulk, first order
+   * - ``pnt.extrapolate_sheath_o2(f)``
+     - Returns the value at the boundary, extrapolating from the bulk, second order
+   * - ``pnt.extrapolate_next_o{1,2}(f)``
+     - Extrapolate into the boundary from the bulk, first or second order
+   * - ``pnt.extrapolate_grad_o{1,2}(f)``
+     - Extrapolate the gradient into the boundary, first or second order
+   * - ``pnt.dirichlet_o{1,2,3}(f, v)``
+     - Apply dirichlet boundary conditions with value ``v`` and given order
+   * - ``pnt.neumann_o{1,2,3}(f, v)``
+     - Applies a gradient of ``v / dy`` boundary condition.
+   * - ``pnt.limitFree(f)``
+     - Extrapolate into the boundary using only monotonic decreasing values.
+       ``f`` needs to be positive.
+   * - ``pnt.dir``
+     - The direction of the boundary.
+
+
+
+
 Boundary regions
 ----------------
 
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 9154f78a4f..373d00d667 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -445,14 +445,42 @@ void Field3D::applyTDerivBoundary() {
   }
 }
 
-void Field3D::setBoundaryTo(const Field3D& f3d) {
+void Field3D::setBoundaryTo(const Field3D& f3d, bool copyParallelSlices) {
 
   checkData(f3d);
 
   allocate(); // Make sure data allocated
 
-  /// Loop over boundary regions
+  if (isFci()) {
+    ASSERT1(f3d.hasParallelSlices());
+    if (copyParallelSlices) {
+      splitParallelSlices();
+      for (int i = 0; i < fieldmesh->ystart; ++i) {
+        yup(i) = f3d.yup(i);
+        ydown(i) = f3d.ydown(i);
+      }
+    } else {
+      // Set yup/ydown using midpoint values from f3d
+      ASSERT1(hasParallelSlices());
+
+      for (auto& region : fieldmesh->getBoundariesPar()) {
+        for (const auto& pnt : *region) {
+          // Interpolate midpoint value in f3d
+          const BoutReal val = pnt.interpolate_sheath_o2(f3d);
+          // Set the same boundary value in this field
+          pnt.dirichlet_o1(*this, val);
+        }
+      }
+    }
+  }
+
+  // Non-FCI.
+  // Transform to field-aligned coordinates?
+  // Loop over boundary regions
   for (const auto& reg : fieldmesh->getBoundaries()) {
+    if (isFci() && reg->by != 0) {
+      continue;
+    }
     /// Loop within each region
     for (reg->first(); !reg->isDone(); reg->next()) {
       for (int z = 0; z < nz; z++) {
diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx
index e9182042de..9af6fd20f9 100644
--- a/src/invert/laplace/invert_laplace.cxx
+++ b/src/invert/laplace/invert_laplace.cxx
@@ -218,10 +218,10 @@ Field3D Laplacian::solve(const Field3D& b, const Field3D& x0) {
 
   // Setting the start and end range of the y-slices
   int ys = localmesh->ystart, ye = localmesh->yend;
-  if (localmesh->hasBndryLowerY() && include_yguards) {
+  if (include_yguards && localmesh->hasBndryLowerY()) {
     ys = 0; // Mesh contains a lower boundary
   }
-  if (localmesh->hasBndryUpperY() && include_yguards) {
+  if (include_yguards && localmesh->hasBndryUpperY()) {
     ye = localmesh->LocalNy - 1; // Contains upper boundary
   }
 
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index ea1b6a41b2..a19e416beb 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -643,23 +643,28 @@ int BoutMesh::load() {
   // Add boundary regions
   addBoundaryRegions();
 
-  // Set cached values
-  {
-    int mybndry = static_cast<int>(!(iterateBndryLowerY().isDone()));
-    int allbndry = 0;
-    mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend));
-    has_boundary_lower_y = static_cast<bool>(allbndry);
-  }
-  {
-    int mybndry = static_cast<int>(!(iterateBndryUpperY().isDone()));
-    int allbndry = 0;
-    mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart));
-    has_boundary_upper_y = static_cast<bool>(allbndry);
-  }
-
   // Initialize default coordinates
   getCoordinates();
 
+  // Set cached values
+  if (isFci()) {
+    has_boundary_lower_y = false;
+    has_boundary_upper_y = false;
+  } else {
+    {
+      int mybndry = static_cast<int>(!(iterateBndryLowerY().isDone()));
+      int allbndry = 0;
+      mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend));
+      has_boundary_lower_y = static_cast<bool>(allbndry);
+    }
+    {
+      int mybndry = static_cast<int>(!(iterateBndryUpperY().isDone()));
+      int allbndry = 0;
+      mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart));
+      has_boundary_upper_y = static_cast<bool>(allbndry);
+    }
+  }
+
   output_info.write(_("\tdone\n"));
 
   return 0;
@@ -2218,9 +2223,9 @@ void BoutMesh::topology() {
     }
 
     for (int i = 0; i < limiter_count; ++i) {
-      int const yind = limiter_yinds[i];
-      int const xstart = limiter_xstarts[i];
-      int const xend = limiter_xends[i];
+      const int yind = limiter_yinds[i];
+      const int xstart = limiter_xstarts[i];
+      const int xend = limiter_xends[i];
       output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n",
                         yind, yind + 1, xstart, xend);
       add_target(yind, xstart, xend);
@@ -2928,6 +2933,9 @@ void BoutMesh::addBoundaryRegions() {
 }
 
 RangeIterator BoutMesh::iterateBndryLowerInnerY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
 
   int xs = 0;
   int xe = LocalNx - 1;
@@ -2963,6 +2971,9 @@ RangeIterator BoutMesh::iterateBndryLowerInnerY() const {
 }
 
 RangeIterator BoutMesh::iterateBndryLowerOuterY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
 
   int xs = 0;
   int xe = LocalNx - 1;
@@ -2997,6 +3008,10 @@ RangeIterator BoutMesh::iterateBndryLowerOuterY() const {
 }
 
 RangeIterator BoutMesh::iterateBndryLowerY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
+
   int xs = 0;
   int xe = LocalNx - 1;
   if ((DDATA_INDEST >= 0) && (DDATA_XSPLIT > xstart)) {
@@ -3026,6 +3041,10 @@ RangeIterator BoutMesh::iterateBndryLowerY() const {
 }
 
 RangeIterator BoutMesh::iterateBndryUpperInnerY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
+
   int xs = 0;
   int xe = LocalNx - 1;
 
@@ -3060,6 +3079,10 @@ RangeIterator BoutMesh::iterateBndryUpperInnerY() const {
 }
 
 RangeIterator BoutMesh::iterateBndryUpperOuterY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
+
   int xs = 0;
   int xe = LocalNx - 1;
 
@@ -3094,6 +3117,10 @@ RangeIterator BoutMesh::iterateBndryUpperOuterY() const {
 }
 
 RangeIterator BoutMesh::iterateBndryUpperY() const {
+  if (this->isFci()) {
+    throw BoutException("FCI should never use this iterator");
+  }
+
   int xs = 0;
   int xe = LocalNx - 1;
   if ((UDATA_INDEST >= 0) && (UDATA_XSPLIT > xstart)) {
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index fa965a335a..01dd6f72f8 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -62,6 +62,13 @@
 #include <string>
 #include <string_view>
 
+namespace {
+template <typename T>
+int sgn(T val) {
+  return (T(0) < val) - (val < T(0));
+}
+} // namespace
+
 using namespace std::string_view_literals;
 
 FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy,
@@ -183,6 +190,8 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy,
   BoutMask to_remove(map_mesh);
   const int xend = map_mesh->xstart
                    + ((map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE()) - 1;
+  // Default to the maximum number of points
+  const int defValid{map_mesh->ystart - 1 + std::abs(offset)};
   // Serial loop because call to BoundaryRegionPar::addPoint
   // (probably?) can't be done in parallel
   BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) {
@@ -252,11 +261,12 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy,
     // need at least 2 points in the domain.
     ASSERT2(map_mesh->xend - map_mesh->xstart >= 2);
     auto boundary = (xt_prime[i] < map_mesh->xstart) ? inner_boundary : outer_boundary;
-    boundary->add_point(x, y, z, x + dx, y + (0.5 * offset_),
-                        z + dz, // Intersection point in local index space
-                        0.5,    // Distance to intersection
-                        1       // Default to that there is a point in the other direction
-    );
+    if (!boundary->contains(x, y, z)) {
+      boundary->add_point(x, y, z, x + dx, y + offset - (sgn(offset) * 0.5),
+                          z + dz, // Intersection point in local index space
+                          std::abs(offset) - 0.5, // Distance to intersection
+                          defValid, offset);
+    }
   }
   region_no_boundary = region_no_boundary.mask(to_remove);
 
@@ -344,7 +354,6 @@ FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool
     field_line_maps.emplace_back(mesh, dy, options, -offset, backward_boundary_xin,
                                  backward_boundary_xout, zperiodic);
   }
-  ASSERT0(mesh.ystart == 1);
   const std::array bndries = {forward_boundary_xin, forward_boundary_xout,
                               backward_boundary_xin, backward_boundary_xout};
   for (const auto& bndry : bndries) {
@@ -352,9 +361,13 @@ FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool
       if (bndry->dir == bndry2->dir) {
         continue;
       }
-      for (bndry->first(); !bndry->isDone(); bndry->next()) {
-        if (bndry2->contains(*bndry)) {
-          bndry->setValid(0);
+      for (auto pnt : *bndry) {
+        for (auto pnt2 : *bndry2) {
+#warning this could likely be done faster
+          if (pnt.ind() == pnt2.ind()) {
+            pnt.setValid(
+                static_cast<signed char>(std::abs((pnt2.offset() - pnt.offset())) - 2));
+          }
         }
       }
     }
diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx
index c71618ab19..4b470da2ad 100644
--- a/src/mesh/parallel/shiftedmetricinterp.cxx
+++ b/src/mesh/parallel/shiftedmetricinterp.cxx
@@ -146,7 +146,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in,
           0.25
               * (1                                                     // dy/2
                  + dy(it.ind, mesh.yend + 1) / dy(it.ind, mesh.yend)), // length
-          yvalid);
+          yvalid, 1);
     }
   }
   auto backward_boundary_xin = std::make_shared<BoundaryRegionPar>(
@@ -162,7 +162,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in,
           0.25
               * (1 // dy/2
                  + dy(it.ind, mesh.ystart - 1) / dy(it.ind, mesh.ystart)),
-          yvalid);
+          yvalid, -1);
     }
   }
   // Create regions for parallel boundary conditions
@@ -179,7 +179,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in,
           0.25
               * (1 // dy/2
                  + dy(it.ind, mesh.yend + 1) / dy(it.ind, mesh.yend)),
-          yvalid);
+          yvalid, 1);
     }
   }
   auto backward_boundary_xout = std::make_shared<BoundaryRegionPar>(
@@ -195,7 +195,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in,
           0.25
               * (dy(it.ind, mesh.ystart - 1) / dy(it.ind, mesh.ystart) // dy/2
                  + 1),
-          yvalid);
+          yvalid, -1);
     }
   }
 
@@ -215,7 +215,7 @@ void ShiftedMetricInterp::checkInputGrid() {
                             "Should be 'orthogonal'.");
     }
   } // else: coordinate_system variable not found in grid input, indicates older input
-    //       file so must rely on the user having ensured the type is correct
+  //       file so must rely on the user having ensured the type is correct
 }
 
 /*!
diff --git a/src/mesh/parallel_boundary_op.cxx b/src/mesh/parallel_boundary_op.cxx
index ebd9852791..644331221e 100644
--- a/src/mesh/parallel_boundary_op.cxx
+++ b/src/mesh/parallel_boundary_op.cxx
@@ -1,21 +1,20 @@
 #include "bout/parallel_boundary_op.hxx"
+#include "bout/bout_types.hxx"
 #include "bout/constants.hxx"
 #include "bout/field_factory.hxx"
 #include "bout/globals.hxx"
 #include "bout/mesh.hxx"
 #include "bout/output.hxx"
+#include "bout/parallel_boundary_region.hxx"
 
-BoutReal BoundaryOpPar::getValue(const BoundaryRegionPar& bndry, BoutReal t) {
-  BoutReal value;
-
+BoutReal BoundaryOpPar::getValue(const BoundaryRegionParIter& bndry, BoutReal t) {
   switch (value_type) {
   case ValueType::GEN:
     return gen_values->generate(bout::generator::Context(
         bndry.s_x(), bndry.s_y(), bndry.s_z(), CELL_CENTRE, bndry.localmesh, t));
   case ValueType::FIELD:
     // FIXME: Interpolate to s_x, s_y, s_z...
-    value = (*field_values)[bndry.ind()];
-    return value;
+    return (*field_values)[bndry.ind()];
   case ValueType::REAL:
     return real_value;
   default:
diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx
index 8e3cfac2f7..0871577b86 100644
--- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx
+++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx
@@ -7,6 +7,8 @@
 #include "bout/output.hxx"
 #include "bout/parallel_boundary_region.hxx"
 
+#include "bout/yboundary_regions.hxx"
+
 #include <fmt/format.h>
 
 #include <vector>
@@ -24,17 +26,39 @@ int main(int argc, char** argv) {
     const auto boundary = static_cast<BoundaryParType>(i);
     const auto boundary_name = toString(boundary);
     mesh->communicate(fields[i]);
-    for (const auto& bndry_par : mesh->getBoundariesPar(boundary)) {
-      output.write("{:s} region\n", boundary_name);
-      for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) {
-        fields[i][bndry_par->ind()] += 1;
-        output.write("{:s} increment\n", boundary_name);
+    for (auto& bndry_par : mesh->getBoundariesPar(static_cast<BoundaryParType>(i))) {
+      output.write("{:s} region\n", toString(static_cast<BoundaryParType>(i)));
+      for (const auto& pnt : *bndry_par) {
+        fields[i][pnt.ind()] += 1;
+        output.write("{:s} increment\n", toString(static_cast<BoundaryParType>(i)));
       }
     }
     output.write("{:s} done\n", boundary_name);
 
     dump[fmt::format("field_{:s}", boundary_name)] = fields[i];
   }
+  for (const auto& name : {"forward_xt_prime", "backward_xt_prime"}) {
+    Field3D tmp;
+    mesh->get(tmp, name);
+    dump[name] = tmp;
+  }
+
+  {
+    Options dummy;
+    YBoundary ybndry;
+    ybndry.init(dummy, mesh);
+
+    std::vector<Field3D> fields((mesh->ystart * 2) + 1, Field3D{0.0});
+    for (auto& field : fields) {
+      field.allocate();
+    }
+    ybndry.iter_pnts(
+        [&](const auto& pnt) { fields[pnt.dir + mesh->ystart][pnt.ind()] += 1; });
+
+    for (int i = -mesh->ystart; i <= mesh->ystart; ++i) {
+      dump[fmt::format("ybndry_{}", i)] = fields[i + mesh->ystart];
+    }
+  }
 
   bout::writeDefaultOutputFile(dump);
 
diff --git a/tests/integrated/test-fci-boundary/runtest b/tests/integrated/test-fci-boundary/runtest
index e749055185..74a8abf741 100755
--- a/tests/integrated/test-fci-boundary/runtest
+++ b/tests/integrated/test-fci-boundary/runtest
@@ -4,9 +4,19 @@
 
 from boututils.run_wrapper import launch_safe
 from boututils.datafile import DataFile
-from boutdata.collect import collect
+from boutdata.collect import collect as _collect
 
-import numpy as np
+from numpy.testing import assert_allclose
+
+
+def collect(*args):
+    return _collect(
+        *args,
+        info=False,
+        path=directory,
+        xguards=False,
+        yguards=False,
+    )
 
 
 nprocs = [1]
@@ -29,27 +39,27 @@ regions = {
 }
 regions = {k: v.astype(int) for k, v in regions.items()}
 
+
 for x in "xout", "xin":
     regions[x] = regions[f"{x}_fwd"] + regions[f"{x}_bwd"]
 for x in "fwd", "bwd":
     regions[x] = regions[f"xin_{x}"] + regions[f"xout_{x}"]
 regions["all"] = regions["xin"] + regions["xout"]
 
+bndrys = {
+    "ybndry_-1": regions["xout_bwd"],
+    "ybndry_0": regions["xout_fwd"] * 0,
+    "ybndry_1": regions["xout_fwd"],
+}
+
 for nproc in nprocs:
     cmd = "./get_par_bndry"
     _, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True)
 
     for k, v in regions.items():
-        data = collect(
-            f"field_{k}",
-            info=False,
-            path=directory,
-            xguards=False,
-            yguards=False,
-        )
-        assert np.allclose(data, v), (
-            f"{k} does not match",
-            np.sum(data),
-            np.sum(v),
-            np.max(data),
-        )
+        data = collect(f"field_{k}")
+        assert_allclose(data, v)
+    for i in range(-1, 2):
+        name = f"ybndry_{i}"
+        data = collect(name)
+        assert_allclose(bndrys[name], data)