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
Hydrology.cc
Go to the documentation of this file.
1// Copyright (C) 2012-2024 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 "pism/hydrology/Hydrology.hh"
20#include "pism/util/error_handling.hh"
21#include "pism/util/io/File.hh"
22#include "pism/util/array/CellType.hh"
23#include "pism/geometry/Geometry.hh"
24
25namespace pism {
26namespace hydrology {
27
28namespace diagnostics {
29
30class TendencyOfWaterMass : public DiagAverageRate<Hydrology>
31{
32public:
34 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass", TOTAL_CHANGE) {
35
36 m_accumulator.metadata()["units"] = "kg";
37
38 m_vars = { { m_sys, "tendency_of_subglacial_water_mass" } };
39 m_vars[0]
40 .long_name("rate of change of the total mass of subglacial water")
41 .units("kg second^-1")
42 .output_units("Gt year^-1");
43 m_vars[0]["cell_methods"] = "time: mean";
44
45 m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
46 m_vars[0]["comment"] = "positive flux corresponds to water gain";
47 }
48
49protected:
51 return model->mass_change();
52 }
53};
54
55/*! @brief Report water input rate from the ice surface into the subglacial water system.
56 */
57class TotalInputRate : public DiagAverageRate<Hydrology>
58{
59public:
61 : DiagAverageRate<Hydrology>(m, "subglacial_water_input_rate_from_surface", RATE) {
62
63 m_accumulator.metadata()["units"] = "m";
64
65 m_vars = { { m_sys, "subglacial_water_input_rate_from_surface" } };
66 m_vars[0]
67 .long_name("water input rate from the ice surface into the subglacial water system")
68 .units("m second^-1")
69 .output_units("m year^-1");
70 m_vars[0]["cell_methods"] = "time: mean";
71 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
72 m_vars[0]["comment"] = "positive flux corresponds to water gain";
73 }
74
75protected:
77 return model->surface_input_rate();
78 }
79};
80
81/*! @brief Report water flux due to input (basal melt and drainage from the surface). */
82class WaterInputFlux : public DiagAverageRate<Hydrology> {
83public:
85 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_due_to_input",
87 m_accumulator.metadata()["units"] = "kg";
88
89 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_due_to_input" } };
90 m_vars[0]
91 .long_name("subglacial water flux due to input")
92 .units("kg second^-1")
93 .output_units("Gt year^-1");
94 m_vars[0]["cell_methods"] = "time: mean";
95 m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
96 m_vars[0]["comment"] = "positive flux corresponds to water gain";
97 }
98
99protected:
101 return model->mass_change_due_to_input();
102 }
103};
104
105/*! @brief Report advective subglacial water flux. */
106class SubglacialWaterFlux : public DiagAverageRate<Hydrology>
107{
108public:
110 : DiagAverageRate<Hydrology>(m, "subglacial_water_flux", RATE),
111 m_flux_magnitude(m_grid, "flux_magnitude"){
112
113 m_accumulator.metadata()["units"] = "m^2";
114
115 m_vars = { { m_sys, "subglacial_water_flux" } };
116 m_vars[0]
117 .long_name("magnitude of the subglacial water flux")
118 .units("m^2 second^-1")
119 .output_units("m^2 year^-1");
120 m_vars[0]["cell_methods"] = "time: mean";
121 m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
122
124 .long_name("magnitude of the subglacial water flux")
125 .units("m^2 s^-1");
126 }
127
128protected:
129 void update_impl(double dt) {
130 compute_magnitude(model->flux(), m_flux_magnitude);
131
133
134 m_interval_length += dt;
135 }
136
138};
139
140
141/*! @brief Report water flux at the grounded margin. */
142class GroundedMarginFlux : public DiagAverageRate<Hydrology> {
143public:
145 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_at_grounded_margins",
146 TOTAL_CHANGE) {
147
148 m_accumulator.metadata()["units"] = "kg";
149
150 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_at_grounded_margins" } };
151 m_vars[0]
152 .long_name("subglacial water flux at grounded ice margins")
153 .units("kg second^-1")
154 .output_units("Gt year^-1");
155 m_vars[0]["cell_methods"] = "time: mean";
156 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
157 m_vars[0]["comment"] = "positive flux corresponds to water gain";
158 }
159
160protected:
162 return model->mass_change_at_grounded_margin();
163 }
164};
165
166/*! @brief Report subglacial water flux at grounding lines. */
167class GroundingLineFlux : public DiagAverageRate<Hydrology> {
168public:
170 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_at_grounding_line",
171 TOTAL_CHANGE) {
172
173 m_accumulator.metadata()["units"] = "kg";
174
175 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_at_grounding_line" } };
176 m_vars[0]
177 .long_name("subglacial water flux at grounding lines")
178 .units("kg second^-1")
179 .output_units("Gt year^-1");
180 m_vars[0]["cell_methods"] = "time: mean";
181
182 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
183 m_vars[0]["comment"] = "positive flux corresponds to water gain";
184 }
185
186protected:
188 return model->mass_change_at_grounding_line();
189 }
190};
191
192/*! @brief Report subglacial water conservation error flux (mass added to preserve non-negativity). */
193class ConservationErrorFlux : public DiagAverageRate<Hydrology> {
194public:
196 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_due_to_conservation_error",
197 TOTAL_CHANGE) {
198
199 m_accumulator.metadata()["units"] = "kg";
200
201 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_due_to_conservation_error" } };
202 m_vars[0]
203 .long_name(
204 "subglacial water flux due to conservation error (mass added to preserve non-negativity)")
205 .units("kg second^-1")
206 .output_units("Gt year^-1");
207 m_vars[0]["cell_methods"] = "time: mean";
208 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
209 m_vars[0]["comment"] = "positive flux corresponds to water gain";
210 }
211
212protected:
214 return model->mass_change_due_to_conservation_error();
215 }
216};
217
218/*! @brief Report subglacial water flux at domain boundary (in regional model configurations). */
219class DomainBoundaryFlux : public DiagAverageRate<Hydrology> {
220public:
222 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_at_domain_boundary",
223 TOTAL_CHANGE) {
224
225 m_accumulator.metadata()["units"] = "kg";
226
227 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_at_domain_boundary" } };
228 m_vars[0]
229 .long_name("subglacial water flux at domain boundary (in regional model configurations)")
230 .units("kg second^-1")
231 .output_units("Gt year^-1");
232 m_vars[0]["cell_methods"] = "time: mean";
233 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
234 m_vars[0]["comment"] = "positive flux corresponds to water gain";
235 }
236
237protected:
239 return model->mass_change_at_domain_boundary();
240 }
241};
242
243/*! @brief Report water flux at the grounded margin. */
245public:
247 : DiagAverageRate<Hydrology>(m, "tendency_of_subglacial_water_mass_due_to_flow",
248 TOTAL_CHANGE) {
249
250 m_accumulator.metadata()["units"] = "kg";
251
252 m_vars = { { m_sys, "tendency_of_subglacial_water_mass_due_to_flow" } };
253 m_vars[0]
254 .long_name("rate of change subglacial water mass due to lateral flow")
255 .units("kg second^-1")
256 .output_units("Gt year^-1");
257 m_vars[0]["cell_methods"] = "time: mean";
258 m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
259 m_vars[0]["comment"] = "positive flux corresponds to water gain";
260 }
261
262protected:
264 return model->mass_change_due_to_lateral_flow();
265 }
266};
267
268} // end of namespace diagnostics
269
271 geometry = nullptr;
272 surface_input_rate = nullptr;
273 basal_melt_rate = nullptr;
274 ice_sliding_speed = nullptr;
275 no_model_mask = nullptr;
276}
277
278Hydrology::Hydrology(std::shared_ptr<const Grid> g)
279 : Component(g),
280 m_Q(m_grid, "water_flux"),
281 m_Wtill(m_grid, "tillwat"),
282 m_W(m_grid, "bwat"),
283 m_Pover(m_grid, "overburden_pressure"),
284 m_surface_input_rate(m_grid, "water_input_rate_from_surface"),
285 m_basal_melt_rate(m_grid, "water_input_rate_due_to_basal_melt"),
286 m_flow_change_incremental(m_grid, "water_thickness_change_due_to_flow"),
287 m_conservation_error_change(m_grid, "conservation_error_change"),
288 m_grounded_margin_change(m_grid, "grounded_margin_change"),
289 m_grounding_line_change(m_grid, "grounding_line_change"),
290 m_input_change(m_grid, "water_mass_change_due_to_input"),
291 m_no_model_mask_change(m_grid, "no_model_mask_change"),
292 m_total_change(m_grid, "water_mass_change"),
293 m_flow_change(m_grid, "water_mass_change_due_to_flow") {
294
296 .long_name("hydrology model workspace for water input rate from the ice surface")
297 .units("m s^-1");
298
300 .long_name("hydrology model workspace for water input rate due to basal melt")
301 .units("m s^-1");
302
303 // *all* Hydrology classes have layer of water stored in till as a state variable
305 .long_name("effective thickness of subglacial water stored in till")
306 .units("m");
307 m_Wtill.metadata()["valid_min"] = { 0.0 };
308
309 m_Pover.metadata(0).long_name("overburden pressure").units("Pa");
310 m_Pover.metadata()["valid_min"] = { 0.0 };
311
312 // needs ghosts in Routing and Distributed
313 m_W.metadata(0)
314 .long_name("thickness of transportable subglacial water layer")
315 .units("m");
316 m_W.metadata()["valid_min"] = { 0.0 };
317
318 m_Q.metadata(0)
319 .long_name("advective subglacial water flux")
320 .units("m^2 s^-1")
321 .output_units("m^2 day^-1");
322 m_Q.set(0.0);
323
324 // storage for water conservation reporting quantities
325 m_total_change.metadata(0).long_name("total change in water mass over one time step").units("kg");
326
328 .long_name(
329 "change in water mass over one time step due to the input (basal melt and surface drainage)")
330 .units("kg");
331
333 .long_name("change in water mass due to lateral flow (over one time step)")
334 .units("kg");
335
337 .long_name("changes in subglacial water thickness at the grounded margin")
338 .units("kg");
339
341 .long_name("changes in subglacial water thickness at the grounding line")
342 .units("kg");
343
345 .long_name(
346 "changes in subglacial water thickness at the edge of the modeling domain (regional models)")
347 .units("kg");
348
350 .long_name(
351 "changes in subglacial water thickness required to preserve non-negativity or keep water thickness within bounds")
352 .units("kg");
353}
354
355void Hydrology::restart(const File &input_file, int record) {
357 this->restart_impl(input_file, record);
358}
359
360void Hydrology::bootstrap(const File &input_file, const array::Scalar &ice_thickness) {
362 this->bootstrap_impl(input_file, ice_thickness);
363}
364
365void Hydrology::init(const array::Scalar &W_till, const array::Scalar &W, const array::Scalar &P) {
367 this->init_impl(W_till, W, P);
368}
369
370void Hydrology::restart_impl(const File &input_file, int record) {
371 m_Wtill.read(input_file, record);
372
373 // whether or not we could initialize from file, we could be asked to regrid from file
374 regrid("Hydrology", m_Wtill);
375}
376
377void Hydrology::bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness) {
378 (void)ice_thickness;
379
380 double tillwat_default = m_config->get_number("bootstrapping.defaults.tillwat");
381 m_Wtill.regrid(input_file, io::Default(tillwat_default));
382
383 // whether or not we could initialize from file, we could be asked to regrid from file
384 regrid("Hydrology", m_Wtill);
385}
386
388 const array::Scalar &P) {
389 (void)W;
390 (void)P;
391 m_Wtill.copy_from(W_till);
392}
393
394void Hydrology::update(double t, double dt, const Inputs &inputs) {
395
396 // reset water thickness changes
397 {
402 m_flow_change.set(0.0);
403 m_input_change.set(0.0);
404 }
405
407
411
413
414 for (auto p = m_grid->points(); p; p.next()) {
415 const int i = p.i(), j = p.j();
416
417 m_total_change(i, j) = m_W(i, j) + m_Wtill(i, j);
418 }
419
420 this->update_impl(t, dt, inputs);
421
422 // compute total change in meters
423 for (auto p = m_grid->points(); p; p.next()) {
424 const int i = p.i(), j = p.j();
425 m_total_change(i, j) = (m_W(i, j) + m_Wtill(i, j)) - m_total_change(i, j);
426 }
427
428 // convert from m to kg
429 // kg = m * (kg / m^3) * m^2
430
431 double water_density = m_config->get_number("constants.fresh_water.density"),
432 kg_per_m = water_density * m_grid->cell_area();
433
434 list.add({ &m_flow_change, &m_input_change });
435 for (auto p = m_grid->points(); p; p.next()) {
436 const int i = p.i(), j = p.j();
437 m_total_change(i, j) *= kg_per_m;
438 m_input_change(i, j) *= kg_per_m;
439 m_flow_change(i, j) *= kg_per_m;
440 }
441}
442
444 using namespace diagnostics;
445 DiagnosticList result = {
446 { "bwat", Diagnostic::wrap(m_W) },
447 { "tillwat", Diagnostic::wrap(m_Wtill) },
448 { "subglacial_water_input_rate", Diagnostic::Ptr(new TotalInputRate(this)) },
449 { "tendency_of_subglacial_water_mass_due_to_input", Diagnostic::Ptr(new WaterInputFlux(this)) },
450 { "tendency_of_subglacial_water_mass_due_to_flow",
451 Diagnostic::Ptr(new TendencyOfWaterMassDueToFlow(this)) },
452 { "tendency_of_subglacial_water_mass_due_to_conservation_error",
453 Diagnostic::Ptr(new ConservationErrorFlux(this)) },
454 { "tendency_of_subglacial_water_mass", Diagnostic::Ptr(new TendencyOfWaterMass(this)) },
455 { "tendency_of_subglacial_water_mass_at_grounded_margins",
456 Diagnostic::Ptr(new GroundedMarginFlux(this)) },
457 { "tendency_of_subglacial_water_mass_at_grounding_line",
458 Diagnostic::Ptr(new GroundingLineFlux(this)) },
459 { "tendency_of_subglacial_water_mass_at_domain_boundary",
460 Diagnostic::Ptr(new DomainBoundaryFlux(this)) },
461 { "subglacial_water_flux_mag", Diagnostic::Ptr(new SubglacialWaterFlux(this)) },
462 };
463
464 return result;
465}
466
467void Hydrology::define_model_state_impl(const File &output) const {
469}
470
471void Hydrology::write_model_state_impl(const File &output) const {
472 m_Wtill.write(output);
473}
474
475//! Update the overburden pressure from ice thickness.
476/*!
477 Uses the standard hydrostatic (shallow) approximation of overburden pressure,
478 \f[ P_0 = \rho_i g H \f]
479*/
481 array::Scalar &result) const {
482 // FIXME issue #15
483
484 const double ice_density = m_config->get_number("constants.ice.density"),
485 standard_gravity = m_config->get_number("constants.standard_gravity");
486
487 array::AccessScope list{ &ice_thickness, &result };
488
489 for (auto p = m_grid->points(); p; p.next()) {
490 const int i = p.i(), j = p.j();
491
492 result(i, j) = ice_density * standard_gravity * ice_thickness(i, j);
493 }
494}
495
497 return m_Pover;
498}
499
500//! Return the effective thickness of the water stored in till.
502 return m_Wtill;
503}
504
505//! Return the effective thickness of the transportable basal water layer.
509
510/*!
511 * Return subglacial water flux (time-average over the time step requested at the time of
512 * the update() call).
513 */
515 return m_Q;
516}
517
521
525
529
533
537
539 return m_total_change;
540}
541
545
549
550/*!
551 Checks @f$ 0 \le W \le W^{max} @f$.
552*/
553void check_bounds(const array::Scalar &W, double W_max) {
554
555 std::string name = W.metadata()["long_name"];
556
557 auto grid = W.grid();
558
559 array::AccessScope list(W);
560 ParallelSection loop(grid->com);
561 try {
562 for (auto p = grid->points(); p; p.next()) {
563 const int i = p.i(), j = p.j();
564
565 if (W(i,j) < 0.0) {
567 "Hydrology: negative %s of %.6f m at (i, j) = (%d, %d)",
568 name.c_str(), W(i, j), i, j);
569 }
570
571 if (W(i,j) > W_max) {
573 "Hydrology: %s of %.6f m exceeds the maximum value %.6f at (i, j) = (%d, %d)",
574 name.c_str(), W(i, j), W_max, i, j);
575 }
576 }
577 } catch (...) {
578 loop.failed();
579 }
580 loop.check();
581}
582
583
584//! Compute the surface water input rate into the basal hydrology layer in the ice-covered
585//! region.
586/*!
587 This method ignores the input rate in the ice-free region.
588
589 @param[in] mask cell type mask
590 @param[in] surface_input_rate surface input rate (kg m-2 s-1); set to NULL to ignore
591 @param[out] result resulting input rate (water thickness per time)
592*/
594 const array::Scalar *surface_input_rate,
595 array::Scalar &result) {
596
597 if (not surface_input_rate) {
598 result.set(0.0);
599 return;
600 }
601
602 array::AccessScope list{surface_input_rate, &mask, &result};
603
604 const double
605 water_density = m_config->get_number("constants.fresh_water.density");
606
607 for (auto p = m_grid->points(); p; p.next()) {
608 const int i = p.i(), j = p.j();
609
610 if (mask.icy(i, j)) {
611 result(i,j) = (*surface_input_rate)(i, j) / water_density;
612 } else {
613 result(i,j) = 0.0;
614 }
615 }
616}
617
618//! Compute the input rate into the basal hydrology layer in the ice-covered
619//! region due to basal melt rate.
620/*!
621 This method ignores the input in the ice-free region.
622
623 @param[in] mask cell type mask
624 @param[in] basal_melt_rate basal melt rate (ice thickness per time)
625 @param[out] result resulting input rate (water thickness per time)
626*/
628 const array::Scalar &basal_melt_rate,
629 array::Scalar &result) {
630
631 array::AccessScope list{&basal_melt_rate, &mask, &result};
632
633 const double
634 ice_density = m_config->get_number("constants.ice.density"),
635 water_density = m_config->get_number("constants.fresh_water.density"),
636 C = ice_density / water_density;
637
638 for (auto p = m_grid->points(); p; p.next()) {
639 const int i = p.i(), j = p.j();
640
641 if (mask.icy(i, j)) {
642 result(i,j) = C * basal_melt_rate(i, j);
643 } else {
644 result(i,j) = 0.0;
645 }
646 }
647}
648
649//! Correct the new water thickness based on boundary requirements.
650/*! At ice free locations and ocean locations we require that water thicknesses (i.e. both
651 the transportable water thickness \f$W\f$ and the till water thickness \f$W_{till}\f$)
652 be zero at the end of each time step. Also we require that any negative water
653 thicknesses be set to zero (i.e. we do projection to enforce lower bound).
654
655 This method should be called once for each thickness field which needs to be processed.
656 This method alters the field water_thickness in-place.
657
658 @param[in] cell_type cell type mask
659 @param[in] no_model_mask (optional) mask of zeros and ones, zero within the modeling
660 domain, one outside
661 @param[in] max_thickness maximum allowed water thickness (use a zero or a negative value
662 to disable)
663 @param[in,out] water_thickness adjusted water thickness (till storage or the transport system)
664 @param[in,out] grounded_margin_change change in water thickness at the grounded margin
665 @param[in,out] grounding_line_change change in water thickness at the grounding line
666 @param[in,out] conservation_error_change change in water thickness due to mass conservation errors
667 @param[in,out] no_model_mask_change change in water thickness outside the modeling domain (regional models)
668*/
670 const array::Scalar *no_model_mask,
671 double max_thickness,
672 double ocean_water_thickness,
673 array::Scalar &water_thickness,
674 array::Scalar &grounded_margin_change,
675 array::Scalar &grounding_line_change,
676 array::Scalar &conservation_error_change,
677 array::Scalar &no_model_mask_change) {
678
679 bool include_floating = m_config->get_flag("hydrology.routing.include_floating_ice");
680
681 array::AccessScope list{&water_thickness, &cell_type,
682 &grounded_margin_change, &grounding_line_change, &conservation_error_change,
683 &no_model_mask_change};
684
685 double
686 fresh_water_density = m_config->get_number("constants.fresh_water.density"),
687 kg_per_m = m_grid->cell_area() * fresh_water_density; // kg m-1
688
689 for (auto p = m_grid->points(); p; p.next()) {
690 const int i = p.i(), j = p.j();
691
692 if (water_thickness(i, j) < 0.0) {
693 conservation_error_change(i, j) += -water_thickness(i, j) * kg_per_m;
694 water_thickness(i, j) = 0.0;
695 }
696
697 if (max_thickness > 0.0 and water_thickness(i, j) > max_thickness) {
698 double excess = water_thickness(i, j) - max_thickness;
699
700 conservation_error_change(i, j) += -excess * kg_per_m;
701 water_thickness(i, j) = max_thickness;
702 }
703
704 if (cell_type.ice_free_land(i, j)) {
705 double grounded_ice_free_max_thickness = 0.0;
706 double excess = water_thickness(i, j) - grounded_ice_free_max_thickness;
707 grounded_margin_change(i, j) += -excess * kg_per_m;
708 water_thickness(i, j) = grounded_ice_free_max_thickness;
709 }
710
711 // This keeps track of water mass changes at the grounding line due to
712 //
713 // - water leaving the system (drainage into the ocean) and
714 //
715 // - water added to the system when the sea level rises and previously grounded areas
716 // come in contact with the ocean.
717 //
718 // If the sea level rises and covers previously ice free land, the till water amount
719 // at that location should change to the maximum till water thickness.
720 //
721 // When the sea level recedes, the till water at that location will be set to zero by
722 // the if block above. All these changes will be accounted for.
723 if ((include_floating and cell_type.ice_free_ocean(i, j)) or
724 (not include_floating and cell_type.ocean(i, j))) {
725
726 double mismatch = water_thickness(i, j) - ocean_water_thickness;
727
728 grounding_line_change(i, j) += -mismatch * kg_per_m;
729 water_thickness(i, j) = ocean_water_thickness;
730 }
731 }
732
733 if (no_model_mask != nullptr) {
734 const array::Scalar &M = *no_model_mask;
735
736 list.add(M);
737
738 for (auto p = m_grid->points(); p; p.next()) {
739 const int i = p.i(), j = p.j();
740
741 if (M(i, j) > 0.0) {
742 no_model_mask_change(i, j) += -water_thickness(i, j) * kg_per_m;
743
744 water_thickness(i, j) = 0.0;
745 }
746 }
747 }
748}
749
750} // end of namespace hydrology
751} // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
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
A class defining a common interface for most PISM sub-models.
Definition Component.hh:118
const Model * 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::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
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
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & output_units(const std::string &input)
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
std::shared_ptr< const Grid > grid() const
Definition Array.cc:131
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
bool ice_free_ocean(int i, int j) const
Definition CellType.hh:58
bool ocean(int i, int j) const
Definition CellType.hh:34
bool ice_free_land(int i, int j) const
Definition CellType.hh:62
bool icy(int i, int j) const
Definition CellType.hh:42
"Cell type" mask. Adds convenience methods to array::Scalar.
Definition CellType.hh:30
array::Scalar m_flow_change
Definition Hydrology.hh:196
const array::Scalar & till_water_thickness() const
Return the effective thickness of the water stored in till.
Definition Hydrology.cc:501
void restart(const File &input_file, int record)
Definition Hydrology.cc:355
Hydrology(std::shared_ptr< const Grid > g)
Definition Hydrology.cc:278
virtual void bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness)
Definition Hydrology.cc:377
array::Scalar m_surface_input_rate
Definition Hydrology.hh:179
array::Scalar m_total_change
Definition Hydrology.hh:195
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition Hydrology.cc:467
void bootstrap(const File &input_file, const array::Scalar &ice_thickness)
Definition Hydrology.cc:360
array::Scalar m_grounding_line_change
Definition Hydrology.hh:192
const array::Scalar & mass_change_at_domain_boundary() const
Definition Hydrology.cc:534
const array::Scalar & surface_input_rate() const
Definition Hydrology.cc:518
virtual void initialization_message() const =0
void enforce_bounds(const array::CellType &cell_type, const array::Scalar *no_model_mask, double max_thickness, double ocean_water_thickness, array::Scalar &water_thickness, array::Scalar &grounded_margin_change, array::Scalar &grounding_line_change, array::Scalar &conservation_error_change, array::Scalar &no_model_mask_change)
Correct the new water thickness based on boundary requirements.
Definition Hydrology.cc:669
const array::Scalar & subglacial_water_thickness() const
Return the effective thickness of the transportable basal water layer.
Definition Hydrology.cc:506
array::Scalar m_basal_melt_rate
Definition Hydrology.hh:182
const array::Scalar & mass_change_due_to_lateral_flow() const
Definition Hydrology.cc:546
virtual void init_impl(const array::Scalar &W_till, const array::Scalar &W, const array::Scalar &P)
Definition Hydrology.cc:387
const array::Scalar & mass_change_due_to_conservation_error() const
Definition Hydrology.cc:530
array::Scalar1 m_W
effective thickness of transportable basal water
Definition Hydrology.hh:173
void compute_overburden_pressure(const array::Scalar &ice_thickness, array::Scalar &result) const
Update the overburden pressure from ice thickness.
Definition Hydrology.cc:480
const array::Scalar & mass_change_at_grounded_margin() const
Definition Hydrology.cc:522
virtual void restart_impl(const File &input_file, int record)
Definition Hydrology.cc:370
array::Scalar m_input_change
Definition Hydrology.hh:193
void update(double t, double dt, const Inputs &inputs)
Definition Hydrology.cc:394
array::Scalar m_Wtill
effective thickness of basal water stored in till
Definition Hydrology.hh:170
const array::Scalar & overburden_pressure() const
Definition Hydrology.cc:496
virtual void update_impl(double t, double dt, const Inputs &inputs)=0
const array::Scalar & mass_change_due_to_input() const
Definition Hydrology.cc:542
array::Scalar m_Pover
overburden pressure
Definition Hydrology.hh:176
virtual std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
Definition Hydrology.cc:443
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition Hydrology.cc:471
const array::Vector & flux() const
Definition Hydrology.cc:514
void compute_basal_melt_rate(const array::CellType &mask, const array::Scalar &basal_melt_rate, array::Scalar &result)
Definition Hydrology.cc:627
array::Scalar m_grounded_margin_change
Definition Hydrology.hh:191
const array::Scalar & mass_change_at_grounding_line() const
Definition Hydrology.cc:526
array::Scalar m_no_model_mask_change
Definition Hydrology.hh:194
void compute_surface_input_rate(const array::CellType &mask, const array::Scalar *surface_input_rate, array::Scalar &result)
Definition Hydrology.cc:593
void init(const array::Scalar &W_till, const array::Scalar &W, const array::Scalar &P)
Definition Hydrology.cc:365
const array::Scalar & mass_change() const
Definition Hydrology.cc:538
array::Scalar m_conservation_error_change
Definition Hydrology.hh:190
The PISM subglacial hydrology model interface.
Definition Hydrology.hh:109
const Geometry * geometry
Definition Hydrology.hh:37
const array::Scalar * ice_sliding_speed
Definition Hydrology.hh:41
const array::Scalar * surface_input_rate
Definition Hydrology.hh:39
const array::Scalar1 * no_model_mask
Definition Hydrology.hh:35
const array::Scalar * basal_melt_rate
Definition Hydrology.hh:40
Report subglacial water conservation error flux (mass added to preserve non-negativity).
Definition Hydrology.cc:193
Report subglacial water flux at domain boundary (in regional model configurations).
Definition Hydrology.cc:219
Report water flux at the grounded margin.
Definition Hydrology.cc:142
Report subglacial water flux at grounding lines.
Definition Hydrology.cc:167
Report advective subglacial water flux.
Definition Hydrology.cc:107
Report water flux at the grounded margin.
Definition Hydrology.cc:244
Report water input rate from the ice surface into the subglacial water system.
Definition Hydrology.cc:58
Report water flux due to input (basal melt and drainage from the surface).
Definition Hydrology.cc:82
#define PISM_ERROR_LOCATION
void check_bounds(const array::Scalar &W, double W_max)
Definition Hydrology.cc:553
@ PISM_DOUBLE
Definition IO_Flags.hh:52
static const double g
Definition exactTestP.cc:36
std::map< std::string, Diagnostic::Ptr > DiagnosticList