diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index 8bf6b59f9..fed44209c 100644 --- a/include/libcloudph++/lgrngn/opts.hpp +++ b/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 { diff --git a/models/kinematic_2D/cases/icmw8_case1.hpp b/models/kinematic_2D/cases/icmw8_case1.hpp index 3f44a4f79..76133a3fb 100644 --- a/models/kinematic_2D/cases/icmw8_case1.hpp +++ b/models/kinematic_2D/cases/icmw8_case1.hpp @@ -46,6 +46,7 @@ namespace config //aerosol chemical composition parameters (needed for activation) // for lgrngn: quantity kappa; // CCN-derived value from Table 1 in Petters and Kreidenweis 2007 + quantity rd_insol; // insoluble dry radius // for blk_2m: quantity chem_b; //ammonium sulphate //chem_b = 1.33; // sodium chloride // for lagrangian simulations with aq. chemistry diff --git a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp index 29c964985..039d2144c 100644 --- a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp +++ b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp @@ -27,7 +27,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common // member fields std::unique_ptr> 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 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 make_arrinfo( @@ -129,6 +164,10 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common { 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 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 bool async = true; libcloudphxx::lgrngn::opts_t cloudph_opts; libcloudphxx::lgrngn::opts_init_t cloudph_opts_init; - outmom_t out_dry, out_wet; + outmom_t out_dry, out_wet, out_ice; }; protected: diff --git a/models/kinematic_2D/src/opts_common.hpp b/models/kinematic_2D/src/opts_common.hpp index c07501357..09c288f06 100644 --- a/models/kinematic_2D/src/opts_common.hpp +++ b/models/kinematic_2D/src/opts_common.hpp @@ -60,6 +60,7 @@ void setopts_common( ("n1_stp", po::value()->default_value(60e6), "n1_stp") ("n2_stp", po::value()->default_value(40e6), "n2_stp") ("kappa", po::value()->default_value(.61), "kappa") + ("rd_insol", po::value()->default_value(0), "rd_insol") ("chem_b", po::value()->default_value(.55), "chem_b") ("chem_rho", po::value()->default_value(1.8e3), "chem_rho") ("tau_rlx", po::value()->default_value(300), "tau_rlx") @@ -91,6 +92,7 @@ void setopts_common( setup.n1_stp = vm["n1_stp"].as() / si::cubic_metres; setup.n2_stp = vm["n2_stp"].as() / si::cubic_metres; setup.kappa = vm["kappa"].as(); + setup.rd_insol = vm["rd_insol"].as() * si::metres; setup.chem_b = vm["chem_b"].as(); setup.chem_rho = vm["chem_rho"].as() * si::kilograms / si::cubic_metres; setup.tau_rlx = vm["tau_rlx"].as() * si::seconds; diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index 924921d24..f699ab78e 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -76,7 +76,8 @@ void parse_moms( std::vector> min_maxnum; outmom_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 &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 @@ -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()->default_value(rt_params.cloudph_opts.sedi) , "particle sedimentation (1=on, 0=off)") ("cond", po::value()->default_value(rt_params.cloudph_opts.cond) , "condensational growth (1=on, 0=off)") ("coal", po::value()->default_value(rt_params.cloudph_opts.coal) , "collisional growth (1=on, 0=off)") + ("coal_switch", po::value()->default_value(rt_params.cloudph_opts_init.coal_switch) , "enable collisions (1=on, 0=off)") ("rcyc", po::value()->default_value(rt_params.cloudph_opts.rcyc) , "recycling of droplets (1=on, 0=off)") ("chem_dsl", po::value()->default_value(rt_params.cloudph_opts.chem_dsl) , "dissolving trace gases (1=on, 0=off)") ("chem_dsc", po::value()->default_value(rt_params.cloudph_opts.chem_dsc) , "dissociation (1=on, 0=off)") ("chem_rct", po::value()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry (1=on, 0=off)") ("chem_switch", po::value()->default_value(rt_params.cloudph_opts_init.chem_switch) , "aqueous chemistry (1=on, 0=off)") + ("ice_switch", po::value()->default_value(rt_params.cloudph_opts_init.ice_switch) , "enable ice (1=on, 0=off)") + ("ice_nucl", po::value()->default_value(rt_params.cloudph_opts.ice_nucl) , "ice nucleation (1=on, 0=off)") + ("time_dep_ice_nucl", po::value()->default_value(rt_params.cloudph_opts_init.time_dep_ice_nucl) , "time dependent ice nucleation (1=on, 0=off)") + ("depo", po::value()->default_value(rt_params.cloudph_opts.depo) , "ice depositional growth (1=on, 0=off)") // free parameters ("sstp_cond", po::value()->default_value(rt_params.cloudph_opts_init.sstp_cond), "no. of substeps for condensation") ("sstp_coal", po::value()->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()->default_value("0:1|0"), "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)") ("out_wet", po::value()->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()->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()->default_value("0:1|0"), "wet radius ranges for output of H+ and S_VI)") ("out_chem", po::value()->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{setup.kappa, config::real_t(0.)}, // kappa, rd_insol + libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, setup.rd_insol / si::meters}, // kappa, rd_insol std::make_shared> (setup) ); @@ -306,6 +314,11 @@ void setopts_micro( rt_params.cloudph_opts.sedi = vm["sedi"].as(); rt_params.cloudph_opts.cond = vm["cond"].as(); rt_params.cloudph_opts.coal = vm["coal"].as(); + rt_params.cloudph_opts_init.coal_switch = vm["coal_switch"].as(); + rt_params.cloudph_opts_init.ice_switch = vm["ice_switch"].as(); + rt_params.cloudph_opts.ice_nucl = vm["ice_nucl"].as(); + rt_params.cloudph_opts_init.time_dep_ice_nucl = vm["time_dep_ice_nucl"].as(); + rt_params.cloudph_opts.depo = vm["depo"].as(); rt_params.cloudph_opts.rcyc = vm["rcyc"].as(); rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as();