tflux_balance.hh - 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
---
tflux_balance.hh (17491B)
---
1 /* Copyright (C) 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 /*
21 * This source code is included in diagnostics.cc and so it does not need any #include
22 * directives here.
23 */
24
25 namespace pism {
26 namespace diagnostics {
27
28 enum AmountKind {AMOUNT, MASS};
29
30 //! @brief Computes tendency_of_ice_amount, the ice amount rate of change.
31 class TendencyOfIceAmount : public Diag<IceModel>
32 {
33 public:
34 TendencyOfIceAmount(const IceModel *m, AmountKind kind)
35 : Diag<IceModel>(m),
36 m_kind(kind),
37 m_last_amount(m_grid, "last_ice_amount", WITHOUT_GHOSTS),
38 m_interval_length(0.0) {
39
40 std::string
41 name = "tendency_of_ice_amount",
42 long_name = "rate of change of the ice amount",
43 standard_name = "",
44 internal_units = "kg m-2 second-1",
45 external_units = "kg m-2 year-1";
46 if (kind == MASS) {
47 name = "tendency_of_ice_mass";
48 long_name = "rate of change of the ice mass";
49 internal_units = "kg second-1";
50 external_units = "Gt year-1" ;
51 }
52
53 // set metadata:
54 m_vars = {SpatialVariableMetadata(m_sys, name)};
55
56 set_attrs(long_name, standard_name, internal_units, external_units, 0);
57
58 auto large_number = to_internal(1e6);
59
60 m_vars[0].set_numbers("valid_range", {-large_number, large_number});
61 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
62 m_vars[0].set_string("cell_methods", "time: mean");
63
64 auto units = internal_units + " second";
65 m_last_amount.set_attrs("internal",
66 "ice amount at the time of the last report of " + name,
67 units, units, "", 0);
68 }
69 protected:
70 IceModelVec::Ptr compute_impl() const {
71
72 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "", WITHOUT_GHOSTS));
73 result->metadata() = m_vars[0];
74
75 if (m_interval_length > 0.0) {
76 double ice_density = m_config->get_number("constants.ice.density");
77
78 auto cell_area = m_grid->cell_area();
79
80 const IceModelVec2S& thickness = model->geometry().ice_thickness;
81 const IceModelVec2S& area_specific_volume = model->geometry().ice_area_specific_volume;
82
83 IceModelVec::AccessList list{result.get(),
84 &thickness, &area_specific_volume, &m_last_amount};
85
86 for (Points p(*m_grid); p; p.next()) {
87 const int i = p.i(), j = p.j();
88
89 // m * (kg / m^3) = kg / m^2
90 double amount = (thickness(i, j) + area_specific_volume(i, j)) * ice_density;
91
92 (*result)(i, j) = (amount - m_last_amount(i, j)) / m_interval_length;
93
94 if (m_kind == MASS) {
95 // kg / m^2 * m^2 = kg
96 (*result)(i, j) *= cell_area;
97 }
98 }
99 } else {
100 result->set(m_fill_value);
101 }
102
103 return result;
104 }
105
106 void reset_impl() {
107 m_interval_length = 0.0;
108
109 const IceModelVec2S& thickness = model->geometry().ice_thickness;
110 const IceModelVec2S& area_specific_volume = model->geometry().ice_area_specific_volume;
111
112 double ice_density = m_config->get_number("constants.ice.density");
113
114 IceModelVec::AccessList list{&m_last_amount, &thickness, &area_specific_volume};
115
116 for (Points p(*m_grid); p; p.next()) {
117 const int i = p.i(), j = p.j();
118
119 // m * (kg / m^3) = kg / m^2
120 m_last_amount(i, j) = (thickness(i, j) + area_specific_volume(i, j)) * ice_density;
121 }
122 }
123
124 void update_impl(double dt) {
125 m_interval_length += dt;
126 }
127
128 protected:
129 AmountKind m_kind;
130 IceModelVec2S m_last_amount;
131 double m_interval_length;
132 };
133
134 //! @brief Computes tendency_of_ice_amount_due_to_flow, the rate of change of ice amount due to
135 //! flow.
136 /*! @brief Report rate of change of ice amount due to flow. */
137 class TendencyOfIceAmountDueToFlow : public DiagAverageRate<IceModel>
138 {
139 public:
140 TendencyOfIceAmountDueToFlow(const IceModel *m, AmountKind kind)
141 : DiagAverageRate<IceModel>(m,
142 kind == AMOUNT
143 ? "tendency_of_ice_amount_due_to_flow"
144 : "tendency_of_ice_mass_due_to_flow", TOTAL_CHANGE),
145 m_kind(kind) {
146
147 std::string
148 name = "tendency_of_ice_amount_due_to_flow",
149 long_name = "rate of change of ice amount due to flow",
150 standard_name = "",
151 accumulator_units = "kg m-2",
152 internal_units = "kg m-2 second-1",
153 external_units = "kg m-2 year-1";
154
155 if (kind == MASS) {
156 name = "tendency_of_ice_mass_due_to_flow";
157 long_name = "rate of change of ice mass due to flow";
158 standard_name = "";
159 accumulator_units = "kg";
160 internal_units = "kg second-1";
161 external_units = "Gt year-1";
162 }
163
164 m_factor = m_config->get_number("constants.ice.density");
165
166 m_vars = {SpatialVariableMetadata(m_sys, name)};
167 m_accumulator.metadata().set_string("units", accumulator_units);
168
169 set_attrs(long_name, standard_name, internal_units, external_units, 0);
170 m_vars[0].set_string("cell_methods", "time: mean");
171 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
172 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
173 }
174
175 protected:
176 void update_impl(double dt) {
177 const IceModelVec2S
178 &dH = model->geometry_evolution().thickness_change_due_to_flow(),
179 &dV = model->geometry_evolution().area_specific_volume_change_due_to_flow();
180
181 auto cell_area = m_grid->cell_area();
182
183 IceModelVec::AccessList list{&m_accumulator, &dH, &dV};
184
185 for (Points p(*m_grid); p; p.next()) {
186 const int i = p.i(), j = p.j();
187
188 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
189
190 m_accumulator(i, j) += C * (dH(i, j) + dV(i, j));
191 }
192
193 m_interval_length += dt;
194 }
195 AmountKind m_kind;
196 };
197
198 /*! @brief Report surface mass balance flux, averaged over the reporting interval */
199 class SurfaceFlux : public DiagAverageRate<IceModel>
200 {
201 public:
202 SurfaceFlux(const IceModel *m, AmountKind kind)
203 : DiagAverageRate<IceModel>(m,
204 kind == AMOUNT
205 ? "tendency_of_ice_amount_due_to_surface_mass_flux"
206 : "tendency_of_ice_mass_due_to_surface_mass_flux",
207 TOTAL_CHANGE),
208 m_kind(kind) {
209 m_factor = m_config->get_number("constants.ice.density");
210
211 auto ismip6 = m_config->get_flag("output.ISMIP6");
212
213 std::string
214 name = ismip6 ? "acabf" : "tendency_of_ice_amount_due_to_surface_mass_flux",
215 accumulator_units = "kg m-2",
216 long_name = "average surface mass flux over reporting interval",
217 standard_name = "land_ice_surface_specific_mass_balance_flux",
218 internal_units = "kg m-2 s-1",
219 external_units = "kg m-2 year-1";
220 if (kind == MASS) {
221 name = "tendency_of_ice_mass_due_to_surface_mass_flux",
222 accumulator_units = "kg",
223 long_name = "average surface mass flux over reporting interval",
224 standard_name = "",
225 internal_units = "kg second-1",
226 external_units = "Gt year-1";
227 }
228
229 m_vars = {SpatialVariableMetadata(m_sys, name)};
230 m_accumulator.metadata().set_string("units", accumulator_units);
231
232 set_attrs(long_name, standard_name, internal_units, external_units, 0);
233
234 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
235 m_vars[0].set_string("cell_methods", "time: mean");
236 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
237 }
238
239 protected:
240 void update_impl(double dt) {
241 const IceModelVec2S
242 &SMB = model->geometry_evolution().top_surface_mass_balance();
243
244 auto cell_area = m_grid->cell_area();
245
246 IceModelVec::AccessList list{&m_accumulator, &SMB};
247
248 for (Points p(*m_grid); p; p.next()) {
249 const int i = p.i(), j = p.j();
250
251 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
252
253 m_accumulator(i, j) += C * SMB(i, j);
254 }
255
256 m_interval_length += dt;
257 }
258 AmountKind m_kind;
259 };
260
261 /*! @brief Report basal mass balance flux, averaged over the reporting interval */
262 class BasalFlux : public DiagAverageRate<IceModel>
263 {
264 public:
265 BasalFlux(const IceModel *m, AmountKind kind)
266 : DiagAverageRate<IceModel>(m,
267 kind == AMOUNT
268 ? "tendency_of_ice_amount_due_to_basal_mass_flux"
269 : "tendency_of_ice_mass_due_to_basal_mass_flux",
270 TOTAL_CHANGE),
271 m_kind(kind) {
272 m_factor = m_config->get_number("constants.ice.density");
273
274 std::string
275 name = "tendency_of_ice_amount_due_to_basal_mass_flux",
276 accumulator_units = "kg m-2",
277 long_name = "average basal mass flux over reporting interval",
278 standard_name = "",
279 internal_units = "kg m-2 second-1",
280 external_units = "kg m-2 year-1";
281 if (kind == MASS) {
282 name = "tendency_of_ice_mass_due_to_basal_mass_flux",
283 accumulator_units = "kg",
284 long_name = "average basal mass flux over reporting interval",
285 standard_name = "tendency_of_land_ice_mass_due_to_basal_mass_balance",
286 internal_units = "kg second-1",
287 external_units = "Gt year-1";
288 }
289 m_vars = {SpatialVariableMetadata(m_sys, name)};
290 m_accumulator.metadata().set_string("units", accumulator_units);
291
292 set_attrs(long_name, standard_name, internal_units, external_units, 0);
293
294 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
295 m_vars[0].set_string("cell_methods", "time: mean");
296 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
297 }
298
299 protected:
300 void update_impl(double dt) {
301 const IceModelVec2S
302 &BMB = model->geometry_evolution().bottom_surface_mass_balance();
303
304 auto cell_area = m_grid->cell_area();
305
306 IceModelVec::AccessList list{&m_accumulator, &BMB};
307
308 for (Points p(*m_grid); p; p.next()) {
309 const int i = p.i(), j = p.j();
310
311 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
312
313 m_accumulator(i, j) += C * BMB(i, j);
314 }
315
316 m_interval_length += dt;
317 }
318 AmountKind m_kind;
319 };
320
321 class ConservationErrorFlux : public DiagAverageRate<IceModel>
322 {
323 public:
324 ConservationErrorFlux(const IceModel *m, AmountKind kind)
325 : DiagAverageRate<IceModel>(m,
326 kind == AMOUNT
327 ? "tendency_of_ice_amount_due_to_conservation_error"
328 : "tendency_of_ice_mass_due_to_conservation_error" ,
329 TOTAL_CHANGE),
330 m_kind(kind) {
331 m_factor = m_config->get_number("constants.ice.density");
332
333 std::string
334 name = "tendency_of_ice_amount_due_to_conservation_error",
335 accumulator_units = "kg m-2",
336 long_name = "average mass conservation error flux over reporting interval",
337 standard_name = "",
338 internal_units = "kg m-2 second-1",
339 external_units = "kg m-2 year-1";
340 if (kind == MASS) {
341 name = "tendency_of_ice_mass_due_to_conservation_error",
342 accumulator_units = "kg",
343 long_name = "average mass conservation error flux over reporting interval",
344 standard_name = "",
345 internal_units = "kg second-1",
346 external_units = "Gt year-1";
347 }
348
349 m_vars = {SpatialVariableMetadata(m_sys, name)};
350 m_accumulator.metadata().set_string("units", accumulator_units);
351
352 set_attrs(long_name, standard_name, internal_units, external_units, 0);
353
354 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
355 m_vars[0].set_string("cell_methods", "time: mean");
356 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
357 }
358
359 protected:
360 void update_impl(double dt) {
361 const IceModelVec2S
362 &error = model->geometry_evolution().conservation_error();
363
364 IceModelVec::AccessList list{&m_accumulator, &error};
365
366 auto cell_area = m_grid->cell_area();
367
368 for (Points p(*m_grid); p; p.next()) {
369 const int i = p.i(), j = p.j();
370
371 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
372
373 m_accumulator(i, j) += C * error(i, j);
374 }
375
376 m_interval_length += dt;
377 }
378 AmountKind m_kind;
379 };
380
381 /*! @brief Report discharge (calving and frontal melt) flux. */
382 class DischargeFlux : public DiagAverageRate<IceModel>
383 {
384 public:
385 DischargeFlux(const IceModel *m, AmountKind kind)
386 : DiagAverageRate<IceModel>(m,
387 kind == AMOUNT
388 ? "tendency_of_ice_amount_due_to_discharge"
389 : "tendency_of_ice_mass_due_to_discharge",
390 TOTAL_CHANGE),
391 m_kind(kind) {
392
393 m_factor = m_config->get_number("constants.ice.density");
394
395 auto ismip6 = m_config->get_flag("output.ISMIP6");
396
397 std::string
398 name = ismip6 ? "lifmassbf" : "tendency_of_ice_amount_due_to_discharge",
399 long_name = "discharge flux (calving, frontal melt, forced retreat)",
400 accumulator_units = "kg m-2",
401 standard_name = "land_ice_specific_mass_flux_due_to_calving_and_ice_front_melting",
402 internal_units = "kg m-2 s-1",
403 external_units = "kg m-2 year-1";
404 if (kind == MASS) {
405 name = "tendency_of_ice_mass_due_to_discharge";
406 long_name = "discharge flux (calving, frontal melt, forced retreat)";
407 accumulator_units = "kg";
408 standard_name = "";
409 internal_units = "kg second-1";
410 external_units = "Gt year-1";
411 }
412
413 m_vars = {SpatialVariableMetadata(m_sys, name)};
414 m_accumulator.metadata().set_string("units", accumulator_units);
415
416 set_attrs(long_name, standard_name, internal_units, external_units, 0);
417 m_vars[0].set_string("cell_methods", "time: mean");
418
419 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
420 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
421 }
422
423 protected:
424 void update_impl(double dt) {
425 const IceModelVec2S &calving = model->calving();
426 const IceModelVec2S &frontal_melt = model->frontal_melt();
427 const IceModelVec2S &forced_retreat = model->forced_retreat();
428
429 IceModelVec::AccessList list{&m_accumulator, &calving, &frontal_melt, &forced_retreat};
430
431 auto cell_area = m_grid->cell_area();
432
433 for (Points p(*m_grid); p; p.next()) {
434 const int i = p.i(), j = p.j();
435
436 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
437
438 m_accumulator(i, j) += C * (calving(i, j) + frontal_melt(i, j) + forced_retreat(i, j));
439 }
440
441 m_interval_length += dt;
442 }
443 AmountKind m_kind;
444 };
445
446 /*! @brief Report discharge (calving and frontal melt) flux. */
447 class CalvingFlux : public DiagAverageRate<IceModel>
448 {
449 public:
450 CalvingFlux(const IceModel *m, AmountKind kind)
451 : DiagAverageRate<IceModel>(m,
452 kind == AMOUNT
453 ? "tendency_of_ice_amount_due_to_calving"
454 : "tendency_of_ice_mass_due_to_calving",
455 TOTAL_CHANGE),
456 m_kind(kind) {
457
458 m_factor = m_config->get_number("constants.ice.density");
459
460 auto ismip6 = m_config->get_flag("output.ISMIP6");
461
462 std::string
463 name = ismip6 ? "licalvf" : "tendency_of_ice_amount_due_to_calving",
464 long_name = "calving flux",
465 accumulator_units = "kg m-2",
466 standard_name = "land_ice_specific_mass_flux_due_to_calving",
467 internal_units = "kg m-2 s-1",
468 external_units = "kg m-2 year-1";
469 if (kind == MASS) {
470 name = "tendency_of_ice_mass_due_to_calving";
471 long_name = "calving flux";
472 accumulator_units = "kg";
473 standard_name = "";
474 internal_units = "kg second-1";
475 external_units = "Gt year-1";
476 }
477
478 m_vars = {SpatialVariableMetadata(m_sys, name)};
479 m_accumulator.metadata().set_string("units", accumulator_units);
480
481 set_attrs(long_name, standard_name, internal_units, external_units, 0);
482 m_vars[0].set_string("cell_methods", "time: mean");
483
484 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
485 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
486 }
487
488 protected:
489 void update_impl(double dt) {
490 const IceModelVec2S &calving = model->calving();
491
492 IceModelVec::AccessList list{&m_accumulator, &calving};
493
494 auto cell_area = m_grid->cell_area();
495
496 for (Points p(*m_grid); p; p.next()) {
497 const int i = p.i(), j = p.j();
498
499 double C = m_factor * (m_kind == AMOUNT ? 1.0 : cell_area);
500
501 m_accumulator(i, j) += C * calving(i, j);
502 }
503
504 m_interval_length += dt;
505 }
506 AmountKind m_kind;
507 };
508
509
510 } // end of namespace diagnostics
511 } // end of namespace pism