diff --git a/models/kinematic_2D/cases/icmw8_case1.hpp b/models/kinematic_2D/cases/icmw8_case1.hpp index 3f44a4f79..cb7e18b54 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 soluble_fraction; // soluble volume fraction // 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..16d9feeb8 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("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 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..1323ba202 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") + ("soluble_fraction", po::value()->default_value(1),"soluble_fraction") ("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.soluble_fraction = vm["soluble_fraction"].as(); 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; @@ -101,4 +103,7 @@ void setopts_common( setup.CO2_g_0 = vm["CO2_g_0"].as(); setup.HNO3_g_0 = vm["HNO3_g_0"].as(); setup.NH3_g_0 = vm["NH3_g_0"].as(); + + if (setup.soluble_fraction < real_t(0) || setup.soluble_fraction > real_t(1)) + throw std::runtime_error("Soluble fraction must be between [0,1]"); } diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index 2c5c37d59..3c701980d 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -76,7 +76,9 @@ 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 : + 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 &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 @@ -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()->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), must be 0 when ice_switch=1") ("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 +263,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 +306,7 @@ void setopts_micro( */ rt_params.cloudph_opts_init.dry_distros.emplace( - libcloudphxx::lgrngn::kappa_soluble_fraction_t{setup.kappa, config::real_t(1.)}, // kappa, soluble_fraction + libcloudphxx::lgrngn::kappa_soluble_fraction_t{setup.kappa, setup.soluble_fraction}, std::make_shared> (setup) ); @@ -306,6 +315,27 @@ 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(); + + 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(); rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as(); diff --git a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_const.h5 b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_const.h5 index cb2d45082..ccfe453e8 100644 Binary files a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_const.h5 and b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_const.h5 differ diff --git a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000000.h5 b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000000.h5 index 1c55b71b9..1bc149d3f 100644 Binary files a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000000.h5 and b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000000.h5 differ diff --git a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000020.h5 b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000020.h5 index 03eb66376..e5d505c8a 100644 Binary files a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000020.h5 and b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/refdata/travis_out_lgrngn/travis_timestep0000000020.h5 differ