tlocalMassBalance.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch) (HTM) git clone git://src.adamsgaard.dk/pism (DIR) Log (DIR) Files (DIR) Refs (DIR) LICENSE --- tlocalMassBalance.cc (17051B) --- 1 // Copyright (C) 2009, 2010, 2011, 2013, 2014, 2015, 2016, 2017, 2018 Ed Bueler and Constantine Khroulev and Andy Aschwanden 2 // 3 // This file is part of PISM. 4 // 5 // PISM is free software; you can redistribute it and/or modify it under the 6 // terms of the GNU General Public License as published by the Free Software 7 // Foundation; either version 3 of the License, or (at your option) any later 8 // version. 9 // 10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY 11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 13 // details. 14 // 15 // You should have received a copy of the GNU General Public License 16 // along with PISM; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 #include <cassert> 20 #include <ctime> // for time(), used to initialize random number gen 21 #include <gsl/gsl_rng.h> 22 #include <gsl/gsl_randist.h> 23 #include <gsl/gsl_math.h> // M_PI 24 #include <cmath> // for erfc() in CalovGreveIntegrand() 25 #include <algorithm> 26 27 #include "pism/util/pism_utilities.hh" 28 #include "pism/util/ConfigInterface.hh" 29 #include "localMassBalance.hh" 30 #include "pism/util/IceGrid.hh" 31 32 namespace pism { 33 namespace surface { 34 35 LocalMassBalance::Changes::Changes() { 36 firn_depth = 0.0; 37 snow_depth = 0.0; 38 melt = 0.0; 39 runoff = 0.0; 40 smb = 0.0; 41 } 42 43 LocalMassBalance::LocalMassBalance(Config::ConstPtr myconfig, units::System::Ptr system) 44 : m_config(myconfig), m_unit_system(system), 45 m_seconds_per_day(86400) { 46 // empty 47 } 48 49 LocalMassBalance::~LocalMassBalance() { 50 // empty 51 } 52 53 std::string LocalMassBalance::method() const { 54 return m_method; 55 } 56 57 PDDMassBalance::PDDMassBalance(Config::ConstPtr config, units::System::Ptr system) 58 : LocalMassBalance(config, system) { 59 precip_as_snow = m_config->get_flag("surface.pdd.interpret_precip_as_snow"); 60 Tmin = m_config->get_number("surface.pdd.air_temp_all_precip_as_snow"); 61 Tmax = m_config->get_number("surface.pdd.air_temp_all_precip_as_rain"); 62 pdd_threshold_temp = m_config->get_number("surface.pdd.positive_threshold_temp"); 63 refreeze_ice_melt = m_config->get_flag("surface.pdd.refreeze_ice_melt"); 64 65 m_method = "an expectation integral"; 66 } 67 68 69 /*! \brief Compute the number of points for temperature and 70 precipitation time-series. 71 */ 72 unsigned int PDDMassBalance::get_timeseries_length(double dt) { 73 const unsigned int NperYear = static_cast<unsigned int>(m_config->get_number("surface.pdd.max_evals_per_year")); 74 const double dt_years = units::convert(m_unit_system, dt, "seconds", "years"); 75 76 return std::max(1U, static_cast<unsigned int>(ceil(NperYear * dt_years))); 77 } 78 79 80 //! Compute the integrand in integral (6) in [\ref CalovGreve05]. 81 /*! 82 The integral is 83 \f[\mathrm{PDD} = \int_{t_0}^{t_0+\mathtt{dt}} dt\, 84 \bigg[\frac{\sigma}{\sqrt{2\pi}}\,\exp\left(-\frac{T_{ac}(t)^2}{2\sigma^2}\right) 85 + \frac{T_{ac}(t)}{2}\,\mathrm{erfc} 86 \left(-\frac{T_{ac}(t)}{\sqrt{2}\,\sigma}\right)\bigg] \f] 87 This procedure computes the quantity in square brackets. The value \f$T_{ac}(t)\f$ 88 in the above integral is in degrees C. Here we think of the argument `TacC` 89 as temperature in Celsius, but really it is the temperature above a threshold 90 at which it is "positive". 91 92 This integral is used for the expected number of positive degree days. The user can choose 93 \f$\sigma\f$ by option `-pdd_std_dev`. Note that the integral is over a time interval of 94 length `dt` instead of a whole year as stated in \ref CalovGreve05 . If `sigma` is zero, 95 return the positive part of `TacC`. 96 */ 97 double PDDMassBalance::CalovGreveIntegrand(double sigma, double TacC) { 98 99 if (sigma == 0) { 100 return std::max(TacC, 0.0); 101 } else { 102 const double Z = TacC / (sqrt(2.0) * sigma); 103 return (sigma / sqrt(2.0 * M_PI)) * exp(-Z*Z) + (TacC / 2.0) * erfc(-Z); 104 } 105 } 106 107 108 //! Compute the expected number of positive degree days from the input temperature time-series. 109 /** 110 * Use the rectangle method for simplicity. 111 * 112 * @param S standard deviation for air temperature excursions 113 * @param dt_series length of the step for the time-series 114 * @param T air temperature (array of length N) 115 * @param N length of the T array 116 * @param[out] PDDs pointer to a pre-allocated array with N-1 elements 117 */ 118 void PDDMassBalance::get_PDDs(double dt_series, 119 const std::vector<double> &S, 120 const std::vector<double> &T, 121 std::vector<double> &PDDs) { 122 assert(S.size() == T.size() and T.size() == PDDs.size()); 123 assert(dt_series > 0.0); 124 125 const double h_days = dt_series / m_seconds_per_day; 126 const size_t N = S.size(); 127 128 for (unsigned int k = 0; k < N; ++k) { 129 PDDs[k] = h_days * CalovGreveIntegrand(S[k], T[k] - pdd_threshold_temp); 130 } 131 } 132 133 134 //! \brief Extract snow accumulation from mixed (snow and rain) 135 //! precipitation using the temperature time-series. 136 /** Uses the temperature time-series to determine whether the 137 * precipitation is snow or rain. Rain is removed entirely from the 138 * surface mass balance, and will not be included in the computed 139 * runoff, which is meltwater runoff. There is an allowed linear 140 * transition for Tmin below which all precipitation is interpreted as 141 * snow, and Tmax above which all precipitation is rain (see, e.g. 142 * [\ref Hock2005b]). 143 * 144 * Sets P[i] to the *solid* (snow) accumulation *rate*. 145 * 146 * @param[in,out] P precipitation rate (array of length N) 147 * @param[in] T air temperature (array of length N) 148 * @param[in] N array length 149 */ 150 void PDDMassBalance::get_snow_accumulation(const std::vector<double> &T, 151 std::vector<double> &P) { 152 153 assert(T.size() == P.size()); 154 const size_t N = T.size(); 155 156 // Following \ref Hock2005b we employ a linear transition from Tmin to Tmax 157 for (unsigned int i = 0; i < N; i++) { 158 // do not allow negative precipitation 159 if (P[i] < 0.0) { 160 P[i] = 0.0; 161 continue; 162 } 163 164 if (precip_as_snow || T[i] <= Tmin) { // T <= Tmin, all precip is snow 165 // no change 166 } else if (T[i] < Tmax) { // linear transition from Tmin to Tmax 167 P[i] *= (Tmax - T[i]) / (Tmax - Tmin); 168 } else { // T >= Tmax, all precip is rain -- ignore it 169 P[i] = 0.0; 170 } 171 } 172 173 } 174 175 176 //! \brief Compute the surface mass balance at a location from the number of positive 177 //! degree days and the accumulation amount in a time interval. 178 /*! 179 * This is a PDD scheme. The input parameter `ddf.snow` is a rate of 180 * melting per positive degree day for snow. 181 * 182 * `accumulation` has units "meter / second". 183 * 184 * - a fraction of the melted snow and ice refreezes, conceptualized 185 * as superimposed ice, and this is controlled by parameter \c 186 * ddf.refreeze_fraction 187 * 188 * - the excess number of PDDs is used to melt both the ice that came 189 * from refreeze and then any ice which is already present. 190 * 191 * Ice melts at a constant rate per positive degree day, controlled by 192 * parameter `ddf.ice`. 193 * 194 * The scheme here came from EISMINT-Greenland [\ref RitzEISMINT], but 195 * is influenced by R. Hock (personal communication). 196 */ 197 PDDMassBalance::Changes PDDMassBalance::step(const DegreeDayFactors &ddf, 198 double PDDs, 199 double thickness, 200 double old_firn_depth, 201 double old_snow_depth, 202 double accumulation) { 203 double 204 firn_depth = old_firn_depth, 205 snow_depth = old_snow_depth, 206 max_snow_melted = PDDs * ddf.snow, 207 firn_melted = 0.0, 208 snow_melted = 0.0, 209 excess_pdds = 0.0; 210 211 assert(thickness >= 0); 212 213 // snow depth cannot exceed total thickness 214 snow_depth = std::min(snow_depth, thickness); 215 216 assert(snow_depth >= 0); 217 218 // firn depth cannot exceed thickness - snow_depth 219 firn_depth = std::min(firn_depth, thickness - snow_depth); 220 221 assert(firn_depth >= 0); 222 223 double ice_thickness = thickness - snow_depth - firn_depth; 224 225 assert(ice_thickness >= 0); 226 227 snow_depth += accumulation; 228 229 if (PDDs <= 0.0) { // The "no melt" case. 230 snow_melted = 0.0; 231 firn_melted = 0.0, 232 excess_pdds = 0.0; 233 } else if (max_snow_melted <= snow_depth) { 234 // Some of the snow melted and some is left; in any case, all of 235 // the energy available for melt, namely all of the positive 236 // degree days (PDDs) were used up in melting snow. 237 snow_melted = max_snow_melted; 238 firn_melted = 0.0; 239 excess_pdds = 0.0; 240 } else if (max_snow_melted <= firn_depth + snow_depth) { 241 // All of the snow is melted but some firn is left; in any case, all of 242 // the energy available for melt, namely all of the positive 243 // degree days (PDDs) were used up in melting snow and firn. 244 snow_melted = snow_depth; 245 firn_melted = max_snow_melted - snow_melted; 246 excess_pdds = 0.0; 247 } else { 248 // All (firn_depth and snow_depth meters) of snow and firn melted. Excess_pdds is the 249 // positive degree days available to melt ice. 250 firn_melted = firn_depth; 251 snow_melted = snow_depth; 252 excess_pdds = PDDs - ((firn_melted + snow_melted) / ddf.snow); // units: K day 253 } 254 255 double 256 ice_melted = std::min(excess_pdds * ddf.ice, ice_thickness), 257 melt = snow_melted + firn_melted + ice_melted, 258 ice_created_by_refreeze = 0.0; 259 260 if (refreeze_ice_melt) { 261 ice_created_by_refreeze = melt * ddf.refreeze_fraction; 262 } else { 263 // Should this only be snow melted? 264 ice_created_by_refreeze = (firn_melted + snow_melted) * ddf.refreeze_fraction; 265 } 266 267 const double runoff = melt - ice_created_by_refreeze; 268 269 snow_depth = std::max(snow_depth - snow_melted, 0.0); 270 firn_depth = std::max(firn_depth - firn_melted, 0.0); 271 272 // FIXME: need to add snow that hasn't melted, is this correct? 273 // firn_depth += (snow_depth - snow_melted); 274 // Turn firn into ice at X times accumulation 275 // firn_depth -= accumulation * m_config->get_number("surface.pdd.firn_compaction_to_accumulation_ratio"); 276 277 const double smb = accumulation - runoff; 278 279 Changes result; 280 // Ensure that we never generate negative ice thicknesses. As far as I can tell the code 281 // above guarantees that thickness + smb >= *in exact arithmetic*. The check below 282 // should make sure that we don't get bitten by rounding errors. 283 result.smb = thickness + smb >= 0 ? smb : -thickness; 284 result.firn_depth = firn_depth - old_firn_depth; 285 result.snow_depth = snow_depth - old_snow_depth; 286 result.melt = melt; 287 result.runoff = runoff; 288 289 assert(thickness + result.smb >= 0); 290 291 return result; 292 } 293 294 295 /*! 296 Initializes the random number generator (RNG). The RNG is GSL's recommended default, 297 which seems to be "mt19937" and is DIEHARD (whatever that means ...). Seed with 298 wall clock time in seconds in non-repeatable case, and with 0 in repeatable case. 299 */ 300 PDDrandMassBalance::PDDrandMassBalance(Config::ConstPtr config, units::System::Ptr system, 301 Kind kind) 302 : PDDMassBalance(config, system) { 303 pddRandGen = gsl_rng_alloc(gsl_rng_default); // so pddRandGen != NULL now 304 gsl_rng_set(pddRandGen, kind == REPEATABLE ? 0 : time(0)); 305 306 m_method = (kind == NOT_REPEATABLE 307 ? "simulation of a random process" 308 : "repeatable simulation of a random process"); 309 } 310 311 312 PDDrandMassBalance::~PDDrandMassBalance() { 313 if (pddRandGen != NULL) { 314 gsl_rng_free(pddRandGen); 315 pddRandGen = NULL; 316 } 317 } 318 319 320 /*! We need to compute simulated random temperature each actual \e 321 day, or at least as close as we can reasonably get. Output `N` is 322 number of days or number of days plus one. 323 324 Thus this method ignores 325 `config.get_number("surface.pdd.max_evals_per_year")`, which is 326 used in the base class PDDMassBalance. 327 328 Implementation of get_PDDs() requires returned N >= 2, so we 329 guarantee that. 330 */ 331 unsigned int PDDrandMassBalance::get_timeseries_length(double dt) { 332 return std::max(static_cast<size_t>(ceil(dt / m_seconds_per_day)), (size_t)2); 333 } 334 335 /** 336 * Computes 337 * \f[ 338 * \text{PDD} = \sum_{i=0}^{N-1} h_{\text{days}} \cdot \text{max}(T_i-T_{\text{threshold}}, 0). 339 * \f] 340 * 341 * @param S \f$\sigma\f$ (standard deviation for daily temperature excursions) 342 * @param dt_series time-series step, in seconds 343 * @param T air temperature 344 * @param N number of points in the temperature time-series, each corresponds to a sub-interval 345 * @param PDDs pointer to a pre-allocated array of length N 346 */ 347 void PDDrandMassBalance::get_PDDs(double dt_series, 348 const std::vector<double> &S, 349 const std::vector<double> &T, 350 std::vector<double> &PDDs) { 351 assert(S.size() == T.size() and T.size() == PDDs.size()); 352 assert(dt_series > 0.0); 353 354 const double h_days = dt_series / m_seconds_per_day; 355 const size_t N = S.size(); 356 357 for (unsigned int k = 0; k < N; ++k) { 358 // average temperature in k-th interval 359 double T_k = T[k] + gsl_ran_gaussian(pddRandGen, S[k]); // add random: N(0,sigma) 360 361 if (T_k > pdd_threshold_temp) { 362 PDDs[k] = h_days * (T_k - pdd_threshold_temp); 363 } 364 } 365 } 366 367 368 FaustoGrevePDDObject::FaustoGrevePDDObject(IceGrid::ConstPtr g) 369 : m_grid(g), m_config(g->ctx()->config()) { 370 371 m_beta_ice_w = m_config->get_number("surface.pdd.fausto.beta_ice_w"); 372 m_beta_snow_w = m_config->get_number("surface.pdd.fausto.beta_snow_w"); 373 374 m_T_c = m_config->get_number("surface.pdd.fausto.T_c"); 375 m_T_w = m_config->get_number("surface.pdd.fausto.T_w"); 376 m_beta_ice_c = m_config->get_number("surface.pdd.fausto.beta_ice_c"); 377 m_beta_snow_c = m_config->get_number("surface.pdd.fausto.beta_snow_c"); 378 379 m_fresh_water_density = m_config->get_number("constants.fresh_water.density"); 380 m_ice_density = m_config->get_number("constants.ice.density"); 381 m_pdd_fausto_latitude_beta_w = m_config->get_number("surface.pdd.fausto.latitude_beta_w"); 382 m_refreeze_fraction = m_config->get_number("surface.pdd.refreeze"); 383 384 385 m_temp_mj.create(m_grid, "temp_mj_faustogreve", WITHOUT_GHOSTS); 386 m_temp_mj.set_attrs("internal", 387 "mean July air temp from Fausto et al (2009) parameterization", 388 "K", "K", "", 0); 389 } 390 391 FaustoGrevePDDObject::~FaustoGrevePDDObject() { 392 // empty 393 } 394 395 LocalMassBalance::DegreeDayFactors FaustoGrevePDDObject::degree_day_factors(int i, int j, 396 double latitude) { 397 398 LocalMassBalance::DegreeDayFactors ddf; 399 ddf.refreeze_fraction = m_refreeze_fraction; 400 401 IceModelVec::AccessList list(m_temp_mj); 402 const double T_mj = m_temp_mj(i,j); 403 404 if (latitude < m_pdd_fausto_latitude_beta_w) { // case latitude < 72 deg N 405 ddf.ice = m_beta_ice_w; 406 ddf.snow = m_beta_snow_w; 407 } else { // case > 72 deg N 408 if (T_mj >= m_T_w) { 409 ddf.ice = m_beta_ice_w; 410 ddf.snow = m_beta_snow_w; 411 } else if (T_mj <= m_T_c) { 412 ddf.ice = m_beta_ice_c; 413 ddf.snow = m_beta_snow_c; 414 } else { // middle case T_c < T_mj < T_w 415 const double 416 lam_i = pow((m_T_w - T_mj) / (m_T_w - m_T_c) , 3.0), 417 lam_s = (T_mj - m_T_c) / (m_T_w - m_T_c); 418 ddf.ice = m_beta_ice_w + (m_beta_ice_c - m_beta_ice_w) * lam_i; 419 ddf.snow = m_beta_snow_w + (m_beta_snow_c - m_beta_snow_w) * lam_s; 420 } 421 } 422 423 // degree-day factors in \ref Faustoetal2009 are water-equivalent 424 // thickness per degree day; ice-equivalent thickness melted per degree 425 // day is slightly larger; for example, iwfactor = 1000/910 426 const double iwfactor = m_fresh_water_density / m_ice_density; 427 ddf.snow *= iwfactor; 428 ddf.ice *= iwfactor; 429 430 return ddf; 431 } 432 433 434 //! Updates mean July near-surface air temperature. 435 /*! 436 Unfortunately this duplicates code in SeaRISEGreenland::update(); 437 */ 438 void FaustoGrevePDDObject::update_temp_mj(const IceModelVec2S &surfelev, 439 const IceModelVec2S &lat, 440 const IceModelVec2S &lon) { 441 const double 442 d_mj = m_config->get_number("atmosphere.fausto_air_temp.d_mj"), // K 443 gamma_mj = m_config->get_number("atmosphere.fausto_air_temp.gamma_mj"), // K m-1 444 c_mj = m_config->get_number("atmosphere.fausto_air_temp.c_mj"), // K (degN)-1 445 kappa_mj = m_config->get_number("atmosphere.fausto_air_temp.kappa_mj"); // K (degW)-1 446 447 const IceModelVec2S 448 &h = surfelev, 449 &lat_degN = lat, 450 &lon_degE = lon; 451 452 IceModelVec::AccessList list{&h, &lat_degN, &lon_degE, &m_temp_mj}; 453 454 for (Points p(*m_grid); p; p.next()) { 455 const int i = p.i(), j = p.j(); 456 m_temp_mj(i,j) = d_mj + gamma_mj * h(i,j) + c_mj * lat_degN(i,j) + kappa_mj * (-lon_degE(i,j)); 457 } 458 } 459 460 } // end of namespace surface 461 } // end of namespace pism