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

EAMxx: Adds aerosols heterogeneous freezing calculations in P3 microphysics #6947

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

EAMxx: Adds aerosols heterogeneous freezing calculations in P3 microphysics #6947

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
ede97b6
Adds heterogeneous freezing from aerosols in P3 microphysics
singhbalwinder Nov 28, 2024
9be599e
Passes all tests on Compy
singhbalwinder Feb 3, 2025
3a76876
Turn of this feature by default, all P3 tests pass on PM-GPUs
singhbalwinder Feb 4, 2025
dcf1f5f
Fixes p3_test by adding new args
singhbalwinder Feb 5, 2025
f98ac90
Removes debug statements
singhbalwinder Feb 5, 2025
85b6e6e
Partially remove use_hetfrz_classnuc from arg lists and use runtime_…
singhbalwinder Feb 5, 2025
dee8f35
use_hetfrz_classnuc is not controlled only by runtime_options
singhbalwinder Feb 5, 2025
5808e72
Modified CNT function name to more readable name
singhbalwinder Feb 5, 2025
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
1 change: 1 addition & 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 @@ -192,6 +192,7 @@ be lost if SCREAM_HACK_XML is not enabled.
<!-- P3 microphysics -->
<p3 inherit="atm_proc_base">
<do_prescribed_ccn>true</do_prescribed_ccn>
<use_hetfrz_classnuc type="logical" doc="Switch to turn on heterogeneous freezing due to prognostic aerosols">false</use_hetfrz_classnuc>
<do_prescribed_ccn COMPSET=".*SCREAM.*noAero">false</do_prescribed_ccn>
<do_predict_nc>true</do_predict_nc>
<do_predict_nc COMPSET=".*SCREAM.*noAero">false</do_predict_nc>
Expand Down
1 change: 1 addition & 0 deletions components/eamxx/src/physics/p3/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ if (NOT EAMXX_ENABLE_GPU OR Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE OR Kokkos
eti/p3_back_to_cell_average.cpp
eti/p3_cloud_rain_acc.cpp
eti/p3_calc_rime_density.cpp
eti/p3_CNT_couple_impl.cpp
eti/p3_cldliq_imm_freezing.cpp
eti/p3_rain_imm_freezing.cpp
eti/p3_droplet_self_coll.cpp
Expand Down
6 changes: 6 additions & 0 deletions components/eamxx/src/physics/p3/disp/p3_main_impl_disp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,11 @@ ::p3_main_internal_disp(
auto nevapr = diagnostic_outputs.nevapr;
auto qv_prev = diagnostic_inputs.qv_prev;
auto t_prev = diagnostic_inputs.t_prev;
// Inputs for the heteogeneous freezing
auto hetfrz_immersion_nucleation_tend = diagnostic_inputs.hetfrz_immersion_nucleation_tend;
auto hetfrz_contact_nucleation_tend = diagnostic_inputs.hetfrz_contact_nucleation_tend;
auto hetfrz_deposition_nucleation_tend = diagnostic_inputs.hetfrz_deposition_nucleation_tend;

auto liq_ice_exchange = history_only.liq_ice_exchange;
auto vap_liq_exchange = history_only.vap_liq_exchange;
auto vap_ice_exchange = history_only.vap_ice_exchange;
Expand Down Expand Up @@ -252,6 +257,7 @@ ::p3_main_internal_disp(

p3_main_part2_disp(
nj, nk, runtime_options.max_total_ni, infrastructure.predictNc, infrastructure.prescribedCCN, infrastructure.dt, inv_dt,
runtime_options.use_hetfrz_classnuc, hetfrz_immersion_nucleation_tend, hetfrz_contact_nucleation_tend, hetfrz_deposition_nucleation_tend,
lookup_tables.dnu_table_vals, lookup_tables.ice_table_vals, lookup_tables.collect_table_vals,
lookup_tables.revap_table_vals, pres, dpres, dz, nc_nuceat_tend, inv_exner,
exner, inv_cld_frac_l, inv_cld_frac_i, inv_cld_frac_r, ni_activated, inv_qc_relvar, cld_frac_i,
Expand Down
6 changes: 6 additions & 0 deletions components/eamxx/src/physics/p3/disp/p3_main_impl_part2_disp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,10 @@ ::p3_main_part2_disp(
const bool& do_prescribed_CCN,
const Scalar& dt,
const Scalar& inv_dt,
const bool& use_hetfrz_classnuc,
const uview_2d<const Spack>& hetfrz_immersion_nucleation_tend,
const uview_2d<const Spack>& hetfrz_contact_nucleation_tend,
const uview_2d<const Spack>& hetfrz_deposition_nucleation_tend,
const view_dnu_table& dnu_table_vals,
const view_ice_table& ice_table_vals,
const view_collect_table& collect_table_vals,
Expand Down Expand Up @@ -111,6 +115,8 @@ ::p3_main_part2_disp(
// main k-loop (for processes):
p3_main_part2(
team, nk_pack, max_total_ni, predictNc, do_prescribed_CCN, dt, inv_dt,
use_hetfrz_classnuc, ekat::subview(hetfrz_immersion_nucleation_tend, i),
ekat::subview(hetfrz_contact_nucleation_tend, i),ekat::subview(hetfrz_deposition_nucleation_tend, i),
dnu_table_vals, ice_table_vals, collect_table_vals, revap_table_vals,
ekat::subview(pres, i), ekat::subview(dpres, i), ekat::subview(dz, i), ekat::subview(nc_nuceat_tend, i), ekat::subview(inv_exner, i),
ekat::subview(exner, i), ekat::subview(inv_cld_frac_l, i), ekat::subview(inv_cld_frac_i, i), ekat::subview(inv_cld_frac_r, i),
Expand Down
23 changes: 23 additions & 0 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 @@ -91,6 +91,23 @@ void P3Microphysics::set_grids(const std::shared_ptr<const GridsManager> grids_m
add_field<Updated> ("qv_prev_micro_step", scalar3d_layout_mid, kg/kg, grid_name, ps);
add_field<Updated> ("T_prev_micro_step", scalar3d_layout_mid, K, grid_name, ps);

// Input from MAM4xx-ACI for heterogeneous freezing calculations
constexpr auto cm = m / 100;

// units of number mixing ratios of tracers
constexpr auto frz_unit = 1 / (cm * cm * cm * s);
// heterogeneous freezing by immersion nucleation [cm^-3 s^-1]
add_field<Required>("hetfrz_immersion_nucleation_tend", scalar3d_layout_mid,
frz_unit, grid_name, ps);

// heterogeneous freezing by contact nucleation [cm^-3 s^-1]
add_field<Required>("hetfrz_contact_nucleation_tend", scalar3d_layout_mid, frz_unit,
grid_name, ps);

// heterogeneous freezing by deposition nucleation [cm^-3 s^-1]
add_field<Required>("hetfrz_deposition_nucleation_tend", scalar3d_layout_mid,
frz_unit, grid_name, ps);
Comment on lines +100 to +109
Copy link
Contributor

Choose a reason for hiding this comment

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

@singhbalwinder this is the last potential important comment, otherwise, I think the PR is mergeable:

What happens if MAM is inactive? What happens to these variables? Maybe we should make these hidden behind if-else? IF so, make sure to do the same below (around 333):

  // Inputs for the heteogeneous freezing
  diag_inputs.hetfrz_immersion_nucleation_tend  = get_field_in("hetfrz_immersion_nucleation_tend").get_view<const Pack**>();
  diag_inputs.hetfrz_contact_nucleation_tend    = get_field_in("hetfrz_contact_nucleation_tend").get_view<const Pack**>();
  diag_inputs.hetfrz_deposition_nucleation_tend = get_field_in("hetfrz_deposition_nucleation_tend").get_view<const Pack**>();

where in the unused/uneeded case, you can use the buffer.unused or something like that

Copy link
Contributor Author

Choose a reason for hiding this comment

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

When MAM is inactive, values for these variables will be picked up from the namelist and will stay the same for the entire simulation. The simulation will run as usual with no impact from these variables.

I can hide it behind an if/else if it is not the desired behavior. Is there a way for P3 to know if MAM4 is active or not? Previously, when I discussed this, the design principle did not allow this, as each parameterization should be able to run independently without the knowledge of other parameterizations/processes. Otherwise, it makes the logic complex (e.g., to decide if different combinations of processes are valid or not). If it has changed, please let me know, and I will add if/else blocks.

Copy link
Contributor

@mahf708 mahf708 Feb 3, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

If they are defined, then there's no need to worry (I wasn't aware they would be defined; this is a little surprising to me tbh)

For hiding them, I would use the boolean you're adding here (hetfrz). It doesn't matter though, whether these statements are hidden or not, won't make a difference as far as I could tell. I was worried a CIME case would error out if these variables are here, but MAM is inactive. I guess we will pick this up in testing soon enough :)

Copy link
Contributor

@mahf708 mahf708 Feb 5, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

I think I would hide these behind the boolean you introduced in this PR. For the other branch of the if, set the buffer.unused or any other dummy values

 FAIL:
!m_add_time_dim
/__w/E3SM/E3SM/components/eamxx/src/share/io/scorpio_output.cpp:477
Error! Time-dependent output field 'hetfrz_contact_nucleation_tend' has not been initialized yet
.

 FAIL:
!m_add_time_dim
/__w/E3SM/E3SM/components/eamxx/src/share/io/scorpio_output.cpp:477
Error! Time-dependent output field 'hetfrz_contact_nucleation_tend' has not been initialized yet
.

 FAIL:
!m_add_time_dim
/__w/E3SM/E3SM/components/eamxx/src/share/io/scorpio_output.cpp:477
Error! Time-dependent output field 'hetfrz_contact_nucleation_tend' has not been initialized yet
.

 FAIL:
!m_add_time_dim
/__w/E3SM/E3SM/components/eamxx/src/share/io/scorpio_output.cpp:477
Error! Time-dependent output field 'hetfrz_contact_nucleation_tend' has not been initialized yet
.


// Diagnostic Outputs: (all fields are just outputs w.r.t. P3)
add_field<Updated>("precip_liq_surf_mass", scalar2d_layout, kg/m2, grid_name, "ACCUMULATED");
add_field<Updated>("precip_ice_surf_mass", scalar2d_layout, kg/m2, grid_name, "ACCUMULATED");
Expand Down Expand Up @@ -313,6 +330,12 @@ void P3Microphysics::initialize_impl (const RunType /* run_type */)
diag_inputs.cld_frac_r = p3_preproc.cld_frac_r;
diag_inputs.dz = p3_preproc.dz;
diag_inputs.inv_exner = p3_preproc.inv_exner;

// Inputs for the heteogeneous freezing
diag_inputs.hetfrz_immersion_nucleation_tend = get_field_in("hetfrz_immersion_nucleation_tend").get_view<const Pack**>();
diag_inputs.hetfrz_contact_nucleation_tend = get_field_in("hetfrz_contact_nucleation_tend").get_view<const Pack**>();
diag_inputs.hetfrz_deposition_nucleation_tend = get_field_in("hetfrz_deposition_nucleation_tend").get_view<const Pack**>();

// --Diagnostic Outputs
diag_outputs.diag_eff_radius_qc = get_field_out("eff_radius_qc").get_view<Pack**>();
diag_outputs.diag_eff_radius_qi = get_field_out("eff_radius_qi").get_view<Pack**>();
Expand Down
14 changes: 14 additions & 0 deletions components/eamxx/src/physics/p3/eti/p3_CNT_couple_impl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
#include "p3_CNT_couple_impl.hpp"

namespace scream {
namespace p3 {

/*
* Explicit instantiation for doing conservation functions on Reals using the
* default device.
*/

template struct Functions<Real,DefaultDevice>;

} // namespace p3
} // namespace scream
62 changes: 62 additions & 0 deletions components/eamxx/src/physics/p3/impl/p3_CNT_couple_impl.hpp
View file
Open in desktop
singhbalwinder marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
#ifndef P3_CNT_COUPLE_IMPL_HPP
#define P3_CNT_COUPLE_IMPL_HPP

#include "p3_functions.hpp" // for ETI only but harmless for GPU

namespace scream {
namespace p3 {

/*
*/

template <typename S, typename D>
KOKKOS_FUNCTION
void Functions<S,D>
::CNT_couple(
const Spack& frzimm, const Spack& frzcnt,
const Spack& frzdep, const Spack& rho,
const Spack& qc_incld, const Spack& nc_incld,
const int Iflag,
Spack& ncheti_cnt, Spack& qcheti_cnt,
Spack& nicnt, Spack& qicnt,
Spack& ninuc_cnt, Spack& qinuc_cnt)
{
constexpr Scalar pi = C::Pi;
constexpr Scalar rho_h2o = C::RHO_H2O;

// TODO: Verify if qsmall can be unified with other "small" numeric literals
constexpr Scalar qsmall = 1.0e-18;
constexpr Scalar piov3 = pi/3.0;

// TODO: Verify if 1.0e-18 can be unified with other "small" numeric literals
constexpr Scalar mi0 = 4.0*piov3*900.0*1.0e-18; // BAD_CONSTANT!

const Spack Zero(0.0);
// minimum mass of new crystal due to freezing of cloud droplets done
// externally (kg)

const Scalar mi0l_min = (4.0/3.0)*pi*rho_h2o*(4.0e-6)*(4.0e-6)*(4.0e-6);
Spack mi0l = qc_incld/ekat::max(nc_incld,1.0e6/rho);
mi0l = ekat::max(mi0l_min, mi0l);

const auto mask = qc_incld > qsmall;
switch (Iflag) {
case 1: // cloud droplet immersion freezing
ncheti_cnt.set(mask, frzimm*1.0e6/rho /* frzimm input is in [#/cm3] */ , Zero);
Copy link
Contributor

Choose a reason for hiding this comment

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

In all these "set" calls, how often do you expect the mask to be true/false? If the mask could often be ALL false (not sometimes, often), then you may consider using if statements, to avoid computing the packs for the true case for nothing (e.g., in the 1st line we have to compute frzimm*1e6/rho regardless of whether we need it or not).

Note: this nano-opt makes sense only if you expect mask to be often false. I assume that's not the case, since qsmall is very small. But I don't know how qc_incld is computed, so maybe it's often 0?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is a good question. I am not sure about that. @kaizhangpnl or @AaronDonahue might know if mask can often be false or not.

qcheti_cnt.set(mask, ncheti_cnt*mi0l, Zero);
break;
case 2: // deposition freezing / contact freezing
nicnt.set(mask, frzcnt*1.0e6/rho, Zero);
qicnt.set(mask, nicnt*mi0l, Zero);
ninuc_cnt.set(mask, frzdep*1.0e6/rho, Zero);
qinuc_cnt.set(mask, ninuc_cnt*mi0, Zero);
break;
default:
singhbalwinder marked this conversation as resolved.
Show resolved Hide resolved
EKAT_KERNEL_ERROR_MSG("Error! Unhandled case in switch statement for Iflag in p3_CNT_couple_impl.hpp .\n");
break;
}
}
} // namespace p3
} // namespace scream

#endif
10 changes: 9 additions & 1 deletion components/eamxx/src/physics/p3/impl/p3_back_to_cell_average_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ ::back_to_cell_average(
Spack& nr_collect_tend, Spack& ni_selfcollect_tend, Spack& qv2qi_vapdep_tend,
Spack& nr2ni_immers_freeze_tend, Spack& ni_sublim_tend, Spack& qv2qi_nucleat_tend,
Spack& ni_nucleat_tend, Spack& qc2qi_berg_tend,
const Smask& context)
Spack& ncheti_cnt, Spack& qcheti_cnt, Spack& nicnt, Spack& qicnt, Spack& ninuc_cnt,
Spack& qinuc_cnt, const Smask& context)
{
Spack ir_cldm, il_cldm, lr_cldm;
ir_cldm = min(cld_frac_i,cld_frac_r); // Intersection of ICE and RAIN cloud
Expand Down Expand Up @@ -70,6 +71,13 @@ ::back_to_cell_average(
ni_sublim_tend.set(context, ni_sublim_tend * cld_frac_i); // Number change due to sublimation of ice
qc2qi_berg_tend.set(context, qc2qi_berg_tend * il_cldm); // Bergeron process

ncheti_cnt.set(context,ncheti_cnt*cld_frac_l);
qcheti_cnt.set(context, qcheti_cnt*cld_frac_l);
nicnt.set(context, nicnt*cld_frac_l);
qicnt.set(context, qicnt*cld_frac_l);
ninuc_cnt.set(context, ninuc_cnt*cld_frac_l);
qinuc_cnt.set(context, qinuc_cnt*cld_frac_l);

// AaronDonahue: These variables are related to aerosol activation and their usage will be changed in a later PR.
//qv2qi_nucleat_tend = qv2qi_nucleat_tend; // Deposition and condensation-freezing nucleation, already cell-averaged
//ni_nucleat_tend = ni_nucleat_tend; // Number change due to deposition and condensation-freezing, already cell-averaged
Expand Down
33 changes: 27 additions & 6 deletions components/eamxx/src/physics/p3/impl/p3_conservation_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,16 @@ void Functions<S,D>
::cloud_water_conservation(const Spack& qc, const Scalar dt,
Spack& qc2qr_autoconv_tend, Spack& qc2qr_accret_tend, Spack &qc2qi_collect_tend, Spack& qc2qi_hetero_freeze_tend,
Spack& qc2qr_ice_shed_tend, Spack& qc2qi_berg_tend, Spack& qi2qv_sublim_tend, Spack& qv2qi_vapdep_tend,
const Smask& context)
Spack& qcheti_cnt, Spack& qicnt, const bool& use_hetfrz_classnuc, const Smask& context)
{
const auto sinks = (qc2qr_autoconv_tend+qc2qr_accret_tend+qc2qi_collect_tend+qc2qi_hetero_freeze_tend+qc2qr_ice_shed_tend+qc2qi_berg_tend)*dt; // Sinks of cloud water

Spack sinks;
if(use_hetfrz_classnuc){
sinks = (qc2qr_autoconv_tend+qc2qr_accret_tend+qc2qi_collect_tend+qcheti_cnt+qc2qr_ice_shed_tend+qc2qi_berg_tend)*dt; // Sinks of cloud water
}
else{
sinks = (qc2qr_autoconv_tend+qc2qr_accret_tend+qc2qi_collect_tend+qc2qi_hetero_freeze_tend+qc2qr_ice_shed_tend+qc2qi_berg_tend)*dt; // Sinks of cloud water
}
const auto sources = qc; // Source of cloud water
Spack ratio;

Expand All @@ -28,7 +35,13 @@ ::cloud_water_conservation(const Spack& qc, const Scalar dt,
qc2qr_autoconv_tend.set(enforce_conservation, qc2qr_autoconv_tend*ratio);
qc2qr_accret_tend.set(enforce_conservation, qc2qr_accret_tend*ratio);
qc2qi_collect_tend.set(enforce_conservation, qc2qi_collect_tend*ratio);
qc2qi_hetero_freeze_tend.set(enforce_conservation, qc2qi_hetero_freeze_tend*ratio);
if(use_hetfrz_classnuc){
qcheti_cnt.set(enforce_conservation, qcheti_cnt*ratio);
qicnt.set(enforce_conservation, qicnt*ratio);
}
else{
qc2qi_hetero_freeze_tend.set(enforce_conservation, qc2qi_hetero_freeze_tend*ratio);
}
qc2qr_ice_shed_tend.set(enforce_conservation, qc2qr_ice_shed_tend*ratio);
qc2qi_berg_tend.set(enforce_conservation, qc2qi_berg_tend*ratio);
}
Expand Down Expand Up @@ -78,13 +91,21 @@ void Functions<S,D>
::ice_water_conservation(
const Spack& qi,const Spack& qv2qi_vapdep_tend,const Spack& qv2qi_nucleat_tend,const Spack& qc2qi_berg_tend,
const Spack &qr2qi_collect_tend,const Spack &qc2qi_collect_tend,const Spack& qr2qi_immers_freeze_tend,
const Spack& qc2qi_hetero_freeze_tend,const Scalar dt,
Spack& qi2qv_sublim_tend, Spack& qi2qr_melt_tend,
const Spack& qc2qi_hetero_freeze_tend,const Scalar dt, Spack &qinuc_cnt, Spack &qcheti_cnt, Spack &qicnt,
Spack& qi2qv_sublim_tend, Spack& qi2qr_melt_tend, const bool& use_hetfrz_classnuc,
const Smask& context)
{
const auto sinks = (qi2qv_sublim_tend+qi2qr_melt_tend)*dt; // Sinks of ice water
const auto sources = qi + (qv2qi_vapdep_tend+qv2qi_nucleat_tend+qr2qi_collect_tend+qc2qi_collect_tend

Spack sources;
if(use_hetfrz_classnuc){
sources = qi + (qv2qi_vapdep_tend+qv2qi_nucleat_tend+qr2qi_collect_tend+qc2qi_collect_tend
+ qr2qi_immers_freeze_tend+qc2qi_berg_tend+qinuc_cnt+qcheti_cnt+qicnt)*dt; // Sources of ice water
}
else{
sources = qi + (qv2qi_vapdep_tend+qv2qi_nucleat_tend+qr2qi_collect_tend+qc2qi_collect_tend
+ qr2qi_immers_freeze_tend+qc2qi_hetero_freeze_tend+qc2qi_berg_tend)*dt; // Sources of ice water
}
Spack ratio;
constexpr Scalar qtendsmall = C::QTENDSMALL;
Smask enforce_conservation = sinks > sources && sinks >= qtendsmall && context; // determine if conservation corrction is necessary
Expand Down
13 changes: 11 additions & 2 deletions components/eamxx/src/physics/p3/impl/p3_ice_supersat_conservation_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ namespace p3 {

template<typename S, typename D>
KOKKOS_FUNCTION
void Functions<S,D>::ice_supersat_conservation(Spack& qv2qi_vapdep_tend, Spack& qv2qi_nucleat_tend, const Spack& cld_frac_i, const Spack& qv, const Spack& qv_sat_i, const Spack& t_atm, const Real& dt, const Spack& qi2qv_sublim_tend, const Spack& qr2qv_evap_tend, const Smask& context)
void Functions<S,D>::ice_supersat_conservation(Spack& qv2qi_vapdep_tend, Spack& qv2qi_nucleat_tend, Spack& qinuc_cnt, const Spack& cld_frac_i, const Spack& qv, const Spack& qv_sat_i, const Spack& t_atm, const Real& dt, const Spack& qi2qv_sublim_tend, const Spack& qr2qv_evap_tend, const bool& use_hetfrz_classnuc, const Smask& context)
{
constexpr Scalar qsmall = C::QSMALL;
constexpr Scalar cp = C::CP;
Expand All @@ -22,7 +22,13 @@ void Functions<S,D>::ice_supersat_conservation(Spack& qv2qi_vapdep_tend, Spack&
constexpr Scalar latice = C::LatIce;
constexpr Scalar latsublim2 = (latvap+latice)*(latvap+latice);

const auto qv_sink = qv2qi_vapdep_tend + qv2qi_nucleat_tend; // in [kg/kg] cell-avg values
Spack qv_sink;
if(use_hetfrz_classnuc){
qv_sink = qv2qi_vapdep_tend + qv2qi_nucleat_tend + qinuc_cnt; // in [kg/kg] cell-avg values
}
else{
qv_sink = qv2qi_vapdep_tend + qv2qi_nucleat_tend; // in [kg/kg] cell-avg values
}

const auto mask = qv_sink > qsmall && cld_frac_i > 1e-20 && context;
if (mask.any()) {
Expand All @@ -36,6 +42,9 @@ void Functions<S,D>::ice_supersat_conservation(Spack& qv2qi_vapdep_tend, Spack&
const auto sink_gt_avail = qv_sink > qv_avail && mask;
if (sink_gt_avail.any()) {
const auto fract = qv_avail / qv_sink;
if(use_hetfrz_classnuc){
qinuc_cnt.set(sink_gt_avail, qinuc_cnt * fract);
}
qv2qi_nucleat_tend.set(sink_gt_avail, qv2qi_nucleat_tend * fract);
qv2qi_vapdep_tend.set(sink_gt_avail, qv2qi_vapdep_tend * fract);
}
Expand Down
6 changes: 6 additions & 0 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 @@ -201,6 +201,11 @@ ::p3_main_internal(
const auto oqv_prev = ekat::subview(diagnostic_inputs.qv_prev, i);
const auto ot_prev = ekat::subview(diagnostic_inputs.t_prev, i);

// Inputs for the heteogeneous freezing
const auto ohetfrz_immersion_nucleation_tend = ekat::subview(diagnostic_inputs.hetfrz_immersion_nucleation_tend, i);
const auto ohetfrz_contact_nucleation_tend = ekat::subview(diagnostic_inputs.hetfrz_contact_nucleation_tend, i);
const auto ohetfrz_deposition_nucleation_tend = ekat::subview(diagnostic_inputs.hetfrz_deposition_nucleation_tend, i);

// Use Kokkos' scratch pad for allocating 2 bools
// per team to determine early exits
ScratchViewType bools(team.team_scratch(0), 2);
Expand Down Expand Up @@ -243,6 +248,7 @@ ::p3_main_internal(

p3_main_part2(
team, nk_pack, runtime_options.max_total_ni, infrastructure.predictNc, infrastructure.prescribedCCN, infrastructure.dt, inv_dt,
runtime_options.use_hetfrz_classnuc, ohetfrz_immersion_nucleation_tend, ohetfrz_contact_nucleation_tend, ohetfrz_deposition_nucleation_tend,
singhbalwinder marked this conversation as resolved.
Show resolved Hide resolved
lookup_tables.dnu_table_vals, lookup_tables.ice_table_vals, lookup_tables.collect_table_vals, lookup_tables.revap_table_vals, opres, odpres, odz, onc_nuceat_tend, oinv_exner,
exner, inv_cld_frac_l, inv_cld_frac_i, inv_cld_frac_r, oni_activated, oinv_qc_relvar, ocld_frac_i,
ocld_frac_l, ocld_frac_r, oqv_prev, ot_prev, T_atm, rho, inv_rho, qv_sat_l, qv_sat_i, qv_supersat_i, rhofacr, rhofaci, acn,
Expand Down
Loading
Loading