tdiagnostics.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
---
tdiagnostics.cc (105809B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Constantine Khroulev
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 <cassert>
20 #include <algorithm>
21
22 #include "pism/icemodel/IceModel.hh"
23 #include "pism/age/AgeModel.hh"
24 #include "pism/energy/EnergyModel.hh"
25 #include "pism/energy/utilities.hh"
26 #include "pism/geometry/grounded_cell_fraction.hh"
27 #include "pism/geometry/part_grid_threshold_thickness.hh"
28 #include "pism/rheology/FlowLaw.hh"
29 #include "pism/stressbalance/StressBalance.hh"
30 #include "pism/stressbalance/SSB_Modifier.hh"
31 #include "pism/stressbalance/ShallowStressBalance.hh"
32 #include "pism/util/EnthalpyConverter.hh"
33 #include "pism/util/Diagnostic.hh"
34 #include "pism/util/Vars.hh"
35 #include "pism/util/error_handling.hh"
36 #include "pism/util/iceModelVec3Custom.hh"
37 #include "pism/util/pism_utilities.hh"
38 #include "pism/util/projection.hh"
39 #include "pism/earth/BedDef.hh"
40
41 #if (Pism_USE_PROJ==1)
42 #include "pism/util/Proj.hh"
43 #endif
44
45 #include "flux_balance.hh"
46
47 namespace pism {
48
49 // Horrendous names used by InitMIP (and ISMIP6, and CMIP5). Ugh.
50 static const char* land_ice_area_fraction_name = "sftgif";
51 static const char* grounded_ice_sheet_area_fraction_name = "sftgrf";
52 static const char* floating_ice_sheet_area_fraction_name = "sftflf";
53
54 namespace diagnostics {
55
56 enum AreaType {GROUNDED, SHELF, BOTH};
57
58 enum TermType {SMB, BMB, FLOW, ERROR};
59
60 /*! @brief Ocean pressure difference at calving fronts. Used to debug CF boundary conditins. */
61 class IceMarginPressureDifference : public Diag<IceModel>
62 {
63 public:
64 IceMarginPressureDifference(IceModel *m);
65 protected:
66 IceModelVec::Ptr compute_impl() const;
67 };
68
69 IceMarginPressureDifference::IceMarginPressureDifference(IceModel *m)
70 : Diag<IceModel>(m) {
71
72 /* set metadata: */
73 m_vars = {SpatialVariableMetadata(m_sys, "ice_margin_pressure_difference")};
74 m_vars[0].set_number("_FillValue", m_fill_value);
75
76 set_attrs("vertically-integrated pressure difference"
77 " at ice margins (including calving fronts)", "",
78 "Pa m", "Pa m", 0);
79 }
80
81 IceModelVec::Ptr IceMarginPressureDifference::compute_impl() const {
82
83 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "ice_margin_pressure_difference", WITHOUT_GHOSTS));
84 result->metadata(0) = m_vars[0];
85
86 IceModelVec2CellType mask(m_grid, "mask", WITH_GHOSTS);
87
88 auto
89 &H = model->geometry().ice_thickness,
90 &bed = model->geometry().bed_elevation,
91 &sea_level = model->geometry().sea_level_elevation;
92
93 {
94 const double H_threshold = m_config->get_number("stress_balance.ice_free_thickness_standard");
95 GeometryCalculator gc(*m_config);
96 gc.set_icefree_thickness(H_threshold);
97
98 gc.compute_mask(sea_level, bed, H, mask);
99 }
100
101 const double
102 rho_ice = m_config->get_number("constants.ice.density"),
103 rho_ocean = m_config->get_number("constants.sea_water.density"),
104 g = m_config->get_number("constants.standard_gravity");
105
106 const bool dry_mode = m_config->get_flag("ocean.always_grounded");
107
108 IceModelVec::AccessList list{&H, &bed, &mask, &sea_level, result.get()};
109
110 ParallelSection loop(m_grid->com);
111 try {
112 for (Points p(*m_grid); p; p.next()) {
113 const int i = p.i(), j = p.j();
114
115 double delta_p = 0.0;
116 if (mask.grounded_ice(i, j) and grid_edge(*m_grid, i, j)) {
117 delta_p = 0.0;
118 } else if (mask.icy(i, j) and mask.next_to_ice_free_ocean(i, j)) {
119 delta_p = stressbalance::margin_pressure_difference(mask.ocean(i, j), dry_mode,
120 H(i, j), bed(i, j), sea_level(i, j),
121 rho_ice, rho_ocean, g);
122 } else {
123 delta_p = m_fill_value;
124 }
125
126 (*result)(i, j) = delta_p;
127 }
128 } catch (...) {
129 loop.failed();
130 }
131 loop.check();
132
133
134 return result;
135 }
136
137 /*! @brief Report average basal mass balance flux over the reporting interval (grounded or floating
138 areas) */
139 class BMBSplit : public DiagAverageRate<IceModel>
140 {
141 public:
142 BMBSplit(const IceModel *m, AreaType flag)
143 : DiagAverageRate<IceModel>(m,
144 flag == GROUNDED
145 ? "basal_mass_flux_grounded"
146 : "basal_mass_flux_floating",
147 TOTAL_CHANGE), m_kind(flag) {
148 assert(flag != BOTH);
149
150 auto ismip6 = m_config->get_flag("output.ISMIP6");
151
152 std::string name, description, standard_name;
153 if (m_kind == GROUNDED) {
154 name = ismip6 ? "libmassbfgr" : "basal_mass_flux_grounded";
155 description = "average basal mass flux over the reporting interval (grounded areas)";
156 standard_name = ismip6 ? "land_ice_basal_specific_mass_balance_flux" : "";
157 } else {
158 name = ismip6 ? "libmassbffl" : "basal_mass_flux_floating";
159 description = "average basal mass flux over the reporting interval (floating areas)";
160 standard_name = ismip6 ? "land_ice_basal_specific_mass_balance_flux" : "";
161 }
162
163 m_vars = {SpatialVariableMetadata(m_sys, name)};
164 m_accumulator.metadata().set_string("units", "kg m-2");
165
166 set_attrs(description, standard_name, "kg m-2 s-1", "kg m-2 year-1", 0);
167 m_vars[0].set_string("cell_methods", "time: mean");
168
169 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
170 m_vars[0].set_string("comment", "positive flux corresponds to ice gain");
171 }
172
173 protected:
174 AreaType m_kind;
175 void update_impl(double dt) {
176 const IceModelVec2S &input = model->geometry_evolution().bottom_surface_mass_balance();
177 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
178
179 double ice_density = m_config->get_number("constants.ice.density");
180
181 // the accumulator has the units of kg/m^2, computed as
182 //
183 // accumulator += BMB (m) * ice_density (kg / m^3)
184
185 IceModelVec::AccessList list{&input, &cell_type, &m_accumulator};
186
187 for (Points p(*m_grid); p; p.next()) {
188 const int i = p.i(), j = p.j();
189
190 if (m_kind == GROUNDED and cell_type.grounded(i, j)) {
191 m_accumulator(i, j) += input(i, j) * ice_density;
192 } else if (m_kind == SHELF and cell_type.ocean(i, j)) {
193 m_accumulator(i, j) += input(i, j) * ice_density;
194 } else {
195 m_accumulator(i, j) = 0.0;
196 }
197 }
198
199 m_interval_length += dt;
200 }
201 };
202
203 //! \brief Computes vertically-averaged ice hardness.
204 class HardnessAverage : public Diag<IceModel>
205 {
206 public:
207 HardnessAverage(const IceModel *m);
208 protected:
209 virtual IceModelVec::Ptr compute_impl() const;
210 };
211
212 HardnessAverage::HardnessAverage(const IceModel *m)
213 : Diag<IceModel>(m) {
214
215 // set metadata:
216 m_vars = {SpatialVariableMetadata(m_sys, "hardav")};
217
218 // choice to use SSA power; see #285
219 const double power = 1.0 / m_config->get_number("stress_balance.ssa.Glen_exponent");
220 auto unitstr = pism::printf("Pa s%f", power);
221
222 set_attrs("vertical average of ice hardness", "",
223 unitstr, unitstr, 0);
224
225 m_vars[0].set_number("valid_min", 0);
226 m_vars[0].set_number("_FillValue", m_fill_value);
227 }
228
229 //! \brief Computes vertically-averaged ice hardness.
230 IceModelVec::Ptr HardnessAverage::compute_impl() const {
231
232 const rheology::FlowLaw *flow_law = model->stress_balance()->shallow()->flow_law().get();
233 if (flow_law == NULL) {
234 flow_law = model->stress_balance()->modifier()->flow_law().get();
235 if (flow_law == NULL) {
236 throw RuntimeError(PISM_ERROR_LOCATION,
237 "Can't compute vertically-averaged hardness: no flow law is used.");
238 }
239 }
240
241 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "hardav", WITHOUT_GHOSTS));
242 result->metadata() = m_vars[0];
243
244 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
245
246 const IceModelVec3& ice_enthalpy = model->energy_balance_model()->enthalpy();
247 const IceModelVec2S& ice_thickness = model->geometry().ice_thickness;
248
249 IceModelVec::AccessList list{&cell_type, &ice_enthalpy, &ice_thickness, result.get()};
250 ParallelSection loop(m_grid->com);
251 try {
252 for (Points p(*m_grid); p; p.next()) {
253 const int i = p.i(), j = p.j();
254
255 const double *Eij = ice_enthalpy.get_column(i,j);
256 const double H = ice_thickness(i,j);
257 if (cell_type.icy(i, j)) {
258 (*result)(i,j) = rheology::averaged_hardness(*flow_law,
259 H, m_grid->kBelowHeight(H),
260 &(m_grid->z()[0]), Eij);
261 } else { // put negative value below valid range
262 (*result)(i,j) = m_fill_value;
263 }
264 }
265 } catch (...) {
266 loop.failed();
267 }
268 loop.check();
269
270 return result;
271 }
272
273 //! \brief Computes a diagnostic field filled with processor rank values.
274 class Rank : public Diag<IceModel>
275 {
276 public:
277 Rank(const IceModel *m);
278 protected:
279 virtual IceModelVec::Ptr compute_impl() const;
280 };
281
282 Rank::Rank(const IceModel *m)
283 : Diag<IceModel>(m) {
284 m_vars = {SpatialVariableMetadata(m_sys, "rank")};
285 set_attrs("processor rank", "", "1", "", 0);
286 m_vars[0].set_time_independent(true);
287 m_vars[0].set_output_type(PISM_INT);
288 }
289
290 IceModelVec::Ptr Rank::compute_impl() const {
291
292 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "rank", WITHOUT_GHOSTS));
293 result->metadata() = m_vars[0];
294
295 IceModelVec::AccessList list{result.get()};
296
297 for (Points p(*m_grid); p; p.next()) {
298 (*result)(p.i(),p.j()) = m_grid->rank();
299 }
300
301 return result;
302 }
303
304 //! \brief Computes CTS, CTS = E/E_s(p).
305 class CTS : public Diag<IceModel>
306 {
307 public:
308 CTS(const IceModel *m);
309 protected:
310 virtual IceModelVec::Ptr compute_impl() const;
311 };
312
313 CTS::CTS(const IceModel *m)
314 : Diag<IceModel>(m) {
315
316 // set metadata:
317 m_vars = {SpatialVariableMetadata(m_sys, "cts", m_grid->z())};
318
319 set_attrs("cts = E/E_s(p), so cold-temperate transition surface is at cts = 1", "",
320 "1", "1", 0);
321 }
322
323 IceModelVec::Ptr CTS::compute_impl() const {
324
325 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "cts", WITHOUT_GHOSTS));
326 result->metadata() = m_vars[0];
327
328 energy::compute_cts(model->energy_balance_model()->enthalpy(),
329 model->geometry().ice_thickness, *result);
330
331 return result;
332 }
333
334 //! \brief Computes ice temperature from enthalpy.
335 class Temperature : public Diag<IceModel>
336 {
337 public:
338 Temperature(const IceModel *m);
339 protected:
340 virtual IceModelVec::Ptr compute_impl() const;
341 };
342
343 Temperature::Temperature(const IceModel *m)
344 : Diag<IceModel>(m) {
345
346 // set metadata:
347 m_vars = {SpatialVariableMetadata(m_sys, "temp", m_grid->z())};
348
349 set_attrs("ice temperature", "land_ice_temperature", "K", "K", 0);
350 m_vars[0].set_number("valid_min", 0);
351 }
352
353 IceModelVec::Ptr Temperature::compute_impl() const {
354
355 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "temp", WITHOUT_GHOSTS));
356 result->metadata() = m_vars[0];
357
358 const IceModelVec2S &thickness = model->geometry().ice_thickness;
359 const IceModelVec3 &enthalpy = model->energy_balance_model()->enthalpy();
360
361 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
362
363 double *Tij;
364 const double *Enthij; // columns of these values
365
366 IceModelVec::AccessList list{result.get(), &enthalpy, &thickness};
367
368 ParallelSection loop(m_grid->com);
369 try {
370 for (Points p(*m_grid); p; p.next()) {
371 const int i = p.i(), j = p.j();
372
373 Tij = result->get_column(i,j);
374 Enthij = enthalpy.get_column(i,j);
375 for (unsigned int k=0; k <m_grid->Mz(); ++k) {
376 const double depth = thickness(i,j) - m_grid->z(k);
377 Tij[k] = EC->temperature(Enthij[k], EC->pressure(depth));
378 }
379 }
380 } catch (...) {
381 loop.failed();
382 }
383 loop.check();
384
385 return result;
386 }
387
388 //! \brief Compute the pressure-adjusted temperature in degrees C corresponding
389 //! to ice temperature.
390 class TemperaturePA : public Diag<IceModel>
391 {
392 public:
393 TemperaturePA(const IceModel *m);
394 protected:
395 virtual IceModelVec::Ptr compute_impl() const;
396 };
397
398
399 TemperaturePA::TemperaturePA(const IceModel *m)
400 : Diag<IceModel>(m) {
401
402 // set metadata:
403 m_vars = {SpatialVariableMetadata(m_sys, "temp_pa", m_grid->z())};
404
405 set_attrs("pressure-adjusted ice temperature (degrees above pressure-melting point)", "",
406 "deg_C", "deg_C", 0);
407 m_vars[0].set_number("valid_max", 0);
408 }
409
410 IceModelVec::Ptr TemperaturePA::compute_impl() const {
411 bool cold_mode = m_config->get_flag("energy.temperature_based");
412 double melting_point_temp = m_config->get_number("constants.fresh_water.melting_point_temperature");
413
414 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "temp_pa", WITHOUT_GHOSTS));
415 result->metadata() = m_vars[0];
416
417 const IceModelVec2S &thickness = model->geometry().ice_thickness;
418 const IceModelVec3 &enthalpy = model->energy_balance_model()->enthalpy();
419
420 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
421
422 double *Tij;
423 const double *Enthij; // columns of these values
424
425 IceModelVec::AccessList list{result.get(), &enthalpy, &thickness};
426
427 ParallelSection loop(m_grid->com);
428 try {
429 for (Points pt(*m_grid); pt; pt.next()) {
430 const int i = pt.i(), j = pt.j();
431
432 Tij = result->get_column(i,j);
433 Enthij = enthalpy.get_column(i,j);
434 for (unsigned int k=0; k < m_grid->Mz(); ++k) {
435 const double depth = thickness(i,j) - m_grid->z(k),
436 p = EC->pressure(depth);
437 Tij[k] = EC->pressure_adjusted_temperature(Enthij[k], p);
438
439 if (cold_mode) { // if ice is temperate then its pressure-adjusted temp
440 // is 273.15
441 if (EC->is_temperate_relaxed(Enthij[k],p) && (thickness(i,j) > 0)) {
442 Tij[k] = melting_point_temp;
443 }
444 }
445
446 }
447 }
448 } catch (...) {
449 loop.failed();
450 }
451 loop.check();
452
453 result->shift(-melting_point_temp);
454
455 return result;
456 }
457
458 //! \brief Computes basal values of the pressure-adjusted temperature.
459 class TemperaturePABasal : public Diag<IceModel>
460 {
461 public:
462 TemperaturePABasal(const IceModel *m);
463 protected:
464 virtual IceModelVec::Ptr compute_impl() const;
465 };
466
467 TemperaturePABasal::TemperaturePABasal(const IceModel *m)
468 : Diag<IceModel>(m) {
469
470 // set metadata:
471 m_vars = {SpatialVariableMetadata(m_sys, "temppabase")};
472
473 set_attrs("pressure-adjusted ice temperature at the base of ice", "",
474 "Celsius", "Celsius", 0);
475 }
476
477 IceModelVec::Ptr TemperaturePABasal::compute_impl() const {
478
479 bool cold_mode = m_config->get_flag("energy.temperature_based");
480 double melting_point_temp = m_config->get_number("constants.fresh_water.melting_point_temperature");
481
482 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "temp_pa_base", WITHOUT_GHOSTS));
483 result->metadata() = m_vars[0];
484
485 const IceModelVec2S &thickness = model->geometry().ice_thickness;
486 const IceModelVec3 &enthalpy = model->energy_balance_model()->enthalpy();
487
488 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
489
490 const double *Enthij; // columns of these values
491
492 IceModelVec::AccessList list{result.get(), &enthalpy, &thickness};
493
494 ParallelSection loop(m_grid->com);
495 try {
496 for (Points pt(*m_grid); pt; pt.next()) {
497 const int i = pt.i(), j = pt.j();
498
499 Enthij = enthalpy.get_column(i,j);
500
501 const double depth = thickness(i,j),
502 p = EC->pressure(depth);
503 (*result)(i,j) = EC->pressure_adjusted_temperature(Enthij[0], p);
504
505 if (cold_mode) { // if ice is temperate then its pressure-adjusted temp
506 // is 273.15
507 if (EC->is_temperate_relaxed(Enthij[0],p) && (thickness(i,j) > 0)) {
508 (*result)(i,j) = melting_point_temp;
509 }
510 }
511 }
512 } catch (...) {
513 loop.failed();
514 }
515 loop.check();
516
517 result->shift(-melting_point_temp);
518
519 return result;
520 }
521
522 //! \brief Computes surface values of ice enthalpy.
523 class IceEnthalpySurface : public Diag<IceModel>
524 {
525 public:
526 IceEnthalpySurface(const IceModel *m);
527 protected:
528 virtual IceModelVec::Ptr compute_impl() const;
529 };
530
531 IceEnthalpySurface::IceEnthalpySurface(const IceModel *m)
532 : Diag<IceModel>(m) {
533
534 // set metadata:
535 m_vars = {SpatialVariableMetadata(m_sys, "enthalpysurf")};
536
537 set_attrs("ice enthalpy at 1m below the ice surface", "",
538 "J kg-1", "J kg-1", 0);
539 m_vars[0].set_number("_FillValue", m_fill_value);
540 }
541
542 IceModelVec::Ptr IceEnthalpySurface::compute_impl() const {
543
544 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "enthalpysurf", WITHOUT_GHOSTS));
545 result->metadata() = m_vars[0];
546
547 // compute levels corresponding to 1 m below the ice surface:
548
549 const IceModelVec3& ice_enthalpy = model->energy_balance_model()->enthalpy();
550 const IceModelVec2S& ice_thickness = model->geometry().ice_thickness;
551
552 IceModelVec::AccessList list{&ice_thickness, result.get()};
553
554 for (Points p(*m_grid); p; p.next()) {
555 const int i = p.i(), j = p.j();
556
557 (*result)(i,j) = std::max(ice_thickness(i,j) - 1.0, 0.0);
558 }
559
560 ice_enthalpy.getSurfaceValues(*result, *result); // z=0 slice
561
562 for (Points p(*m_grid); p; p.next()) {
563 const int i = p.i(), j = p.j();
564
565 if (ice_thickness(i,j) <= 1.0) {
566 (*result)(i,j) = m_fill_value;
567 }
568 }
569
570 return result;
571 }
572
573 //! \brief Computes enthalpy at the base of the ice.
574 class IceEnthalpyBasal : public Diag<IceModel>
575 {
576 public:
577 IceEnthalpyBasal(const IceModel *m);
578 protected:
579 virtual IceModelVec::Ptr compute_impl() const;
580 };
581
582 IceEnthalpyBasal::IceEnthalpyBasal(const IceModel *m)
583 : Diag<IceModel>(m) {
584
585 // set metadata:
586 m_vars = {SpatialVariableMetadata(m_sys, "enthalpybase")};
587
588 set_attrs("ice enthalpy at the base of ice", "",
589 "J kg-1", "J kg-1", 0);
590 m_vars[0].set_number("_FillValue", m_fill_value);
591 }
592
593 IceModelVec::Ptr IceEnthalpyBasal::compute_impl() const {
594
595 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "enthalpybase", WITHOUT_GHOSTS));
596 result->metadata() = m_vars[0];
597
598 model->energy_balance_model()->enthalpy().getHorSlice(*result, 0.0); // z=0 slice
599
600 result->mask_by(model->geometry().ice_thickness, m_fill_value);
601
602 return result;
603 }
604
605 //! \brief Computes ice temperature at the base of the ice.
606 class TemperatureBasal : public Diag<IceModel>
607 {
608 public:
609 TemperatureBasal(const IceModel *m, AreaType flag);
610 private:
611 IceModelVec::Ptr compute_impl() const;
612
613 AreaType m_area_type;
614 };
615
616 TemperatureBasal::TemperatureBasal(const IceModel *m, AreaType area_type)
617 : Diag<IceModel>(m), m_area_type(area_type) {
618
619 std::string name, long_name, standard_name;
620 switch (area_type) {
621 case GROUNDED:
622 name = "litempbotgr";
623 long_name = "ice temperature at the bottom surface of grounded ice";
624 standard_name = "temperature_at_base_of_ice_sheet_model";
625 break;
626 case SHELF:
627 name = "litempbotfl";
628 long_name = "ice temperature at the bottom surface of floating ice";
629 standard_name = "temperature_at_base_of_ice_sheet_model";
630 break;
631 case BOTH:
632 name = "tempbase";
633 long_name = "ice temperature at the base of ice";
634 standard_name = "land_ice_basal_temperature";
635 break;
636 }
637 // set metadata:
638 m_vars = {SpatialVariableMetadata(m_sys, name)};
639
640 set_attrs(long_name, standard_name, "K", "K", 0);
641 m_vars[0].set_number("_FillValue", m_fill_value);
642 }
643
644 IceModelVec::Ptr TemperatureBasal::compute_impl() const {
645
646 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "basal_temperature", WITHOUT_GHOSTS));
647 result->metadata(0) = m_vars[0];
648
649 const IceModelVec2S &thickness = model->geometry().ice_thickness;
650
651 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
652
653 model->energy_balance_model()->enthalpy().getHorSlice(*result, 0.0); // z=0 (basal) slice
654 // Now result contains basal enthalpy.
655
656 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
657
658 IceModelVec::AccessList list{&cell_type, result.get(), &thickness};
659
660 ParallelSection loop(m_grid->com);
661 try {
662 for (Points p(*m_grid); p; p.next()) {
663 const int i = p.i(), j = p.j();
664
665 double
666 depth = thickness(i, j),
667 pressure = EC->pressure(depth),
668 T = EC->temperature((*result)(i, j), pressure);
669
670 if ((m_area_type == BOTH and cell_type.icy(i, j)) or
671 (m_area_type == GROUNDED and cell_type.grounded_ice(i, j)) or
672 (m_area_type == SHELF and cell_type.floating_ice(i, j))) {
673 (*result)(i, j) = T;
674 } else {
675 (*result)(i,j) = m_fill_value;
676 }
677 }
678 } catch (...) {
679 loop.failed();
680 }
681 loop.check();
682
683 return result;
684 }
685
686 //! \brief Computes ice temperature at the surface of the ice.
687 class TemperatureSurface : public Diag<IceModel>
688 {
689 public:
690 TemperatureSurface(const IceModel *m);
691 protected:
692 virtual IceModelVec::Ptr compute_impl() const;
693 };
694
695 TemperatureSurface::TemperatureSurface(const IceModel *m)
696 : Diag<IceModel>(m) {
697
698 // set metadata:
699 m_vars = {SpatialVariableMetadata(m_sys, "tempsurf")};
700
701 set_attrs("ice temperature at 1m below the ice surface",
702 "temperature_at_ground_level_in_snow_or_firn", // InitMIP "standard" name
703 "K", "K", 0);
704 m_vars[0].set_number("_FillValue", m_fill_value);
705 }
706
707 IceModelVec::Ptr TemperatureSurface::compute_impl() const {
708
709 const IceModelVec2S &thickness = model->geometry().ice_thickness;
710
711 IceModelVec::Ptr enth = IceEnthalpySurface(model).compute();
712 IceModelVec2S::Ptr result = IceModelVec2S::To2DScalar(enth);
713
714 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
715
716 // result contains surface enthalpy; note that it is allocated by
717 // IceEnthalpySurface::compute().
718
719 IceModelVec::AccessList list{result.get(), &thickness};
720
721 double depth = 1.0,
722 pressure = EC->pressure(depth);
723 ParallelSection loop(m_grid->com);
724 try {
725 for (Points p(*m_grid); p; p.next()) {
726 const int i = p.i(), j = p.j();
727
728 if (thickness(i,j) > 1) {
729 (*result)(i,j) = EC->temperature((*result)(i,j), pressure);
730 } else {
731 (*result)(i,j) = m_fill_value;
732 }
733 }
734 } catch (...) {
735 loop.failed();
736 }
737 loop.check();
738
739 result->metadata(0) = m_vars[0];
740 return result;
741 }
742
743 //! \brief Computes the liquid water fraction.
744 class LiquidFraction : public Diag<IceModel>
745 {
746 public:
747 LiquidFraction(const IceModel *m);
748 protected:
749 virtual IceModelVec::Ptr compute_impl() const;
750 };
751
752 LiquidFraction::LiquidFraction(const IceModel *m)
753 : Diag<IceModel>(m) {
754
755 // set metadata:
756 m_vars = {SpatialVariableMetadata(m_sys, "liqfrac", m_grid->z())};
757
758 set_attrs("liquid water fraction in ice (between 0 and 1)", "",
759 "1", "1", 0);
760 m_vars[0].set_numbers("valid_range", {0.0, 1.0});
761 }
762
763 IceModelVec::Ptr LiquidFraction::compute_impl() const {
764
765 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "liqfrac", WITHOUT_GHOSTS));
766 result->metadata(0) = m_vars[0];
767
768 bool cold_mode = m_config->get_flag("energy.temperature_based");
769
770 if (cold_mode) {
771 result->set(0.0);
772 } else {
773 energy::compute_liquid_water_fraction(model->energy_balance_model()->enthalpy(),
774 model->geometry().ice_thickness,
775 *result);
776 }
777
778 return result;
779 }
780
781 //! \brief Computes the total thickness of temperate ice in a column.
782 class TemperateIceThickness : public Diag<IceModel>
783 {
784 public:
785 TemperateIceThickness(const IceModel *m);
786 protected:
787 virtual IceModelVec::Ptr compute_impl() const;
788 };
789
790 TemperateIceThickness::TemperateIceThickness(const IceModel *m)
791 : Diag<IceModel>(m) {
792
793 // set metadata:
794 m_vars = {SpatialVariableMetadata(m_sys,
795 "tempicethk")};
796
797 set_attrs("temperate ice thickness (total column content)", "",
798 "m", "m", 0);
799 m_vars[0].set_number("_FillValue", m_fill_value);
800 }
801
802 IceModelVec::Ptr TemperateIceThickness::compute_impl() const {
803
804 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "tempicethk", WITHOUT_GHOSTS));
805 result->metadata(0) = m_vars[0];
806
807 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
808 const IceModelVec3& ice_enthalpy = model->energy_balance_model()->enthalpy();
809 const IceModelVec2S& ice_thickness = model->geometry().ice_thickness;
810
811 IceModelVec::AccessList list{&cell_type, result.get(), &ice_enthalpy, &ice_thickness};
812
813 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
814
815 ParallelSection loop(m_grid->com);
816 try {
817 for (Points p(*m_grid); p; p.next()) {
818 const int i = p.i(), j = p.j();
819
820 if (cell_type.icy(i, j)) {
821 const double *Enth = ice_enthalpy.get_column(i,j);
822 double H_temperate = 0.0;
823 const double H = ice_thickness(i,j);
824 const unsigned int ks = m_grid->kBelowHeight(H);
825
826 for (unsigned int k=0; k<ks; ++k) { // FIXME issue #15
827 double pressure = EC->pressure(H - m_grid->z(k));
828
829 if (EC->is_temperate_relaxed(Enth[k], pressure)) {
830 H_temperate += m_grid->z(k+1) - m_grid->z(k);
831 }
832 }
833
834 double pressure = EC->pressure(H - m_grid->z(ks));
835 if (EC->is_temperate_relaxed(Enth[ks], pressure)) {
836 H_temperate += H - m_grid->z(ks);
837 }
838
839 (*result)(i,j) = H_temperate;
840 } else {
841 // ice-free
842 (*result)(i,j) = m_fill_value;
843 }
844 }
845 } catch (...) {
846 loop.failed();
847 }
848 loop.check();
849
850 return result;
851 }
852
853 //! \brief Computes the thickness of the basal layer of temperate ice.
854 class TemperateIceThicknessBasal : public Diag<IceModel>
855 {
856 public:
857 TemperateIceThicknessBasal(const IceModel *m);
858 protected:
859 virtual IceModelVec::Ptr compute_impl() const;
860 };
861
862 TemperateIceThicknessBasal::TemperateIceThicknessBasal(const IceModel *m)
863 : Diag<IceModel>(m) {
864
865 // set metadata:
866 m_vars = {SpatialVariableMetadata(m_sys,
867 "tempicethk_basal")};
868
869 set_attrs("thickness of the basal layer of temperate ice", "",
870 "m", "m", 0);
871 m_vars[0].set_number("_FillValue", m_fill_value);
872 }
873
874 /*!
875 * Uses linear interpolation to go beyond vertical grid resolution.
876 */
877 IceModelVec::Ptr TemperateIceThicknessBasal::compute_impl() const {
878
879 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "tempicethk_basal", WITHOUT_GHOSTS));
880 result->metadata(0) = m_vars[0];
881
882 EnthalpyConverter::Ptr EC = model->ctx()->enthalpy_converter();
883
884 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
885 const IceModelVec3& ice_enthalpy = model->energy_balance_model()->enthalpy();
886 const IceModelVec2S& ice_thickness = model->geometry().ice_thickness;
887
888 IceModelVec::AccessList list{&cell_type, result.get(), &ice_thickness, &ice_enthalpy};
889
890 ParallelSection loop(m_grid->com);
891 try {
892 for (Points p(*m_grid); p; p.next()) {
893 const int i = p.i(), j = p.j();
894
895 double H = ice_thickness(i,j);
896
897 // if we have no ice, go on to the next grid point (this cell will be
898 // marked as "missing" later)
899 if (cell_type.ice_free(i, j)) {
900 (*result)(i,j) = m_fill_value;
901 continue;
902 }
903
904 const double *Enth = ice_enthalpy.get_column(i,j);
905
906 unsigned int ks = m_grid->kBelowHeight(H);
907
908 unsigned int k = 0;
909 double pressure = EC->pressure(H - m_grid->z(k));
910 while (k <= ks) { // FIXME issue #15
911 pressure = EC->pressure(H - m_grid->z(k));
912
913 if (EC->is_temperate_relaxed(Enth[k],pressure)) {
914 k++;
915 } else {
916 break;
917 }
918 }
919 // after this loop 'pressure' is equal to the pressure at the first level
920 // that is cold
921
922 // no temperate ice at all; go to the next grid point
923 if (k == 0) {
924 (*result)(i,j) = 0.0;
925 continue;
926 }
927
928 // the whole column is temperate (except, possibly, some ice between
929 // z(ks) and the total thickness; we ignore it)
930 if (k == ks + 1) {
931 (*result)(i,j) = m_grid->z(ks);
932 continue;
933 }
934
935 double
936 pressure_0 = EC->pressure(H - m_grid->z(k-1)),
937 dz = m_grid->z(k) - m_grid->z(k-1),
938 slope1 = (Enth[k] - Enth[k-1]) / dz,
939 slope2 = (EC->enthalpy_cts(pressure) - EC->enthalpy_cts(pressure_0)) / dz;
940
941 if (slope1 != slope2) {
942 (*result)(i,j) = m_grid->z(k-1) +
943 (EC->enthalpy_cts(pressure_0) - Enth[k-1]) / (slope1 - slope2);
944
945 // check if the resulting thickness is valid:
946 (*result)(i,j) = std::max((*result)(i,j), m_grid->z(k-1));
947 (*result)(i,j) = std::min((*result)(i,j), m_grid->z(k));
948 } else {
949 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Linear interpolation of the thickness of"
950 " the basal temperate layer failed:\n"
951 "(i=%d, j=%d, k=%d, ks=%d)\n",
952 i, j, k, ks);
953 }
954 }
955 } catch (...) {
956 loop.failed();
957 }
958 loop.check();
959
960
961 return result;
962 }
963
964 namespace scalar {
965
966 //! \brief Computes the total ice volume in glacierized areas.
967 class IceVolumeGlacierized : public TSDiag<TSSnapshotDiagnostic, IceModel>
968 {
969 public:
970 IceVolumeGlacierized(IceModel *m)
971 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_glacierized") {
972
973 set_units("m3", "m3");
974 m_ts.variable().set_string("long_name", "volume of the ice in glacierized areas");
975 m_ts.variable().set_number("valid_min", 0.0);
976 }
977 double compute() {
978 return ice_volume(model->geometry(),
979 m_config->get_number("output.ice_free_thickness_standard"));
980 }
981 };
982
983 //! \brief Computes the total ice volume.
984 class IceVolume : public TSDiag<TSSnapshotDiagnostic, IceModel>
985 {
986 public:
987 IceVolume(IceModel *m)
988 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume") {
989
990 set_units("m3", "m3");
991 m_ts.variable().set_string("long_name", "volume of the ice, including seasonal cover");
992 m_ts.variable().set_number("valid_min", 0.0);
993 }
994
995 double compute() {
996 return ice_volume(model->geometry(), 0.0);
997 }
998 };
999
1000 //! \brief Computes the total ice volume which is relevant for sea-level
1001 class SeaLevelRisePotential : public TSDiag<TSSnapshotDiagnostic, IceModel>
1002 {
1003 public:
1004 SeaLevelRisePotential(const IceModel *m)
1005 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "sea_level_rise_potential") {
1006
1007 set_units("m", "m");
1008 m_ts.variable().set_string("long_name", "the sea level rise that would result if all the ice were melted");
1009 m_ts.variable().set_number("valid_min", 0.0);
1010 }
1011
1012 double compute() {
1013 return sea_level_rise_potential(model->geometry(),
1014 m_config->get_number("output.ice_free_thickness_standard"));
1015 }
1016 };
1017
1018 //! \brief Computes the rate of change of the total ice volume in glacierized areas.
1019 class IceVolumeRateOfChangeGlacierized : public TSDiag<TSRateDiagnostic, IceModel>
1020 {
1021 public:
1022 IceVolumeRateOfChangeGlacierized(IceModel *m)
1023 : TSDiag<TSRateDiagnostic, IceModel>(m, "tendency_of_ice_volume_glacierized") {
1024
1025 set_units("m3 s-1", "m3 year-1");
1026 m_ts.variable().set_string("long_name", "rate of change of the ice volume in glacierized areas");
1027 }
1028
1029 double compute() {
1030 return ice_volume(model->geometry(),
1031 m_config->get_number("output.ice_free_thickness_standard"));
1032 }
1033 };
1034
1035 //! \brief Computes the rate of change of the total ice volume.
1036 class IceVolumeRateOfChange : public TSDiag<TSRateDiagnostic, IceModel>
1037 {
1038 public:
1039 IceVolumeRateOfChange(IceModel *m)
1040 : TSDiag<TSRateDiagnostic, IceModel>(m, "tendency_of_ice_volume") {
1041
1042 set_units("m3 s-1", "m3 year-1");
1043 m_ts.variable().set_string("long_name",
1044 "rate of change of the ice volume, including seasonal cover");
1045 }
1046
1047 double compute() {
1048 return ice_volume(model->geometry(), 0.0);
1049 }
1050 };
1051
1052 //! \brief Computes the total ice area.
1053 class IceAreaGlacierized : public TSDiag<TSSnapshotDiagnostic, IceModel>
1054 {
1055 public:
1056 IceAreaGlacierized(IceModel *m)
1057 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_area_glacierized") {
1058
1059 set_units("m2", "m2");
1060 m_ts.variable().set_string("long_name", "glacierized area");
1061 m_ts.variable().set_number("valid_min", 0.0);
1062 }
1063
1064 double compute() {
1065 return ice_area(model->geometry(), m_config->get_number("output.ice_free_thickness_standard"));
1066 }
1067 };
1068
1069 //! \brief Computes the total mass of the ice not displacing sea water.
1070 class IceMassNotDisplacingSeaWater : public TSDiag<TSSnapshotDiagnostic, IceModel>
1071 {
1072 public:
1073 IceMassNotDisplacingSeaWater(const IceModel *m)
1074 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "limnsw") {
1075
1076 set_units("kg", "kg");
1077 m_ts.variable().set_string("long_name", "mass of the ice not displacing sea water");
1078 m_ts.variable().set_string("standard_name", "land_ice_mass_not_displacing_sea_water");
1079 m_ts.variable().set_number("valid_min", 0.0);
1080 }
1081
1082 double compute() {
1083
1084 const double
1085 thickness_standard = m_config->get_number("output.ice_free_thickness_standard"),
1086 ice_density = m_config->get_number("constants.ice.density"),
1087 ice_volume = ice_volume_not_displacing_seawater(model->geometry(),
1088 thickness_standard),
1089 ice_mass = ice_volume * ice_density;
1090
1091 return ice_mass;
1092 }
1093 };
1094
1095 //! \brief Computes the total ice mass in glacierized areas.
1096 class IceMassGlacierized : public TSDiag<TSSnapshotDiagnostic, IceModel>
1097 {
1098 public:
1099 IceMassGlacierized(IceModel *m)
1100 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_mass_glacierized") {
1101
1102 set_units("kg", "kg");
1103 m_ts.variable().set_string("long_name", "mass of the ice in glacierized areas");
1104 m_ts.variable().set_number("valid_min", 0.0);
1105 }
1106
1107 double compute() {
1108 double
1109 ice_density = m_config->get_number("constants.ice.density"),
1110 thickness_standard = m_config->get_number("output.ice_free_thickness_standard");
1111 return ice_volume(model->geometry(), thickness_standard) * ice_density;
1112 }
1113 };
1114
1115 //! \brief Computes the total ice mass.
1116 class IceMass : public TSDiag<TSSnapshotDiagnostic, IceModel>
1117 {
1118 public:
1119 IceMass(IceModel *m)
1120 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_mass") {
1121
1122 if (m_config->get_flag("output.ISMIP6")) {
1123 m_ts.variable().set_name("lim");
1124 }
1125
1126 set_units("kg", "kg");
1127 m_ts.variable().set_string("long_name", "mass of the ice, including seasonal cover");
1128 m_ts.variable().set_string("standard_name", "land_ice_mass");
1129 m_ts.variable().set_number("valid_min", 0.0);
1130 }
1131
1132 double compute() {
1133 return (ice_volume(model->geometry(), 0.0) *
1134 m_config->get_number("constants.ice.density"));
1135 }
1136 };
1137
1138 //! \brief Computes the rate of change of the total ice mass in glacierized areas.
1139 class IceMassRateOfChangeGlacierized : public TSDiag<TSRateDiagnostic, IceModel>
1140 {
1141 public:
1142 IceMassRateOfChangeGlacierized(IceModel *m)
1143 : TSDiag<TSRateDiagnostic, IceModel>(m, "tendency_of_ice_mass_glacierized") {
1144
1145 set_units("kg s-1", "Gt year-1");
1146 m_ts.variable().set_string("long_name", "rate of change of the ice mass in glacierized areas");
1147 }
1148
1149 double compute() {
1150 double
1151 ice_density = m_config->get_number("constants.ice.density"),
1152 thickness_threshold = m_config->get_number("output.ice_free_thickness_standard");
1153 return ice_volume(model->geometry(), thickness_threshold) * ice_density;
1154 }
1155 };
1156
1157 //! \brief Computes the rate of change of the total ice mass due to flow (influx due to
1158 //! prescribed constant-in-time ice thickness).
1159 /*!
1160 * This is the change in mass resulting from prescribing (fixing) ice thickness.
1161 */
1162 class IceMassRateOfChangeDueToFlow : public TSDiag<TSFluxDiagnostic, IceModel>
1163 {
1164 public:
1165 IceMassRateOfChangeDueToFlow(IceModel *m)
1166 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_flow") {
1167
1168 set_units("kg s-1", "Gt year-1");
1169 m_ts.variable().set_string("long_name", "rate of change of the mass of ice due to flow"
1170 " (i.e. prescribed ice thickness)");
1171 }
1172
1173 double compute() {
1174
1175 const double
1176 ice_density = m_config->get_number("constants.ice.density");
1177
1178 const IceModelVec2S
1179 &dH = model->geometry_evolution().thickness_change_due_to_flow(),
1180 &dV = model->geometry_evolution().area_specific_volume_change_due_to_flow();
1181
1182 auto cell_area = m_grid->cell_area();
1183
1184 IceModelVec::AccessList list{&dH, &dV};
1185
1186 double volume_change = 0.0;
1187 for (Points p(*m_grid); p; p.next()) {
1188 const int i = p.i(), j = p.j();
1189 // m * m^2 = m^3
1190 volume_change += (dH(i, j) + dV(i, j)) * cell_area;
1191 }
1192
1193 // (kg/m^3) * m^3 = kg
1194 return ice_density * GlobalSum(m_grid->com, volume_change);
1195 }
1196 };
1197
1198 //! \brief Computes the rate of change of the total ice mass.
1199 class IceMassRateOfChange : public TSDiag<TSRateDiagnostic, IceModel>
1200 {
1201 public:
1202 IceMassRateOfChange(IceModel *m)
1203 : TSDiag<TSRateDiagnostic, IceModel>(m, "tendency_of_ice_mass") {
1204
1205 set_units("kg s-1", "Gt year-1");
1206 m_ts.variable().set_string("long_name",
1207 "rate of change of the mass of ice, including seasonal cover");
1208 }
1209
1210 double compute() {
1211 const double ice_density = m_config->get_number("constants.ice.density");
1212 return ice_volume(model->geometry(), 0.0) * ice_density;
1213 }
1214 };
1215
1216
1217 //! \brief Computes the total volume of the temperate ice in glacierized areas.
1218 class IceVolumeGlacierizedTemperate : public TSDiag<TSSnapshotDiagnostic, IceModel>
1219 {
1220 public:
1221 IceVolumeGlacierizedTemperate(IceModel *m)
1222 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_glacierized_temperate") {
1223
1224 set_units("m3", "m3");
1225 m_ts.variable().set_string("long_name", "volume of temperate ice in glacierized areas");
1226 m_ts.variable().set_number("valid_min", 0.0);
1227 }
1228
1229 double compute() {
1230 return model->ice_volume_temperate(m_config->get_number("output.ice_free_thickness_standard"));
1231 }
1232 };
1233
1234 //! \brief Computes the total volume of the temperate ice.
1235 class IceVolumeTemperate : public TSDiag<TSSnapshotDiagnostic, IceModel>
1236 {
1237 public:
1238 IceVolumeTemperate(IceModel *m)
1239 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_temperate") {
1240
1241 set_units("m3", "m3");
1242 m_ts.variable().set_string("long_name", "volume of temperate ice, including seasonal cover");
1243 m_ts.variable().set_number("valid_min", 0.0);
1244 }
1245
1246 double compute() {
1247 return model->ice_volume_temperate(0.0);
1248 }
1249 };
1250
1251 //! \brief Computes the total volume of the cold ice in glacierized areas.
1252 class IceVolumeGlacierizedCold : public TSDiag<TSSnapshotDiagnostic, IceModel>
1253 {
1254 public:
1255 IceVolumeGlacierizedCold(IceModel *m)
1256 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_glacierized_cold") {
1257
1258 set_units("m3", "m3");
1259 m_ts.variable().set_string("long_name", "volume of cold ice in glacierized areas");
1260 m_ts.variable().set_number("valid_min", 0.0);
1261 }
1262
1263 double compute() {
1264 return model->ice_volume_cold(m_config->get_number("output.ice_free_thickness_standard"));
1265 }
1266 };
1267
1268 //! \brief Computes the total volume of the cold ice.
1269 class IceVolumeCold : public TSDiag<TSSnapshotDiagnostic, IceModel>
1270 {
1271 public:
1272 IceVolumeCold(IceModel *m)
1273 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_cold") {
1274
1275 set_units("m3", "m3");
1276 m_ts.variable().set_string("long_name", "volume of cold ice, including seasonal cover");
1277 m_ts.variable().set_number("valid_min", 0.0);
1278 }
1279
1280 double compute() {
1281 return model->ice_volume_cold(0.0);
1282 }
1283 };
1284
1285 //! \brief Computes the total area of the temperate ice.
1286 class IceAreaGlacierizedTemperateBase : public TSDiag<TSSnapshotDiagnostic, IceModel>
1287 {
1288 public:
1289 IceAreaGlacierizedTemperateBase(IceModel *m)
1290 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_area_glacierized_temperate_base") {
1291
1292 set_units("m2", "m2");
1293 m_ts.variable().set_string("long_name", "glacierized area where basal ice is temperate");
1294 m_ts.variable().set_number("valid_min", 0.0);
1295 }
1296
1297 double compute() {
1298 return model->ice_area_temperate(m_config->get_number("output.ice_free_thickness_standard"));
1299 }
1300 };
1301
1302 //! \brief Computes the total area of the cold ice.
1303 class IceAreaGlacierizedColdBase : public TSDiag<TSSnapshotDiagnostic, IceModel>
1304 {
1305 public:
1306 IceAreaGlacierizedColdBase(IceModel *m)
1307 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_area_glacierized_cold_base") {
1308
1309 set_units("m2", "m2");
1310 m_ts.variable().set_string("long_name", "glacierized area where basal ice is cold");
1311 m_ts.variable().set_number("valid_min", 0.0);
1312 }
1313
1314 double compute() {
1315 return model->ice_area_cold(m_config->get_number("output.ice_free_thickness_standard"));
1316 }
1317 };
1318
1319 //! \brief Computes the total ice enthalpy in glacierized areas.
1320 class IceEnthalpyGlacierized : public TSDiag<TSSnapshotDiagnostic, IceModel>
1321 {
1322 public:
1323 IceEnthalpyGlacierized(IceModel *m)
1324 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_enthalpy_glacierized") {
1325
1326 set_units("J", "J");
1327 m_ts.variable().set_string("long_name", "enthalpy of the ice in glacierized areas");
1328 m_ts.variable().set_number("valid_min", 0.0);
1329 }
1330
1331 double compute() {
1332 return energy::total_ice_enthalpy(m_config->get_number("output.ice_free_thickness_standard"),
1333 model->energy_balance_model()->enthalpy(),
1334 model->geometry().ice_thickness);
1335 }
1336 };
1337
1338 //! \brief Computes the total ice enthalpy.
1339 class IceEnthalpy : public TSDiag<TSSnapshotDiagnostic, IceModel>
1340 {
1341 public:
1342 IceEnthalpy(IceModel *m)
1343 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_enthalpy") {
1344
1345 set_units("J", "J");
1346 m_ts.variable().set_string("long_name", "enthalpy of the ice, including seasonal cover");
1347 m_ts.variable().set_number("valid_min", 0.0);
1348 }
1349
1350 double compute() {
1351 return energy::total_ice_enthalpy(0.0,
1352 model->energy_balance_model()->enthalpy(),
1353 model->geometry().ice_thickness);
1354 }
1355 };
1356
1357 //! \brief Computes the total grounded ice area.
1358 class IceAreaGlacierizedGrounded : public TSDiag<TSSnapshotDiagnostic, IceModel>
1359 {
1360 public:
1361 IceAreaGlacierizedGrounded(IceModel *m)
1362 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_area_glacierized_grounded") {
1363
1364 if (m_config->get_flag("output.ISMIP6")) {
1365 m_ts.variable().set_name("iareagr");
1366 }
1367
1368 set_units("m2", "m2");
1369 m_ts.variable().set_string("long_name", "area of grounded ice in glacierized areas");
1370 m_ts.variable().set_string("standard_name", "grounded_ice_sheet_area");
1371 m_ts.variable().set_number("valid_min", 0.0);
1372 }
1373
1374 double compute() {
1375 return ice_area_grounded(model->geometry(),
1376 m_config->get_number("output.ice_free_thickness_standard"));
1377 }
1378 };
1379
1380 //! \brief Computes the total floating ice area.
1381 class IceAreaGlacierizedShelf : public TSDiag<TSSnapshotDiagnostic, IceModel>
1382 {
1383 public:
1384 IceAreaGlacierizedShelf(IceModel *m)
1385 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_area_glacierized_floating") {
1386
1387 if (m_config->get_flag("output.ISMIP6")) {
1388 m_ts.variable().set_name("iareafl");
1389 }
1390
1391 set_units("m2", "m2");
1392 m_ts.variable().set_string("long_name", "area of ice shelves in glacierized areas");
1393 m_ts.variable().set_string("standard_name", "floating_ice_shelf_area");
1394 m_ts.variable().set_number("valid_min", 0.0);
1395 }
1396
1397 double compute() {
1398 return ice_area_floating(model->geometry(),
1399 m_config->get_number("output.ice_free_thickness_standard"));
1400 }
1401 };
1402
1403 //! \brief Computes the total grounded ice volume.
1404 class IceVolumeGlacierizedGrounded : public TSDiag<TSSnapshotDiagnostic, IceModel>
1405 {
1406 public:
1407 IceVolumeGlacierizedGrounded(IceModel *m)
1408 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_glacierized_grounded") {
1409
1410 set_units("m3", "m3");
1411 m_ts.variable().set_string("long_name", "volume of grounded ice in glacierized areas");
1412 m_ts.variable().set_number("valid_min", 0.0);
1413 }
1414
1415 double compute() {
1416 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
1417
1418 const IceModelVec2S &ice_thickness = model->geometry().ice_thickness;
1419
1420 const double
1421 thickness_threshold = m_config->get_number("output.ice_free_thickness_standard"),
1422 cell_area = m_grid->cell_area();
1423
1424 IceModelVec::AccessList list{&ice_thickness, &cell_type};
1425
1426 double volume = 0.0;
1427 for (Points p(*m_grid); p; p.next()) {
1428 const int i = p.i(), j = p.j();
1429
1430 const double H = ice_thickness(i, j);
1431
1432 if (cell_type.grounded(i, j) and H >= thickness_threshold) {
1433 volume += cell_area * H;
1434 }
1435 }
1436
1437 return GlobalSum(m_grid->com, volume);
1438 }
1439 };
1440
1441 //! \brief Computes the total floating ice volume.
1442 class IceVolumeGlacierizedShelf : public TSDiag<TSSnapshotDiagnostic, IceModel>
1443 {
1444 public:
1445 IceVolumeGlacierizedShelf(IceModel *m)
1446 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "ice_volume_glacierized_floating") {
1447
1448 set_units("m3", "m3");
1449 m_ts.variable().set_string("long_name", "volume of ice shelves in glacierized areas");
1450 m_ts.variable().set_number("valid_min", 0.0);
1451 }
1452
1453 double compute() {
1454 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
1455
1456 const IceModelVec2S &ice_thickness = model->geometry().ice_thickness;
1457
1458 const double
1459 thickness_threshold = m_config->get_number("output.ice_free_thickness_standard"),
1460 cell_area = m_grid->cell_area();
1461
1462 IceModelVec::AccessList list{&ice_thickness, &cell_type};
1463
1464 double volume = 0.0;
1465 for (Points p(*m_grid); p; p.next()) {
1466 const int i = p.i(), j = p.j();
1467
1468 const double H = ice_thickness(i, j);
1469
1470 if (cell_type.ocean(i, j) and H >= thickness_threshold) {
1471 volume += cell_area * H;
1472 }
1473 }
1474
1475 return GlobalSum(m_grid->com, volume);
1476 }
1477 };
1478
1479 //! \brief Reports the mass continuity time step.
1480 class TimeStepLength : public TSDiag<TSSnapshotDiagnostic, IceModel>
1481 {
1482 public:
1483 TimeStepLength(const IceModel *m)
1484 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "dt") {
1485
1486 set_units("second", "year");
1487 m_ts.variable().set_string("long_name", "mass continuity time step");
1488 m_ts.variable().set_number("valid_min", 0.0);
1489 }
1490
1491 double compute() {
1492 return model->dt();
1493 }
1494 };
1495
1496 //! \brief Reports maximum diffusivity.
1497 class MaxDiffusivity : public TSDiag<TSSnapshotDiagnostic, IceModel>
1498 {
1499 public:
1500 MaxDiffusivity(const IceModel *m)
1501 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "max_diffusivity") {
1502
1503 set_units("m2 s-1", "m2 s-1");
1504 m_ts.variable().set_string("long_name", "maximum diffusivity");
1505 m_ts.variable().set_number("valid_min", 0.0);
1506 }
1507
1508 double compute() {
1509 return model->stress_balance()->max_diffusivity();
1510 }
1511 };
1512
1513 //! \brief Reports the maximum horizontal absolute velocity component over the grid.
1514 /*!
1515 * This is the value used by the adaptive time-stepping code in the CFL condition
1516 * for horizontal advection (i.e. in energy and mass conservation time steps).
1517 *
1518 * This is not the maximum horizontal speed, but rather the maximum of components.
1519 *
1520 * Note that this picks up the value computed during the time-step taken at a
1521 * reporting time. (It is not the "average over the reporting interval computed using
1522 * differencing in time", as other rate-of-change diagnostics.)
1523 */
1524 class MaxHorizontalVelocity : public TSDiag<TSSnapshotDiagnostic, IceModel>
1525 {
1526 public:
1527 MaxHorizontalVelocity(const IceModel *m)
1528 : TSDiag<TSSnapshotDiagnostic, IceModel>(m, "max_hor_vel") {
1529
1530 set_units("m second-1", "m year-1");
1531 m_ts.variable().set_string("long_name",
1532 "maximum abs component of horizontal ice velocity"
1533 " over grid in last time step during time-series reporting interval");
1534 m_ts.variable().set_number("valid_min", 0.0);
1535 }
1536
1537 double compute() {
1538 CFLData cfl = model->stress_balance()->max_timestep_cfl_3d();
1539
1540 return std::max(cfl.u_max, cfl.v_max);
1541 }
1542 };
1543
1544 /*!
1545 * Return total mass change due to one of the terms in the mass continuity equation.
1546 *
1547 * Possible terms are
1548 *
1549 * - SMB: surface mass balance
1550 * - BMB: basal mass balance
1551 * - FLOW: ice flow
1552 * - ERROR: numerical flux needed to preserve non-negativity of thickness
1553 *
1554 * This computation can be restricted to grounded and floating areas
1555 * using the `area` argument.
1556 *
1557 * - BOTH: include all contributions
1558 * - GROUNDED: include grounded areas only
1559 * - SHELF: include floating areas only
1560 *
1561 * When computing mass changes due to flow it is important to remember
1562 * that ice mass in a cell can be represented by its thickness *or* an
1563 * "area specific volume". Transferring mass from one representation
1564 * to the other does not change the mass in a cell. This explains the
1565 * special case used when `term == FLOW`. (Note that surface and basal
1566 * mass balances do not affect the area specific volume field.)
1567 */
1568 double mass_change(const IceModel *model, TermType term, AreaType area) {
1569 const IceGrid &grid = *model->grid();
1570 const Config &config = *grid.ctx()->config();
1571
1572 const double
1573 ice_density = config.get_number("constants.ice.density"),
1574 cell_area = grid.cell_area();
1575
1576 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
1577
1578 const IceModelVec2S *thickness_change = nullptr;
1579
1580 switch (term) {
1581 case FLOW:
1582 thickness_change = &model->geometry_evolution().thickness_change_due_to_flow();
1583 break;
1584 case SMB:
1585 thickness_change = &model->geometry_evolution().top_surface_mass_balance();
1586 break;
1587 case BMB:
1588 thickness_change = &model->geometry_evolution().bottom_surface_mass_balance();
1589 break;
1590 case ERROR:
1591 thickness_change = &model->geometry_evolution().conservation_error();
1592 break;
1593 default:
1594 // can't happen
1595 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid term type");
1596 }
1597
1598 const IceModelVec2S &dV_flow = model->geometry_evolution().area_specific_volume_change_due_to_flow();
1599
1600 IceModelVec::AccessList list{&cell_type, thickness_change};
1601
1602 if (term == FLOW) {
1603 list.add(dV_flow);
1604 }
1605
1606 double volume_change = 0.0;
1607 for (Points p(grid); p; p.next()) {
1608 const int i = p.i(), j = p.j();
1609
1610 if ((area == BOTH) or
1611 (area == GROUNDED and cell_type.grounded(i, j)) or
1612 (area == SHELF and cell_type.ocean(i, j))) {
1613
1614 double dV = term == FLOW ? dV_flow(i, j) : 0.0;
1615
1616 // m^3 = m^2 * m
1617 volume_change += cell_area * ((*thickness_change)(i, j) + dV);
1618 }
1619 }
1620
1621 // (kg / m^3) * m^3 = kg
1622 return ice_density * GlobalSum(grid.com, volume_change);
1623 }
1624
1625 //! \brief Reports the total bottom surface ice flux.
1626 class IceMassFluxBasal : public TSDiag<TSFluxDiagnostic, IceModel>
1627 {
1628 public:
1629 IceMassFluxBasal(const IceModel *m)
1630 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_basal_mass_flux") {
1631
1632 if (m_config->get_flag("output.ISMIP6")) {
1633 m_ts.variable().set_name("tendlibmassbf");
1634 }
1635
1636 set_units("kg s-1", "Gt year-1");
1637 m_ts.variable().set_string("long_name", "total over ice domain of bottom surface ice mass flux");
1638 m_ts.variable().set_string("standard_name", "tendency_of_land_ice_mass_due_to_basal_mass_balance");
1639 m_ts.variable().set_string("comment", "positive means ice gain");
1640 }
1641
1642 double compute() {
1643 return mass_change(model, BMB, BOTH);
1644 }
1645 };
1646
1647 //! \brief Reports the total top surface ice flux.
1648 class IceMassFluxSurface : public TSDiag<TSFluxDiagnostic, IceModel>
1649 {
1650 public:
1651 IceMassFluxSurface(const IceModel *m)
1652 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_surface_mass_flux") {
1653
1654 if (m_config->get_flag("output.ISMIP6")) {
1655 m_ts.variable().set_name("tendacabf");
1656 }
1657
1658 set_units("kg s-1", "Gt year-1");
1659 m_ts.variable().set_string("long_name", "total over ice domain of top surface ice mass flux");
1660 m_ts.variable().set_string("standard_name", "tendency_of_land_ice_mass_due_to_surface_mass_balance");
1661 m_ts.variable().set_string("comment", "positive means ice gain");
1662 }
1663
1664 double compute() {
1665 return mass_change(model, SMB, BOTH);
1666 }
1667 };
1668
1669 //! \brief Reports the total basal ice flux over the grounded region.
1670 class IceMassFluxBasalGrounded : public TSDiag<TSFluxDiagnostic, IceModel>
1671 {
1672 public:
1673 IceMassFluxBasalGrounded(const IceModel *m)
1674 : TSDiag<TSFluxDiagnostic, IceModel>(m, "basal_mass_flux_grounded") {
1675
1676 set_units("kg s-1", "Gt year-1");
1677 m_ts.variable().set_string("long_name", "total over grounded ice domain of basal mass flux");
1678 m_ts.variable().set_string("standard_name", "tendency_of_land_ice_mass_due_to_basal_mass_balance");
1679 m_ts.variable().set_string("comment", "positive means ice gain");
1680 }
1681
1682 double compute() {
1683 return mass_change(model, BMB, GROUNDED);
1684 }
1685 };
1686
1687 //! \brief Reports the total sub-shelf ice flux.
1688 class IceMassFluxBasalFloating : public TSDiag<TSFluxDiagnostic, IceModel>
1689 {
1690 public:
1691 IceMassFluxBasalFloating(const IceModel *m)
1692 : TSDiag<TSFluxDiagnostic, IceModel>(m, "basal_mass_flux_floating") {
1693
1694 if (m_config->get_flag("output.ISMIP6")) {
1695 m_ts.variable().set_name("tendlibmassbffl");
1696 }
1697
1698 set_units("kg s-1", "Gt year-1");
1699 m_ts.variable().set_string("long_name", "total sub-shelf ice flux");
1700 m_ts.variable().set_string("standard_name", "tendency_of_land_ice_mass_due_to_basal_mass_balance");
1701 m_ts.variable().set_string("comment", "positive means ice gain");
1702 }
1703
1704 double compute() {
1705 return mass_change(model, BMB, SHELF);
1706 }
1707 };
1708
1709 //! \brief Reports the total numerical mass flux needed to preserve
1710 //! non-negativity of ice thickness.
1711 class IceMassFluxConservationError : public TSDiag<TSFluxDiagnostic, IceModel>
1712 {
1713 public:
1714 IceMassFluxConservationError(const IceModel *m)
1715 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_conservation_error") {
1716
1717 set_units("kg s-1", "Gt year-1");
1718 m_ts.variable().set_string("long_name", "total numerical flux needed to preserve non-negativity"
1719 " of ice thickness");
1720 m_ts.variable().set_string("comment", "positive means ice gain");
1721 }
1722
1723 double compute() {
1724 return mass_change(model, ERROR, BOTH);
1725 }
1726 };
1727
1728 //! \brief Reports the total discharge flux.
1729 class IceMassFluxDischarge : public TSDiag<TSFluxDiagnostic, IceModel>
1730 {
1731 public:
1732 IceMassFluxDischarge(const IceModel *m)
1733 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_discharge") {
1734
1735 if (m_config->get_flag("output.ISMIP6")) {
1736 m_ts.variable().set_name("tendlifmassbf");
1737 }
1738
1739 set_units("kg s-1", "Gt year-1");
1740 m_ts.variable().set_string("long_name",
1741 "discharge flux (frontal melt, calving, forced retreat)");
1742 m_ts.variable().set_string("standard_name",
1743 "tendency_of_land_ice_mass_due_to_calving_and_ice_front_melting");
1744 m_ts.variable().set_string("comment", "positive means ice gain");
1745 }
1746
1747 double compute() {
1748 const double ice_density = m_config->get_number("constants.ice.density");
1749
1750 const IceModelVec2S &calving = model->calving();
1751 const IceModelVec2S &frontal_melt = model->frontal_melt();
1752 const IceModelVec2S &forced_retreat = model->forced_retreat();
1753
1754 auto cell_area = m_grid->cell_area();
1755
1756 double volume_change = 0.0;
1757
1758 IceModelVec::AccessList list{&calving, &frontal_melt, &forced_retreat};
1759
1760 for (Points p(*m_grid); p; p.next()) {
1761 const int i = p.i(), j = p.j();
1762 // m^2 * m = m^3
1763 volume_change += cell_area * (calving(i, j) + frontal_melt(i, j) + forced_retreat(i, j));
1764 }
1765
1766 // (kg/m^3) * m^3 = kg
1767 return ice_density * GlobalSum(m_grid->com, volume_change);
1768 }
1769 };
1770
1771 //! \brief Reports the total calving flux.
1772 class IceMassFluxCalving : public TSDiag<TSFluxDiagnostic, IceModel>
1773 {
1774 public:
1775 IceMassFluxCalving(const IceModel *m)
1776 : TSDiag<TSFluxDiagnostic, IceModel>(m, "tendency_of_ice_mass_due_to_calving") {
1777
1778 if (m_config->get_flag("output.ISMIP6")) {
1779 m_ts.variable().set_name("tendlicalvf");
1780 }
1781
1782 set_units("kg s-1", "Gt year-1");
1783 m_ts.variable().set_string("long_name", "calving flux");
1784 m_ts.variable().set_string("standard_name", "tendency_of_land_ice_mass_due_to_calving");
1785 m_ts.variable().set_string("comment", "positive means ice gain");
1786 }
1787
1788 double compute() {
1789 const double ice_density = m_config->get_number("constants.ice.density");
1790
1791 const IceModelVec2S &calving = model->calving();
1792
1793 auto cell_area = m_grid->cell_area();
1794
1795 double volume_change = 0.0;
1796
1797 IceModelVec::AccessList list{&calving};
1798
1799 for (Points p(*m_grid); p; p.next()) {
1800 const int i = p.i(), j = p.j();
1801 // m^2 * m = m^3
1802 volume_change += cell_area * calving(i, j);
1803 }
1804
1805 // (kg/m^3) * m^3 = kg
1806 return ice_density * GlobalSum(m_grid->com, volume_change);
1807 }
1808 };
1809
1810 //! @brief Reports the total flux across the grounding line.
1811 class IceMassFluxAtGroundingLine : public TSDiag<TSFluxDiagnostic, IceModel>
1812 {
1813 public:
1814 IceMassFluxAtGroundingLine(const IceModel *m)
1815 : TSDiag<TSFluxDiagnostic, IceModel>(m, "grounding_line_flux") {
1816
1817 if (m_config->get_flag("output.ISMIP6")) {
1818 m_ts.variable().set_name("tendligroundf");
1819 m_ts.variable().set_string("standard_name",
1820 "tendency_of_grounded_ice_mass");
1821 }
1822
1823 set_units("kg s-1", "Gt year-1");
1824 m_ts.variable().set_string("long_name",
1825 "total ice flux across the grounding line");
1826 m_ts.variable().set_string("comment", "negative flux corresponds to ice loss into the ocean");
1827 }
1828
1829 double compute() {
1830 return total_grounding_line_flux(model->geometry().cell_type,
1831 model->geometry_evolution().flux_staggered(),
1832 model->dt());
1833 }
1834 };
1835
1836 } // end of namespace scalar
1837
1838
1839 //! \brief Computes dHdt, the ice thickness rate of change.
1840 class ThicknessRateOfChange : public Diag<IceModel>
1841 {
1842 public:
1843 ThicknessRateOfChange(const IceModel *m)
1844 : Diag<IceModel>(m),
1845 m_last_thickness(m_grid, "last_ice_thickness", WITHOUT_GHOSTS),
1846 m_interval_length(0.0) {
1847
1848 auto ismip6 = m_config->get_flag("output.ISMIP6");
1849
1850 // set metadata:
1851 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "dlithkdt" : "dHdt")};
1852
1853 set_attrs("ice thickness rate of change",
1854 "tendency_of_land_ice_thickness",
1855 "m s-1", "m year-1", 0);
1856
1857 auto large_number = to_internal(1e6);
1858
1859 m_vars[0].set_numbers("valid_range", {-large_number, large_number});
1860 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
1861 m_vars[0].set_string("cell_methods", "time: mean");
1862
1863 m_last_thickness.set_attrs("internal",
1864 "ice thickness at the time of the last report of dHdt",
1865 "m", "m", "land_ice_thickness", 0);
1866 }
1867 protected:
1868 IceModelVec::Ptr compute_impl() const {
1869
1870 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "dHdt", WITHOUT_GHOSTS));
1871 result->metadata() = m_vars[0];
1872
1873 if (m_interval_length > 0.0) {
1874 model->geometry().ice_thickness.add(-1.0, m_last_thickness, *result);
1875 result->scale(1.0 / m_interval_length);
1876 } else {
1877 result->set(m_fill_value);
1878 }
1879
1880 return result;
1881 }
1882
1883 void reset_impl() {
1884 m_interval_length = 0.0;
1885 m_last_thickness.copy_from(model->geometry().ice_thickness);
1886 }
1887
1888 void update_impl(double dt) {
1889 m_interval_length += dt;
1890 }
1891
1892 protected:
1893 IceModelVec2S m_last_thickness;
1894 double m_interval_length;
1895 };
1896
1897 //! \brief Computes latitude and longitude bounds.
1898 class LatLonBounds : public Diag<IceModel>
1899 {
1900 public:
1901 LatLonBounds(const IceModel *m,
1902 const std::string &var_name,
1903 const std::string &proj_string);
1904 protected:
1905 virtual IceModelVec::Ptr compute_impl() const;
1906 protected:
1907 std::string m_var_name, m_proj_string;
1908 };
1909
1910 LatLonBounds::LatLonBounds(const IceModel *m,
1911 const std::string &var_name,
1912 const std::string &proj_string)
1913 : Diag<IceModel>(m) {
1914 assert(var_name == "lat" || var_name == "lon");
1915 m_var_name = var_name;
1916
1917 // set metadata:
1918 std::vector<double> levels(4);
1919 for (int k = 0; k < 4; ++k) {
1920 levels[k] = k;
1921 }
1922
1923 m_vars = {SpatialVariableMetadata(m_sys, m_var_name + "_bnds", levels)};
1924 m_vars[0].get_z().set_name("nv4");
1925 m_vars[0].get_z().clear_all_strings();
1926 m_vars[0].get_z().clear_all_doubles();
1927 m_vars[0].set_time_independent(true);
1928
1929 if (m_var_name == "lon") {
1930 set_attrs("longitude bounds", "", "degree_east", "degree_east", 0);
1931 m_vars[0].set_numbers("valid_range", {-180, 180});
1932 } else {
1933 set_attrs("latitude bounds", "", "degree_north", "degree_north", 0);
1934 m_vars[0].set_numbers("valid_range", {-90, 90});
1935 }
1936
1937 m_proj_string = proj_string;
1938
1939 #if (Pism_USE_PROJ==1)
1940 // create the transformation from the provided projection to lat,lon to check if
1941 // proj_string is valid.
1942 Proj crs(m_proj_string, "EPSG:4326");
1943 #endif
1944 // If PISM_USE_PROJ is not 1 we don't need to check validity of m_proj_string: this diagnostic
1945 // will not be available and so this code will not run.
1946 }
1947
1948 IceModelVec::Ptr LatLonBounds::compute_impl() const {
1949 std::map<std::string,std::string> attrs;
1950 std::vector<double> indices(4);
1951
1952 IceModelVec3Custom::Ptr result(new IceModelVec3Custom(m_grid, m_var_name + "_bnds", "nv4",
1953 indices, attrs));
1954 result->metadata(0) = m_vars[0];
1955
1956 bool latitude = true;
1957 if (m_var_name == "lon") {
1958 latitude = false;
1959 }
1960
1961 if (latitude) {
1962 compute_lat_bounds(m_proj_string, *result);
1963 } else {
1964 compute_lon_bounds(m_proj_string, *result);
1965 }
1966
1967 return result;
1968 }
1969
1970 class IceAreaFraction : public Diag<IceModel>
1971 {
1972 public:
1973 IceAreaFraction(const IceModel *m);
1974 protected:
1975 virtual IceModelVec::Ptr compute_impl() const;
1976 };
1977
1978 IceAreaFraction::IceAreaFraction(const IceModel *m)
1979 : Diag<IceModel>(m) {
1980 m_vars = {SpatialVariableMetadata(m_sys, land_ice_area_fraction_name)};
1981 set_attrs("fraction of a grid cell covered by ice (grounded or floating)",
1982 "land_ice_area_fraction", // InitMIP "standard" name
1983 "1", "1", 0);
1984 }
1985
1986 IceModelVec::Ptr IceAreaFraction::compute_impl() const {
1987
1988 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, land_ice_area_fraction_name, WITHOUT_GHOSTS));
1989 result->metadata(0) = m_vars[0];
1990
1991 const IceModelVec2S
1992 &thickness = model->geometry().ice_thickness,
1993 &surface_elevation = model->geometry().ice_surface_elevation,
1994 &bed_topography = model->geometry().bed_elevation;
1995
1996 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
1997
1998 IceModelVec::AccessList list{&thickness, &surface_elevation, &bed_topography, &cell_type,
1999 result.get()};
2000
2001 const bool do_part_grid = m_config->get_flag("geometry.part_grid.enabled");
2002 const IceModelVec2S &Href = model->geometry().ice_area_specific_volume;;
2003 if (do_part_grid) {
2004 list.add(Href);
2005 }
2006
2007 ParallelSection loop(m_grid->com);
2008 try {
2009 for (Points p(*m_grid); p; p.next()) {
2010 const int i = p.i(), j = p.j();
2011
2012 if (cell_type.icy(i, j)) {
2013 // an "icy" cell: the area fraction is one
2014 (*result)(i, j) = 1.0;
2015 } else if (cell_type.ice_free_ocean(i, j)) {
2016 // an ice-free ocean cell may be "partially-filled", in which case we need to compute its
2017 // ice area fraction by dividing Href by the threshold thickness.
2018
2019 double H_reference = do_part_grid ? Href(i, j) : 0.0;
2020
2021 if (H_reference > 0.0) {
2022 const double H_threshold = part_grid_threshold_thickness(cell_type.int_star(i, j),
2023 thickness.star(i, j),
2024 surface_elevation.star(i, j),
2025 bed_topography(i,j));
2026 // protect from a division by zero
2027 if (H_threshold > 0.0) {
2028 (*result)(i, j) = std::min(H_reference / H_threshold, 1.0);
2029 } else {
2030 (*result)(i, j) = 1.0;
2031 }
2032 } else {
2033 // H_reference is zero
2034 (*result)(i, j) = 0.0;
2035 }
2036 } else {
2037 // an ice-free-ground cell: the area fraction is zero
2038 (*result)(i, j) = 0.0;
2039 }
2040 } // end of the loop over grid points
2041 } catch (...) {
2042 loop.failed();
2043 }
2044 loop.check();
2045
2046 return result;
2047 }
2048
2049 class IceAreaFractionGrounded : public Diag<IceModel>
2050 {
2051 public:
2052 IceAreaFractionGrounded(const IceModel *m);
2053 protected:
2054 virtual IceModelVec::Ptr compute_impl() const;
2055 };
2056
2057 IceAreaFractionGrounded::IceAreaFractionGrounded(const IceModel *m)
2058 : Diag<IceModel>(m) {
2059 m_vars = {SpatialVariableMetadata(m_sys, grounded_ice_sheet_area_fraction_name)};
2060 set_attrs("fraction of a grid cell covered by grounded ice",
2061 "grounded_ice_sheet_area_fraction", // InitMIP "standard" name
2062 "1", "1", 0);
2063 }
2064
2065 IceModelVec::Ptr IceAreaFractionGrounded::compute_impl() const {
2066 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, grounded_ice_sheet_area_fraction_name, WITHOUT_GHOSTS));
2067 result->metadata() = m_vars[0];
2068
2069 const double
2070 ice_density = m_config->get_number("constants.ice.density"),
2071 ocean_density = m_config->get_number("constants.sea_water.density");
2072
2073 auto
2074 &ice_thickness = model->geometry().ice_thickness,
2075 &sea_level = model->geometry().sea_level_elevation,
2076 &bed_topography = model->geometry().bed_elevation;
2077
2078 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
2079
2080 compute_grounded_cell_fraction(ice_density, ocean_density, sea_level,
2081 ice_thickness, bed_topography,
2082 *result);
2083
2084 // All grounded areas have the grounded cell fraction of one, so now we make sure that ice-free
2085 // areas get the value of 0 (they are grounded but not covered by a grounded ice sheet).
2086
2087 IceModelVec::AccessList list{&cell_type, result.get()};
2088
2089 ParallelSection loop(m_grid->com);
2090 try {
2091 for (Points p(*m_grid); p; p.next()) {
2092 const int i = p.i(), j = p.j();
2093 if (cell_type.ice_free(i, j)) {
2094 (*result)(i, j) = 0.0;
2095 }
2096 }
2097 } catch (...) {
2098 loop.failed();
2099 }
2100 loop.check();
2101
2102 return result;
2103 }
2104
2105 class IceAreaFractionFloating : public Diag<IceModel>
2106 {
2107 public:
2108 IceAreaFractionFloating(const IceModel *m);
2109 protected:
2110 virtual IceModelVec::Ptr compute_impl() const;
2111 };
2112
2113 IceAreaFractionFloating::IceAreaFractionFloating(const IceModel *m)
2114 : Diag<IceModel>(m) {
2115 m_vars = {SpatialVariableMetadata(m_sys, floating_ice_sheet_area_fraction_name)};
2116 set_attrs("fraction of a grid cell covered by floating ice",
2117 "floating_ice_shelf_area_fraction",
2118 "1", "1", 0);
2119 }
2120
2121 IceModelVec::Ptr IceAreaFractionFloating::compute_impl() const {
2122
2123 IceAreaFraction land_ice_area_fraction(model);
2124 IceModelVec::Ptr ice_area_fraction = land_ice_area_fraction.compute();
2125
2126 IceAreaFractionGrounded grounded_ice_sheet_area_fraction(model);
2127 IceModelVec::Ptr grounded_area_fraction = grounded_ice_sheet_area_fraction.compute();
2128
2129 IceModelVec::Ptr result = ice_area_fraction;
2130 result->metadata() = m_vars[0];
2131
2132 // Floating area fraction is total area fraction minus grounded area fraction.
2133 result->add(-1.0, *grounded_area_fraction);
2134
2135 return result;
2136 }
2137
2138 //! \brief Computes the 2D height above flotation.
2139 class HeightAboveFloatation : public Diag<IceModel>
2140 {
2141 public:
2142 HeightAboveFloatation(const IceModel *m);
2143 protected:
2144 virtual IceModelVec::Ptr compute_impl() const;
2145 };
2146
2147 HeightAboveFloatation::HeightAboveFloatation(const IceModel *m)
2148 : Diag<IceModel>(m) {
2149
2150 // set metadata:
2151 m_vars = {SpatialVariableMetadata(m_sys, "height_above_flotation")};
2152
2153 set_attrs("ice thickness in excess of the maximum floating ice thickness",
2154 "", "m", "m", 0);
2155 m_vars[0].set_number("_FillValue", m_fill_value);
2156 m_vars[0].set_string("comment",
2157 "shows how close to floatation the ice is at a given location");
2158 }
2159
2160 IceModelVec::Ptr HeightAboveFloatation::compute_impl() const {
2161
2162 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "height_above_flotation", WITHOUT_GHOSTS));
2163 result->metadata(0) = m_vars[0];
2164
2165 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
2166
2167 const double
2168 ice_density = m_config->get_number("constants.ice.density"),
2169 ocean_density = m_config->get_number("constants.sea_water.density");
2170
2171 auto
2172 &sea_level = model->geometry().sea_level_elevation,
2173 &ice_thickness = model->geometry().ice_thickness,
2174 &bed_topography = model->geometry().bed_elevation;
2175
2176 IceModelVec::AccessList list{&cell_type, result.get(), &ice_thickness, &bed_topography, &sea_level};
2177
2178 ParallelSection loop(m_grid->com);
2179 try {
2180 for (Points p(*m_grid); p; p.next()) {
2181 const int i = p.i(), j = p.j();
2182
2183 const double
2184 thickness = ice_thickness(i, j),
2185 bed = bed_topography(i, j),
2186 ocean_depth = sea_level(i, j) - bed;
2187
2188 if (cell_type.icy(i, j) and ocean_depth > 0.0) {
2189 const double max_floating_thickness = ocean_depth * (ocean_density / ice_density);
2190 (*result)(i, j) = thickness - max_floating_thickness;
2191 } else {
2192 (*result)(i, j) = m_fill_value;
2193 }
2194 }
2195 } catch (...) {
2196 loop.failed();
2197 }
2198 loop.check();
2199
2200 return result;
2201 }
2202
2203 //! \brief Computes the mass per cell.
2204 class IceMass : public Diag<IceModel>
2205 {
2206 public:
2207 IceMass(const IceModel *m);
2208 protected:
2209 virtual IceModelVec::Ptr compute_impl() const;
2210 };
2211
2212 IceMass::IceMass(const IceModel *m)
2213 : Diag<IceModel>(m) {
2214
2215 // set metadata:
2216 m_vars = {SpatialVariableMetadata(m_sys, "ice_mass")};
2217
2218 set_attrs("mass per cell",
2219 "", // no standard name
2220 "kg", "kg", 0);
2221 m_vars[0].set_number("_FillValue", m_fill_value);
2222 }
2223
2224 IceModelVec::Ptr IceMass::compute_impl() const {
2225
2226 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "ice_mass", WITHOUT_GHOSTS));
2227 result->metadata(0) = m_vars[0];
2228
2229 const IceModelVec2CellType &cell_type = model->geometry().cell_type;
2230
2231 const double
2232 ice_density = m_config->get_number("constants.ice.density");
2233
2234 const IceModelVec2S
2235 &ice_thickness = model->geometry().ice_thickness;
2236
2237 auto cell_area = m_grid->cell_area();
2238
2239 IceModelVec::AccessList list{&cell_type, result.get(), &ice_thickness};
2240
2241 ParallelSection loop(m_grid->com);
2242 try {
2243 for (Points p(*m_grid); p; p.next()) {
2244 const int i = p.i(), j = p.j();
2245
2246 // count all ice, including cells which have so little they
2247 // are considered "ice-free"
2248 if (ice_thickness(i, j) > 0.0) {
2249 (*result)(i,j) = ice_density * ice_thickness(i, j) * cell_area;
2250 } else {
2251 (*result)(i,j) = m_fill_value;
2252 }
2253 } // end of loop over grid points
2254
2255 } catch (...) {
2256 loop.failed();
2257 }
2258 loop.check();
2259
2260 // Add the mass of ice in Href:
2261 if (m_config->get_flag("geometry.part_grid.enabled")) {
2262 const IceModelVec2S &Href = model->geometry().ice_area_specific_volume;
2263 list.add(Href);
2264 for (Points p(*m_grid); p; p.next()) {
2265 const int i = p.i(), j = p.j();
2266
2267 if (ice_thickness(i, j) <= 0.0 and Href(i, j) > 0.0) {
2268 (*result)(i, j) = ice_density * Href(i, j) * cell_area;
2269 }
2270 }
2271 }
2272
2273 return result;
2274 }
2275
2276 /*! @brief Sea-level adjusted bed topography (zero at sea level). */
2277 class BedTopographySeaLevelAdjusted : public Diag<IceModel>
2278 {
2279 public:
2280 BedTopographySeaLevelAdjusted(const IceModel *m);
2281 protected:
2282 IceModelVec::Ptr compute_impl() const;
2283 };
2284
2285 BedTopographySeaLevelAdjusted::BedTopographySeaLevelAdjusted(const IceModel *m)
2286 : Diag<IceModel>(m) {
2287
2288 /* set metadata: */
2289 m_vars = {SpatialVariableMetadata(m_sys, "topg_sl_adjusted")};
2290
2291 set_attrs("sea-level adjusted bed topography (zero is at sea level)", "",
2292 "meters", "meters", 0);
2293 }
2294
2295 IceModelVec::Ptr BedTopographySeaLevelAdjusted::compute_impl() const {
2296
2297 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "topg_sl_adjusted", WITHOUT_GHOSTS));
2298 result->metadata(0) = m_vars[0];
2299
2300 auto
2301 &bed = model->geometry().bed_elevation,
2302 &sea_level = model->geometry().sea_level_elevation;
2303
2304 IceModelVec::AccessList list{&bed, &sea_level, result.get()};
2305
2306 for (Points p(*m_grid); p; p.next()) {
2307 const int i = p.i(), j = p.j();
2308
2309 (*result)(i, j) = bed(i, j) - sea_level(i, j);
2310 }
2311
2312 return result;
2313 }
2314
2315 /*! @brief Ice hardness computed using the SIA flow law. */
2316 class IceHardness : public Diag<IceModel>
2317 {
2318 public:
2319 IceHardness(const IceModel *m);
2320 protected:
2321 IceModelVec::Ptr compute_impl() const;
2322 };
2323
2324 IceHardness::IceHardness(const IceModel *m)
2325 : Diag<IceModel>(m) {
2326
2327 /* set metadata: */
2328 m_vars = {SpatialVariableMetadata(m_sys, "hardness", m_grid->z())};
2329
2330 const double power = 1.0 / m_config->get_number("stress_balance.sia.Glen_exponent");
2331 auto unitstr = pism::printf("Pa s%f", power);
2332
2333 set_attrs("ice hardness computed using the SIA flow law", "",
2334 unitstr, unitstr, 0);
2335 }
2336
2337 IceModelVec::Ptr IceHardness::compute_impl() const {
2338
2339 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "hardness", WITHOUT_GHOSTS));
2340 result->metadata(0) = m_vars[0];
2341
2342 EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
2343
2344 const IceModelVec3 &ice_enthalpy = model->energy_balance_model()->enthalpy();
2345 const IceModelVec2S &ice_thickness = model->geometry().ice_thickness;
2346
2347 const rheology::FlowLaw &flow_law = *model->stress_balance()->modifier()->flow_law();
2348
2349 IceModelVec::AccessList list{&ice_enthalpy, &ice_thickness, result.get()};
2350
2351 const unsigned int Mz = m_grid->Mz();
2352
2353 ParallelSection loop(m_grid->com);
2354 try {
2355 for (Points p(*m_grid); p; p.next()) {
2356 const int i = p.i(), j = p.j();
2357 const double *E = ice_enthalpy.get_column(i, j);
2358 const double H = ice_thickness(i, j);
2359
2360 double *hardness = result->get_column(i, j);
2361
2362 for (unsigned int k = 0; k < Mz; ++k) {
2363 const double depth = H - m_grid->z(k);
2364
2365 // EC->pressure() handles negative depths correctly
2366 const double pressure = EC->pressure(depth);
2367
2368 hardness[k] = flow_law.hardness(E[k], pressure);
2369 }
2370 }
2371 } catch (...) {
2372 loop.failed();
2373 }
2374 loop.check();
2375
2376 return result;
2377 }
2378
2379 /*! @brief Effective viscosity of ice (3D). */
2380 class IceViscosity : public Diag<IceModel>
2381 {
2382 public:
2383 IceViscosity(IceModel *m);
2384 protected:
2385 IceModelVec::Ptr compute_impl() const;
2386 };
2387
2388 IceViscosity::IceViscosity(IceModel *m)
2389 : Diag<IceModel>(m) {
2390
2391 /* set metadata: */
2392 m_vars = {SpatialVariableMetadata(m_sys, "effective_viscosity", m_grid->z())};
2393
2394 set_attrs("effective viscosity of ice", "",
2395 "Pascal second", "kPascal second", 0);
2396 m_vars[0].set_number("valid_min", 0);
2397 m_vars[0].set_number("_FillValue", m_fill_value);
2398 }
2399
2400 static inline double square(double x) {
2401 return x * x;
2402 }
2403
2404 IceModelVec::Ptr IceViscosity::compute_impl() const {
2405
2406 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "effective_viscosity", WITHOUT_GHOSTS));
2407 result->metadata(0) = m_vars[0];
2408
2409 IceModelVec3 W(m_grid, "wvel", WITH_GHOSTS);
2410
2411 using mask::ice_free;
2412
2413 EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
2414
2415 const rheology::FlowLaw &flow_law = *model->stress_balance()->modifier()->flow_law();
2416
2417 const IceModelVec2S &ice_thickness = model->geometry().ice_thickness;
2418
2419 const IceModelVec3
2420 &ice_enthalpy = model->energy_balance_model()->enthalpy(),
2421 &U = model->stress_balance()->velocity_u(),
2422 &V = model->stress_balance()->velocity_v(),
2423 &W_without_ghosts = model->stress_balance()->velocity_w();
2424
2425 W_without_ghosts.update_ghosts(W);
2426
2427 const unsigned int Mz = m_grid->Mz();
2428 const double
2429 dx = m_grid->dx(),
2430 dy = m_grid->dy();
2431 const std::vector<double> &z = m_grid->z();
2432
2433 const IceModelVec2CellType &mask = model->geometry().cell_type;
2434
2435 IceModelVec::AccessList list{&U, &V, &W, &ice_enthalpy, &ice_thickness, &mask, result.get()};
2436
2437 ParallelSection loop(m_grid->com);
2438 try {
2439 for (Points p(*m_grid); p; p.next()) {
2440 const int i = p.i(), j = p.j();
2441
2442 const double *E = ice_enthalpy.get_column(i, j);
2443 const double H = ice_thickness(i, j);
2444
2445 const double
2446 *u = U.get_column(i, j),
2447 *u_n = U.get_column(i, j + 1),
2448 *u_e = U.get_column(i + 1, j),
2449 *u_s = U.get_column(i, j - 1),
2450 *u_w = U.get_column(i - 1, j);
2451
2452 const double
2453 *v = V.get_column(i, j),
2454 *v_n = V.get_column(i, j + 1),
2455 *v_e = V.get_column(i + 1, j),
2456 *v_s = V.get_column(i, j - 1),
2457 *v_w = V.get_column(i - 1, j);
2458
2459 const double
2460 *w = W.get_column(i, j),
2461 *w_n = W.get_column(i, j + 1),
2462 *w_e = W.get_column(i + 1, j),
2463 *w_s = W.get_column(i, j - 1),
2464 *w_w = W.get_column(i - 1, j);
2465
2466 StarStencil<int> m = mask.int_star(i, j);
2467 const unsigned int
2468 east = ice_free(m.e) ? 0 : 1,
2469 west = ice_free(m.w) ? 0 : 1,
2470 south = ice_free(m.s) ? 0 : 1,
2471 north = ice_free(m.n) ? 0 : 1;
2472
2473 double *viscosity = result->get_column(i, j);
2474
2475 if (ice_free(m.ij)) {
2476 result->set_column(i, j, m_fill_value);
2477 continue;
2478 }
2479
2480 for (unsigned int k = 0; k < Mz; ++k) {
2481 const double depth = H - z[k];
2482
2483 if (depth < 0.0) {
2484 viscosity[k] = m_fill_value;
2485 continue;
2486 }
2487
2488 // EC->pressure() handles negative depths correctly
2489 const double pressure = EC->pressure(depth);
2490
2491 const double hardness = flow_law.hardness(E[k], pressure);
2492
2493 double u_x = 0.0, v_x = 0.0, w_x = 0.0;
2494 if (west + east > 0) {
2495 const double D = 1.0 / (dx * (west + east));
2496 u_x = D * (west * (u[k] - u_w[k]) + east * (u_e[k] - u[k]));
2497 v_x = D * (west * (v[k] - v_w[k]) + east * (v_e[k] - v[k]));
2498 w_x = D * (west * (w[k] - w_w[k]) + east * (w_e[k] - w[k]));
2499 }
2500
2501 double u_y = 0.0, v_y = 0.0, w_y = 0.0;
2502 if (south + north > 0) {
2503 const double D = 1.0 / (dy * (south + north));
2504 u_y = D * (south * (u[k] - u_s[k]) + north * (u_n[k] - u[k]));
2505 v_y = D * (south * (v[k] - v_s[k]) + north * (v_n[k] - v[k]));
2506 w_y = D * (south * (w[k] - w_s[k]) + north * (w_n[k] - w[k]));
2507 }
2508
2509 double
2510 u_z = 0.0,
2511 v_z = 0.0,
2512 w_z = 0.0;
2513
2514 if (k == 0) {
2515 const double dz = z[1] - z[0];
2516 u_z = (u[1] - u[0]) / dz;
2517 v_z = (v[1] - v[0]) / dz;
2518 w_z = (w[1] - w[0]) / dz;
2519 } else if (k == Mz - 1) {
2520 const double dz = z[Mz - 1] - z[Mz - 2];
2521 u_z = (u[Mz - 1] - u[Mz - 2]) / dz;
2522 v_z = (v[Mz - 1] - v[Mz - 2]) / dz;
2523 w_z = (w[Mz - 1] - w[Mz - 2]) / dz;
2524 } else {
2525 const double
2526 dz_p = z[k + 1] - z[k],
2527 dz_m = z[k] - z[k - 1];
2528 u_z = 0.5 * ((u[k + 1] - u[k]) / dz_p + (u[k] - u[k - 1]) / dz_m);
2529 v_z = 0.5 * ((v[k + 1] - v[k]) / dz_p + (v[k] - v[k - 1]) / dz_m);
2530 w_z = 0.5 * ((w[k + 1] - w[k]) / dz_p + (w[k] - w[k - 1]) / dz_m);
2531 }
2532
2533 // These should be "epsilon dot", but that's just too long.
2534 const double
2535 eps_xx = u_x,
2536 eps_yy = v_y,
2537 eps_zz = w_z,
2538 eps_xy = 0.5 * (u_y + v_x),
2539 eps_xz = 0.5 * (u_z + w_x),
2540 eps_yz = 0.5 * (v_z + w_y);
2541
2542 // The second invariant of the 3D strain rate tensor; see equation 4.8 in [@ref
2543 // GreveBlatter2009]. Unlike secondInvariant_2D(), this code does not make assumptions about
2544 // the input velocity field: we do not ignore w_x and w_y and do not assume that u_z and v_z
2545 // are zero.
2546 const double
2547 gamma = (square(eps_xx) + square(eps_yy) + square(eps_zz) +
2548 2.0 * (square(eps_xy) + square(eps_xz) + square(eps_yz)));
2549
2550 double nu = 0.0;
2551 // Note: in PISM gamma has an extra factor of 1/2; compare to
2552 flow_law.effective_viscosity(hardness, 0.5 * gamma, &nu, NULL);
2553
2554 viscosity[k] = nu;
2555 }
2556 }
2557 } catch (...) {
2558 loop.failed();
2559 }
2560 loop.check();
2561
2562 return result;
2563 }
2564
2565 /*! @brief Report ice thickness */
2566 class IceThickness : public Diag<IceModel>
2567 {
2568 public:
2569 IceThickness(const IceModel *m)
2570 : Diag<IceModel>(m) {
2571
2572 auto ismip6 = m_config->get_flag("output.ISMIP6");
2573
2574 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "lithk" : "thk")};
2575
2576 set_attrs("land ice thickness", "land_ice_thickness",
2577 "m", "m", 0);
2578 m_vars[0].set_number("valid_min", 0.0);
2579 }
2580
2581 protected:
2582 IceModelVec::Ptr compute_impl() const {
2583
2584 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "thk", WITHOUT_GHOSTS));
2585 result->metadata(0) = m_vars[0];
2586
2587 result->copy_from(model->geometry().ice_thickness);
2588
2589 return result;
2590 }
2591 };
2592
2593 /*! @brief Report ice top surface elevation */
2594 class IceBottomSurfaceElevation : public Diag<IceModel>
2595 {
2596 public:
2597 IceBottomSurfaceElevation(const IceModel *m)
2598 : Diag<IceModel>(m) {
2599
2600 auto ismip6 = m_config->get_flag("output.ISMIP6");
2601
2602 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "base" : "ice_base_elevation")};
2603
2604 set_attrs("ice bottom surface elevation", "", // no standard name
2605 "m", "m", 0);
2606 }
2607
2608 protected:
2609 IceModelVec::Ptr compute_impl() const {
2610
2611 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "ice_base_elevation", WITHOUT_GHOSTS));
2612 result->metadata(0) = m_vars[0];
2613
2614 ice_bottom_surface(model->geometry(), *result);
2615
2616 return result;
2617 }
2618 };
2619
2620 /*! @brief Report ice top surface elevation */
2621 class IceSurfaceElevation : public Diag<IceModel>
2622 {
2623 public:
2624 IceSurfaceElevation(const IceModel *m)
2625 : Diag<IceModel>(m) {
2626
2627 auto ismip6 = m_config->get_flag("output.ISMIP6");
2628
2629 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "orog" : "usurf")};
2630
2631 set_attrs("ice top surface elevation", "surface_altitude",
2632 "m", "m", 0);
2633 }
2634
2635 protected:
2636 IceModelVec::Ptr compute_impl() const {
2637
2638 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "usurf", WITHOUT_GHOSTS));
2639 result->metadata(0) = m_vars[0];
2640
2641 result->copy_from(model->geometry().ice_surface_elevation);
2642
2643 return result;
2644 }
2645 };
2646
2647 /*! @brief Report grounding line flux. */
2648 class GroundingLineFlux : public DiagAverageRate<IceModel>
2649 {
2650 public:
2651 GroundingLineFlux(const IceModel *m)
2652 : DiagAverageRate<IceModel>(m, "grounding_line_flux", RATE) {
2653
2654 auto ismip6 = m_config->get_flag("output.ISMIP6");
2655
2656 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "ligroundf" : "grounding_line_flux")};
2657 m_accumulator.metadata().set_string("units", "kg m-2");
2658
2659 set_attrs("grounding line flux", "",
2660 "kg m-2 second-1", "kg m-2 year-1", 0);
2661 m_vars[0].set_string("cell_methods", "time: mean");
2662
2663 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
2664 m_vars[0].set_string("comment",
2665 "Positive flux corresponds to mass moving from the ocean to"
2666 " an icy grounded area. This convention makes it easier to compare"
2667 " grounding line flux to the total discharge into the ocean");
2668 }
2669
2670 protected:
2671 void update_impl(double dt) {
2672 grounding_line_flux(model->geometry().cell_type,
2673 model->geometry_evolution().flux_staggered(),
2674 dt,
2675 ADD_VALUES,
2676 m_accumulator);
2677
2678 m_interval_length += dt;
2679 }
2680 };
2681
2682
2683 } // end of namespace diagnostics
2684
2685 void IceModel::init_diagnostics() {
2686
2687 using namespace diagnostics;
2688
2689 typedef Diagnostic d;
2690 typedef Diagnostic::Ptr f; // "f" for "field"
2691 m_diagnostics = {
2692 // geometry
2693 {"cell_grounded_fraction", d::wrap(m_geometry.cell_grounded_fraction)},
2694 {"height_above_flotation", f(new HeightAboveFloatation(this))},
2695 {"ice_area_specific_volume", d::wrap(m_geometry.ice_area_specific_volume)},
2696 {"ice_mass", f(new IceMass(this))},
2697 {"lat", d::wrap(m_geometry.latitude)},
2698 {"lon", d::wrap(m_geometry.longitude)},
2699 {"mask", d::wrap(m_geometry.cell_type)},
2700 {"thk", f(new IceThickness(this))},
2701 {"topg_sl_adjusted", f(new BedTopographySeaLevelAdjusted(this))},
2702 {"usurf", f(new IceSurfaceElevation(this))},
2703 {"ice_base_elevation", f(new IceBottomSurfaceElevation(this))},
2704 {floating_ice_sheet_area_fraction_name, f(new IceAreaFractionFloating(this))},
2705 {grounded_ice_sheet_area_fraction_name, f(new IceAreaFractionGrounded(this))},
2706 {land_ice_area_fraction_name, f(new IceAreaFraction(this))},
2707
2708 // temperature, enthalpy, and liquid water fraction
2709 {"enthalpybase", f(new IceEnthalpyBasal(this))},
2710 {"enthalpysurf", f(new IceEnthalpySurface(this))},
2711 {"bedtoptemp", d::wrap(m_bedtoptemp)},
2712 {"cts", f(new CTS(this))},
2713 {"liqfrac", f(new LiquidFraction(this))},
2714 {"temp", f(new Temperature(this))},
2715 {"temp_pa", f(new TemperaturePA(this))},
2716 {"tempbase", f(new TemperatureBasal(this, BOTH))},
2717 {"temppabase", f(new TemperaturePABasal(this))},
2718 {"tempsurf", f(new TemperatureSurface(this))},
2719
2720 // rheology-related stuff
2721 {"tempicethk", f(new TemperateIceThickness(this))},
2722 {"tempicethk_basal", f(new TemperateIceThicknessBasal(this))},
2723 {"hardav", f(new HardnessAverage(this))},
2724 {"hardness", f(new IceHardness(this))},
2725 {"effective_viscosity", f(new IceViscosity(this))},
2726
2727 // boundary conditions
2728 {"ssa_bc_mask", d::wrap(m_ssa_dirichlet_bc_mask)},
2729 {"ssa_bc_vel", d::wrap(m_ssa_dirichlet_bc_values)},
2730 {"ice_margin_pressure_difference", f(new IceMarginPressureDifference(this))},
2731
2732 // balancing the books
2733 // tendency_of_ice_amount = (tendency_of_ice_amount_due_to_flow +
2734 // tendency_of_ice_amount_due_to_conservation_error +
2735 // tendency_of_ice_amount_due_to_surface_mass_balance +
2736 // tendency_of_ice_amount_due_to_basal_mass_balance +
2737 // tendency_of_ice_amount_due_to_discharge)
2738 {"tendency_of_ice_amount", f(new TendencyOfIceAmount(this, AMOUNT))},
2739 {"tendency_of_ice_amount_due_to_flow", f(new TendencyOfIceAmountDueToFlow(this, AMOUNT))},
2740 {"tendency_of_ice_amount_due_to_conservation_error", f(new ConservationErrorFlux(this, AMOUNT))},
2741 {"tendency_of_ice_amount_due_to_surface_mass_flux", f(new SurfaceFlux(this, AMOUNT))},
2742 {"tendency_of_ice_amount_due_to_basal_mass_flux", f(new BasalFlux(this, AMOUNT))},
2743 {"tendency_of_ice_amount_due_to_discharge", f(new DischargeFlux(this, AMOUNT))},
2744 {"tendency_of_ice_amount_due_to_calving", f(new CalvingFlux(this, AMOUNT))},
2745
2746 // same, in terms of mass
2747 // tendency_of_ice_mass = (tendency_of_ice_mass_due_to_flow +
2748 // tendency_of_ice_mass_due_to_conservation_error +
2749 // tendency_of_ice_mass_due_to_surface_mass_flux +
2750 // tendency_of_ice_mass_due_to_basal_mass_balance +
2751 // tendency_of_ice_mass_due_to_discharge)
2752 {"tendency_of_ice_mass", f(new TendencyOfIceAmount(this, MASS))},
2753 {"tendency_of_ice_mass_due_to_flow", f(new TendencyOfIceAmountDueToFlow(this, MASS))},
2754 {"tendency_of_ice_mass_due_to_conservation_error", f(new ConservationErrorFlux(this, MASS))},
2755 {"tendency_of_ice_mass_due_to_surface_mass_flux", f(new SurfaceFlux(this, MASS))},
2756 {"tendency_of_ice_mass_due_to_basal_mass_flux", f(new BasalFlux(this, MASS))},
2757 {"tendency_of_ice_mass_due_to_discharge", f(new DischargeFlux(this, MASS))},
2758 {"tendency_of_ice_mass_due_to_calving", f(new CalvingFlux(this, MASS))},
2759
2760 // other rates and fluxes
2761 {"basal_mass_flux_grounded", f(new BMBSplit(this, GROUNDED))},
2762 {"basal_mass_flux_floating", f(new BMBSplit(this, SHELF))},
2763 {"dHdt", f(new ThicknessRateOfChange(this))},
2764 {"bmelt", d::wrap(m_basal_melt_rate)},
2765 {"grounding_line_flux", f(new GroundingLineFlux(this))},
2766
2767 // misc
2768 {"rank", f(new Rank(this))},
2769 };
2770
2771 #if (Pism_USE_PROJ==1)
2772 std::string proj = m_grid->get_mapping_info().proj;
2773 if (not proj.empty()) {
2774 m_diagnostics["lat_bnds"] = f(new LatLonBounds(this, "lat", proj));
2775 m_diagnostics["lon_bnds"] = f(new LatLonBounds(this, "lon", proj));
2776 }
2777 #endif
2778
2779 // add ISMIP6 variable names
2780 if (m_config->get_flag("output.ISMIP6")) {
2781 m_diagnostics["base"] = m_diagnostics["ice_base_elevation"];
2782 m_diagnostics["lithk"] = m_diagnostics["thk"];
2783 m_diagnostics["dlithkdt"] = m_diagnostics["dHdt"];
2784 m_diagnostics["orog"] = m_diagnostics["usurf"];
2785 m_diagnostics["acabf"] = m_diagnostics["tendency_of_ice_amount_due_to_surface_mass_flux"];
2786 m_diagnostics["libmassbfgr"] = m_diagnostics["basal_mass_flux_grounded"];
2787 m_diagnostics["libmassbffl"] = m_diagnostics["basal_mass_flux_floating"];
2788 m_diagnostics["lifmassbf"] = m_diagnostics["tendency_of_ice_amount_due_to_discharge"];
2789 m_diagnostics["licalvf"] = m_diagnostics["tendency_of_ice_amount_due_to_calving"];
2790 m_diagnostics["litempbotgr"] = f(new TemperatureBasal(this, GROUNDED));
2791 m_diagnostics["litempbotfl"] = f(new TemperatureBasal(this, SHELF));
2792 m_diagnostics["ligroundf"] = m_diagnostics["grounding_line_flux"];
2793 }
2794
2795 typedef TSDiagnostic::Ptr s; // "s" for "scalar"
2796 m_ts_diagnostics = {
2797 // area
2798 {"ice_area_glacierized", s(new scalar::IceAreaGlacierized(this))},
2799 {"ice_area_glacierized_cold_base", s(new scalar::IceAreaGlacierizedColdBase(this))},
2800 {"ice_area_glacierized_grounded", s(new scalar::IceAreaGlacierizedGrounded(this))},
2801 {"ice_area_glacierized_floating", s(new scalar::IceAreaGlacierizedShelf(this))},
2802 {"ice_area_glacierized_temperate_base", s(new scalar::IceAreaGlacierizedTemperateBase(this))},
2803 // mass
2804 {"ice_mass_glacierized", s(new scalar::IceMassGlacierized(this))},
2805 {"ice_mass", s(new scalar::IceMass(this))},
2806 {"tendency_of_ice_mass_glacierized", s(new scalar::IceMassRateOfChangeGlacierized(this))},
2807 {"limnsw", s(new scalar::IceMassNotDisplacingSeaWater(this))},
2808 // volume
2809 {"ice_volume_glacierized", s(new scalar::IceVolumeGlacierized(this))},
2810 {"ice_volume_glacierized_cold", s(new scalar::IceVolumeGlacierizedCold(this))},
2811 {"ice_volume_glacierized_grounded", s(new scalar::IceVolumeGlacierizedGrounded(this))},
2812 {"ice_volume_glacierized_floating", s(new scalar::IceVolumeGlacierizedShelf(this))},
2813 {"ice_volume_glacierized_temperate", s(new scalar::IceVolumeGlacierizedTemperate(this))},
2814 {"ice_volume", s(new scalar::IceVolume(this))},
2815 {"ice_volume_cold", s(new scalar::IceVolumeCold(this))},
2816 {"ice_volume_temperate", s(new scalar::IceVolumeTemperate(this))},
2817 {"tendency_of_ice_volume_glacierized", s(new scalar::IceVolumeRateOfChangeGlacierized(this))},
2818 {"tendency_of_ice_volume", s(new scalar::IceVolumeRateOfChange(this))},
2819 {"sea_level_rise_potential", s(new scalar::SeaLevelRisePotential(this))},
2820 // energy
2821 {"ice_enthalpy_glacierized", s(new scalar::IceEnthalpyGlacierized(this))},
2822 {"ice_enthalpy", s(new scalar::IceEnthalpy(this))},
2823 // time-stepping
2824 {"max_diffusivity", s(new scalar::MaxDiffusivity(this))},
2825 {"max_hor_vel", s(new scalar::MaxHorizontalVelocity(this))},
2826 {"dt", s(new scalar::TimeStepLength(this))},
2827 // balancing the books
2828 {"tendency_of_ice_mass", s(new scalar::IceMassRateOfChange(this))},
2829 {"tendency_of_ice_mass_due_to_flow", s(new scalar::IceMassRateOfChangeDueToFlow(this))},
2830 {"tendency_of_ice_mass_due_to_conservation_error", s(new scalar::IceMassFluxConservationError(this))},
2831 {"tendency_of_ice_mass_due_to_basal_mass_flux", s(new scalar::IceMassFluxBasal(this))},
2832 {"tendency_of_ice_mass_due_to_surface_mass_flux", s(new scalar::IceMassFluxSurface(this))},
2833 {"tendency_of_ice_mass_due_to_discharge", s(new scalar::IceMassFluxDischarge(this))},
2834 {"tendency_of_ice_mass_due_to_calving", s(new scalar::IceMassFluxCalving(this))},
2835 // other fluxes
2836 {"basal_mass_flux_grounded", s(new scalar::IceMassFluxBasalGrounded(this))},
2837 {"basal_mass_flux_floating", s(new scalar::IceMassFluxBasalFloating(this))},
2838 {"grounding_line_flux", s(new scalar::IceMassFluxAtGroundingLine(this))},
2839 };
2840
2841 // add ISMIP6 variable names
2842 if (m_config->get_flag("output.ISMIP6")) {
2843 m_ts_diagnostics["iareafl"] = m_ts_diagnostics["ice_area_glacierized_floating"];
2844 m_ts_diagnostics["iareagr"] = m_ts_diagnostics["ice_area_glacierized_grounded"];
2845 m_ts_diagnostics["lim"] = m_ts_diagnostics["ice_mass"];
2846 m_ts_diagnostics["tendacabf"] = m_ts_diagnostics["tendency_of_ice_mass_due_to_surface_mass_flux"];
2847 m_ts_diagnostics["tendlibmassbf"] = m_ts_diagnostics["tendency_of_ice_mass_due_to_basal_mass_flux"];
2848 m_ts_diagnostics["tendlibmassbffl"] = m_ts_diagnostics["basal_mass_flux_floating"];
2849 m_ts_diagnostics["tendlicalvf"] = m_ts_diagnostics["tendency_of_ice_mass_due_to_calving"];
2850 m_ts_diagnostics["tendlifmassbf"] = m_ts_diagnostics["tendency_of_ice_mass_due_to_discharge"];
2851 m_ts_diagnostics["tendligroundf"] = m_ts_diagnostics["grounding_line_flux"];
2852 }
2853
2854 // get diagnostics from submodels
2855 for (auto m : m_submodels) {
2856 m_diagnostics = pism::combine(m_diagnostics, m.second->diagnostics());
2857 m_ts_diagnostics = pism::combine(m_ts_diagnostics, m.second->ts_diagnostics());
2858 }
2859 }
2860
2861 typedef std::map<std::string, std::vector<VariableMetadata>> Metadata;
2862
2863 static void print_diagnostics(const Logger &log, const Metadata &list) {
2864 for (const auto &d : list) {
2865 const std::string &name = d.first;
2866 log.message(1, " Name: %s\n", name.c_str());
2867
2868 for (const auto &v : d.second) {
2869
2870 std::string
2871 var_name = v.get_name(),
2872 units = v.get_string("units"),
2873 glaciological_units = v.get_string("glaciological_units"),
2874 long_name = v.get_string("long_name"),
2875 comment = v.get_string("comment");
2876
2877 if (not glaciological_units.empty()) {
2878 units = glaciological_units;
2879 }
2880
2881 log.message(1, " %s [%s]\n", var_name.c_str(), units.c_str());
2882 log.message(1, " %s\n", long_name.c_str());
2883 if (not comment.empty()) {
2884 log.message(1, " %s\n", comment.c_str());
2885 }
2886 }
2887 log.message(1, "\n");
2888 }
2889 }
2890
2891 static void print_diagnostics_json(const Logger &log, const Metadata &list) {
2892 log.message(1, "{\n");
2893 bool first_diagnostic = true;
2894 for (const auto &d : list) {
2895
2896 if (not first_diagnostic) {
2897 log.message(1, ",\n");
2898 } else {
2899 first_diagnostic = false;
2900 }
2901
2902 log.message(1, "\"%s\" : [\n", d.first.c_str());
2903
2904 bool first_variable = true;
2905 for (const auto &variable : d.second) {
2906
2907 std::string
2908 var_name = variable.get_name(),
2909 units = variable.get_string("units"),
2910 glaciological_units = variable.get_string("glaciological_units"),
2911 long_name = variable.get_string("long_name"),
2912 standard_name = variable.get_string("standard_name"),
2913 comment = variable.get_string("comment");
2914
2915 if (not glaciological_units.empty()) {
2916 units = glaciological_units;
2917 }
2918
2919 if (not first_variable) {
2920 log.message(1, ",\n");
2921 } else {
2922 first_variable = false;
2923 }
2924
2925 log.message(1, "[\"%s\", \"%s\", \"%s\", \"%s\", \"%s\"]",
2926 var_name.c_str(), units.c_str(), long_name.c_str(), standard_name.c_str(), comment.c_str());
2927 }
2928 log.message(1, "]");
2929 }
2930 log.message(1, "}\n");
2931 }
2932
2933 /*!
2934 * Return metadata of 2D and 3D diagnostics.
2935 */
2936 static Metadata diag_metadata(const std::map<std::string,Diagnostic::Ptr> &diags) {
2937 Metadata result;
2938
2939 for (auto f : diags) {
2940 Diagnostic::Ptr diag = f.second;
2941
2942 for (unsigned int k = 0; k < diag->n_variables(); ++k) {
2943 result[f.first].push_back(diag->metadata(k));
2944 }
2945 }
2946
2947 return result;
2948 }
2949
2950 /*!
2951 * Return metadata of scalar diagnostics.
2952 */
2953 static Metadata ts_diag_metadata(const std::map<std::string,TSDiagnostic::Ptr> &ts_diags) {
2954 Metadata result;
2955
2956 for (auto d : ts_diags) {
2957 // always one variable per diagnostic
2958 result[d.first] = {d.second->metadata()};
2959 }
2960
2961 return result;
2962 }
2963
2964 void IceModel::list_diagnostics_json() const {
2965
2966 m_log->message(1, "{\n");
2967
2968 m_log->message(1, "\"spatial\" :\n");
2969 print_diagnostics_json(*m_log, diag_metadata(m_diagnostics));
2970
2971 m_log->message(1, ",\n"); // separator
2972
2973 m_log->message(1, "\"scalar\" :\n");
2974 print_diagnostics_json(*m_log, ts_diag_metadata(m_ts_diagnostics));
2975
2976 m_log->message(1, "}\n");
2977 }
2978
2979 void IceModel::list_diagnostics() const {
2980
2981 m_log->message(1, "\n");
2982 m_log->message(1, "======== Available 2D and 3D diagnostics ========\n");
2983
2984 print_diagnostics(*m_log, diag_metadata(m_diagnostics));
2985
2986 // scalar time-series
2987 m_log->message(1, "======== Available time-series ========\n");
2988
2989 print_diagnostics(*m_log, ts_diag_metadata(m_ts_diagnostics));
2990 }
2991
2992 /*!
2993 Computes fraction of the base which is melted.
2994
2995 Communication occurs here.
2996 */
2997 double IceModel::compute_temperate_base_fraction(double total_ice_area) {
2998
2999 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
3000
3001 auto cell_area = m_grid->cell_area();
3002
3003 double result = 0.0, meltarea = 0.0;
3004
3005 const IceModelVec3 &enthalpy = m_energy_model->enthalpy();
3006
3007 IceModelVec::AccessList list{&enthalpy, &m_geometry.cell_type, &m_geometry.ice_thickness};
3008 ParallelSection loop(m_grid->com);
3009 try {
3010 for (Points p(*m_grid); p; p.next()) {
3011 const int i = p.i(), j = p.j();
3012
3013 if (m_geometry.cell_type.icy(i, j)) {
3014 const double
3015 E_basal = enthalpy.get_column(i, j)[0],
3016 pressure = EC->pressure(m_geometry.ice_thickness(i,j)); // FIXME issue #15
3017 // accumulate area of base which is at melt point
3018 if (EC->is_temperate_relaxed(E_basal, pressure)) {
3019 meltarea += cell_area;
3020 }
3021 }
3022 }
3023 } catch (...) {
3024 loop.failed();
3025 }
3026 loop.check();
3027
3028 // convert from m2 to km2
3029 meltarea = units::convert(m_sys, meltarea, "m2", "km2");
3030 // communication
3031 result = GlobalSum(m_grid->com, meltarea);
3032
3033 // normalize fraction correctly
3034 if (total_ice_area > 0.0) {
3035 result = result / total_ice_area;
3036 } else {
3037 result = 0.0;
3038 }
3039 return result;
3040 }
3041
3042
3043 /*!
3044 Computes fraction of the ice which is as old as the start of the run (original).
3045 Communication occurs here.
3046 */
3047 double IceModel::compute_original_ice_fraction(double total_ice_volume) {
3048
3049 double result = -1.0; // result value if not age.enabled
3050
3051 if (not m_age_model) {
3052 return result; // leave now
3053 }
3054
3055 const double a = m_grid->cell_area() * 1e-3 * 1e-3, // area unit (km^2)
3056 currtime = m_time->current(); // seconds
3057
3058 const IceModelVec3 &ice_age = m_age_model->age();
3059
3060 IceModelVec::AccessList list{&m_geometry.cell_type, &m_geometry.ice_thickness, &ice_age};
3061
3062 const double one_year = units::convert(m_sys, 1.0, "year", "seconds");
3063 double original_ice_volume = 0.0;
3064
3065 // compute local original volume
3066 ParallelSection loop(m_grid->com);
3067 try {
3068 for (Points p(*m_grid); p; p.next()) {
3069 const int i = p.i(), j = p.j();
3070
3071 if (m_geometry.cell_type.icy(i, j)) {
3072 // accumulate volume of ice which is original
3073 const double *age = ice_age.get_column(i, j);
3074 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i,j));
3075 for (int k = 1; k <= ks; k++) {
3076 // ice in segment is original if it is as old as one year less than current time
3077 if (0.5 * (age[k - 1] + age[k]) > currtime - one_year) {
3078 original_ice_volume += a * 1.0e-3 * (m_grid->z(k) - m_grid->z(k - 1));
3079 }
3080 }
3081 }
3082 }
3083 } catch (...) {
3084 loop.failed();
3085 }
3086 loop.check();
3087
3088
3089 // communicate to turn into global original fraction
3090 result = GlobalSum(m_grid->com, original_ice_volume);
3091
3092 // normalize fraction correctly
3093 if (total_ice_volume > 0.0) {
3094 result = result / total_ice_volume;
3095 } else {
3096 result = 0.0;
3097 }
3098 return result;
3099 }
3100
3101
3102 //! Computes the temperate ice volume, in m^3.
3103 double IceModel::ice_volume_temperate(double thickness_threshold) const {
3104
3105 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
3106
3107 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
3108
3109 auto cell_area = m_grid->cell_area();
3110
3111 double volume = 0.0;
3112
3113 IceModelVec::AccessList list{&m_geometry.ice_thickness, &ice_enthalpy};
3114 ParallelSection loop(m_grid->com);
3115 try {
3116 for (Points p(*m_grid); p; p.next()) {
3117 const int i = p.i(), j = p.j();
3118
3119 if (m_geometry.ice_thickness(i,j) >= thickness_threshold) {
3120 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i,j));
3121 const double *Enth = ice_enthalpy.get_column(i,j);
3122
3123 for (int k = 0; k < ks; ++k) {
3124 if (EC->is_temperate_relaxed(Enth[k],EC->pressure(m_geometry.ice_thickness(i,j)))) { // FIXME issue #15
3125 volume += (m_grid->z(k + 1) - m_grid->z(k)) * cell_area;
3126 }
3127 }
3128
3129 if (EC->is_temperate_relaxed(Enth[ks],EC->pressure(m_geometry.ice_thickness(i,j)))) { // FIXME issue #15
3130 volume += (m_geometry.ice_thickness(i,j) - m_grid->z(ks)) * cell_area;
3131 }
3132 }
3133 }
3134 } catch (...) {
3135 loop.failed();
3136 }
3137 loop.check();
3138
3139
3140 return GlobalSum(m_grid->com, volume);
3141 }
3142
3143 //! Computes the cold ice volume, in m^3.
3144 double IceModel::ice_volume_cold(double thickness_threshold) const {
3145
3146 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
3147
3148 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
3149
3150 double volume = 0.0;
3151
3152 auto cell_area = m_grid->cell_area();
3153
3154 IceModelVec::AccessList list{&m_geometry.ice_thickness, &ice_enthalpy};
3155
3156 ParallelSection loop(m_grid->com);
3157 try {
3158 for (Points p(*m_grid); p; p.next()) {
3159 const int i = p.i(), j = p.j();
3160
3161 const double thickness = m_geometry.ice_thickness(i, j);
3162
3163 // count all ice, including cells which have so little they
3164 // are considered "ice-free"
3165 if (thickness >= thickness_threshold) {
3166 const int ks = m_grid->kBelowHeight(thickness);
3167 const double *Enth = ice_enthalpy.get_column(i, j);
3168
3169 for (int k=0; k<ks; ++k) {
3170 if (not EC->is_temperate_relaxed(Enth[k], EC->pressure(thickness))) { // FIXME issue #15
3171 volume += (m_grid->z(k+1) - m_grid->z(k)) * cell_area;
3172 }
3173 }
3174
3175 if (not EC->is_temperate_relaxed(Enth[ks], EC->pressure(thickness))) { // FIXME issue #15
3176 volume += (thickness - m_grid->z(ks)) * cell_area;
3177 }
3178 }
3179 }
3180 } catch (...) {
3181 loop.failed();
3182 }
3183 loop.check();
3184
3185
3186 return GlobalSum(m_grid->com, volume);
3187 }
3188
3189 //! Computes area of basal ice which is temperate, in m^2.
3190 double IceModel::ice_area_temperate(double thickness_threshold) const {
3191
3192 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
3193
3194 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
3195
3196 double area = 0.0;
3197
3198 auto cell_area = m_grid->cell_area();
3199
3200 IceModelVec::AccessList list{&m_geometry.ice_thickness, &ice_enthalpy};
3201 ParallelSection loop(m_grid->com);
3202 try {
3203 for (Points p(*m_grid); p; p.next()) {
3204 const int i = p.i(), j = p.j();
3205
3206 const double
3207 thickness = m_geometry.ice_thickness(i, j),
3208 basal_enthalpy = ice_enthalpy.get_column(i, j)[0];
3209
3210 if (thickness >= thickness_threshold and
3211 EC->is_temperate_relaxed(basal_enthalpy, EC->pressure(thickness))) { // FIXME issue #15
3212 area += cell_area;
3213 }
3214 }
3215 } catch (...) {
3216 loop.failed();
3217 }
3218 loop.check();
3219
3220 return GlobalSum(m_grid->com, area);
3221 }
3222
3223 //! Computes area of basal ice which is cold, in m^2.
3224 double IceModel::ice_area_cold(double thickness_threshold) const {
3225
3226 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
3227
3228 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
3229
3230 double area = 0.0;
3231
3232 auto cell_area = m_grid->cell_area();
3233
3234 IceModelVec::AccessList list{&ice_enthalpy, &m_geometry.ice_thickness};
3235 ParallelSection loop(m_grid->com);
3236 try {
3237 for (Points p(*m_grid); p; p.next()) {
3238 const int i = p.i(), j = p.j();
3239
3240 const double
3241 thickness = m_geometry.ice_thickness(i, j),
3242 basal_enthalpy = ice_enthalpy.get_column(i, j)[0];
3243
3244 if (thickness >= thickness_threshold and
3245 not EC->is_temperate_relaxed(basal_enthalpy, EC->pressure(thickness))) { // FIXME issue #15
3246 area += cell_area;
3247 }
3248 }
3249 } catch (...) {
3250 loop.failed();
3251 }
3252 loop.check();
3253
3254 return GlobalSum(m_grid->com, area);
3255 }
3256
3257 } // end of namespace pism