This compared well with Kessler. (erf-model#1634)

Co-authored-by: Aaron Lattanzi <[email protected]>

Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler

@@ -0,0 +1,81 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step  = 2000
+stop_time = 3600.0
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent = 20000.0 400.0  10000.0
+amr.n_cell           = 200     4      100
+geometry.is_periodic = 0 1 0
+xlo.type = "SlipWall"
+xhi.type = "SlipWall"    
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt = 0.5
+erf.fixed_mri_dt_ratio = 4
+#erf.no_substepping = 1
+#erf.fixed_dt = 0.1
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v              = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+erf.check_int       = 100       # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+erf.plot_int_1      = 100        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qrain
+
+# SOLVER CHOICES
+erf.use_gravity          = true
+erf.use_coriolis         = false
+
+erf.dycore_horiz_adv_type    = "Upwind_3rd"
+erf.dycore_vert_adv_type     = "Upwind_3rd"
+erf.dryscal_horiz_adv_type   = "Upwind_3rd"
+erf.dryscal_vert_adv_type    = "Upwind_3rd"
+erf.moistscal_horiz_adv_type = "Upwind_3rd"
+erf.moistscal_vert_adv_type  = "Upwind_3rd"       
+
+# PHYSICS OPTIONS
+erf.les_type        = "None"
+erf.pbl_type        = "None"
+erf.moisture_model  = "Kessler"
+erf.buoyancy_type   = 1
+erf.use_moist_background = true
+
+erf.molec_diff_type  = "ConstantAlpha"
+erf.rho0_trans       = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T          = 0.0 # [m^2/s]
+erf.alpha_C          = 0.0
+
+# INITIAL CONDITIONS
+#erf.init_type = "input_sounding"
+#erf.input_sounding_file = "BF02_moist_sounding"
+#erf.init_sounding_ideal = true
+
+# PROBLEM PARAMETERS (optional)
+# warm bubble input
+prob.x_c    = 10000.0
+prob.z_c    =  2000.0
+prob.x_r    =  2000.0
+prob.z_r    =  2000.0
+prob.T_0    =   300.0
+
+prob.do_moist_bubble = true
+prob.theta_pert  = 2.0
+prob.qt_init     = 0.02
+prob.eq_pot_temp = 320.0

Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce

@@ -0,0 +1,81 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step  = 2000
+stop_time = 3600.0
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent = 20000.0 400.0  10000.0
+amr.n_cell           = 200     4      100
+geometry.is_periodic = 0 1 0
+xlo.type = "SlipWall"
+xhi.type = "SlipWall"    
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt = 0.5
+erf.fixed_mri_dt_ratio = 4
+#erf.no_substepping = 1
+#erf.fixed_dt = 0.1
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v              = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+erf.check_int       = 100       # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+erf.plot_int_1      = 100        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qi qp qrain qsnow qgraup
+
+# SOLVER CHOICES
+erf.use_gravity          = true
+erf.use_coriolis         = false
+
+erf.dycore_horiz_adv_type    = "Upwind_3rd"
+erf.dycore_vert_adv_type     = "Upwind_3rd"
+erf.dryscal_horiz_adv_type   = "Upwind_3rd"
+erf.dryscal_vert_adv_type    = "Upwind_3rd"
+erf.moistscal_horiz_adv_type = "Upwind_3rd"
+erf.moistscal_vert_adv_type  = "Upwind_3rd"       
+
+# PHYSICS OPTIONS
+erf.les_type        = "None"
+erf.pbl_type        = "None"
+erf.moisture_model  = "SAM_NoIce"
+erf.buoyancy_type   = 1
+erf.use_moist_background = true
+
+erf.molec_diff_type  = "ConstantAlpha"
+erf.rho0_trans       = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T          = 0.0 # [m^2/s]
+erf.alpha_C          = 0.0
+
+# INITIAL CONDITIONS
+#erf.init_type = "input_sounding"
+#erf.input_sounding_file = "BF02_moist_sounding"
+#erf.init_sounding_ideal = true
+
+# PROBLEM PARAMETERS (optional)
+# warm bubble input
+prob.x_c    = 10000.0
+prob.z_c    =  2000.0
+prob.x_r    =  2000.0
+prob.z_r    =  2000.0
+prob.T_0    =   300.0
+
+prob.do_moist_bubble = true
+prob.theta_pert  = 2.0
+prob.qt_init     = 0.02
+prob.eq_pot_temp = 320.0

Exec/RegTests/Bubble/prob.cpp

@@ -364,7 +364,8 @@ Problem::init_custom_pert(
 
             if(sc.moisture_type == MoistureType::SAM) {
                 moisture_type = 1;
-            } else if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
+            } else if(sc.moisture_type == MoistureType::SAM_NoIce ||
+                      sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
                 moisture_type = 2;
             }
 

Source/DataStructs/DataStruct.H

@@ -33,7 +33,7 @@ enum struct MoistureModelType{
 };
 
 enum struct MoistureType {
-    Kessler, SAM, SAM_NoPrecip_NoIce, Kessler_NoRain, None
+    Kessler, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler_NoRain, None
 };
 
 enum struct WindFarmType {
@@ -95,6 +95,8 @@ struct SolverChoice {
         pp.query("moisture_model", moisture_model_string);
         if (moisture_model_string == "SAM") {
             moisture_type = MoistureType::SAM;
+        } else if (moisture_model_string == "SAM_NoIce") {
+            moisture_type = MoistureType::SAM_NoIce;
         } else if (moisture_model_string == "SAM_NoPrecip_NoIce") {
             moisture_type = MoistureType::SAM_NoPrecip_NoIce;
         } else if (moisture_model_string == "Kessler") {
@@ -110,6 +112,7 @@ struct SolverChoice {
         if (moisture_type != MoistureType::None) {
             if (moisture_type == MoistureType::Kessler_NoRain ||
                 moisture_type == MoistureType::SAM ||
+                moisture_type == MoistureType::SAM_NoIce ||
                 moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
                 buoyancy_type = 1; // asserted in make buoyancy
             } else {

Source/ERF.cpp

@@ -1451,7 +1451,7 @@ ERF::ReadParameters ()
         lsm.SetModel<NullSurf>();
         Print() << "Null land surface model!\n";
     } else {
-        Abort("Dont know this moisture_type!") ;
+        Abort("Dont know this LandSurfaceType!") ;
     }
 
     if (verbose > 0) {

Source/IO/Plotfile.cpp

@@ -77,11 +77,13 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     for (int i = 0; i < derived_names.size(); ++i) {
         if ( containerHasElement(plot_var_names, derived_names[i]) ) {
             if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) {
-                if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "qi" &&
-                                                                           derived_names[i] != "qsnow" &&
-                                                                           derived_names[i] != "qgraup" &&
-                                                                           derived_names[i] != "snow_accum" &&
-                                                                           derived_names[i] != "graup_accum") )
+                if ( (solverChoice.moisture_type == MoistureType::SAM ||
+                      solverChoice.moisture_type == MoistureType::SAM_NoIce) ||
+                     (derived_names[i] != "qi" &&
+                      derived_names[i] != "qsnow" &&
+                      derived_names[i] != "qgraup" &&
+                      derived_names[i] != "snow_accum" &&
+                      derived_names[i] != "graup_accum") )
                 {
                     tmp_plot_names.push_back(derived_names[i]);
                 } // moisture_type

Source/Microphysics/EulerianMicrophysics.H

@@ -27,11 +27,12 @@ public:
     {
         AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Eulerian );
         m_moist_model.resize(nlev);
-        if (a_model_type == MoistureType::SAM or
+        if (a_model_type == MoistureType::SAM ||
+            a_model_type == MoistureType::SAM_NoIce ||
             a_model_type == MoistureType::SAM_NoPrecip_NoIce) {
             SetModel<SAM>();
             amrex::Print() << "SAM moisture model!\n";
-        } else if (a_model_type == MoistureType::Kessler or
+        } else if (a_model_type == MoistureType::Kessler ||
                    a_model_type == MoistureType::Kessler_NoRain) {
             SetModel<Kessler>();
             amrex::Print() << "Kessler moisture model!\n";

Source/Microphysics/Kessler/Kessler.cpp

@@ -144,12 +144,12 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
 
                     autor = 0.0;
                     if (qcc > qcw0) {
-                        autor = 0.001;
+                        autor = alphaelq;
                     }
 
                     accrr = 0.0;
                     accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875);
-                    dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001));
+                    dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - qcw0));
 
                     // If the amount of change is more than the amount of qc present, then dq = qc
                     dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k));

Source/Microphysics/Microphysics.H

@@ -58,6 +58,7 @@ public:
     static MoistureModelType modelType (const MoistureType a_moisture_type)
     {
         if (    (a_moisture_type == MoistureType::SAM)
+             || (a_moisture_type == MoistureType::SAM_NoIce)
              || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce)
              || (a_moisture_type == MoistureType::Kessler)
              || (a_moisture_type == MoistureType::Kessler_NoRain)

Source/Microphysics/SAM/Cloud_SAM.cpp

@@ -8,7 +8,9 @@ using namespace amrex;
 /**
  * Split cloud components according to saturation pressures; source theta from latent heat.
  */
-void SAM::Cloud (const SolverChoice& solverChoice) {
+void
+SAM::Cloud (const SolverChoice& sc)
+{
 
     constexpr Real an = 1.0/(tbgmax-tbgmin);
     constexpr Real bn = tbgmin*an;
@@ -18,13 +20,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
     Real fac_fus  = m_fac_fus;
     Real rdOcp    = m_rdOcp;
 
-    SolverChoice sc = solverChoice;
-
     int SAM_moisture_type = 1;
-
-    if(solverChoice.moisture_type == MoistureType::SAM) {
-        SAM_moisture_type = 1;
-    } else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
+    if (sc.moisture_type == MoistureType::SAM_NoIce ||
+        sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
         SAM_moisture_type = 2;
     }
 
@@ -128,18 +126,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
                                                   qv_array  , qcl_array  , qci_array,
                                                   qn_array  , qt_array);
 
-                // Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
-                // This would be inconsistent with updating the pressure as part of the iteration above.
-                // Empirically based on the moist bubble rise case, getting the correct theta here
-                // depends on using the old (unchanged by the phase changes) pressure.
+                // Update theta
                 theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
-                // Update pressure to be consistent with updated theta_array
-                pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));
-
-                // Pressure unit conversion
-                pres_array(i,j,k) *= 0.01;
-
             //
             // We cannot blindly relax to qsat, but we can convert qc/qi -> qv.
             // The concept here is that if we put all the moisture into qv and modify
@@ -163,18 +152,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
                 // Update temperature (endothermic since we evap/sublime)
                 tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi;
 
-                // Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
-                // This would be inconsistent with updating the pressure as part of the iteration above.
-                // Empirically based on the moist bubble rise case, getting the correct theta here
-                // depends on using the old (unchanged by the phase changes) pressure.
+                // Update theta
                 theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
-                // Update pressure to be consistent with updated theta_array
-                pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));
-
-                // Pressure unit conversion
-                pres_array(i,j,k) *= 0.01;
-
                 // Verify assumption that qv > qsat does not occur
                 erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
                 erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
@@ -189,18 +169,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
                                                       qv_array  , qcl_array  , qci_array,
                                                       qn_array  , qt_array);
 
-                    // Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
-                    // This would be inconsistent with updating the pressure as part of the iteration above.
-                    // Empirically based on the moist bubble rise case, getting the correct theta here
-                    // depends on using the old (unchanged by the phase changes) pressure.
+                    // Update theta
                     theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
-                    // Update pressure to be consistent with updated theta_array
-                    pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));
-
-                    // Pressure unit conversion
-                    pres_array(i,j,k) *= 0.01;
-
                 }
             }
         });

Source/Microphysics/SAM/IceFall.cpp

@@ -7,7 +7,11 @@ using namespace amrex;
 /**
  * Sedimentation of cloud ice (A32)
  */
-void SAM::IceFall () {
+void SAM::IceFall (const SolverChoice& sc) {
+
+    if(sc.moisture_type == MoistureType::SAM_NoIce ||
+       sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce)
+      return;
 
     Real dz   = m_geom.CellSize(2);
     Real dtn  = dt;

Source/Microphysics/SAM/Init_SAM.cpp

@@ -18,7 +18,8 @@ using namespace amrex;
  * @param[in] geom Geometry associated with these MultiFabs and grids
  * @param[in] dt_advance Timestep for the advance
  */
-void SAM::Init (const MultiFab& cons_in,
+void
+SAM::Init (const MultiFab& cons_in,
                 const BoxArray& grids,
                 const Geometry& geom,
                 const Real& dt_advance,
@@ -83,7 +84,8 @@ void SAM::Init (const MultiFab& cons_in,
  *
  * @param[in] cons_in Conserved variables input
  */
-void SAM::Copy_State_to_Micro (const MultiFab& cons_in)
+void
+SAM::Copy_State_to_Micro (const MultiFab& cons_in)
 {
     // Get the temperature, density, theta, qt and qp from input
     for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {