From 391920e1d18d4f9153c072f10989558150696a7e Mon Sep 17 00:00:00 2001 From: Aaron Lattanzi Date: Wed, 29 May 2024 10:12:11 -0700 Subject: [PATCH] This compared well with Kessler. --- .../Bubble/inputs_BF02_moist_bubble_Kessler | 81 +++++++++++++++++++ .../Bubble/inputs_BF02_moist_bubble_SAM_NoIce | 81 +++++++++++++++++++ Exec/RegTests/Bubble/prob.cpp | 3 +- Source/DataStructs/DataStruct.H | 5 +- Source/ERF.cpp | 2 +- Source/IO/Plotfile.cpp | 12 +-- Source/Microphysics/EulerianMicrophysics.H | 5 +- Source/Microphysics/Kessler/Kessler.cpp | 4 +- Source/Microphysics/Microphysics.H | 1 + Source/Microphysics/SAM/Cloud_SAM.cpp | 47 +++-------- Source/Microphysics/SAM/IceFall.cpp | 6 +- Source/Microphysics/SAM/Init_SAM.cpp | 6 +- Source/Microphysics/SAM/Precip.cpp | 41 ++++++---- Source/Microphysics/SAM/PrecipFall.cpp | 40 +++++++-- Source/Microphysics/SAM/SAM.H | 14 ++-- Source/Microphysics/SAM/Update_SAM.cpp | 3 +- 16 files changed, 268 insertions(+), 83 deletions(-) create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler new file mode 100644 index 000000000..56e38b8fc --- /dev/null +++ b/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 diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce new file mode 100644 index 000000000..7faf745ed --- /dev/null +++ b/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 diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index b6382081b..911f18aa8 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/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; } diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 69a5c1b02..a3351c756 100644 --- a/Source/DataStructs/DataStruct.H +++ b/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 { diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c0f1effb5..68c267430 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1451,7 +1451,7 @@ ERF::ReadParameters () lsm.SetModel(); Print() << "Null land surface model!\n"; } else { - Abort("Dont know this moisture_type!") ; + Abort("Dont know this LandSurfaceType!") ; } if (verbose > 0) { diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index ada718f94..c9e0ff65e 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -77,11 +77,13 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector 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 diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H index 4e1ac1af2..ba5f72735 100644 --- a/Source/Microphysics/EulerianMicrophysics.H +++ b/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(); 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(); amrex::Print() << "Kessler moisture model!\n"; diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 7e8abbf96..f00d8759b 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/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)); diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index ea96f0e79..dd3bc87b4 100644 --- a/Source/Microphysics/Microphysics.H +++ b/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) diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 7783cb5bf..d1c6fcdba 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/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; @@ -20,13 +22,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { Real tol = 1.0e-4; - 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; } @@ -162,7 +160,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) { domn = an; } } else if (SAM_moisture_type == 2) { - omn = 1.0; + omn = 1.0; domn = 0.0; } @@ -184,12 +182,6 @@ void SAM::Cloud (const SolverChoice& solverChoice) { dtabs = -fff/dfff; tabs = tabs+dtabs; - // For now at least we perform this iteration at constant pressure - // For the moist bubble case, the results are indistinguisable - // between running with this used vs commented out - // pres = rho_array(i,j,k) * R_d * tabs - // * (1.0 + R_v/R_d * qsat) * 0.01; - // Update iteration niter = niter+1; } while(std::abs(dtabs) > tol && niter < 20); @@ -212,23 +204,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // Update temperature tabs_array(i,j,k) = tabs; - // 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)); - - // This was used in the earlier implmentation when we updated theta using this new pressure - // pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - // * (1.0 + R_v/R_d * 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 assume the temperature change due to this conversion @@ -251,15 +229,8 @@ 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 pressure - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - - // Update theta from temperature - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - - // Pressure unit conversion - pres_array(i,j,k) *= 0.01; + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); } }); } // mfi diff --git a/Source/Microphysics/SAM/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp index 120b7c725..4951c70f9 100644 --- a/Source/Microphysics/SAM/IceFall.cpp +++ b/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; diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 1d593044f..69abdf31f 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/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) { diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp index 296cc776c..6582c5290 100644 --- a/Source/Microphysics/SAM/Precip.cpp +++ b/Source/Microphysics/SAM/Precip.cpp @@ -6,7 +6,11 @@ using namespace amrex; /** * Autoconversion (A30), Accretion (A28), Evaporation (A24) */ -void SAM::Precip () { +void +SAM::Precip (const SolverChoice& sc) +{ + + if (sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return; Real powr1 = (3.0 + b_rain) / 4.0; Real powr2 = (5.0 + b_rain) / 8.0; @@ -37,6 +41,11 @@ void SAM::Precip () { Real dtn = dt; + int SAM_moisture_type = 1; + if (sc.moisture_type == MoistureType::SAM_NoIce) { + SAM_moisture_type = 2; + } + // get the temperature, dentisy, theta, qt and qp from input for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); @@ -77,9 +86,15 @@ void SAM::Precip () { // Work to be done for autoc/accr or evap if (qn_array(i,j,k)+qp_array(i,j,k) > 0.0) { - omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); - omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); - omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + if (SAM_moisture_type == 2) { + omn = 1.0; + omp = 1.0; + omg = 0.0; + } else { + omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + } qcc = qcl_array(i,j,k); qii = qci_array(i,j,k); @@ -167,12 +182,11 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta (endothermic fusion & exothermic melting) + // Update temperature tabs_array(i,j,k) += fac_fus * ( dqca * (1.0 - omp) - dqia * omp ); - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) *= 0.01; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); } //================================================== @@ -211,12 +225,11 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta (endothermic) + // Update temperature tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg); - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); } } }); diff --git a/Source/Microphysics/SAM/PrecipFall.cpp b/Source/Microphysics/SAM/PrecipFall.cpp index facf45966..d72277cc7 100644 --- a/Source/Microphysics/SAM/PrecipFall.cpp +++ b/Source/Microphysics/SAM/PrecipFall.cpp @@ -11,8 +11,11 @@ using namespace amrex; * * @param[in] hydro_type Type selection for the precipitation advection hydrodynamics scheme (0-3) */ -void SAM::PrecipFall () +void +SAM::PrecipFall (const SolverChoice& sc) { + if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return; + Real rho_0 = 1.29; Real gamr3 = erf_gammafff(4.0+b_rain); @@ -49,6 +52,11 @@ void SAM::PrecipFall () MultiFab fz; fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow); + int SAM_moisture_type = 1; + if (sc.moisture_type == MoistureType::SAM_NoIce) { + SAM_moisture_type = 2; + } + // Add sedimentation of precipitation field to the vert. vel. for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qp_array = qp->array(mfi); @@ -80,8 +88,14 @@ void SAM::PrecipFall () 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 omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + 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; @@ -99,8 +113,14 @@ void SAM::PrecipFall () fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg); if(k==k_lo){ - 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 omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + } rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*1000.0; // Divide by rho_water and convert to mm snow_accum_array(i,j,k) = snow_accum_array(i,j,k) + rho_avg*(1.0-omp)*(1.0-omg)*qp_avg*vrain*dtn/rhos*1000.0; // Divide by rho_snow and convert to mm graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(1.0-omp)*(omg)*qp_avg*vrain*dtn/rhog*1000.0; // Divide by rho_graupel and convert to mm @@ -133,8 +153,14 @@ void SAM::PrecipFall () // Precipitating sedimentation (A19) //================================================== Real dqp = dJinv * (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)); + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + 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)); diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 3b4d9f280..e7fe9706e 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -66,13 +66,13 @@ public: void Cloud (const SolverChoice& sc); // ice physics - void IceFall (); + void IceFall (const SolverChoice& sc); // precip - void Precip (); + void Precip (const SolverChoice& sc); // precip fall - void PrecipFall (); + void PrecipFall (const SolverChoice& sc); // Set up for first time void @@ -126,11 +126,9 @@ public: dt = dt_advance; this->Cloud(sc); - if(sc.moisture_type == MoistureType::SAM) { - this->IceFall(); - this->Precip(); - this->PrecipFall(); - } + this->IceFall(sc); + this->Precip(sc); + this->PrecipFall(sc); } amrex::MultiFab* diff --git a/Source/Microphysics/SAM/Update_SAM.cpp b/Source/Microphysics/SAM/Update_SAM.cpp index 3ddf97a29..ce4bb29ba 100644 --- a/Source/Microphysics/SAM/Update_SAM.cpp +++ b/Source/Microphysics/SAM/Update_SAM.cpp @@ -11,7 +11,8 @@ using namespace amrex; * @param[out] cons Conserved variables * @param[out] qmoist: qv, qc, qi, qr, qs, qg */ -void SAM::Copy_Micro_to_State (MultiFab& cons) +void +SAM::Copy_Micro_to_State (MultiFab& cons) { // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) {