Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DEBMSimple.cc
Go to the documentation of this file.
1// Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2022, 2023, 2024, 2025 PISM Authors
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 <algorithm> // std::min
20
21#include "pism/coupler/surface/DEBMSimple.hh"
22
23#include "pism/coupler/AtmosphereModel.hh"
24#include "pism/util/Grid.hh"
25#include "pism/util/Mask.hh"
26#include "pism/util/Time.hh"
27#include "pism/util/io/File.hh"
28
29#include "pism/coupler/util/options.hh"
30#include "pism/geometry/Geometry.hh"
31#include "pism/util/array/CellType.hh"
32#include "pism/util/error_handling.hh"
33#include "pism/util/pism_utilities.hh"
34#include "pism/util/Vars.hh"
35#include "pism/util/array/Forcing.hh"
36
37namespace pism {
38namespace surface {
39
40///// PISM surface model implementing a dEBM-Simple scheme.
41
42DEBMSimple::DEBMSimple(std::shared_ptr<const Grid> g, std::shared_ptr<atmosphere::AtmosphereModel> input)
43 : SurfaceModel(g, std::move(input)),
44 m_model(*g->ctx()),
45 m_mass_flux(m_grid, "climatic_mass_balance"),
46 m_snow_depth(m_grid, "snow_depth"),
47 m_temperature_driven_melt(m_grid, "debm_temperature_driven_melt_flux"),
48 m_insolation_driven_melt(m_grid, "debm_insolation_driven_melt_flux"),
49 m_offset_melt(m_grid, "debm_offset_melt_flux"),
50 m_surface_albedo(m_grid, "surface_albedo"),
51 m_transmissivity(m_grid, "atmosphere_transmissivity") {
52
53 m_sd_use_param = m_config->get_flag("surface.debm_simple.std_dev.param.enabled");
54 m_sd_param_a = m_config->get_number("surface.debm_simple.std_dev.param.a");
55 m_sd_param_b = m_config->get_number("surface.debm_simple.std_dev.param.b");
56
57 m_precip_as_snow = m_config->get_flag("surface.debm_simple.interpret_precip_as_snow");
58 m_Tmax = m_config->get_number("surface.debm_simple.air_temp_all_precip_as_rain");
59 m_Tmin = m_config->get_number("surface.debm_simple.air_temp_all_precip_as_snow");
60
61 if (m_Tmax <= m_Tmin) {
63 "surface.debm_simple.air_temp_all_precip_as_rain has to exceed "
64 "surface.debm_simple.air_temp_all_precip_as_snow");
65 }
66
67 // note: does not need to be calendar-aware
68 m_year_length = units::convert(g->ctx()->unit_system(), 1.0, "years", "seconds");
69
70 m_n_per_year = static_cast<unsigned int>(m_config->get_number("surface.debm_simple.max_evals_per_year"));
71
72 auto albedo_input = m_config->get_string("surface.debm_simple.albedo_input.file");
73 if (not albedo_input.empty()) {
74
75 m_log->message(2, " Surface albedo is read in from %s...", albedo_input.c_str());
76
77 File file(m_grid->com, albedo_input, io::PISM_GUESS, io::PISM_READONLY);
78
79 int buffer_size = static_cast<int>(m_config->get_number("input.forcing.buffer_size"));
80 bool periodic = m_config->get_flag("surface.debm_simple.albedo_input.periodic");
81 m_input_albedo = std::make_shared<array::Forcing>(m_grid,
82 file,
83 "surface_albedo",
84 "surface_albedo",
85 buffer_size,
86 periodic,
87 LINEAR);
88 } else {
89 m_input_albedo = nullptr;
90 }
91
92 // initialize the spatially-variable air temperature standard deviation
93
94 auto air_temp_sd = m_config->get_string("surface.debm_simple.std_dev.file");
95 if (not air_temp_sd.empty()) {
96 m_log->message(2, " Reading standard deviation of near-surface air temperature from '%s'...\n",
97 air_temp_sd.c_str());
98
99 int buffer_size = static_cast<int>(m_config->get_number("input.forcing.buffer_size"));
100
102
103 bool periodic = m_config->get_flag("surface.debm_simple.std_dev.periodic");
104 m_air_temp_sd = std::make_shared<array::Forcing>(m_grid, file, "air_temp_sd",
105 "", // no standard name
106 buffer_size, periodic, LINEAR);
107 } else {
108 double temp_std_dev = m_config->get_number("surface.debm_simple.std_dev");
109
110 m_air_temp_sd = array::Forcing::Constant(m_grid, "air_temp_sd", temp_std_dev);
111 m_log->message(2, " Using constant standard deviation of near-surface air temperature.\n");
112 }
113
114 m_air_temp_sd->metadata(0)
115 .long_name("standard deviation of near-surface air temperature")
116 .units("kelvin");
117
119 .long_name("instantaneous surface mass balance (accumulation/ablation) rate")
120 .units("kg m^-2 s^-1")
121 .standard_name("land_ice_surface_specific_mass_balance_flux");
122 m_mass_flux.metadata().set_string("comment", "positive values correspond to ice gain");
123
124 // diagnostic fields:
125
126 {
131
133 .long_name("temperature-driven melt in dEBM-simple")
134 .units("kg m^-2");
136
138 .long_name("insolation-driven melt in dEBM-simple")
139 .units("kg m^-2");
141
143 .long_name("offset melt in dEBM-simple")
144 .units("kg m^-2");
145 m_offset_melt.set(0.0);
146 }
147
149 .long_name("snow cover depth (set to zero once a year)")
150 .units("m");
151 m_snow_depth.set(0.0);
152
154 .long_name("surface_albedo")
155 .units("1")
156 .standard_name("surface_albedo");
158
160 .long_name("atmosphere_transmissivity")
161 .units("1");
162
164}
165
166void DEBMSimple::init_impl(const Geometry &geometry) {
167
168 // call the default implementation (not the interface method init())
169 SurfaceModel::init_impl(geometry);
170
171 {
172 m_log->message(2,
173 "* Initializing dEBM-simple, the diurnal Energy Balance Model (simple version).\n"
174 " Inputs: precipitation and 2m air temperature from an atmosphere model.\n"
175 " Outputs: SMB and ice upper surface temperature.\n");
176 }
177
178 // initializing the model state
180
181 auto default_albedo = m_config->get_number("surface.debm_simple.albedo_max");
182 if (input.type == INIT_RESTART) {
183 m_snow_depth.read(input.filename, input.record);
184 m_surface_albedo.read(input.filename, input.record);
185 } else if (input.type == INIT_BOOTSTRAP) {
187 m_surface_albedo.regrid(input.filename, io::Default(default_albedo));
188 } else {
189 m_snow_depth.set(0.0);
190 m_surface_albedo.set(default_albedo);
191 }
192
193 {
194 regrid("dEBM-simple surface model", m_snow_depth);
195 regrid("dEBM-simple surface model", m_surface_albedo);
196 }
197
198 if ((bool)m_input_albedo) {
199 auto filename = m_config->get_string("surface.debm_simple.albedo_input.file");
200 bool periodic = m_config->get_flag("surface.debm_simple.albedo_input.periodic");
201 m_input_albedo->init(filename, periodic);
202 }
203
204 // finish up
205 {
207
208 m_accumulation->set(0.0);
209 m_melt->set(0.0);
210 m_runoff->set(0.0);
211 }
212}
213
215 return m_atmosphere->max_timestep(my_t);
216}
217
219 // compute the time corresponding to the beginning of the next balance year
220 double balance_year_start_day = m_config->get_number("surface.mass_balance_year_start_day"),
221 one_day = units::convert(m_sys, 1.0, "days", "seconds"),
222 year_start = this->time().calendar_year_start(time),
223 balance_year_start = year_start + (balance_year_start_day - 1.0) * one_day;
224
225 if (balance_year_start > time) {
226 return balance_year_start;
227 }
228 return this->time().increment_date(balance_year_start, 1);
229}
230
231/** @brief Extracts snow accumulation from mixed (snow and rain) precipitation using a
232 * temperature threshold with a linear transition.
233 *
234 * Rain is removed entirely from the surface mass balance, and will not be included in
235 * the computed runoff, which is meltwater runoff.
236 *
237 * There is an linear transition for Tmin below which all precipitation is interpreted as
238 * snow, and Tmax above which all precipitation is rain (see, e.g. [\ref Hock2005b]).
239 *
240 * Returns the *solid* (snow) accumulation *rate*.
241 *
242 * @param[in] T air temperature
243 * @param[in] P precipitation rate
244 */
245double DEBMSimple::snow_accumulation(double T, double P) const {
246
247 // do not allow negative precipitation
248 if (P < 0.0) {
249 return 0.0;
250 }
251
252 if (m_precip_as_snow or T <= m_Tmin) {
253 // T <= Tmin, all precip is snow
254 return P;
255 }
256
257 if (T < m_Tmax) { // linear transition from Tmin to Tmax
258 return P * (m_Tmax - T) / (m_Tmax - m_Tmin);
259 }
260
261 // T >= Tmax, all precip is rain -- ignore it
262 return 0.0;
263}
264
265
266void DEBMSimple::update_impl(const Geometry &geometry, double t, double dt) {
267
268 const double melting_point = 273.15;
269
270 // update to ensure that temperature and precipitation time series are correct:
271 m_atmosphere->update(geometry, t, dt);
272
273 // Use near-surface air temperature as the top-of-the-ice temperature:
274 m_temperature->copy_from(m_atmosphere->air_temperature());
275
276 // Set up air temperature and precipitation time series
277 int N = static_cast<int>(timeseries_length(dt));
278
279 const double dtseries = dt / N;
280 std::vector<double> ts(N), T(N), S(N), P(N), Alb(N);
281 std::vector<DEBMSimpleOrbitalParameters> orbital(N);
282
283 for (int k = 0; k < N; ++k) {
284 ts[k] = t + k * dtseries;
285
286 // pre-compute orbital parameters which depend on time and *not* on the map-plane
287 // location
288 orbital[k] = m_model.orbital_parameters(ts[k]);
289 }
290
291 // update standard deviation time series
292 m_air_temp_sd->update(t, dt);
293 m_air_temp_sd->init_interpolation(ts);
294
295 const auto &mask = geometry.cell_type;
296 const auto &H = geometry.ice_thickness;
297 const auto &surface_altitude = geometry.ice_surface_elevation;
298
300 { &mask,
301 &H,
302 m_air_temp_sd.get(),
305 m_accumulation.get(),
306 m_melt.get(),
307 m_runoff.get(),
313 &geometry.latitude,
314 &surface_altitude
315 };
316
317 if ((bool)m_input_albedo) {
318 m_input_albedo->update(t, dt);
319 m_input_albedo->init_interpolation(ts);
320 list.add(*m_input_albedo);
321 }
322
323 double
324 ice_density = m_config->get_number("constants.ice.density"),
325 sigmalapserate = m_config->get_number("surface.pdd.std_dev.lapse_lat_rate"),
326 sigmabaselat = m_config->get_number("surface.pdd.std_dev.lapse_lat_base");
327
328 m_atmosphere->init_timeseries(ts);
329 m_atmosphere->begin_pointwise_access();
330
331 ParallelSection loop(m_grid->com);
332 try {
333 for (auto p = m_grid->points(); p; p.next()) {
334 const int i = p.i(), j = p.j();
335
336 double latitude = geometry.latitude(i, j);
337
338 // Get the temperature time series from an atmosphere model and its modifiers
339 m_atmosphere->temp_time_series(i, j, T);
340
341 if (mask.ice_free_ocean(i, j)) {
342 // ignore precipitation over ice-free ocean
343 for (int k = 0; k < N; ++k) {
344 P[k] = 0.0;
345 }
346 } else {
347 // elsewhere, get precipitation from the atmosphere model
348 m_atmosphere->precip_time_series(i, j, P);
349
350 // Use temperature time series to remove rainfall from precipitation and convert to
351 // m/s ice equivalent.
352 for (int k = 0; k < N; ++k) {
353 P[k] = snow_accumulation(T[k], // air temperature (input)
354 P[k] / ice_density); // precipitation rate (input, gets overwritten)
355 }
356 }
357
358 if ((bool)m_input_albedo) {
359 m_input_albedo->interp(i, j, Alb);
360 }
361
362 // standard deviation of daily variability of air temperature
363 {
364 // interpolate temperature standard deviation time series
365 //
366 // Note: this works when m_air_temp_sd is constant in time.
367 m_air_temp_sd->interp(i, j, S);
368
369 if (sigmalapserate != 0.0) {
370 // apply standard deviation lapse rate on top of prescribed values
371 for (int k = 0; k < N; ++k) {
372 S[k] += sigmalapserate * (latitude - sigmabaselat);
373 }
374 (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
375 } else if (m_sd_use_param and mask.icy(i, j)) {
376 // apply standard deviation parameterization over ice if in use
377 for (int k = 0; k < N; ++k) {
378 S[k] = std::max(m_sd_param_a * (T[k] - melting_point) + m_sd_param_b, 0.0);
379 }
380 (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
381 }
382 }
383
384 {
385 double next_snow_depth_reset = m_next_balance_year_start;
386
387 // make copies of firn and snow depth values at this point to avoid accessing 2D
388 // fields in the inner loop
389 double
390 ice_thickness = H(i, j),
391 snow = m_snow_depth(i, j),
392 surfelev = surface_altitude(i, j),
393 albedo = m_surface_albedo(i, j);
394
395 auto cell_type = static_cast<MaskValue>(mask.as_int(i, j));
396
397 double
398 A = 0.0, // accumulation
399 M = 0.0, // melt
400 R = 0.0, // runoff
401 SMB = 0.0, // resulting mass balance
402 Mi = 0.0, // insolation melt contribution
403 Mt = 0.0, // temperature melt contribution
404 Mc = 0.0, // offset melt contribution
405 Al = 0.0; // albedo
406
407 // beginning of the loop over small time steps:
408 for (int k = 0; k < N; ++k) {
409
410 if (ts[k] >= next_snow_depth_reset) {
411 snow = 0.0;
412 while (next_snow_depth_reset <= ts[k]) {
413 next_snow_depth_reset = time().increment_date(next_snow_depth_reset, 1);
414 }
415 }
416
417 auto accumulation = P[k] * dtseries;
418
419 DEBMSimpleMelt melt_info{};
420 if (not mask::ice_free_ocean(cell_type)) {
421
422 melt_info = m_model.melt(orbital[k].solar_declination,
423 orbital[k].distance_factor,
424 dtseries,
425 S[k],
426 T[k],
427 surfelev,
428 latitude,
429 (bool)m_input_albedo ? Alb[k] : albedo);
430 }
431
432 auto changes = m_model.step(ice_thickness,
433 melt_info.total_melt,
434 snow,
436
437 if ((bool) m_input_albedo) {
438 albedo = Alb[k];
439 } else {
440 albedo = m_model.albedo(changes.melt / dtseries, cell_type);
441 }
442
443 // update ice thickness
444 ice_thickness += changes.smb;
445 assert(ice_thickness >= 0);
446 // update snow depth
447 snow += changes.snow_depth;
448 assert(snow >= 0);
449 // update total accumulation, melt, and runoff
450 {
451 A += accumulation;
452 M += changes.melt;
453 Mt += melt_info.temperature_melt;
454 Mi += melt_info.insolation_melt;
455 Mc += melt_info.offset_melt;
456 R += changes.runoff;
457 SMB += changes.smb;
458 Al += albedo;
459 }
460 } // end of the time-stepping loop
461
462 // set firn and snow depths
463 m_snow_depth(i, j) = snow;
464 m_surface_albedo(i, j) = Al / N;
466
467 // set melt terms at this point, converting
468 // from "meters, ice equivalent" to "kg / m^2"
469 m_temperature_driven_melt(i, j) = Mt * ice_density;
470 m_insolation_driven_melt(i, j) = Mi * ice_density;
471 m_offset_melt(i, j) = Mc * ice_density;
472
473 // set total accumulation, melt, and runoff, and SMB at this point, converting
474 // from "meters, ice equivalent" to "kg / m^2"
475 {
476 (*m_accumulation)(i, j) = A * ice_density;
477 (*m_melt)(i, j) = M * ice_density;
478 (*m_runoff)(i, j) = R * ice_density;
479 // m_mass_flux (unlike m_accumulation, m_melt, and m_runoff), is a
480 // rate. m * (kg / m^3) / second = kg / m^2 / second
481 m_mass_flux(i, j) = SMB * ice_density / dt;
482 }
483 }
484
485 if (mask.ice_free_ocean(i, j)) {
486 m_snow_depth(i, j) = 0.0; // snow over the ocean does not stick
487 }
488 }
489 } catch (...) {
490 loop.failed();
491 }
492 loop.check();
493
494 m_atmosphere->end_pointwise_access();
495
498}
499
500
502 return m_mass_flux;
503}
504
508
512
514 return *m_melt;
515}
516
518 return *m_runoff;
519}
520
522 return m_snow_depth;
523}
524
526 return *m_air_temp_sd;
527}
528
532
536
538 return m_offset_melt;
539}
540
544
548
554
555void DEBMSimple::write_model_state_impl(const File &output) const {
557 m_snow_depth.write(output);
558 m_surface_albedo.write(output);
559}
560
564
565namespace diagnostics {
566
567/*! @brief Report mean top of atmosphere insolation */
568class DEBMSInsolation : public Diag<DEBMSimple>
569{
570public:
572 m_vars = { { m_sys, "insolation" } };
573 m_vars[0]
574 .long_name(
575 "mean top of atmosphere insolation during the period when the sun is above the critical angle Phi")
576 .units("W m^-2");
577 }
578
579protected:
580 std::shared_ptr<array::Array> compute_impl() const {
581
582 auto result = allocate<array::Scalar>("insolation");
583
584 const auto *latitude = m_grid->variables().get_2d_scalar("latitude");
585 auto ctx = m_grid->ctx();
586
587 {
588 const auto& M = model->pointwise_model();
589
590 auto orbital = M.orbital_parameters(ctx->time()->current());
591
592 array::AccessScope list{latitude, result.get()};
593
594 for (auto p = m_grid->points(); p; p.next()) {
595 const int i = p.i(), j = p.j();
596
597 (*result)(i, j) = M.insolation_diagnostic(orbital.solar_declination,
598 orbital.distance_factor,
599 (*latitude)(i, j));
600 }
601 }
602
603 return result;
604 }
605};
606
608
609/*! @brief Report surface insolation melt, averaged over the reporting interval */
610class DEBMSInsolationMelt : public DiagAverageRate<DEBMSimple> {
611public:
614 kind == AMOUNT
615 ? "debm_insolation_driven_melt_flux"
616 : "debm_insolation_driven_melt_rate",
618 m_kind(kind),
619 m_melt_mass(m_grid, "debm_insolation_driven_melt_mass") {
620
621 std::string
622 name = "debm_insolation_driven_melt_flux",
623 long_name = "surface insolation melt, averaged over the reporting interval",
624 accumulator_units = "kg m^-2",
625 internal_units = "kg m^-2 second^-1",
626 external_units = "kg m^-2 year^-1";
627 if (kind == MASS) {
628 name = "debm_insolation_driven_melt_rate";
629 accumulator_units = "kg";
630 internal_units = "kg second^-1";
631 external_units = "Gt year^-1";
632 }
633
634 m_accumulator.metadata().units(accumulator_units);
635
636 m_vars = { { m_sys, name } };
637 m_vars[0]
638 .long_name(long_name)
639 .units(internal_units)
640 .output_units(external_units);
641 m_vars[0].set_string("cell_methods", "time: mean");
642
643 double fill_value = units::convert(m_sys, m_fill_value, external_units, internal_units);
644 m_vars[0].set_number("_FillValue", fill_value);
645 }
646
647protected:
649 const auto &melt_amount = model->insolation_driven_melt();
650
651 if (m_kind == MASS) {
652 m_melt_mass.copy_from(melt_amount);
653 m_melt_mass.scale(m_grid->cell_area());
654 return m_melt_mass;
655 }
656
657 return melt_amount;
658 }
659
660private:
663};
664
665/*! @brief Report surface temperature melt, averaged over the reporting interval */
666class DEBMSTemperatureMelt : public DiagAverageRate<DEBMSimple> {
667public:
670 kind == AMOUNT
671 ? "debm_temperature_driven_melt_flux"
672 : "debm_temperature_driven_melt_rate",
674 m_kind(kind),
675 m_melt_mass(m_grid, "temperature_melt_mass") {
676
677 std::string
678 name = "debm_temperature_driven_melt_flux",
679 long_name = "temperature-driven melt, averaged over the reporting interval",
680 accumulator_units = "kg m^-2",
681 internal_units = "kg m^-2 second^-1",
682 external_units = "kg m^-2 year^-1";
683 if (kind == MASS) {
684 name = "debm_temperature_driven_melt_rate";
685 accumulator_units = "kg";
686 internal_units = "kg second^-1";
687 external_units = "Gt year^-1";
688 }
689
690 m_accumulator.metadata().units(accumulator_units);
691
692 m_vars = { { m_sys, name } };
693 m_vars[0]
694 .long_name(long_name)
695 .units(internal_units)
696 .output_units(external_units);
697 m_vars[0].set_string("cell_methods", "time: mean");
698 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
699 }
700
701protected:
703 const auto &melt_amount = model->temperature_driven_melt();
704
705 if (m_kind == MASS) {
706 m_melt_mass.copy_from(melt_amount);
707 m_melt_mass.scale(m_grid->cell_area());
708 return m_melt_mass;
709 }
710
711 return melt_amount;
712 }
713
714private:
717};
718
719/*! @brief Report surface backround melt, averaged over the reporting interval */
720class DEBMSBackroundMelt : public DiagAverageRate<DEBMSimple> {
721public:
724 kind == AMOUNT
725 ? "debm_offset_melt_flux"
726 : "debm_offset_melt_rate",
728 m_kind(kind),
729 m_melt_mass(m_grid, "backround_melt_mass") {
730
731 std::string name = "debm_offset_melt_flux",
732 long_name = "offset melt, averaged over the reporting interval",
733 accumulator_units = "kg m^-2",
734 internal_units = "kg m^-2 second^-1",
735 external_units = "kg m^-2 year^-1";
736
737 if (kind == MASS) {
738 name = "debm_offset_melt_rate";
739 accumulator_units = "kg";
740 internal_units = "kg second^-1";
741 external_units = "Gt year^-1";
742 }
743 m_accumulator.metadata().units(accumulator_units);
744
745 m_vars = { { m_sys, name } };
746 m_vars[0]
747 .long_name(long_name)
748 .units(internal_units)
749 .output_units(external_units);
750 m_vars[0].set_string("cell_methods", "time: mean");
751 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
752 }
753
754protected:
756 const auto &melt_amount = model->offset_melt();
757
758 if (m_kind == MASS) {
759 m_melt_mass.copy_from(melt_amount);
760 m_melt_mass.scale(m_grid->cell_area());
761 return m_melt_mass;
762 }
763
764 return melt_amount;
765 }
766
767private:
770};
771
772} // end of namespace diagnostics
773
774/*! @brief The number of points for temperature and precipitation time-series.
775 */
776unsigned int DEBMSimple::timeseries_length(double dt) const {
777 double dt_years = dt / m_year_length;
778
779 return std::max(1U, static_cast<unsigned int>(ceil(m_n_per_year * dt_years)));
780}
781
783 using namespace diagnostics;
784
785 DiagnosticList result = {
786 { "debm_insolation_driven_melt_flux", Diagnostic::Ptr(new DEBMSInsolationMelt(this, AMOUNT)) },
787 { "debm_insolation_driven_melt_rate", Diagnostic::Ptr(new DEBMSInsolationMelt(this, MASS)) },
788 { "debm_temperature_driven_melt_flux", Diagnostic::Ptr(new DEBMSTemperatureMelt(this, AMOUNT)) },
789 { "debm_temperature_driven_melt_rate", Diagnostic::Ptr(new DEBMSTemperatureMelt(this, MASS)) },
790 { "debm_offset_melt_flux", Diagnostic::Ptr(new DEBMSBackroundMelt(this, AMOUNT)) },
791 { "debm_offset_melt_rate", Diagnostic::Ptr(new DEBMSBackroundMelt(this, MASS)) },
792 { "air_temp_sd", Diagnostic::wrap(*m_air_temp_sd) },
793 { "snow_depth", Diagnostic::wrap(m_snow_depth) },
794 { "surface_albedo", Diagnostic::wrap(m_surface_albedo) },
795 { "atmosphere_transmissivity", Diagnostic::wrap(m_transmissivity) },
796 { "insolation", Diagnostic::Ptr(new DEBMSInsolation(this)) }
797 };
798
800
801 return result;
802}
803
804} // end of namespace surface
805} // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition Component.hh:160
const Time & time() const
Definition Component.cc:109
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition Component.hh:156
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition Component.cc:159
DiagnosticList diagnostics() const
Definition Component.cc:89
const DEBMSimple * model
A template derived from Diagnostic, adding a "Model".
static Ptr wrap(const T &input)
double m_fill_value
fill value (used often enough to justify storing it)
const units::System::Ptr m_sys
the unit system
double to_internal(double x) const
Definition Diagnostic.cc:59
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
std::shared_ptr< Diagnostic > Ptr
Definition Diagnostic.hh:65
std::shared_ptr< const Grid > m_grid
the grid
High-level PISM I/O class.
Definition File.hh:55
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar latitude
Definition Geometry.hh:41
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
double increment_date(double T, double years) const
Definition Time.cc:841
double calendar_year_start(double T) const
Returns the model time in seconds corresponding to the beginning of the year T falls into.
Definition Time.cc:832
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & standard_name(const std::string &input)
VariableMetadata & set_number(const std::string &name, double value)
Set a scalar attribute to a single (scalar) value.
VariableMetadata & output_units(const std::string &input)
VariableMetadata & set_string(const std::string &name, const std::string &value)
Set a string attribute.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void copy_from(const Array2D< T > &source)
Definition Array2D.hh:73
void add(double alpha, const Array2D< T > &x)
Definition Array2D.hh:65
void read(const std::string &filename, unsigned int time)
Definition Array.cc:731
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition Array.cc:463
void write(const std::string &filename) const
Definition Array.cc:722
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition Array.cc:224
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
static std::shared_ptr< Forcing > Constant(std::shared_ptr< const Grid > grid, const std::string &short_name, double value)
Definition Forcing.cc:148
double albedo(double melt_rate, MaskValue cell_type) const
DEBMSimpleMelt melt(double declination, double distance_factor, double dt, double T_std_deviation, double T, double surface_elevation, double lat, double albedo) const
double atmosphere_transmissivity(double elevation) const
DEBMSimpleOrbitalParameters orbital_parameters(double time) const
DEBMSimpleChanges step(double ice_thickness, double max_melt, double snow_depth, double accumulation) const
Compute the surface mass balance at a location from the amount of melted snow and the solid accumulat...
A dEBM-simple implementation.
array::Scalar m_insolation_driven_melt
total insolation melt during the last time step
const array::Scalar & atmosphere_transmissivity() const
virtual const array::Scalar & mass_flux_impl() const
std::shared_ptr< array::Scalar > m_melt
total melt during the last time step
virtual void update_impl(const Geometry &geometry, double t, double dt)
double snow_accumulation(double T, double P) const
Extracts snow accumulation from mixed (snow and rain) precipitation using a temperature threshold wit...
const array::Scalar & insolation_driven_melt() const
const array::Scalar & surface_albedo() const
array::Scalar m_surface_albedo
albedo field
array::Scalar m_mass_flux
cached surface mass balance rate
Definition DEBMSimple.hh:87
std::shared_ptr< array::Scalar > m_runoff
total runoff during the last time step
array::Scalar m_temperature_driven_melt
total temperature melt during the last time step
const array::Scalar & air_temp_sd() const
array::Scalar m_snow_depth
snow depth (reset once a year)
Definition DEBMSimple.hh:92
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
std::shared_ptr< array::Forcing > m_air_temp_sd
standard deviation of the daily variability of the air temperature
Definition DEBMSimple.hh:95
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
double m_Tmax
the temperature above which all precipitation is rain
const array::Scalar & runoff_impl() const
unsigned int m_n_per_year
number of small time steps per year
array::Scalar m_transmissivity
transmissivity field
double compute_next_balance_year_start(double time)
std::shared_ptr< array::Scalar > m_accumulation
total accumulation during the last time step
Definition DEBMSimple.hh:98
unsigned int timeseries_length(double dt) const
The number of points for temperature and precipitation time-series.
DEBMSimplePointwise m_model
Definition DEBMSimple.hh:82
virtual DiagnosticList diagnostics_impl() const
const array::Scalar & offset_melt() const
const array::Scalar & accumulation_impl() const
const array::Scalar & snow_depth() const
const array::Scalar & temperature_driven_melt() const
virtual MaxTimestep max_timestep_impl(double t) const
const array::Scalar & melt_impl() const
array::Scalar m_offset_melt
total offset_melt during the last timestep
DEBMSimple(std::shared_ptr< const Grid > g, std::shared_ptr< atmosphere::AtmosphereModel > input)
Definition DEBMSimple.cc:42
virtual void init_impl(const Geometry &geometry)
double m_Tmin
the temperature below which all precipitation is snow
std::shared_ptr< array::Scalar > m_temperature
Definition DEBMSimple.hh:89
const DEBMSimplePointwise & pointwise_model() const
std::shared_ptr< array::Forcing > m_input_albedo
if albedo is given as input field
virtual const array::Scalar & temperature_impl() const
bool m_precip_as_snow
interpret all the precipitation as snow (no rain)
A class implementing a temperature-index (positive degree-day) scheme to compute melt and runoff,...
Definition DEBMSimple.hh:39
static std::shared_ptr< array::Scalar > allocate_runoff(std::shared_ptr< const Grid > grid)
std::shared_ptr< atmosphere::AtmosphereModel > m_atmosphere
virtual DiagnosticList diagnostics_impl() const
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
static std::shared_ptr< array::Scalar > allocate_temperature(std::shared_ptr< const Grid > grid)
const array::Scalar & accumulation() const
Returns accumulation.
static std::shared_ptr< array::Scalar > allocate_accumulation(std::shared_ptr< const Grid > grid)
static std::shared_ptr< array::Scalar > allocate_melt(std::shared_ptr< const Grid > grid)
virtual void init_impl(const Geometry &geometry)
The interface of PISM's surface models.
DEBMSBackroundMelt(const DEBMSimple *m, AmountKind kind)
Report surface backround melt, averaged over the reporting interval.
DEBMSInsolationMelt(const DEBMSimple *m, AmountKind kind)
Report surface insolation melt, averaged over the reporting interval.
std::shared_ptr< array::Array > compute_impl() const
Report mean top of atmosphere insolation.
DEBMSTemperatureMelt(const DEBMSimple *m, AmountKind kind)
Report surface temperature melt, averaged over the reporting interval.
#define PISM_ERROR_LOCATION
@ PISM_GUESS
Definition IO_Flags.hh:56
@ PISM_READONLY
open an existing file for reading only
Definition IO_Flags.hh:68
@ PISM_DOUBLE
Definition IO_Flags.hh:52
bool ice_free_ocean(int M)
Definition Mask.hh:61
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition Units.cc:70
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition Component.cc:43
static const double g
Definition exactTestP.cc:36
@ INIT_BOOTSTRAP
Definition Component.hh:56
@ INIT_RESTART
Definition Component.hh:56
std::map< std::string, Diagnostic::Ptr > DiagnosticList
static const double k
Definition exactTestP.cc:42
MaskValue
Definition Mask.hh:30
T combine(const T &a, const T &b)
InitializationType type
initialization type
Definition Component.hh:61
std::string filename
name of the input file (if applicable)
Definition Component.hh:63
unsigned int record
index of the record to re-start from
Definition Component.hh:65
static double S(unsigned n)
Definition test_cube.c:58