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