PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Distributed.cc
Go to the documentation of this file.
1 // Copyright (C) 2012-2019, 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 <algorithm> // std::min, std::max
20 
21 #include "Distributed.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/io/File.hh"
26 #include "pism/util/pism_options.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/IceModelVec2CellType.hh"
29 #include "pism/geometry/Geometry.hh"
30 
31 namespace pism {
32 namespace hydrology {
33 
35  : Routing(g),
36  m_P(m_grid, "bwp", WITH_GHOSTS, 1),
37  m_Pnew(m_grid, "Pnew_internal", WITHOUT_GHOSTS) {
38 
39  // additional variables beyond hydrology::Routing
40  m_P.set_attrs("model_state",
41  "pressure of transportable water in subglacial layer",
42  "Pa", "Pa", "", 0);
43  m_P.metadata()["valid_min"] = {0.0};
44 
45  m_Pnew.set_attrs("internal",
46  "new transportable subglacial water pressure during update",
47  "Pa", "Pa", "", 0);
48  m_Pnew.metadata()["valid_min"] = {0.0};
49 }
50 
52  m_log->message(2,
53  "* Initializing the distributed, linked-cavities subglacial hydrology model...\n");
54 }
55 
56 void Distributed::restart_impl(const File &input_file, int record) {
57  Routing::restart_impl(input_file, record);
58 
59  m_P.read(input_file, record);
60 
61  regrid("Hydrology", m_P);
62 }
63 
64 void Distributed::bootstrap_impl(const File &input_file,
65  const IceModelVec2S &ice_thickness) {
66  Routing::bootstrap_impl(input_file, ice_thickness);
67 
68  double bwp_default = m_config->get_number("bootstrapping.defaults.bwp");
69  m_P.regrid(input_file, OPTIONAL, bwp_default);
70 
71  regrid("Hydrology", m_P);
72 
73  bool init_P_from_steady = m_config->get_flag("hydrology.distributed.init_p_from_steady");
74 
75  if (init_P_from_steady) { // if so, just overwrite -i or -bootstrap value of P=bwp
76  m_log->message(2,
77  " initializing P from P(W) formula which applies in steady state\n");
78 
79  compute_overburden_pressure(ice_thickness, m_Pover);
80 
81  IceModelVec2S sliding_speed(m_grid, "velbase_mag", WITHOUT_GHOSTS);
82  sliding_speed.set_attrs("internal", "basal sliding speed",
83  "m s-1", "m s-1", "", 0);
84 
85  std::string filename = m_config->get_string("hydrology.distributed.sliding_speed_file");
86 
87  if (filename.empty()) {
89  "hydrology.distributed.sliding_speed_file is not set");
90  }
91 
92  sliding_speed.regrid(filename, CRITICAL);
93 
94  P_from_W_steady(m_W, m_Pover, sliding_speed,
95  m_P);
96  }
97 }
98 
100  const IceModelVec2S &W,
101  const IceModelVec2S &P) {
102  Routing::init_impl(W_till, W, P);
103 
104  m_P.copy_from(P);
105 }
106 
107 void Distributed::define_model_state_impl(const File &output) const {
109  m_P.define(output);
110 }
111 
112 void Distributed::write_model_state_impl(const File &output) const {
114  m_P.write(output);
115 }
116 
117 std::map<std::string, TSDiagnostic::Ptr> Distributed::ts_diagnostics_impl() const {
118  std::map<std::string, TSDiagnostic::Ptr> result = {
119  // FIXME: add mass-conservation time-series diagnostics
120  };
121  return result;
122 }
123 
124 //! Copies the P state variable which is the modeled water pressure.
126  return m_P;
127 }
128 
129 //! Check bounds on P and fail with message if not satisfied. Optionally, enforces the
130 //! upper bound instead of checking it.
131 /*!
132  The bounds are \f$0 \le P \le P_o\f$ where \f$P_o\f$ is the overburden pressure.
133 */
135  const IceModelVec2S &P_o,
136  bool enforce_upper) {
137 
138  IceModelVec::AccessList list{&P, &P_o};
139 
140  ParallelSection loop(m_grid->com);
141  try {
142  for (Points p(*m_grid); p; p.next()) {
143  const int i = p.i(), j = p.j();
144 
145  if (P(i,j) < 0.0) {
147  "negative subglacial water pressure\n"
148  "P(%d, %d) = %.6f Pa",
149  i, j, P(i, j));
150  }
151 
152  if (enforce_upper) {
153  P(i,j) = std::min(P(i,j), P_o(i,j));
154  } else if (P(i,j) > P_o(i,j) + 0.001) {
156  "subglacial water pressure P = %.16f Pa exceeds\n"
157  "overburden pressure Po = %.16f Pa at (i,j)=(%d,%d)",
158  P(i, j), P_o(i, j), i, j);
159  }
160  }
161  } catch (...) {
162  loop.failed();
163  }
164  loop.check();
165 
166 }
167 
168 
169 //! Compute functional relationship P(W) which applies only in steady state.
170 /*!
171  In steady state in this model, water pressure is determined by a balance of
172  cavitation (opening) caused by sliding and creep closure.
173 
174  This will be used in initialization when P is otherwise unknown, and
175  in verification and/or reporting. It is not used during time-dependent
176  model runs. To be more complete, \f$P = P(W,P_o,|v_b|)\f$.
177 */
179  const IceModelVec2S &P_overburden,
180  const IceModelVec2S &sliding_speed,
181  IceModelVec2S &result) {
182 
183  const double
184  ice_softness = m_config->get_number("flow_law.isothermal_Glen.ice_softness"),
185  creep_closure_coefficient = m_config->get_number("hydrology.creep_closure_coefficient"),
186  cavitation_opening_coefficient = m_config->get_number("hydrology.cavitation_opening_coefficient"),
187  Glen_exponent = m_config->get_number("stress_balance.sia.Glen_exponent"),
188  Wr = m_config->get_number("hydrology.roughness_scale");
189 
190  const double CC = cavitation_opening_coefficient / (creep_closure_coefficient * ice_softness);
191 
192  IceModelVec::AccessList list{&W, &P_overburden, &sliding_speed, &result};
193 
194  for (Points p(*m_grid); p; p.next()) {
195  const int i = p.i(), j = p.j();
196 
197  double sb = pow(CC * sliding_speed(i, j), 1.0 / Glen_exponent);
198  if (W(i, j) == 0.0) {
199  // see P(W) formula in steady state; note P(W) is continuous (in steady
200  // state); these facts imply:
201  if (sb > 0.0) {
202  // no water + cavitation = underpressure
203  result(i, j) = 0.0;
204  } else {
205  // no water + no cavitation = creep repressurizes = overburden
206  result(i, j) = P_overburden(i, j);
207  }
208  } else {
209  double Wratio = std::max(0.0, Wr - W(i, j)) / W(i, j);
210  // in cases where steady state is actually possible this will come out positive, but
211  // otherwise we should get underpressure P=0, and that is what it yields
212  result(i, j) = std::max(0.0, P_overburden(i, j) - sb * pow(Wratio, 1.0 / Glen_exponent));
213  }
214  } // end of the loop over grid points
215 }
216 
217 double Distributed::max_timestep_P_diff(double phi0, double dt_diff_w) const {
218  return 2.0 * phi0 * dt_diff_w;
219 }
220 
221 void Distributed::update_P(double dt,
222  const IceModelVec2CellType &cell_type,
223  const IceModelVec2S &sliding_speed,
224  const IceModelVec2S &surface_input_rate,
225  const IceModelVec2S &basal_melt_rate,
226  const IceModelVec2S &P_overburden,
227  const IceModelVec2S &Wtill,
228  const IceModelVec2S &Wtill_new,
229  const IceModelVec2S &P,
230  const IceModelVec2S &W,
231  const IceModelVec2Stag &Ws,
232  const IceModelVec2Stag &K,
233  const IceModelVec2Stag &Q,
234  IceModelVec2S &P_new) const {
235 
236  const double
237  n = m_config->get_number("stress_balance.sia.Glen_exponent"),
238  A = m_config->get_number("flow_law.isothermal_Glen.ice_softness"),
239  c1 = m_config->get_number("hydrology.cavitation_opening_coefficient"),
240  c2 = m_config->get_number("hydrology.creep_closure_coefficient"),
241  Wr = m_config->get_number("hydrology.roughness_scale"),
242  phi0 = m_config->get_number("hydrology.regularizing_porosity");
243 
244  // update Pnew from time step
245  const double
246  CC = (m_rg * dt) / phi0,
247  wux = 1.0 / (m_dx * m_dx),
248  wuy = 1.0 / (m_dy * m_dy);
249 
250  IceModelVec::AccessList list{&P, &W, &Wtill, &Wtill_new, &sliding_speed, &Ws,
251  &K, &Q, &surface_input_rate, &basal_melt_rate,
252  &cell_type, &P_overburden, &P_new};
253 
254  for (Points p(*m_grid); p; p.next()) {
255  const int i = p.i(), j = p.j();
256 
257  auto w = W.star(i, j);
258  double P_o = P_overburden(i, j);
259 
260  if (cell_type.ice_free_land(i, j)) {
261  P_new(i, j) = 0.0;
262  } else if (cell_type.ocean(i, j)) {
263  P_new(i, j) = P_o;
264  } else if (w.ij <= 0.0) {
265  P_new(i, j) = P_o;
266  } else {
267  auto q = Q.star(i, j);
268  auto k = K.star(i, j);
269  auto ws = Ws.star(i, j);
270 
271  double
272  Open = c1 * sliding_speed(i, j) * std::max(0.0, Wr - w.ij),
273  Close = c2 * A * pow(P_o - P(i, j), n) * w.ij;
274 
275  // compute the flux divergence the same way as in update_W()
276  const double divadflux = (q.e - q.w) / m_dx + (q.n - q.s) / m_dy;
277  const double
278  De = m_rg * k.e * ws.e,
279  Dw = m_rg * k.w * ws.w,
280  Dn = m_rg * k.n * ws.n,
281  Ds = m_rg * k.s * ws.s;
282 
283  double diffW = (wux * (De * (w.e - w.ij) - Dw * (w.ij - w.w)) +
284  wuy * (Dn * (w.n - w.ij) - Ds * (w.ij - w.s)));
285 
286  double divflux = -divadflux + diffW;
287 
288  // pressure update equation
289  double Wtill_change = Wtill_new(i, j) - Wtill(i, j);
290  double total_input = surface_input_rate(i, j) + basal_melt_rate(i, j);
291  double ZZ = Close - Open + total_input - Wtill_change / dt;
292 
293  P_new(i, j) = P(i, j) + CC * (divflux + ZZ);
294 
295  // projection to enforce 0 <= P <= P_o
296  P_new(i, j) = clip(P_new(i, j), 0.0, P_o);
297  }
298  } // end of the loop over grid points
299 }
300 
301 
302 //! Update the model state variables W,P by running the subglacial hydrology model.
303 /*!
304  Runs the hydrology model from time t to time t + dt. Here [t,dt]
305  is generally on the order of months to years. This hydrology model will take its
306  own shorter time steps, perhaps hours to weeks.
307 */
308 void Distributed::update_impl(double t, double dt, const Inputs& inputs) {
309 
311 
312  double
313  ht = t,
314  hdt = 0.0;
315 
316  const double
317  t_final = t + dt,
318  dt_max = m_config->get_number("hydrology.maximum_time_step", "seconds"),
319  phi0 = m_config->get_number("hydrology.regularizing_porosity"),
320  tillwat_max = m_config->get_number("hydrology.tillwat_max");
321 
322  m_Qstag_average.set(0.0);
323 
324  // make sure W,P have valid ghosts before starting hydrology steps
325  m_W.update_ghosts();
326  m_P.update_ghosts();
327 
328  unsigned int step_counter = 0;
329  for (; ht < t_final; ht += hdt) {
330  step_counter++;
331 
332 #if (Pism_DEBUG==1)
333  double huge_number = 1e6;
334  check_bounds(m_W, huge_number);
335  check_bounds(m_Wtill, tillwat_max);
336 #endif
337 
338  // note that ice dynamics can change overburden pressure, so we can only check P
339  // bounds if thk has not changed; if thk could have just changed, such as in the first
340  // time through the current loop, we enforce them
341  bool enforce_upper = (step_counter == 1);
342  check_P_bounds(m_P, m_Pover, enforce_upper);
343 
345  inputs.geometry->cell_type,
346  m_Wstag);
347 
348  double maxKW = 0.0;
352  m_Kstag, maxKW);
353 
357  m_Kstag,
358  inputs.no_model_mask,
359  m_Vstag);
360 
361  // to get Q, W needs valid ghosts
363 
365 
366  {
367  const double
368  dt_cfl = max_timestep_W_cfl(),
369  dt_diff_w = max_timestep_W_diff(maxKW),
370  dt_diff_p = max_timestep_P_diff(phi0, dt_diff_w);
371 
372  hdt = std::min(t_final - ht, dt_max);
373  hdt = std::min(hdt, dt_cfl);
374  hdt = std::min(hdt, dt_diff_w);
375  hdt = std::min(hdt, dt_diff_p);
376  }
377 
378  m_log->message(3, " hydrology step %05d, dt = %f s\n", step_counter, hdt);
379 
380  // update Wtillnew from Wtill and input_rate
381  update_Wtill(hdt,
382  m_Wtill,
385  m_Wtillnew);
386  // remove water in ice-free areas and account for changes
388  inputs.no_model_mask,
389  0.0, // do not limit maximum thickness
390  tillwat_max, // till water thickness under the ocean
391  m_Wtillnew,
396 
397  update_P(hdt,
398  inputs.geometry->cell_type,
399  *inputs.ice_sliding_speed,
402  m_Pover,
405  m_W, m_Wstag,
406  m_Kstag, m_Qstag,
407  m_Pnew);
408 
409  // update Wnew from W, Wtill, Wtillnew, Wstag, Q, input_rate
410  update_W(hdt,
413  m_W, m_Wstag,
415  m_Kstag, m_Qstag,
416  m_Wnew);
417  // remove water in ice-free areas and account for changes
419  inputs.no_model_mask,
420  0.0, // do not limit maximum thickness
421  0.0, // transportable water layer thickness under the ocean
422  m_Wnew,
427 
428  // transfer new into old
432  } // end of the time-stepping loop
433 
435  m_config->get_flag("hydrology.routing.include_floating_ice"),
436  m_Q);
437  m_Q.scale(1.0 / dt);
438 
439  m_log->message(2,
440  " took %d hydrology sub-steps with average dt = %.6f years (%.6f s)\n",
441  step_counter,
442  units::convert(m_sys, dt/step_counter, "seconds", "years"),
443  dt/step_counter);
444 }
445 
446 } // end of namespace hydrology
447 } // 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
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 ice_free_land(int i, int j) const
bool ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
void copy_from(const IceModelVec2S &source)
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.
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
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 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
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
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
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Distributed.cc:112
void P_from_W_steady(const IceModelVec2S &W, const IceModelVec2S &P_overburden, const IceModelVec2S &sliding_speed, IceModelVec2S &result)
Compute functional relationship P(W) which applies only in steady state.
Definition: Distributed.cc:178
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Distributed.cc:107
virtual void bootstrap_impl(const File &input_file, const IceModelVec2S &ice_thickness)
Definition: Distributed.cc:64
void check_P_bounds(IceModelVec2S &P, const IceModelVec2S &P_o, bool enforce_upper)
Definition: Distributed.cc:134
void update_P(double dt, const IceModelVec2CellType &cell_type, const IceModelVec2S &sliding_speed, const IceModelVec2S &surface_input_rate, const IceModelVec2S &basal_melt_rate, const IceModelVec2S &P_overburden, const IceModelVec2S &Wtill, const IceModelVec2S &Wtill_new, const IceModelVec2S &P, const IceModelVec2S &W, const IceModelVec2Stag &Ws, const IceModelVec2Stag &K, const IceModelVec2Stag &Q, IceModelVec2S &P_new) const
Definition: Distributed.cc:221
const IceModelVec2S & subglacial_water_pressure() const
Copies the P state variable which is the modeled water pressure.
Definition: Distributed.cc:125
std::map< std::string, TSDiagnostic::Ptr > ts_diagnostics_impl() const
Definition: Distributed.cc:117
void initialization_message() const
Definition: Distributed.cc:51
Distributed(IceGrid::ConstPtr g)
Definition: Distributed.cc:34
void update_impl(double t, double dt, const Inputs &inputs)
Update the model state variables W,P by running the subglacial hydrology model.
Definition: Distributed.cc:308
virtual double max_timestep_P_diff(double phi0, double dt_diff_w) const
Definition: Distributed.cc:217
virtual void init_impl(const IceModelVec2S &W_till, const IceModelVec2S &W, const IceModelVec2S &P)
Definition: Distributed.cc:99
virtual void restart_impl(const File &input_file, int record)
Definition: Distributed.cc:56
IceModelVec2S m_surface_input_rate
Definition: Hydrology.hh:184
IceModelVec2S m_no_model_mask_change
Definition: Hydrology.hh:199
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
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
void compute_overburden_pressure(const IceModelVec2S &ice_thickness, IceModelVec2S &result) const
Update the overburden pressure from ice thickness.
Definition: Hydrology.cc:483
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 Geometry * geometry
Definition: Hydrology.hh:42
const IceModelVec2Int * no_model_mask
Definition: Hydrology.hh:40
const IceModelVec2S * ice_sliding_speed
Definition: Hydrology.hh:46
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
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
double max_timestep_W_cfl() const
Definition: Routing.cc:716
IceModelVec2S m_Wtillnew
Definition: Routing.hh:127
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
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
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Routing.cc:375
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
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
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 g
Definition: exactTestP.cc:39
T clip(T x, T a, T b)
static const double c2
Definition: exactTestP.cc:48
static const double k
Definition: exactTestP.cc:45
@ OPTIONAL
Definition: IO_Flags.hh:70
@ CRITICAL
Definition: IO_Flags.hh:70
static const double Wr
Definition: exactTestP.cc:46
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double c1
Definition: exactTestP.cc:47
void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result)
Definition: Geometry.cc:223
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