EAMxx: add ability to use separate cloud fractions in p3

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>
+      <use_separate_ice_liq_frac type="logical" doc="use separate ice and liquid cloud fractions from shoc">false</use_separate_ice_liq_frac>
     </p3>
 
     <!-- SHOC macrophysics -->

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

@@ -62,6 +62,8 @@ void P3Microphysics::set_grids(const std::shared_ptr<const GridsManager> grids_m
 
   // These variables are needed by the interface, but not actually passed to p3_main.
   add_field<Required>("cldfrac_tot", scalar3d_layout_mid, nondim, grid_name, ps);
+  add_field<Required>("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name, ps);
+  add_field<Required>("cldfrac_ice", scalar3d_layout_mid, nondim, grid_name, ps);
 
 //should we use one pressure only, wet/full?
   add_field<Required>("p_mid",       scalar3d_layout_mid, Pa,     grid_name, ps);
@@ -254,7 +256,9 @@ void P3Microphysics::initialize_impl (const RunType /* run_type */)
   const  auto& pseudo_density = get_field_in("pseudo_density").get_view<const Pack**>();
   const  auto& pseudo_density_dry = get_field_in("pseudo_density_dry").get_view<const Pack**>();
   const  auto& T_atm          = get_field_out("T_mid").get_view<Pack**>();
-  const  auto& cld_frac_t     = get_field_in("cldfrac_tot").get_view<const Pack**>();
+  const  auto& cld_frac_t_in  = get_field_in("cldfrac_tot").get_view<const Pack**>();
+  const  auto& cld_frac_l_in  = get_field_in("cldfrac_liq").get_view<const Pack**>();
+  const  auto& cld_frac_i_in  = get_field_in("cldfrac_ice").get_view<const Pack**>();
   const  auto& qv             = get_field_out("qv").get_view<Pack**>();
   const  auto& qc             = get_field_out("qc").get_view<Pack**>();
   const  auto& nc             = get_field_out("nc").get_view<Pack**>();
@@ -278,7 +282,7 @@ void P3Microphysics::initialize_impl (const RunType /* run_type */)
 
   // -- Set values for the pre-amble structure
   p3_preproc.set_variables(m_num_cols,nk_pack,pmid,pmid_dry,pseudo_density,pseudo_density_dry,
-                        T_atm,cld_frac_t,
+                        T_atm,cld_frac_t_in,cld_frac_l_in,cld_frac_i_in,
                         qv, qc, nc, qr, nr, qi, qm, ni, bm, qv_prev,
                         inv_exner, th_atm, cld_frac_l, cld_frac_i, cld_frac_r, dz, runtime_options);
   // --Prognostic State Variables:

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

@@ -69,7 +69,9 @@ class P3Microphysics : public AtmosphereProcess
         // The ipack slice of input variables used more than once
         const Spack& pmid_pack(pmid(icol,ipack));
         const Spack& T_atm_pack(T_atm(icol,ipack));
-        const Spack& cld_frac_t_pack(cld_frac_t(icol,ipack));
+        const Spack& cld_frac_t_in_pack(cld_frac_t_in(icol,ipack));
+        const Spack& cld_frac_l_in_pack(cld_frac_l_in(icol,ipack));
+        const Spack& cld_frac_i_in_pack(cld_frac_i_in(icol,ipack));
         const Spack& pseudo_density_pack(pseudo_density(icol,ipack));
         const Spack& pseudo_density_dry_pack(pseudo_density_dry(icol,ipack));
 
@@ -106,9 +108,14 @@ class P3Microphysics : public AtmosphereProcess
         // Cloud fraction
         // Set minimum cloud fraction - avoids division by zero
         // Alternatively set fraction to 1 everywhere to disable subgrid effects
-        cld_frac_l(icol,ipack) = runtime_opts.set_cld_frac_l_to_one ? 1 : ekat::max(cld_frac_t_pack,mincld);
-        cld_frac_i(icol,ipack) = runtime_opts.set_cld_frac_i_to_one ? 1 : ekat::max(cld_frac_t_pack,mincld);
-        cld_frac_r(icol,ipack) = runtime_opts.set_cld_frac_r_to_one ? 1 : ekat::max(cld_frac_t_pack,mincld);
+        if (runtime_opts.use_separate_ice_liq_frac){
+          cld_frac_l(icol,ipack) = runtime_opts.set_cld_frac_l_to_one ? 1 : ekat::max(cld_frac_l_in_pack,mincld);
+          cld_frac_i(icol,ipack) = runtime_opts.set_cld_frac_i_to_one ? 1 : ekat::max(cld_frac_i_in_pack,mincld);
+        } else {
+          cld_frac_l(icol,ipack) = runtime_opts.set_cld_frac_l_to_one ? 1 : ekat::max(cld_frac_t_in_pack,mincld);
+          cld_frac_i(icol,ipack) = runtime_opts.set_cld_frac_i_to_one ? 1 : ekat::max(cld_frac_t_in_pack,mincld);
+        }
+        cld_frac_r(icol,ipack) = runtime_opts.set_cld_frac_r_to_one ? 1 : ekat::max(cld_frac_t_in_pack,mincld);
 
         // update rain cloud fraction given neighboring levels using max-overlap approach.
         if ( !runtime_opts.set_cld_frac_r_to_one ) {
@@ -121,8 +128,8 @@ class P3Microphysics : public AtmosphereProcess
             Int ipack_m1 = (lev - 1) / Spack::n;
             Int ivec_m1  = (lev - 1) % Spack::n;
             if (lev != 0) { /* Not applicable at the very top layer */
-              cld_frac_r(icol,ipack)[ivec] = cld_frac_t(icol,ipack_m1)[ivec_m1]>cld_frac_r(icol,ipack)[ivec] ?
-                                                cld_frac_t(icol,ipack_m1)[ivec_m1] :
+              cld_frac_r(icol,ipack)[ivec] = cld_frac_t_in(icol,ipack_m1)[ivec_m1]>cld_frac_r(icol,ipack)[ivec] ?
+                                                cld_frac_t_in(icol,ipack_m1)[ivec_m1] :
                                                 cld_frac_r(icol,ipack)[ivec];
             }
           }
@@ -138,7 +145,9 @@ class P3Microphysics : public AtmosphereProcess
     view_2d_const pseudo_density;
     view_2d_const pseudo_density_dry;
     view_2d       T_atm;
-    view_2d_const cld_frac_t;
+    view_2d_const cld_frac_t_in;
+    view_2d_const cld_frac_l_in;
+    view_2d_const cld_frac_i_in;
     view_2d       qv;
     view_2d       qc;
     view_2d       nc;
@@ -162,7 +171,8 @@ class P3Microphysics : public AtmosphereProcess
            const view_2d_const& pmid_, const view_2d_const& pmid_dry_,
            const view_2d_const& pseudo_density_,
            const view_2d_const& pseudo_density_dry_, const view_2d& T_atm_,
-           const view_2d_const& cld_frac_t_, const view_2d& qv_, const view_2d& qc_,
+           const view_2d_const& cld_frac_t_in_, const view_2d_const& cld_frac_l_in_, const view_2d_const& cld_frac_i_in_,
+           const view_2d& qv_, const view_2d& qc_,
            const view_2d& nc_, const view_2d& qr_, const view_2d& nr_, const view_2d& qi_,
            const view_2d& qm_, const view_2d& ni_, const view_2d& bm_, const view_2d& qv_prev_,
            const view_2d& inv_exner_, const view_2d& th_atm_, const view_2d& cld_frac_l_,
@@ -178,7 +188,9 @@ class P3Microphysics : public AtmosphereProcess
       pseudo_density = pseudo_density_;
       pseudo_density_dry = pseudo_density_dry_;
       T_atm          = T_atm_;
-      cld_frac_t     = cld_frac_t_;
+      cld_frac_t_in     = cld_frac_t_in_;
+      cld_frac_l_in     = cld_frac_l_in_;
+      cld_frac_i_in     = cld_frac_i_in_;
       // OUT
       qv             = qv_;
       qc             = qc_;

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

@@ -65,7 +65,11 @@ ::back_to_cell_average(
   nc2ni_immers_freeze_tend.set(context, nc2ni_immers_freeze_tend * cld_frac_l);   // Number change associated with freexzing of cld drops
   nr_collect_tend.set(context, nr_collect_tend * ir_cldm);  // Rain number change due to collection from ice
   ni_selfcollect_tend.set(context, ni_selfcollect_tend * cld_frac_i);    // Ice self collection
-  qv2qi_vapdep_tend.set(context, qv2qi_vapdep_tend * cld_frac_i);    // Vapor deposition to ice phase
+  if ((cld_frac_i == cld_frac_l).all()) {
+    qv2qi_vapdep_tend.set(context, qv2qi_vapdep_tend * cld_frac_i);    // Vapor deposition to ice phase
+  } else {
+    qv2qi_vapdep_tend.set(context, qv2qi_vapdep_tend * (cld_frac_i-il_cldm));    // Vapor deposition to ice phase
+  }
   nr2ni_immers_freeze_tend.set(context, nr2ni_immers_freeze_tend * cld_frac_r);   // Change in number due to immersion freezing of rain
   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

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

@@ -7,6 +7,54 @@
 namespace scream {
 namespace p3 {
 
+template<typename S, typename D>
+KOKKOS_FUNCTION
+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 Spack& cld_frac_l, const Spack& cld_frac_i,
+  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
+  const auto sources = qc; // Source of cloud water
+  const auto il_cldm = min(cld_frac_i,cld_frac_l); 
+  Spack ratio;
+
+  constexpr Scalar qtendsmall = C::QTENDSMALL;
+  Smask enforce_conservation  = sinks > sources && sinks >= qtendsmall && context;  // determine if  conservation corrction is necessary
+  Smask nothing_todo = !enforce_conservation && context;
+
+  if (enforce_conservation.any()){
+    ratio.set(enforce_conservation, sources/sinks);
+    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);
+    qc2qr_ice_shed_tend.set(enforce_conservation, qc2qr_ice_shed_tend*ratio);
+    qc2qi_berg_tend.set(enforce_conservation, qc2qi_berg_tend*ratio);
+  }
+
+  if(nothing_todo.any()){
+    ratio.set(nothing_todo, 1); // If not limiting sinks on qc then most likely did not run out of qc
+  }
+
+  //PMC: ratio is also frac of step w/ liq. thus we apply qc2qi_berg_tend for
+  //"ratio" of timestep and vapor deposition and sublimation  for the
+  //remaining frac of the timestep.  Only limit if there will be cloud
+  // water to begin with.
+  // in instances where ratio < 1 and we have qc, qidep needs to take over
+   //but this is an in addition to the qidep we computed outside the mixed
+   //phase cloud. qidep*(1._rtype-ratio)*(il_cldm/(cld_frac_i-il_cldm)) is the additional
+   // vapor depositional growth rate that takes place within the mixed phase cloud
+   // after qc is depleted
+  enforce_conservation = sources > qtendsmall && context;
+  if (enforce_conservation.any()){
+    qv2qi_vapdep_tend.set(enforce_conservation, qv2qi_vapdep_tend + qv2qi_vapdep_tend*(1-ratio)*(il_cldm/(cld_frac_i-il_cldm)));
+    qi2qv_sublim_tend.set(enforce_conservation, qi2qv_sublim_tend*(1-ratio));
+  }
+}
+
 template<typename S, typename D>
 KOKKOS_FUNCTION
 void Functions<S,D>

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

@@ -105,6 +105,7 @@ ::p3_main_part2(
   constexpr Scalar latice       = C::LatIce;
 
   const bool do_ice_production = runtime_options.do_ice_production;
+  const bool use_separate_ice_liq_frac = runtime_options.use_separate_ice_liq_frac;
 
   team.team_barrier();
   hydrometeorsPresent = false;
@@ -411,10 +412,19 @@ ::p3_main_part2(
     // check qv because it is typically much greater than zero so seldom goes negative (and if it does
     // catastrophic failure is appropriate)]
 
-    // cloud
-    cloud_water_conservation(
-      qc(k), dt,
-      qc2qr_autoconv_tend, qc2qr_accret_tend, qc2qi_collect_tend, qc2qi_hetero_freeze_tend, qc2qr_ice_shed_tend, qc2qi_berg_tend, qi2qv_sublim_tend, qv2qi_vapdep_tend, not_skip_all);
+    if (use_separate_ice_liq_frac) {
+      // cloud
+      cloud_water_conservation(
+        qc(k), dt,
+        qc2qr_autoconv_tend, qc2qr_accret_tend, qc2qi_collect_tend, qc2qi_hetero_freeze_tend, qc2qr_ice_shed_tend, qc2qi_berg_tend, qi2qv_sublim_tend, qv2qi_vapdep_tend,
+        cld_frac_l(k), cld_frac_i(k), not_skip_all);
+    } else {
+      // cloud
+      cloud_water_conservation(
+        qc(k), dt,
+        qc2qr_autoconv_tend, qc2qr_accret_tend, qc2qi_collect_tend, qc2qi_hetero_freeze_tend, qc2qr_ice_shed_tend, qc2qi_berg_tend, qi2qv_sublim_tend, qv2qi_vapdep_tend,
+        not_skip_all);
+    }
 
     // rain
     rain_water_conservation(

components/eamxx/src/physics/p3/p3_functions.hpp

@@ -134,6 +134,7 @@ struct Functions
     bool set_cld_frac_l_to_one = false;
     bool set_cld_frac_i_to_one = false;
     bool set_cld_frac_r_to_one = false;
+    bool use_separate_ice_liq_frac = false;
 
     void load_runtime_options_from_file(ekat::ParameterList& params) {
       max_total_ni = params.get<double>("max_total_ni", max_total_ni);
@@ -159,6 +160,7 @@ struct Functions
       set_cld_frac_l_to_one = params.get<bool>("set_cld_frac_l_to_one", set_cld_frac_l_to_one);
       set_cld_frac_i_to_one = params.get<bool>("set_cld_frac_i_to_one", set_cld_frac_i_to_one);
       set_cld_frac_r_to_one = params.get<bool>("set_cld_frac_r_to_one", set_cld_frac_r_to_one);
+      use_separate_ice_liq_frac = params.get<bool>("use_separate_ice_liq_frac", use_separate_ice_liq_frac);
     }
 
   };
@@ -688,6 +690,13 @@ struct Functions
     const Int& kbot, const Int& ktop, const Int& kdir,
     bool& log_present);
 
+  KOKKOS_FUNCTION
+  static void 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 Spack& cld_frac_l, const Spack& cld_frac_i,
+    const Smask& context = Smask(true) );
+
   KOKKOS_FUNCTION
   static void 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,