EAMxx: add p3 process rates as outputs

components/eamxx/cime_config/namelist_defaults_scream.xml

@@ -227,6 +227,7 @@ be lost if SCREAM_HACK_XML is not enabled.
       <set_cld_frac_l_to_one type="logical" doc="set P3 input liquid cloud fraction to 1 everywhere">false</set_cld_frac_l_to_one>
       <set_cld_frac_r_to_one type="logical" doc="set P3 input rain cloud fraction to 1 everywhere"  >false</set_cld_frac_r_to_one>
       <set_cld_frac_i_to_one type="logical" doc="set P3 input ice cloud fraction to 1 everywhere"   >false</set_cld_frac_i_to_one>
+      <p3_extra_diags type="logical" doc="Extra P3 diagnostics">false</p3_extra_diags>
     </p3>
 
     <!-- SHOC macrophysics -->

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

@@ -166,6 +166,18 @@ ::p3_main_internal_disp(
   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;
+  auto P3_qr2qv_evap           = history_only.P3_qr2qv_evap;
+  auto P3_qi2qv_sublim         = history_only.P3_qi2qv_sublim;
+  auto P3_qc2qr_accret         = history_only.P3_qc2qr_accret;
+  auto P3_qc2qr_autoconv       = history_only.P3_qc2qr_autoconv;
+  auto P3_qv2qi_vapdep         = history_only.P3_qv2qi_vapdep;
+  auto P3_qc2qi_berg           = history_only.P3_qc2qi_berg;
+  auto P3_qc2qr_ice_shed       = history_only.P3_qc2qr_ice_shed;
+  auto P3_qc2qi_collect        = history_only.P3_qc2qi_collect;
+  auto P3_qr2qi_collect        = history_only.P3_qr2qi_collect;
+  auto P3_qc2qi_hetero_freeze  = history_only.P3_qc2qi_hetero_freeze;
+  auto P3_qr2qi_immers_freeze  = history_only.P3_qr2qi_immers_freeze;
+  auto P3_qi2qr_melt           = history_only.P3_qi2qr_melt;
   auto mu_r                    = temporaries.mu_r;
   auto T_atm                   = temporaries.T_atm;
   auto lamr                    = temporaries.lamr;
@@ -260,6 +272,8 @@ ::p3_main_internal_disp(
       nr_incld, ni_incld, bm_incld, mu_c, nu, lamc, cdist, cdist1, cdistr,
       mu_r, lamr, logn0r, qv2qi_depos_tend, precip_total_tend, nevapr, qr_evap_tend,
       vap_liq_exchange, vap_ice_exchange, liq_ice_exchange,
+      P3_qr2qv_evap, P3_qi2qv_sublim, P3_qc2qr_accret, P3_qc2qr_autoconv, P3_qv2qi_vapdep, P3_qc2qi_berg,
+      P3_qc2qr_ice_shed, P3_qc2qi_collect, P3_qr2qi_collect, P3_qc2qi_hetero_freeze, P3_qr2qi_immers_freeze, P3_qi2qr_melt,
       pratot, prctot, nucleationPossible, hydrometeorsPresent, runtime_options);
 
   //NOTE: At this point, it is possible to have negative (but small) nc, nr, ni.  This is not

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

@@ -86,6 +86,18 @@ ::p3_main_part2_disp(
   const uview_2d<Spack>& vap_liq_exchange,
   const uview_2d<Spack>& vap_ice_exchange,
   const uview_2d<Spack>& liq_ice_exchange,
+  const uview_2d<Spack>& P3_qr2qv_evap,
+  const uview_2d<Spack>& P3_qi2qv_sublim,
+  const uview_2d<Spack>& P3_qc2qr_accret,
+  const uview_2d<Spack>& P3_qc2qr_autoconv,
+  const uview_2d<Spack>& P3_qv2qi_vapdep,
+  const uview_2d<Spack>& P3_qc2qi_berg,
+  const uview_2d<Spack>& P3_qc2qr_ice_shed,
+  const uview_2d<Spack>& P3_qc2qi_collect,
+  const uview_2d<Spack>& P3_qr2qi_collect,
+  const uview_2d<Spack>& P3_qc2qi_hetero_freeze,
+  const uview_2d<Spack>& P3_qr2qi_immers_freeze,
+  const uview_2d<Spack>& P3_qi2qr_melt,
   const uview_2d<Spack>& pratot,
   const uview_2d<Spack>& prctot,
   const uview_1d<bool>& nucleationPossible,
@@ -124,7 +136,9 @@ ::p3_main_part2_disp(
       ekat::subview(ni_incld, i), ekat::subview(bm_incld, i), ekat::subview(mu_c, i), ekat::subview(nu, i), ekat::subview(lamc, i), ekat::subview(cdist, i),
       ekat::subview(cdist1, i), ekat::subview(cdistr, i), ekat::subview(mu_r, i), ekat::subview(lamr, i), ekat::subview(logn0r, i),
       ekat::subview(qv2qi_depos_tend, i), ekat::subview(precip_total_tend, i),
-      ekat::subview(nevapr, i), ekat::subview(qr_evap_tend, i), ekat::subview(vap_liq_exchange, i), ekat::subview(vap_ice_exchange, i), ekat::subview(liq_ice_exchange, i),
+      ekat::subview(nevapr, i), ekat::subview(qr_evap_tend, i), ekat::subview(vap_liq_exchange, i), ekat::subview(vap_ice_exchange, i), ekat::subview(liq_ice_exchange, i), ekat::subview(P3_qr2qv_evap, i), ekat::subview(P3_qi2qv_sublim, i), ekat::subview(P3_qc2qr_accret,i),  ekat::subview(P3_qc2qr_autoconv,i), ekat::subview(P3_qv2qi_vapdep,i), ekat::subview(P3_qc2qi_berg,i),
+      ekat::subview(P3_qc2qr_ice_shed,i), ekat::subview(P3_qc2qi_collect,i), ekat::subview(P3_qr2qi_collect,i), ekat::subview(P3_qc2qi_hetero_freeze,i), 
+      ekat::subview(P3_qr2qi_immers_freeze,i), ekat::subview(P3_qi2qr_melt,i),
       ekat::subview(pratot, i), ekat::subview(prctot, i), hydrometeorsPresent(i), nk, runtime_options);
 
     if (!hydrometeorsPresent(i)) return;

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

@@ -35,6 +35,9 @@ void P3Microphysics::set_grids(const std::shared_ptr<const GridsManager> grids_m
   m_num_cols = m_grid->get_num_local_dofs(); // Number of columns on this rank
   m_num_levs = m_grid->get_num_vertical_levels();  // Number of levels per column
 
+  // Gather runtime options from file
+  runtime_options.load_runtime_options_from_file(m_params);
+
   // --Infrastructure
   // dt is passed as an argument to run_impl
   infrastructure.it  = 0;
@@ -108,6 +111,24 @@ void P3Microphysics::set_grids(const std::shared_ptr<const GridsManager> grids_m
   add_field<Computed>("micro_vap_liq_exchange", scalar3d_layout_mid, kg/kg,  grid_name, ps);
   add_field<Computed>("micro_vap_ice_exchange", scalar3d_layout_mid, kg/kg,  grid_name, ps);
   add_field<Computed>("rainfrac",               scalar3d_layout_mid, nondim, grid_name, ps);
+  // P3 process rate diagnostics, use convention P3_<process_rate>
+  if (runtime_options.p3_extra_diags){
+    add_field<Computed>("P3_qr2qv_evap", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qi2qv_sublim", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qr_accret", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qr_autoconv", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qv2qi_vapdep", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qi_berg", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qr_ice_shed", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qi_collect", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qr2qi_collect", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc2qi_hetero_freeze", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qr2qi_immers_freeze", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qi2qr_melt", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qr_sed", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qc_sed", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+    add_field<Computed>("P3_qi_sed", scalar3d_layout_mid, kg/kg/s,  grid_name, ps);
+  }
 
   // Boundary flux fields for energy and mass conservation checks
   if (has_column_conservation_check()) {
@@ -217,8 +238,6 @@ void P3Microphysics::init_buffers(const ATMBufferManager &buffer_manager)
 // =========================================================================================
 void P3Microphysics::initialize_impl (const RunType /* run_type */)
 {
-  // Gather runtime options from file
-  runtime_options.load_runtime_options_from_file(m_params);
 
   // Set property checks for fields in this process
   add_invariant_check<FieldWithinIntervalCheck>(get_field_out("T_mid"),m_grid,100.0,500.0,false);
@@ -332,6 +351,45 @@ void P3Microphysics::initialize_impl (const RunType /* run_type */)
   history_only.liq_ice_exchange = get_field_out("micro_liq_ice_exchange").get_view<Pack**>();
   history_only.vap_liq_exchange = get_field_out("micro_vap_liq_exchange").get_view<Pack**>();
   history_only.vap_ice_exchange = get_field_out("micro_vap_ice_exchange").get_view<Pack**>();
+  if(runtime_options.p3_extra_diags) {
+    // if we are doing extra diagnostics, assign the fields to the history only struct
+    history_only.P3_qr2qv_evap   = get_field_out("P3_qr2qv_evap").get_view<Pack**>();
+    history_only.P3_qi2qv_sublim = get_field_out("P3_qi2qv_sublim").get_view<Pack**>();
+    history_only.P3_qc2qr_accret = get_field_out("P3_qc2qr_accret").get_view<Pack**>();
+    history_only.P3_qc2qr_autoconv = get_field_out("P3_qc2qr_autoconv").get_view<Pack**>();
+    history_only.P3_qv2qi_vapdep = get_field_out("P3_qv2qi_vapdep").get_view<Pack**>();
+    history_only.P3_qc2qi_berg = get_field_out("P3_qc2qi_berg").get_view<Pack**>();
+    history_only.P3_qc2qr_ice_shed = get_field_out("P3_qc2qr_ice_shed").get_view<Pack**>();
+    history_only.P3_qc2qi_collect = get_field_out("P3_qc2qi_collect").get_view<Pack**>();
+    history_only.P3_qr2qi_collect = get_field_out("P3_qr2qi_collect").get_view<Pack**>();
+    history_only.P3_qc2qi_hetero_freeze = get_field_out("P3_qc2qi_hetero_freeze").get_view<Pack**>();
+    history_only.P3_qr2qi_immers_freeze = get_field_out("P3_qr2qi_immers_freeze").get_view<Pack**>();
+    history_only.P3_qi2qr_melt = get_field_out("P3_qi2qr_melt").get_view<Pack**>();
+    history_only.P3_qr_sed = get_field_out("P3_qr_sed").get_view<Pack**>();
+    history_only.P3_qc_sed = get_field_out("P3_qc_sed").get_view<Pack**>();
+    history_only.P3_qi_sed = get_field_out("P3_qi_sed").get_view<Pack**>();
+  } else {
+    // if not, let's just use the buffer for the unused? field
+    // TODO: check if this is actually okay and doesn't have uintended consequences
+    // if we are not outputing these fields, we really don't care about their values
+    // but would this have side effects or memory issues? idk
+    // TODO: maybe just use more buffer and assign stuff?
+    history_only.P3_qr2qv_evap   = m_buffer.unused;
+    history_only.P3_qi2qv_sublim = m_buffer.unused;
+    history_only.P3_qc2qr_accret = m_buffer.unused;
+    history_only.P3_qc2qr_autoconv = m_buffer.unused;
+    history_only.P3_qv2qi_vapdep = m_buffer.unused;
+    history_only.P3_qc2qi_berg = m_buffer.unused;
+    history_only.P3_qc2qr_ice_shed = m_buffer.unused;
+    history_only.P3_qc2qi_collect = m_buffer.unused;
+    history_only.P3_qr2qi_collect = m_buffer.unused;
+    history_only.P3_qc2qi_hetero_freeze = m_buffer.unused;
+    history_only.P3_qr2qi_immers_freeze = m_buffer.unused;
+    history_only.P3_qi2qr_melt = m_buffer.unused;
+    history_only.P3_qr_sed = m_buffer.unused;
+    history_only.P3_qc_sed = m_buffer.unused;
+    history_only.P3_qi_sed = m_buffer.unused;
+  }
 #ifdef SCREAM_P3_SMALL_KERNELS
   // Temporaries
   temporaries.mu_r                    = m_buffer.mu_r;

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

@@ -28,6 +28,25 @@ void P3Microphysics::run_impl (const double dt)
   get_field_out("micro_vap_liq_exchange").deep_copy(0.0);
   get_field_out("micro_vap_ice_exchange").deep_copy(0.0);
 
+  // Optional extra p3 diags
+  if (runtime_options.p3_extra_diags) {
+    get_field_out("P3_qr2qv_evap").deep_copy(0.0);
+    get_field_out("P3_qi2qv_sublim").deep_copy(0.0);
+    get_field_out("P3_qc2qr_accret").deep_copy(0.0);
+    get_field_out("P3_qc2qr_autoconv").deep_copy(0.0);
+    get_field_out("P3_qv2qi_vapdep").deep_copy(0.0);
+    get_field_out("P3_qc2qi_berg").deep_copy(0.0);
+    get_field_out("P3_qc2qr_ice_shed").deep_copy(0.0);
+    get_field_out("P3_qc2qi_collect").deep_copy(0.0);
+    get_field_out("P3_qr2qi_collect").deep_copy(0.0);
+    get_field_out("P3_qc2qi_hetero_freeze").deep_copy(0.0);
+    get_field_out("P3_qr2qi_immers_freeze").deep_copy(0.0);
+    get_field_out("P3_qi2qr_melt").deep_copy(0.0);
+    get_field_out("P3_qr_sed").deep_copy(0.0);
+    get_field_out("P3_qc_sed").deep_copy(0.0);
+    get_field_out("P3_qi_sed").deep_copy(0.0);
+  }
+
   P3F::p3_main(runtime_options, prog_state, diag_inputs, diag_outputs, infrastructure,
                history_only, lookup_tables,
 #ifdef SCREAM_P3_SMALL_KERNELS

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

@@ -198,6 +198,21 @@ ::p3_main_internal(
     const auto oliq_ice_exchange   = ekat::subview(history_only.liq_ice_exchange, i);
     const auto ovap_liq_exchange   = ekat::subview(history_only.vap_liq_exchange, i);
     const auto ovap_ice_exchange   = ekat::subview(history_only.vap_ice_exchange, i);
+    const auto oP3_qr2qv_evap      = ekat::subview(history_only.P3_qr2qv_evap, i);
+    const auto oP3_qi2qv_sublim    = ekat::subview(history_only.P3_qi2qv_sublim, i);
+    const auto oP3_qc2qr_accret    = ekat::subview(history_only.P3_qc2qr_accret,i);
+    const auto oP3_qc2qr_autoconv  = ekat::subview(history_only.P3_qc2qr_autoconv,i);
+    const auto oP3_qv2qi_vapdep    = ekat::subview(history_only.P3_qv2qi_vapdep,i);
+    const auto oP3_qc2qi_berg      = ekat::subview(history_only.P3_qc2qi_berg,i);
+    const auto oP3_qc2qr_ice_shed  = ekat::subview(history_only.P3_qc2qr_ice_shed,i);
+    const auto oP3_qc2qi_collect   = ekat::subview(history_only.P3_qc2qi_collect,i);
+    const auto oP3_qr2qi_collect   = ekat::subview(history_only.P3_qr2qi_collect,i);
+    const auto oP3_qc2qi_hetero_freeze = ekat::subview(history_only.P3_qc2qi_hetero_freeze,i);
+    const auto oP3_qr2qi_immers_freeze = ekat::subview(history_only.P3_qr2qi_immers_freeze,i);
+    const auto oP3_qi2qr_melt      =  ekat::subview(history_only.P3_qi2qr_melt,i);
+    const auto oP3_qr_sed          = ekat::subview(history_only.P3_qr_sed, i);
+    const auto oP3_qc_sed          = ekat::subview(history_only.P3_qc_sed, i);
+    const auto oP3_qi_sed          = ekat::subview(history_only.P3_qi_sed, i);
     const auto oqv_prev            = ekat::subview(diagnostic_inputs.qv_prev, i);
     const auto ot_prev             = ekat::subview(diagnostic_inputs.t_prev, i);
 
@@ -250,7 +265,8 @@ ::p3_main_internal(
       qc_incld, qr_incld, qi_incld, qm_incld, nc_incld,
       nr_incld, ni_incld, bm_incld, mu_c, nu, lamc, cdist, cdist1, cdistr,
       mu_r, lamr, logn0r, oqv2qi_depos_tend, oprecip_total_tend, onevapr, qr_evap_tend,
-      ovap_liq_exchange, ovap_ice_exchange, oliq_ice_exchange,
+      ovap_liq_exchange, ovap_ice_exchange, oliq_ice_exchange,oP3_qr2qv_evap,oP3_qi2qv_sublim,oP3_qc2qr_accret, oP3_qc2qr_autoconv, oP3_qv2qi_vapdep,
+      oP3_qc2qi_berg, oP3_qc2qr_ice_shed, oP3_qc2qi_collect, oP3_qr2qi_collect, oP3_qc2qi_hetero_freeze, oP3_qr2qi_immers_freeze, oP3_qi2qr_melt,
       pratot, prctot, hydrometeorsPresent, nk, runtime_options);
 
     //NOTE: At this point, it is possible to have negative (but small) nc, nr, ni.  This is not
@@ -271,21 +287,21 @@ ::p3_main_internal(
     cloud_sedimentation(
       qc_incld, rho, inv_rho, ocld_frac_l, acn, inv_dz, lookup_tables.dnu_table_vals, team, workspace,
       nk, ktop, kbot, kdir, infrastructure.dt, inv_dt, infrastructure.predictNc,
-      oqc, onc, nc_incld, mu_c, lamc, qtend_ignore, ntend_ignore,
+      oqc, onc, nc_incld, mu_c, lamc, oP3_qc_sed, ntend_ignore,
       diagnostic_outputs.precip_liq_surf(i));
 
     // Rain sedimentation:  (adaptive substepping)
     rain_sedimentation(
       rho, inv_rho, rhofacr, ocld_frac_r, inv_dz, qr_incld, team, workspace,
       lookup_tables.vn_table_vals, lookup_tables.vm_table_vals, nk, ktop, kbot, kdir, infrastructure.dt, inv_dt, oqr,
-      onr, nr_incld, mu_r, lamr, oprecip_liq_flux, qtend_ignore, ntend_ignore,
+      onr, nr_incld, mu_r, lamr, oprecip_liq_flux, oP3_qr_sed, ntend_ignore,
       diagnostic_outputs.precip_liq_surf(i), runtime_options);
 
     // Ice sedimentation:  (adaptive substepping)
     ice_sedimentation(
       rho, inv_rho, rhofaci, ocld_frac_i, inv_dz, team, workspace, nk, ktop, kbot,
       kdir, infrastructure.dt, inv_dt, oqi, qi_incld, oni, ni_incld,
-      oqm, qm_incld, obm, bm_incld, qtend_ignore, ntend_ignore,
+      oqm, qm_incld, obm, bm_incld, oP3_qi_sed, ntend_ignore,
       lookup_tables.ice_table_vals, diagnostic_outputs.precip_ice_surf(i), runtime_options);
 
     // homogeneous freezing of cloud and rain

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

@@ -89,6 +89,18 @@ ::p3_main_part2(
   const uview_1d<Spack>& vap_liq_exchange,
   const uview_1d<Spack>& vap_ice_exchange,
   const uview_1d<Spack>& liq_ice_exchange,
+  const uview_1d<Spack>& P3_qr2qv_evap,
+  const uview_1d<Spack>& P3_qi2qv_sublim,
+  const uview_1d<Spack>& P3_qc2qr_accret,
+  const uview_1d<Spack>& P3_qc2qr_autoconv,
+  const uview_1d<Spack>& P3_qv2qi_vapdep,
+  const uview_1d<Spack>& P3_qc2qi_berg,
+  const uview_1d<Spack>& P3_qc2qr_ice_shed,
+  const uview_1d<Spack>& P3_qc2qi_collect,
+  const uview_1d<Spack>& P3_qr2qi_collect,
+  const uview_1d<Spack>& P3_qc2qi_hetero_freeze,
+  const uview_1d<Spack>& P3_qr2qi_immers_freeze,
+  const uview_1d<Spack>& P3_qi2qr_melt,
   const uview_1d<Spack>& pratot,
   const uview_1d<Spack>& prctot,
   bool& hydrometeorsPresent, const Int& nk,
@@ -469,6 +481,20 @@ ::p3_main_part2(
     vap_liq_exchange(k).set(not_skip_all, -qr2qv_evap_tend);
     liq_ice_exchange(k).set(not_skip_all, qc2qi_hetero_freeze_tend + qr2qi_immers_freeze_tend - qi2qr_melt_tend + qc2qi_berg_tend + qc2qi_collect_tend + qr2qi_collect_tend);
 
+   // set tendencies
+    P3_qr2qv_evap(k).set(not_skip_all, qr2qv_evap_tend);
+    P3_qi2qv_sublim(k).set(not_skip_all, qi2qv_sublim_tend);
+    P3_qc2qr_accret(k).set(not_skip_all, qc2qr_accret_tend);
+    P3_qc2qr_autoconv(k).set(not_skip_all, qc2qr_autoconv_tend);
+    P3_qv2qi_vapdep(k).set(not_skip_all, qv2qi_vapdep_tend);
+    P3_qc2qi_berg(k).set(not_skip_all, qc2qi_berg_tend);
+    P3_qc2qr_ice_shed(k).set(not_skip_all, qc2qr_ice_shed_tend);
+    P3_qc2qi_collect(k).set(not_skip_all, qc2qi_collect_tend);
+    P3_qr2qi_collect(k).set(not_skip_all, qr2qi_collect_tend);
+    P3_qc2qi_hetero_freeze(k).set(not_skip_all, qc2qi_hetero_freeze_tend);
+    P3_qr2qi_immers_freeze(k).set(not_skip_all, qr2qi_immers_freeze_tend);
+    P3_qi2qr_melt(k).set(not_skip_all, qi2qr_melt_tend);
+
     // clipping for small hydrometeor values
     const auto qc_small    = qc(k) < qsmall    && not_skip_all;
     const auto qr_small    = qr(k) < qsmall    && not_skip_all;