From 12f67132b9f9b8083675dec923d766fdb103c590 Mon Sep 17 00:00:00 2001 From: Chris Marsh Date: Tue, 2 May 2023 13:04:15 -0600 Subject: [PATCH] Add eqn 14 and 13 from Essery 199 to pomli param. Don't remove out of triangles that did not saltate. --- src/modules/PBSM3D.cpp | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/src/modules/PBSM3D.cpp b/src/modules/PBSM3D.cpp index 23583219..eef014ec 100644 --- a/src/modules/PBSM3D.cpp +++ b/src/modules/PBSM3D.cpp @@ -845,12 +845,16 @@ void PBSM3D::run(mesh& domain) if (use_PomLi_probability) // Pomeroy and Li 2000 upscaled // probability { - // Essery, Li, and Pomeroy 1999 + // 1.Essery, R., Li, L. & Pomeroy, J. A distributed model of blowing snow over complex terrain. Hydrological Processes 13, 2423–2438 (1999). // Probability of blowing snow double A = (*face)["p_snow_hours"_s]; // hours since last snowfall double u_mean = 11.2 + 0.365 * T + 0.00706 * T * T + 0.9 * log(A); // eqn 10 T -> air temp, degC double delta = 0.145 * T + 0.00196 * T * T + 4.3; // eqn 11 - double Pu10 = 1.0 / (1.0 + exp((sqrt(M_PI) * (u_mean - u10)) / delta)); // eqn 12 + + double z0v = (d.N * d.dv * height_diff) / 2.0; // eqn 14 + double us = u10 / sqrt((1+340.0*z0v)); // eqn 13 + + double Pu10 = 1.0 / (1.0 + exp((sqrt(M_PI) * (u_mean - us)) / delta)); // eqn 12 (*face)["blowingsnow_probability"_s] = Pu10; // decrease the saltation by the probability amount @@ -1715,11 +1719,19 @@ void PBSM3D::run(mesh& domain) double swe = (*face)["swe"_s]; // mm --> kg/m^2 swe = is_nan(swe) ? 0 : swe; // handle the first timestep where swe won't have been // updated if we override the module order - if( mass < 0 && std::fabs(mass) > swe ) - { - (*face)["pbsm_more_than_avail"_s] = 1; - mass = -swe; - } + if( mass < 0 && std::fabs(mass) > swe ) + { + (*face)["pbsm_more_than_avail"_s] = 1; + mass = -swe; + } + + // if this triangle did not saltate, it is not a candidate for removal + auto& d = face->get_module_data(ID); + if (mass < 0 && !d.saltation) + { + mass = 0; + } + (*face)["drift_mass"_s] = mass; (*face)["sum_drift"_s] += mass; }