tNullTransport.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tNullTransport.cc (6822B)
---
1 // Copyright (C) 2012-2019 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 "NullTransport.hh"
20 #include "pism/util/error_handling.hh"
21 #include "pism/util/MaxTimestep.hh"
22 #include "pism/util/IceModelVec2CellType.hh"
23 #include "pism/util/pism_utilities.hh" // clip
24 #include "pism/geometry/Geometry.hh"
25
26 namespace pism {
27 namespace hydrology {
28
29 NullTransport::NullTransport(IceGrid::ConstPtr g)
30 : Hydrology(g) {
31 m_diffuse_tillwat = m_config->get_flag("hydrology.null_diffuse_till_water");
32 m_diffusion_time = m_config->get_number("hydrology.null_diffusion_time", "seconds");
33 m_diffusion_distance = m_config->get_number("hydrology.null_diffusion_distance", "meters");
34 m_tillwat_max = m_config->get_number("hydrology.tillwat_max", "meters");
35 m_tillwat_decay_rate = m_config->get_number("hydrology.tillwat_decay_rate", "m / second");
36
37 if (m_tillwat_max < 0.0) {
38 throw RuntimeError(PISM_ERROR_LOCATION,
39 "hydrology::NullTransport: hydrology.tillwat_max is negative.\n"
40 "This is not allowed.");
41 }
42
43 if (m_diffuse_tillwat) {
44 m_Wtill_old.create(m_grid, "Wtill_old", WITH_GHOSTS);
45 }
46 }
47
48 NullTransport::~NullTransport() {
49 // empty
50 }
51
52 void NullTransport::initialization_message() const {
53 m_log->message(2,
54 "* Initializing the null-transport (till only) subglacial hydrology model ...\n");
55
56 if (m_diffuse_tillwat) {
57 m_log->message(2,
58 " [using lateral diffusion of stored till water as in Bueler and Brown, 2009]\n");
59 }
60 }
61
62 void NullTransport::restart_impl(const File &input_file, int record) {
63 Hydrology::restart_impl(input_file, record);
64 }
65
66 void NullTransport::bootstrap_impl(const File &input_file,
67 const IceModelVec2S &ice_thickness) {
68 Hydrology::bootstrap_impl(input_file, ice_thickness);
69 }
70
71 void NullTransport::init_impl(const IceModelVec2S &W_till,
72 const IceModelVec2S &W,
73 const IceModelVec2S &P) {
74 Hydrology::init_impl(W_till, W, P);
75 }
76
77 MaxTimestep NullTransport::max_timestep_impl(double t) const {
78 (void) t;
79 if (m_diffuse_tillwat) {
80 const double
81 dx2 = m_grid->dx() * m_grid->dx(),
82 dy2 = m_grid->dy() * m_grid->dy(),
83 L = m_diffusion_distance,
84 T = m_diffusion_time,
85 K = L * L / (2.0 * T);
86
87 return MaxTimestep(dx2 * dy2 / (2.0 * K * (dx2 + dy2)), "null-transport hydrology");
88 } else {
89 return MaxTimestep("null-transport hydrology");
90 }
91 }
92
93 //! Update the till water thickness by simply integrating the melt input.
94 /*!
95 Does a step of the trivial integration
96
97 \f[ \frac{\partial W_{till}}{\partial t} = \frac{m}{\rho_w} - C\f]
98
99 where \f$C=\f$`hydrology_tillwat_decay_rate`. Enforces bounds
100 \f$0 \le W_{till} \le W_{till}^{max}\f$ where the upper bound is
101 `hydrology_tillwat_max`. Here \f$m/\rho_w\f$ is `total_input`.
102
103 Uses the current mass-continuity timestep `dt`. (Compare
104 hydrology::Routing::update_Wtill() which will generally be taking time steps
105 determined by the evolving transportable water layer in that model.)
106
107 There is no attempt to report on conservation errors because this
108 hydrology::NullTransport model does not conserve water.
109
110 There is no tranportable water thickness variable and no interaction with it.
111 */
112 void NullTransport::update_impl(double t, double dt, const Inputs& inputs) {
113 (void) t;
114
115 bool add_surface_input = m_config->get_flag("hydrology.add_water_input_to_till_storage");
116
117 // no transportable basal water
118 m_W.set(0.0);
119
120 m_input_change.add(dt, m_basal_melt_rate);
121 if (add_surface_input) {
122 m_input_change.add(dt, m_surface_input_rate);
123 }
124
125 const double
126 water_density = m_config->get_number("constants.fresh_water.density"),
127 kg_per_m = m_grid->cell_area() * water_density; // kg m-1
128
129 const IceModelVec2CellType &cell_type = inputs.geometry->cell_type;
130
131 IceModelVec::AccessList list{&cell_type, &m_Wtill, &m_basal_melt_rate,
132 &m_conservation_error_change};
133
134 if (add_surface_input) {
135 list.add(m_surface_input_rate);
136 }
137
138 for (Points p(*m_grid); p; p.next()) {
139 const int i = p.i(), j = p.j();
140
141 double
142 W_old = m_Wtill(i, j),
143 dW_input = dt * m_basal_melt_rate(i, j);
144
145 if (add_surface_input) {
146 dW_input += dt * m_surface_input_rate(i, j);
147 }
148
149 if (W_old < 0.0) {
150 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
151 "negative subglacial water thickness of %f m at (%d, %d)",
152 W_old, i, j);
153 }
154
155 // decay rate in areas under grounded ice
156 double dW_decay = dt * (- m_tillwat_decay_rate);
157
158 if (not cell_type.grounded_ice(i, j)) {
159 dW_decay = 0.0;
160 }
161
162 m_Wtill(i, j) = (W_old + dW_input) + dW_decay;
163
164 // dW_decay is a "conservation error"
165 m_conservation_error_change(i, j) += dW_decay * kg_per_m;
166 }
167
168 if (m_diffuse_tillwat) {
169 diffuse_till_water(dt);
170 }
171
172 // remove water in ice-free areas and account for changes
173 enforce_bounds(inputs.geometry->cell_type,
174 inputs.no_model_mask,
175 m_tillwat_max,
176 m_Wtill,
177 m_grounded_margin_change,
178 m_grounding_line_change,
179 m_conservation_error_change,
180 m_no_model_mask_change);
181 }
182
183 void NullTransport::diffuse_till_water(double dt) {
184 // note: this call updates ghosts of m_Wtill_old
185 m_Wtill_old.copy_from(m_Wtill);
186
187 const double
188 dx = m_grid->dx(),
189 dy = m_grid->dy(),
190 L = m_diffusion_distance,
191 T = m_diffusion_time,
192 K = L * L / (2.0 * T),
193 Rx = K * dt / (dx * dx),
194 Ry = K * dt / (dy * dy);
195
196 IceModelVec::AccessList list{&m_Wtill, &m_Wtill_old, &m_flow_change};
197 for (Points p(*m_grid); p; p.next()) {
198 const int i = p.i(), j = p.j();
199
200 auto W = m_Wtill_old.star(i, j);
201
202 double dWtill = (Rx * (W.w - 2.0 * W.ij + W.e) + Ry * (W.s - 2.0 * W.ij + W.n));
203
204 m_Wtill(i, j) = W.ij + dWtill;
205
206 m_flow_change(i, j) = dWtill;
207 }
208 }
209
210 } // end of namespace hydrology
211 } // end of namespace pism