Skip to content

Ice in super-droplets #196

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AgnieszkaMakulska wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Ice in super-droplets #196

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
552d25f
kappa_rd_insol
AgnieszkaMakulska Dec 8, 2025
04dc35d
ice opts and diag
AgnieszkaMakulska Dec 8, 2025
c787a0e
opts
AgnieszkaMakulska Dec 8, 2025
81d5f1d
ice spectrum
AgnieszkaMakulska Dec 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions include/UWLCM/output_bins.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,4 +24,15 @@ namespace
ret.push_back(setup::real_t(1e-6 * pow(10, -3 + i * .2)) * si::metres);
return ret;
}

std::vector<quantity<si::length, setup::real_t> > bins_ice()
{
std::vector<quantity<si::length, setup::real_t> > ret;

// ice polar and equatorial radius bins: .001 ... .01 ... 1 mm (30 bins in total)
for (int i = 0; i < 30; ++i)
ret.push_back(setup::real_t(1e-6 * pow(10, -3 + i * .2)) * si::metres);
return ret;
}

};
1 change: 1 addition & 0 deletions src/cases/CasesCommon.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@ namespace cases

// hygroscopicity kappa of the aerosol
quantity<si::dimensionless, real_t> kappa = .61; // defaults to ammonium sulphate; CCN-derived value from Table 1 in Petters and Kreidenweis 2007
quantity<si::length, real_t> rd_insol = 0 * si::metres;

real_t div_LS = 0.; // large-scale wind divergence (same as ForceParameters::D), 0. to turn off large-scale subsidence of SDs, TODO: add a process switch in libcloudph++ like for coal/cond/etc

Expand Down
2 changes: 1 addition & 1 deletion src/detail/user_params.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ struct user_params_t
quantity<si::length, setup::real_t> mean_rd1, mean_rd2;
quantity<si::dimensionless, setup::real_t> sdev_rd1, sdev_rd2;
quantity<power_typeof_helper<si::length, static_rational<-3>>::type, setup::real_t> n1_stp, n2_stp;
quantity<si::dimensionless, setup::real_t> kappa1, kappa2;
quantity<si::dimensionless, setup::real_t> kappa1, kappa2, rd_insol1, rd_insol2;
quantity<si::dimensionless, setup::real_t> case_n_stp_multiplier;

bool relax_th_rv,
Expand Down
75 changes: 49 additions & 26 deletions src/opts/opts_lgrngn.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,9 @@ void setopts_micro(
("chem_rct", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry (1=on, 0=off)")
("dev_count", po::value<int>()->default_value(0), "no. of CUDA devices")
("dev_id", po::value<int>()->default_value(-1), "CUDA backend - id of device to be used")
("ice_switch", po::value<bool>()->default_value(rt_params.cloudph_opts_init.ice_switch) , "enable ice microphysics (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)")
// free parameters
("exact_sstp_cond", po::value<bool>()->default_value(rt_params.cloudph_opts_init.exact_sstp_cond), "exact(per-particle) logic for substeps for condensation")
("diag_incloud_time", po::value<bool>()->default_value(rt_params.cloudph_opts_init.diag_incloud_time), "diagnose incloud time of droplets")
Expand All @@ -80,6 +83,7 @@ void setopts_micro(
//
("out_dry", po::value<std::string>()->default_value(""), "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)")
("out_wet", po::value<std::string>()->default_value(""), "wet radius ranges and moment numbers (r1:r2|n1,n2...;...)")
("out_ice", po::value<std::string>()->default_value(""), "ice radius (polar and equatorial) ranges and moment numbers (r1:r2|n1,n2...;...)")
("gccn", po::value<setup::real_t>()->default_value(0) , "concentration of giant aerosols = gccn * VOCALS observations")
// ("unit_test", po::value<bool>()->default_value(false) , "very low number concentration for unit tests")
("adve_scheme", po::value<std::string>()->default_value("euler") , "one of: euler, implicit, pred_corr")
Expand All @@ -89,6 +93,7 @@ void setopts_micro(
("ReL", po::value<setup::real_t>()->default_value(100) , "taylor-microscale reynolds number (onishi kernel)")
("out_dry_spec", po::value<bool>()->default_value(false), "enable output for plotting dry spectrum")
("out_wet_spec", po::value<bool>()->default_value(false), "enable output for plotting wet spectrum")
("out_ice_spec", po::value<bool>()->default_value(false), "enable output for plotting ice spectrum")
("rd_min", po::value<setup::real_t>()->default_value(rt_params.cloudph_opts_init.rd_min), "minimum dry radius of initialized droplets [m] (negative means automatic detection)")
("rd_max", po::value<setup::real_t>()->default_value(rt_params.cloudph_opts_init.rd_max), "maximum dry radius of initialized droplets [m] (negative means automatic detection); sd_conc_large_tail==true may result in initialization of even larger droplets")
("relax_ccn", po::value<bool>()->default_value(false) , "add CCN if per-level mean of CCN concentration is lower than (case-specific) desired concentration")
Expand Down Expand Up @@ -143,7 +148,7 @@ void setopts_micro(
}
if(user_params.n1_stp*si::cubic_metres >= 0) {
rt_params.cloudph_opts_init.dry_distros.emplace(
user_params.kappa1,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{user_params.kappa1, user_params.rd_insol1},
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
user_params.mean_rd1,
thrust_real_t(1.0e-6) * si::metres,
Expand All @@ -156,7 +161,7 @@ void setopts_micro(
}
if(user_params.n2_stp*si::cubic_metres >= 0) {
rt_params.cloudph_opts_init.dry_distros.emplace(
user_params.kappa2,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{user_params.kappa2, user_params.rd_insol2},
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
thrust_real_t(1.0e-6) * si::metres,
user_params.mean_rd2,
Expand All @@ -169,7 +174,7 @@ void setopts_micro(
}
if(user_params.n1_stp*si::cubic_metres < 0 && user_params.n2_stp*si::cubic_metres < 0) {
rt_params.cloudph_opts_init.dry_distros.emplace(
case_ptr->kappa,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{case_ptr->kappa, case_ptr->rd_insol},
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
case_ptr->mean_rd1,
case_ptr->mean_rd2,
Expand Down Expand Up @@ -235,7 +240,7 @@ void setopts_micro(

if(user_params.n1_stp*si::cubic_metres >= 0) {
rt_params.cloudph_opts_init.rlx_dry_distros.emplace(
user_params.kappa1,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{user_params.kappa1, user_params.rd_insol1},
std::make_tuple(
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
user_params.mean_rd1,
Expand All @@ -252,7 +257,7 @@ void setopts_micro(
}
if(user_params.n2_stp*si::cubic_metres >= 0) {
rt_params.cloudph_opts_init.rlx_dry_distros.emplace(
user_params.kappa2,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{user_params.kappa2, user_params.rd_insol2},
std::make_tuple(
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
thrust_real_t(1.0e-6) * si::metres,
Expand All @@ -270,7 +275,7 @@ void setopts_micro(

if(user_params.n1_stp*si::cubic_metres < 0 && user_params.n2_stp*si::cubic_metres < 0) {
rt_params.cloudph_opts_init.rlx_dry_distros.emplace(
case_ptr->kappa,
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{case_ptr->kappa, case_ptr->rd_insol},
std::make_tuple(
std::make_shared<setup::log_dry_radii<thrust_real_t>> (
case_ptr->mean_rd1,
Expand Down Expand Up @@ -305,11 +310,10 @@ void setopts_micro(
rt_params.cloudph_opts_init.src_z1 = case_ptr->gccn_max_height / si::meters;// 700;
// rt_params.cloudph_opts_init.src_z1 = 200;

rt_params.cloudph_opts_init.src_sd_conc = 38;

/*
rt_params.cloudph_opts_init.src_dry_sizes.emplace(
1.28, // kappa
rt_params.cloudph_opts.src_dry_sizes.emplace(
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{thrust_real_t(1.28), thrust_real_t(0.)},
std::map<setup::real_t, std::pair<setup::real_t, int> > {
{0.8e-6, {rt_params.gccn / rt_params.dt * 111800, 1}},
{1.0e-6, {rt_params.gccn / rt_params.dt * 68490, 1}},
Expand Down Expand Up @@ -353,20 +357,24 @@ void setopts_micro(
);
*/

rt_params.cloudph_opts_init.src_dry_distros.emplace(
1.28, // kappa
std::make_shared<setup::log_dry_radii_gccn<thrust_real_t>> (
log(0.8e-6), // minimum radius
log(10e-6), // maximum radius
rt_params.gccn / rt_params.dt // concenctration multiplier
)
rt_params.cloudph_opts.src_dry_distros.emplace(
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{thrust_real_t(1.28), thrust_real_t(0.)},
std::make_tuple(
std::make_shared<setup::log_dry_radii_gccn<thrust_real_t>> (
log(0.8e-6), // minimum radius
log(10e-6), // maximum radius
rt_params.gccn / rt_params.dt // concenctration multiplier
),
38, // number of SD created per cell (src_sd_conc)
1 // timestep interval at which source will be applied
)
);

// GCCN relaxation stuff
if(rt_params.user_params.relax_ccn)
{
rt_params.cloudph_opts_init.rlx_dry_distros.emplace(
1.28, // kappa
libcloudphxx::lgrngn::kappa_rd_insol_t<thrust_real_t>{thrust_real_t(1.28), thrust_real_t(0.)},
std::make_tuple(
std::make_shared<setup::log_dry_radii_gccn<thrust_real_t>> (
log(0.8e-6), // minimum radius
Expand Down Expand Up @@ -405,6 +413,7 @@ 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.ice_nucl = vm["ice_nucl"].as<bool>();

rt_params.cloudph_opts.rcyc = vm["rcyc"].as<bool>();
rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as<bool>();
Expand All @@ -425,6 +434,9 @@ void setopts_micro(
rt_params.cloudph_opts_init.rng_seed_init = user_params.rng_seed_init;
rt_params.cloudph_opts_init.rng_seed_init_switch = true;

rt_params.cloudph_opts_init.ice_switch = vm["ice_switch"].as<bool>();
rt_params.cloudph_opts_init.time_dep_ice_nucl = vm["time_dep_ice_nucl"].as<bool>();

// coalescence kernel choice
std::string kernel_str = vm["coal_kernel"].as<std::string>();

Expand Down Expand Up @@ -454,9 +466,9 @@ void setopts_micro(
rt_params.cloudph_opts_init.subs_switch = rt_params.subsidence == subs_t::local || rt_params.subsidence == subs_t::mean ? true : false;
rt_params.cloudph_opts.subs = rt_params.subsidence == subs_t::local || rt_params.subsidence == subs_t::mean ? true : false;

// parsing --out_dry and --out_wet options values
// parsing --out_dry and --out_wet --out_ice options values
// the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..."
for (auto &opt : std::set<std::string>({"out_dry", "out_wet"}))
for (auto &opt : std::set<std::string>({"out_dry", "out_wet", "out_ice"}))
{
namespace qi = boost::spirit::qi;
namespace phoenix = boost::phoenix;
Expand All @@ -466,10 +478,11 @@ void setopts_micro(
auto last = val.end();

std::vector<std::pair<std::string, std::string>> min_maxnum;
outmom_t<thrust_real_t> &moms =
opt == "out_dry"
? rt_params.out_dry
: rt_params.out_wet;
outmom_t<thrust_real_t>* moms_ptr = nullptr;
if (opt == "out_dry") moms_ptr = &rt_params.out_dry;
else if (opt == "out_wet") moms_ptr = &rt_params.out_wet;
else moms_ptr = &rt_params.out_ice;
outmom_t<thrust_real_t>& moms = *moms_ptr;

const bool result = qi::phrase_parse(first, last,
*(
Expand Down Expand Up @@ -515,12 +528,22 @@ void setopts_micro(
}
}

for (auto &opt : std::set<std::string>({"out_dry_spec", "out_wet_spec"}))
for (auto &opt : std::set<std::string>({"out_dry_spec", "out_wet_spec", "out_ice_spec"}))
{
if(vm[opt].as<bool>())
{
auto left_edges = opt == "out_dry_spec" ? bins_dry() : bins_wet();
auto &out = opt == "out_dry_spec" ? rt_params.out_dry : rt_params.out_wet;

std::vector<thrust_real_t> left_edges;
if (opt == "out_dry_spec") left_edges = bins_dry();
else if (opt == "out_wet_spec") left_edges = bins_wet();
else left_edges = bins_ice();

outmom_t<thrust_real_t>* out_ptr = nullptr;
if (opt == "out_dry_spec") out_ptr = &rt_params.out_dry;
else if (opt == "out_wet_spec") out_ptr = &rt_params.out_wet;
else out_ptr = &rt_params.out_ice;
auto &out = *out_ptr;

for (int i = 0; i < left_edges.size()-1; ++i)
{
out.push_back(outmom_t<thrust_real_t>::value_type(
Expand Down
59 changes: 58 additions & 1 deletion src/solvers/lgrngn/diag_lgrngn.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,12 @@ void slvr_lgrngn<ct_params_t>::diag()
*/

// recording precipitation rate per grid cel
prtcls->diag_all();
prtcls->diag_water();
prtcls->diag_precip_rate();
this->record_aux("precip_rate", prtcls->outbuf());
prtcls->diag_ice();
prtcls->diag_precip_rate_ice();
this->record_aux("precip_rate_ice_mass", prtcls->outbuf());

// recording 0th mom of rw of rd>=0.8um
// prtcls->diag_dry_rng(0.7999e-6, 1);
Expand Down Expand Up @@ -67,6 +70,8 @@ void slvr_lgrngn<ct_params_t>::diag()
// prtcls->diag_wet_mom(0);
// this->record_aux("non_gccn_rw_mom0", prtcls->outbuf());

prtcls->diag_water();

// recording 0th mom of rw of activated drops
prtcls->diag_rw_ge_rc();
prtcls->diag_wet_mom(0);
Expand Down Expand Up @@ -203,6 +208,26 @@ void slvr_lgrngn<ct_params_t>::diag()
// prtcls->diag_wet_mom(0);
// this->record_aux("bigrain_gccn_rw_mom0", prtcls->outbuf());

if (params.cloudph_opts_init.ice_switch)
{
prtcls->diag_ice();

// recording ice mixing ratio
prtcls->diag_ice_mass();
this->record_aux("r_i", prtcls->outbuf());

// recording 0th moment of ice
prtcls->diag_ice_a_mom(0);
this->record_aux("ice_mom0", prtcls->outbuf());

// recording 1st moment of ice_a
prtcls->diag_ice_a_mom(1);
this->record_aux("ice_a_mom1", prtcls->outbuf());

// recording 1st moment of ice_c
prtcls->diag_ice_c_mom(1);
this->record_aux("ice_c_mom1", prtcls->outbuf());
}
// recording requested statistical moments
if ((this->timestep ) % static_cast<int>(params.outfreq_spec) == 0)
{
Expand All @@ -221,6 +246,7 @@ void slvr_lgrngn<ct_params_t>::diag()
}

// wet
prtcls->diag_water();
rng_num = 0;
for (auto &rng_moms : params.out_wet)
{
Expand All @@ -233,6 +259,37 @@ void slvr_lgrngn<ct_params_t>::diag()
}
rng_num++;
}

if (params.cloudph_opts_init.ice_switch)
{
prtcls->diag_ice();
// ice_a
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
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++;
}
}
}

diag_prev_step = 1;
Expand Down
24 changes: 15 additions & 9 deletions src/solvers/lgrngn/hook_ante_loop_lgrngn.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,8 +42,12 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
n_sd_from_dry_sizes += rcm.second.second;

// src_dry_distros is used only in the first step to add GCCN below some level
int n_sd_from_src_dry_distros = params.cloudph_opts_init.src_sd_conc * params.cloudph_opts_init.src_z1 / params.cloudph_opts_init.z1 + 0.5;

//int n_sd_from_src_dry_distros = params.cloudph_opts.src_sd_conc * params.cloudph_opts_init.src_z1 / params.cloudph_opts_init.z1 + 0.5;
int n_sd_from_src_dry_distros = 0;
for (auto const& kv : params.cloudph_opts.src_dry_distros)
n_sd_from_src_dry_distros += std::get<1>(kv.second); // src_sd_conc
n_sd_from_src_dry_distros *= params.cloudph_opts_init.src_z1 / params.cloudph_opts_init.z1;

const int n_sd_per_cell = params.cloudph_opts_init.sd_conc + n_sd_from_dry_sizes + n_sd_from_src_dry_distros;

if(parent_t::n_dims == 2) // 2D
Expand Down Expand Up @@ -160,13 +164,16 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
this->record_aux_const("subs", "lgrngn", params.cloudph_opts.subs);
this->record_aux_const("cond", "lgrngn", params.cloudph_opts.cond);
this->record_aux_const("coal", "lgrngn", params.flag_coal); // cloudph_opts.coal could be 0 here due to spinup
this->record_aux_const("ice_nucl", "lgrngn", params.cloudph_opts.ice_nucl);
this->record_aux_const("rcyc", "lgrngn", params.cloudph_opts.rcyc);
this->record_aux_const("src", "lgrngn", params.cloudph_opts.src);
this->record_aux_const("rlx", "lgrngn", params.cloudph_opts.rlx);
this->record_aux_const("gccn", "lgrngn", params.gccn);
this->record_aux_const("turb_adve", "lgrngn", params.cloudph_opts.turb_adve);
this->record_aux_const("turb_cond", "lgrngn", params.cloudph_opts.turb_cond);
this->record_aux_const("turb_coal", "lgrngn", params.cloudph_opts.turb_coal);
this->record_aux_const("turb_coal", "lgrngn", params.cloudph_opts.turb_coal);
this->record_aux_const("ice_switch", "lgrngn", params.cloudph_opts_init.ice_switch);
this->record_aux_const("time_dep_ice_nucl", "lgrngn", params.cloudph_opts_init.time_dep_ice_nucl);
this->record_aux_const("chem_switch", "lgrngn", params.cloudph_opts_init.chem_switch);
this->record_aux_const("coal_switch", "lgrngn", params.cloudph_opts_init.coal_switch);
this->record_aux_const("sedi_switch", "lgrngn", params.cloudph_opts_init.sedi_switch);
Expand All @@ -179,10 +186,8 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
this->record_aux_const("chem_dsc", "lgrngn", params.cloudph_opts.chem_dsc);
this->record_aux_const("chem_rct", "lgrngn", params.cloudph_opts.chem_rct);
this->record_aux_const("chem_rho", "lgrngn", params.cloudph_opts_init.chem_rho);
this->record_aux_const("opts_init RH_max", "lgrngn", params.cloudph_opts_init.RH_max);
this->record_aux_const("supstp_src", "lgrngn", params.cloudph_opts_init.supstp_src);
this->record_aux_const("supstp_rlx", "lgrngn", params.cloudph_opts_init.supstp_rlx);
this->record_aux_const("src_sd_conc", "lgrngn", params.cloudph_opts_init.src_sd_conc);
this->record_aux_const("opts_init RH_max", "lgrngn", params.cloudph_opts_init.RH_max);
this->record_aux_const("supstp_rlx", "lgrngn", params.cloudph_opts_init.supstp_rlx);
this->record_aux_const("src_z1", "lgrngn", params.cloudph_opts_init.src_z1);
this->record_aux_const("rlx_bins", "lgrngn", params.cloudph_opts_init.rlx_bins);
this->record_aux_const("rlx_sd_per_bin", "lgrngn", params.cloudph_opts_init.rlx_sd_per_bin);
Expand All @@ -201,8 +206,9 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
else
this->record_aux_const("aerosol_conc_factor", "uniform");
this->record_aux_const("outfreq_spec", "lgrngn", params.outfreq_spec);
// store left and right edges of bins for wet and dry radii moments
for (auto &out : std::set<std::pair<outmom_t<real_t>, std::string>>({{params.out_dry, "out_dry"}, {params.out_wet, "out_wet"}}))
// store left and right edges of bins for wet, dry and ice radii moments
for (auto &out : std::set<std::pair<outmom_t<real_t>, std::string>>({{params.out_dry, "out_dry"}, {params.out_wet, "out_wet"},
{params.out_ice, "out_ice_a"}, {params.out_ice, "out_ice_c"}}))
{
if(!out.first.empty())
{
Expand Down
Loading
Loading