EAMxx: Adds aerosols heterogeneous freezing calculations in P3 microphysics

components/eamxx/cime_config/namelist_defaults_scream.xml

@@ -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>

components/eamxx/src/physics/p3/CMakeLists.txt

@@ -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_ice_classical_nucleation_impl.cpp
     eti/p3_cldliq_imm_freezing.cpp
     eti/p3_rain_imm_freezing.cpp
     eti/p3_droplet_self_coll.cpp

components/eamxx/src/physics/p3/disp/p3_main_impl_disp.cpp

@@ -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;
@@ -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,
+      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,

components/eamxx/src/physics/p3/disp/p3_main_impl_part2_disp.cpp

@@ -23,6 +23,9 @@ ::p3_main_part2_disp(
   const bool& do_prescribed_CCN,
   const Scalar& dt,
   const Scalar& inv_dt,
+  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,
@@ -111,6 +114,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,
+      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),

components/eamxx/src/physics/p3/eamxx_p3_process_interface.cpp

@@ -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);
+
   // 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");
@@ -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**>();

components/eamxx/src/physics/p3/eti/p3_ice_classical_nucleation_impl.cpp

@@ -0,0 +1,14 @@
+#include "p3_ice_classical_nucleation_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

components/eamxx/src/physics/p3/impl/p3_back_to_cell_average_impl.hpp

@@ -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
@@ -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

components/eamxx/src/physics/p3/impl/p3_conservation_impl.hpp

@@ -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;
 
@@ -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);
   }
@@ -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

components/eamxx/src/physics/p3/impl/p3_ice_classical_nucleation_impl.hpp

@@ -0,0 +1,62 @@
+#ifndef P3_ICE_CLASSICAL_NUCLEATION_IMPL_HPP
+#define P3_ICE_CLASSICAL_NUCLEATION_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>
+::ice_classical_nucleation(
+  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);
+      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:
+      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

components/eamxx/src/physics/p3/impl/p3_ice_supersat_conservation_impl.hpp

@@ -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;
@@ -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()) {
@@ -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);
     }

components/eamxx/src/physics/p3/impl/p3_main_impl.hpp

@@ -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);
@@ -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,
+      ohetfrz_immersion_nucleation_tend, ohetfrz_contact_nucleation_tend, ohetfrz_deposition_nucleation_tend,
       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,