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