diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
index 3e7349b88..865db5f64 100644
--- a/Source/DataStructs/ERF_DataStruct.H
+++ b/Source/DataStructs/ERF_DataStruct.H
@@ -29,6 +29,10 @@ AMREX_ENUM(SubsteppingType,
     None, Explicit, Implicit
 );
 
+AMREX_ENUM(MeshType,
+    ConstantDz, StretchedDz, VariableDz
+);
+
 AMREX_ENUM(TerrainType,
     None, Static, Moving
 );
@@ -87,6 +91,9 @@ struct SolverChoice {
                                       "The grid stretching ratio must be greater than 1");
         }
         if (grid_stretching_ratio >= 1) {
+            if (mesh_type == MeshType::ConstantDz) {
+                mesh_type    = MeshType::StretchedDz;
+            }
             if (terrain_type != TerrainType::Static) {
                 amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl;
                 terrain_type = TerrainType::Static;
@@ -148,10 +155,20 @@ struct SolverChoice {
 
         // Is the terrain none, static or moving?
         pp.query_enum_case_insensitive("terrain_type",terrain_type);
+
+        if (terrain_type != TerrainType::None) {
+            mesh_type = MeshType::VariableDz;
+        }
+
         int n_zlevels = pp.countval("terrain_z_levels");
-        if (n_zlevels > 0 and terrain_type == TerrainType::None)
+        if (n_zlevels > 0)
         {
-            terrain_type = TerrainType::Static;
+            if (terrain_type == TerrainType::None) {
+                terrain_type = TerrainType::Static;
+            }
+            if (mesh_type == MeshType::ConstantDz) {
+                mesh_type    = MeshType::StretchedDz;
+            }
         }
 
         // Use lagged_delta_rt in the fast integrator?
@@ -434,6 +451,17 @@ struct SolverChoice {
             amrex::Print() << "    None" << std::endl;
         }
 
+        amrex::Print() << "   Mesh Type: " << std::endl;
+        if (mesh_type == MeshType::ConstantDz) {
+            amrex::Print() << "    ConstantDz" << std::endl;
+        } else if (mesh_type == MeshType::StretchedDz) {
+            amrex::Print() << "    StretchedDz" << std::endl;
+        } else if (mesh_type == MeshType::VariableDz) {
+            amrex::Print() << "    VariableDz" << std::endl;
+        } else {
+            amrex::Abort("No mesh_type set!");
+        }
+
         amrex::Print() << "ABL Driver Type: " << std::endl;
         if (abl_driver_type == ABLDriverType::None) {
             amrex::Print() << "    None" << std::endl;
@@ -533,6 +561,9 @@ struct SolverChoice {
     inline static
     TerrainType terrain_type = TerrainType::None;
 
+    inline static
+    MeshType mesh_type = MeshType::ConstantDz;
+
     static
     void set_flat_terrain_flag ()
     {
diff --git a/Source/ERF.H b/Source/ERF.H
index 89163e184..2523005a5 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -826,8 +826,10 @@ private:
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q1fx1_lev, SFS_q1fx2_lev, SFS_q1fx3_lev;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q2fx3_lev;
 
-    // Terrain / grid stretching
+    // Grid stretching
     amrex::Vector<amrex::Vector<amrex::Real>> zlevels_stag; // nominal height levels
+
+    // Terrain data structures
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_cc;
 
@@ -850,6 +852,7 @@ private:
 
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_t_rk;
 
+    // Map scale factors
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_m;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_u;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_v;
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index 3c759adcf..38c87fea4 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -143,7 +143,7 @@ ERF::ERF_shared ()
     const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
     const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
 
-    // This is only used when we have flat terrain and stretched z levels.
+    // This is only used when we have mesh_type == MeshType::StretchedDz
     stretched_dz_h.resize(nlevs_max);
     stretched_dz_d.resize(nlevs_max);
 
@@ -159,7 +159,8 @@ ERF::ERF_shared ()
                  solverChoice.zsurf,
                  solverChoice.dz0);
 
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type == MeshType::StretchedDz ||
+        SolverChoice::mesh_type == MeshType::VariableDz) {
         int nz = geom[0].Domain().length(2) + 1; // staggered
         if (std::fabs(zlevels_stag[0][nz-1]-geom[0].ProbHi(2)) > 1.0e-4) {
             Print() << "Note: prob_hi[2]=" << geom[0].ProbHi(2)
@@ -448,7 +449,6 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
 
     if (solverChoice.coupling_type == CouplingType::TwoWay)
     {
-        bool use_terrain = (SolverChoice::terrain_type != TerrainType::None);
         int ncomp = vars_new[0][Vars::cons].nComp();
         for (int lev = finest_level-1; lev >= 0; lev--)
         {
@@ -459,16 +459,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
                 const Box& bx = mfi.tilebox();
                 const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
                 const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
-                if (use_terrain) {
-                    const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
+                if (solverChoice.mesh_type == MeshType::ConstantDz) {
                     ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
-                        cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
+                        cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
                     });
                 } else {
+                    const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
                     ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
-                        cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
+                        cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
                     });
                 }
             } // mfi
@@ -482,16 +482,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
                 const Box& bx = mfi.tilebox();
                 const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
                 const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
-                if (use_terrain) {
-                    const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
+                if (solverChoice.mesh_type == MeshType::ConstantDz) {
                     ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
-                        cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k);
+                        cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
                     });
                 } else {
+                    const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
                     ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
-                        cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
+                        cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k);
                     });
                 }
             } // mfi
@@ -633,16 +633,16 @@ void
 ERF::InitData_post ()
 {
     if (restart_chkfile.empty()) {
-        if (SolverChoice::terrain_type != TerrainType::None) {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz) {
             if (init_type == InitType::Ideal) {
-                Abort("We do not currently support init_type = ideal with terrain");
+                Abort("We do not currently support init_type = ideal with non-constant dz");
             }
         }
 
         //
         // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one
         //
-        if (SolverChoice::terrain_type != TerrainType::None) {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz) {
             for (int crse_lev = finest_level-1; crse_lev >= 0; crse_lev--) {
                 average_down(  *detJ_cc[crse_lev+1],   *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev));
                 average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev));
@@ -704,8 +704,8 @@ ERF::InitData_post ()
         {
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 0,
                 "Thin immersed body with refinement not currently supported.");
-            if (SolverChoice::terrain_type != TerrainType::None) {
-                amrex::Print() << "NOTE: Thin immersed body with terrain has not been tested." << std::endl;
+            if (SolverChoice::mesh_type != MeshType::ConstantDz) {
+                amrex::Print() << "NOTE: Thin immersed body with non-constant dz has not been tested." << std::endl;
             }
         }
 
@@ -782,8 +782,8 @@ ERF::InitData_post ()
                                                  h_v_geos[lev], d_v_geos[lev],
                                                  geom[lev], z_phys_cc[lev]);
             } else {
-                if (SolverChoice::terrain_type != TerrainType::None && !SolverChoice::terrain_is_flat) {
-                    amrex::Print() << "Note: 1-D geostrophic wind profile input is only defined for grid stretching, not real terrain" << std::endl;
+                if (SolverChoice::mesh_type == MeshType::VariableDz) {
+                    amrex::Print() << "Note: 1-D geostrophic wind profile input is not defined for real terrain" << std::endl;
                 }
                 init_geo_wind_profile(solverChoice.abl_geo_wind_table,
                                       h_u_geos[lev], d_u_geos[lev],
@@ -875,7 +875,6 @@ ERF::InitData_post ()
     //
     if (restart_chkfile == "")
     {
-        // Note -- this projection is only defined for no terrain
         if (solverChoice.project_initial_velocity) {
             Real dummy_dt = 1.0;
             for (int lev = 0; lev <= finest_level; ++lev)
@@ -917,7 +916,7 @@ ERF::InitData_post ()
     for (int lev = 0; lev <= finest_level; ++lev)
     {
         dz_min[lev] = geom[lev].CellSize(2);
-        if ( SolverChoice::terrain_type != TerrainType::None ) {
+        if ( SolverChoice::mesh_type != MeshType::ConstantDz ) {
             dz_min[lev] *= (*detJ_cc[lev]).min(0);
         }
     }
diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp
index c0874e153..40a67bbb5 100644
--- a/Source/ERF_MakeNewArrays.cpp
+++ b/Source/ERF_MakeNewArrays.cpp
@@ -42,7 +42,8 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     // ********************************************************************************************
     // Allocate terrain arrays
     // ********************************************************************************************
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type == MeshType::StretchedDz ||
+        SolverChoice::mesh_type == MeshType::VariableDz) {
         z_phys_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
 
         if (solverChoice.terrain_type == TerrainType::Moving)
@@ -452,7 +453,8 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
 void
 ERF::init_zphys (int lev, Real time)
 {
-    if (SolverChoice::terrain_type != TerrainType::None)
+    if (SolverChoice::mesh_type == MeshType::StretchedDz ||
+        SolverChoice::mesh_type == MeshType::VariableDz)
     {
         if (init_type != InitType::Real && init_type != InitType::Metgrid)
         {
@@ -490,7 +492,7 @@ ERF::init_zphys (int lev, Real time)
 void
 ERF::remake_zphys (int lev, std::unique_ptr<MultiFab>& temp_zphys_nd)
 {
-    if (lev > 0 && SolverChoice::terrain_type != TerrainType::None)
+    if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz)
     {
         //
         // First interpolate from coarser level
@@ -517,7 +519,8 @@ ERF::remake_zphys (int lev, std::unique_ptr<MultiFab>& temp_zphys_nd)
 void
 ERF::update_terrain_arrays (int lev)
 {
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type == MeshType::StretchedDz ||
+        SolverChoice::mesh_type == MeshType::VariableDz) {
         make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
         make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);
         make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
@@ -549,7 +552,7 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf)
 void
 ERF::make_physbcs (int lev)
 {
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type == MeshType::VariableDz) {
         AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr);
     }
 
diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp
index c524b2434..fdcfd4f4b 100644
--- a/Source/ERF_MakeNewLevel.cpp
+++ b/Source/ERF_MakeNewLevel.cpp
@@ -210,7 +210,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
     //
     // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one
     //
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type != MeshType::ConstantDz) {
         for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) {
             average_down(  *detJ_cc[crse_lev+1],   *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev));
             average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev));
@@ -339,7 +339,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     //
     // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one
     //
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type != MeshType::ConstantDz) {
         for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) {
             average_down(  *detJ_cc[crse_lev+1],   *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev));
             average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev));
diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp
index e6d839911..db63f88eb 100644
--- a/Source/ERF_Tagging.cpp
+++ b/Source/ERF_Tagging.cpp
@@ -161,7 +161,7 @@ ERF::refinement_criteria_setup ()
                     jlo = static_cast<int>((rbox_lo[1] - plo[1])/dx[1]);
                     ihi = static_cast<int>((rbox_hi[0] - plo[0])/dx[0]-1);
                     jhi = static_cast<int>((rbox_hi[1] - plo[1])/dx[1]-1);
-                    if (SolverChoice::terrain_type != TerrainType::None) {
+                    if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                         // Search for k indices corresponding to nominal grid
                         // AGL heights
                         const Box& domain = geom[lev_for_box].Domain();
diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp
index 72c3a9da1..d2294277f 100644
--- a/Source/IO/ERF_Checkpoint.cpp
+++ b/Source/IO/ERF_Checkpoint.cpp
@@ -145,7 +145,7 @@ ERF::WriteCheckpointFile () const
         MultiFab::Copy(base,base_state[lev],0,0,ncomp_base,ng_base);
         VisMF::Write(base, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState"));
 
-        if (SolverChoice::terrain_type != TerrainType::None)  {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz)  {
             // Note that we also write the ghost cells of z_phys_nd
             IntVect ng = z_phys_nd[lev]->nGrowVect();
             MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng);
@@ -443,7 +443,7 @@ ERF::ReadCheckpointFile ()
         }
         base_state[lev].FillBoundary(geom[lev].periodicity());
 
-        if (SolverChoice::terrain_type != TerrainType::None)  {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz)  {
            // Note that we also read the ghost cells of z_phys_nd
            IntVect ng = z_phys_nd[lev]->nGrowVect();
            MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng);
diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp
index 4b9ab1084..6378a9bbd 100644
--- a/Source/IO/ERF_Plotfile.cpp
+++ b/Source/IO/ERF_Plotfile.cpp
@@ -72,7 +72,8 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     //
     for (int i = 0; i < derived_names.size(); ++i) {
         if ( containerHasElement(plot_var_names, derived_names[i]) ) {
-            if (SolverChoice::terrain_type != TerrainType::None || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) {
+            if ( SolverChoice::mesh_type != MeshType::ConstantDz ||
+                (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) {
                 if ( (solverChoice.moisture_type == MoistureType::SAM ||
                       solverChoice.moisture_type == MoistureType::SAM_NoIce) ||
                      (derived_names[i] != "qi" &&
@@ -225,7 +226,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 
     // Vector of MultiFabs for nodal data
     Vector<MultiFab> mf_nd(finest_level+1);
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if ( SolverChoice::mesh_type != MeshType::ConstantDz) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes();
             mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0);
@@ -578,7 +579,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
                 const Array4<Real>& derdat = mf[lev].array(mfi);
                 const Array4<Real> & p_arr  = pres.array(mfi);
 
-                if (SolverChoice::terrain_type != TerrainType::None) {
+                if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                     const Array4<Real const>& z_nd = z_phys_nd[lev]->const_array(mfi);
 
                     ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -672,7 +673,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
                 const Array4<Real>& derdat = mf[lev].array(mfi);
                 const Array4<Real> & p_arr  = pres.array(mfi);
 
-                if (SolverChoice::terrain_type != TerrainType::None) {
+                if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                     const Array4<Real const>& z_nd = z_phys_nd[lev]->const_array(mfi);
 
                     ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -843,7 +844,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
             mf_comp += 1;
         } // pres_hse_y
 
-        if (SolverChoice::terrain_type != TerrainType::None) {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz) {
             if (containerHasElement(plot_var_names, "z_phys"))
             {
                 MultiFab::Copy(mf[lev],*z_phys_cc[lev],0,mf_comp,1,0);
@@ -1350,7 +1351,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 #endif
 
     // Fill terrain distortion MF
-    if (SolverChoice::terrain_type != TerrainType::None) {
+    if (SolverChoice::mesh_type != MeshType::ConstantDz) {
         for (int lev(0); lev <= finest_level; ++lev) {
             MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0);
             Real dz = Geom()[lev].CellSizeArray()[2];
@@ -1399,7 +1400,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
         if (plotfile_type == PlotFileType::Amrex)
         {
             Print() << "Writing native plotfile " << plotfilename << "\n";
-            if (SolverChoice::terrain_type != TerrainType::None) {
+            if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                 WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
                                                    GetVecOfConstPtrs(mf),
                                                    GetVecOfConstPtrs(mf_nd),
@@ -1510,7 +1511,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
                 }
 
                Print() << "Writing plotfile " << plotfilename << "\n";
-               if (SolverChoice::terrain_type != TerrainType::None) {
+               if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                    WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
                                                       GetVecOfConstPtrs(mf2),
                                                       GetVecOfConstPtrs(mf_nd),
@@ -1523,7 +1524,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
                }
 
             } else {
-                if (SolverChoice::terrain_type != TerrainType::None) {
+               if (SolverChoice::mesh_type != MeshType::ConstantDz) {
                     WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
                                                        GetVecOfConstPtrs(mf),
                                                        GetVecOfConstPtrs(mf_nd),
diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp
index d839640ab..a8ef59498 100644
--- a/Source/IO/ERF_WriteScalarProfiles.cpp
+++ b/Source/IO/ERF_WriteScalarProfiles.cpp
@@ -406,7 +406,7 @@ ERF::volWgtSumMF (int lev,
     auto const& dx = geom[lev].CellSizeArray();
     Real cell_vol = dx[0]*dx[1]*dx[2];
     volume.setVal(cell_vol);
-    if (solverChoice.terrain_type != TerrainType::None) {
+    if (solverChoice.mesh_type != MeshType::ConstantDz) {
         MultiFab::Multiply(volume, *detJ_cc[lev], 0, 0, 1, 0);
     }
     sum = MultiFab::Dot(tmp, 0, volume, 0, 1, 0, local);
diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp
index 6af99e5c7..c1ad1e081 100644
--- a/Source/Initialization/ERF_Init1D.cpp
+++ b/Source/Initialization/ERF_Init1D.cpp
@@ -95,7 +95,7 @@ ERF::initHSE (int lev)
 
         std::unique_ptr<MultiFab> new_z_phys_cc;
         std::unique_ptr<MultiFab> new_z_phys_nd;
-        if (solverChoice.terrain_type != TerrainType::None) {
+        if (solverChoice.mesh_type != MeshType::ConstantDz) {
             new_z_phys_cc = std::make_unique<MultiFab>(ba_new,dm_new,1,1);
             new_z_phys_cc->ParallelCopy(*z_phys_cc[lev],0,0,1,1,1);
 
@@ -151,7 +151,7 @@ ERF::erf_enforce_hse (int lev,
                       std::unique_ptr<MultiFab>& z_cc)
 {
     Real l_gravity = solverChoice.gravity;
-    bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);
+    bool l_use_terrain = (solverChoice.mesh_type != MeshType::ConstantDz);
 
     const auto geomdata = geom[lev].data();
     const Real dz = geomdata.CellSize(2);
diff --git a/Source/LinearSolvers/ERF_ComputeDivergence.cpp b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
index 292b0a62d..fbbeb133e 100644
--- a/Source/LinearSolvers/ERF_ComputeDivergence.cpp
+++ b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
@@ -11,8 +11,6 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
 {
     BL_PROFILE("ERF::compute_divergence()");
 
-    bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);
-
     auto dxInv = geom_at_lev.InvCellSizeArray();
 
     // ****************************************************************************
@@ -23,7 +21,11 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
     bool already_on_centroids = true;
     EB_computeDivergence(rhs, rho0_u_const, geom_at_lev, already_on_centroids);
 #else
-    if (l_use_terrain && SolverChoice::terrain_is_flat)
+    if (SolverChoice::mesh_type == MeshType::ConstantDz)
+    {
+        computeDivergence(rhs, rho0_u_const, geom_at_lev);
+    }
+    else
     {
         for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
@@ -33,7 +35,7 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
             const Array4<Real const>& rho0w_arr = rho0_u_const[2]->const_array(mfi);
             const Array4<Real      >&  rhs_arr = rhs.array(mfi);
 
-            if (SolverChoice::terrain_is_flat) { // flat terrain
+            if (SolverChoice::mesh_type == MeshType::StretchedDz) {
                 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];
@@ -41,7 +43,7 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
                                      +(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
+            } else {
 
                 //
                 // Note we compute the divergence using "rho0w" == Omega
@@ -60,11 +62,6 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
                 });
             }
         } // mfi
-
-    }
-    else // no terrain
-    {
-        computeDivergence(rhs, rho0_u_const, geom_at_lev);
     }
 #endif
 }
diff --git a/Source/LinearSolvers/ERF_PoissonSolve.cpp b/Source/LinearSolvers/ERF_PoissonSolve.cpp
index 5aff21d30..7c9662e9b 100644
--- a/Source/LinearSolvers/ERF_PoissonSolve.cpp
+++ b/Source/LinearSolvers/ERF_PoissonSolve.cpp
@@ -16,8 +16,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     auto const dom_hi = ubound(geom[lev].Domain());
 #endif
 
-    bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None;
-
     // Make sure the solver only sees the levels over which we are solving
     Vector<BoxArray>            ba_tmp;   ba_tmp.push_back(mom_mf[Vars::cons].boxArray());
     Vector<DistributionMapping> dm_tmp;   dm_tmp.push_back(mom_mf[Vars::cons].DistributionMap());
@@ -44,7 +42,8 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // ****************************************************************************
     //
 #ifndef ERF_USE_EB
-    if (l_use_terrain && !SolverChoice::terrain_is_flat) {
+    if (solverChoice.mesh_type == MeshType::VariableDz)
+    {
         for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
             const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
@@ -133,7 +132,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // ****************************************************************************
     // No terrain or grid stretching
     // ****************************************************************************
-    if (!l_use_terrain) {
+    if (solverChoice.mesh_type == MeshType::ConstantDz) {
 #ifdef ERF_USE_FFT
         if (use_fft) {
             if (boxes_make_rectangle) {
@@ -156,7 +155,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // ****************************************************************************
     // Grid stretching (flat terrain)
     // ****************************************************************************
-    else if (l_use_terrain && SolverChoice::terrain_is_flat) {
+    else if (solverChoice.mesh_type == MeshType::StretchedDz) {
 #ifndef ERF_USE_FFT
         amrex::Abort("Rebuild with USE_FFT = TRUE so you can use the FFT solver");
 #else
@@ -174,7 +173,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // ****************************************************************************
     // General terrain
     // ****************************************************************************
-    else if (l_use_terrain && !SolverChoice::terrain_is_flat) {
+    else if (solverChoice.mesh_type == MeshType::VariableDz) {
 #ifdef ERF_USE_FFT
         if (use_fft)
         {
@@ -243,7 +242,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // Now convert the rho0w MultiFab back to holding (rho0w) rather than Omega
     // ****************************************************************************
     //
-    if (l_use_terrain && !solverChoice.terrain_is_flat)
+    if (solverChoice.mesh_type == MeshType::VariableDz)
     {
         for (MFIter mfi(mom_mf[Vars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
diff --git a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp
index 100fdf75b..107f6c0b7 100644
--- a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp
+++ b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp
@@ -20,7 +20,7 @@ bool ERF::projection_has_dirichlet (Array<LinOpBCType,AMREX_SPACEDIM> bcs) const
 void ERF::project_velocities_tb (int lev, Real l_dt, Vector<MultiFab>& vmf, MultiFab& pmf)
 {
     BL_PROFILE("ERF::project_velocities_tb()");
-    AMREX_ALWAYS_ASSERT(solverChoice.terrain_type == TerrainType::None);
+    AMREX_ALWAYS_ASSERT(solverChoice.mesh_type == MeshType::ConstantDz);
 
     // Make sure the solver only sees the levels over which we are solving
     Vector<BoxArray>            ba_tmp;   ba_tmp.push_back(vmf[Vars::cons].boxArray());
diff --git a/Source/LinearSolvers/ERF_SolveWithFFT.cpp b/Source/LinearSolvers/ERF_SolveWithFFT.cpp
index e499a7f07..22228ea40 100644
--- a/Source/LinearSolvers/ERF_SolveWithFFT.cpp
+++ b/Source/LinearSolvers/ERF_SolveWithFFT.cpp
@@ -12,8 +12,6 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array<MultiFab,
 {
     BL_PROFILE("ERF::solve_with_fft()");
 
-    bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None;
-
     auto const dom_lo = lbound(geom[lev].Domain());
     auto const dom_hi = ubound(geom[lev].Domain());
 
@@ -29,10 +27,9 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array<MultiFab,
     // FFT solve
     // ****************************************************************************
     //
-    // No terrain or stretched grids
     // This calls the full 3D FFT solver with bc's set through bc_fft
     //
-    if (!l_use_terrain)
+    if (solverChoice.mesh_type == MeshType::ConstantDz)
     {
         if (mg_verbose > 0) {
             amrex::Print() << "Using the 3D FFT solver..." << std::endl;
@@ -49,7 +46,8 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array<MultiFab,
     // This calls the hybrid 2D FFT solver + tridiagonal in z with lateral bc's set through bc_fft
     // and Neumann at top and bottom z-boundaries
     //
-    } else if (l_use_terrain && SolverChoice::terrain_is_flat)
+    }
+    else if (solverChoice.mesh_type == MeshType::StretchedDz)
     {
         if (mg_verbose > 0) {
             amrex::Print() << "Using the hybrid FFT solver..." << std::endl;
@@ -89,7 +87,7 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array<MultiFab,
 
         Box const& zbx = mfi.nodaltilebox(2);
         Array4<Real> const& fz_arr  = fluxes[2].array(mfi);
-        if (l_use_terrain && SolverChoice::terrain_is_flat) {
+        if (solverChoice.mesh_type == MeshType::StretchedDz) {
             Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
             ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp
index e5ef6b0bc..f724406d7 100644
--- a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp
+++ b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp
@@ -14,7 +14,7 @@ using namespace amrex;
  * @param[in]  S_stage_prim primitive variables (i.e. conserved variables divided by density) at the last stage
  * @param[in]  pi_stage Exner function at the last stage
  * @param[in]  geom   Container for geometric information
- * @param[in]  terrain_type  Are we using terrain-fitted coordinates
+ * @param[in]  mesh_type     Do we have constant dz?
  * @param[in]  gravity       Magnitude of gravity
  * @param[in]  c_p           Coefficient at constant pressure
  * @param[in]  r0            Reference (hydrostatically stratified) density
@@ -30,7 +30,7 @@ void make_fast_coeffs (int /*level*/,
                        const MultiFab& pi_stage,                       // Exner function evaluated at least stage
                        const amrex::Geometry geom,
                        bool l_use_moisture,
-                       TerrainType terrain_type,
+                       MeshType mesh_type,
                        Real gravity, Real c_p,
                        std::unique_ptr<MultiFab>& detJ_cc,
                        const MultiFab* r0, const MultiFab* pi0,
@@ -76,7 +76,7 @@ void make_fast_coeffs (int /*level*/,
         const Array4<const Real> & stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
         const Array4<const Real> & prim       = S_stage_prim.const_array(mfi);
 
-        const Array4<const Real>& detJ        = (terrain_type != TerrainType::None) ? detJ_cc->const_array(mfi) : Array4<const Real>{};
+        const Array4<const Real>& detJ        = (mesh_type != MeshType::ConstantDz) ? detJ_cc->const_array(mfi) : Array4<const Real>{};
 
         const Array4<const Real>& r0_ca       = r0->const_array(mfi);
         const Array4<const Real>& pi0_ca      = pi0->const_array(mfi);
@@ -107,7 +107,7 @@ void make_fast_coeffs (int /*level*/,
         Real halfg = std::abs(0.5 * grav_gpu[2]);
 
         //Note we don't act on the bottom or top boundaries of the domain
-        if (terrain_type != TerrainType::None)
+        if (mesh_type != MeshType::ConstantDz)
         {
             ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp
index 2cb7c9f60..aaee5f3fe 100644
--- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp
+++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp
@@ -110,7 +110,7 @@ void erf_slow_rhs_post (int level, int finest_level,
     const MultiFab* t_mean_mf = nullptr;
     if (most) t_mean_mf = most->get_mac_avg(0,2);
 
-    const bool l_use_terrain      = (solverChoice.terrain_type != TerrainType::None);
+    const bool l_use_terrain      = (solverChoice.mesh_type != MeshType::ConstantDz);
     const bool l_moving_terrain   = (solverChoice.terrain_type == TerrainType::Moving);
     const bool l_reflux = (solverChoice.coupling_type != CouplingType::OneWay);
     if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain);
diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H
index e21fd3c74..a0731419e 100644
--- a/Source/TimeIntegration/ERF_TI_fast_headers.H
+++ b/Source/TimeIntegration/ERF_TI_fast_headers.H
@@ -112,7 +112,7 @@ void make_fast_coeffs (int level,
                        const amrex::MultiFab& pi_stage,
                        const amrex::Geometry geom,
                        const bool use_moisture,
-                       const TerrainType terrain_type,
+                       const MeshType mesh_type,
                        const amrex::Real gravity,
                        const amrex::Real c_p,
                        std::unique_ptr<amrex::MultiFab>& detJ_cc,
diff --git a/Source/TimeIntegration/ERF_TI_substep_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H
index e336631eb..0717fd577 100644
--- a/Source/TimeIntegration/ERF_TI_substep_fun.H
+++ b/Source/TimeIntegration/ERF_TI_substep_fun.H
@@ -93,7 +93,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
             // We have to call this each step since it depends on the substep time now
             // Note we pass in the *old* detJ here
             make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
-                             l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p,
+                             l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p,
                              detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
             if (fast_step == 0) {
@@ -121,12 +121,12 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                                 mapfac_m[level], mapfac_u[level], mapfac_v[level],
                                 fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             }
-        } else if (solverChoice.terrain_type == TerrainType::Static) {
+        } else if (solverChoice.mesh_type != MeshType::ConstantDz) {
             if (fast_step == 0) {
 
                 // If this is the first substep we make the coefficients since they are based only on stage data
                 make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
-                                 l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p,
+                                 l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p,
                                  detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
                 // If this is the first substep we pass in S_old as the previous step's solution
@@ -152,7 +152,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
 
                 // If this is the first substep we make the coefficients since they are based only on stage data
                 make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
-                                 l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p,
+                                 l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p,
                                  detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
                 // If this is the first substep we pass in S_old as the previous step's solution
diff --git a/Source/Utils/ERF_AverageDown.cpp b/Source/Utils/ERF_AverageDown.cpp
index c080d1e1b..1c28e439b 100644
--- a/Source/Utils/ERF_AverageDown.cpp
+++ b/Source/Utils/ERF_AverageDown.cpp
@@ -55,7 +55,7 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT
         const Box& bx = mfi.tilebox();
         const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
         const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
-        if (SolverChoice::terrain_type != TerrainType::None) {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz) {
             const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
             ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
             {
@@ -119,7 +119,7 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT
         const Box& bx = mfi.tilebox();
         const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
         const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
-        if (SolverChoice::terrain_type != TerrainType::None) {
+        if (SolverChoice::mesh_type != MeshType::ConstantDz) {
             const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
             ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
             {