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