Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

p3 auto, accr, scol runtime params #3037

Merged
merged 5 commits into from
Oct 15, 2024
Merged

p3 auto, accr, scol runtime params #3037

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
938631a
accr/auto runtime
mahf708 Oct 11, 2024
90c7b94
avoid nbfb
mahf708 Oct 11, 2024
7c859f0
add selfcollection prefactor
mahf708 Oct 11, 2024
540f8e0
add autoconversion radius to namelist
mahf708 Oct 11, 2024
fed60e8
fixes nbfb, and better document breakup thresh
mahf708 Oct 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions components/eamxx/cime_config/namelist_defaults_scream.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,18 +202,24 @@ be lost if SCREAM_HACK_XML is not enabled.
${DIN_LOC_ROOT}/atm/scream/tables/vn_table_vals.dat8,
${DIN_LOC_ROOT}/atm/scream/tables/vm_table_vals.dat8
</tables>
<p3_autoconversion_prefactor type="real" doc="P3 autoconversion_prefactor (scale factor in autoconversion)">1350.0</p3_autoconversion_prefactor>
<p3_autoconversion_prefactor type="real" doc="Autoconversion linear prefactor in P3">1350.0</p3_autoconversion_prefactor>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Style question: do we need to prepend all param names with p3_, considering we are already in the p3 sublist?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't like the p3_ but I hesitated removing them. Your opinion will nudge me to remove this. Should I remove them? If so, I will do that in the next PR that will be NBFB to address the two TODO items.

<p3_autoconversion_qc_exponent type="real" doc="Autoconversion qc exponent in P3">2.47</p3_autoconversion_qc_exponent>
<p3_autoconversion_nc_exponent type="real" doc="Autoconversion nc exponent in P3">1.79</p3_autoconversion_nc_exponent>
<p3_autoconversion_radius type="real" doc="Autoconversion droplet radius (meter) in P3">25.0e-6</p3_autoconversion_radius>
<p3_accretion_prefactor type="real" doc="Accretion linear prefactor in P3">67.0</p3_accretion_prefactor>
<p3_accretion_qc_exponent type="real" doc="Accretion qc exponent in P3">1.15</p3_accretion_qc_exponent>
<p3_accretion_qr_exponent type="real" doc="Accretion qr exponent in P3">1.15</p3_accretion_qr_exponent>
<p3_rain_selfcollection_prefactor type="real" doc="Rain selfcollection prefactor in P3">5.78</p3_rain_selfcollection_prefactor>
<p3_rain_selfcollection_breakup_diameter type="real" doc="Rain selfcollection breakup diameter (meter) in P3">0.00028</p3_rain_selfcollection_breakup_diameter>
<p3_mu_r_constant type="real" doc="P3 mu_r_constant (rain shape parameter in gamma drop-size distribution)">1.0</p3_mu_r_constant>
<p3_spa_to_nc type="real" doc="P3 spa_to_nc (scaling factor for turning CCN into nc in SPA)">1.0</p3_spa_to_nc>
<p3_k_accretion type="real" doc="P3 k_accretion (scaling factor on accretion)">67.0</p3_k_accretion>
<p3_eci type="real" doc="P3 eci (liquid/ice collision/collection coefficient)">0.5</p3_eci>
<p3_eri type="real" doc="P3 eri (ice/rain collision/collection coefficient)">1.0</p3_eri>
<p3_rho_rime_min type="real" doc="P3 rho_rime_min (riming density maximum)">50.0</p3_rho_rime_min>
<p3_rho_rime_max type="real" doc="P3 rho_rime_max (riming density minimum)">900.0</p3_rho_rime_max>
<p3_a_imm type="real" doc="P3 a_imm (drop freezing exponent )">0.65</p3_a_imm>
<p3_dep_nucleation_exponent type="real" doc="P3 dep_nucleation_exponent (deposition nucleation)">0.304</p3_dep_nucleation_exponent>
<p3_ice_sed_knob type="real" doc="P3 ice_sed_knob (ice fall speed)">1.0</p3_ice_sed_knob>
<p3_d_breakup_cutoff type="real" doc="P3 d_breakup_cutoff (rain self collection and breakup)">0.00028</p3_d_breakup_cutoff>
<p3_do_ice_production type="logical" doc="Flag to turn on ice production processes in P3 (loss processes unaffected by this flag)">true</p3_do_ice_production>
</p3>

Expand Down
10 changes: 8 additions & 2 deletions components/eamxx/src/physics/p3/eamxx_p3_process_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,17 +222,23 @@ void P3Microphysics::initialize_impl (const RunType /* run_type */)
// Gather runtime options
runtime_options.max_total_ni = m_params.get<double>("max_total_ni", 740.0e3);
runtime_options.p3_autoconversion_prefactor = m_params.get<double>("p3_autoconversion_prefactor", 1350.0);
runtime_options.p3_autoconversion_qc_exponent = m_params.get<double>("p3_autoconversion_qc_exponent", 2.47);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another style question: should we maybe move all the params loading into a P3Runtime::load_params kind of method? May help keep the process interface a tad leaner.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd like to do this in the next PR as well. Thanks!

runtime_options.p3_autoconversion_nc_exponent = m_params.get<double>("p3_autoconversion_nc_exponent", 1.79);
runtime_options.p3_autoconversion_radius = m_params.get<double>("p3_autoconversion_radius", 25.0e-6);
runtime_options.p3_accretion_prefactor = m_params.get<double>("p3_accretion_prefactor", 67.0);
runtime_options.p3_accretion_qc_exponent = m_params.get<double>("p3_accretion_qc_exponent", 1.15);
runtime_options.p3_accretion_qr_exponent = m_params.get<double>("p3_accretion_qr_exponent", 1.15);
runtime_options.p3_rain_selfcollection_prefactor = m_params.get<double>("p3_rain_selfcollection_prefactor", 5.78);
runtime_options.p3_rain_selfcollection_breakup_diameter = m_params.get<double>("p3_rain_selfcollection_breakup_diameter", 0.00028);
runtime_options.p3_mu_r_constant = m_params.get<double>("p3_mu_r_constant", 1.0);
runtime_options.p3_spa_to_nc = m_params.get<double>("p3_spa_to_nc", 1.0);
runtime_options.p3_k_accretion = m_params.get<double>("p3_k_accretion", 67.0);
runtime_options.p3_eci = m_params.get<double>("p3_eci", 0.5);
runtime_options.p3_eri = m_params.get<double>("p3_eri", 1.0);
runtime_options.p3_rho_rime_min = m_params.get<double>("p3_rho_rime_min", 50.0);
runtime_options.p3_rho_rime_max = m_params.get<double>("p3_rho_rime_max", 900.0);
runtime_options.p3_a_imm = m_params.get<double>("p3_a_imm", 0.65);
runtime_options.p3_dep_nucleation_exponent = m_params.get<double>("p3_dep_nucleation_exponent", 0.304);
runtime_options.p3_ice_sed_knob = m_params.get<double>("p3_ice_sed_knob", 1.0);
runtime_options.p3_d_breakup_cutoff = m_params.get<double>("p3_d_breakup_cutoff", 0.00028);
runtime_options.p3_do_ice_production = m_params.get<bool>("p3_do_ice_production", true);

// Set property checks for fields in this process
Expand Down
21 changes: 15 additions & 6 deletions components/eamxx/src/physics/p3/impl/p3_autoconversion_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,20 +19,29 @@ ::cloud_water_autoconversion(

// Khroutdinov and Kogan (2000)
const auto qc_not_small = qc_incld >= 1e-8 && context;
constexpr Scalar CONS3 = C::CONS3;

const Scalar p3_autoconversion_prefactor = runtime_options.p3_autoconversion_prefactor;
const Scalar p3_autoconversion_qc_exponent = runtime_options.p3_autoconversion_qc_exponent;
const Scalar p3_autoconversion_nc_exponent = runtime_options.p3_autoconversion_nc_exponent;
const Scalar p3_autoconversion_radius = runtime_options.p3_autoconversion_radius;

if(qc_not_small.any()){
// TODO: correct this later (by keeping commented-out def) once BFB reqs are satisfied
const Scalar CONS3 = C::CONS3; // sp(1.0) / (C::CONS2 * pow(p3_autoconversion_radius, sp(3.0)));

if(qc_not_small.any()) {
Spack sgs_var_coef;
// sgs_var_coef = subgrid_variance_scaling(inv_qc_relvar, sp(2.47) );
sgs_var_coef = 1;

qc2qr_autoconv_tend.set(qc_not_small,
sgs_var_coef*p3_autoconversion_prefactor*pow(qc_incld,sp(2.47))*pow(nc_incld*sp(1.e-6)*rho,sp(-1.79)));
qc2qr_autoconv_tend.set(
qc_not_small,
sgs_var_coef * p3_autoconversion_prefactor *
pow(qc_incld, p3_autoconversion_qc_exponent) *
pow(nc_incld * sp(1.e-6) * rho, -p3_autoconversion_nc_exponent));
// note: ncautr is change in Nr; nc2nr_autoconv_tend is change in Nc
ncautr.set(qc_not_small, qc2qr_autoconv_tend*CONS3);
nc2nr_autoconv_tend.set(qc_not_small, qc2qr_autoconv_tend*nc_incld/qc_incld);
ncautr.set(qc_not_small, qc2qr_autoconv_tend * CONS3);
nc2nr_autoconv_tend.set(qc_not_small,
qc2qr_autoconv_tend * nc_incld / qc_incld);
}

nc2nr_autoconv_tend.set(qc2qr_autoconv_tend == 0 && context, 0);
Expand Down
17 changes: 12 additions & 5 deletions components/eamxx/src/physics/p3/impl/p3_cloud_rain_acc_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,18 +25,25 @@ ::cloud_rain_accretion(
{
constexpr Scalar qsmall = C::QSMALL;

const Scalar p3_k_accretion = runtime_options.p3_k_accretion;;
const Scalar p3_accretion_prefactor = runtime_options.p3_accretion_prefactor;
const Scalar p3_accretion_qc_exponent = runtime_options.p3_accretion_qc_exponent;
const Scalar p3_accretion_qr_exponent = runtime_options.p3_accretion_qr_exponent;

Spack sgs_var_coef;
// sgs_var_coef = subgrid_variance_scaling(inv_qc_relvar, sp(1.15) );
sgs_var_coef = 1;

const auto qr_and_qc_not_small = (qr_incld >= qsmall) && (qc_incld >= qsmall) && context;
if (qr_and_qc_not_small.any()) {
if(qr_and_qc_not_small.any()) {
// Khroutdinov and Kogan (2000)
qc2qr_accret_tend.set(qr_and_qc_not_small,
sgs_var_coef * sp(p3_k_accretion) * pow(qc_incld * qr_incld, sp(1.15)));
nc_accret_tend.set(qr_and_qc_not_small, qc2qr_accret_tend * nc_incld / qc_incld);
qc2qr_accret_tend.set(
qr_and_qc_not_small,
sgs_var_coef * sp(p3_accretion_prefactor) *
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's no need to use sp for accretion prefactor. sp should only be used with actual numeric constants, to avoid implicit promotion to double inside the calculation. Here, p3_accretion_prefactor already has type Scalar.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Likewise, I will fix this in the PR making switching things on

// TODO: split this pow into two in a later PR due to NBFB behavior
// pow(qc_incld, p3_accretion_qc_exponent) * pow(qr_incld, p3_accretion_qr_exponent));
pow(qc_incld * qr_incld, p3_accretion_qr_exponent));
nc_accret_tend.set(qr_and_qc_not_small,
qc2qr_accret_tend * nc_incld / qc_incld);

qc2qr_accret_tend.set(nc_accret_tend == 0 && context, 0);
nc_accret_tend.set(qc2qr_accret_tend == 0 && context, 0);
Expand Down
28 changes: 19 additions & 9 deletions components/eamxx/src/physics/p3/impl/p3_rain_self_collection_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,23 +22,33 @@ ::rain_self_collection(
constexpr Scalar rho_h2o = C::RHO_H2O;
constexpr Scalar pi = C::Pi;

const Scalar p3_d_breakup_cutoff = runtime_options.p3_d_breakup_cutoff;
const Scalar p3_rain_selfcollection_breakup_diameter = runtime_options.p3_rain_selfcollection_breakup_diameter;
const Scalar p3_rain_selfcollection_prefactor = runtime_options.p3_rain_selfcollection_prefactor;

const auto qr_incld_not_small = qr_incld >= qsmall && context;

if (qr_incld_not_small.any()) {

const auto dum2 = cbrt((qr_incld)/(pi*rho_h2o*nr_incld));
if(qr_incld_not_small.any()) {
// use mass-mean diameter (do this by using
// the old version of lambda w/o mu dependence)
// note there should be a factor of 6^(1/3), but we
// want to keep breakup threshold consistent so 'dum'
// is expressed in terms of lambda rather than mass-mean D
const auto dum2 = cbrt((qr_incld) / (pi * rho_h2o * nr_incld));

Spack dum;
const auto dum2_lt_dum1 = dum2 < p3_d_breakup_cutoff && qr_incld_not_small;
const auto dum2_gt_dum1 = dum2 >= p3_d_breakup_cutoff && qr_incld_not_small;
const auto dum2_lt_dum1 =
dum2 < p3_rain_selfcollection_breakup_diameter && qr_incld_not_small;
const auto dum2_gt_dum1 =
dum2 >= p3_rain_selfcollection_breakup_diameter && qr_incld_not_small;
dum.set(dum2_lt_dum1, 1);
if (dum2_gt_dum1.any()) {
dum.set(dum2_gt_dum1, 2 - exp(2300 * (dum2-p3_d_breakup_cutoff)));
if(dum2_gt_dum1.any()) {
dum.set(dum2_gt_dum1,
2 - exp(2300 * (dum2 - p3_rain_selfcollection_breakup_diameter)));
}

nr_selfcollect_tend.set(qr_incld_not_small, dum*sp(5.78)*nr_incld*qr_incld*rho);
nr_selfcollect_tend.set(
qr_incld_not_small,
dum * p3_rain_selfcollection_prefactor * nr_incld * qr_incld * rho);
}
}

Expand Down
10 changes: 8 additions & 2 deletions components/eamxx/src/physics/p3/p3_functions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,17 +111,23 @@ struct Functions
struct P3Runtime {
Scalar max_total_ni = 740.0e3;
Scalar p3_autoconversion_prefactor = 1350.0;
Scalar p3_autoconversion_qc_exponent = 2.47;
Scalar p3_autoconversion_nc_exponent = 1.79;
Scalar p3_autoconversion_radius = 25.0e-6;
Scalar p3_accretion_prefactor = 67.0;
Scalar p3_accretion_qc_exponent = 1.15;
Scalar p3_accretion_qr_exponent = 1.15;
Scalar p3_rain_selfcollection_prefactor = 5.78;
Scalar p3_rain_selfcollection_breakup_diameter = 0.00028;
Scalar p3_mu_r_constant = 1.0;
Scalar p3_spa_to_nc = 1.0;
Scalar p3_k_accretion = 67.0;
Scalar p3_eci = 0.5;
Scalar p3_eri = 1.0;
Scalar p3_rho_rime_min = 50.0;
Scalar p3_rho_rime_max = 900.0;
Scalar p3_a_imm = 0.65;
Scalar p3_dep_nucleation_exponent = 0.304;
Scalar p3_ice_sed_knob = 1.0;
Scalar p3_d_breakup_cutoff = 0.00028;
bool p3_do_ice_production = true;
};

Expand Down