diff --git a/include/bout/boundary_standard.hxx b/include/bout/boundary_standard.hxx
index b1116e159f..3cb4b23803 100644
--- a/include/bout/boundary_standard.hxx
+++ b/include/bout/boundary_standard.hxx
@@ -33,6 +33,30 @@ private:
   BoutReal val;
 };
 
+/// Dirichlet (set to zero) boundary condition
+class BoundaryDirichlet_O1 : public BoundaryOp {
+public:
+  BoundaryDirichlet_O1() : gen(nullptr) {}
+  BoundaryDirichlet_O1(BoundaryRegion* region, std::shared_ptr<FieldGenerator> g)
+      : BoundaryOp(region), gen(std::move(g)) {}
+
+  using BoundaryOp::clone;
+  BoundaryOp* clone(BoundaryRegion* region, const std::list<std::string>& args) override;
+
+  using BoundaryOp::apply;
+  void apply(Field2D& f) override;
+  void apply(Field2D& f, BoutReal t) override;
+  void apply(Field3D& f) override;
+  void apply(Field3D& f, BoutReal t) override;
+
+  using BoundaryOp::apply_ddt;
+  void apply_ddt(Field2D& f) override;
+  void apply_ddt(Field3D& f) override;
+
+private:
+  std::shared_ptr<FieldGenerator> gen; // Generator
+};
+
 /// Dirichlet (set to zero) boundary condition
 class BoundaryDirichlet : public BoundaryOp {
 public:
@@ -163,6 +187,30 @@ public:
   void apply(Field3D& f) override;
 };
 
+/// Neumann (zero-gradient) boundary condition, using 1st order on boundary
+class BoundaryNeumann_O1 : public BoundaryOp {
+public:
+  BoundaryNeumann_O1() : gen(nullptr) {}
+  BoundaryNeumann_O1(BoundaryRegion* region, std::shared_ptr<FieldGenerator> g)
+      : BoundaryOp(region), gen(std::move(g)) {}
+
+  using BoundaryOp::clone;
+  BoundaryOp* clone(BoundaryRegion* region, const std::list<std::string>& args) override;
+
+  using BoundaryOp::apply;
+  void apply(Field2D& f) override;
+  void apply(Field2D& f, BoutReal t) override;
+  void apply(Field3D& f) override;
+  void apply(Field3D& f, BoutReal t) override;
+
+  using BoundaryOp::apply_ddt;
+  void apply_ddt(Field2D& f) override;
+  void apply_ddt(Field3D& f) override;
+
+private:
+  std::shared_ptr<FieldGenerator> gen;
+};
+
 /// Neumann boundary condition set half way between guard cell and grid cell at 2nd order accuracy
 class BoundaryNeumann_2ndOrder : public BoundaryOp {
 public:
diff --git a/src/mesh/boundary_factory.cxx b/src/mesh/boundary_factory.cxx
index f0c41cf8a1..3841a95b11 100644
--- a/src/mesh/boundary_factory.cxx
+++ b/src/mesh/boundary_factory.cxx
@@ -21,15 +21,18 @@ using std::string;
 BoundaryFactory* BoundaryFactory::instance = nullptr;
 
 BoundaryFactory::BoundaryFactory() {
-  add(new BoundaryDirichlet(), "dirichlet");
-  add(new BoundaryDirichlet(), "dirichlet_o2"); // Synonym for "dirichlet"
+  add(new BoundaryDirichlet(), "dirichlet");       // Default
+  add(new BoundaryDirichlet_O1(), "dirichlet_o1"); // Old implementation in v3
+  add(new BoundaryDirichlet(), "dirichlet_o2");    // Synonym for "dirichlet"
   add(new BoundaryDirichlet_O3(), "dirichlet_o3");
   add(new BoundaryDirichlet_O4(), "dirichlet_o4");
   add(new BoundaryDirichlet_4thOrder(), "dirichlet_4thorder");
-  add(new BoundaryNeumann(), "neumann");
-  add(new BoundaryNeumann(), "neumann_O2"); // Synonym for "neumann"
+
+  add(new BoundaryNeumann(), "neumann");       // Default
+  add(new BoundaryNeumann_O1(), "neumann_o1"); // Old implementation in v3
+  add(new BoundaryNeumann(), "neumann_o2");    // Synonym for "neumann"
   add(new BoundaryNeumann_4thOrder(), "neumann_4thorder");
-  add(new BoundaryNeumann_O4(), "neumann_O4");
+  add(new BoundaryNeumann_O4(), "neumann_o4");
   add(new BoundaryNeumannPar(), "neumannpar");
   add(new BoundaryNeumann_NonOrthogonal(), "neumann_nonorthogonal");
   add(new BoundaryRobin(), "robin");
diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx
index 141cab0a43..491b4a5b71 100644
--- a/src/mesh/boundary_standard.cxx
+++ b/src/mesh/boundary_standard.cxx
@@ -117,6 +117,133 @@ void verifyNumPoints(BoundaryRegion*, int) {}
 
 ///////////////////////////////////////////////////////////////
 
+BoundaryOp* BoundaryDirichlet_O1::clone(BoundaryRegion* region,
+                                        const std::list<std::string>& args) {
+  verifyNumPoints(region, 1);
+
+  std::shared_ptr<FieldGenerator> newgen;
+  if (!args.empty()) {
+    // First argument should be an expression
+    newgen = FieldFactory::get()->parse(args.front());
+  }
+  return new BoundaryDirichlet_O1(region, newgen);
+}
+
+void BoundaryDirichlet_O1::apply(Field2D& f) { BoundaryDirichlet_O1::apply(f, 0.); }
+
+void BoundaryDirichlet_O1::apply(Field2D& f, BoutReal t) {
+  // Set (at 1st order) the value at the grid cell to the guard cells.
+
+  Mesh* mesh = bndry->localmesh;
+  ASSERT1(mesh == f.getMesh());
+  bndry->first();
+
+  // Decide which generator to use
+  std::shared_ptr<FieldGenerator> fg = gen;
+  if (!fg) {
+    fg = f.getBndryGenerator(bndry->location);
+  }
+
+  BoutReal val = 0.0;
+
+  // Check for staggered grids
+
+  CELL_LOC loc = f.getLocation();
+  if (mesh->StaggerGrids) {
+    // Staggered
+    throw BoutException("dirichlet_o1 BC is not implementated for staggered grids.");
+
+  } else {
+    // Non-staggered, standard case
+    for (; !bndry->isDone(); bndry->next1d()) {
+
+      if (fg) {
+        val = fg->generate(Context(bndry, loc, t, mesh));
+      }
+      f(bndry->x, bndry->y) = val;
+
+      // Need to set second guard cell, as may be used for interpolation or upwinding derivatives
+      // This is not very efficient. Both boundary cells can be treated in one loop.
+      for (int i = 1; i < bndry->width; i++) {
+        int xi = bndry->x + i * bndry->bx;
+        int yi = bndry->y + i * bndry->by;
+        f(xi, yi) = val;
+      }
+    }
+  }
+}
+
+void BoundaryDirichlet_O1::apply(Field3D& f) { BoundaryDirichlet_O1::apply(f, 0.); }
+
+void BoundaryDirichlet_O1::apply(Field3D& f, BoutReal t) {
+  // Set (at 1st order) the value at the grid cell to the guard cells.
+
+  Mesh* mesh = bndry->localmesh;
+  ASSERT1(mesh == f.getMesh());
+  bndry->first();
+
+  // Decide which generator to use
+  std::shared_ptr<FieldGenerator> fg = gen;
+  if (!fg) {
+    fg = f.getBndryGenerator(bndry->location);
+  }
+
+  BoutReal val = 0.0;
+
+  // Check for staggered grids
+
+  CELL_LOC loc = f.getLocation();
+  if (mesh->StaggerGrids) {
+    // Staggered.
+    throw BoutException("dirichlet_o1 BC is not implementated for staggered grids.");
+
+  } else {
+    // Standard (non-staggered) case
+    for (; !bndry->isDone(); bndry->next1d()) {
+      for (int zk = mesh->zstart; zk <= mesh->zend; zk++) {
+        if (fg) {
+          val = fg->generate(Context(bndry, zk, loc, t, mesh));
+        }
+        f(bndry->x, bndry->y, zk) = val;
+      }
+
+      // This is not very efficient. Both boundary cells can be treated in one loop.
+      for (int i = 1; i < bndry->width; i++) {
+        // Set any other guard cells using the values on the cells
+        int xi = bndry->x + i * bndry->bx;
+        int yi = bndry->y + i * bndry->by;
+        for (int zk = mesh->zstart; zk <= mesh->zend; zk++) {
+          if (fg) {
+            val = fg->generate(Context(bndry, zk, loc, t, mesh));
+          }
+          f(xi, yi, zk) = val;
+        }
+      }
+    }
+  }
+}
+
+void BoundaryDirichlet_O1::apply_ddt(Field2D& f) {
+  Field2D* dt = f.timeDeriv();
+  for (bndry->first(); !bndry->isDone(); bndry->next()) {
+    (*dt)(bndry->x, bndry->y) = 0.; // Set time derivative to zero
+  }
+}
+
+void BoundaryDirichlet_O1::apply_ddt(Field3D& f) {
+  Mesh* mesh = bndry->localmesh;
+  ASSERT1(mesh == f.getMesh());
+  Field3D* dt = f.timeDeriv();
+
+  for (bndry->first(); !bndry->isDone(); bndry->next()) {
+    for (int z = mesh->zstart; z <= mesh->zend; z++) {
+      (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero
+    }
+  }
+}
+
+///////////////////////////////////////////////////////////////
+
 BoundaryOp* BoundaryDirichlet::clone(BoundaryRegion* region,
                                      const std::list<std::string>& args) {
   verifyNumPoints(region, 1);
@@ -1714,6 +1841,138 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) {
 
   ///////////////////////////////////////////////////////////////
 
+  BoundaryOp* BoundaryNeumann_O1::clone(BoundaryRegion * region,
+                                        const std::list<std::string>& args) {
+    verifyNumPoints(region, 1);
+    std::shared_ptr<FieldGenerator> newgen = nullptr;
+    if (!args.empty()) {
+      // First argument should be an expression
+      newgen = FieldFactory::get()->parse(args.front());
+    }
+    return new BoundaryNeumann_O1(region, newgen);
+  }
+
+  void BoundaryNeumann_O1::apply(Field2D & f) { BoundaryNeumann_O1::apply(f, 0.); }
+
+  void BoundaryNeumann_O1::apply(Field2D & f, BoutReal t) {
+    // Set (at 1st order) the gradient/value at the grid cell to the guard cells.
+
+#if not(BOUT_USE_METRIC_3D)
+    Mesh* mesh = bndry->localmesh;
+    ASSERT1(mesh == f.getMesh());
+    Coordinates* metric = f.getCoordinates();
+
+    bndry->first();
+
+    // Decide which generator to use
+    std::shared_ptr<FieldGenerator> fg = gen;
+    if (!fg) {
+      fg = f.getBndryGenerator(bndry->location);
+    }
+
+    BoutReal val = 0.0;
+
+    // Check for staggered grids
+
+    CELL_LOC loc = f.getLocation();
+    if (mesh->StaggerGrids) {
+      // Staggered.
+      throw BoutException("neumann_o1 BC is not implementated for staggered grids.");
+
+    } else {
+      // Non-staggered, standard case
+
+      for (bndry->first(); !bndry->isDone(); bndry->next1d()) {
+        BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y)
+                         + bndry->by * metric->dy(bndry->x, bndry->y);
+
+        if (fg) {
+          val = fg->generate(Context(bndry, loc, t, mesh));
+        }
+
+        f(bndry->x, bndry->y) =
+            f(bndry->x - bndry->bx, bndry->y - bndry->by) + delta * val;
+        if (bndry->width == 2) {
+          f(bndry->x + bndry->bx, bndry->y + bndry->by) =
+              f(bndry->x, bndry->y) + delta * val;
+        }
+      }
+    }
+#else
+  throw BoutException("Applying boundary condition 'neumann' to Field2D "
+                      "not compatible with 3D metrics in all cases.");
+#endif
+  }
+
+  void BoundaryNeumann_O1::apply(Field3D & f) { BoundaryNeumann_O1::apply(f, 0.); }
+
+  void BoundaryNeumann_O1::apply(Field3D & f, BoutReal t) {
+    Mesh* mesh = bndry->localmesh;
+    ASSERT1(mesh == f.getMesh());
+    Coordinates* metric = f.getCoordinates();
+
+    bndry->first();
+
+    // Decide which generator to use
+    std::shared_ptr<FieldGenerator> fg = gen;
+    if (!fg) {
+      fg = f.getBndryGenerator(bndry->location);
+    }
+
+    BoutReal val = 0.0;
+
+    // Check for staggered grids
+
+    CELL_LOC loc = f.getLocation();
+    if (mesh->StaggerGrids) {
+      // Staggered.
+      throw BoutException("neumann_o1 BC is not implementated for staggered grids.");
+
+    } else {
+      for (; !bndry->isDone(); bndry->next1d()) {
+#if BOUT_USE_METRIC_3D
+        for (int zk = mesh->zstart; zk <= mesh->zend; zk++) {
+          BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk)
+                           + bndry->by * metric->dy(bndry->x, bndry->y, zk);
+#else
+      BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y)
+                       + bndry->by * metric->dy(bndry->x, bndry->y);
+      for (int zk = mesh->zstart; zk <= mesh->zend; zk++) {
+#endif
+          if (fg) {
+            val = fg->generate(Context(bndry, zk, loc, t, mesh));
+          }
+          f(bndry->x, bndry->y, zk) =
+              f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) + delta * val;
+          if (bndry->width == 2) {
+            f(bndry->x + bndry->bx, bndry->y + bndry->by, zk) =
+                f(bndry->x, bndry->y, zk) + delta * val;
+          }
+        }
+      }
+    }
+  }
+
+  void BoundaryNeumann_O1::apply_ddt(Field2D & f) {
+    Field2D* dt = f.timeDeriv();
+    for (bndry->first(); !bndry->isDone(); bndry->next()) {
+      (*dt)(bndry->x, bndry->y) = 0.; // Set time derivative to zero
+    }
+  }
+
+  void BoundaryNeumann_O1::apply_ddt(Field3D & f) {
+    Mesh* mesh = bndry->localmesh;
+    ASSERT1(mesh == f.getMesh());
+    Field3D* dt = f.timeDeriv();
+    for (bndry->first(); !bndry->isDone(); bndry->next()) {
+      for (int z = mesh->zstart; z <= mesh->zend; z++) {
+        (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero
+      }
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////
+
   BoundaryOp* BoundaryNeumann::clone(BoundaryRegion * region,
                                      const std::list<std::string>& args) {
     verifyNumPoints(region, 1);