Ice in kinematic 2D

include/libcloudph++/lgrngn/opts.hpp

@@ -20,7 +20,7 @@ namespace libcloudphxx
     struct opts_t 
     {
       // process toggling
-      bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal, ice_nucl;
+      bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal, ice_nucl, depo;
 
       // RH limit for drop growth
       real_t RH_max;       
@@ -42,7 +42,7 @@ namespace libcloudphxx
       opts_t() : 
         adve(true), sedi(true), subs(false), cond(true), coal(true), src(false), rlx(false), rcyc(false),
         chem_dsl(false), chem_dsc(false), chem_rct(false),
-        turb_adve(false), turb_cond(false), turb_coal(false), ice_nucl(false),
+        turb_adve(false), turb_cond(false), turb_coal(false), ice_nucl(false), depo(false),
         RH_max(44), // :) (anything greater than 1.1 would be enough
         dt(-1) // negative means that we do not override dt in this step
       {

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::length, real_t> rd_insol; // insoluble dry radius
     // 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("rd_insol", this->setup.rd_insol / si::metres);
+      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")
+    ("rd_insol", po::value<real_t>()->default_value(0),        "rd_insol")
     ("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.rd_insol = vm["rd_insol"].as<real_t>() * si::metres;
   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;

models/kinematic_2D/src/opts_lgrngn.hpp

@@ -76,7 +76,8 @@ 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 :
+        rt_params.out_ice;
 
   const bool result = qi::phrase_parse(first, last, 
     *(
@@ -140,9 +141,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 +200,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 +219,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 +243,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)")
     ("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 +262,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 +305,7 @@ void setopts_micro(
 */
 
   rt_params.cloudph_opts_init.dry_distros.emplace(
-    libcloudphxx::lgrngn::kappa_rd_insol_t<config::real_t>{setup.kappa, config::real_t(0.)}, // kappa, rd_insol
+    libcloudphxx::lgrngn::kappa_rd_insol_t<config::real_t>{setup.kappa, setup.rd_insol / si::meters}, // kappa, rd_insol
     std::make_shared<config::log_dry_radii<config::real_t>> (setup)
   );
 
@@ -306,6 +314,11 @@ 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>();
 
   rt_params.cloudph_opts.rcyc = vm["rcyc"].as<bool>();
   rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as<bool>();