diff --git a/Source/ERF.H b/Source/ERF.H
index d98d2be63..d49fffe32 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -164,8 +164,8 @@ public:
     void WriteMyEBSurface ();
 #endif
 
-    // Compute the divergence of the momenta
-    void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Vector<amrex::MultiFab>& mom_mf,
+    // Compute the divergence -- whether EB, no-terrain, flat terrain or general terrain
+    void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> rho0_u_const,
                              amrex::Geometry const& geom_at_lev);
 
     // Project the velocities to be divergence-free -- this is only relevant if anelastic == 1
diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp
index 01a79c1cd..4b9ab1084 100644
--- a/Source/IO/ERF_Plotfile.cpp
+++ b/Source/IO/ERF_Plotfile.cpp
@@ -1,7 +1,6 @@
 #include <ERF_EOS.H>
 #include <ERF.H>
 #include "AMReX_Interp_3D_C.H"
-#include "AMReX_PlotFileUtil.H"
 #include "ERF_TerrainMetrics.H"
 #include "ERF_Constants.H"
 
@@ -368,7 +367,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
             u[0] = &(vars_new[lev][Vars::xvel]);
             u[1] = &(vars_new[lev][Vars::yvel]);
             u[2] = &(vars_new[lev][Vars::zvel]);
-            computeDivergence(dmf, u, geom[lev]);
+            compute_divergence (lev, dmf, u, geom[lev]);
             mf_comp += 1;
         }
 
diff --git a/Source/LinearSolvers/ERF_ComputeDivergence.cpp b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
index 08102e387..292b0a62d 100644
--- a/Source/LinearSolvers/ERF_ComputeDivergence.cpp
+++ b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
@@ -7,18 +7,13 @@ using namespace amrex;
  * Project the single-level velocity field to enforce incompressibility
  * Note that the level may or may not be level 0.
  */
-void ERF::compute_divergence (int lev, MultiFab& rhs, Vector<MultiFab>& mom_mf, Geometry const& geom_at_lev)
+void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMREX_SPACEDIM> rho0_u_const, Geometry const& geom_at_lev)
 {
     BL_PROFILE("ERF::compute_divergence()");
 
     bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);
 
-    auto dxInv = geom[lev].InvCellSizeArray();
-
-    Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
-    rho0_u_const[0] = &mom_mf[IntVars::xmom];
-    rho0_u_const[1] = &mom_mf[IntVars::ymom];
-    rho0_u_const[2] = &mom_mf[IntVars::zmom];
+    auto dxInv = geom_at_lev.InvCellSizeArray();
 
     // ****************************************************************************
     // Compute divergence which will form RHS
@@ -33,45 +28,37 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Vector<MultiFab>& mom_mf,
         for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
             Box bx = mfi.tilebox();
-            const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
-            const Array4<Real const>& rho0v_arr = mom_mf[IntVars::ymom].const_array(mfi);
-            const Array4<Real const>& rho0w_arr = mom_mf[IntVars::zmom].const_array(mfi);
+            const Array4<Real const>& rho0u_arr = rho0_u_const[0]->const_array(mfi);
+            const Array4<Real const>& rho0v_arr = rho0_u_const[1]->const_array(mfi);
+            const Array4<Real const>& rho0w_arr = rho0_u_const[2]->const_array(mfi);
             const Array4<Real      >&  rhs_arr = rhs.array(mfi);
 
-            Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
-            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
-                Real dz = stretched_dz_d_ptr[k];
-                rhs_arr(i,j,k) =  (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
-                                 +(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
-                                 +(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
-            });
-        } // mfi
-    }
-    else if (l_use_terrain) // terrain is not flat
-    {
-        //
-        // Note we compute the divergence using "rho0w" == Omega
-        //
-        for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
-        {
-            Box bx = mfi.tilebox();
-            const Array4<Real >& rho0u_arr = mom_mf[IntVars::xmom].array(mfi);
-            const Array4<Real >& rho0v_arr = mom_mf[IntVars::ymom].array(mfi);
-            const Array4<Real >& rho0w_arr = mom_mf[IntVars::zmom].array(mfi);
-            const Array4<Real      >&  rhs_arr = rhs.array(mfi);
+            if (SolverChoice::terrain_is_flat) { // flat terrain
+                Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                    Real dz = stretched_dz_d_ptr[k];
+                    rhs_arr(i,j,k) =  (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
+                                     +(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
+                                     +(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
+                });
+            } else { // terrain is not flat
 
-            const Array4<Real const>& ax_arr = ax[lev]->const_array(mfi);
-            const Array4<Real const>& ay_arr = ay[lev]->const_array(mfi);
-            const Array4<Real const>& dJ_arr = detJ_cc[lev]->const_array(mfi);
-            //
-            // az == 1 for terrain-fitted coordinates
-            //
-            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-            {
-                rhs_arr(i,j,k) =   ((ax_arr(i+1,j,k)*rho0u_arr(i+1,j,k) - ax_arr(i,j,k)*rho0u_arr(i,j,k)) * dxInv[0]
-                                   +(ay_arr(i,j+1,k)*rho0v_arr(i,j+1,k) - ay_arr(i,j,k)*rho0v_arr(i,j,k)) * dxInv[1]
-                                   +(                rho0w_arr(i,j,k+1) -               rho0w_arr(i,j,k)) * dxInv[2]) / dJ_arr(i,j,k);
-            });
+                //
+                // Note we compute the divergence using "rho0w" == Omega
+                //
+                const Array4<Real const>& ax_arr = ax[lev]->const_array(mfi);
+                const Array4<Real const>& ay_arr = ay[lev]->const_array(mfi);
+                const Array4<Real const>& dJ_arr = detJ_cc[lev]->const_array(mfi);
+                //
+                // az == 1 for terrain-fitted coordinates
+                //
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    rhs_arr(i,j,k) =   ((ax_arr(i+1,j,k)*rho0u_arr(i+1,j,k) - ax_arr(i,j,k)*rho0u_arr(i,j,k)) * dxInv[0]
+                                       +(ay_arr(i,j+1,k)*rho0v_arr(i,j+1,k) - ay_arr(i,j,k)*rho0v_arr(i,j,k)) * dxInv[1]
+                                       +(                rho0w_arr(i,j,k+1) -               rho0w_arr(i,j,k)) * dxInv[2]) / dJ_arr(i,j,k);
+                });
+            }
         } // mfi
 
     }
diff --git a/Source/LinearSolvers/ERF_PoissonSolve.cpp b/Source/LinearSolvers/ERF_PoissonSolve.cpp
index 34d3ea898..5aff21d30 100644
--- a/Source/LinearSolvers/ERF_PoissonSolve.cpp
+++ b/Source/LinearSolvers/ERF_PoissonSolve.cpp
@@ -72,7 +72,12 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // Compute divergence which will form RHS
     // Note that we replace "rho0w" with the contravariant momentum, Omega
     // ****************************************************************************
-    compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
+    Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
+    rho0_u_const[0] = &mom_mf[IntVars::xmom];
+    rho0_u_const[1] = &mom_mf[IntVars::ymom];
+    rho0_u_const[2] = &mom_mf[IntVars::zmom];
+
+    compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);
 
     Real rhsnorm = rhs[0].norm0();
 
@@ -212,7 +217,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     //
     if (mg_verbose > 0)
     {
-        compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
+        compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);
 
         Print() << "Max/L2 norm of divergence  after solve at level " << lev << " : " << rhs[0].norm0() << " " << rhs[0].norm2() << std::endl;