Ice in kinematic 2D

models/kinematic_2D/cases/icmw8_case1.hpp

@@ -46,6 +46,7 @@ namespace config
     //aerosol chemical composition parameters (needed for activation)
     // for lgrngn:
     quantity<si::dimensionless, real_t> kappa; // CCN-derived value from Table 1 in Petters and Kreidenweis 2007
+    quantity<si::dimensionless, real_t> soluble_fraction; // soluble volume fraction
     // for blk_2m:
     quantity<si::dimensionless, real_t> chem_b; //ammonium sulphate //chem_b = 1.33; // sodium chloride
     // for lagrangian simulations with aq. chemistry

models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp

@@ -27,7 +27,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
   // member fields
   std::unique_ptr<libcloudphxx::lgrngn::particles_proto_t<real_t>> prtcls;
 
-  bool coal, sedi;
+  bool coal, sedi, ice_switch, ice_nucl, time_dep_ice_nucl, depo;
 
   // helper methods
   void diag()
@@ -92,6 +92,41 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
         rng_num++;
       }
     }
+    if (params.cloudph_opts_init.ice_switch)
+    {
+      prtcls->diag_ice_mix_ratio();
+      this->record_aux("ice_mix_ratio", prtcls->outbuf());
+      {
+        // ice_a
+        int rng_num = 0;
+        for (auto &rng_moms : params.out_ice)
+        {
+          auto &rng(rng_moms.first);
+          prtcls->diag_ice_a_rng(rng.first / si::metres, rng.second / si::metres);
+          for (auto &mom : rng_moms.second)
+          {
+            prtcls->diag_ice_a_mom(mom);
+            this->record_aux(aux_name("ice_a", rng_num, mom), prtcls->outbuf());
+          }
+          rng_num++;
+        }
+      }
+      {
+        // ice_c
+        int rng_num = 0;
+        for (auto &rng_moms : params.out_ice)
+        {
+          auto &rng(rng_moms.first);
+          prtcls->diag_ice_c_rng(rng.first / si::metres, rng.second / si::metres);
+          for (auto &mom : rng_moms.second)
+          {
+            prtcls->diag_ice_c_mom(mom);
+            this->record_aux(aux_name("ice_c", rng_num, mom), prtcls->outbuf());
+          }
+          rng_num++;
+        }
+      }
+    }
   } 
 
   libcloudphxx::lgrngn::arrinfo_t<real_t> make_arrinfo(
@@ -129,6 +164,10 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
   {
     coal = params.cloudph_opts.coal;
     sedi = params.cloudph_opts.sedi;
+    ice_switch = params.cloudph_opts_init.ice_switch;
+    ice_nucl = params.cloudph_opts.ice_nucl;
+    time_dep_ice_nucl = params.cloudph_opts_init.time_dep_ice_nucl;
+    depo = params.cloudph_opts.depo;
 
     parent_t::hook_ante_loop(nt); 
 
@@ -151,6 +190,11 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
       this->record_aux_const("n1_stp",   this->setup.n1_stp * si::cubic_metres);
       this->record_aux_const("n2_stp",   this->setup.n2_stp * si::cubic_metres);
       this->record_aux_const("kappa",    this->setup.kappa);
+      this->record_aux_const("soluble_fraction", this->setup.soluble_fraction);
+      this->record_aux_const("ice_switch", ice_switch);
+      this->record_aux_const("ice_nucl", ice_nucl);
+      this->record_aux_const("time_dep_ice_nucl", time_dep_ice_nucl);
+      this->record_aux_const("depo", depo);
 
       assert(params.backend != -1);
       assert(params.dt != 0); 
@@ -302,7 +346,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
     bool async = true;
     libcloudphxx::lgrngn::opts_t<real_t> cloudph_opts;
     libcloudphxx::lgrngn::opts_init_t<real_t> cloudph_opts_init;
-    outmom_t<real_t> out_dry, out_wet;
+    outmom_t<real_t> out_dry, out_wet, out_ice;
   };
 
   protected:

models/kinematic_2D/src/opts_common.hpp

@@ -60,6 +60,7 @@ void setopts_common(
     ("n1_stp",   po::value<real_t>()->default_value(60e6),     "n1_stp")
     ("n2_stp",   po::value<real_t>()->default_value(40e6),     "n2_stp")
     ("kappa",    po::value<real_t>()->default_value(.61),      "kappa")
+    ("soluble_fraction", po::value<real_t>()->default_value(1),"soluble_fraction")
     ("chem_b",   po::value<real_t>()->default_value(.55),      "chem_b")
     ("chem_rho", po::value<real_t>()->default_value(1.8e3),    "chem_rho")
     ("tau_rlx",  po::value<real_t>()->default_value(300),      "tau_rlx")
@@ -91,6 +92,7 @@ void setopts_common(
   setup.n1_stp   = vm["n1_stp"].as<real_t>() / si::cubic_metres;
   setup.n2_stp   = vm["n2_stp"].as<real_t>() / si::cubic_metres;
   setup.kappa    = vm["kappa"].as<real_t>();
+  setup.soluble_fraction = vm["soluble_fraction"].as<real_t>();
   setup.chem_b   = vm["chem_b"].as<real_t>();
   setup.chem_rho = vm["chem_rho"].as<real_t>() * si::kilograms / si::cubic_metres;
   setup.tau_rlx  = vm["tau_rlx"].as<real_t>() * si::seconds;
@@ -101,4 +103,7 @@ void setopts_common(
   setup.CO2_g_0  = vm["CO2_g_0"].as<real_t>();
   setup.HNO3_g_0 = vm["HNO3_g_0"].as<real_t>();
   setup.NH3_g_0  = vm["NH3_g_0"].as<real_t>();
+
+  if (setup.soluble_fraction < real_t(0) || setup.soluble_fraction > real_t(1))
+    throw std::runtime_error("Soluble fraction must be between [0,1]");
 }

models/kinematic_2D/src/opts_lgrngn.hpp

@@ -76,7 +76,9 @@ void parse_moms(
   std::vector<std::pair<std::string, std::string>> min_maxnum;
   outmom_t<config::real_t> &moms =
     opt == "out_dry" ? rt_params.out_dry :
-      rt_params.out_wet;
+      opt == "out_wet" ? rt_params.out_wet :
+        opt == "out_ice" ? rt_params.out_ice :
+          throw std::runtime_error("parse_moms(): unknown opt \"" + opt + "\"");
 
   const bool result = qi::phrase_parse(first, last, 
     *(
@@ -140,9 +142,10 @@ void parse_moms(
 
   outmom_t<config::real_t> &moms =
     opt == "out_dry" ? rt_params.out_dry :
-      opt == "out_wet" ? rt_params.out_wet : 
-         opt == "out_chem" ? rt_params.out_chem :
-           rt_params.out_wet_pH;
+      opt == "out_wet" ? rt_params.out_wet :
+        opt == "out_ice" ? rt_params.out_ice :
+          opt == "out_chem" ? rt_params.out_chem :
+            rt_params.out_wet_pH;
 
   const bool result = qi::phrase_parse(first, last, 
     *(
@@ -198,7 +201,7 @@ void setopts_micro_chem(
     {solver_t::ix::th, {"th", "[K]"}},
     {solver_t::ix::rv, {"rv", "[kg kg-1]"}}
   };
-  out_set = {"out_dry", "out_wet"};
+  out_set = {"out_dry", "out_wet", "out_ice"};
 }
 
 template <class solver_t, class ct_params_t>
@@ -217,7 +220,7 @@ void setopts_micro_chem(
     {solver_t::ix::NH3g,  {"NH3g", "[dimesnionless]"}},
     {solver_t::ix::HNO3g, {"HNO3g","[dimensionless]"}}
   };
-  out_set = {"out_dry", "out_wet", "out_chem", "out_wet_pH"};
+  out_set = {"out_dry", "out_wet", "out_ice", "out_chem", "out_wet_pH"};
 }
 
 // simulation and output parameters for micro=lgrngn
@@ -241,11 +244,16 @@ void setopts_micro(
     ("sedi", po::value<bool>()->default_value(rt_params.cloudph_opts.sedi) , "particle sedimentation (1=on, 0=off)")
     ("cond", po::value<bool>()->default_value(rt_params.cloudph_opts.cond) , "condensational growth  (1=on, 0=off)")
     ("coal", po::value<bool>()->default_value(rt_params.cloudph_opts.coal) , "collisional growth     (1=on, 0=off)")
+    ("coal_switch", po::value<bool>()->default_value(rt_params.cloudph_opts_init.coal_switch) , "enable collisions (1=on, 0=off), must be 0 when ice_switch=1")
     ("rcyc", po::value<bool>()->default_value(rt_params.cloudph_opts.rcyc) , "recycling of droplets  (1=on, 0=off)")
     ("chem_dsl", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_dsl) , "dissolving trace gases (1=on, 0=off)")
     ("chem_dsc", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_dsc) , "dissociation           (1=on, 0=off)")
     ("chem_rct", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry      (1=on, 0=off)")
     ("chem_switch", po::value<bool>()->default_value(rt_params.cloudph_opts_init.chem_switch) , "aqueous chemistry (1=on, 0=off)")
+    ("ice_switch", po::value<bool>()->default_value(rt_params.cloudph_opts_init.ice_switch) , "enable ice (1=on, 0=off)")
+    ("ice_nucl", po::value<bool>()->default_value(rt_params.cloudph_opts.ice_nucl) , "ice nucleation (1=on, 0=off)")
+    ("time_dep_ice_nucl", po::value<bool>()->default_value(rt_params.cloudph_opts_init.time_dep_ice_nucl) , "time dependent ice nucleation (1=on, 0=off)")
+    ("depo", po::value<bool>()->default_value(rt_params.cloudph_opts.depo) , "ice depositional growth (1=on, 0=off)")
     // free parameters
     ("sstp_cond", po::value<int>()->default_value(rt_params.cloudph_opts_init.sstp_cond), "no. of substeps for condensation")
     ("sstp_coal", po::value<int>()->default_value(rt_params.cloudph_opts_init.sstp_coal), "no. of substeps for coalescence")
@@ -255,6 +263,7 @@ void setopts_micro(
     // output
     ("out_dry", po::value<std::string>()->default_value("0:1|0"),       "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)")
     ("out_wet", po::value<std::string>()->default_value(".5e-6:25e-6|0,1,2,3;25e-6:1|0,3,6"),  "wet radius ranges and moment numbers (r1:r2|n1,n2...;...)")
+    ("out_ice", po::value<std::string>()->default_value(".1e-6:1|0,1,2,3"),  "ice semi-axis ranges and moment numbers (r1:r2|n1,n2...;...)")
     ("out_wet_pH", po::value<std::string>()->default_value("0:1|0"), "wet radius ranges for output of H+ and S_VI)")
     ("out_chem", po::value<std::string>()->default_value("0:1|0"),   "dry radius ranges for which chem mass is outputted")
     // collision and sedimentation
@@ -297,7 +306,7 @@ void setopts_micro(
 */
 
   rt_params.cloudph_opts_init.dry_distros.emplace(
-    libcloudphxx::lgrngn::kappa_soluble_fraction_t<config::real_t>{setup.kappa, config::real_t(1.)}, // kappa, soluble_fraction
+    libcloudphxx::lgrngn::kappa_soluble_fraction_t<config::real_t>{setup.kappa, setup.soluble_fraction},
     std::make_shared<config::log_dry_radii<config::real_t>> (setup)
   );
 
@@ -306,6 +315,27 @@ void setopts_micro(
   rt_params.cloudph_opts.sedi = vm["sedi"].as<bool>();
   rt_params.cloudph_opts.cond = vm["cond"].as<bool>();
   rt_params.cloudph_opts.coal = vm["coal"].as<bool>();
+  rt_params.cloudph_opts_init.coal_switch = vm["coal_switch"].as<bool>();
+  rt_params.cloudph_opts_init.ice_switch = vm["ice_switch"].as<bool>();
+  rt_params.cloudph_opts.ice_nucl = vm["ice_nucl"].as<bool>();
+  rt_params.cloudph_opts_init.time_dep_ice_nucl = vm["time_dep_ice_nucl"].as<bool>();
+  rt_params.cloudph_opts.depo = vm["depo"].as<bool>();
+
+  if (rt_params.cloudph_opts_init.ice_switch && rt_params.cloudph_opts_init.coal_switch)
+    throw std::runtime_error(
+      "Invalid options: ice_switch and coal_switch cannot be enabled at the same time"
+    );
+  if (!rt_params.cloudph_opts_init.ice_switch)
+  {
+    if (rt_params.cloudph_opts.ice_nucl)
+      throw std::runtime_error(
+        "Invalid options: ice nucleation requires ice_switch=1"
+      );
+    if (rt_params.cloudph_opts.depo)
+      throw std::runtime_error(
+        "Invalid options: deposition requires ice_switch=1"
+      );
+  }
 
   rt_params.cloudph_opts.rcyc = vm["rcyc"].as<bool>();
   rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as<bool>();