diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H
index f557403aa..ecac975fc 100644
--- a/Source/BoundaryConditions/ABLMost.H
+++ b/Source/BoundaryConditions/ABLMost.H
@@ -32,7 +32,8 @@ public:
     explicit ABLMost (const amrex::Vector<amrex::Geometry>& geom,
                       amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
-                      amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd)
+                      amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
+                      int ng)
       : m_geom(geom), m_ma(geom,vars_old,Theta_prim,z_phys_nd)
     {
         amrex::ParmParse pp("erf");
@@ -88,8 +89,7 @@ public:
         {
             // Z0 heights
             //--------------------------------------------------------
-            // TODO: generalize the "3" for the number of ghost cells
-            amrex::Box bx = amrex::grow(m_geom[lev].Domain(),3);
+            amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng);
             bx.setSmall(2,0);
             bx.setBig(2,0);
             z_0[lev].resize(bx,1);
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index d49de1435..bcb358dd8 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -532,7 +532,8 @@ ERF::InitData ()
     //       WritePlotFile calls FillPatch in order to compute gradients
     if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST)
     {
-        m_most = std::make_unique<ABLMost>(geom,vars_old,Theta_prim,z_phys_nd);
+        int ng_for_most = ComputeGhostCells(solverChoice.advChoice,solverChoice.use_NumDiff)+1;
+        m_most = std::make_unique<ABLMost>(geom,vars_old,Theta_prim,z_phys_nd,ng_for_most);
 
         // We now configure ABLMost params here so that we can print the averages at t=0
         // Note we don't fill ghost cells here because this is just for diagnostics