diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx
index cdb0a067ed..d4d2ffc859 100644
--- a/include/bout/fv_ops.hxx
+++ b/include/bout/fv_ops.hxx
@@ -10,8 +10,12 @@
 #include "bout/field3d.hxx"
 #include "bout/mesh.hxx"
 #include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner)
+#include "bout/region.hxx"
 #include "bout/vector2d.hxx"
 
+#include <map>
+#include <memory>
+
 namespace FV {
 /// Vorticity:
 /// \f[
@@ -185,5 +189,46 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 template <typename CellEdges = MC>
 Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f);
 
+/// Div ( a Grad_perp(f) )
+///
+/// This version uses pre-computed coefficient.  They are generated by
+/// zoidberg during grid generation.
+
+/// @param[in]  a            The diffusion coefficient
+/// @param[in]  f            The field to diffuse
+/// @param[out] flow_xlow    Flow at the lower X cell boundary
+///                          Already includes area factor * flux
+///                          This field is optional
+/// @param[out] flow_zlow    Flow at the lower Z cell boundary
+///                          Already includes area factor * flux
+///                          This field is optional
+/// @param[in]  upwinding    Use upwinding to to calculate the fluxes
+
+class Div_a_Grad_perp_FCI {
+public:
+  Field3D operator()(const Field3D& a, const Field3D& f, Field3D& low_xlow,
+                     Field3D& flow_zlow, bool upwinding) const;
+  Field3D operator()(const Field3D& a, const Field3D& f, bool upwinding) const;
+  Div_a_Grad_perp_FCI(Mesh& mesh);
+
+  static std::shared_ptr<Div_a_Grad_perp_FCI> create(Mesh* mesh, BoutReal rho_s0);
+
+private:
+  template <bool extra, bool upwinding>
+  Field3D operator()(const Field3D& a, const Field3D& f, Field3D* low_xlow,
+                     Field3D* flow_zlow) const;
+  Field3D fac_XX;
+  Field3D fac_XZ;
+  Field3D fac_ZX;
+  Field3D fac_ZZ;
+  Field3D volume;
+
+  template <bool upwinding>
+  BoutReal xflux(const Field3D& a, const Field3D& f, const Ind3D& i) const;
+  template <bool upwinding>
+  BoutReal zflux(const Field3D& a, const Field3D& f, const Ind3D& i) const;
+  static std::map<Mesh*, std::weak_ptr<Div_a_Grad_perp_FCI>> op_cache;
+};
+
 } // namespace FV
 #endif // BOUT_FV_OPS_H
diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx
index 66ce16e2ed..53a3290607 100644
--- a/src/mesh/boundary_standard.cxx
+++ b/src/mesh/boundary_standard.cxx
@@ -1857,10 +1857,10 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) {
 
   void BoundaryNeumann_O1::apply(Field2D & f) { BoundaryNeumann_O1::apply(f, 0.); }
 
+#if not(BOUT_USE_METRIC_3D)
   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());
     const Coordinates* metric = f.getCoordinates();
@@ -1901,6 +1901,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) {
     }
 
 #else
+void BoundaryNeumann_O1::apply([[maybe_unused]] Field2D& f, [[maybe_unused]] BoutReal t) {
   throw BoutException("Applying boundary condition 'neumann' to Field2D "
                       "not compatible with 3D metrics in all cases.");
 #endif
diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx
index afe8088523..c111ab585d 100644
--- a/src/mesh/fv_ops.cxx
+++ b/src/mesh/fv_ops.cxx
@@ -10,10 +10,15 @@
 #include "bout/field3d.hxx"
 #include "bout/globals.hxx"
 #include "bout/msg_stack.hxx"
+#include "bout/output.hxx"
 #include "bout/region.hxx"
 #include "bout/utils.hxx"
 #include <bout/yboundary_regions.hxx>
 
+#include <map>
+#include <math.h>
+#include <memory>
+
 namespace {
 template <class T>
 struct Slices {
@@ -678,4 +683,119 @@ template Field3D Div_par_fvv_heating<WENO3>(const Field3D& f_in, const Field3D&
                                             Field3D& flow_ylow, bool fixflux);
 template Field3D Div_a_Grad_perp_limit<WENO3>(const Field3D& a, const Field3D& g,
                                               const Field3D& f);
+
+std::map<Mesh*, std::weak_ptr<Div_a_Grad_perp_FCI>> Div_a_Grad_perp_FCI::op_cache;
+std::shared_ptr<Div_a_Grad_perp_FCI> Div_a_Grad_perp_FCI::create(Mesh* mesh,
+                                                                 BoutReal rho_s0) {
+  if (auto search = Div_a_Grad_perp_FCI::op_cache.find(mesh);
+      search != Div_a_Grad_perp_FCI::op_cache.end()) {
+    if (std::shared_ptr<Div_a_Grad_perp_FCI> spt = search->second.lock()) {
+      return spt;
+    }
+    output_warn.write("\tPreformance Warning: DIV_A_GRAD_PERP needs to be re-created!\n");
+  }
+  auto dagp = std::make_shared<Div_a_Grad_perp_FCI>(*mesh);
+  dagp->volume /= rho_s0 * rho_s0;
+  Div_a_Grad_perp_FCI::op_cache[mesh] = dagp;
+  return dagp;
+}
+
+Div_a_Grad_perp_FCI::Div_a_Grad_perp_FCI(Mesh& mesh)
+    : fac_XX(&mesh), fac_XZ(&mesh), fac_ZX(&mesh), fac_ZZ(&mesh), volume(&mesh) {
+  ASSERT0(mesh.get(fac_XX, "dagp_fv_XX", 0.0, false) == 0);
+  ASSERT0(mesh.get(fac_XZ, "dagp_fv_XZ", 0.0, false) == 0);
+  ASSERT0(mesh.get(fac_ZX, "dagp_fv_ZX", 0.0, false) == 0);
+  ASSERT0(mesh.get(fac_ZZ, "dagp_fv_ZZ", 0.0, false) == 0);
+  ASSERT0(volume.hasParallelSlices() == false);
+  ASSERT0(mesh.get(volume, "dagp_fv_volume", 0.0, false) == 0);
+  ASSERT0(volume.hasParallelSlices() == false);
+  volume.setRegion("RGN_NOBNDRY");
+  if (!mesh.hasRegion3D("RGN_dapg_fv_xbndry")) {
+    mesh.addRegion("RGN_dapg_fv_xbndry",
+                   Region<>(mesh.xstart - 1, mesh.xstart - 1, mesh.ystart, mesh.yend,
+                            mesh.zstart, mesh.zend, mesh.LocalNy, mesh.LocalNz));
+  }
+}
+Field3D Div_a_Grad_perp_FCI::operator()(const Field3D& a, const Field3D& f,
+                                        Field3D& flow_xlow, Field3D& flow_zlow,
+                                        bool upwinding) const {
+  if (upwinding) {
+    return operator()<true, true>(a, f, &flow_xlow, &flow_zlow);
+  }
+  return operator()<true, false>(a, f, &flow_xlow, &flow_zlow);
+}
+Field3D Div_a_Grad_perp_FCI::operator()(const Field3D& a, const Field3D& f,
+                                        bool upwinding) const {
+  if (upwinding) {
+    return operator()<false, true>(a, f, nullptr, nullptr);
+  }
+  return operator()<false, false>(a, f, nullptr, nullptr);
+}
+template <bool extra, bool upwinding>
+Field3D Div_a_Grad_perp_FCI::operator()(const Field3D& a, const Field3D& f,
+                                        Field3D* flow_xlow, Field3D* flow_zlow) const {
+  auto result{zeroFrom(f)};
+  ASSERT1_FIELDS_COMPATIBLE(a, f);
+  if constexpr (extra) {
+    flow_xlow->allocate();
+    flow_zlow->allocate();
+  }
+  BOUT_FOR(i, f.getRegion("RGN_dapg_fv_xbndry")) {
+    const auto xf = xflux<upwinding>(a, f, i);
+    result[i.xp()] = xf;
+    if constexpr (extra) {
+      (*flow_xlow)[i.xp()] = xf;
+    }
+  }
+  BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) {
+    const auto xf = xflux<upwinding>(a, f, i);
+    result[i.xp()] = xf;
+    result[i] -= xf;
+    if constexpr (extra) {
+      (*flow_xlow)[i.xp()] = xf;
+    }
+    const auto zf = zflux<upwinding>(a, f, i);
+    result[i.zp()] += zf;
+    result[i] -= zf;
+    if constexpr (extra) {
+      (*flow_zlow)[i.zp()] = zf;
+    }
+  }
+  result.setRegion("RGN_NOBNDRY");
+  if constexpr (extra) {
+    flow_xlow->setRegion("RGN_NOBNDRY");
+    flow_zlow->setRegion("RGN_NOBNDRY");
+    *flow_xlow /= volume;
+    *flow_zlow /= volume;
+  }
+  return result / volume;
+}
+template <bool upwinding>
+inline BoutReal Div_a_Grad_perp_FCI::xflux(const Field3D& a, const Field3D& f,
+                                           const Ind3D& i) const {
+  const auto ixp = i.xp();
+  const auto dx = f[ixp] - f[i];
+  BoutReal av = BoutNaN;
+  if constexpr (upwinding) {
+    av = dx > 0 ? a[ixp] : a[i];
+  } else {
+    av = 0.5 * (a[i] + a[ixp]);
+  };
+  const auto dz = 0.5 * (f[i.zp()] - f[i.zm()] + f[i.zp().xp()] - f[i.zm().xp()]);
+  return -((fac_XX[i] * dx) + (fac_XZ[i] * dz)) * av;
+}
+template <bool upwinding>
+inline BoutReal Div_a_Grad_perp_FCI::zflux(const Field3D& a, const Field3D& f,
+                                           const Ind3D& i) const {
+  const auto izp = i.zp();
+  const auto dz = f[izp] - f[i];
+  BoutReal av = BoutNaN;
+  if constexpr (upwinding) {
+    av = dz > 0 ? a[izp] : a[i];
+  } else {
+    av = 0.5 * (a[i] + a[izp]);
+  }
+  const auto dx = 0.5 * (f[i.xp()] - f[i.xm()] + f[izp.xp()] - f[izp.xm()]);
+  return -((fac_ZX[i] * dx) + (fac_ZZ[i] * dz)) * av;
+}
 } // Namespace FV
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 4384440c74..0c0ddf706c 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -97,6 +97,7 @@ set(serial_tests_source
     ./mesh/test_boutmesh.cxx
     ./mesh/test_coordinates.cxx
     ./mesh/test_coordinates_accessor.cxx
+    ./mesh/test_fv_ops.cxx
     ./mesh/test_interpolation.cxx
     ./mesh/test_invert3x3.cxx
     ./mesh/test_mesh.cxx
diff --git a/tests/unit/mesh/test_fv_ops.cxx b/tests/unit/mesh/test_fv_ops.cxx
new file mode 100644
index 0000000000..2aecd707c4
--- /dev/null
+++ b/tests/unit/mesh/test_fv_ops.cxx
@@ -0,0 +1,35 @@
+#include "gtest/gtest.h"
+#include <cmath>
+
+#include "bout/build_defines.hxx"
+#include "bout/constants.hxx"
+#include "bout/coordinates.hxx"
+#include "bout/fv_ops.hxx"
+#include "bout/mesh.hxx"
+#include "bout/output.hxx"
+
+#include "fake_mesh.hxx"
+#include "fake_mesh_fixture.hxx"
+#include "test_extras.hxx"
+
+using bout::globals::mesh;
+
+using FVOpsTest = FakeMeshFixture;
+
+TEST_F(FVOpsTest, Div_a_Grad_perp_FCI_no_coeffs) {
+  EXPECT_THROW(FV::Div_a_Grad_perp_FCI::create(mesh, 1), BoutException);
+}
+
+TEST_F(FVOpsTest, Div_a_Grad_perp_FCI_coeffs) {
+  static_cast<FakeMesh*>(bout::globals::mesh)
+      ->setGridDataSource(new FakeGridDataSource({{"dagp_fv_XX", 1.0},
+                                                  {"dagp_fv_XZ", 1.0},
+                                                  {"dagp_fv_ZX", 1.0},
+                                                  {"dagp_fv_ZZ", 1.0},
+                                                  {"dagp_fv_volume", 1.0}}));
+
+  auto op = FV::Div_a_Grad_perp_FCI::create(mesh, 1);
+  Field3D a{1.0}, f{1.0};
+  (*op)(a, f, true);
+  (*op)(a, f, false);
+}