diff --git a/Build/cmake_with_eb.sh b/Build/cmake_with_eb.sh
deleted file mode 100755
index 386491f48..000000000
--- a/Build/cmake_with_eb.sh
+++ /dev/null
@@ -1,19 +0,0 @@
-#!/bin/bash
-
-# Example CMake config script for an OSX laptop with OpenMPI
-
-cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \
-      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
-      -DCMAKE_C_COMPILER:STRING=mpicc \
-      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
-      -DMPIEXEC_PREFLAGS:STRING=--oversubscribe \
-      -DCMAKE_BUILD_TYPE:STRING=Release \
-      -DAMREX_EB:BOOL=ON \
-      -DERF_DIM:STRING=3 \
-      -DERF_ENABLE_MPI:BOOL=ON \
-      -DERF_ENABLE_TESTS:BOOL=ON \
-      -DERF_ENABLE_FCOMPARE:BOOL=ON \
-      -DERF_ENABLE_DOCUMENTATION:BOOL=OFF \
-      -DERF_ENABLE_EB:BOOL=ON \
-      -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \
-      .. && make -j8
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
index 113125f33..aa6c6ab12 100644
--- a/CMake/BuildERFExe.cmake
+++ b/CMake/BuildERFExe.cmake
@@ -35,19 +35,6 @@ function(build_erf_lib erf_lib_name)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_PARTICLES)
   endif()
 
-  if(ERF_ENABLE_EB)
-    target_sources(${erf_lib_name} PRIVATE
-                   ${SRC_DIR}/EB/ERF_InitEB.cpp
-                   ${SRC_DIR}/EB/ERF_EBBox.cpp
-                   ${SRC_DIR}/EB/ERF_EBCylinder.cpp
-                   ${SRC_DIR}/EB/ERF_EBRegular.cpp
-                   ${SRC_DIR}/EB/ERF_InitEB.cpp
-                   ${SRC_DIR}/EB/ERF_WriteEBSurface.cpp 
-                   ${SRC_DIR}/LinearSolvers/ERF_SolveWithEBMLMG.cpp)
-    target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/EB>)
-    target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_EB)
-  endif()
-
   if(ERF_ENABLE_FFT)
     target_sources(${erf_lib_name} PRIVATE
                    ${SRC_DIR}/LinearSolvers/ERF_SolveWithFFT.cpp)
@@ -138,6 +125,11 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Diffusion/ERF_ComputeStrain_N.cpp
        ${SRC_DIR}/Diffusion/ERF_ComputeStrain_T.cpp
        ${SRC_DIR}/Diffusion/ERF_ComputeTurbulentViscosity.cpp
+       ${SRC_DIR}/EB/ERF_InitEB.cpp
+       ${SRC_DIR}/EB/ERF_EBBox.cpp
+       ${SRC_DIR}/EB/ERF_EBRegular.cpp
+       ${SRC_DIR}/EB/ERF_InitEB.cpp
+       ${SRC_DIR}/EB/ERF_WriteEBSurface.cpp 
        ${SRC_DIR}/Initialization/ERF_InitBCs.cpp
        ${SRC_DIR}/Initialization/ERF_InitCustom.cpp
        ${SRC_DIR}/Initialization/ERF_InitFromHSE.cpp
@@ -192,6 +184,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve_tb.cpp
        ${SRC_DIR}/LinearSolvers/ERF_PoissonWallDist.cpp
        ${SRC_DIR}/LinearSolvers/ERF_ComputeDivergence.cpp 
+       ${SRC_DIR}/LinearSolvers/ERF_SolveWithEBMLMG.cpp
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithGMRES.cpp
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithMLMG.cpp
        ${SRC_DIR}/LinearSolvers/ERF_TerrainPoisson.cpp
@@ -247,6 +240,7 @@ endif()
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/BoundaryConditions>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/DataStructs>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Diffusion>)
+  target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/EB>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Initialization>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/IO>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/LinearSolvers>)
diff --git a/CMake/SetAmrexOptions.cmake b/CMake/SetAmrexOptions.cmake
index ad1c504dc..4a62fd03f 100644
--- a/CMake/SetAmrexOptions.cmake
+++ b/CMake/SetAmrexOptions.cmake
@@ -26,6 +26,7 @@ set(AMReX_CUDA ${ERF_ENABLE_CUDA})
 set(AMReX_ACC OFF)
 set(AMReX_PLOTFILE_TOOLS ${ERF_ENABLE_FCOMPARE})
 set(AMReX_FORTRAN OFF)
+set(AMReX_EB ON)
 
 set(AMReX_LINEAR_SOLVERS ON)
 set(AMReX_LINEAR_SOLVERS_EM OFF)
@@ -41,11 +42,6 @@ if(ERF_ENABLE_FFT)
   set(AMReX_FFT ON)
 endif()
 
-set(AMReX_EB OFF)
-if(ERF_ENABLE_EB)
-  set(AMReX_EB ON)
-endif()
-
 if(ERF_ENABLE_CUDA)
   set(AMReX_GPU_BACKEND CUDA CACHE STRING "AMReX GPU type" FORCE)
   set(AMReX_CUDA_WARN_CAPTURE_THIS OFF)
diff --git a/Exec/DevTests/EB_Test/ERF_Prob.cpp b/Exec/DevTests/EB_Test/ERF_Prob.cpp
index e15bd4222..91c78998b 100644
--- a/Exec/DevTests/EB_Test/ERF_Prob.cpp
+++ b/Exec/DevTests/EB_Test/ERF_Prob.cpp
@@ -169,14 +169,17 @@ Problem::init_custom_pert(
     dxInv[2] = 1. / dx[2];
 
     // Set the z-velocity from impenetrable condition
-    ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-    {
-#ifdef ERF_USE_EB
-        z_vel_pert(i, j, k) = 0.0;
-#else
-        z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv);
-#endif
-    });
+    if (sc.terrain_type == TerrainType::EB) {
+        ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            z_vel_pert(i, j, k) = 0.0;
+        });
+    } else {
+        ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv);
+        });
+    }
 
     amrex::Gpu::streamSynchronize();
 }
diff --git a/Exec/DevTests/EB_Test/GNUmakefile b/Exec/DevTests/EB_Test/GNUmakefile
index 0ad0b3ebf..e9b9c8d3c 100644
--- a/Exec/DevTests/EB_Test/GNUmakefile
+++ b/Exec/DevTests/EB_Test/GNUmakefile
@@ -24,8 +24,6 @@ DEBUG = FALSE
 TEST = TRUE
 USE_ASSERTION = TRUE
 
-USE_EB = TRUE
-
 # GNU Make
 Bpack := ./Make.package
 Blocs := .
diff --git a/Exec/Make.ERF b/Exec/Make.ERF
index 32a208695..f9ce7d172 100644
--- a/Exec/Make.ERF
+++ b/Exec/Make.ERF
@@ -2,66 +2,59 @@ AMREX_HOME      ?= $(ERF_HOME)/Submodules/AMReX
 
 BL_NO_FORT = TRUE
 
+USE_EB = TRUE
+
 include $(AMREX_HOME)/Tools/GNUMake/Make.defs
 
 EBASE = ERF
 
 ERF_SOURCE_DIR = $(ERF_HOME)/Source
 include $(ERF_SOURCE_DIR)/Make.package
-
 VPATH_LOCATIONS   += $(ERF_SOURCE_DIR)
 INCLUDE_LOCATIONS += $(ERF_SOURCE_DIR)
 
 ERF_BC_DIR = $(ERF_SOURCE_DIR)/BoundaryConditions
 include $(ERF_BC_DIR)/Make.package
-
 VPATH_LOCATIONS   += $(ERF_BC_DIR)
 INCLUDE_LOCATIONS += $(ERF_BC_DIR)
 
 ERF_ADVECTION_DIR = $(ERF_SOURCE_DIR)/Advection
 include $(ERF_ADVECTION_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_ADVECTION_DIR)
+INCLUDE_LOCATIONS += $(ERF_ADVECTION_DIR)
 
 ERF_DIFFUSION_DIR = $(ERF_SOURCE_DIR)/Diffusion
 include $(ERF_DIFFUSION_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_DIFFUSION_DIR)
+INCLUDE_LOCATIONS += $(ERF_DIFFUSION_DIR)
 
 ERF_PBL_DIR = $(ERF_SOURCE_DIR)/PBL
 include $(ERF_PBL_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_PBL_DIR)
+INCLUDE_LOCATIONS += $(ERF_PBL_DIR)
 
 ERF_INIT_DIR = $(ERF_SOURCE_DIR)/Initialization
 include $(ERF_INIT_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_INIT_DIR)
+INCLUDE_LOCATIONS += $(ERF_INIT_DIR)
 
 ERF_DATA_DIR = $(ERF_SOURCE_DIR)/DataStructs
 include $(ERF_DATA_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_DATA_DIR)
+INCLUDE_LOCATIONS += $(ERF_DATA_DIR)
 
 ERF_SOLVERS_DIR = $(ERF_SOURCE_DIR)/LinearSolvers
 include $(ERF_SOLVERS_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_SOLVERS_DIR)
+INCLUDE_LOCATIONS += $(ERF_SOLVERS_DIR)
 
 ERF_UTIL_DIR = $(ERF_SOURCE_DIR)/Utils
 include $(ERF_UTIL_DIR)/Make.package
-
-ERF_MULTIBLOCK_DIR = $(ERF_SOURCE_DIR)/MultiBlock
-include $(ERF_MULTIBLOCK_DIR)/Make.package
-
-VPATH_LOCATIONS   += $(ERF_ADVECTION_DIR)
-INCLUDE_LOCATIONS += $(ERF_ADVECTION_DIR)
-
-VPATH_LOCATIONS   += $(ERF_DIFFUSION_DIR)
-INCLUDE_LOCATIONS += $(ERF_DIFFUSION_DIR)
-
-VPATH_LOCATIONS   += $(ERF_PBL_DIR)
-INCLUDE_LOCATIONS += $(ERF_PBL_DIR)
-
-VPATH_LOCATIONS   += $(ERF_INIT_DIR)
-INCLUDE_LOCATIONS += $(ERF_INIT_DIR)
-
-VPATH_LOCATIONS   += $(ERF_DATA_DIR)
-INCLUDE_LOCATIONS += $(ERF_DATA_DIR)
-
 VPATH_LOCATIONS   += $(ERF_UTIL_DIR)
 INCLUDE_LOCATIONS += $(ERF_UTIL_DIR)
 
-VPATH_LOCATIONS   += $(ERF_SOLVERS_DIR)
-INCLUDE_LOCATIONS += $(ERF_SOLVERS_DIR)
+ERF_MULTIBLOCK_DIR = $(ERF_SOURCE_DIR)/MultiBlock
+include $(ERF_MULTIBLOCK_DIR)/Make.package
 
 ifeq ($(USE_MULTIBLOCK),TRUE)
 VPATH_LOCATIONS   += $(ERF_MULTIBLOCK_DIR)
@@ -75,60 +68,30 @@ VPATH_LOCATIONS   += $(ERF_PARTICLES_DIR)
 INCLUDE_LOCATIONS += $(ERF_PARTICLES_DIR)
 endif
 
-ifeq ($(USE_EB),TRUE)
+include $(ERF_PROBLEM_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_PROBLEM_DIR)
+INCLUDE_LOCATIONS += $(ERF_PROBLEM_DIR)
+
 ERF_EB_DIR = $(ERF_SOURCE_DIR)/EB
 include $(ERF_EB_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_EB_DIR)
 INCLUDE_LOCATIONS += $(ERF_EB_DIR)
-endif
 
 ERF_SRCTERMS_DIR = $(ERF_SOURCE_DIR)/SourceTerms
 include $(ERF_SRCTERMS_DIR)/Make.package
-
 VPATH_LOCATIONS   += $(ERF_SRCTERMS_DIR)
 INCLUDE_LOCATIONS += $(ERF_SRCTERMS_DIR)
 
 ERF_TIMEINT_DIR = $(ERF_SOURCE_DIR)/TimeIntegration
 include $(ERF_TIMEINT_DIR)/Make.package
-
 VPATH_LOCATIONS   += $(ERF_TIMEINT_DIR)
 INCLUDE_LOCATIONS += $(ERF_TIMEINT_DIR)
 
 ERF_IO_DIR = $(ERF_SOURCE_DIR)/IO
 include $(ERF_IO_DIR)/Make.package
-
 VPATH_LOCATIONS   += $(ERF_IO_DIR)
 INCLUDE_LOCATIONS += $(ERF_IO_DIR)
 
-include $(ERF_PROBLEM_DIR)/Make.package
-
-VPATH_LOCATIONS   += $(ERF_PROBLEM_DIR)
-INCLUDE_LOCATIONS += $(ERF_PROBLEM_DIR)
-
-include $(AMREX_HOME)/Src/Base/Make.package
-
-AMReXdirs             := Base Boundary AmrCore
-
-ifeq ($(USE_PARTICLES),TRUE)
-AMReXdirs             += Particle
-endif
-
-ifeq ($(USE_EB),TRUE)
-AMReXdirs             += EB
-endif
-
-USE_LINEAR_SOLVERS_INCFLO = TRUE
-USE_LINEAR_SOLVERS_EM     = FALSE
-AMReXdirs                += LinearSolvers
-
-ifeq ($(USE_FFT),TRUE)
-AMReXdirs                += FFT
-endif
-
-AMReXpack             += $(foreach dir, $(AMReXdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
-
-include $(AMReXpack)
-
 ERF_MOISTURE_DIR = $(ERF_SOURCE_DIR)/Microphysics
 include $(ERF_MOISTURE_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_MOISTURE_DIR)
@@ -238,6 +201,30 @@ ERF_LSM_MM5_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel/MM5
 include $(ERF_LSM_MM5_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_LSM_MM5_DIR)
 INCLUDE_LOCATIONS += $(ERF_LSM_MM5_DIR)
+
+include $(AMREX_HOME)/Src/Base/Make.package
+
+AMReXdirs             := Base Boundary AmrCore
+
+ifeq ($(USE_PARTICLES),TRUE)
+AMReXdirs             += Particle
+endif
+
+ifeq ($(USE_EB),TRUE)
+AMReXdirs             += EB
+endif
+
+USE_LINEAR_SOLVERS_INCFLO = TRUE
+USE_LINEAR_SOLVERS_EM     = FALSE
+AMReXdirs                += LinearSolvers
+
+ifeq ($(USE_FFT),TRUE)
+AMReXdirs                += FFT
+endif
+
+AMReXpack             += $(foreach dir, $(AMReXdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
+
+include $(AMReXpack)
         
 ifeq ($(COMPUTE_ERROR), TRUE)
   DEFINES += -DERF_COMPUTE_ERROR
@@ -247,10 +234,6 @@ ifeq ($(USE_PARTICLES), TRUE)
   DEFINES += -DERF_USE_PARTICLES
 endif
 
-ifeq ($(USE_EB), TRUE)
-  DEFINES += -DERF_USE_EB
-endif
-
 ifeq ($(USE_MULTIBLOCK), TRUE)
   DEFINES += -DERF_USE_MULTIBLOCK
 endif
diff --git a/Source/Advection/ERF_Advection.H b/Source/Advection/ERF_Advection.H
index d780b7238..311747476 100644
--- a/Source/Advection/ERF_Advection.H
+++ b/Source/Advection/ERF_Advection.H
@@ -5,6 +5,7 @@
 #include <AMReX_MultiFab.H>
 #include <AMReX_BCRec.H>
 #include <AMReX_YAFluxRegister.H>
+
 #include <ERF_DataStruct.H>
 #include <ERF_IndexDefines.H>
 #include <ERF_ABLMost.H>
@@ -74,7 +75,8 @@ void AdvectionSrcForMom (const amrex::Box& bx,
                          const amrex::Array4<const amrex::Real>& mf_v,
                          const AdvType horiz_adv_type, const AdvType vert_adv_type,
                          const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac,
-                         const bool use_terrain, const int lo_z_face, const int hi_z_face,
+                         TerrainType& terrain_type,
+                         const int lo_z_face, const int hi_z_face,
                          const amrex::Box& domain,
                          const amrex::BCRec* bc_ptr_h);
 
diff --git a/Source/Advection/ERF_AdvectionSrcForMom.cpp b/Source/Advection/ERF_AdvectionSrcForMom.cpp
index ce4cf326e..e9929a4a1 100644
--- a/Source/Advection/ERF_AdvectionSrcForMom.cpp
+++ b/Source/Advection/ERF_AdvectionSrcForMom.cpp
@@ -28,14 +28,13 @@ using namespace amrex;
  * @param[in] ax   Area fraction of x-faces
  * @param[in] ay   Area fraction of y-faces
  * @param[in] az   Area fraction of z-faces
- * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
+ * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain_fitted_coords is false)
  * @param[in] cellSizeInv inverse of the mesh spacing
  * @param[in] mf_m map factor at cell centers
  * @param[in] mf_u map factor at x-faces
  * @param[in] mf_v map factor at y-faces
  * @param[in] horiz_adv_type sets the spatial order to be used for lateral derivatives
  * @param[in] vert_adv_type  sets the spatial order to be used for vertical derivatives
- * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms)
  */
 void
 AdvectionSrcForMom (const Box& bx,
@@ -63,7 +62,7 @@ AdvectionSrcForMom (const Box& bx,
                     const AdvType vert_adv_type,
                     const Real horiz_upw_frac,
                     const Real vert_upw_frac,
-                    const bool use_terrain,
+                    TerrainType& terrain_type,
                     const int lo_z_face, const int hi_z_face,
                     const Box& domain,
                     const BCRec* bc_ptr_h)
@@ -82,6 +81,9 @@ AdvectionSrcForMom (const Box& bx,
     const Array4<Real>& mf_u_inv = mf_u_invFAB.array();
     const Array4<Real>& mf_v_inv = mf_v_invFAB.array();
 
+    const bool use_terrain_fitted_coords = ( terrain_type == TerrainType::StaticFittedMesh ||
+                                             terrain_type == TerrainType::MovingFittedMesh);
+
     ParallelFor(box2d_u, box2d_v,
     [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
     {
@@ -92,23 +94,25 @@ AdvectionSrcForMom (const Box& bx,
         mf_v_inv(i,j,0) = 1. / mf_v(i,j,0);
     });
 
-#ifdef ERF_USE_EB
-    amrex::ignore_unused(use_terrain);
-    ParallelFor(bxx, bxy, bxz,
-    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-    {
-        rho_u_rhs(i, j, k) = 0.0;
-    },
-    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-    {
-        rho_v_rhs(i, j, k) = 0.0;
-    },
-    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    if ( terrain_type == TerrainType::EB ||
+         terrain_type == TerrainType::ImmersedForcing)
     {
-        rho_w_rhs(i, j, k) = 0.0;
-    });
-#else
-    if (!use_terrain) {
+        amrex::ignore_unused(use_terrain_fitted_coords);
+        ParallelFor(bxx, bxy, bxz,
+        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            rho_u_rhs(i, j, k) = 0.0;
+        },
+        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            rho_v_rhs(i, j, k) = 0.0;
+        },
+        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            rho_w_rhs(i, j, k) = 0.0;
+        });
+    } else {
+        if ( !use_terrain_fitted_coords) {
         // Inline with 2nd order for efficiency
         if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd)
         {
@@ -215,9 +219,9 @@ AdvectionSrcForMom (const Box& bx,
                 AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
             }
         }
-    } // end of use_terrain == false
+    } // end of use_terrain_fitted_coords == false
     else
-    { // now do use_terrain = true
+    { // now do use_terrain_fitted_coords = true
         // Inline with 2nd order for efficiency
         if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd)
         {
@@ -354,7 +358,8 @@ AdvectionSrcForMom (const Box& bx,
                     AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
             }
         } // higher order
-    } // terrain
+        } // terrain
+    } // Not EB
 
     // Open bc will be imposed upon all vars (we only access cons here for simplicity)
     const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open);
@@ -440,6 +445,5 @@ AdvectionSrcForMom (const Box& bx,
                                            ax, ay, az, detJ, cellSizeInv,
                                            domhi_z);
     }
-#endif
 }
 
diff --git a/Source/BoundaryConditions/ERF_BoundaryConditionsZvel.cpp b/Source/BoundaryConditions/ERF_BoundaryConditionsZvel.cpp
index ca16c6a35..39dc79b04 100644
--- a/Source/BoundaryConditions/ERF_BoundaryConditionsZvel.cpp
+++ b/Source/BoundaryConditions/ERF_BoundaryConditionsZvel.cpp
@@ -156,7 +156,7 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4<Real      >& dest_a
  * @param[in] bccomp_u  index into m_domain_bcs_type corresponding to u
  * @param[in] bccomp_v  index into m_domain_bcs_type corresponding to v
  * @param[in] bccomp_w  index into m_domain_bcs_type corresponding to w
- * @param[in] terrain_type if Moving then the terrain is moving; otherwise fixed
+ * @param[in] terrain_type if MovingFittedMesh then the terrain is moving; otherwise fixed
  */
 
 void ERFPhysBCFunct_w::impose_vertical_zvel_bcs (const Array4<Real>& dest_arr,
@@ -193,7 +193,7 @@ void ERFPhysBCFunct_w::impose_vertical_zvel_bcs (const Array4<Real>& dest_arr,
     const BCRec* bc_ptr_w_h = bcrs_w.data();
 
     bool l_use_terrain = (m_z_phys_nd != nullptr);
-    bool l_moving_terrain = (terrain_type == TerrainType::Moving);
+    bool l_moving_terrain = (terrain_type == TerrainType::MovingFittedMesh);
 
     GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>,1> l_bc_extdir_vals_d;
 
diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
index 13270cdd0..12ddb0195 100644
--- a/Source/DataStructs/ERF_DataStruct.H
+++ b/Source/DataStructs/ERF_DataStruct.H
@@ -34,7 +34,7 @@ AMREX_ENUM(MeshType,
 );
 
 AMREX_ENUM(TerrainType,
-    None, Static, Moving
+    None, StaticFittedMesh, MovingFittedMesh, EB, ImmersedForcing
 );
 
 AMREX_ENUM(MoistureModelType,
@@ -94,9 +94,9 @@ struct SolverChoice {
             if (mesh_type == MeshType::ConstantDz) {
                 mesh_type    = MeshType::StretchedDz;
             }
-            if (terrain_type != TerrainType::Static) {
+            if (terrain_type != TerrainType::StaticFittedMesh) {
                 amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl;
-                terrain_type = TerrainType::Static;
+                terrain_type = TerrainType::StaticFittedMesh;
             }
             pp.query("zsurface", zsurf);
             if (zsurf != 0.0) {
@@ -158,7 +158,17 @@ struct SolverChoice {
         pp.query_enum_case_insensitive("land_surface_model",lsm_type);
 
         // Is the terrain none, static or moving?
-        pp.query_enum_case_insensitive("terrain_type",terrain_type);
+        std::string terrain_type_temp = "";
+        pp.query("terrain_type", terrain_type_temp);
+        if (terrain_type_temp == "Moving") {
+            amrex::Warning("erf.terrain_type = Moving is deprecated; please replace Moving by MovingFittedMesh");
+            terrain_type = TerrainType::MovingFittedMesh;
+        } else if (terrain_type_temp == "Static") {
+            amrex::Warning("erf.terrain_type = Static is deprecated; please replace Static by StaticFittedMesh");
+            terrain_type = TerrainType::StaticFittedMesh;
+        } else {
+            pp.query_enum_case_insensitive("terrain_type",terrain_type);
+        }
 
         if (terrain_type != TerrainType::None) {
             mesh_type = MeshType::VariableDz;
@@ -168,7 +178,7 @@ struct SolverChoice {
         if (n_zlevels > 0)
         {
             if (terrain_type == TerrainType::None) {
-                terrain_type = TerrainType::Static;
+                terrain_type = TerrainType::StaticFittedMesh;
             }
             if (mesh_type == MeshType::ConstantDz) {
                 mesh_type    = MeshType::StretchedDz;
@@ -403,7 +413,7 @@ struct SolverChoice {
             amrex::Abort("buoyancy_type must be 1, 2, 3 or 4");
         }
 
-        if (!use_lagged_delta_rt && !(terrain_type == TerrainType::Moving)) {
+        if (!use_lagged_delta_rt && !(terrain_type == TerrainType::MovingFittedMesh)) {
             amrex::Error("Can't turn off lagged_delta_rt when terrain not moving");
         }
 
@@ -448,11 +458,10 @@ struct SolverChoice {
         amrex::Print() << "use_coriolis                : " << use_coriolis << std::endl;
         amrex::Print() << "use_gravity                 : " << use_gravity << std::endl;
 
-        amrex::Print() << "Terrain Type: " << std::endl;
-        if (terrain_type == TerrainType::Static) {
-            amrex::Print() << "    Static" << std::endl;
-        } else if (terrain_type == TerrainType::Moving) {
-            amrex::Print() << "    Moving" << std::endl;
+        if (terrain_type == TerrainType::StaticFittedMesh) {
+            amrex::Print() << "Terrain Type: StaticFittedMesh" << std::endl;
+        } else if (terrain_type == TerrainType::MovingFittedMesh) {
+            amrex::Print() << "Terrain Type: MovingFittedMesh" << std::endl;
         } else {
             amrex::Print() << "    None" << std::endl;
         }
@@ -695,8 +704,5 @@ struct SolverChoice {
 
     // Use forest canopy model?
     bool do_forest_drag {false};
-
-    // Use immersed forcing representation of terrain?
-    bool do_terrain_drag {false};
 };
 #endif
diff --git a/Source/EB/ERF_EBCylinder.cpp b/Source/EB/ERF_EBCylinder.cpp
deleted file mode 100644
index 4811e7eb9..000000000
--- a/Source/EB/ERF_EBCylinder.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-#include <AMReX_EB2.H>
-#include <AMReX_EB2_IF.H>
-#include <AMReX_ParmParse.H>
-
-#include <algorithm>
-#include <ERF.H>
-
-using namespace amrex;
-
-/********************************************************************************
- *                                                                              *
- * Function to create a simple cylinder EB.                                     *
- *                                                                              *
- ********************************************************************************/
-void ERF::make_eb_cylinder()
-{
-    // Initialise cylinder parameters
-    bool inside = true;
-    Real radius = 0.0002;
-    int direction = 0;
-    Vector<Real> centervec(3);
-    Real rotation  = 0;
-    int rotation_axe  = 0;
-
-    // Get cylinder information from inputs file.                               *
-    ParmParse pp("cylinder");
-
-    pp.query("internal_flow", inside);
-    pp.query("radius", radius);
-    pp.query("direction", direction);
-    pp.query("rotation",   rotation);
-    pp.query("rotation_axe",   rotation_axe);
-    pp.getarr("center", centervec, 0, 3);
-    Array<Real, AMREX_SPACEDIM> center = {AMREX_D_DECL(centervec[0], centervec[1], centervec[2])};
-
-    rotation = (rotation/180.)*M_PI;
-
-    // Print info about cylinder
-    amrex::Print() << " " << std::endl;
-    amrex::Print() << " Internal Flow: " << inside << std::endl;
-    amrex::Print() << " Radius:    " << radius << std::endl;
-    amrex::Print() << " Direction: " << direction << std::endl;
-    amrex::Print() << " Rotation angle(rad): " << rotation << std::endl;
-    amrex::Print() << " Rotation axe: " << rotation_axe << std::endl;
-#if (AMREX_SPACDEIM == 3)
-    amrex::Print() << " Center:    " << center[0] << ", " << center[1] << ", " << center[2]
-                   << std::endl;
-#else
-    amrex::Print() << " Center:    " << center[0] << ", " << center[1] << std::endl;
-#endif
-
-    // Build the Cylinder implficit function representing the curved walls
-    EB2::CylinderIF my_cyl(radius, direction, center, inside);
-
-    auto my_cyl_rot = EB2::rotate(my_cyl, rotation, rotation_axe);
-
-    // Generate GeometryShop
-    auto gshop = EB2::makeShop(my_cyl_rot);
-
-    // Build index space
-    int max_level_here = 0;
-    int max_coarsening_level = 100;
-    EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level);
-}
diff --git a/Source/EB/ERF_InitEB.cpp b/Source/EB/ERF_InitEB.cpp
index 7b878c72d..532879241 100644
--- a/Source/EB/ERF_InitEB.cpp
+++ b/Source/EB/ERF_InitEB.cpp
@@ -11,8 +11,8 @@ using namespace amrex;
 void ERF::MakeEBGeometry()
 {
    /******************************************************************************
-   * ERF.geometry=<string> specifies the EB geometry. <string> can be one of    *
-   * box, cylinder, or terrain */
+   * ERF.geometry=<string> specifies the EB geometry. <string> can be either     *
+   * box or terrain */
 
     ParmParse pp("eb2");
 
@@ -32,11 +32,7 @@ void ERF::MakeEBGeometry()
         max_coarsening_level = 0;
     }
 
-    if(geom_type == "cylinder") {
-        amrex::Print() << "\n Building cylinder geometry." << std::endl;
-        make_eb_cylinder();
-
-    } else if (geom_type == "terrain") {
+    if (geom_type == "terrain") {
         amrex::Print() << "\n Building EB geometry based on idealized terrain." << std::endl;
         Real dummy_time = 0.0;
         Box bx(surroundingNodes(Geom(0).Domain())); bx.grow(2);
diff --git a/Source/EB/ERF_Redistribute.cpp b/Source/EB/ERF_Redistribute.cpp
index aa98941c1..0fd9b25eb 100644
--- a/Source/EB/ERF_Redistribute.cpp
+++ b/Source/EB/ERF_Redistribute.cpp
@@ -1,7 +1,5 @@
 #include <AMReX_Config.H>
 
-#ifdef ERF_USE_EB
-
 #include <ERF.H>
 #include "AMReX_EB_Redistribution.H"
 
@@ -97,4 +95,3 @@ ERF::redistribute_term ( MFIter const& mfi, int lev,
         });
     }
 }
-#endif
diff --git a/Source/EB/ERF_TerrainIF.H b/Source/EB/ERF_TerrainIF.H
index cbcfbcef6..bd1222be4 100644
--- a/Source/EB/ERF_TerrainIF.H
+++ b/Source/EB/ERF_TerrainIF.H
@@ -19,20 +19,20 @@ public:
           m_geom(a_geom)
         {
             amrex::Print() << " EB type = Terrain " << std::endl;
+
+            dx = m_geom.CellSizeArray()[0];
+            dy = m_geom.CellSizeArray()[1];
+
+            terr_arr = m_terr.const_array();
         }
 
     AMREX_GPU_HOST_DEVICE inline
     amrex::Real operator() (AMREX_D_DECL(amrex::Real x, amrex::Real y, amrex::Real z))
         const noexcept
     {
-        amrex::Real dx = m_geom.CellSizeArray()[0];
-        amrex::Real dy = m_geom.CellSizeArray()[1];
-
         const int i = static_cast<int>(x / dx);
         const int j = static_cast<int>(y / dy);
 
-        amrex::Array4<amrex::Real const> const& terr_arr = m_terr.const_array();
-
         return -(z - terr_arr(i,j,0));
     }
 
@@ -42,8 +42,11 @@ public:
     }
 
 protected:
-    amrex::FArrayBox&  m_terr;
-    amrex::Geometry&   m_geom;
+    amrex::FArrayBox& m_terr;
+    amrex::Array4<amrex::Real const> terr_arr;
+
+    amrex::Geometry&  m_geom;
+    amrex::Real dx, dy;
 };
 
 #endif
diff --git a/Source/EB/Make.package b/Source/EB/Make.package
index 8bae36093..deddc170b 100644
--- a/Source/EB/Make.package
+++ b/Source/EB/Make.package
@@ -2,7 +2,6 @@
 CEXE_sources += main.cpp
 CEXE_sources += ERF_InitEB.cpp
 CEXE_sources += ERF_EBBox.cpp
-CEXE_sources += ERF_EBCylinder.cpp
 CEXE_sources += ERF_EBRegular.cpp
 CEXE_sources += ERF_Redistribute.cpp
 CEXE_sources += ERF_WriteEBSurface.cpp
diff --git a/Source/ERF.H b/Source/ERF.H
index 4a6d87af4..03f55334c 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -21,6 +21,8 @@
 #include <AMReX_PhysBCFunct.H>
 #include <AMReX_YAFluxRegister.H>
 #include <AMReX_ErrorList.H>
+#include <AMReX_EBFabFactory.H>
+#include <AMReX_EBMultiFabUtil.H>
 
 #ifdef ERF_USE_FFT
 #include <AMReX_FFT_Poisson.H>
@@ -161,9 +163,7 @@ public:
     // Initialize multilevel data after MultiBlock
     void InitData_post ();
 
-#ifdef ERF_USE_EB
     void WriteMyEBSurface ();
-#endif
 
     // 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,
@@ -184,10 +184,8 @@ public:
     void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p,
                          amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes);
 #endif
-#ifdef ERF_USE_EB
     void solve_with_EB_mlmg (int lev, amrex::Vector<amrex::MultiFab>& rhs, amrex::Vector<amrex::MultiFab>& p,
                              amrex::Vector<amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>>& fluxes);
-#endif
 
     void solve_with_gmres (int lev, amrex::Vector<amrex::MultiFab>& rhs, amrex::Vector<amrex::MultiFab>& p,
                            amrex::Vector<amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>>& fluxes);
@@ -474,10 +472,8 @@ public:
                        const amrex::Real& time);
 #endif
 
-#ifdef ERF_USE_EB
     void MakeEBGeometry ();
     void make_eb_box ();
-    void make_eb_cylinder ();
     void make_eb_regular ();
     void redistribute_term ( int lev,
                         amrex::MultiFab& result,
@@ -491,7 +487,6 @@ public:
                         amrex::MultiFab const& state,
                         amrex::BCRec const* bc,
                         amrex::Real const dt);
-#endif
 
     // more flexible version of AverageDown() that lets you average down across multiple levels
     void AverageDownTo (int crse_lev, int scomp, int ncomp); // NOLINT
@@ -994,11 +989,9 @@ private:
                                                     "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat",
                                                     "rain_accum", "snow_accum", "graup_accum",
                                                     // Terrain IB mask
-                                                    "terrain_IB_mask"
-#ifdef ERF_USE_EB
+                                                    "terrain_IB_mask",
                                                     // EB variables
-                                                    ,"volfrac",
-#endif
+                                                    "volfrac"
 #ifdef ERF_COMPUTE_ERROR
                                                     // error vars
                                                     ,"xvel_err", "yvel_err", "zvel_err", "pp_err"
@@ -1397,9 +1390,8 @@ private:
     //! The filename of the ith samplelinelog file.
     [[nodiscard]] std::string SampleLineLogName (int i) const noexcept { return samplelinelogname[i]; }
 
-#ifdef ERF_USE_EB
-
-    amrex::Vector<std::unique_ptr<amrex::EBFArrayBoxFactory> > m_factory;
+    // amrex::Vector<std::unique_ptr<amrex::EBFArrayBoxFactory> > m_factory;
+    amrex::Vector<std::unique_ptr<amrex::FabFactory<amrex::FArrayBox> > > m_factory;
 
     [[nodiscard]] amrex::FabFactory<amrex::FArrayBox> const&
     Factory (int lev) const noexcept { return *m_factory[lev]; }
@@ -1417,9 +1409,6 @@ private:
 
     [[nodiscard]] static int nghost_eb_full ()
     { return 4; }
-#else
-    amrex::Vector<std::unique_ptr<amrex::FabFactory<amrex::FArrayBox> > > m_factory;
-#endif
 
 #ifdef ERF_USE_FFT
     amrex::Vector<std::unique_ptr<amrex::FFT::Poisson<amrex::MultiFab>>> m_3D_poisson;
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index f7a5acb4c..6decc4a2c 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -361,10 +361,12 @@ ERF::ERF_shared ()
     // We will create each of these in MakeNewLevel.../RemakeLevel
     m_factory.resize(max_level+1);
 
-#ifdef ERF_USE_EB
     // This is needed before initializing level MultiFabs
-    MakeEBGeometry();
-#endif
+    if ( solverChoice.terrain_type == TerrainType::EB ||
+         solverChoice.terrain_type == TerrainType::ImmersedForcing)
+    {
+        MakeEBGeometry();
+    }
 }
 
 ERF::~ERF () = default;
@@ -587,7 +589,7 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
     }
 
     // Moving terrain
-    if ( solverChoice.terrain_type == TerrainType::Moving )
+    if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
     {
       for (int lev = finest_level; lev >= 0; lev--)
       {
@@ -937,7 +939,7 @@ ERF::InitData_post ()
         base_state[lev].FillBoundary(geom[lev].periodicity());
 
         // For moving terrain only
-        if (solverChoice.terrain_type == TerrainType::Moving) {
+        if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
             MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
             base_state_new[lev].FillBoundary(geom[lev].periodicity());
         }
@@ -1174,11 +1176,13 @@ ERF::InitData_post ()
     pp.query("do_line_sampling",do_line); pp.query("do_plane_sampling",do_plane);
     if (do_line || do_plane) { data_sampler = std::make_unique<SampleData>(do_line, do_plane); }
 
-#ifdef ERF_USE_EB
-    bool write_eb_surface = false;
-    pp.query("write_eb_surface", write_eb_surface);
-    if (write_eb_surface) WriteMyEBSurface();
-#endif
+    if ( solverChoice.terrain_type == TerrainType::EB ||
+         solverChoice.terrain_type == TerrainType::ImmersedForcing)
+    {
+        bool write_eb_surface = false;
+        pp.query("write_eb_surface", write_eb_surface);
+        if (write_eb_surface) WriteMyEBSurface();
+    }
 
 }
 
@@ -1584,24 +1588,6 @@ ERF::ReadParameters ()
         pp.query("cf_width", cf_width);
         pp.query("cf_set_width", cf_set_width);
 
-        // Query the canopy model file name
-        std::string forestfile;
-        solverChoice.do_forest_drag = pp.query("forest_file", forestfile);
-        if (solverChoice.do_forest_drag) {
-            for (int lev = 0; lev <= max_level; ++lev) {
-                m_forest_drag[lev] = std::make_unique<ForestDrag>(forestfile);
-            }
-        }
-
-        //Query the terrain file name
-        std::string terrainfile;
-        solverChoice.do_terrain_drag = pp.query("terrain_file", terrainfile);
-        if (solverChoice.do_terrain_drag) {
-            for (int lev = 0; lev <= max_level; ++lev) {
-                m_terrain_drag[lev] = std::make_unique<TerrainDrag>(terrainfile);
-            }
-        }
-
         // AmrMesh iterate on grids?
         bool iterate(true);
         pp_amr.query("iterate_grids",iterate);
@@ -1618,9 +1604,27 @@ ERF::ReadParameters ()
 
     solverChoice.init_params(max_level);
 
+    // Query the canopy model file name
+    std::string forestfile;
+    solverChoice.do_forest_drag = pp.query("forest_file", forestfile);
+    if (solverChoice.do_forest_drag) {
+        for (int lev = 0; lev <= max_level; ++lev) {
+            m_forest_drag[lev] = std::make_unique<ForestDrag>(forestfile);
+        }
+    }
+
+    // Query the terrain file name (*after* reading in solverChoice inputs)
+    std::string terrainfile;
+    pp.query("terrain_file", terrainfile);
+    if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+        for (int lev = 0; lev <= max_level; ++lev) {
+            m_terrain_drag[lev] = std::make_unique<TerrainDrag>(terrainfile);
+        }
+    }
+
     // No moving terrain with init real (we must do this after init_params
     //    because that is where we set terrain_type
-    if (init_type == InitType::Real && solverChoice.terrain_type == TerrainType::Moving) {
+    if (init_type == InitType::Real && solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
         Abort("Moving terrain is not supported with init real");
     }
 
diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp
index 5c3d91ad8..dc0996042 100644
--- a/Source/ERF_MakeNewArrays.cpp
+++ b/Source/ERF_MakeNewArrays.cpp
@@ -34,7 +34,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     tmp_base_state.define(ba,dm,BaseState::num_comps,3);
     tmp_base_state.setVal(0.);
 
-    if (solverChoice.terrain_type == TerrainType::Moving) {
+    if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
         base_state_new[lev].define(ba,dm,BaseState::num_comps,base_state[lev].nGrowVect());
         base_state_new[lev].setVal(0.);
     }
@@ -46,7 +46,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
         SolverChoice::mesh_type == MeshType::VariableDz) {
         z_phys_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
 
-        if (solverChoice.terrain_type == TerrainType::Moving)
+        if (solverChoice.terrain_type == TerrainType::MovingFittedMesh)
         {
             detJ_cc_new[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
             detJ_cc_src[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
@@ -69,7 +69,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
         int ngrow = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_num_diff) + 2;
         tmp_zphys_nd = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
 
-        if (solverChoice.terrain_type == TerrainType::Moving) {
+        if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
             z_phys_nd_new[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
             z_phys_nd_src[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
         }
@@ -580,5 +580,5 @@ ERF::make_physbcs (int lev)
                                                             solverChoice.terrain_type, z_phys_nd[lev],
                                                             use_real_bcs, zvel_bc_data[lev].data());
     physbcs_base[lev] = std::make_unique<ERFPhysBCFunct_base> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
-                                                               (solverChoice.terrain_type == TerrainType::Moving));
+                                                               (solverChoice.terrain_type == TerrainType::MovingFittedMesh));
 }
diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp
index addc31964..8876df4b8 100644
--- a/Source/ERF_MakeNewLevel.cpp
+++ b/Source/ERF_MakeNewLevel.cpp
@@ -53,15 +53,17 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
 
     if (lev == 0) init_bcs();
 
-#ifdef ERF_USE_EB
-    m_factory[lev] = makeEBFabFactory(geom[lev], grids[lev], dmap[lev],
-                                      {nghost_eb_basic(),
-                                       nghost_eb_volume(),
-                                       nghost_eb_full()},
-                                       EBSupport::full);
-#else
-    m_factory[lev] = std::make_unique<FArrayBoxFactory>();
-#endif
+    if ( solverChoice.terrain_type == TerrainType::EB ||
+         solverChoice.terrain_type == TerrainType::ImmersedForcing)
+    {
+        m_factory[lev] = makeEBFabFactory(geom[lev], grids[lev], dmap[lev],
+                                          {nghost_eb_basic(),
+                                           nghost_eb_volume(),
+                                           nghost_eb_full()},
+                                           EBSupport::full);
+    } else {
+        // m_factory[lev] = std::make_unique<FabFactory<FArrayBox>>();
+    }
 
     auto& lev_new = vars_new[lev];
     auto& lev_old = vars_old[lev];
@@ -138,12 +140,20 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
     #ifdef ERF_USE_WINDFARM
         init_windfarm(lev);
     #endif
+
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    if (solverChoice.do_forest_drag) {
+        m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
 
-    if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    // ********************************************************************************************
+    // Build the data structures for immersed forcing representation of terrain
+    // ********************************************************************************************
+    if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+        m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
 
     //********************************************************************************************
     // Create wall distance field for RANS model (depends upon z_phys)
@@ -229,9 +239,17 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    if (solverChoice.do_forest_drag) {
+        m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
+
+    // ********************************************************************************************
+    // Build the data structures for immersed forcing representation of terrain
+    // ********************************************************************************************
+    if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+        m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
 
-    if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -314,7 +332,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     }
 
     AMREX_ALWAYS_ASSERT(lev > 0);
-    AMREX_ALWAYS_ASSERT(solverChoice.terrain_type != TerrainType::Moving);
+    AMREX_ALWAYS_ASSERT(solverChoice.terrain_type != TerrainType::MovingFittedMesh);
 
     BoxArray            ba_old(vars_new[lev][Vars::cons].boxArray());
     DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap());
@@ -359,9 +377,17 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    if (solverChoice.do_forest_drag) {
+        m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
+
+    // ********************************************************************************************
+    // Build the data structures for immersed forcing representation of terrain
+    // ********************************************************************************************
+    if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+        m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get());
+    }
 
-    if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
     // *****************************************************************************************************
     // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch
     // *****************************************************************************************************
diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp
index 5ce963f46..51c96ccf4 100644
--- a/Source/IO/ERF_Plotfile.cpp
+++ b/Source/IO/ERF_Plotfile.cpp
@@ -1207,12 +1207,16 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
         }
 #endif
 
-#ifdef ERF_USE_EB
         if (containerHasElement(plot_var_names, "volfrac")) {
-            MultiFab::Copy(mf[lev], EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
+            if ( solverChoice.terrain_type == TerrainType::EB ||
+                 solverChoice.terrain_type == TerrainType::ImmersedForcing)
+            {
+                MultiFab::Copy(mf[lev], EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
+            } else {
+                mf[lev].setVal(1.0, mf_comp, 1, 0);
+            }
             mf_comp += 1;
         }
-#endif
 
 #ifdef ERF_COMPUTE_ERROR
         // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some
@@ -1389,11 +1393,12 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 #endif
     }
 
-#ifdef ERF_USE_EB
-    for (int lev = 0; lev <= finest_level; ++lev) {
-        EB_set_covered(mf[lev], 0.0);
+    if (solverChoice.terrain_type == TerrainType::EB)
+    {
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            EB_set_covered(mf[lev], 0.0);
+        }
     }
-#endif
 
     // Fill terrain distortion MF
     if (SolverChoice::mesh_type != MeshType::ConstantDz) {
@@ -1402,7 +1407,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
             Real dz = Geom()[lev].CellSizeArray()[2];
             for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                 const Box& bx = mfi.tilebox();
-                Array4<      Real> mf_arr = mf_nd[lev].array(mfi);
+                Array4<Real> mf_arr = mf_nd[lev].array(mfi);
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                     mf_arr(i,j,k,2) -= k * dz;
                 });
diff --git a/Source/LinearSolvers/ERF_ComputeDivergence.cpp b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
index fbbeb133e..36f9179d6 100644
--- a/Source/LinearSolvers/ERF_ComputeDivergence.cpp
+++ b/Source/LinearSolvers/ERF_ComputeDivergence.cpp
@@ -17,11 +17,12 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
     // Compute divergence which will form RHS
     // Note that we replace "rho0w" with the contravariant momentum, Omega
     // ****************************************************************************
-#ifdef ERF_USE_EB
-    bool already_on_centroids = true;
-    EB_computeDivergence(rhs, rho0_u_const, geom_at_lev, already_on_centroids);
-#else
-    if (SolverChoice::mesh_type == MeshType::ConstantDz)
+    if (solverChoice.terrain_type == TerrainType::EB)
+    {
+        bool already_on_centroids = true;
+        EB_computeDivergence(rhs, rho0_u_const, geom_at_lev, already_on_centroids);
+    }
+    else if (SolverChoice::mesh_type == MeshType::ConstantDz)
     {
         computeDivergence(rhs, rho0_u_const, geom_at_lev);
     }
@@ -63,5 +64,4 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array<MultiFab const*,AMRE
             }
         } // mfi
     }
-#endif
 }
diff --git a/Source/LinearSolvers/ERF_PoissonSolve.cpp b/Source/LinearSolvers/ERF_PoissonSolve.cpp
index 7c9662e9b..df055c8f3 100644
--- a/Source/LinearSolvers/ERF_PoissonSolve.cpp
+++ b/Source/LinearSolvers/ERF_PoissonSolve.cpp
@@ -11,10 +11,8 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
 {
     BL_PROFILE("ERF::project_velocities()");
 
-#ifndef ERF_USE_EB
     auto const dom_lo = lbound(geom[lev].Domain());
     auto const dom_hi = ubound(geom[lev].Domain());
-#endif
 
     // 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());
@@ -26,13 +24,14 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     Vector<MultiFab> rhs;
     Vector<MultiFab> phi;
 
-#ifdef ERF_USE_EB
-    rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
-    phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1, MFInfo(), Factory(lev));
-#else
-    rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0);
-    phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1);
-#endif
+    if (solverChoice.terrain_type == TerrainType::EB)
+    {
+        rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
+        phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1, MFInfo(), Factory(lev));
+    } else {
+        rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0);
+        phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1);
+    }
 
     auto dxInv = geom[lev].InvCellSizeArray();
 
@@ -41,7 +40,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // Now convert the rho0w MultiFab to hold Omega rather than rhow
     // ****************************************************************************
     //
-#ifndef ERF_USE_EB
     if (solverChoice.mesh_type == MeshType::VariableDz)
     {
         for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
@@ -65,7 +63,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
             });
         } // mfi
     }
-#endif
 
     // ****************************************************************************
     // Compute divergence which will form RHS
@@ -106,11 +103,11 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     Vector<Array<MultiFab,AMREX_SPACEDIM> > fluxes;
     fluxes.resize(1);
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-#ifdef ERF_USE_EB
-        fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
-#else
-        fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0);
-#endif
+        if (solverChoice.terrain_type == TerrainType::EB) {
+            fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
+        } else {
+            fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0);
+        }
     }
 
     // ****************************************************************************
@@ -120,9 +117,9 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // ****************************************************************************
     // EB
     // ****************************************************************************
-#ifdef ERF_USE_EB
-    solve_with_EB_mlmg(lev, rhs, phi, fluxes);
-#else
+    if (solverChoice.terrain_type == TerrainType::EB) {
+        solve_with_EB_mlmg(lev, rhs, phi, fluxes);
+    } else {
 
 #ifdef ERF_USE_FFT
         Box my_region(ba_tmp[0].minimalBox());
@@ -189,7 +186,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
 #endif
     } // general terrain
 
-#endif // not EB
+    } // not EB
 
     // ****************************************************************************
     // Print time in solve
diff --git a/Source/LinearSolvers/ERF_PoissonWallDist.cpp b/Source/LinearSolvers/ERF_PoissonWallDist.cpp
index 2555e9204..e94b87e0b 100644
--- a/Source/LinearSolvers/ERF_PoissonWallDist.cpp
+++ b/Source/LinearSolvers/ERF_PoissonWallDist.cpp
@@ -54,12 +54,12 @@ void ERF::poisson_wall_dist (int lev)
     Vector<MultiFab> rhs;
     Vector<MultiFab> phi;
 
-#ifdef ERF_USE_EB
-    Error("Wall dist calc not implemented for EB";
-#else
-    rhs.resize(1);   rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0);
-    phi.resize(1);   phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1);
-#endif
+    if (solverChoice.terrain_type == TerrainType::EB) {
+        amrex::Error("Wall dist calc not implemented for EB");
+    } else {
+        rhs.resize(1);   rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0);
+        phi.resize(1);   phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1);
+    }
 
     rhs[0].setVal(-1.0);
 
@@ -173,7 +173,7 @@ void ERF::poisson_wall_dist (int lev)
     }
 
     // ****************************************************************************
-    // Setup BCs, with solid domain boundaries being dirichlet
+    // Setup BCs, with solid domain boundaries being Dirichlet
     // We assume that the zlo boundary corresponds to the land surface
     // ****************************************************************************
     amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> bc3d_lo, bc3d_hi;
@@ -227,7 +227,7 @@ void ERF::poisson_wall_dist (int lev)
     const Real abstol = solverChoice.poisson_abstol;
 
     Real sigma = 1.0;
-    MLNodeLaplacian mlpoisson(geom_tmp, ba_tmp, dm_tmp, info, {}, sigma);
+    MLNodeLaplacian mlpoisson(geom_tmp, ba_tmp, dm_tmp, info, {m_factory[lev].get()}, sigma);
 
     mlpoisson.setDomainBC(bc3d_lo, bc3d_hi);
 
diff --git a/Source/LinearSolvers/ERF_SolveWithEBMLMG.cpp b/Source/LinearSolvers/ERF_SolveWithEBMLMG.cpp
index 87c7e4fc7..0e5af2aa2 100644
--- a/Source/LinearSolvers/ERF_SolveWithEBMLMG.cpp
+++ b/Source/LinearSolvers/ERF_SolveWithEBMLMG.cpp
@@ -49,7 +49,7 @@ void ERF::solve_with_EB_mlmg (int lev, Vector<MultiFab>& rhs, Vector<MultiFab>&
     // Multigrid solve
     // ****************************************************************************
 
-    MLEBABecLap mleb (geom_tmp, ba_tmp, dm_tmp, info, {m_factory[lev].get()});
+    MLEBABecLap mleb (geom_tmp, ba_tmp, dm_tmp, info, {&EBFactory(lev)});
 
     mleb.setMaxOrder(2);
     mleb.setDomainBC(bclo, bchi);
diff --git a/Source/LinearSolvers/Make.package b/Source/LinearSolvers/Make.package
index a499733e3..28c47468b 100644
--- a/Source/LinearSolvers/Make.package
+++ b/Source/LinearSolvers/Make.package
@@ -6,6 +6,7 @@ CEXE_sources += ERF_ComputeDivergence.cpp
 CEXE_sources += ERF_PoissonSolve.cpp
 CEXE_sources += ERF_PoissonSolve_tb.cpp
 CEXE_sources += ERF_PoissonWallDist.cpp
+CEXE_sources += ERF_SolveWithEBMLMG.cpp
 CEXE_sources += ERF_SolveWithGMRES.cpp
 CEXE_sources += ERF_SolveWithMLMG.cpp
 
@@ -13,7 +14,3 @@ ifeq ($(USE_FFT), TRUE)
 CEXE_headers += ERF_FFTUtils.H
 CEXE_sources += ERF_SolveWithFFT.cpp
 endif
-
-ifeq ($(USE_EB), TRUE)
-CEXE_sources += ERF_SolveWithEBMLMG.cpp
-endif
diff --git a/Source/SourceTerms/ERF_MakeMomSources.cpp b/Source/SourceTerms/ERF_MakeMomSources.cpp
index b0981c485..344ee67ae 100644
--- a/Source/SourceTerms/ERF_MakeMomSources.cpp
+++ b/Source/SourceTerms/ERF_MakeMomSources.cpp
@@ -91,18 +91,19 @@ void make_mom_sources (int level,
     //   10. Immersed Forcing
     // *****************************************************************************
   //const bool l_use_ndiff       = solverChoice.use_num_diff;
-    const bool l_use_zphys       = (solverChoice.mesh_type != MeshType::ConstantDz);
-    const bool l_do_forest_drag  = solverChoice.do_forest_drag;
-    const bool l_do_terrain_drag = solverChoice.do_terrain_drag;
+    const bool l_use_zphys       = (solverChoice.terrain_type == TerrainType::StaticFittedMesh ||
+                                    solverChoice.terrain_type == TerrainType::MovingFittedMesh);
 
-    // Check if terrain and immersed terrain clash
-    if(l_use_zphys && l_do_terrain_drag){
-        amrex::Error(" Cannot use immersed forcing for terrain with terrain-fitted coordinates");
-    }
-    if(l_do_forest_drag && l_do_terrain_drag){
-        amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
+    if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+        if (solverChoice.do_forest_drag) {
+            amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
+        }
+        if (solverChoice.mesh_type != MeshType::ConstantDz) {
+            amrex::Error(" Currently immersed forcing for terrain is only valid for constant dz");
+        }
     }
 
+
     // *****************************************************************************
     // Data for Coriolis forcing
     // *****************************************************************************
@@ -507,7 +508,7 @@ void make_mom_sources (int level,
         // *****************************************************************************
         // 9. Add CANOPY source terms
         // *****************************************************************************
-        if (l_do_forest_drag) {
+        if (solverChoice.do_forest_drag) {
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 const Real ux = u(i, j, k);
@@ -545,7 +546,7 @@ void make_mom_sources (int level,
         // *****************************************************************************
         // 10. Add Immersed source terms
         // *****************************************************************************
-        if (l_do_terrain_drag) {
+        if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
             const Real drag_coefficient=10.0/dz;
             const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H
index 4992e80cd..244b74f5d 100644
--- a/Source/SourceTerms/ERF_SrcHeaders.H
+++ b/Source/SourceTerms/ERF_SrcHeaders.H
@@ -2,14 +2,12 @@
 #define ERF_INTEGRATION_H_
 
 #include <AMReX_MultiFab.H>
+#include <AMReX_EBMultiFabUtil.H>
+
 #include "ERF_DataStruct.H"
 #include "ERF_InputSoundingData.H"
 #include "ERF_TurbPertStruct.H"
 
-#ifdef ERF_USE_EB
-#include <AMReX_EBMultiFabUtil.H>
-#endif
-
 /**
  * Function for computing the buoyancy term to be used in the evolution
  * equation for the z-component of momentum in the slow integrator.  There
diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp
index bd43780be..141596047 100644
--- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp
+++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp
@@ -78,30 +78,80 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const
     int l_implicit_substepping = (solverChoice.substepping_type[level] == SubsteppingType::Implicit);
     int l_anelastic      = solverChoice.anelastic[level];
 
-#ifdef ERF_USE_EB
-    EBFArrayBoxFactory ebfact = EBFactory(level);
-    const MultiFab& detJ = ebfact.getVolFrac();
-#endif
+    Real estdt_comp_inv;
 
-#ifdef ERF_USE_EB
-    Real estdt_comp_inv = ReduceMax(S_new, ccvel, detJ, 0,
-       [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                  Array4<Real const> const& s,
-                                  Array4<Real const> const& u,
-                                  Array4<Real const> const& vf) -> Real
-#else
-    Real estdt_comp_inv = ReduceMax(S_new, ccvel, 0,
+    if (solverChoice.terrain_type == TerrainType::EB)
+    {
+        EBFArrayBoxFactory ebfact = EBFactory(level);
+        const MultiFab& detJ = ebfact.getVolFrac();
+
+        estdt_comp_inv = ReduceMax(S_new, ccvel, detJ, 0,
+        [=] AMREX_GPU_HOST_DEVICE (Box const& b,
+                                   Array4<Real const> const& s,
+                                   Array4<Real const> const& u,
+                                   Array4<Real const> const& vf) -> Real
+        {
+           Real new_comp_dt = -1.e100;
+           amrex::Loop(b, [=,&new_comp_dt] (int i, int j, int k) noexcept
+           {
+               if (vf(i,j,k) > 0.)
+               {
+                   const Real rho      = s(i, j, k, Rho_comp);
+                   const Real rhotheta = s(i, j, k, RhoTheta_comp);
+
+                   // NOTE: even when moisture is present,
+                   //       we only use the partial pressure of the dry air
+                   //       to compute the soundspeed
+                   Real pressure = getPgivenRTh(rhotheta);
+                   Real c = std::sqrt(Gamma * pressure / rho);
+
+                   // If we are doing implicit acoustic substepping, then the z-direction does not contribute
+                   //    to the computation of the time step
+                   if (l_implicit_substepping) {
+                       if ((nxc > 1) && (nyc==1)) {
+                           // 2-D in x-z
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), new_comp_dt);
+                       } else if ((nyc > 1) && (nxc==1)) {
+                           // 2-D in y-z
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt);
+                       } else {
+                           // 3-D or SCM
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
+                                                    ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt);
+                       }
+
+                   // If we are not doing implicit acoustic substepping, then the z-direction contributes
+                   //    to the computation of the time step
+                   } else {
+                       if (nxc > 1 && nyc > 1) {
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
+                                                    ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
+                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
+                       } else if (nxc > 1) {
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
+                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
+                       } else if (nyc > 1) {
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
+                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
+                       } else {
+                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
+                       }
+
+                   }
+               }
+           });
+           return new_comp_dt;
+       });
+
+    } else {
+       estdt_comp_inv = ReduceMax(S_new, ccvel, 0,
        [=] AMREX_GPU_HOST_DEVICE (Box const& b,
                                   Array4<Real const> const& s,
                                   Array4<Real const> const& u) -> Real
-#endif
        {
            Real new_comp_dt = -1.e100;
            amrex::Loop(b, [=,&new_comp_dt] (int i, int j, int k) noexcept
            {
-#ifdef ERF_USE_EB
-               if (vf(i,j,k) > 0.)
-#endif
                {
                    const Real rho      = s(i, j, k, Rho_comp);
                    const Real rhotheta = s(i, j, k, RhoTheta_comp);
@@ -149,6 +199,7 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const
            });
            return new_comp_dt;
        });
+    } // not EB
 
     ParallelDescriptor::ReduceRealMax(estdt_comp_inv);
     estdt_comp = cfl / estdt_comp_inv;
diff --git a/Source/TimeIntegration/ERF_MakeTauTerms.cpp b/Source/TimeIntegration/ERF_MakeTauTerms.cpp
index ecb122eaf..63f16f9ff 100644
--- a/Source/TimeIntegration/ERF_MakeTauTerms.cpp
+++ b/Source/TimeIntegration/ERF_MakeTauTerms.cpp
@@ -37,7 +37,7 @@ void erf_make_tau_terms (int level, int nrk,
     TurbChoice tc = solverChoice.turbChoice[level];
 
     const bool    l_use_terrain    = (solverChoice.terrain_type != TerrainType::None);
-    const bool    l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving);
+    const bool    l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
     if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain);
 
 
diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp
index 2a32c4d65..269740ade 100644
--- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp
+++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp
@@ -80,9 +80,7 @@ void erf_slow_rhs_post (int level, int finest_level,
                         std::unique_ptr<MultiFab>& mapfac_m,
                         std::unique_ptr<MultiFab>& mapfac_u,
                         std::unique_ptr<MultiFab>& mapfac_v,
-#ifdef ERF_USE_EB
                         amrex::EBFArrayBoxFactory const& ebfact,
-#endif
 #if defined(ERF_USE_NETCDF)
                         const bool& moist_set_rhs_bool,
                         const Real& bdy_time_interval,
@@ -111,7 +109,7 @@ void erf_slow_rhs_post (int level, int finest_level,
     if (most) t_mean_mf = most->get_mac_avg(level,2);
 
     const bool l_use_terrain      = (solverChoice.mesh_type != MeshType::ConstantDz);
-    const bool l_moving_terrain   = (solverChoice.terrain_type == TerrainType::Moving);
+    const bool l_moving_terrain   = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
     const bool l_reflux = (solverChoice.coupling_type != CouplingType::OneWay);
     if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain);
 
@@ -320,17 +318,21 @@ void erf_slow_rhs_post (int level, int finest_level,
         // Define updates in the RHS of continuity, temperature, and scalar equations
         // **************************************************************************
         // Metric terms
-#ifdef ERF_USE_EB
-        auto const& ax_arr   = ebfact.getAreaFrac()[0]->const_array(mfi);
-        auto const& ay_arr   = ebfact.getAreaFrac()[1]->const_array(mfi);
-        auto const& az_arr   = ebfact.getAreaFrac()[2]->const_array(mfi);
-        const auto& detJ_arr = ebfact.getVolFrac().const_array(mfi);
-#else
-        auto const& ax_arr   = ax->const_array(mfi);
-        auto const& ay_arr   = ay->const_array(mfi);
-        auto const& az_arr   = az->const_array(mfi);
-        auto const& detJ_arr = detJ->const_array(mfi);
-#endif
+        Array4<const Real> ax_arr;
+        Array4<const Real> ay_arr;
+        Array4<const Real> az_arr;
+        Array4<const Real> detJ_arr;
+        if (solverChoice.terrain_type == TerrainType::EB) {
+            ax_arr   = ebfact.getAreaFrac()[0]->const_array(mfi);
+            ay_arr   = ebfact.getAreaFrac()[1]->const_array(mfi);
+            az_arr   = ebfact.getAreaFrac()[2]->const_array(mfi);
+            detJ_arr = ebfact.getVolFrac().const_array(mfi);
+        } else {
+            ax_arr   = ax->const_array(mfi);
+            ay_arr   = ay->const_array(mfi);
+            az_arr   = az->const_array(mfi);
+            detJ_arr = detJ->const_array(mfi);
+        }
 
         AdvType horiz_adv_type, vert_adv_type;
         Real    horiz_upw_frac, vert_upw_frac;
diff --git a/Source/TimeIntegration/ERF_SlowRhsPre.cpp b/Source/TimeIntegration/ERF_SlowRhsPre.cpp
index 2168315c8..c28cf84d9 100644
--- a/Source/TimeIntegration/ERF_SlowRhsPre.cpp
+++ b/Source/TimeIntegration/ERF_SlowRhsPre.cpp
@@ -109,18 +109,12 @@ void erf_slow_rhs_pre (int level, int finest_level,
                        std::unique_ptr<MultiFab>& mapfac_m,
                        std::unique_ptr<MultiFab>& mapfac_u,
                        std::unique_ptr<MultiFab>& mapfac_v,
-#ifdef ERF_USE_EB
                        EBFArrayBoxFactory const& ebfact,
-#endif
                        YAFluxRegister* fr_as_crse,
                        YAFluxRegister* fr_as_fine)
 {
     BL_PROFILE_REGION("erf_slow_rhs_pre()");
 
-#ifdef ERF_USE_EB
-    amrex::ignore_unused(ax,ay,az,detJ);
-#endif
-
     const BCRec* bc_ptr_d = domain_bcs_type_d.data();
     const BCRec* bc_ptr_h = domain_bcs_type_h.data();
 
@@ -139,7 +133,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
     const Real    l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac;
     const Real    l_vert_upw_frac  = solverChoice.advChoice.dycore_vert_upw_frac;
     const bool    l_use_terrain    = (solverChoice.terrain_type != TerrainType::None);
-    const bool    l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving);
+    const bool    l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
     if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain);
 
     const bool l_use_mono_adv   = solverChoice.use_mono_adv;
@@ -447,17 +441,22 @@ void erf_slow_rhs_pre (int level, int finest_level,
         // *****************************************************************************
         // Define updates in the RHS of continuity and potential temperature equations
         // *****************************************************************************
-#ifdef ERF_USE_EB
-        auto const& ax_arr   = ebfact.getAreaFrac()[0]->const_array(mfi);
-        auto const& ay_arr   = ebfact.getAreaFrac()[1]->const_array(mfi);
-        auto const& az_arr   = ebfact.getAreaFrac()[2]->const_array(mfi);
-        const auto& detJ_arr = ebfact.getVolFrac().const_array(mfi);
-#else
-        auto const& ax_arr   = ax->const_array(mfi);
-        auto const& ay_arr   = ay->const_array(mfi);
-        auto const& az_arr   = az->const_array(mfi);
-        auto const& detJ_arr = detJ->const_array(mfi);
-#endif
+        Array4<const Real> ax_arr;
+        Array4<const Real> ay_arr;
+        Array4<const Real> az_arr;
+        Array4<const Real> detJ_arr;
+        if (solverChoice.terrain_type == TerrainType::EB)
+        {
+            ax_arr   = ebfact.getAreaFrac()[0]->const_array(mfi);
+            ay_arr   = ebfact.getAreaFrac()[1]->const_array(mfi);
+            az_arr   = ebfact.getAreaFrac()[2]->const_array(mfi);
+            detJ_arr = ebfact.getVolFrac().const_array(mfi);
+        } else {
+            ax_arr   = ax->const_array(mfi);
+            ay_arr   = ay->const_array(mfi);
+            az_arr   = az->const_array(mfi);
+            detJ_arr = detJ->const_array(mfi);
+        }
 
         AdvectionSrcForRho(bx, cell_rhs,
                            rho_u, rho_v, omega_arr,      // these are being used to build the fluxes
@@ -561,7 +560,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
                            dxInv, mf_m, mf_u, mf_v,
                            l_horiz_adv_type, l_vert_adv_type,
                            l_horiz_upw_frac, l_vert_upw_frac,
-                           l_use_terrain, lo_z_face, hi_z_face,
+                           solverChoice.terrain_type, lo_z_face, hi_z_face,
                            domain, bc_ptr_h);
 
         if (l_use_diff) {
diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H
index a0731419e..1ffc0fdae 100644
--- a/Source/TimeIntegration/ERF_TI_fast_headers.H
+++ b/Source/TimeIntegration/ERF_TI_fast_headers.H
@@ -4,17 +4,14 @@
 #include <AMReX_MultiFab.H>
 #include <AMReX_BCRec.H>
 #include <AMReX_YAFluxRegister.H>
+#include <AMReX_EBMultiFabUtil.H>
+
 #include "ERF_DataStruct.H"
 #include "ERF_IndexDefines.H"
 #include <ERF_TerrainMetrics.H>
-
 #include <ERF_TileNoZ.H>
 #include <ERF_ProbCommon.H>
 
-#ifdef ERF_USE_EB
-#include <AMReX_EBMultiFabUtil.H>
-#endif
-
 /**
  * Function for computing the fast RHS with no terrain
  *
diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H
index 074c0c2c7..1ecd4a8a2 100644
--- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H
+++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H
@@ -55,7 +55,7 @@
                 Array4<Real>* fslow = fslow_d.dataPtr();
 
                 // Moving terrain
-                if ( solverChoice.terrain_type == TerrainType::Moving )
+                if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
                 {
                     const Array4<const Real>& dJ_old =     detJ_cc[level]->const_array(mfi);
                     const Array4<const Real>& dJ_new = detJ_cc_new[level]->const_array(mfi);
@@ -106,59 +106,57 @@
                     });
 
                 } else { // Fixed or no terrain
-#ifdef ERF_USE_EB
                     const Array4<const Real>& dJ_old =     detJ_cc[level]->const_array(mfi);
-#endif
                     ParallelFor(bx, ncomp_fast[IntVars::cons],
                     [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
                         const int n = scomp_fast[IntVars::cons] + nn;
-#ifdef ERF_USE_EB
                         if (dJ_old(i,j,k) > 0.0) {
-#endif
                             ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
                                ( fslow[IntVars::cons](i,j,k,n) );
-#ifdef ERF_USE_EB
                         }
-#endif
                     });
 
-                    // Commenting out the update is a total HACK while developing the EB capability
-                    ParallelFor(tbx, tby, tbz,
-                    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
-#ifdef ERF_USE_EB
-                        ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
-#else
-                        ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
-                                                  + slow_dt * fslow[IntVars::xmom](i,j,k);
-#endif
-                    },
-                    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
-#ifdef ERF_USE_EB
-                        ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
-#else
-                        ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
-                                                  + slow_dt * fslow[IntVars::ymom](i,j,k);
-#endif
-                    },
-                    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
-#ifdef ERF_USE_EB
-                        ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
-#else
-                        ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
-                                                  + slow_dt * fslow[IntVars::zmom](i,j,k);
-#endif
-                    });
-                }
+                    // Commenting out the update is a HACK while developing the EB capability
+                    if (solverChoice.terrain_type == TerrainType::EB)
+                    {
+                        ParallelFor(tbx, tby, tbz,
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
+                        },
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
+                        },
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
+                        });
+
+                    } else {
+                        ParallelFor(tbx, tby, tbz,
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
+                                                      + slow_dt * fslow[IntVars::xmom](i,j,k);
+                        },
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
+                                                      + slow_dt * fslow[IntVars::ymom](i,j,k);
+                        },
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                            ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
+                                                      + slow_dt * fslow[IntVars::zmom](i,j,k);
+                        });
+                    } // no EB
+                } // not moving terrain
             } // mfi
-#ifdef ERF_USE_EB
-            // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF
-            // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level));
-            // dUdt_tmp.setVal(0.);
-            // dUdt_tmp.FillBoundary(fine_geom.periodicity());
-            // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both?
-            // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp,
-            //                     S_sum[IntVars::cons], bc_ptr_h, slow_dt);
-#endif
+
+            if (solverChoice.terrain_type == TerrainType::EB) {
+                // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF
+                // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level));
+                // dUdt_tmp.setVal(0.);
+                // dUdt_tmp.FillBoundary(fine_geom.periodicity());
+                // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both?
+                // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp,
+                //                     S_sum[IntVars::cons], bc_ptr_h, slow_dt);
+            } // EB
         } // omp
 
         // Even if we update all the conserved variables we don't need
diff --git a/Source/TimeIntegration/ERF_TI_slow_headers.H b/Source/TimeIntegration/ERF_TI_slow_headers.H
index 467b065f5..4f1a7bcd5 100644
--- a/Source/TimeIntegration/ERF_TI_slow_headers.H
+++ b/Source/TimeIntegration/ERF_TI_slow_headers.H
@@ -4,10 +4,12 @@
 #include <AMReX_MultiFab.H>
 #include <AMReX_BCRec.H>
 #include <AMReX_YAFluxRegister.H>
+#include <AMReX_MultiCutFab.H>
+#include <AMReX_EBMultiFabUtil.H>
+
 #include "ERF_DataStruct.H"
 #include "ERF_IndexDefines.H"
 #include "ERF_ABLMost.H"
-
 #include <ERF_Advection.H>
 #include <ERF_Diffusion.H>
 #include <ERF_EddyViscosity.H>
@@ -15,10 +17,6 @@
 #include <ERF_TerrainMetrics.H>
 #include <ERF_TileNoZ.H>
 
-#ifdef ERF_USE_EB
-#include <AMReX_MultiCutFab.H>
-#include <AMReX_EBMultiFabUtil.H>
-#endif
 
 void erf_make_tau_terms (int level, int nrk,
                          const amrex::Vector<amrex::BCRec>& domain_bcs_type,
@@ -100,9 +98,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk,
                       std::unique_ptr<amrex::MultiFab>& mapfac_m,
                       std::unique_ptr<amrex::MultiFab>& mapfac_u,
                       std::unique_ptr<amrex::MultiFab>& mapfac_v,
-#ifdef ERF_USE_EB
                       amrex::EBFArrayBoxFactory const& ebfact,
-#endif
                       amrex::YAFluxRegister* fr_as_crse,
                       amrex::YAFluxRegister* fr_as_fine);
 
@@ -147,9 +143,7 @@ void erf_slow_rhs_post (int level, int finest_level, int nrk,
                        std::unique_ptr<amrex::MultiFab>& mapfac_m,
                        std::unique_ptr<amrex::MultiFab>& mapfac_u,
                        std::unique_ptr<amrex::MultiFab>& mapfac_v,
-#ifdef ERF_USE_EB
                        amrex::EBFArrayBoxFactory const& ebfact,
-#endif
 #if defined(ERF_USE_NETCDF)
                        const bool& moist_zero,
                        const amrex::Real& bdy_time_interval,
diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H
index 5a6b09a6f..8e21174ca 100644
--- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H
+++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H
@@ -65,10 +65,12 @@
         if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
         // Immersed Forcing
         MultiFab* terrain_blank = nullptr;
-        if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
+        if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+            terrain_blank = m_terrain_drag[level]->get_terrain_blank_field();
+        }
 
         // Moving terrain
-        if ( solverChoice.terrain_type == TerrainType::Moving )
+        if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
         {
             // Note that the "old" and "new" metric terms correspond to
             // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source
@@ -151,9 +153,7 @@
                              z_phys_nd_src[level], ax_src[level], ay_src[level], az_src[level], detJ_cc_src[level], p0_new,
                              pp_inc[level],
                              mapfac_m[level], mapfac_u[level], mapfac_v[level],
-#ifdef ERF_USE_EB
                              EBFactory(level),
-#endif
                              fr_as_crse, fr_as_fine);
 
             add_thin_body_sources(xmom_src, ymom_src, zmom_src,
@@ -260,9 +260,7 @@
                              z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0,
                              pp_inc[level],
                              mapfac_m[level], mapfac_u[level], mapfac_v[level],
-#ifdef ERF_USE_EB
                              EBFactory(level),
-#endif
                              fr_as_crse, fr_as_fine);
 
             add_thin_body_sources(xmom_src, ymom_src, zmom_src,
@@ -349,7 +347,7 @@
         }
 
         // Moving terrain
-        if ( solverChoice.terrain_type == TerrainType::Moving ) {
+        if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh ) {
             erf_slow_rhs_post(level, finest_level, nrk, slow_dt, n_qstate,
                               S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
                               xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
@@ -357,9 +355,7 @@
                               fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
                               z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc_new[level],
                               mapfac_m[level], mapfac_u[level], mapfac_v[level],
-#ifdef ERF_USE_EB
                               EBFactory(level),
-#endif
 #if defined(ERF_USE_NETCDF)
                               moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time,
                               real_width, real_set_width,
@@ -374,9 +370,7 @@
                               fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
                               z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc[level],
                               mapfac_m[level], mapfac_u[level], mapfac_v[level],
-#ifdef ERF_USE_EB
                               EBFactory(level),
-#endif
 #if defined(ERF_USE_NETCDF)
                               moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time,
                               real_width, real_set_width,
@@ -438,7 +432,9 @@
         if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
         // Immersed Forcing
         MultiFab* terrain_blank = nullptr;
-        if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
+        if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+            terrain_blank = m_terrain_drag[level]->get_terrain_blank_field();
+        }
 
         make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level],
 #if defined(ERF_USE_RRTMGP)
@@ -473,9 +469,7 @@
                          z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0,
                          pp_inc[level],
                          mapfac_m[level], mapfac_u[level], mapfac_v[level],
-#ifdef ERF_USE_EB
                          EBFactory(level),
-#endif
                          fr_as_crse, fr_as_fine);
 
          add_thin_body_sources(xmom_src, ymom_src, zmom_src,
diff --git a/Source/TimeIntegration/ERF_TI_substep_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H
index e2df0963d..f6695a364 100644
--- a/Source/TimeIntegration/ERF_TI_substep_fun.H
+++ b/Source/TimeIntegration/ERF_TI_substep_fun.H
@@ -48,7 +48,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
 
         // Moving terrain
         std::unique_ptr<MultiFab> z_t_pert;
-        if ( solverChoice.terrain_type == TerrainType::Moving )
+        if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
         {
             // Make "old" fast geom -- store in z_phys_nd for convenience
             if (verbose) Print() << "Making geometry at start of substep time: " << old_substep_time << std::endl;