tutilities.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
---
tutilities.cc (16942B)
---
1 /* Copyright (C) 2016, 2017, 2018, 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
20 #include "utilities.hh"
21
22 #include "pism/util/IceGrid.hh"
23 #include "pism/util/iceModelVec.hh"
24 #include "pism/util/Logger.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/EnthalpyConverter.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "bootstrapping.hh"
29
30 namespace pism {
31 namespace energy {
32
33 //! Compute ice enthalpy from temperature temperature by assuming the ice has zero liquid fraction.
34 /*!
35 First this method makes sure the temperatures is at most the pressure-melting
36 value, before computing the enthalpy for that temperature, using zero liquid
37 fraction.
38
39 Because of how EnthalpyConverter::pressure() works, the energy
40 content in the air is set to the value that ice would have if it a chunk of it
41 occupied the air; the atmosphere actually has much lower energy content. It is
42 done this way for regularity (i.e. dEnth/dz computations).
43 */
44 void compute_enthalpy_cold(const IceModelVec3 &temperature,
45 const IceModelVec2S &ice_thickness,
46 IceModelVec3 &result) {
47
48 IceGrid::ConstPtr grid = result.grid();
49 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
50
51 IceModelVec::AccessList list{&temperature, &result, &ice_thickness};
52
53 const unsigned int Mz = grid->Mz();
54 const std::vector<double> &z = grid->z();
55
56 for (Points p(*grid); p; p.next()) {
57 const int i = p.i(), j = p.j();
58
59 const double *Tij = temperature.get_column(i,j);
60 double *Enthij = result.get_column(i,j);
61
62 for (unsigned int k = 0; k < Mz; ++k) {
63 const double depth = ice_thickness(i, j) - z[k]; // FIXME issue #15
64 Enthij[k] = EC->enthalpy_permissive(Tij[k], 0.0, EC->pressure(depth));
65 }
66 }
67
68 result.inc_state_counter();
69
70 result.update_ghosts();
71 }
72
73 void compute_temperature(const IceModelVec3 &enthalpy,
74 const IceModelVec2S &ice_thickness,
75 IceModelVec3 &result) {
76
77 IceGrid::ConstPtr grid = result.grid();
78 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
79
80 IceModelVec::AccessList list{&enthalpy, &ice_thickness, &result};
81
82 const unsigned int Mz = grid->Mz();
83 const std::vector<double> &z = grid->z();
84
85 for (Points p(*grid); p; p.next()) {
86 const int i = p.i(), j = p.j();
87
88 const double
89 *E = enthalpy.get_column(i, j),
90 H = ice_thickness(i, j);
91 double *T = result.get_column(i, j);
92
93 for (unsigned int k = 0; k < Mz; ++k) {
94 const double depth = H - z[k]; // FIXME issue #15
95 T[k] = EC->temperature(E[k], EC->pressure(depth));
96 }
97 }
98
99 result.inc_state_counter();
100
101 result.update_ghosts();
102 }
103
104 //! Compute `result` (enthalpy) from `temperature` and liquid fraction.
105 void compute_enthalpy(const IceModelVec3 &temperature,
106 const IceModelVec3 &liquid_water_fraction,
107 const IceModelVec2S &ice_thickness,
108 IceModelVec3 &result) {
109
110 IceGrid::ConstPtr grid = result.grid();
111 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
112
113 IceModelVec::AccessList list{&temperature, &liquid_water_fraction, &ice_thickness, &result};
114
115 const unsigned int Mz = grid->Mz();
116 const std::vector<double> &z = grid->z();
117
118 for (Points p(*grid); p; p.next()) {
119 const int i = p.i(), j = p.j();
120
121 const double *T = temperature.get_column(i,j);
122 const double *omega = liquid_water_fraction.get_column(i,j);
123 double *E = result.get_column(i,j);
124
125 for (unsigned int k = 0; k < Mz; ++k) {
126 const double depth = ice_thickness(i,j) - z[k]; // FIXME issue #15
127 E[k] = EC->enthalpy_permissive(T[k], omega[k], EC->pressure(depth));
128 }
129 }
130
131 result.update_ghosts();
132
133 result.inc_state_counter();
134 }
135
136 //! Compute the liquid fraction corresponding to enthalpy and ice_thickness.
137 void compute_liquid_water_fraction(const IceModelVec3 &enthalpy,
138 const IceModelVec2S &ice_thickness,
139 IceModelVec3 &result) {
140
141 IceGrid::ConstPtr grid = result.grid();
142
143 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
144
145 result.set_name("liqfrac");
146 result.metadata(0).set_name("liqfrac");
147 result.set_attrs("diagnostic",
148 "liquid water fraction in ice (between 0 and 1)",
149 "1", "1", "", 0);
150
151 IceModelVec::AccessList list{&result, &enthalpy, &ice_thickness};
152
153 ParallelSection loop(grid->com);
154 try {
155 for (Points p(*grid); p; p.next()) {
156 const int i = p.i(), j = p.j();
157
158 const double *Enthij = enthalpy.get_column(i,j);
159 double *omegaij = result.get_column(i,j);
160
161 for (unsigned int k=0; k < grid->Mz(); ++k) {
162 const double depth = ice_thickness(i,j) - grid->z(k); // FIXME issue #15
163 omegaij[k] = EC->water_fraction(Enthij[k],EC->pressure(depth));
164 }
165 }
166 } catch (...) {
167 loop.failed();
168 }
169 loop.check();
170
171 result.inc_state_counter();
172 }
173
174 //! @brief Compute the CTS field, CTS = E/E_s(p), from `ice_enthalpy` and `ice_thickness`, and put
175 //! in `result`.
176 /*!
177 * The actual cold-temperate transition surface (CTS) is the level set CTS = 1.
178 *
179 * Does not communicate ghosts for IceModelVec3 result.
180 */
181 void compute_cts(const IceModelVec3 &ice_enthalpy,
182 const IceModelVec2S &ice_thickness,
183 IceModelVec3 &result) {
184
185 IceGrid::ConstPtr grid = result.grid();
186 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
187
188 result.set_name("cts");
189 result.metadata(0).set_name("cts");
190 result.set_attrs("diagnostic",
191 "cts = E/E_s(p), so cold-temperate transition surface is at cts = 1",
192 "1", "1", "", 0);
193
194 IceModelVec::AccessList list{&ice_enthalpy, &ice_thickness, &result};
195
196 const unsigned int Mz = grid->Mz();
197 const std::vector<double> &z = grid->z();
198
199 for (Points p(*grid); p; p.next()) {
200 const int i = p.i(), j = p.j();
201
202 double *CTS = result.get_column(i,j);
203 const double *enthalpy = ice_enthalpy.get_column(i,j);
204
205 for (unsigned int k = 0; k < Mz; ++k) {
206 const double depth = ice_thickness(i,j) - z[k]; // FIXME issue #15
207 CTS[k] = enthalpy[k] / EC->enthalpy_cts(EC->pressure(depth));
208 }
209 }
210
211 result.inc_state_counter();
212 }
213
214 //! Computes the total ice enthalpy in J.
215 /*!
216 Units of the specific enthalpy field are J kg-1. We integrate
217 \f$E(t,x,y,z)\f$ over the entire ice fluid region \f$\Omega(t)\f$, multiplying
218 by the density to get units of energy:
219 \f[ E_{\text{total}}(t) = \int_{\Omega(t)} E(t,x,y,z) \rho_i \,dx\,dy\,dz. \f]
220 */
221 double total_ice_enthalpy(double thickness_threshold,
222 const IceModelVec3 &ice_enthalpy,
223 const IceModelVec2S &ice_thickness) {
224 double enthalpy_sum = 0.0;
225
226 IceGrid::ConstPtr grid = ice_enthalpy.grid();
227 Config::ConstPtr config = grid->ctx()->config();
228
229 auto cell_area = grid->cell_area();
230
231 const std::vector<double> &z = grid->z();
232
233 IceModelVec::AccessList list{&ice_enthalpy, &ice_thickness};
234 ParallelSection loop(grid->com);
235 try {
236 for (Points p(*grid); p; p.next()) {
237 const int i = p.i(), j = p.j();
238
239 const double H = ice_thickness(i, j);
240
241 if (H >= thickness_threshold) {
242 const int ks = grid->kBelowHeight(H);
243
244 const double
245 *E = ice_enthalpy.get_column(i, j);
246
247 for (int k = 0; k < ks; ++k) {
248 enthalpy_sum += cell_area * E[k] * (z[k+1] - z[k]);
249 }
250 enthalpy_sum += cell_area * E[ks] * (H - z[ks]);
251 }
252 }
253 } catch (...) {
254 loop.failed();
255 }
256 loop.check();
257
258 enthalpy_sum *= config->get_number("constants.ice.density");
259
260 return GlobalSum(grid->com, enthalpy_sum);
261 }
262
263 //! Create a temperature field within the ice from provided ice thickness, surface temperature, surface mass balance, and geothermal flux.
264 /*!
265 In bootstrapping we need to determine initial values for the temperature within
266 the ice (and the bedrock). There are various data available at bootstrapping,
267 but not the 3D temperature field needed as initial values for the temperature. Here
268 we take a "guess" based on an assumption of steady state and a simple model of
269 the vertical velocity in the column. The rule is certainly heuristic but it
270 seems to work well anyway.
271
272 The result is *not* the temperature field which is in steady state with the ice
273 dynamics. Spinup is most-definitely needed in many applications. Such spinup
274 usually starts from the temperature field computed by this procedure and then
275 runs for a long time (e.g. \f$10^4\f$ to \f$10^6\f$ years), possibly with fixed
276 geometry, to get closer to thermomechanically-coupled equilibrium.
277
278 Consider a horizontal grid point. Suppose the surface temperature
279 \f$T_s\f$, surface mass balance \f$m\f$, and geothermal flux \f$g\f$ are given at that location.
280 Within the column denote the temperature by \f$T(z)\f$ at height \f$z\f$ above
281 the base of the ice. Suppose the column of ice has height \f$H\f$, the ice
282 thickness.
283
284 There are two alternative bootstrap methods determined by the configuration parameter
285 `config.get_number("bootstrapping.temperature_heuristic"))`. Allowed values are `"smb"` and
286 `"quartic_guess"`.
287
288 1. If the `smb` method is chosen, which is the default, and if \f$m>0\f$,
289 then the method sets the ice
290 temperature to the solution of the steady problem [\ref Paterson]
291 \f[\rho_i c w \frac{\partial T}{\partial z} = k_i \frac{\partial^2 T}{\partial z^2} \qquad \text{with boundary conditions} \qquad T(H) = T_s \quad \text{and} \quad \frac{\partial T}{\partial z}(0) = - \frac{g}{k_i}, \f]
292 where the vertical velocity is linear between the surface value \f$w=-m\f$ and
293 a velocity of zero at the base:
294 \f[w(z) = - m z / H.\f]
295 (Note that because \f$m>0\f$, this vertical velocity is downward.)
296 This is a two-point boundary value problem for a linear ODE. In fact, if
297 \f$K = k_i / (\rho_i c)\f$ then we can write the ODE as
298 \f[K T'' + \frac{m z}{H} T' = 0.\f]
299 Then let
300 \f[C_0 = \frac{g \sqrt{\pi H K}}{k_i \sqrt{2 m}}, \qquad \gamma_0 = \sqrt{\frac{mH}{2K}}.\f]
301 (Note \f$\gamma_0\f$ is, up to a constant, the square root of the Peclet number
302 [\ref Paterson]; compare [\ref vanderWeletal2013].) The solution to the
303 two-point boundary value problem is then
304 \f[T(z) = T_s + C_0 \left(\operatorname{erf}(\gamma_0) - \operatorname{erf}\left(\gamma_0 \frac{z}{H}\right)\right).\f]
305 If `usesmb` is true and \f$m \le 0\f$, then the velocity in the column, relative
306 to the base, is taken to be zero. Thus the solution is
307 \f[ T(z) = \frac{g}{k_i} \left( H - z \right) + T_s, \f]
308 a straight line whose slope is determined by the geothermal flux and whose value
309 at the ice surface is the surface temperature, \f$T(H) = T_s\f$.
310 2. If the `quartic_guess` method is chosen, the "quartic guess" formula which was in older
311 versions of PISM is used. Namely, within the ice we set
312 \f[T(z) = T_s + \alpha (H-z)^2 + \beta (H-z)^4\f]
313 where \f$\alpha,\beta\f$ are chosen so that
314 \f[\frac{\partial T}{\partial z}\Big|_{z=0} = - \frac{g}{k_i} \qquad \text{and} \qquad \frac{\partial T}{\partial z}\Big|_{z=H/4} = - \frac{g}{2 k_i}.\f]
315 The purpose of the second condition is that when ice is advecting downward then
316 the temperature gradient is much larger in roughly the bottom quarter of the
317 ice column. However, without the surface mass balance, much less the solution
318 of the stress balance equations, we cannot estimate the vertical velocity, so
319 we make such a rough guess.
320
321 In either case the temperature within the ice is not allowed to exceed the
322 pressure-melting temperature.
323
324 We set \f$T(z)=T_s\f$ above the top of the ice.
325
326 This method determines \f$T(0)\f$, the ice temperature at the ice base. This
327 temperature is used by BedThermalUnit::bootstrap() to determine a
328 bootstrap temperature profile in the bedrock.
329 */
330 void bootstrap_ice_temperature(const IceModelVec2S &ice_thickness,
331 const IceModelVec2S &ice_surface_temp,
332 const IceModelVec2S &surface_mass_balance,
333 const IceModelVec2S &basal_heat_flux,
334 IceModelVec3 &result) {
335
336 IceGrid::ConstPtr grid = result.grid();
337 Context::ConstPtr ctx = grid->ctx();
338 Config::ConstPtr config = ctx->config();
339 Logger::ConstPtr log = ctx->log();
340 EnthalpyConverter::Ptr EC = ctx->enthalpy_converter();
341
342 const bool use_smb = config->get_string("bootstrapping.temperature_heuristic") == "smb";
343
344 if (use_smb) {
345 log->message(2,
346 " - Filling 3D ice temperatures using surface temperature"
347 " (and mass balance for velocity estimate)...\n");
348
349 } else {
350 log->message(2,
351 " - Filling 3D ice temperatures using surface temperature"
352 " (and a quartic guess without SMB)...\n");
353 }
354
355 const double
356 ice_k = config->get_number("constants.ice.thermal_conductivity"),
357 ice_density = config->get_number("constants.ice.density"),
358 ice_c = config->get_number("constants.ice.specific_heat_capacity"),
359 K = ice_k / (ice_density * ice_c),
360 T_min = config->get_number("energy.minimum_allowed_temperature"),
361 T_melting = config->get_number("constants.fresh_water.melting_point_temperature",
362 "Kelvin");
363
364 IceModelVec::AccessList list{&ice_surface_temp, &surface_mass_balance,
365 &ice_thickness, &basal_heat_flux, &result};
366
367 ParallelSection loop(grid->com);
368 try {
369 for (Points p(*grid); p; p.next()) {
370 const int i = p.i(), j = p.j();
371
372 const double
373 T_surface = std::min(ice_surface_temp(i, j), T_melting),
374 H = ice_thickness(i, j),
375 G = basal_heat_flux(i, j);
376
377 if (G < 0.0) {
378 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
379 "geothermal flux G(%d,%d) = %f < 0.0 %s",
380 i, j, G,
381 basal_heat_flux.metadata().get_string("units").c_str());
382 }
383
384 if (T_surface < T_min) {
385 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
386 "T_surface(%d,%d) = %f < T_min = %f Kelvin",
387 i, j, T_surface, T_min);
388 }
389
390 const unsigned int ks = grid->kBelowHeight(H);
391
392 double *T = result.get_column(i, j);
393
394 // within ice
395 if (use_smb) { // method 1: includes surface mass balance in estimate
396
397 // Convert SMB from "kg m-2 s-1" to "m second-1".
398 const double SMB = surface_mass_balance(i, j) / ice_density;
399
400 for (unsigned int k = 0; k < ks; k++) {
401 const double z = grid->z(k);
402 T[k] = ice_temperature_guess_smb(EC, H, z, T_surface, G, ice_k, K, SMB);
403 }
404
405 } else { // method 2: a quartic guess; does not use SMB
406
407 for (unsigned int k = 0; k < ks; k++) {
408 const double z = grid->z(k);
409 T[k] = ice_temperature_guess(EC, H, z, T_surface, G, ice_k);
410 }
411
412 }
413
414 // Make sure that resulting temperatures are not too low.
415 for (unsigned int k = 0; k < ks; k++) {
416 if (T[k] < T_min) {
417 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
418 "T(%d,%d,%d) = %f < T_min = %f Kelvin",
419 i, j, k, T[k], T_min);
420 }
421 }
422
423 // above ice
424 for (unsigned int k = ks; k < grid->Mz(); k++) {
425 T[k] = T_surface;
426 }
427 }
428 } catch (...) {
429 loop.failed();
430 }
431 loop.check();
432
433 result.update_ghosts();
434 }
435
436 void bootstrap_ice_enthalpy(const IceModelVec2S &ice_thickness,
437 const IceModelVec2S &ice_surface_temp,
438 const IceModelVec2S &surface_mass_balance,
439 const IceModelVec2S &basal_heat_flux,
440 IceModelVec3 &result) {
441
442 bootstrap_ice_temperature(ice_thickness, ice_surface_temp,
443 surface_mass_balance, basal_heat_flux,
444 result);
445
446 compute_enthalpy_cold(result, ice_thickness, result);
447 }
448
449 } // end of namespace energy
450 } // end of namespace pism