tTemperatureModel.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
---
tTemperatureModel.cc (18463B)
---
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 "TemperatureModel.hh"
21 #include "pism/energy/tempSystem.hh"
22 #include "pism/util/pism_utilities.hh"
23 #include "pism/energy/utilities.hh"
24 #include "pism/util/IceModelVec2CellType.hh"
25 #include "pism/util/Vars.hh"
26 #include "pism/util/io/File.hh"
27
28 namespace pism {
29 namespace energy {
30
31 TemperatureModel::TemperatureModel(IceGrid::ConstPtr grid,
32 stressbalance::StressBalance *stress_balance)
33 : EnergyModel(grid, stress_balance) {
34
35 m_ice_temperature.create(m_grid, "temp", WITH_GHOSTS);
36 m_ice_temperature.set_attrs("model_state",
37 "ice temperature",
38 "K", "K", "land_ice_temperature", 0);
39 m_ice_temperature.metadata().set_number("valid_min", 0.0);
40 }
41
42 const IceModelVec3 & TemperatureModel::temperature() const {
43 return m_ice_temperature;
44 }
45
46 void TemperatureModel::restart_impl(const File &input_file, int record) {
47
48 m_log->message(2, "* Restarting the temperature-based energy balance model from %s...\n",
49 input_file.filename().c_str());
50
51 m_basal_melt_rate.read(input_file, record);
52
53 const IceModelVec2S &ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
54
55 if (input_file.find_variable(m_ice_temperature.metadata().get_name())) {
56 m_ice_temperature.read(input_file, record);
57 } else {
58 init_enthalpy(input_file, false, record);
59 compute_temperature(m_ice_enthalpy, ice_thickness, m_ice_temperature);
60 }
61
62 regrid("Temperature-based energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
63 regrid("Temperature-based energy balance model", m_ice_temperature, REGRID_WITHOUT_REGRID_VARS);
64
65 compute_enthalpy_cold(m_ice_temperature, ice_thickness, m_ice_enthalpy);
66 }
67
68 void TemperatureModel::bootstrap_impl(const File &input_file,
69 const IceModelVec2S &ice_thickness,
70 const IceModelVec2S &surface_temperature,
71 const IceModelVec2S &climatic_mass_balance,
72 const IceModelVec2S &basal_heat_flux) {
73
74 m_log->message(2, "* Bootstrapping the temperature-based energy balance model from %s...\n",
75 input_file.filename().c_str());
76
77 m_basal_melt_rate.regrid(input_file, OPTIONAL,
78 m_config->get_number("bootstrapping.defaults.bmelt"));
79 regrid("Temperature-based energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
80
81 int temp_revision = m_ice_temperature.state_counter();
82 regrid("Temperature-based energy balance model", m_ice_temperature, REGRID_WITHOUT_REGRID_VARS);
83
84 if (temp_revision == m_ice_temperature.state_counter()) {
85 bootstrap_ice_temperature(ice_thickness, surface_temperature, climatic_mass_balance,
86 basal_heat_flux, m_ice_temperature);
87 }
88
89 compute_enthalpy_cold(m_ice_temperature, ice_thickness, m_ice_enthalpy);
90 }
91
92 void TemperatureModel::initialize_impl(const IceModelVec2S &basal_melt_rate,
93 const IceModelVec2S &ice_thickness,
94 const IceModelVec2S &surface_temperature,
95 const IceModelVec2S &climatic_mass_balance,
96 const IceModelVec2S &basal_heat_flux) {
97
98 m_log->message(2, "* Bootstrapping the temperature-based energy balance model...\n");
99
100 m_basal_melt_rate.copy_from(basal_melt_rate);
101 regrid("Temperature-based energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
102
103 int temp_revision = m_ice_temperature.state_counter();
104 regrid("Temperature-based energy balance model", m_ice_temperature, REGRID_WITHOUT_REGRID_VARS);
105
106 if (temp_revision == m_ice_temperature.state_counter()) {
107 bootstrap_ice_temperature(ice_thickness, surface_temperature, climatic_mass_balance,
108 basal_heat_flux, m_ice_temperature);
109 }
110
111 compute_enthalpy_cold(m_ice_temperature, ice_thickness, m_ice_enthalpy);
112 }
113
114 //! Takes a semi-implicit time-step for the temperature equation.
115 /*!
116 This method should be kept because it is worth having alternative physics, and
117 so that older results can be reproduced. In particular, this method is
118 documented by papers [\ref BBL,\ref BBssasliding]. The new browser page
119 \ref bombproofenth essentially documents the cold-ice-BOMBPROOF method here, but
120 the newer enthalpy-based method is slightly different and (we hope) a superior
121 implementation of the conservation of energy principle.
122
123 The conservation of energy equation written in terms of temperature is
124 \f[ \rho c_p(T) \frac{dT}{dt} = k \frac{\partial^2 T}{\partial z^2} + \Sigma,\f]
125 where \f$T(t,x,y,z)\f$ is the temperature of the ice. This equation is the shallow approximation
126 of the full 3D conservation of energy. Note \f$dT/dt\f$ stands for the material
127 derivative, so advection is included. Here \f$\rho\f$ is the density of ice,
128 \f$c_p\f$ is its specific heat, and \f$k\f$ is its conductivity. Also \f$\Sigma\f$ is the volume
129 strain heating (with SI units of \f$J/(\text{s} \text{m}^3) = \text{W}\,\text{m}^{-3}\f$).
130
131 We handle horizontal advection explicitly by first-order upwinding. We handle
132 vertical advection implicitly by centered differencing when possible, and retreat to
133 implicit first-order upwinding when necessary. There is a CFL condition
134 for the horizontal explicit upwinding [\ref MortonMayers]. We report
135 any CFL violations, but they are designed to not occur.
136
137 The vertical conduction term is always handled implicitly (%i.e. by backward Euler).
138
139 We work from the bottom of the column upward in building the system to solve
140 (in the semi-implicit time-stepping scheme). The excess energy above pressure melting
141 is converted to melt-water, and that a fraction of this melt water is transported to
142 the ice base according to the scheme in excessToFromBasalMeltLayer().
143
144 The method uses equally-spaced calculation but the columnSystemCtx
145 methods coarse_to_fine(), fine_to_coarse() interpolate
146 back-and-forth from this equally-spaced computational grid to the
147 (usually) non-equally spaced storage grid.
148
149 An instance of tempSystemCtx is used to solve the tridiagonal system set-up here.
150
151 In this procedure two scalar fields are modified: basal_melt_rate and m_work.
152 But basal_melt_rate will never need to communicate ghosted values (horizontal stencil
153 neighbors). The ghosted values for m_ice_temperature are updated from the values in m_work in the
154 communication done by energy_step().
155
156 The (older) scheme cold-ice-BOMBPROOF, implemented here, is very reliable, but there is
157 still an extreme and rare fjord situation which causes trouble. For example, it
158 occurs in one column of ice in one fjord perhaps only once
159 in a 200ka simulation of the whole sheet, in my (ELB) experience modeling the Greenland
160 ice sheet. It causes the discretized advection bulge to give temperatures below that
161 of the coldest ice anywhere, a continuum impossibility. So as a final protection
162 there is a "bulge limiter" which sets the temperature to the surface temperature of
163 the column minus the bulge maximum (15 K) if it is below that level. The number of
164 times this occurs is reported as a "BPbulge" percentage.
165 */
166 void TemperatureModel::update_impl(double t, double dt, const Inputs &inputs) {
167 // current time does not matter here
168 (void) t;
169
170 using mask::ocean;
171
172 Logger log(MPI_COMM_SELF, m_log->get_threshold());
173
174 const double
175 ice_density = m_config->get_number("constants.ice.density"),
176 ice_c = m_config->get_number("constants.ice.specific_heat_capacity"),
177 L = m_config->get_number("constants.fresh_water.latent_heat_of_fusion"),
178 melting_point_temp = m_config->get_number("constants.fresh_water.melting_point_temperature"),
179 beta_CC_grad = m_config->get_number("constants.ice.beta_Clausius_Clapeyron") * ice_density * m_config->get_number("constants.standard_gravity");
180
181 const bool allow_above_melting = m_config->get_flag("energy.allow_temperature_above_melting");
182
183 // this is bulge limit constant in K; is max amount by which ice
184 // or bedrock can be lower than surface temperature
185 const double bulge_max = m_config->get_number("energy.enthalpy.cold_bulge_max") / ice_c;
186
187 inputs.check();
188 const IceModelVec3
189 &strain_heating3 = *inputs.volumetric_heating_rate,
190 &u3 = *inputs.u3,
191 &v3 = *inputs.v3,
192 &w3 = *inputs.w3;
193
194 const IceModelVec2CellType &cell_type = *inputs.cell_type;
195
196 const IceModelVec2S
197 &basal_frictional_heating = *inputs.basal_frictional_heating,
198 &basal_heat_flux = *inputs.basal_heat_flux,
199 &ice_thickness = *inputs.ice_thickness,
200 &shelf_base_temp = *inputs.shelf_base_temp,
201 &ice_surface_temp = *inputs.surface_temp,
202 &till_water_thickness = *inputs.till_water_thickness;
203
204 IceModelVec::AccessList list{&ice_surface_temp, &shelf_base_temp, &ice_thickness,
205 &cell_type, &basal_heat_flux, &till_water_thickness, &basal_frictional_heating,
206 &u3, &v3, &w3, &strain_heating3, &m_basal_melt_rate, &m_ice_temperature, &m_work};
207
208 energy::tempSystemCtx system(m_grid->z(), "temperature",
209 m_grid->dx(), m_grid->dy(), dt,
210 *m_config,
211 m_ice_temperature, u3, v3, w3, strain_heating3);
212
213 double dz = system.dz();
214 const std::vector<double>& z_fine = system.z();
215 size_t Mz_fine = z_fine.size();
216 std::vector<double> x(Mz_fine);// space for solution of system
217 std::vector<double> Tnew(Mz_fine); // post-processed solution
218
219 // counts unreasonably low temperature values; deprecated?
220 unsigned int maxLowTempCount = m_config->get_number("energy.max_low_temperature_count");
221 const double T_minimum = m_config->get_number("energy.minimum_allowed_temperature");
222
223 double margin_threshold = m_config->get_number("energy.margin_ice_thickness_limit");
224
225 ParallelSection loop(m_grid->com);
226 try {
227 for (Points p(*m_grid); p; p.next()) {
228 const int i = p.i(), j = p.j();
229
230 MaskValue mask = static_cast<MaskValue>(cell_type.as_int(i,j));
231
232 const double H = ice_thickness(i, j);
233 const double T_surface = ice_surface_temp(i, j);
234
235 system.initThisColumn(i, j,
236 marginal(ice_thickness, i, j, margin_threshold),
237 mask, H);
238
239 const int ks = system.ks();
240
241 if (ks > 0) { // if there are enough points in ice to bother ...
242
243 if (system.lambda() < 1.0) {
244 m_stats.reduced_accuracy_counter += 1; // count columns with lambda < 1
245 }
246
247 // set boundary values for tridiagonal system
248 system.setSurfaceBoundaryValuesThisColumn(T_surface);
249 system.setBasalBoundaryValuesThisColumn(basal_heat_flux(i,j),
250 shelf_base_temp(i,j),
251 basal_frictional_heating(i,j));
252
253 // solve the system for this column; melting not addressed yet
254 system.solveThisColumn(x);
255 } // end of "if there are enough points in ice to bother ..."
256
257 // prepare for melting/refreezing
258 double bwatnew = till_water_thickness(i,j);
259
260 // insert solution for generic ice segments
261 for (int k=1; k <= ks; k++) {
262 if (allow_above_melting) { // in the ice
263 Tnew[k] = x[k];
264 } else {
265 const double
266 Tpmp = melting_point_temp - beta_CC_grad * (H - z_fine[k]); // FIXME issue #15
267 if (x[k] > Tpmp) {
268 Tnew[k] = Tpmp;
269 double Texcess = x[k] - Tpmp; // always positive
270 column_drainage(ice_density, ice_c, L, z_fine[k], dz, &Texcess, &bwatnew);
271 // Texcess will always come back zero here; ignore it
272 } else {
273 Tnew[k] = x[k];
274 }
275 }
276 if (Tnew[k] < T_minimum) {
277 log.message(1,
278 " [[too low (<200) ice segment temp T = %f at %d, %d, %d;"
279 " proc %d; mask=%d; w=%f m year-1]]\n",
280 Tnew[k], i, j, k, m_grid->rank(), mask,
281 units::convert(m_sys, system.w(k), "m second-1", "m year-1"));
282
283 m_stats.low_temperature_counter++;
284 }
285 if (Tnew[k] < T_surface - bulge_max) {
286 Tnew[k] = T_surface - bulge_max;
287 m_stats.bulge_counter += 1;
288 }
289 }
290
291 // insert solution for ice base segment
292 if (ks > 0) {
293 if (allow_above_melting == true) { // ice/rock interface
294 Tnew[0] = x[0];
295 } else { // compute diff between x[k0] and Tpmp; melt or refreeze as appropriate
296 const double Tpmp = melting_point_temp - beta_CC_grad * H; // FIXME issue #15
297 double Texcess = x[0] - Tpmp; // positive or negative
298 if (ocean(mask)) {
299 // when floating, only half a segment has had its temperature raised
300 // above Tpmp
301 column_drainage(ice_density, ice_c, L, 0.0, dz/2.0, &Texcess, &bwatnew);
302 } else {
303 column_drainage(ice_density, ice_c, L, 0.0, dz, &Texcess, &bwatnew);
304 }
305 Tnew[0] = Tpmp + Texcess;
306 if (Tnew[0] > (Tpmp + 0.00001)) {
307 throw RuntimeError(PISM_ERROR_LOCATION, "updated temperature came out above Tpmp");
308 }
309 }
310 if (Tnew[0] < T_minimum) {
311 log.message(1,
312 " [[too low (<200) ice/bedrock segment temp T = %f at %d,%d;"
313 " proc %d; mask=%d; w=%f]]\n",
314 Tnew[0],i,j,m_grid->rank(), mask,
315 units::convert(m_sys, system.w(0), "m second-1", "m year-1"));
316
317 m_stats.low_temperature_counter++;
318 }
319 if (Tnew[0] < T_surface - bulge_max) {
320 Tnew[0] = T_surface - bulge_max;
321 m_stats.bulge_counter += 1;
322 }
323 }
324
325 // set to air temp above ice
326 for (unsigned int k = ks; k < Mz_fine; k++) {
327 Tnew[k] = T_surface;
328 }
329
330 // transfer column into m_work; communication later
331 system.fine_to_coarse(Tnew, i, j, m_work);
332
333 // basal_melt_rate(i,j) is rate of mass loss at bottom of ice
334 if (ocean(mask)) {
335 m_basal_melt_rate(i,j) = 0.0;
336 } else {
337 // basalMeltRate is rate of change of bwat; can be negative
338 // (subglacial water freezes-on); note this rate is calculated
339 // *before* limiting or other nontrivial modelling of bwat,
340 // which is Hydrology's job
341 m_basal_melt_rate(i,j) = (bwatnew - till_water_thickness(i,j)) / dt;
342 } // end of the grounded case
343 }
344 } catch (...) {
345 loop.failed();
346 }
347 loop.check();
348
349 m_stats.low_temperature_counter = GlobalSum(m_grid->com, m_stats.low_temperature_counter);
350 if (m_stats.low_temperature_counter > maxLowTempCount) {
351 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "too many low temps: %d",
352 m_stats.low_temperature_counter);
353 }
354
355 // copy to m_ice_temperature, updating ghosts
356 m_work.update_ghosts(m_ice_temperature);
357
358 // Set ice enthalpy in place. EnergyModel::update will scatter ghosts
359 compute_enthalpy_cold(m_work, ice_thickness, m_work);
360 }
361
362 void TemperatureModel::define_model_state_impl(const File &output) const {
363 m_ice_temperature.define(output);
364 m_basal_melt_rate.define(output);
365 // ice enthalpy is not a part of the model state
366 }
367
368 void TemperatureModel::write_model_state_impl(const File &output) const {
369 m_ice_temperature.write(output);
370 m_basal_melt_rate.write(output);
371 // ice enthalpy is not a part of the model state
372 }
373
374 //! Compute the melt water which should go to the base if \f$T\f$ is above pressure-melting.
375 void TemperatureModel::column_drainage(const double rho, const double c, const double L,
376 const double z, const double dz,
377 double *Texcess, double *bwat) const {
378
379 const double
380 darea = m_grid->cell_area(),
381 dvol = darea * dz,
382 dE = rho * c * (*Texcess) * dvol,
383 massmelted = dE / L;
384
385 if (*Texcess >= 0.0) {
386 // T is at or above pressure-melting temp, so temp needs to be set to
387 // pressure-melting; a fraction of excess energy is turned into melt water at base
388 // note massmelted is POSITIVE!
389 const double FRACTION_TO_BASE
390 = (z < 100.0) ? 0.2 * (100.0 - z) / 100.0 : 0.0;
391 // note: ice-equiv thickness:
392 *bwat += (FRACTION_TO_BASE * massmelted) / (rho * darea);
393 *Texcess = 0.0;
394 } else { // neither Texcess nor bwat needs to change if Texcess < 0.0
395 // Texcess negative; only refreeze (i.e. reduce bwat) if at base and bwat > 0.0
396 // note ONLY CALLED IF AT BASE! note massmelted is NEGATIVE!
397 if (z > 0.00001) {
398 throw RuntimeError(PISM_ERROR_LOCATION, "excessToBasalMeltLayer() called with z not at base and negative Texcess");
399 }
400 if (*bwat > 0.0) {
401 const double thicknessToFreezeOn = - massmelted / (rho * darea);
402 if (thicknessToFreezeOn <= *bwat) { // the water *is* available to freeze on
403 *bwat -= thicknessToFreezeOn;
404 *Texcess = 0.0;
405 } else { // only refreeze bwat thickness of water; update Texcess
406 *bwat = 0.0;
407 const double dTemp = L * (*bwat) / (c * dz);
408 *Texcess += dTemp;
409 }
410 }
411 // note: if *bwat == 0 and Texcess < 0.0 then Texcess unmolested; temp will go down
412 }
413 }
414
415 } // end of namespace energy
416 } // end of namespace