PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Routing.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 <cassert>
20 
21 #include "Routing.hh"
22 #include "pism/util/IceModelVec2CellType.hh"
23 #include "pism/util/Mask.hh"
24 #include "pism/util/MaxTimestep.hh"
25 
26 #include "pism/util/error_handling.hh"
27 
28 #include "pism/util/pism_options.hh"
29 #include "pism/util/pism_utilities.hh"
30 #include "pism/util/Vars.hh"
31 #include "pism/geometry/Geometry.hh"
32 #include "pism/util/Profiling.hh"
33 #include "pism/util/Context.hh"
34 
35 namespace pism {
36 namespace hydrology {
37 
38 namespace diagnostics {
39 
40 //! \brief Reports the pressure of the transportable water in the subglacial layer.
41 class BasalWaterPressure : public Diag<Routing>
42 {
43 public:
45  : Diag<Routing>(m) {
47  set_attrs("pressure of transportable water in subglacial layer", "", "Pa", "Pa", 0);
48  }
49 
50 protected:
51  virtual IceModelVec::Ptr compute_impl() const {
53  result->metadata() = m_vars[0];
54  result->copy_from(model->subglacial_water_pressure());
55  return result;
56  }
57 };
58 
59 
60 //! \brief Reports the pressure of the transportable water in the subglacial layer as a
61 //! fraction of the overburden pressure.
62 class RelativeBasalWaterPressure : public Diag<Routing>
63 {
64 public:
66  : Diag<Routing>(m) {
67  m_vars = {SpatialVariableMetadata(m_sys, "bwprel")};
68  set_attrs("pressure of transportable water in subglacial layer"
69  " as fraction of the overburden pressure", "",
70  "", "", 0);
71  m_vars[0]["_FillValue"] = {m_fill_value};
72  }
73 
74 protected:
75  virtual IceModelVec::Ptr compute_impl() const {
76  double fill_value = m_fill_value;
77 
78  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "bwprel", WITHOUT_GHOSTS));
79  result->metadata(0) = m_vars[0];
80 
81  const IceModelVec2S
82  &P = model->subglacial_water_pressure(),
83  &Po = model->overburden_pressure();
84 
85  IceModelVec::AccessList list{result.get(), &Po, &P};
86  for (Points p(*m_grid); p; p.next()) {
87  const int i = p.i(), j = p.j();
88 
89  if (Po(i,j) > 0.0) {
90  (*result)(i,j) = P(i, j) / Po(i,j);
91  } else {
92  (*result)(i,j) = fill_value;
93  }
94  }
95 
96  return result;
97  }
98 };
99 
100 
101 //! \brief Reports the effective pressure of the transportable water in the subglacial
102 //! layer, that is, the overburden pressure minus the pressure.
103 class EffectiveBasalWaterPressure : public Diag<Routing>
104 {
105 public:
107  : Diag<Routing>(m) {
108  m_vars = {SpatialVariableMetadata(m_sys, "effbwp")};
109  set_attrs("effective pressure of transportable water in subglacial layer"
110  " (overburden pressure minus water pressure)",
111  "", "Pa", "Pa", 0);
112  }
113 
114 protected:
115  virtual IceModelVec::Ptr compute_impl() const {
116 
117  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "effbwp", WITHOUT_GHOSTS));
118  result->metadata() = m_vars[0];
119 
120  const IceModelVec2S
121  &P = model->subglacial_water_pressure(),
122  &Po = model->overburden_pressure();
123 
124  IceModelVec::AccessList list{&Po, &P, result.get()};
125 
126  for (Points p(*m_grid); p; p.next()) {
127  const int i = p.i(), j = p.j();
128 
129  (*result)(i, j) = Po(i, j) - P(i, j);
130  }
131 
132  return result;
133  }
134 };
135 
136 
137 //! \brief Report the wall melt rate from dissipation of the potential energy of the
138 //! transportable water.
139 class WallMelt : public Diag<Routing>
140 {
141 public:
142  WallMelt(const Routing *m)
143  : Diag<Routing>(m) {
144  m_vars = {SpatialVariableMetadata(m_sys, "wallmelt")};
145  set_attrs("wall melt into subglacial hydrology layer"
146  " from (turbulent) dissipation of energy in transportable water",
147  "", "m s-1", "m year-1", 0);
148  }
149 
150 protected:
151  virtual IceModelVec::Ptr compute_impl() const {
152  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "wallmelt", WITHOUT_GHOSTS));
153  result->metadata() = m_vars[0];
154 
155  const IceModelVec2S &bed_elevation = *m_grid->variables().get_2d_scalar("bedrock_altitude");
156 
157  wall_melt(*model, bed_elevation, *result);
158  return result;
159  }
160 };
161 
162 //! @brief Diagnostically reports the staggered-grid components of the velocity of the
163 //! water in the subglacial layer.
164 class BasalWaterVelocity : public Diag<Routing>
165 {
166 public:
168  : Diag<Routing>(m) {
169 
170  // set metadata:
171  m_vars = {SpatialVariableMetadata(m_sys, "bwatvel[0]"),
172  SpatialVariableMetadata(m_sys, "bwatvel[1]")};
173 
174  set_attrs("velocity of water in subglacial layer, i-offset", "",
175  "m s-1", "m s-1", 0);
176  set_attrs("velocity of water in subglacial layer, j-offset", "",
177  "m s-1", "m s-1", 1);
178  }
179 protected:
180  virtual IceModelVec::Ptr compute_impl() const {
182  result->metadata(0) = m_vars[0];
183  result->metadata(1) = m_vars[1];
184 
185  result->copy_from(model->velocity_staggered());
186 
187  return result;
188  }
189 };
190 
191 //! Compute the hydraulic potential.
192 /*!
193  Computes \f$\psi = P + \rho_w g (b + W)\f$.
194 */
196  const IceModelVec2S &P,
197  const IceModelVec2S &sea_level,
198  const IceModelVec2S &bed,
199  const IceModelVec2S &ice_thickness,
200  IceModelVec2S &result) {
201 
202  IceGrid::ConstPtr grid = result.grid();
203 
204  Config::ConstPtr config = grid->ctx()->config();
205 
206  double
207  ice_density = config->get_number("constants.ice.density"),
208  sea_water_density = config->get_number("constants.sea_water.density"),
209  C = ice_density / sea_water_density,
210  rg = (config->get_number("constants.fresh_water.density") *
211  config->get_number("constants.standard_gravity"));
212 
213  IceModelVec::AccessList list{&P, &W, &sea_level, &ice_thickness, &bed, &result};
214 
215  for (Points p(*grid); p; p.next()) {
216  const int i = p.i(), j = p.j();
217 
218  double b = std::max(bed(i, j), sea_level(i, j) - C * ice_thickness(i, j));
219 
220  result(i, j) = P(i, j) + rg * (b + W(i, j));
221  }
222 }
223 
224 /*! @brief Report hydraulic potential in the subglacial hydrology system */
225 class HydraulicPotential : public Diag<Routing>
226 {
227 public:
229  : Diag<Routing>(m) {
230 
231  m_vars = {SpatialVariableMetadata(m_sys, "hydraulic_potential")};
232 
233  set_attrs("hydraulic potential in the subglacial hydrology system", "",
234  "Pa", "Pa", 0);
235  }
236 
237 protected:
239 
240  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "hydraulic_potential", WITHOUT_GHOSTS));
241  result->metadata(0) = m_vars[0];
242 
243  const IceModelVec2S &sea_level = *m_grid->variables().get_2d_scalar("sea_level");
244  const IceModelVec2S &bed_elevation = *m_grid->variables().get_2d_scalar("bedrock_altitude");
245  const IceModelVec2S &ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
246 
247  hydraulic_potential(model->subglacial_water_thickness(),
248  model->subglacial_water_pressure(),
249  sea_level,
250  bed_elevation,
251  ice_thickness,
252  *result);
253 
254  return result;
255  }
256 };
257 
258 } // end of namespace diagnostics
259 
261  : Hydrology(grid),
262  m_Qstag(grid, "advection_flux", WITH_GHOSTS, 1),
263  m_Qstag_average(grid, "cumulative_advection_flux", WITH_GHOSTS, 1),
264  m_Vstag(grid, "water_velocity", WITHOUT_GHOSTS),
265  m_Wstag(grid, "W_staggered", WITH_GHOSTS, 1),
266  m_Kstag(grid, "K_staggered", WITH_GHOSTS, 1),
267  m_Wnew(grid, "W_new", WITHOUT_GHOSTS),
268  m_Wtillnew(grid, "Wtill_new", WITHOUT_GHOSTS),
269  m_R(grid, "potential_workspace", WITH_GHOSTS, 1), /* box stencil used */
270  m_dx(grid->dx()),
271  m_dy(grid->dy()),
272  m_bottom_surface(grid, "ice_bottom_surface_elevation", WITH_GHOSTS) {
273 
274  m_W.metadata()["pism_intent"] = "model_state";
275 
276  m_rg = (m_config->get_number("constants.fresh_water.density") *
277  m_config->get_number("constants.standard_gravity"));
278 
279  m_Qstag.set_attrs("internal",
280  "cell face-centered (staggered) components of advective subglacial water flux",
281  "m2 s-1", "m2 s-1", "", 0);
282 
283  m_Qstag_average.set_attrs("internal",
284  "average (over time) advection flux on the staggered grid",
285  "m2 s-1", "m2 s-1", "", 0);
286 
287  m_Vstag.set_attrs("internal",
288  "cell face-centered (staggered) components of water velocity"
289  " in subglacial water layer",
290  "m s-1", "m s-1", "", 0);
291 
292  // auxiliary variables which NEED ghosts
293  m_Wstag.set_attrs("internal",
294  "cell face-centered (staggered) values of water layer thickness",
295  "m", "m", "", 0);
296  m_Wstag.metadata()["valid_min"] = {0.0};
297 
298  m_Kstag.set_attrs("internal",
299  "cell face-centered (staggered) values of nonlinear conductivity",
300  "", "", "", 0);
301  m_Kstag.metadata()["valid_min"] = {0.0};
302 
303  m_R.set_attrs("internal",
304  "work space for modeled subglacial water hydraulic potential",
305  "Pa", "Pa", "", 0);
306 
307  // temporaries during update; do not need ghosts
308  m_Wnew.set_attrs("internal",
309  "new thickness of transportable subglacial water layer during update",
310  "m", "m", "", 0);
311  m_Wnew.metadata()["valid_min"] = {0.0};
312 
313  m_Wtillnew.set_attrs("internal",
314  "new thickness of till (subglacial) water layer during update",
315  "m", "m", "", 0);
316  m_Wtillnew.metadata()["valid_min"] = {0.0};
317 
318  {
319  double alpha = m_config->get_number("hydrology.thickness_power_in_flux");
320  if (alpha < 1.0) {
322  "alpha = %f < 1 which is not allowed", alpha);
323  }
324 
325  if (m_config->get_number("hydrology.tillwat_max") < 0.0) {
327  "hydrology::Routing: hydrology.tillwat_max is negative.\n"
328  "This is not allowed.");
329  }
330  }
331 }
332 
334  m_log->message(2,
335  "* Initializing the routing subglacial hydrology model ...\n");
336 
337  if (m_config->get_flag("hydrology.routing.include_floating_ice")) {
338  m_log->message(2, " ... routing subglacial water under grounded and floating ice.\n");
339  } else {
340  m_log->message(2, " ... routing subglacial water under grounded ice only.\n");
341  }
342 }
343 
344 void Routing::restart_impl(const File &input_file, int record) {
345  Hydrology::restart_impl(input_file, record);
346 
347  m_W.read(input_file, record);
348 
349  regrid("Hydrology", m_W);
350 }
351 
352 void Routing::bootstrap_impl(const File &input_file,
353  const IceModelVec2S &ice_thickness) {
354  Hydrology::bootstrap_impl(input_file, ice_thickness);
355 
356  double bwat_default = m_config->get_number("bootstrapping.defaults.bwat");
357  m_W.regrid(input_file, OPTIONAL, bwat_default);
358 
359  regrid("Hydrology", m_W);
360 }
361 
363  const IceModelVec2S &W,
364  const IceModelVec2S &P) {
365  Hydrology::init_impl(W_till, W, P);
366 
367  m_W.copy_from(W);
368 }
369 
370 void Routing::define_model_state_impl(const File &output) const {
372  m_W.define(output);
373 }
374 
375 void Routing::write_model_state_impl(const File &output) const {
377  m_W.write(output);
378 }
379 
380 //! Returns the (trivial) overburden pressure as the pressure of the transportable water,
381 //! because this is the model.
383  return m_Pover;
384 }
385 
387  return m_Vstag;
388 }
389 
390 
391 //! Average the regular grid water thickness to values at the center of cell edges.
392 /*! Uses mask values to avoid averaging using water thickness values from
393  either ice-free or floating areas. */
395  const IceModelVec2CellType &mask,
396  IceModelVec2Stag &result) {
397 
398  bool include_floating = m_config->get_flag("hydrology.routing.include_floating_ice");
399 
400  IceModelVec::AccessList list{ &mask, &W, &result };
401 
402  for (Points p(*m_grid); p; p.next()) {
403  const int i = p.i(), j = p.j();
404 
405  if (include_floating) {
406  // east
407  if (mask.icy(i, j)) {
408  result(i, j, 0) = mask.icy(i + 1, j) ? 0.5 * (W(i, j) + W(i + 1, j)) : W(i, j);
409  } else {
410  result(i, j, 0) = mask.icy(i + 1, j) ? W(i + 1, j) : 0.0;
411  }
412  // north
413  if (mask.icy(i, j)) {
414  result(i, j, 1) = mask.icy(i, j + 1) ? 0.5 * (W(i, j) + W(i, j + 1)) : W(i, j);
415  } else {
416  result(i, j, 1) = mask.icy(i, j + 1) ? W(i, j + 1) : 0.0;
417  }
418  } else {
419  // east
420  if (mask.grounded_ice(i, j)) {
421  result(i, j, 0) = mask.grounded_ice(i + 1, j) ? 0.5 * (W(i, j) + W(i + 1, j)) : W(i, j);
422  } else {
423  result(i, j, 0) = mask.grounded_ice(i + 1, j) ? W(i + 1, j) : 0.0;
424  }
425  // north
426  if (mask.grounded_ice(i, j)) {
427  result(i, j, 1) = mask.grounded_ice(i, j + 1) ? 0.5 * (W(i, j) + W(i, j + 1)) : W(i, j);
428  } else {
429  result(i, j, 1) = mask.grounded_ice(i, j + 1) ? W(i, j + 1) : 0.0;
430  }
431  }
432  }
433 
434  result.update_ghosts();
435 }
436 
437 
438 //! Compute the nonlinear conductivity at the center of cell edges.
439 /*!
440  Computes
441 
442  \f[ K = K(W, \nabla P, \nabla b) = k W^{\alpha-1} |\nabla R|^{\beta-2} \f]
443 
444  on the staggered grid, where \f$R = P+\rho_w g b\f$. We denote
445 
446  \f[ \Pi = |\nabla R|^2 \f]
447 
448  internally; this is computed on a staggered grid by a Mahaffy-like ([@ref Mahaffy])
449  scheme. This requires \f$R\f$ to be defined on a box stencil of width 1.
450 
451  Also returns the maximum over all staggered points of \f$ K W \f$.
452 */
454  const IceModelVec2S &P,
455  const IceModelVec2S &bed_elevation,
456  IceModelVec2Stag &result,
457  double &KW_max) const {
458  const double
459  k = m_config->get_number("hydrology.hydraulic_conductivity"),
460  alpha = m_config->get_number("hydrology.thickness_power_in_flux"),
461  beta = m_config->get_number("hydrology.gradient_power_in_flux"),
462  betapow = (beta - 2.0) / 2.0;
463 
464  IceModelVec::AccessList list({&result, &W});
465 
466  KW_max = 0.0;
467 
468  if (beta != 2.0) {
469  // Put the squared norm of the gradient of the simplified hydrolic potential (Pi) in
470  // "result"
471  //
472  // FIXME: we don't need to re-compute this during every hydrology time step: the
473  // simplified hydrolic potential does not depend on the water amount and can be
474  // computed *once* in update_impl(), before entering the time-stepping loop
475  {
476  // R <-- P + rhow g b
477  P.add(m_rg, bed_elevation, m_R); // yes, it updates ghosts
478 
479  list.add(m_R);
480  for (Points p(*m_grid); p; p.next()) {
481  const int i = p.i(), j = p.j();
482 
483  double dRdx, dRdy;
484  dRdx = (m_R(i + 1, j) - m_R(i, j)) / m_dx;
485  dRdy = (m_R(i + 1, j + 1) + m_R(i, j + 1) - m_R(i + 1, j - 1) - m_R(i, j - 1)) / (4.0 * m_dy);
486  result(i, j, 0) = dRdx * dRdx + dRdy * dRdy;
487 
488  dRdx = (m_R(i + 1, j + 1) + m_R(i + 1, j) - m_R(i - 1, j + 1) - m_R(i - 1, j)) / (4.0 * m_dx);
489  dRdy = (m_R(i, j + 1) - m_R(i, j)) / m_dy;
490  result(i, j, 1) = dRdx * dRdx + dRdy * dRdy;
491  }
492  }
493 
494  // We regularize negative power |\grad psi|^{beta-2} by adding eps because large
495  // head gradient might be 10^7 Pa per 10^4 m or 10^3 Pa/m.
496  const double eps = beta < 2.0 ? 1.0 : 0.0;
497 
498  for (Points p(*m_grid); p; p.next()) {
499  const int i = p.i(), j = p.j();
500 
501  for (int o = 0; o < 2; ++o) {
502  const double Pi = result(i, j, o);
503 
504  // FIXME: same as Pi above: we don't need to re-compute this each time we make a
505  // short hydrology time step
506  double B = pow(Pi + eps * eps, betapow);
507 
508  result(i, j, o) = k * pow(W(i, j, o), alpha - 1.0) * B;
509 
510  KW_max = std::max(KW_max, result(i, j, o) * W(i, j, o));
511  }
512  }
513  } else {
514  for (Points p(*m_grid); p; p.next()) {
515  const int i = p.i(), j = p.j();
516 
517  for (int o = 0; o < 2; ++o) {
518  result(i, j, o) = k * pow(W(i, j, o), alpha - 1.0);
519 
520  KW_max = std::max(KW_max, result(i, j, o) * W(i, j, o));
521  }
522  }
523  }
524 
525  KW_max = GlobalMax(m_grid->com, KW_max);
526 
527  result.update_ghosts();
528 }
529 
530 
531 //! Compute the wall melt rate which comes from (turbulent) dissipation of flow energy.
532 /*!
533  This code fills `result` with
534  \f[ \frac{m_{wall}}{\rho_w} = - \frac{1}{L \rho_w} \mathbf{q} \cdot \nabla \psi = \left(\frac{k}{L \rho_w}\right) W^\alpha |\nabla R|^\beta \f]
535  where \f$R = P+\rho_w g b\f$.
536 
537  Note that conductivity_staggered() computes the related quantity
538  \f$K = k W^{\alpha-1} |\nabla R|^{\beta-2}\f$ on the staggered grid, but
539  contriving to reuse that code would be inefficient because of the
540  staggered-versus-regular change.
541 
542  At the current state of the code, this is a diagnostic calculation only.
543 */
544 void wall_melt(const Routing &model,
545  const IceModelVec2S &bed_elevation,
546  IceModelVec2S &result) {
547 
548  IceGrid::ConstPtr grid = result.grid();
549 
550  Config::ConstPtr config = grid->ctx()->config();
551 
552  const double
553  k = config->get_number("hydrology.hydraulic_conductivity"),
554  L = config->get_number("constants.fresh_water.latent_heat_of_fusion"),
555  alpha = config->get_number("hydrology.thickness_power_in_flux"),
556  beta = config->get_number("hydrology.gradient_power_in_flux"),
557  g = config->get_number("constants.standard_gravity"),
558  rhow = config->get_number("constants.fresh_water.density"),
559  rg = rhow * g,
560  CC = k / (L * rhow);
561 
562  // FIXME: could be scaled with overall factor hydrology_coefficient_wall_melt ?
563  if (alpha < 1.0) {
565  "alpha = %f < 1 which is not allowed", alpha);
566  }
567 
568  IceModelVec2S R(grid, "R", WITH_GHOSTS);
569 
570  // R <-- P + rhow g b
571  model.subglacial_water_pressure().add(rg, bed_elevation, R);
572  // yes, it updates ghosts
573 
574  IceModelVec2S W(grid, "W", WITH_GHOSTS);
576 
577  IceModelVec::AccessList list{&R, &W, &result};
578 
579  double dx = grid->dx();
580  double dy = grid->dy();
581 
582  for (Points p(*grid); p; p.next()) {
583  const int i = p.i(), j = p.j();
584  double dRdx, dRdy;
585 
586  if (W(i, j) > 0.0) {
587  dRdx = 0.0;
588  if (W(i + 1, j) > 0.0) {
589  dRdx = (R(i + 1, j) - R(i, j)) / (2.0 * dx);
590  }
591  if (W(i - 1, j) > 0.0) {
592  dRdx += (R(i, j) - R(i - 1, j)) / (2.0 * dx);
593  }
594  dRdy = 0.0;
595  if (W(i, j + 1) > 0.0) {
596  dRdy = (R(i, j + 1) - R(i, j)) / (2.0 * dy);
597  }
598  if (W(i, j - 1) > 0.0) {
599  dRdy += (R(i, j) - R(i, j - 1)) / (2.0 * dy);
600  }
601  result(i, j) = CC * pow(W(i, j), alpha) * pow(dRdx * dRdx + dRdy * dRdy, beta/2.0);
602  } else {
603  result(i, j) = 0.0;
604  }
605  }
606 }
607 
608 
609 //! Get the advection velocity V at the center of cell edges.
610 /*!
611  Computes the advection velocity @f$\mathbf{V}@f$ on the staggered
612  (edge-centered) grid. If V = (u, v) in components then we have
613  <code> result(i, j, 0) = u(i+1/2, j) </code> and
614  <code> result(i, j, 1) = v(i, j+1/2) </code>
615 
616  The advection velocity is given by the formula
617 
618  @f[ \mathbf{V} = - K \left(\nabla P + \rho_w g \nabla b\right) @f]
619 
620  where @f$\mathbf{V}@f$ is the water velocity, @f$P@f$ is the water
621  pressure, and @f$b@f$ is the bedrock elevation.
622 
623  If the corresponding staggered grid value of the water thickness is zero then that
624  component of V is set to zero. This does not change the flux value (which would be zero
625  anyway) but it does provide the correct max velocity in the CFL calculation. We assume
626  bed has valid ghosts.
627 */
629  const IceModelVec2S &pressure,
630  const IceModelVec2S &bed,
631  const IceModelVec2Stag &K,
632  const IceModelVec2Int *no_model_mask,
633  IceModelVec2Stag &result) const {
634  IceModelVec2S &P = m_R;
635  P.copy_from(pressure); // yes, it updates ghosts
636 
637  IceModelVec::AccessList list{&P, &W, &K, &bed, &result};
638 
639  for (Points p(*m_grid); p; p.next()) {
640  const int i = p.i(), j = p.j();
641 
642  if (W(i, j, 0) > 0.0) {
643  double
644  P_x = (P(i + 1, j) - P(i, j)) / m_dx,
645  b_x = (bed(i + 1, j) - bed(i, j)) / m_dx;
646  result(i, j, 0) = - K(i, j, 0) * (P_x + m_rg * b_x);
647  } else {
648  result(i, j, 0) = 0.0;
649  }
650 
651  if (W(i, j, 1) > 0.0) {
652  double
653  P_y = (P(i, j + 1) - P(i, j)) / m_dy,
654  b_y = (bed(i, j + 1) - bed(i, j)) / m_dy;
655  result(i, j, 1) = - K(i, j, 1) * (P_y + m_rg * b_y);
656  } else {
657  result(i, j, 1) = 0.0;
658  }
659  }
660 
661  if (no_model_mask) {
662  list.add(*no_model_mask);
663 
664  for (Points p(*m_grid); p; p.next()) {
665  const int i = p.i(), j = p.j();
666 
667  auto M = no_model_mask->star(i, j);
668 
669  if (M.ij or M.e) {
670  result(i, j, 0) = 0.0;
671  }
672 
673  if (M.ij or M.n) {
674  result(i, j, 1) = 0.0;
675  }
676  }
677  }
678 }
679 
680 
681 //! Compute Q = V W at edge-centers (staggered grid) by first-order upwinding.
682 /*!
683  The field W must have valid ghost values, but V does not need them.
684 
685  FIXME: This could be re-implemented using the Koren (1993) flux-limiter.
686 */
688  const IceModelVec2S &W,
689  IceModelVec2Stag &result) const {
690  IceModelVec::AccessList list{&W, &V, &result};
691 
692  assert(W.stencil_width() >= 1);
693 
694  for (Points p(*m_grid); p; p.next()) {
695  const int i = p.i(), j = p.j();
696 
697  result(i, j, 0) = V(i, j, 0) * (V(i, j, 0) >= 0.0 ? W(i, j) : W(i + 1, j));
698  result(i, j, 1) = V(i, j, 1) * (V(i, j, 1) >= 0.0 ? W(i, j) : W(i, j + 1));
699  }
700 
701  result.update_ghosts();
702 }
703 
704 /*!
705  * See equation (51) in Bueler and van Pelt.
706  */
707 double Routing::max_timestep_W_diff(double KW_max) const {
708  double D_max = m_rg * KW_max;
709  double result = 1.0 / (m_dx * m_dx) + 1.0 / (m_dy * m_dy);
710  return 0.25 / (D_max * result);
711 }
712 
713 /*!
714  * See equation (50) in Bueler and van Pelt.
715  */
717  // V could be zero if P is constant and bed is flat
718  auto tmp = absmax(m_Vstag);
719 
720  // add a safety margin
721  double alpha = 0.95;
722  double eps = 1e-6;
723 
724  return alpha * 0.5 / (tmp[0]/m_dx + tmp[1]/m_dy + eps);
725 }
726 
727 
728 //! The computation of Wtillnew, called by update().
729 /*!
730  Does a step of the trivial integration
731  \f[ \frac{\partial W_{till}}{\partial t} = \frac{m}{\rho_w} - C\f]
732 
733  where \f$C=\f$`hydrology_tillwat_decay_rate`. Enforces bounds
734  \f$0 \le W_{till} \le W_{till}^{max}\f$ where the upper bound is
735  `hydrology_tillwat_max`. Here \f$m/\rho_w\f$ is `total_input`.
736 
737  Compare hydrology::NullTransport::update_impl().
738 
739  The current code is not quite "code duplication" because the code here:
740 
741  1. computes `Wtill_new` instead of updating `Wtill` in place;
742  2. uses time steps determined by the rest of the hydrology::Routing model;
743  3. does not check mask because the enforce_bounds() call addresses that.
744 
745  Otherwise this is the same physical model with the same configurable parameters.
746 */
747 void Routing::update_Wtill(double dt,
748  const IceModelVec2S &Wtill,
749  const IceModelVec2S &surface_input_rate,
750  const IceModelVec2S &basal_melt_rate,
751  IceModelVec2S &Wtill_new) {
752  const double
753  tillwat_max = m_config->get_number("hydrology.tillwat_max"),
754  C = m_config->get_number("hydrology.tillwat_decay_rate", "m / second");
755 
756  IceModelVec::AccessList list{&Wtill, &Wtill_new, &basal_melt_rate};
757 
758  bool add_surface_input = m_config->get_flag("hydrology.add_water_input_to_till_storage");
759  if (add_surface_input) {
760  list.add(surface_input_rate);
761  }
762 
763  for (Points p(*m_grid); p; p.next()) {
764  const int i = p.i(), j = p.j();
765 
766  double input_rate = basal_melt_rate(i, j);
767  if (add_surface_input) {
768  input_rate += surface_input_rate(i, j);
769  }
770 
771  Wtill_new(i, j) = clip(Wtill(i, j) + dt * (input_rate - C),
772  0.0, tillwat_max);
773  }
774 }
775 
777  const IceModelVec2S &W,
778  const IceModelVec2Stag &Wstag,
779  const IceModelVec2Stag &K,
780  const IceModelVec2Stag &Q,
781  IceModelVec2S &result) {
782  const double
783  wux = 1.0 / (m_dx * m_dx),
784  wuy = 1.0 / (m_dy * m_dy);
785 
786  IceModelVec::AccessList list{&W, &Wstag, &K, &Q, &result};
787 
788  for (Points p(*m_grid); p; p.next()) {
789  const int i = p.i(), j = p.j();
790 
791  auto q = Q.star(i, j);
792  const double divQ = (q.e - q.w) / m_dx + (q.n - q.s) / m_dy;
793 
794  auto k = K.star(i, j);
795  auto ws = Wstag.star(i, j);
796 
797  const double
798  De = m_rg * k.e * ws.e,
799  Dw = m_rg * k.w * ws.w,
800  Dn = m_rg * k.n * ws.n,
801  Ds = m_rg * k.s * ws.s;
802 
803  auto w = W.star(i, j);
804  const double diffW = (wux * (De * (w.e - w.ij) - Dw * (w.ij - w.w)) +
805  wuy * (Dn * (w.n - w.ij) - Ds * (w.ij - w.s)));
806 
807  result(i, j) = dt * (- divQ + diffW);
808  }
809 }
810 
811 
812 //! The computation of Wnew, called by update().
813 void Routing::update_W(double dt,
814  const IceModelVec2S &surface_input_rate,
815  const IceModelVec2S &basal_melt_rate,
816  const IceModelVec2S &W,
817  const IceModelVec2Stag &Wstag,
818  const IceModelVec2S &Wtill,
819  const IceModelVec2S &Wtill_new,
820  const IceModelVec2Stag &K,
821  const IceModelVec2Stag &Q,
822  IceModelVec2S &W_new) {
823 
825 
826  IceModelVec::AccessList list{&W, &Wtill, &Wtill_new, &surface_input_rate,
827  &basal_melt_rate, &m_flow_change_incremental, &W_new};
828 
829  for (Points p(*m_grid); p; p.next()) {
830  const int i = p.i(), j = p.j();
831 
832  double input_rate = surface_input_rate(i, j) + basal_melt_rate(i, j);
833 
834  double Wtill_change = Wtill_new(i, j) - Wtill(i, j);
835  W_new(i, j) = (W(i, j) + (dt * input_rate - Wtill_change) + m_flow_change_incremental(i, j));
836  }
837 
840  m_input_change.add(dt, basal_melt_rate);
841 }
842 
843 //! Update the model state variables W and Wtill by applying the subglacial hydrology model equations.
844 /*!
845  Runs the hydrology model from time t to time t + dt. Here [t, dt]
846  is generally on the order of months to years. This hydrology model will take its
847  own shorter time steps, perhaps hours to weeks.
848 
849  To update W = `bwat` we call update_W(), and to update Wtill = `tillwat` we
850  call update_Wtill().
851 */
852 void Routing::update_impl(double t, double dt, const Inputs& inputs) {
853 
855 
856  double
857  ht = t,
858  hdt = 0.0;
859 
860  const double
861  t_final = t + dt,
862  dt_max = m_config->get_number("hydrology.maximum_time_step", "seconds"),
863  tillwat_max = m_config->get_number("hydrology.tillwat_max");
864 
865  m_Qstag_average.set(0.0);
866 
867  // make sure W has valid ghosts before starting hydrology steps
868  m_W.update_ghosts();
869 
870  unsigned int step_counter = 0;
871  for (; ht < t_final; ht += hdt) {
872  step_counter++;
873 
874 #if (Pism_DEBUG==1)
875  double huge_number = 1e6;
876  check_bounds(m_W, huge_number);
877  check_bounds(m_Wtill, tillwat_max);
878 #endif
879 
880  // updates ghosts of m_Wstag
882  inputs.geometry->cell_type,
883  m_Wstag);
884 
885  double maxKW = 0.0;
886  // updates ghosts of m_Kstag
887  m_grid->ctx()->profiling().begin("routing_conductivity");
891  m_Kstag, maxKW);
892  m_grid->ctx()->profiling().end("routing_conductivity");
893 
894  // ghosts of m_Vstag are not updated
895  m_grid->ctx()->profiling().begin("routing_velocity");
899  m_Kstag,
900  inputs.no_model_mask,
901  m_Vstag);
902  m_grid->ctx()->profiling().end("routing_velocity");
903 
904  // to get Q, W needs valid ghosts (ghosts of m_Vstag are not used)
905  // updates ghosts of m_Qstag
906  m_grid->ctx()->profiling().begin("routing_flux");
908  m_grid->ctx()->profiling().end("routing_flux");
909 
911 
912  {
913  const double
914  dt_cfl = max_timestep_W_cfl(),
915  dt_diff_w = max_timestep_W_diff(maxKW);
916 
917  hdt = std::min(t_final - ht, dt_max);
918  hdt = std::min(hdt, dt_cfl);
919  hdt = std::min(hdt, dt_diff_w);
920  }
921 
922  m_log->message(3, " hydrology step %05d, dt = %f s\n", step_counter, hdt);
923 
924  // update Wtillnew from Wtill and input_rate
925  {
926  m_grid->ctx()->profiling().begin("routing_Wtill");
927  update_Wtill(hdt,
928  m_Wtill,
931  m_Wtillnew);
932  // remove water in ice-free areas and account for changes
934  inputs.no_model_mask,
935  0.0, // do not limit maximum thickness
936  tillwat_max, // till water thickness under the ocean
937  m_Wtillnew,
942  m_grid->ctx()->profiling().end("routing_Wtill");
943  }
944 
945  // update Wnew from W, Wtill, Wtillnew, Wstag, Q, input_rate
946  // uses ghosts of m_W, m_Wstag, m_Qstag, m_Kstag
947  {
948  m_grid->ctx()->profiling().begin("routing_W");
949  update_W(hdt,
952  m_W, m_Wstag,
954  m_Kstag, m_Qstag,
955  m_Wnew);
956  // remove water in ice-free areas and account for changes
958  inputs.no_model_mask,
959  0.0, // do not limit maximum thickness
960  0.0, // transportable water thickness under the ocean
961  m_Wnew,
966 
967  // transfer new into old (updates ghosts of m_W)
969  m_grid->ctx()->profiling().end("routing_W");
970  }
971 
972  // m_Wtill has no ghosts
974  } // end of the time-stepping loop
975 
977  m_config->get_flag("hydrology.routing.include_floating_ice"),
978  m_Q);
979  m_Q.scale(1.0 / dt);
980 
981  m_log->message(2,
982  " took %d hydrology sub-steps with average dt = %.6f years (%.3f s or %.3f hours)\n",
983  step_counter,
984  units::convert(m_sys, dt / step_counter, "seconds", "years"),
985  dt / step_counter,
986  (dt / step_counter) / 3600.0);
987 }
988 
989 std::map<std::string, Diagnostic::Ptr> Routing::diagnostics_impl() const {
990  using namespace diagnostics;
991 
992  DiagnosticList result = {
993  {"bwatvel", Diagnostic::Ptr(new BasalWaterVelocity(this))},
994  {"bwp", Diagnostic::Ptr(new BasalWaterPressure(this))},
995  {"bwprel", Diagnostic::Ptr(new RelativeBasalWaterPressure(this))},
996  {"effbwp", Diagnostic::Ptr(new EffectiveBasalWaterPressure(this))},
997  {"wallmelt", Diagnostic::Ptr(new WallMelt(this))},
998  {"hydraulic_potential", Diagnostic::Ptr(new HydraulicPotential(this))},
999  };
1000  return combine(result, Hydrology::diagnostics_impl());
1001 }
1002 
1003 std::map<std::string, TSDiagnostic::Ptr> Routing::ts_diagnostics_impl() const {
1004  std::map<std::string, TSDiagnostic::Ptr> result = {
1005  // FIXME: add mass-conservation diagnostics
1006  };
1007  return result;
1008 }
1009 
1010 } // end of namespace hydrology
1011 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:140
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
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
std::shared_ptr< const Config > ConstPtr
const Routing * model
Definition: Diagnostic.hh:166
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:161
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
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
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool icy(int i, int j) const
bool grounded_ice(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
stencils::Star< int > star(int i, int j) const
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)
std::shared_ptr< IceModelVec2S > Ptr
Definition: iceModelVec.hh:341
stencils::Star< double > star(int i, int j) const
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
std::shared_ptr< IceModelVec2Stag > Ptr
Definition: iceModelVec.hh:454
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
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
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
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 add(double alpha, const IceModelVec &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition: iceModelVec.cc:234
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
IceModelVec2S m_flow_change
Definition: Hydrology.hh:201
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_flow_change_incremental
Definition: Hydrology.hh:190
IceModelVec2S m_basal_melt_rate
Definition: Hydrology.hh:187
IceModelVec2S m_W
effective thickness of transportable basal water
Definition: Hydrology.hh:178
const IceModelVec2S & surface_input_rate() const
Definition: Hydrology.cc:522
IceModelVec2S m_grounded_margin_change
Definition: Hydrology.hh:196
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
IceModelVec2S m_Pover
overburden pressure
Definition: Hydrology.hh:181
IceModelVec2S m_Wtill
effective thickness of basal water stored in till
Definition: Hydrology.hh:175
IceModelVec2S m_grounding_line_change
Definition: Hydrology.hh:197
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
const IceModelVec2S & subglacial_water_thickness() const
Return the effective thickness of the transportable basal water layer.
Definition: Hydrology.cc:510
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
The PISM subglacial hydrology model interface.
Definition: Hydrology.hh:114
const Geometry * geometry
Definition: Hydrology.hh:42
const IceModelVec2Int * no_model_mask
Definition: Hydrology.hh:40
void water_thickness_staggered(const IceModelVec2S &W, const IceModelVec2CellType &mask, IceModelVec2Stag &result)
Average the regular grid water thickness to values at the center of cell edges.
Definition: Routing.cc:394
IceModelVec2S m_Wnew
Definition: Routing.hh:127
virtual void init_impl(const IceModelVec2S &W_till, const IceModelVec2S &W, const IceModelVec2S &P)
Definition: Routing.cc:362
virtual void initialization_message() const
Definition: Routing.cc:333
void W_change_due_to_flow(double dt, const IceModelVec2S &W, const IceModelVec2Stag &Wstag, const IceModelVec2Stag &K, const IceModelVec2Stag &Q, IceModelVec2S &result)
Definition: Routing.cc:776
virtual std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
Definition: Routing.cc:989
double max_timestep_W_diff(double KW_max) const
Definition: Routing.cc:707
void compute_conductivity(const IceModelVec2Stag &W, const IceModelVec2S &P, const IceModelVec2S &bed, IceModelVec2Stag &result, double &maxKW) const
Compute the nonlinear conductivity at the center of cell edges.
Definition: Routing.cc:453
IceModelVec2Stag m_Vstag
Definition: Routing.hh:118
IceModelVec2S m_R
Definition: Routing.hh:130
double max_timestep_W_cfl() const
Definition: Routing.cc:716
IceModelVec2S m_Wtillnew
Definition: Routing.hh:127
const IceModelVec2S & subglacial_water_pressure() const
Definition: Routing.cc:382
IceModelVec2S m_bottom_surface
Definition: Routing.hh:135
IceModelVec2Stag m_Kstag
Definition: Routing.hh:124
void compute_velocity(const IceModelVec2Stag &W, const IceModelVec2S &P, const IceModelVec2S &bed, const IceModelVec2Stag &K, const IceModelVec2Int *no_model_mask, IceModelVec2Stag &result) const
Get the advection velocity V at the center of cell edges.
Definition: Routing.cc:628
Routing(IceGrid::ConstPtr g)
Definition: Routing.cc:260
virtual void update_impl(double t, double dt, const Inputs &inputs)
Update the model state variables W and Wtill by applying the subglacial hydrology model equations.
Definition: Routing.cc:852
virtual void bootstrap_impl(const File &input_file, const IceModelVec2S &ice_thickness)
Definition: Routing.cc:352
IceModelVec2Stag m_Qstag_average
Definition: Routing.hh:115
void advective_fluxes(const IceModelVec2Stag &V, const IceModelVec2S &W, IceModelVec2Stag &result) const
Compute Q = V W at edge-centers (staggered grid) by first-order upwinding.
Definition: Routing.cc:687
void update_W(double dt, const IceModelVec2S &surface_input_rate, const IceModelVec2S &basal_melt_rate, const IceModelVec2S &W, const IceModelVec2Stag &Wstag, const IceModelVec2S &Wtill, const IceModelVec2S &Wtill_new, const IceModelVec2Stag &K, const IceModelVec2Stag &Q, IceModelVec2S &W_new)
The computation of Wnew, called by update().
Definition: Routing.cc:813
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Routing.cc:370
IceModelVec2Stag m_Qstag
Definition: Routing.hh:113
const IceModelVec2Stag & velocity_staggered() const
Definition: Routing.cc:386
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Routing.cc:375
virtual std::map< std::string, TSDiagnostic::Ptr > ts_diagnostics_impl() const
Definition: Routing.cc:1003
virtual void restart_impl(const File &input_file, int record)
Definition: Routing.cc:344
void update_Wtill(double dt, const IceModelVec2S &Wtill, const IceModelVec2S &surface_input_rate, const IceModelVec2S &basal_melt_rate, IceModelVec2S &Wtill_new)
The computation of Wtillnew, called by update().
Definition: Routing.cc:747
IceModelVec2Stag m_Wstag
Definition: Routing.hh:121
A subglacial hydrology model which assumes water pressure equals overburden pressure.
Definition: Routing.hh:81
virtual IceModelVec::Ptr compute_impl() const
Definition: Routing.cc:51
Reports the pressure of the transportable water in the subglacial layer.
Definition: Routing.cc:42
virtual IceModelVec::Ptr compute_impl() const
Definition: Routing.cc:180
Diagnostically reports the staggered-grid components of the velocity of the water in the subglacial l...
Definition: Routing.cc:165
virtual IceModelVec::Ptr compute_impl() const
Definition: Routing.cc:115
Reports the effective pressure of the transportable water in the subglacial layer,...
Definition: Routing.cc:104
Report hydraulic potential in the subglacial hydrology system.
Definition: Routing.cc:226
virtual IceModelVec::Ptr compute_impl() const
Definition: Routing.cc:75
Reports the pressure of the transportable water in the subglacial layer as a fraction of the overburd...
Definition: Routing.cc:63
virtual IceModelVec::Ptr compute_impl() const
Definition: Routing.cc:151
Report the wall melt rate from dissipation of the potential energy of the transportable water.
Definition: Routing.cc:140
#define PISM_ERROR_LOCATION
static const double L
Definition: exactTestL.cc:40
void hydraulic_potential(const IceModelVec2S &W, const IceModelVec2S &P, const IceModelVec2S &sea_level, const IceModelVec2S &bed, const IceModelVec2S &ice_thickness, IceModelVec2S &result)
Compute the hydraulic potential.
Definition: Routing.cc:195
void wall_melt(const Routing &model, const IceModelVec2S &bed_elevation, IceModelVec2S &result)
Compute the wall melt rate which comes from (turbulent) dissipation of flow energy.
Definition: Routing.cc:544
static double K(double psi_x, double psi_y, double speed, double epsilon)
void check_bounds(const IceModelVec2S &W, double W_max)
Definition: Hydrology.cc:557
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:72
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double rhow
Definition: exactTestP.cc:41
static const double g
Definition: exactTestP.cc:39
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
T clip(T x, T a, T b)
double absmax(const IceModelVec2S &input)
Finds maximum over all the absolute values in an IceModelVec2S object. Ignores ghosts.
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
@ OPTIONAL
Definition: IO_Flags.hh:70
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result)
Definition: Geometry.cc:223
T combine(const T &a, const T &b)
void staggered_to_regular(const IceModelVec2CellType &cell_type, const IceModelVec2Stag &input, bool include_floating_ice, IceModelVec2S &result)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
const int C[]
Definition: ssafd_code.cc:44