diff --git a/CHANGES b/CHANGES
index bb63b3a5d..d8429d896 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,22 @@
+# 24.03
+    -- AMReX submodule set to 24.03 release hash (d737886)
+
+    -- Plot rain accumulation (#1478)
+
+    -- Remove FastEddy completely (#1483)
+
+    -- Rain accumulation capability (#1474)
+
+    -- Remove the FastEddy model since it coincides with Kessler (#1473)
+
+    -- Input_sounding with moisture new HSE balance (#1469)
+
+    -- SAM microphysics re-write for ERF dycore (#1461/#1447)
+
+    -- Two-way coupled AMR (#1458)
+
+    -- Numerical diffusion bug fix (#1451)
+
 # 24.02
     -- AMReX submodule set to 24.02 release hash (52393d3)
 
diff --git a/Docs/sphinx_doc/BestPractices.rst b/Docs/sphinx_doc/BestPractices.rst
new file mode 100644
index 000000000..8702a2851
--- /dev/null
+++ b/Docs/sphinx_doc/BestPractices.rst
@@ -0,0 +1,63 @@
+ .. role:: cpp(code)
+    :language: c++
+
+.. _GettingStarted:
+
+Best Practices
+==============
+
+Please note this section is a work in progress and aims to offer guidelines
+but, in general, your mileage may vary.
+
+Large-Eddy Simulations
+----------------------
+
+* Advection Scheme
+
+  - A WRF-like configuration is generally robust, but may under-resolve
+    turbulence features.
+
+    .. code-block:: python
+
+       erf.dycore_horiz_adv_type  = "Upwind_5th"
+       erf.dycore_vert_adv_type   = "Upwind_3rd"
+       erf.dryscal_horiz_adv_type = "Upwind_5th"
+       erf.dryscal_vert_adv_type  = "Upwind_3rd"
+
+  - Centered difference schemes will generally give some non-physical
+    numerical noise, clearly visible in the free atmosphere, but may also
+    better resolve turbulence features. With ``Centered_2nd``, the simulation
+    may remain numerically stable but without any upwinding or numerical
+    diffusion, these results should be carefully interpreted.
+
+  - For higher-order central differencing alone (i.e., without any added
+    upwinding), at least 5% numerical diffusion should be included to stabilize
+    the solution; this was tested with ``Centered_6th``. Note that this does not
+    necessarily kill the numerical noise and is only for numerical stability.
+    These options are identical to WRF's ``diff_6th_opt`` (default: off) and
+    ``diff_6th_factor`` (default: 12%) options.
+
+    .. code-block:: python
+
+       erf.use_NumDiff  = true
+       erf.NumDiffCoeff = 0.05
+
+* Time Integration
+
+  - Split timestepping offers some computational cost savings but still does
+    not allow you to run with an incompressible time-step size.
+  - The acoustic CFL should be less than 0.5, with 4--6 substeps. E.g.,:
+
+    .. code-block:: python
+
+       erf.fixed_dt           = 0.06  # slow timestep
+
+       # These are equivalent and result in a fixed fast timestep size
+       #   if dx=10, speed of sound ~ 350 m/s
+       erf.fixed_fast_dt      = 0.01  # ==> CFL~0.35
+       #erf.fixed_mri_dt_ratio = 6
+
+       # Alternatively, let ERF chose the fast timestep
+       #erf.cfl                = 0.5
+
+
diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
index 8c455531a..acea95839 100644
--- a/Docs/sphinx_doc/Inputs.rst
+++ b/Docs/sphinx_doc/Inputs.rst
@@ -1066,7 +1066,10 @@ Examples of Usage
 Moisture
 ========
 
-ERF has several different moisture models.
+ERF has several different moisture models. The models that are currently implemented
+are Eulerian models; however, ERF has the capability for Lagrangian models when
+compiled with particles.
+
 The following run-time options control how the full moisture model is used.
 
 List of Parameters
diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst
index 1be6e6abe..178d0b7fe 100644
--- a/Docs/sphinx_doc/index.rst
+++ b/Docs/sphinx_doc/index.rst
@@ -39,6 +39,7 @@ In addition to this documentation, there is API documentation for ERF generated
 
    GettingStarted.rst
    Inputs.rst
+   BestPractices.rst
 
 .. toctree::
    :caption: THEORY
diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H
index a87161180..b41ebb0f4 100644
--- a/Source/DataStructs/DataStruct.H
+++ b/Source/DataStructs/DataStruct.H
@@ -28,6 +28,10 @@ enum struct TerrainType {
     Static, Moving
 };
 
+enum struct MoistureModelType{
+    Eulerian, Lagrangian, Undefined
+};
+
 enum struct MoistureType {
     Kessler, SAM, Kessler_NoRain, None
 };
diff --git a/Source/ERF.H b/Source/ERF.H
index 71c1e7b29..6b5808485 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -43,7 +43,8 @@
 #include "ParticleData.H"
 #endif
 
-#include "Microphysics.H"
+#include "EulerianMicrophysics.H"
+#include "LagrangianMicrophysics.H"
 #include "LandSurface.H"
 
 #ifdef ERF_USE_RRTMGP
@@ -374,6 +375,9 @@ private:
 
     void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf);
 
+    // Initialize the microphysics object
+    void initializeMicrophysics (const int&);
+
     // Compute a vector of new MultiFabs by copying from valid region and filling ghost cells
     //
     // NOTE: FillPatch takes in an empty MF, and returns cell-centered + velocities (not momenta)
@@ -552,7 +556,7 @@ private:
     amrex::Vector<amrex::MultiFab> rW_old;
     amrex::Vector<amrex::MultiFab> rW_new;
 
-    Microphysics micro;
+    std::unique_ptr<Microphysics> micro;
     amrex::Vector<amrex::Vector<amrex::MultiFab*>> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg
 
 //#if defined(ERF_USE_WINDFARM)
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index c6dcf39eb..a6672e829 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -111,10 +111,6 @@ ERF::ERF ()
 
     int nlevs_max = max_level + 1;
 
-    // NOTE: size micro before readparams (chooses the model at all levels)
-    micro.ReSize(nlevs_max);
-    qmoist.resize(nlevs_max);
-
 #ifdef ERF_USE_WINDFARM
     if(solverChoice.windfarm_type == WindFarmType::Fitch){
         Nturb.resize(nlevs_max);
@@ -132,6 +128,8 @@ ERF::ERF ()
     lsm_flux.resize(nlevs_max);
 
     ReadParameters();
+    initializeMicrophysics(nlevs_max);
+
     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);
 
@@ -609,11 +607,23 @@ ERF::InitData ()
         if (solverChoice.moisture_type != MoistureType::None)
         {
             for (int lev = 0; lev <= finest_level; lev++) {
-                FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
+                if (qmoist[lev].size() > 0) FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
             }
         }
     }
 
+#ifdef ERF_USE_PARTICLES
+    /* If using a Lagrangian microphysics model, its particle container has now been
+       constructed and initialized (calls to micro->Init). So, add its pointer to
+       ERF::particleData and remove its name from list of unallocated particle containers. */
+    if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+        const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
+        const auto& pc_ptr( dynamic_cast<LagrangianMicrophysics&>(*micro).getParticleContainer() );
+        particleData.pushBack(pc_name, pc_ptr);
+        particleData.getNamesUnalloc().remove(pc_name);
+    }
+#endif
+
     if (input_bndry_planes) {
         // Read the "time.dat" file to know what data is available
         m_r2d->read_time_file();
@@ -801,7 +811,7 @@ ERF::InitData ()
 
     // Update micro vars before first plot file
     if (solverChoice.moisture_type != MoistureType::None) {
-        for (int lev = 0; lev <= finest_level; ++lev) micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
+        for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
     }
 
     // check for additional plotting variables that are available after particle containers
@@ -912,6 +922,34 @@ ERF::InitData ()
     BL_PROFILE_VAR_STOP(InitData);
 }
 
+// Initialize microphysics object
+void
+ERF::initializeMicrophysics (const int& a_nlevsmax /*!< number of AMR levels */)
+{
+    if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Eulerian) {
+
+        micro = std::make_unique<EulerianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+
+    } else if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+#ifdef ERF_USE_PARTICLES
+
+        micro = std::make_unique<LagrangianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+        /* Lagrangian microphysics models will have a particle container; it needs to be added
+           to ERF::particleData */
+        const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
+        /* The particle container has not yet been constructed and initialized, so just add
+           its name here for now (so that functions to set plotting variables can see it). */
+        particleData.addName( pc_name );
+
+#else
+        amrex::Abort("Lagrangian microphysics can be used when compiled with ERF_USE_PARTICLES");
+#endif
+    }
+
+    qmoist.resize(a_nlevsmax);
+    return;
+}
+
 void
 ERF::restart ()
 {
@@ -1227,22 +1265,6 @@ ERF::ReadParameters ()
 
     solverChoice.init_params(max_level);
 
-    // What type of moisture model to use
-    // NOTE: Must be checked after init_params
-    if (solverChoice.moisture_type == MoistureType::SAM) {
-        micro.SetModel<SAM>();
-        amrex::Print() << "SAM moisture model!\n";
-    } else if (solverChoice.moisture_type == MoistureType::Kessler or
-               solverChoice.moisture_type == MoistureType::Kessler_NoRain) {
-        micro.SetModel<Kessler>();
-        amrex::Print() << "Kessler moisture model!\n";
-    } else if (solverChoice.moisture_type == MoistureType::None) {
-        micro.SetModel<NullMoist>();
-        amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n";
-    } else {
-        amrex::Abort("Dont know this moisture_type!") ;
-    }
-
     // What type of land surface model to use
     // NOTE: Must be checked after init_params
     if (solverChoice.lsm_type == LandSurfaceType::SLM) {
@@ -1307,12 +1329,13 @@ ERF::MakeHorizontalAverages ()
             auto  fab_arr = mf.array(mfi);
             auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
             auto const   qv_arr = qmoist[lev][0]->const_array(mfi);
+            int ncomp = vars_new[lev][Vars::cons].nComp();
 
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                 Real dens = cons_arr(i, j, k, Rho_comp);
                 fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
-                fab_arr(i, j, k, 3) = cons_arr(i, j, k, RhoQ1_comp) / dens;
-                fab_arr(i, j, k, 4) = cons_arr(i, j, k, RhoQ2_comp) / dens;
+                fab_arr(i, j, k, 3) = (ncomp > RhoQ1_comp ? cons_arr(i, j, k, RhoQ1_comp) / dens : 0.0);
+                fab_arr(i, j, k, 4) = (ncomp > RhoQ2_comp ? cons_arr(i, j, k, RhoQ2_comp) / dens : 0.0);
             });
         }
 
@@ -1604,10 +1627,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
 
     int nlevs_max = max_level + 1;
 
-    // NOTE: size micro before readparams (chooses the model at all levels)
-    micro.ReSize(nlevs_max);
-    qmoist.resize(nlevs_max);
-
 #ifdef ERF_USE_WINDFARM
     if(solverChoice.windfarm_type == WindFarmType::Fitch){
         Nturb.resize(nlevs_max);
@@ -1625,6 +1644,8 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
     lsm_flux.resize(nlevs_max);
 
     ReadParameters();
+    initializeMicrophysics(nlevs_max);
+
     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);
 
diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp
index 6a559ba6b..9671d52ac 100644
--- a/Source/ERF_make_new_level.cpp
+++ b/Source/ERF_make_new_level.cpp
@@ -44,7 +44,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
 
     init_stuff(lev, ba, dm);
 
-    int n_qstate   = micro.Get_Qstate_Size();
+    int n_qstate   = micro->Get_Qstate_Size();
     int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);
 
     lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state);
@@ -62,15 +62,15 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size  = micro.Get_Qmoist_Size();
+    int q_size  = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     //********************************************************************************************
@@ -225,15 +225,15 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size  = micro.Get_Qmoist_Size();
+    int q_size  = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     init_stuff(lev, ba, dm);
@@ -340,15 +340,15 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size = micro.Get_Qmoist_Size();
+    int q_size = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     init_stuff(lev, ba, dm);
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
index 1a40d414b..8d58f3c3b 100644
--- a/Source/IO/Plotfile.cpp
+++ b/Source/IO/Plotfile.cpp
@@ -44,7 +44,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     // since they may be in any order in the input list
     Vector<std::string> tmp_plot_names;
 
-    int n_qstate   = micro.Get_Qstate_Size();
+    int n_qstate   = micro->Get_Qstate_Size();
     int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);
 
     for (int i = 0; i < ncomp_cons; ++i) {
@@ -164,7 +164,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
     bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
     for (int lev = 0; lev <= finest_level; ++lev) {
         for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-            qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+            qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
         }
     }
 
@@ -260,10 +260,11 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 const Box& bx = mfi.tilebox();
                 const Array4<Real      >& derdat = mf[lev].array(mfi);
                 const Array4<Real const>&  S_arr = vars_new[lev][Vars::cons].const_array(mfi);
+                const int ncomp = vars_new[lev][Vars::cons].nComp();
 
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
-                    Real qv_for_p = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0;
+                    Real qv_for_p = (use_moisture && (ncomp > RhoQ1_comp)) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0;
                     const Real rhotheta = S_arr(i,j,k,RhoTheta_comp);
                     derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p);
                 });
@@ -281,10 +282,11 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 const Array4<Real>& derdat = mf[lev].array(mfi);
                 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
                 const Array4<Real const>& S_arr = vars_new[lev][Vars::cons].const_array(mfi);
+                const int ncomp = vars_new[lev][Vars::cons].nComp();
 
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                 {
-                    Real qv_for_p = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0;
+                    Real qv_for_p = (use_moisture && (ncomp > RhoQ1_comp)) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0;
                     const Real rhotheta = S_arr(i,j,k,RhoTheta_comp);
                     derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k);
                 });
@@ -319,9 +321,10 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 const Box& bx = mfi.tilebox();
                 const Array4<Real>& derdat  = mf[lev].array(mfi);
                 const Array4<Real const>& S_arr = vars_new[lev][Vars::cons].const_array(mfi);
+                const int ncomp = vars_new[lev][Vars::cons].nComp();
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                    Real qv = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0.0;
-                    Real qc = (use_moisture) ? S_arr(i,j,k,RhoQ2_comp)/S_arr(i,j,k,Rho_comp) : 0.0;
+                    Real qv = (use_moisture && (ncomp > RhoQ1_comp)) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0.0;
+                    Real qc = (use_moisture && (ncomp > RhoQ2_comp)) ? S_arr(i,j,k,RhoQ2_comp)/S_arr(i,j,k,Rho_comp) : 0.0;
                     Real T = getTgivenRandRTh(S_arr(i,j,k,Rho_comp), S_arr(i,j,k,RhoTheta_comp), qv);
                     Real pressure = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp), qv);
                     Real fac = Cp_d + Cp_l*(qv + qc);
@@ -683,7 +686,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
 
         // NOTE: Protect against accessing non-existent data
         if (use_moisture) {
-            int n_qstate   = micro.Get_Qstate_Size();
+            int n_qstate   = micro->Get_Qstate_Size();
 
             // Non-precipitating components
             //--------------------------------------------------------------------------
diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp
index e598f304e..6bdf506c7 100644
--- a/Source/Initialization/ERF_init_custom.cpp
+++ b/Source/Initialization/ERF_init_custom.cpp
@@ -100,7 +100,7 @@ ERF::init_custom (int lev)
     }
 
     if (solverChoice.moisture_type != MoistureType::None) {
-        int qstate_size = micro.Get_Qstate_Size();
+        int qstate_size = micro->Get_Qstate_Size();
         MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp,    RhoQ1_comp,    1, cons_pert.nGrow());
         MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ2_comp,    RhoQ2_comp,    1, cons_pert.nGrow());
         for (int q_offset(2); q_offset<qstate_size; ++q_offset) {
diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp
index 6b460d49c..8573f2248 100644
--- a/Source/Initialization/ERF_init_from_wrfinput.cpp
+++ b/Source/Initialization/ERF_init_from_wrfinput.cpp
@@ -150,7 +150,7 @@ ERF::init_from_wrfinput (int lev)
     }
 
     auto& lev_new = vars_new[lev];
-     int n_qstate = micro.Get_Qstate_Size();
+    int n_qstate = micro->Get_Qstate_Size();
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
@@ -320,11 +320,15 @@ init_state_from_wrfinput (int lev,
 
         if (moisture_type != MoistureType::None)
         {
-            state_fab.template copy<RunOn::Device>(NC_QVAPOR_fab[idx], 0, RhoQ1_comp, 1);
-            state_fab.template mult<RunOn::Device>(NC_rho_fab[idx]   , 0, RhoQ1_comp, 1);
+            if (n_qstate >= 1) {
+              state_fab.template copy<RunOn::Device>(NC_QVAPOR_fab[idx], 0, RhoQ1_comp, 1);
+              state_fab.template mult<RunOn::Device>(NC_rho_fab[idx]   , 0, RhoQ1_comp, 1);
+            }
 
-            state_fab.template copy<RunOn::Device>(NC_QCLOUD_fab[idx], 0, RhoQ2_comp, 1);
-            state_fab.template mult<RunOn::Device>(NC_rho_fab[idx]   , 0, RhoQ2_comp, 1);
+            if (n_qstate >= 2) {
+              state_fab.template copy<RunOn::Device>(NC_QCLOUD_fab[idx], 0, RhoQ2_comp, 1);
+              state_fab.template mult<RunOn::Device>(NC_rho_fab[idx]   , 0, RhoQ2_comp, 1);
+            }
 
             if (n_qstate >= 3) {
                 state_fab.template copy<RunOn::Device>(NC_QRAIN_fab[idx], 0, RhoQ3_comp, 1);
diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H
new file mode 100644
index 000000000..c15e5e1e5
--- /dev/null
+++ b/Source/Microphysics/EulerianMicrophysics.H
@@ -0,0 +1,125 @@
+/*! @file EulerianMicrophysics.H
+ *  \brief Containes the Eulerian microphysics class
+ */
+
+#ifndef EULERIANMICROPHYSICS_H
+#define EULERIANMICROPHYSICS_H
+
+#include "NullMoist.H"
+#include "SAM.H"
+#include "Kessler.H"
+#include "Microphysics.H"
+
+/*! \brief Eulerian microphysics interface */
+class EulerianMicrophysics : public Microphysics {
+
+public:
+
+    /*! \brief Null constructor */
+    EulerianMicrophysics () { }
+
+    /*! \brief default destructor */
+    ~EulerianMicrophysics () = default;
+
+    /*! \brief Constructor: create the moisture model */
+    EulerianMicrophysics ( const int& nlev, /*!< Number of AMR levels */
+                           const MoistureType& a_model_type /*!< moisture model */)
+    {
+        AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Eulerian );
+        m_moist_model.resize(nlev);
+        if (a_model_type == MoistureType::SAM) {
+            SetModel<SAM>();
+            amrex::Print() << "SAM moisture model!\n";
+        } else if (a_model_type == MoistureType::Kessler or
+                   a_model_type == MoistureType::Kessler_NoRain) {
+            SetModel<Kessler>();
+            amrex::Print() << "Kessler moisture model!\n";
+        } else if (a_model_type == MoistureType::None) {
+            SetModel<NullMoist>();
+            amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n";
+        } else {
+            amrex::Abort("EulerianMicrophysics: Dont know this moisture_type!") ;
+        }
+    }
+
+    /*! \brief Define the moisture model */
+    void Define (const int& lev, /*!< AMR level */
+                 SolverChoice& sc /*!< Solver choice object */) override
+    {
+        m_moist_model[lev]->Define(sc);
+    }
+
+    /*! \brief Initialize the moisture model */
+    void Init ( const int& lev, /*!< AMR level */
+                const amrex::MultiFab& cons_in, /*!< Conserved state variables */
+                const amrex::BoxArray& grids,   /*!< Grids */
+                const amrex::Geometry& geom,    /*!< Geometry */
+                const amrex::Real& dt_advance   /*!< Time step */ ) override
+    {
+        m_moist_model[lev]->Init(cons_in, grids, geom, dt_advance);
+    }
+
+    /*! \brief Advance the moisture model for one time step */
+    void Advance (  const int& lev, /*!< AMR level */
+                    const amrex::Real& dt_advance, /*!< Time step */
+                    const SolverChoice &solverChoice, /*!< Solver choice object */
+                    amrex::Vector<amrex::Vector<amrex::MultiFab>>&, /*!< Dycore state variables */
+                    const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& /*!< terrain */) override
+    {
+        m_moist_model[lev]->Advance(dt_advance, solverChoice);
+    }
+
+    /*! \brief compute diagnostics */
+    void Diagnose (const int& lev /*!< AMR level */) override
+    {
+        m_moist_model[lev]->Diagnose();
+    }
+
+    /*! \brief update microphysics variables from ERF state variables */
+    void Update_Micro_Vars_Lev (const int& lev, /*! AMR level */
+                                amrex::MultiFab& cons_in /*!< Conserved state variables */) override
+    {
+        m_moist_model[lev]->Update_Micro_Vars(cons_in);
+    }
+
+    /*! \brief update ERF state variables from microphysics variables */
+    void Update_State_Vars_Lev (const int& lev, /*!< AMR level */
+                                amrex::MultiFab& cons_in /*!< Conserved state variables */) override
+    {
+        m_moist_model[lev]->Update_State_Vars(cons_in);
+    }
+
+    /*! \brief get pointer to a moisture variable */
+    amrex::MultiFab* Get_Qmoist_Ptr (const int& lev, /*!< AMR level */
+                                     const int& varIdx /*!< moisture variable index */) override
+    {
+        return m_moist_model[lev]->Qmoist_Ptr(varIdx);
+    }
+
+    /*! \brief get the number of moisture model variables */
+    int Get_Qmoist_Size (const int& /* lev */) override
+    {
+        return m_moist_model[0]->Qmoist_Size();
+    }
+
+    /*! \brief get the number of moisture-model-related conserved state variables */
+    int Get_Qstate_Size () override
+    {
+        return m_moist_model[0]->Qstate_Size();
+    }
+
+protected:
+
+    /*! \brief Create and set the specified moisture model */
+    template<class NewMoistModel>
+    void SetModel ()
+    {
+        for (int lev(0); lev<m_moist_model.size(); ++lev) {
+            m_moist_model[lev] = std::make_unique<NewMoistModel>();
+        }
+    }
+
+private:
+    amrex::Vector<std::unique_ptr<NullMoist>> m_moist_model; /*!< moisture model */
+};
+#endif
diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H
index 916f5084b..3d033bc5d 100644
--- a/Source/Microphysics/Kessler/Kessler.H
+++ b/Source/Microphysics/Kessler/Kessler.H
@@ -105,7 +105,8 @@ public:
 
     // wrapper to advance micro vars
     void
-    Advance (const amrex::Real& dt_advance, const SolverChoice &solverChoice) override
+    Advance (const amrex::Real& dt_advance,
+             const SolverChoice& solverChoice) override
     {
         dt = dt_advance;
 
diff --git a/Source/Microphysics/LagrangianMicrophysics.H b/Source/Microphysics/LagrangianMicrophysics.H
new file mode 100644
index 000000000..3f08452b3
--- /dev/null
+++ b/Source/Microphysics/LagrangianMicrophysics.H
@@ -0,0 +1,141 @@
+/*! @file LagrangianMicrophysics.H
+ *  \brief Containes the Lagrangian microphysics class
+ */
+
+#ifndef LAGRANGIANMICROPHYSICS_H
+#define LAGRANGIANMICROPHYSICS_H
+
+#ifdef ERF_USE_PARTICLES
+
+#include <utility>
+#include <string>
+
+#include "NullMoistLagrangian.H"
+#include "Microphysics.H"
+
+/* forward declaration */
+class ERFPC;
+
+/*! \brief Eulerian microphysics interface
+ *
+ * One key difference from #EulerianMicrophysics is that only the base AMR
+ * level has the moisture model. Thus, for higher AMR levels, the microphysics
+ * interface need not do anything. The number of conserved state variables for
+ * moisture (RhoQ1, RhoQ2, ...) are the same for all AMR levels; however, for
+ * levels other than the base, they just contain and evolve zeros. */
+class LagrangianMicrophysics : public Microphysics {
+
+public:
+
+    /*! \brief Null constructor */
+    LagrangianMicrophysics () { }
+
+    /*! \brief default destructor */
+    ~LagrangianMicrophysics () = default;
+
+    /*! \brief Constructor: create the moisture model */
+    LagrangianMicrophysics (  const int& /* nlev */, /*!< Number of AMR levels */
+                              const MoistureType& a_model_type /*!< moisture model */ )
+    {
+        AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Lagrangian );
+        amrex::Abort("No Lagrangian moisture model implemented yet!") ;
+    }
+
+    /*! \brief Define the moisture model */
+    void Define (const int& lev, /*!< AMR level */
+                 SolverChoice& sc /*!< Solver choice object */) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Define(sc);
+    }
+
+    /*! \brief Initialize the moisture model */
+    void Init ( const int& lev, /*!< AMR level */
+                const amrex::MultiFab& cons_in, /*!< Conserved state variables */
+                const amrex::BoxArray& grids,   /*!< Grids */
+                const amrex::Geometry& geom,    /*!< Geometry */
+                const amrex::Real& dt_advance   /*!< Time step */ ) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Init(cons_in, grids, geom, dt_advance);
+    }
+
+    /*! \brief Advance the moisture model for one time step */
+    void Advance (  const int& lev, /*!< AMR level */
+                    const amrex::Real& dt_advance, /*!< Time step */
+                    const SolverChoice &solverChoice, /*!< Solver choice object */
+                    amrex::Vector<amrex::Vector<amrex::MultiFab>>& a_vars, /*!< Dycore state variables */
+                    const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& a_z/*!< terrain */) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Advance(dt_advance, a_vars, a_z);
+    }
+
+    /*! \brief compute diagnostics */
+    void Diagnose (const int& lev /*!< AMR level */) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Diagnose();
+    }
+
+    /*! \brief update microphysics variables from ERF state variables */
+    void Update_Micro_Vars_Lev (const int& lev, /*! AMR level */
+                                amrex::MultiFab& cons_in /*!< Conserved state variables */) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Update_Micro_Vars(cons_in);
+    }
+
+    /*! \brief update ERF state variables from microphysics variables */
+    void Update_State_Vars_Lev (const int& lev, /*!< AMR level */
+                                amrex::MultiFab& cons_in /*!< Conserved state variables */) override
+    {
+        if (lev > 0) return;
+        m_moist_model->Update_State_Vars(cons_in);
+    }
+
+    /*! \brief get pointer to a moisture variable */
+    amrex::MultiFab* Get_Qmoist_Ptr (const int& lev, /*!< AMR level */
+                                     const int& varIdx /*!< moisture variable index */) override
+    {
+        return (lev > 0 ? nullptr : m_moist_model->Qmoist_Ptr(varIdx));
+    }
+
+    /*! \brief get the number of moisture model variables */
+    int Get_Qmoist_Size (const int& a_lev /*!< AMR level */) override
+    {
+        return (a_lev > 0 ? 0 : m_moist_model->Qmoist_Size());
+    }
+
+    /*! \brief get the number of moisture-model-related conserved state variables */
+    int Get_Qstate_Size () override
+    {
+        return m_moist_model->Qstate_Size();
+    }
+
+    /*! \brief get the particle container from the moisture model */
+    inline ERFPC* getParticleContainer() const
+    {
+        return m_moist_model->getParticleContainer();
+    }
+
+    /*! \brief get the name of the moisture model's particle container */
+    inline const std::string& getName() const
+    {
+        return m_moist_model->getName();
+    }
+
+protected:
+
+    /*! \brief Create and set the specified moisture model */
+    template<class NewMoistModel>
+    void SetModel ()
+    {
+        m_moist_model = std::make_unique<NewMoistModel>();
+    }
+
+    std::unique_ptr<NullMoistLagrangian> m_moist_model; /*!< moisture model */
+};
+
+#endif
+#endif
diff --git a/Source/Microphysics/Make.package b/Source/Microphysics/Make.package
index 36572d9e5..d2b85471d 100644
--- a/Source/Microphysics/Make.package
+++ b/Source/Microphysics/Make.package
@@ -1,2 +1,4 @@
 CEXE_headers += Microphysics.H
+CEXE_headers += EulerianMicrophysics.H
+CEXE_headers += LagrangianMicrophysics.H
 
diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H
index 36ade90c5..aa2c28b63 100644
--- a/Source/Microphysics/Microphysics.H
+++ b/Source/Microphysics/Microphysics.H
@@ -1,82 +1,73 @@
+/*! @file Microphysics.H
+ *  \brief Containes the base class for microphysics
+ */
+
 #ifndef MICROPHYSICS_H
 #define MICROPHYSICS_H
 
-#include <NullMoist.H>
-#include <SAM.H>
-#include <Kessler.H>
-#include <ERF.H>
+#include "DataStruct.H"
 
+/*! \brief Base class for microphysics interface */
 class Microphysics {
 
 public:
 
+    /*! \brief Null constructor */
     Microphysics () { }
 
-    ~Microphysics () = default;
+    /*! \brief default destructor */
+    virtual ~Microphysics () = default;
 
-    void
-    ReSize (const int& nlev) { m_moist_model.resize(nlev); }
+    /*! \brief define the microphysics object */
+    virtual void Define (const int&, SolverChoice&) = 0;
 
-    template<class NewMoistModel>
-    void
-    SetModel ()
-    {
-        for (int lev(0); lev<m_moist_model.size(); ++lev) {
-            m_moist_model[lev] = std::make_unique<NewMoistModel>();
-        }
-    }
+    /*! \brief initialize the microphysics object */
+    virtual void Init (const int&,
+                       const amrex::MultiFab&,
+                       const amrex::BoxArray&,
+                       const amrex::Geometry&,
+                       const amrex::Real&) = 0;
 
-    void
-    Define (const int& lev,
-            SolverChoice& sc)
-    {
-        m_moist_model[lev]->Define(sc);
-    }
+    /*! \brief advance microphysics for one time step */
+    virtual void Advance (  const int&,
+                            const amrex::Real&,
+                            const SolverChoice&,
+                            amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
+                            const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& ) = 0;
 
-    void
-    Init (const int& lev,
-          const amrex::MultiFab& cons_in,
-          const amrex::BoxArray& grids,
-          const amrex::Geometry& geom,
-          const amrex::Real& dt_advance)
-    {
-        m_moist_model[lev]->Init(cons_in, grids, geom, dt_advance);
-    }
+    /*! \brief compute diagnostics */
+    virtual void Diagnose (const int&) = 0;
 
-    void
-    Advance (const int& lev, const amrex::Real& dt_advance, const SolverChoice &solverChoice)
-    {
-        m_moist_model[lev]->Advance(dt_advance, solverChoice);
-    }
+    /*! \brief update microphysics variables from ERF state variables */
+    virtual void Update_Micro_Vars_Lev (const int&, amrex::MultiFab&) = 0;
 
-    void
-    Diagnose (const int& lev)
-    {
-        m_moist_model[lev]->Diagnose();
-    }
+    /*! \brief update ERF state variables from microphysics variables */
+    virtual void Update_State_Vars_Lev (const int&, amrex::MultiFab&) = 0;
 
-    void
-    Update_Micro_Vars_Lev (const int& lev, amrex::MultiFab& cons_in)
-    {
-        m_moist_model[lev]->Update_Micro_Vars(cons_in);
-    }
-
-    void
-    Update_State_Vars_Lev (const int& lev, amrex::MultiFab& cons_in)
-    {
-        m_moist_model[lev]->Update_State_Vars(cons_in);
-    }
+    /*! \brief get pointer to a moisture variable */
+    virtual amrex::MultiFab* Get_Qmoist_Ptr (const int&, const int&) = 0;
 
-    amrex::MultiFab*
-    Get_Qmoist_Ptr (const int& lev, const int& varIdx) { return m_moist_model[lev]->Qmoist_Ptr(varIdx); }
+    /*! \brief get the number of moisture model variables */
+    virtual int Get_Qmoist_Size (const int&) = 0;
 
-    int
-    Get_Qmoist_Size () { return m_moist_model[0]->Qmoist_Size(); }
+    /*! \brief get the number of moisture-model-related conserved state variables */
+    virtual int Get_Qstate_Size () = 0;
 
-    int
-    Get_Qstate_Size () { return m_moist_model[0]->Qstate_Size(); }
+    /*! \brief query if a specified moisture model is Eulerian or Lagrangian */
+    static MoistureModelType modelType(const MoistureType a_moisture_type)
+    {
+        if (    (a_moisture_type == MoistureType::SAM)
+             || (a_moisture_type == MoistureType::Kessler)
+             || (a_moisture_type == MoistureType::Kessler_NoRain)
+             || (a_moisture_type == MoistureType::None) ) {
+            return MoistureModelType::Eulerian;
+        } else {
+            amrex::Abort("Dont know this moisture_type!") ;
+            return MoistureModelType::Undefined;
+        }
+    }
 
 private:
-    amrex::Vector<std::unique_ptr<NullMoist>> m_moist_model;
+
 };
 #endif
diff --git a/Source/Microphysics/Null/Make.package b/Source/Microphysics/Null/Make.package
index d528e5844..d918b542d 100644
--- a/Source/Microphysics/Null/Make.package
+++ b/Source/Microphysics/Null/Make.package
@@ -1,2 +1,3 @@
 CEXE_headers += NullMoist.H
+CEXE_headers += NullMoistLagrangian.H
 
diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H
index 2ba068543..d1a98f325 100644
--- a/Source/Microphysics/Null/NullMoist.H
+++ b/Source/Microphysics/Null/NullMoist.H
@@ -7,7 +7,7 @@
 
 class NullMoist {
 
- public:
+public:
     NullMoist () {}
 
     virtual ~NullMoist () = default;
@@ -59,8 +59,9 @@ class NullMoist {
     int
     Qstate_Size () { return NullMoist::m_qstate_size; }
 
- private:
+private:
     int m_qmoist_size = 1;
     int m_qstate_size = 0;
 };
+
 #endif
diff --git a/Source/Microphysics/Null/NullMoistLagrangian.H b/Source/Microphysics/Null/NullMoistLagrangian.H
new file mode 100644
index 000000000..ed4bfb062
--- /dev/null
+++ b/Source/Microphysics/Null/NullMoistLagrangian.H
@@ -0,0 +1,54 @@
+/*! @file NullMoistLagrangian.H
+ *  \brief Contains the Lagrangian moisture model base class */
+
+#ifndef NULLMOISTLAGRANGIAN_H
+#define NULLMOISTLAGRANGIAN_H
+
+#ifdef ERF_USE_PARTICLES
+
+#include "NullMoist.H"
+
+/* forward declaration */
+class ERFPC;
+
+/*! \brief Base class for a Lagrangian moisture model
+ *
+ * Extends #NullMoist for a Lagrangian model */
+class NullMoistLagrangian : public NullMoist {
+
+public:
+
+    /*! \brief Null constructor */
+    NullMoistLagrangian() {}
+
+    /*! \brief Default destructor */
+    virtual ~NullMoistLagrangian() = default;
+
+    /*! \brief get the particle container */
+    virtual ERFPC* getParticleContainer()
+    {
+        return nullptr;
+    }
+
+    /*! \brief get the name */
+    virtual const std::string& getName() const
+    {
+        return m_name;
+    }
+
+    using NullMoist::Advance;
+
+    /*! \brief advance the moisture model by one time step */
+    virtual void
+    Advance (const amrex::Real&, /* dt */
+             amrex::Vector<amrex::Vector<amrex::MultiFab>>&, /* state variables */
+             const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& /* terrain */) { }
+
+protected:
+
+private:
+    const std::string m_name = "null"; /*!< name of the moisture model */
+};
+
+#endif
+#endif
diff --git a/Source/Microphysics/SAM/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp
index cbb57b214..80b5297a5 100644
--- a/Source/Microphysics/SAM/IceFall.cpp
+++ b/Source/Microphysics/SAM/IceFall.cpp
@@ -5,15 +5,18 @@
 using namespace amrex;
 
 /**
- * Sedimentation of cloud ice (A32/33); sources enthalpy
+ * Sedimentation of cloud ice (A32)
  */
 void SAM::IceFall () {
 
-    Real dz  = m_geom.CellSize(2);
-    Real dtn = dt;
-    int  nz  = nlev;
+    Real dz   = m_geom.CellSize(2);
+    Real dtn  = dt;
+    Real coef = dtn/dz;
+
+    auto domain = m_geom.Domain();
+    int k_lo = domain.smallEnd(2);
+    int k_hi = domain.bigEnd(2);
 
-    int kmax, kmin;
     auto qcl   = mic_fab_vars[MicVar::qcl];
     auto qci   = mic_fab_vars[MicVar::qci];
     auto qn    = mic_fab_vars[MicVar::qn];
@@ -28,97 +31,56 @@ void SAM::IceFall () {
     fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ng);
     fz.setVal(0.);
 
-    kmin = nz;
-    kmax = -1;
+    for (MFIter mfi(fz, TileNoZ()); mfi.isValid(); ++mfi) {
+        auto qci_array = qci->array(mfi);
+        auto rho_array = rho->array(mfi);
+        auto fz_array  = fz.array(mfi);
 
-    { // calculate maximum and minium ice fall vertical region
-        auto const& qcl_arrays  = qcl->const_arrays();
-        auto const& qci_arrays  = qci->const_arrays();
-        auto const& tabs_arrays = tabs->const_arrays();
+        const auto& box3d  = mfi.tilebox();
 
-        GpuTuple<int, int> k_max_min = ParReduce(TypeList<ReduceOpMax, ReduceOpMin>{},
-                                                 TypeList<int, int>{},
-                                                 *qcl, IntVect::TheZeroVector(),
-                                                 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-                                                 -> GpuTuple<int, int>
-                                                 {
-                                                     bool is_positive       = qcl_arrays[box_no](i,j,k)+qci_arrays[box_no](i,j,k) > 0.0;
-                                                     bool smaller_than_zero = tabs_arrays[box_no](i,j,k) < 273.15;
-                                                     int mkmin = is_positive && smaller_than_zero ? k : nz-1;
-                                                     int mkmax = is_positive && smaller_than_zero ? k : 0;
-                                                     return {mkmax, mkmin};
-                                                 });
-        kmax = amrex::get<0>(k_max_min);
-        kmin = amrex::get<1>(k_max_min);
+        ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            Real rho_avg, qci_avg;
+            if (k==k_lo) {
+                rho_avg = rho_array(i,j,k);
+                qci_avg = qci_array(i,j,k);
+            } else if (k==k_hi+1) {
+                rho_avg = rho_array(i,j,k-1);
+                qci_avg = qci_array(i,j,k-1);
+            } else {
+                rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
+                qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k));
+            }
+            Real vt_ice = min( 0.4 , 8.66 * pow( (max(0.,qci_avg)+1.e-10) , 0.24) );
+            fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
+        });
     }
 
-    for ( amrex::MFIter mfi(*tabs, TileNoZ()); mfi.isValid(); ++mfi) {
+    for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) {
         auto qci_array   = qci->array(mfi);
         auto qn_array    = qn->array(mfi);
         auto qt_array    = qt->array(mfi);
         auto rho_array   = rho->array(mfi);
         auto fz_array    = fz.array(mfi);
 
-        const auto& gbox3d = mfi.tilebox(IntVect(0),IntVect(0,0,1));
         const auto& box3d  = mfi.tilebox();
 
-        ParallelFor(gbox3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-        {
-            if (k >= std::max(0,kmin-1) && k <= kmax+1 ) {
-                // Set up indices for x-y planes above and below current plane.
-                int kc = std::min(k+1, nz-1);
-                int kb = std::max(k-1, 0);
-
-                // CFL number based on grid spacing interpolated to interface i,j,k-1/2
-                Real coef = dtn/dz;
-
-                // Compute cloud ice density in this cell and the ones above/below.
-                // Since cloud ice is falling, the above cell is u(icrm,upwind),
-                // this cell is c (center) and the one below is d (downwind).
-                Real qiu = rho_array(i,j,kc)*qci_array(i,j,kc);
-                Real qic = rho_array(i,j,k )*qci_array(i,j,k );
-                Real qid = rho_array(i,j,kb)*qci_array(i,j,kb);
-
-                // Ice sedimentation velocity depends on ice content. The fiting is
-                // based on the data by Heymsfield (JAS,2003). -Marat
-                Real vt_ice = min( 0.4 , 8.66 * pow( (max(0.,qic)+1.e-10) , 0.24) );   // Heymsfield (JAS, 2003, p.2607)
-
-                // Use MC flux limiter in computation of flux correction.
-                // (MC = monotonized centered difference).
-                Real tmp_phi;
-                if ( std::abs(qic-qid) < 1.0e-25 ) {  // when qic, and qid is very small, qic_qid can still be zero
-                    // even if qic is not equal to qid. so add a fix here +++mhwang
-                    tmp_phi = 0.;
-                } else {
-                    Real tmp_theta = (qiu-qic)/(qic-qid+1.0e-20);
-                    tmp_phi = max(0., min(0.5*(1.+tmp_theta), min(2., 2.*tmp_theta)));
-                }
-
-                // Compute limited flux.
-                // Since falling cloud ice is a 1D advection problem, this
-                // flux-limited advection scheme is monotonic.
-                fz_array(i,j,k) = -vt_ice*(qic - 0.5*(1.-coef*vt_ice)*tmp_phi*(qic-qid));
-            }
-        });
-
         ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
-            if ( k >= std::max(0,kmin) && k <= kmax ) {
-                //==================================================
-                // Cloud ice sedimentation (A32)
-                //==================================================
-                Real coef = dtn/dz;
-                Real dqi  = std::max(-qci_array(i,j,k),coef*(fz_array(i,j,k)-fz_array(i,j,k+1)));
-
-                // Add this increment to both non-precipitating and total water.
-                qci_array(i,j,k) += dqi;
-                 qn_array(i,j,k) += dqi;
-                 qt_array(i,j,k) += dqi;
-
-                // NOTE: Sedimentation does not affect the potential temperature,
-                //       but it does affect the liquid/ice static energy.
-                //       No source to Theta occurs here.
-            }
+            //==================================================
+            // Cloud ice sedimentation (A32)
+            //==================================================
+            Real dqi  = (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
+            dqi = std::max(-qci_array(i,j,k), dqi);
+
+            // Add this increment to both non-precipitating and total water.
+            qci_array(i,j,k) += dqi;
+            qn_array(i,j,k) += dqi;
+            qt_array(i,j,k) += dqi;
+
+            // NOTE: Sedimentation does not affect the potential temperature,
+            //       but it does affect the liquid/ice static energy.
+            //       No source to Theta occurs here.
         });
     }
 }
diff --git a/Source/Microphysics/SAM/PrecipFall.cpp b/Source/Microphysics/SAM/PrecipFall.cpp
index 78bd321b5..9e6a764a4 100644
--- a/Source/Microphysics/SAM/PrecipFall.cpp
+++ b/Source/Microphysics/SAM/PrecipFall.cpp
@@ -11,10 +11,9 @@ using namespace amrex;
  *
  * @param[in] hydro_type Type selection for the precipitation advection hydrodynamics scheme (0-3)
  */
-void SAM::PrecipFall (int hydro_type)
+void SAM::PrecipFall ()
 {
     Real eps = std::numeric_limits<Real>::epsilon();
-    bool constexpr nonos = true;
     Real rho_0 = 1.29;
 
     Real gamr3 = erf_gammafff(4.0+b_rain);
@@ -25,14 +24,18 @@ void SAM::PrecipFall (int hydro_type)
     Real vsnow = (a_snow*gams3/6.0)*pow((PI*rhos*nzeros),-csnow);
     Real vgrau = (a_grau*gamg3/6.0)*pow((PI*rhog*nzerog),-cgrau);
 
-    Real dt_advance = dt;
-    int nz = nlev;
+    auto dz   = m_geom.CellSize(2);
+    Real dtn  = dt;
+    Real coef = dtn/dz;
+
+    auto domain = m_geom.Domain();
+    int k_lo = domain.smallEnd(2);
+    int k_hi = domain.bigEnd(2);
 
     auto qpr   = mic_fab_vars[MicVar::qpr];
     auto qps   = mic_fab_vars[MicVar::qps];
     auto qpg   = mic_fab_vars[MicVar::qpg];
     auto qp    = mic_fab_vars[MicVar::qp];
-    auto omega = mic_fab_vars[MicVar::omega];
     auto rho   = mic_fab_vars[MicVar::rho];
     auto tabs  = mic_fab_vars[MicVar::tabs];
     auto theta = mic_fab_vars[MicVar::theta];
@@ -41,266 +44,81 @@ void SAM::PrecipFall (int hydro_type)
     auto dm    = tabs->DistributionMap();
     auto ngrow = tabs->nGrowVect();
 
-    MultiFab mx;
-    MultiFab mn;
-    MultiFab lfac;
-    MultiFab www;
     MultiFab fz;
-    MultiFab wp;
-    MultiFab tmp_qp;
-    MultiFab prec_cfl_fab;
-
-    mx.define(ba,dm, 1, ngrow);
-    mn.define(ba,dm, 1, ngrow);
-    lfac.define(ba, dm, 1, ngrow);
-    www.define(ba, dm, 1, ngrow);
-    fz.define(ba, dm, 1, ngrow);
-    wp.define(ba, dm, 1, ngrow);
-    tmp_qp.define(ba, dm, 1, ngrow);
-    prec_cfl_fab.define(ba, dm, 1, ngrow);
-
-    TableData<Real, 1> irho;
-    TableData<Real, 1> iwmax;
-    TableData<Real, 1> rhofac;
-
-    irho.resize({zlo},{zhi});
-    iwmax.resize({zlo},{zhi});
-    rhofac.resize({zlo},{zhi});
-
-    auto irho_t   = irho.table();
-    auto iwmax_t  = iwmax.table();
-    auto rhofac_t = rhofac.table();
-    auto rho1d_t  = rho1d.table();
-
-    auto dz = m_geom.CellSize(2);
-
-    // Precompute omega variable
-    for ( MFIter mfi(*(mic_fab_vars[MicVar::omega]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-        auto omega_array = mic_fab_vars[MicVar::omega]->array(mfi);
-        auto tabs_array  = mic_fab_vars[MicVar::tabs]->array(mfi);
-
-        const auto& box3d = mfi.tilebox();
-
-        ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-        {
-            omega_array(i,j,k) = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
-        });
-    }
-
-    ParallelFor(nz, [=] AMREX_GPU_DEVICE (int k) noexcept
-    {
-        rhofac_t(k)  = std::sqrt(rho_0/rho1d_t(k));
-        irho_t(k)    = 1.0/rho1d_t(k);
-        Real wmax    = dz/dt_advance;   // Velocity equivalent to a cfl of 1.0.
-        iwmax_t(k)   = 1.0/wmax;
-    });
+    fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow);
 
     //  Add sedimentation of precipitation field to the vert. vel.
-    for ( MFIter mfi(lfac, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-        auto lfac_array     = lfac.array(mfi);
-        auto omega_array    = omega->array(mfi);
-        auto qp_array       = qp->array(mfi);
-        auto rho_array      = rho->array(mfi);
-        auto tabs_array     = tabs->array(mfi);
-        auto wp_array       = wp.array(mfi);
-        auto www_array      = www.array(mfi);
-        auto fz_array       = fz.array(mfi);
-        auto prec_cfl_array = prec_cfl_fab.array(mfi);
+    for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        auto qp_array   = qp->array(mfi);
+        auto rho_array  = rho->array(mfi);
+        auto tabs_array = tabs->array(mfi);
+        auto fz_array   = fz.array(mfi);
 
         const auto& box3d = mfi.tilebox();
 
-        const Real fac_cond = m_fac_cond;
-        const Real fac_sub  = m_fac_sub;
-        const Real fac_fus  = m_fac_fus;
-
         ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
-            if (hydro_type == 0) {
-                lfac_array(i,j,k) = fac_cond;
+            Real rho_avg, tab_avg, qp_avg;
+            if (k==k_lo) {
+                rho_avg =  rho_array(i,j,k);
+                tab_avg = tabs_array(i,j,k);
+                 qp_avg =   qp_array(i,j,k);
+            } else if (k==k_hi+1) {
+                rho_avg =  rho_array(i,j,k-1);
+                tab_avg = tabs_array(i,j,k-1);
+                 qp_avg =   qp_array(i,j,k-1);
+            } else {
+                rho_avg = 0.5*( rho_array(i,j,k-1) +  rho_array(i,j,k));
+                tab_avg = 0.5*(tabs_array(i,j,k-1) + tabs_array(i,j,k));
+                 qp_avg = 0.5*(  qp_array(i,j,k-1) +   qp_array(i,j,k));
             }
-            else if (hydro_type == 1) {
-                lfac_array(i,j,k) = fac_sub;
-            }
-            else if (hydro_type == 2) {
-                lfac_array(i,j,k) = fac_cond + (1.0-omega_array(i,j,k))*fac_fus;
-            }
-            else if (hydro_type == 3) {
-                lfac_array(i,j,k) = 0.0;
-            }
-            Real Veff = term_vel_qp(qp_array(i,j,k), vrain, vsnow, vgrau,
-                                    rho_array(i,j,k), tabs_array(i,j,k));
-            wp_array(i,j,k) = -Veff * std::sqrt(rho_0/rho_array(i,j,k));
-            prec_cfl_array(i,j,k) = wp_array(i,j,k) * iwmax_t(k);
-            wp_array(i,j,k) *= rho_array(i,j,k) * dt_advance/dz;
-            if (k == 0) {
-                fz_array(i,j,nz-1)   = 0.0;
-                www_array(i,j,nz-1)  = 0.0;
-                lfac_array(i,j,nz-1) = 0.0;
+
+            Real Pprecip = 0.0;
+            if(qp_avg > qp_threshold) {
+                Real omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr));
+                Real omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr));
+                Real qrr = omp*qp_avg;
+                Real qss = (1.0-omp)*(1.0-omg)*qp_avg;
+                Real qgg = (1.0-omp)*(omg)*qp_avg;
+                Pprecip = omp*vrain*std::pow(rho_avg*qrr,1.0+crain)
+                        + (1.0-omp)*( (1.0-omg)*vsnow*std::pow(rho_avg*qss,1.0+csnow)
+                                    +      omg *vgrau*std::pow(rho_avg*qgg,1.0+cgrau) );
             }
+            fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg);
         });
     }
 
-    auto const& cfl_arrays = prec_cfl_fab.const_arrays();
-    Real prec_cfl = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{},
-                              prec_cfl_fab, IntVect(0),
-                              [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-                              ->GpuTuple<Real>
-    {
-        return { cfl_arrays[box_no](i,j,k) };
-    });
-
-    // If maximum CFL due to precipitation velocity is greater than 0.9,
-    // take more than one advection step to maintain stability.
-    int nprec;
-    if (prec_cfl > 0.9) {
-        nprec = static_cast<int>(std::ceil(prec_cfl/0.9));
-        for (MFIter mfi(wp, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-            auto wp_array = wp.array(mfi);
-            const auto& box3d = mfi.tilebox();
+    for (MFIter mfi(*qp, TileNoZ()); mfi.isValid(); ++mfi) {
+        auto qpr_array    = qpr->array(mfi);
+        auto qps_array    = qps->array(mfi);
+        auto qpg_array    = qpg->array(mfi);
+        auto qp_array     = qp->array(mfi);
+        auto rho_array    = rho->array(mfi);
+        auto tabs_array   = tabs->array(mfi);
+        auto fz_array     = fz.array(mfi);
 
-            ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i)
-            {
-                // wp already includes factor of dt, so reduce it by a
-                // factor equal to the number of precipitation steps.
-                wp_array(i,j,k) = wp_array(i,j,k)/Real(nprec);
-            });
-        }
-    } else {
-        nprec = 1;
-    }
-
-#ifdef ERF_FIXED_SUBCYCLE
-    nprec = 4;
-#endif
-
-    for(int iprec = 1; iprec<=nprec; iprec++) {
-        for ( MFIter mfi(tmp_qp, TileNoZ()); mfi.isValid(); ++mfi) {
-            auto qpr_array    = qpr->array(mfi);
-            auto qps_array    = qps->array(mfi);
-            auto qpg_array    = qpg->array(mfi);
-            auto qp_array     = qp->array(mfi);
-            auto rho_array    = rho->array(mfi);
-            auto tabs_array   = tabs->array(mfi);
-            auto tmp_qp_array = tmp_qp.array(mfi);
-            auto mx_array     = mx.array(mfi);
-            auto mn_array     = mn.array(mfi);
-            auto fz_array     = fz.array(mfi);
-            auto wp_array     = wp.array(mfi);
-            auto lfac_array   = lfac.array(mfi);
-            auto www_array    = www.array(mfi);
-
-            const auto& box3d = mfi.tilebox();
-
-            ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                tmp_qp_array(i,j,k) = qp_array(i,j,k); // Temporary array for qp in this column
-            });
-
-            ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                if (nonos) {
-                    int kc=min(nz-1,k+1);
-                    int kb=max(0,k-1);
-                    mx_array(i,j,k) = max(tmp_qp_array(i,j,kb), max(tmp_qp_array(i,j,kc), tmp_qp_array(i,j,k)));
-                    mn_array(i,j,k) = min(tmp_qp_array(i,j,kb), min(tmp_qp_array(i,j,kc), tmp_qp_array(i,j,k)));
-                }
-                // Define upwind precipitation flux
-                fz_array(i,j,k) = tmp_qp_array(i,j,k)*wp_array(i,j,k);
-            });
-
-            ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                int kc = min(k+1, nz-1);
-                Real dqp = (fz_array(i,j,kc)-fz_array(i,j,k)) / rho_array(i,j,k);
-                tmp_qp_array(i,j,k) = tmp_qp_array(i,j,k) + dqp; //Update temporary qp
-            });
-
-            ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                // Also, compute anti-diffusive correction to previous
-                // (upwind) approximation to the flux
-                int kb=max(0,k-1);
-                // The precipitation velocity is a cell-centered quantity,
-                // since it is computed from the cell-centered
-                // precipitation mass fraction.  Therefore, a reformulated
-                // anti-diffusive flux is used here which accounts for
-                // this and results in reduced numerical diffusion.
-                www_array(i,j,k) = 0.5*(1.0+wp_array(i,j,k)/rho_array(i,j,k))*(tmp_qp_array(i,j,kb)*wp_array(i,j,kb) -
-                                                                        tmp_qp_array(i,j,k)*wp_array(i,j,k)); // works for wp(k)<0
-            });
-
-            if (nonos) {
-                ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                {
-                    int kc=min(nz-1,k+1);
-                    int kb=max(0,k-1);
-                    mx_array(i,j,k) = max(tmp_qp_array(i,j,kb),max(tmp_qp_array(i,j,kc), max(tmp_qp_array(i,j,k), mx_array(i,j,k))));
-                    mn_array(i,j,k) = min(tmp_qp_array(i,j,kb),min(tmp_qp_array(i,j,kc), min(tmp_qp_array(i,j,k), mn_array(i,j,k))));
-                    kc = min(nz-1,k+1);
-                    mx_array(i,j,k) = rho_array(i,j,k)*(mx_array(i,j,k)-tmp_qp_array(i,j,k))/(pn(www_array(i,j,kc)) +
-                                                                                        pp(www_array(i,j,k))+eps);
-                    mn_array(i,j,k) = rho_array(i,j,k)*(tmp_qp_array(i,j,k)-mn_array(i,j,k))/(pp(www_array(i,j,kc)) +
-                                                                                        pn(www_array(i,j,k))+eps);
-                });
-
-                ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                {
-                    int kb=max(0,k-1);
-                    // Add limited flux correction to fz(k).
-                    fz_array(i,j,k) = fz_array(i,j,k)
-                                    + pp(www_array(i,j,k))*std::min(1.0,std::min(mx_array(i,j,k), mn_array(i,j,kb)))
-                                    - pn(www_array(i,j,k))*std::min(1.0,std::min(mx_array(i,j,kb),mn_array(i,j,k))); // Anti-diffusive flux
-                });
-            }
-
-            // Update precipitation mass fraction and liquid-ice static
-            // energy using precipitation fluxes computed in this column.
-            ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                int kc=min(k+1, nz-1);
-                // Update precipitation mass fraction.
-                // Note that fz is the total flux, including both the
-                // upwind flux and the anti-diffusive correction.
-
-                //==================================================
-                // Precipitating sedimentation (A19)
-                //==================================================
-                Real dqp = ( fz_array(i,j,kc) - fz_array(i,j,k) ) / rho_array(i,j,k);
-                Real omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
-                Real omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr));
-
-                qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp);
-                qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg));
-                qpg_array(i,j,k) = std::max(0.0, qpg_array(i,j,k) + dqp*(1.0-omp)*omg);
-                 qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
-
-                // NOTE: Sedimentation does not affect the potential temperature,
-                //       but it does affect the liquid/ice static energy.
-                //       No source to Theta occurs here.
-            });
+        const auto& box3d = mfi.tilebox();
 
-            if (iprec < nprec) {
-                // Re-compute precipitation velocity using new value of qp.
-                ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                {
-                    Real tmp = term_vel_qp(qp_array(i,j,k),
-                                           vrain, vsnow, vgrau,
-                                           rho_array(i,j,k), tabs_array(i,j,k));
-                    wp_array(i,j,k) = std::sqrt(rho_0/rho_array(i,j,k))*tmp;
-                    // Decrease precipitation velocity by factor of nprec
-                    wp_array(i,j,k) = -wp_array(i,j,k)*rho_array(i,j,k)*dt_advance/dz/nprec;
-                    // Note: Don't bother checking CFL condition at each
-                    // substep since it's unlikely that the CFL will
-                    // increase very much between substeps when using
-                    // monotonic advection schemes.
-                    if (k == 0) {
-                          fz_array(i,j,nz-1) = 0.0;
-                         www_array(i,j,nz-1) = 0.0;
-                        lfac_array(i,j,nz-1) = 0.0;
-                    }
-                });
-            }
-        } // mfi
-    } // iprec loop
+        // Update precipitation mass fraction and liquid-ice static
+        // energy using precipitation fluxes computed in this column.
+        ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            //==================================================
+            // Precipitating sedimentation (A19)
+            //==================================================
+            Real dqp = (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
+            Real omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
+            Real omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr));
+
+            qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp);
+            qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg));
+            qpg_array(i,j,k) = std::max(0.0, qpg_array(i,j,k) + dqp*(1.0-omp)*omg);
+             qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
+
+            // NOTE: Sedimentation does not affect the potential temperature,
+            //       but it does affect the liquid/ice static energy.
+            //       No source to Theta occurs here.
+        });
+    } // mfi
 }
 
diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H
index b181f4c05..807a85b51 100644
--- a/Source/Microphysics/SAM/SAM.H
+++ b/Source/Microphysics/SAM/SAM.H
@@ -69,7 +69,7 @@ public:
     void Precip ();
 
     // precip fall
-    void PrecipFall (int hydro_type=2);
+    void PrecipFall ();
 
     // diagnose
     void Diagnose () override;
diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp
index 1399cef83..f460effd9 100644
--- a/Source/TimeIntegration/ERF_Advance.cpp
+++ b/Source/TimeIntegration/ERF_Advance.cpp
@@ -65,7 +65,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/
 
     if (solverChoice.moisture_type != MoistureType::None) {
         // TODO: This is only qv
-        FillPatchMoistVars(lev, *(qmoist[lev][0]));
+        if (qmoist[lev].size() > 0) FillPatchMoistVars(lev, *(qmoist[lev][0]));
     }
 
 #if defined(ERF_USE_WINDFARM)
diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp
index 02dee83ba..6ffe270d7 100644
--- a/Source/TimeIntegration/ERF_advance_microphysics.cpp
+++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp
@@ -7,8 +7,8 @@ void ERF::advance_microphysics (int lev,
                                 const Real& dt_advance)
 {
     if (solverChoice.moisture_type != MoistureType::None) {
-        micro.Update_Micro_Vars_Lev(lev, cons);
-        micro.Advance(lev, dt_advance, solverChoice);
-        micro.Update_State_Vars_Lev(lev, cons);
+        micro->Update_Micro_Vars_Lev(lev, cons);
+        micro->Advance(lev, dt_advance, solverChoice, vars_new, z_phys_nd);
+        micro->Update_State_Vars_Lev(lev, cons);
     }
 }
diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H
index ba950629f..3e1d534ac 100644
--- a/Source/TimeIntegration/TI_slow_rhs_fun.H
+++ b/Source/TimeIntegration/TI_slow_rhs_fun.H
@@ -95,7 +95,7 @@
             MultiFab* p0_new = &p_hse_new;
 
             make_buoyancy(S_data, S_prim, buoyancy, fine_geom,
-                          solverChoice, r0_new, micro.Get_Qstate_Size());
+                          solverChoice, r0_new, micro->Get_Qstate_Size());
 
             erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
                              xvel_new, yvel_new, zvel_new,
@@ -188,7 +188,7 @@
 
             // If not moving_terrain
             make_buoyancy(S_data, S_prim, buoyancy,
-                          fine_geom, solverChoice, r0, micro.Get_Qstate_Size());
+                          fine_geom, solverChoice, r0, micro->Get_Qstate_Size());
 
             erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
                              xvel_new, yvel_new, zvel_new,
@@ -352,7 +352,7 @@
 
         // If not moving_terrain
         make_buoyancy(S_data, S_prim, buoyancy,
-                      fine_geom, solverChoice, r0, micro.Get_Qstate_Size());
+                      fine_geom, solverChoice, r0, micro->Get_Qstate_Size());
 
         erf_slow_rhs_inc(level, nrk, slow_dt,
                          S_rhs, S_old, S_data, S_prim, S_scratch,
diff --git a/Submodules/AMReX b/Submodules/AMReX
index 52393d3fa..d737886d5 160000
--- a/Submodules/AMReX
+++ b/Submodules/AMReX
@@ -1 +1 @@
-Subproject commit 52393d3fafc835fa379f0420aa2de4ccdddf155a
+Subproject commit d737886d574d4f1c0cf76323337b666ecd5bb4e0