Skip to content

add runtime options for p3 #3033

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

Closed
wants to merge 6 commits into from
Closed

add runtime options for p3 #3033

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
bba279b
Added nl parameters for autoconversion and accretion exponents as wel…
hassanbeydoun Oct 2, 2024
2fc1726
Merge branch 'master' into cumulus_autoconv
hassanbeydoun Oct 2, 2024
e3b223c
fixed bug in assignment
hassanbeydoun Oct 2, 2024
7895945
added nl parameter for radius of newly formed rain drops from autocon…
hassanbeydoun Oct 4, 2024
699d103
formatting and comments
mahf708 Oct 8, 2024
9749a30
try to maintain bfb with older code
mahf708 Oct 8, 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
6 changes: 6 additions & 0 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 @@ -203,6 +203,11 @@ be lost if SCREAM_HACK_XML is not enabled.
${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_qc_exp type="real" doc="P3 autoconversion qc exponent">2.47</p3_autoconversion_qc_exp>
<p3_autoconversion_nc_exp type="real" doc="P3 autoconversion nc exponent">1.79</p3_autoconversion_nc_exp>
<p3_autoconversion_radius type="real" doc="P3 autoconversion radius of new rain drops">0.000025</p3_autoconversion_radius>
Copy link
Contributor

Choose a reason for hiding this comment

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

It may be a good idea to add the units in the doc string.

Copy link
Contributor Author

@mahf708 mahf708 Oct 8, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

I wanna edit the docs for all p3 constants tbh...

<p3_accretion_qc_exp type="real" doc="P3 accretion qc exponent">1.15</p3_accretion_qc_exp>
<p3_accretion_qr_exp type="real" doc="P3 accretion qr exponent">1.15</p3_accretion_qr_exp>
<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>
Expand All @@ -214,6 +219,7 @@ be lost if SCREAM_HACK_XML is not enabled.
<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 type="logical" doc="P3 does ice process or not">true</p3_do_ice>
</p3>

<!-- SHOC macrophysics -->
Expand Down
28 changes: 22 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 @@ -20,22 +20,38 @@ ::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 = p3constants.p3_autoconversion_prefactor;
const Scalar p3_autoconversion_qc_exp = p3constants.p3_autoconversion_qc_exp;
const Scalar p3_autoconversion_nc_exp = p3constants.p3_autoconversion_nc_exp;
const Scalar p3_autoconversion_radius = p3constants.p3_autoconversion_radius;

//printf(" hey inside AAAAAAAAAAAAAAAA %13.6f \n", p3_autoconversion_factor);
Scalar scale_nc2nr_tend;
constexpr Scalar CONS2 = C::CONS2;
constexpr Scalar CONS3 = C::CONS3;
if (p3_autoconversion_radius == 0.000025) {
scale_nc2nr_tend = CONS3;
} else {
scale_nc2nr_tend = sp(1.)/(CONS2 * pow(p3_autoconversion_radius, 3));
}

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, sp(p3_autoconversion_qc_exp)) *
pow(nc_incld * sp(1.e-6) * rho, sp(-p3_autoconversion_nc_exp)));
// 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 * scale_nc2nr_tend);
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
9 changes: 7 additions & 2 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 @@ -26,6 +26,8 @@ ::cloud_rain_accretion(
constexpr Scalar qsmall = C::QSMALL;

const Scalar p3_k_accretion = p3constants.p3_k_accretion;;
const Scalar p3_accretion_qc_exp = p3constants.p3_accretion_qc_exp;
const Scalar p3_accretion_qr_exp = p3constants.p3_accretion_qr_exp;

Spack sgs_var_coef;
// sgs_var_coef = subgrid_variance_scaling(inv_qc_relvar, sp(1.15) );
Expand All @@ -35,8 +37,11 @@ ::cloud_rain_accretion(
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);
sgs_var_coef * sp(p3_k_accretion) *
pow(qc_incld, sp(p3_accretion_qc_exp)) *
pow(qr_incld, sp(p3_accretion_qr_exp)));
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
10 changes: 6 additions & 4 deletions components/eamxx/src/physics/p3/impl/p3_main_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,8 @@ ::p3_main_internal(
// we do not want to measure init stuff
auto start = std::chrono::steady_clock::now();

const bool p3_do_ice = p3constants.p3_do_ice;

// p3_main loop
Kokkos::parallel_for(
"p3 main loop",
Expand Down Expand Up @@ -288,10 +290,10 @@ ::p3_main_internal(
lookup_tables.ice_table_vals, diagnostic_outputs.precip_ice_surf(i), p3constants);

// homogeneous freezing of cloud and rain
homogeneous_freezing(
T_atm, oinv_exner, team, nk, ktop, kbot, kdir, oqc, onc, oqr, onr, oqi,
oni, oqm, obm, oth);

if(p3_do_ice) {
homogeneous_freezing(T_atm, oinv_exner, team, nk, ktop, kbot, kdir, oqc,
onc, oqr, onr, oqi, oni, oqm, obm, oth);
}
//
// final checks to ensure consistency of mass/number
// and compute diagnostic fields for output
Expand Down
127 changes: 76 additions & 51 deletions components/eamxx/src/physics/p3/impl/p3_main_impl_part2.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,8 @@ ::p3_main_part2(
constexpr Scalar latvap = C::LatVap;
constexpr Scalar latice = C::LatIce;

const bool p3_do_ice = p3constants.p3_do_ice;

team.team_barrier();
hydrometeorsPresent = false;
team.team_barrier();
Expand Down Expand Up @@ -267,60 +269,78 @@ ::p3_main_part2(
// ......................................................................

// collection of droplets
ice_cldliq_collection(
rho(k), T_atm(k), rhofaci(k), table_val_qc2qi_collect, qi_incld(k), qc_incld(k), ni_incld(k), nc_incld(k),
qc2qi_collect_tend, nc_collect_tend, qc2qr_ice_shed_tend, ncshdc, p3constants, not_skip_micro);

if(p3_do_ice) {
ice_cldliq_collection(rho(k), T_atm(k), rhofaci(k),
table_val_qc2qi_collect, qi_incld(k), qc_incld(k),
ni_incld(k), nc_incld(k), qc2qi_collect_tend,
nc_collect_tend, qc2qr_ice_shed_tend, ncshdc,
p3constants, not_skip_micro);
}
// collection of rain
ice_rain_collection(
rho(k), T_atm(k), rhofaci(k), logn0r(k), table_val_nr_collect, table_val_qr2qi_collect, qi_incld(k), ni_incld(k), qr_incld(k),
qr2qi_collect_tend, nr_collect_tend, p3constants, not_skip_micro);

if(p3_do_ice) {
ice_rain_collection(
rho(k), T_atm(k), rhofaci(k), logn0r(k), table_val_nr_collect,
table_val_qr2qi_collect, qi_incld(k), ni_incld(k), qr_incld(k),
qr2qi_collect_tend, nr_collect_tend, p3constants, not_skip_micro);
}
// collection between ice categories

// PMC nCat deleted lots of stuff here.

// self-collection of ice
ice_self_collection(
rho(k), rhofaci(k), table_val_ni_self_collect, eii, qm_incld(k), qi_incld(k), ni_incld(k),
ni_selfcollect_tend, not_skip_micro);

if(p3_do_ice) {
ice_self_collection(rho(k), rhofaci(k), table_val_ni_self_collect, eii,
qm_incld(k), qi_incld(k), ni_incld(k),
ni_selfcollect_tend, not_skip_micro);
}
// melting
ice_melting(
rho(k), T_atm(k), pres(k), rhofaci(k), table_val_qi2qr_melting, table_val_qi2qr_vent_melt, dv, sc, mu, kap, qv(k), qi_incld(k), ni_incld(k),
qi2qr_melt_tend, ni2nr_melt_tend, not_skip_micro);

if(p3_do_ice) {
ice_melting(rho(k), T_atm(k), pres(k), rhofaci(k),
table_val_qi2qr_melting, table_val_qi2qr_vent_melt, dv, sc,
mu, kap, qv(k), qi_incld(k), ni_incld(k), qi2qr_melt_tend,
ni2nr_melt_tend, not_skip_micro);
}
// calculate wet growth
ice_cldliq_wet_growth(
rho(k), T_atm(k), pres(k), rhofaci(k), table_val_qi2qr_melting, table_val_qi2qr_vent_melt,
dv, kap, mu, sc, qv(k), qc_incld(k), qi_incld(k), ni_incld(k), qr_incld(k),
wetgrowth, qr2qi_collect_tend, qc2qi_collect_tend, qc_growth_rate, nr_ice_shed_tend, qc2qr_ice_shed_tend, not_skip_micro);

// calculate total inverse ice relaxation timescale combined for all ice categories
// note 'f1pr' values are normalized, so we need to multiply by N
ice_relaxation_timescale(
rho(k), T_atm(k), rhofaci(k), table_val_qi2qr_melting, table_val_qi2qr_vent_melt, dv, mu, sc, qi_incld(k), ni_incld(k),
epsi, epsi_tot, not_skip_micro);

if(p3_do_ice) {
ice_cldliq_wet_growth(
rho(k), T_atm(k), pres(k), rhofaci(k), table_val_qi2qr_melting,
table_val_qi2qr_vent_melt, dv, kap, mu, sc, qv(k), qc_incld(k),
qi_incld(k), ni_incld(k), qr_incld(k), wetgrowth,
qr2qi_collect_tend, qc2qi_collect_tend, qc_growth_rate,
nr_ice_shed_tend, qc2qr_ice_shed_tend, not_skip_micro);
}
// calculate total inverse ice relaxation timescale combined for all ice
// categories note 'f1pr' values are normalized, so we need to multiply by N
if(p3_do_ice) {
ice_relaxation_timescale(
rho(k), T_atm(k), rhofaci(k), table_val_qi2qr_melting,
table_val_qi2qr_vent_melt, dv, mu, sc, qi_incld(k), ni_incld(k),
epsi, epsi_tot, not_skip_micro);
}
// calculate rime density
calc_rime_density(
T_atm(k), rhofaci(k), table_val_qi_fallspd, acn(k), lamc(k), mu_c(k), qc_incld(k), qc2qi_collect_tend,
vtrmi1, rho_qm_cloud, not_skip_micro);

if(p3_do_ice) {
calc_rime_density(T_atm(k), rhofaci(k), table_val_qi_fallspd, acn(k),
lamc(k), mu_c(k), qc_incld(k), qc2qi_collect_tend,
vtrmi1, rho_qm_cloud, not_skip_micro);
}
// contact and immersion freezing droplets
cldliq_immersion_freezing(
T_atm(k), lamc(k), mu_c(k), cdist1(k), qc_incld(k), inv_qc_relvar(k),
qc2qi_hetero_freeze_tend, nc2ni_immers_freeze_tend, p3constants, not_skip_micro);

if(p3_do_ice) {
cldliq_immersion_freezing(
T_atm(k), lamc(k), mu_c(k), cdist1(k), qc_incld(k),
inv_qc_relvar(k), qc2qi_hetero_freeze_tend,
nc2ni_immers_freeze_tend, p3constants, not_skip_micro);
}
// for future: get rid of log statements below for rain freezing
rain_immersion_freezing(
T_atm(k), lamr(k), mu_r(k), cdistr(k), qr_incld(k),
qr2qi_immers_freeze_tend, nr2ni_immers_freeze_tend, p3constants, not_skip_micro);

if(p3_do_ice) {
rain_immersion_freezing(T_atm(k), lamr(k), mu_r(k), cdistr(k),
qr_incld(k), qr2qi_immers_freeze_tend,
nr2ni_immers_freeze_tend, p3constants,
not_skip_micro);
}
// rime splintering (Hallet-Mossop 1974)
// PMC comment: Morrison and Milbrandt 2015 part 1 and 2016 part 3 both say
// that Hallet-Mossop should be neglected if 1 category to compensate for
// artificial smearing out of ice DSD
// PMC comment: Morrison and Milbrandt 2015 part 1 and 2016 part 3 both
// say that Hallet-Mossop should be neglected if 1 category to compensate
// for artificial smearing out of ice DSD

// ................................................
// condensation/evaporation/deposition/sublimation
Expand All @@ -335,19 +355,24 @@ ::p3_main_part2(
cld_frac_l(k),cld_frac_r(k),qv(k),qv_prev(k),qv_sat_l(k),qv_sat_i(k),
ab,abi,epsr,epsi_tot,T_atm(k),t_prev(k),dqsdt,dt,
qr2qv_evap_tend,nr_evap_tend, not_skip_micro);

ice_deposition_sublimation(
qi_incld(k), ni_incld(k), T_atm(k), qv_sat_l(k), qv_sat_i(k), epsi, abi, qv(k), inv_dt,
qv2qi_vapdep_tend, qi2qv_sublim_tend, ni_sublim_tend, qc2qi_berg_tend, not_skip_micro);
if(p3_do_ice) {
ice_deposition_sublimation(
qi_incld(k), ni_incld(k), T_atm(k), qv_sat_l(k), qv_sat_i(k), epsi,
abi, qv(k), inv_dt, qv2qi_vapdep_tend, qi2qv_sublim_tend,
ni_sublim_tend, qc2qi_berg_tend, not_skip_micro);
}
}

// deposition/condensation-freezing nucleation
ice_nucleation(
T_atm(k), inv_rho(k), ni(k), ni_activated(k), qv_supersat_i(k), inv_dt, predictNc, do_prescribed_CCN,
qv2qi_nucleat_tend, ni_nucleat_tend, p3constants, not_skip_all);

if(p3_do_ice) {
ice_nucleation(T_atm(k), inv_rho(k), ni(k), ni_activated(k),
qv_supersat_i(k), inv_dt, predictNc, do_prescribed_CCN,
qv2qi_nucleat_tend, ni_nucleat_tend, p3constants,
not_skip_all);
}
// cloud water autoconversion
// NOTE: cloud_water_autoconversion must be called before droplet_self_collection
// NOTE: cloud_water_autoconversion must be called before
// droplet_self_collection
cloud_water_autoconversion(
rho(k), qc_incld(k), nc_incld(k), inv_qc_relvar(k),
qc2qr_autoconv_tend, nc2nr_autoconv_tend, ncautr, p3constants, not_skip_all);
Expand Down
48 changes: 48 additions & 0 deletions components/eamxx/src/physics/share/physics_constants.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,11 @@ struct P3_Constants
{
public:
Scalar p3_autoconversion_prefactor = 1350.0;
Scalar p3_autoconversion_qc_exp = 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.

Random question: do we know why the struct P3_Constants is in physics/share rathe than physics/p3? I don't see it used anywhere else. Do we foresee other parts of eamxx to interact with it?

Copy link
Contributor

Choose a reason for hiding this comment

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

I agree, I think we should pull it out of physics/share. I fear that in the future we may see physics/share/physics_constants.hpp as a dumping ground for a bunch of different processes that may or may not be in use.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This was done in #2655. I agree with both of you on not having these here, but I view that as a separate PR

Copy link
Contributor

Choose a reason for hiding this comment

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

Thanks Naser for volunteering to do it an a separate PR... :P

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If you look inside components/eamxx/src/physics/share/physics_constants.hpp, I think it is already the case that it is a dumping ground ...

Copy link
Contributor

Choose a reason for hiding this comment

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

(jk)

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 can, but it will have to wait a bit :D

Scalar p3_autoconversion_nc_exp = 1.79;
Scalar p3_autoconversion_radius = 0.000025;
Scalar p3_accretion_qc_exp = 1.15;
Scalar p3_accretion_qr_exp = 1.15;
Scalar p3_mu_r_constant = 1.0;
Scalar p3_spa_to_nc = 1.0;
Scalar p3_k_accretion = 67.0;
Expand All @@ -138,13 +143,34 @@ struct P3_Constants
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 = true;

void set_p3_from_namelist(ekat::ParameterList &params){

std::string nname = "p3_autoconversion_prefactor";
if(params.isParameter(nname))
p3_autoconversion_prefactor = params.get<double>(nname);

nname = "p3_autoconversion_qc_exp";
if(params.isParameter(nname))
p3_autoconversion_qc_exp = params.get<double>(nname);

nname = "p3_autoconversion_nc_exp";
if(params.isParameter(nname))
p3_autoconversion_nc_exp = params.get<double>(nname);

nname = "p3_autoconversion_radius";
if(params.isParameter(nname))
p3_autoconversion_radius = params.get<double>(nname);

nname = "p3_accretion_qc_exp";
if(params.isParameter(nname))
p3_accretion_qc_exp = params.get<double>(nname);

nname = "p3_accretion_qr_exp";
if(params.isParameter(nname))
p3_accretion_qr_exp = params.get<double>(nname);

nname = "p3_mu_r_constant";
if(params.isParameter(nname))
p3_mu_r_constant = params.get<double>(nname);
Expand Down Expand Up @@ -189,13 +215,32 @@ struct P3_Constants
if(params.isParameter(nname))
p3_d_breakup_cutoff = params.get<double>(nname);

nname = "p3_do_ice";
if(params.isParameter(nname))
p3_do_ice = params.get<bool>(nname);

};

void print_p3constants(std::shared_ptr<ekat::logger::LoggerBase> logger){
logger->info("P3 Constants:");

std::string nname = "p3_autoconversion_prefactor";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_prefactor));

nname = "p3_autoconversion_qc_exp";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_qc_exp));

nname = "p3_autoconversion_nc_exp";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_nc_exp));

nname = "p3_autoconversion_radius";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_radius));

nname = "p3_accretion_qc_exp";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_qc_exp));

nname = "p3_accretion_qr_exp";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_autoconversion_nc_exp));

nname = "p3_mu_r_constant";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_mu_r_constant));
Expand Down Expand Up @@ -230,6 +275,9 @@ struct P3_Constants
nname = "p3_d_breakup_cutoff";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_d_breakup_cutoff));

nname = "p3_do_ice";
logger->info(std::string("P3 ") + nname + std::string(" = ") + std::to_string(p3_do_ice));

logger->info(" ");
};

Expand Down