PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
SurfaceModel.cc
Go to the documentation of this file.
1 // Copyright (C) 2008-2023 Ed Bueler, Constantine Khroulev, Ricarda Winkelmann,
2 // Gudfinna Adalgeirsdottir and Andy Aschwanden
3 //
4 // This file is part of PISM.
5 //
6 // PISM is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the Free Software
8 // Foundation; either version 3 of the License, or (at your option) any later
9 // version.
10 //
11 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 // details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with PISM; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 
20 #include <cassert>
21 #include <gsl/gsl_math.h> // GSL_NAN
22 #include <memory>
23 
24 #include "pism/coupler/SurfaceModel.hh"
25 #include "pism/coupler/AtmosphereModel.hh"
26 #include "pism/util/io/File.hh"
27 #include "pism/util/Grid.hh"
28 #include "pism/util/MaxTimestep.hh"
29 #include "pism/util/pism_utilities.hh"
30 #include "pism/util/Context.hh"
31 
32 namespace pism {
33 namespace surface {
34 
35 std::shared_ptr<array::Scalar> SurfaceModel::allocate_layer_mass(std::shared_ptr<const Grid> grid) {
36  auto result = std::make_shared<array::Scalar>(grid, "surface_layer_mass");
37 
38  result->metadata(0)
39  .long_name("mass held in the surface layer")
40  .units("kg");
41 
42  result->metadata()["valid_min"] = { 0.0 };
43 
44  return result;
45 }
46 
47 std::shared_ptr<array::Scalar> SurfaceModel::allocate_layer_thickness(std::shared_ptr<const Grid> grid) {
48 
49  auto result = std::make_shared<array::Scalar>(grid, "surface_layer_thickness");
50 
51  result->metadata(0)
52  .long_name("thickness of the surface process layer at the top surface of the ice")
53  .units("m");
54 
55  result->metadata()["valid_min"] = { 0.0 };
56 
57  return result;
58 }
59 
60 std::shared_ptr<array::Scalar> SurfaceModel::allocate_liquid_water_fraction(std::shared_ptr<const Grid> grid) {
61 
62  auto result = std::make_shared<array::Scalar>(grid, "ice_surface_liquid_water_fraction");
63 
64  result->metadata(0)
65  .long_name("liquid water fraction of the ice at the top surface")
66  .units("1");
67 
68  result->metadata()["valid_range"] = { 0.0, 1.0 };
69 
70  return result;
71 }
72 
73 std::shared_ptr<array::Scalar> SurfaceModel::allocate_mass_flux(std::shared_ptr<const Grid> grid) {
74 
75  auto result = std::make_shared<array::Scalar>(grid, "climatic_mass_balance");
76 
77  result->metadata(0)
78  .long_name("surface mass balance (accumulation/ablation) rate")
79  .units("kg m-2 second-1")
80  .output_units("kg m-2 year-1")
81  .standard_name("land_ice_surface_specific_mass_balance_flux");
82 
83  auto config = grid->ctx()->config();
84  const double smb_max = config->get_number("surface.given.smb_max", "kg m-2 second-1");
85 
86  result->metadata()["valid_range"] = { -smb_max, smb_max };
87 
88  return result;
89 }
90 
91 std::shared_ptr<array::Scalar>
92 SurfaceModel::allocate_temperature(std::shared_ptr<const Grid> grid) {
93 
94  auto result = std::make_shared<array::Scalar>(grid, "ice_surface_temp");
95 
96  result->metadata(0)
97  .long_name("temperature of the ice at the ice surface but below firn processes")
98  .units("Kelvin");
99 
100  result->metadata()["valid_range"] = { 0.0, 323.15 }; // [0C, 50C]
101 
102  return result;
103 }
104 
105 std::shared_ptr<array::Scalar>
106 SurfaceModel::allocate_accumulation(std::shared_ptr<const Grid> grid) {
107 
108  auto result = std::make_shared<array::Scalar>(grid, "surface_accumulation_flux");
109 
110  result->metadata(0).long_name("surface accumulation (precipitation minus rain)").units("kg m-2");
111 
112  return result;
113 }
114 
115 std::shared_ptr<array::Scalar> SurfaceModel::allocate_melt(std::shared_ptr<const Grid> grid) {
116 
117  auto result = std::make_shared<array::Scalar>(grid, "surface_melt_flux");
118 
119  result->metadata(0).long_name("surface melt").units("kg m-2");
120 
121  return result;
122 }
123 
124 std::shared_ptr<array::Scalar> SurfaceModel::allocate_runoff(std::shared_ptr<const Grid> grid) {
125 
126  auto result = std::make_shared<array::Scalar>(grid, "surface_runoff_flux");
127 
128  result->metadata(0).long_name("surface meltwater runoff").units("kg m-2");
129 
130  return result;
131 }
132 
133 SurfaceModel::SurfaceModel(std::shared_ptr<const Grid> grid)
134  : Component(grid) {
135 
142 
143  // default values
144  m_layer_thickness->set(0.0);
145  m_layer_mass->set(0.0);
146  m_liquid_water_fraction->set(0.0);
147  m_accumulation->set(0.0);
148  m_melt->set(0.0);
149  m_runoff->set(0.0);
150 }
151 
152 SurfaceModel::SurfaceModel(std::shared_ptr<const Grid> g, std::shared_ptr<SurfaceModel> input)
153  : Component(g) {
154  m_input_model = input;
155  // this is a modifier: allocate storage only if necessary (in derived classes)
156 }
157 
158 SurfaceModel::SurfaceModel(std::shared_ptr<const Grid> grid,
159  std::shared_ptr<atmosphere::AtmosphereModel> atmosphere)
160  : SurfaceModel(grid) { // this constructor will allocate storage
161 
162  m_atmosphere = atmosphere;
163 }
164 
165 
166 //! \brief Returns accumulation
167 /*!
168  * Basic surface models currently implemented in PISM do not model accumulation
169  */
171  return accumulation_impl();
172 }
173 
174 //! \brief Returns melt
175 /*!
176  * Basic surface models currently implemented in PISM do not model melt
177  */
179  return melt_impl();
180 }
181 
182 //! \brief Returns runoff
183 /*!
184  * Basic surface models currently implemented in PISM do not model runoff
185  */
187  return runoff_impl();
188 }
189 
191  return mass_flux_impl();
192 }
193 
195  return temperature_impl();
196 }
197 
198 //! \brief Returns the liquid water fraction of the ice at the top ice surface.
199 /*!
200  * Most PISM surface models return 0.
201  */
204 }
205 
206 //! \brief Returns mass held in the surface layer.
207 /*!
208  * Basic surface models currently implemented in PISM do not model the mass of
209  * the surface layer.
210  */
212  return layer_mass_impl();
213 }
214 
215 //! \brief Returns thickness of the surface layer. Could be used to compute surface
216 //! elevation as a sum of elevation of the top surface of the ice and surface layer (firn,
217 //! etc) thickness.
218 /*!
219  * Basic surface models currently implemented in PISM do not model surface
220  * layer thickness.
221  */
223  return layer_thickness_impl();
224 }
225 
227  if (m_input_model) {
228  return m_input_model->accumulation();
229  }
230 
231  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
232 }
233 
235  if (m_input_model) {
236  return m_input_model->melt();
237  }
238 
239  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
240 }
241 
243  if (m_input_model) {
244  return m_input_model->runoff();
245  }
246 
247  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
248 }
249 
251  if (m_input_model) {
252  return m_input_model->mass_flux();
253  }
254 
255  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
256 }
257 
259  if (m_input_model) {
260  return m_input_model->temperature();
261  }
262 
263  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
264 }
265 
267  if (m_input_model) {
268  return m_input_model->liquid_water_fraction();
269  }
270 
271  return *m_liquid_water_fraction;
272 }
273 
275  if (m_input_model) {
276  return m_input_model->layer_mass();
277  }
278 
279  return *m_layer_mass;
280 }
281 
283  if (m_input_model) {
284  return m_input_model->layer_thickness();
285  }
286 
287  return *m_layer_thickness;
288 }
289 
290 void SurfaceModel::init(const Geometry &geometry) {
291  this->init_impl(geometry);
292 }
293 
294 void SurfaceModel::init_impl(const Geometry &geometry) {
295  if (m_atmosphere) {
296  m_atmosphere->init(geometry);
297  }
298 
299  if (m_input_model) {
300  m_input_model->init(geometry);
301  }
302 }
303 
304 void SurfaceModel::update(const Geometry &geometry, double t, double dt) {
305  this->update_impl(geometry, t, dt);
306 }
307 
308 void SurfaceModel::update_impl(const Geometry &geometry, double t, double dt) {
309  if (m_atmosphere) {
310  m_atmosphere->update(geometry, t, dt);
311  }
312 
313  if (m_input_model) {
314  m_input_model->update(geometry, t, dt);
315  }
316 }
317 
318 void SurfaceModel::define_model_state_impl(const File &output) const {
319  if (m_atmosphere) {
320  m_atmosphere->define_model_state(output);
321  }
322 
323  if (m_input_model) {
324  m_input_model->define_model_state(output);
325  }
326 }
327 
328 void SurfaceModel::write_model_state_impl(const File &output) const {
329  if (m_atmosphere) {
330  m_atmosphere->write_model_state(output);
331  }
332 
333  if (m_input_model) {
334  m_input_model->write_model_state(output);
335  }
336 }
337 
339  if (m_atmosphere) {
340  return m_atmosphere->max_timestep(t);
341  }
342 
343  if (m_input_model) {
344  return m_input_model->max_timestep(t);
345  }
346 
347  return MaxTimestep("surface model");
348 }
349 
350 /*!
351  * Use the surface mass balance to compute dummy accumulation.
352  *
353  * This is used by surface models that compute the SMB but do not provide accumulation,
354  * melt, and runoff.
355  *
356  * We assume that the positive part of the SMB is accumulation and the negative part is
357  * runoff. This ensures that outputs of PISM's surface models satisfy "SMB = accumulation
358  * - runoff".
359  */
361 
362  array::AccessScope list{&result, &smb};
363 
364  for (auto p = m_grid->points(); p; p.next()) {
365  const int i = p.i(), j = p.j();
366  result(i,j) = std::max(smb(i,j), 0.0);
367  }
368 }
369 
370 /*!
371  * Use the surface mass balance to compute dummy runoff.
372  *
373  * This is used by surface models that compute the SMB but do not provide accumulation,
374  * melt, and runoff.
375  *
376  * We assume that the positive part of the SMB is accumulation and the negative part is
377  * runoff. This ensures that outputs of PISM's surface models satisfy "SMB = accumulation
378  * - runoff".
379  */
381 
382  array::AccessScope list{&result, &smb};
383 
384  for (auto p = m_grid->points(); p; p.next()) {
385  const int i = p.i(), j = p.j();
386  result(i,j) = std::max(-smb(i,j), 0.0);
387  }
388 }
389 
390 /*!
391  * Use the surface mass balance to compute dummy runoff.
392  *
393  * This is used by surface models that compute the SMB but do not provide accumulation,
394  * melt, and runoff.
395  *
396  * We assume that all melt runs off, i.e. runoff = melt, but treat melt as a "derived"
397  * quantity.
398  */
400  dummy_runoff(smb, result);
401 }
402 
403 namespace diagnostics {
404 
405 // SurfaceModel diagnostics (these don't need to be in the header)
406 
407 /*! @brief Climatic mass balance */
408 class PS_climatic_mass_balance : public Diag<SurfaceModel>
409 {
410 public:
412 protected:
413  std::shared_ptr<array::Array> compute_impl() const;
414 };
415 
416 /*! @brief Ice surface temperature. */
417 class PS_ice_surface_temp : public Diag<SurfaceModel>
418 {
419 public:
421 protected:
422  std::shared_ptr<array::Array> compute_impl() const;
423 };
424 
425 /*! @brief Ice liquid water fraction at the ice surface. */
426 class PS_liquid_water_fraction : public Diag<SurfaceModel>
427 {
428 public:
430 protected:
431  std::shared_ptr<array::Array> compute_impl() const;
432 };
433 
434 /*! @brief Mass of the surface layer (snow and firn). */
435 class PS_layer_mass : public Diag<SurfaceModel>
436 {
437 public:
438  PS_layer_mass(const SurfaceModel *m);
439 protected:
440  std::shared_ptr<array::Array> compute_impl() const;
441 };
442 
443 /*! @brief Surface layer (snow and firn) thickness. */
444 class PS_layer_thickness : public Diag<SurfaceModel>
445 {
446 public:
448 protected:
449  std::shared_ptr<array::Array> compute_impl() const;
450 };
451 
453  : Diag<SurfaceModel>(m) {
454 
455  /* set metadata: */
456  m_vars = { { m_sys, "climatic_mass_balance" } };
457  m_vars[0]
458  .long_name("surface mass balance (accumulation/ablation) rate")
459  .standard_name("land_ice_surface_specific_mass_balance_flux")
460  .units("kg m-2 second-1")
461  .output_units("kg m-2 year-1");
462 }
463 
464 std::shared_ptr<array::Array> PS_climatic_mass_balance::compute_impl() const {
465  auto result = allocate<array::Scalar>("climatic_mass_balance");
466 
467  result->copy_from(model->mass_flux());
468 
469  return result;
470 }
471 
473 
474  auto ismip6 = m_config->get_flag("output.ISMIP6");
475 
476  m_vars = { { m_sys, ismip6 ? "litemptop" : "ice_surface_temp" } };
477  m_vars[0]
478  .long_name("ice temperature at the top ice surface")
479  .standard_name("temperature_at_top_of_ice_sheet_model")
480  .units("K");
481 }
482 
483 std::shared_ptr<array::Array> PS_ice_surface_temp::compute_impl() const {
484  auto result = allocate<array::Scalar>("ice_surface_temp");
485 
486  result->copy_from(model->temperature());
487 
488  return result;
489 }
490 
492  m_vars = { { m_sys, "ice_surface_liquid_water_fraction" } };
493  m_vars[0].long_name("ice liquid water fraction at the ice surface").units("1");
494 }
495 
496 std::shared_ptr<array::Array> PS_liquid_water_fraction::compute_impl() const {
497 
498  auto result = allocate<array::Scalar>("ice_surface_liquid_water_fraction");
499 
500  result->copy_from(model->liquid_water_fraction());
501 
502  return result;
503 }
504 
506  m_vars = { { m_sys, "surface_layer_mass" } };
507  m_vars[0].long_name("mass of the surface layer (snow and firn)").units("kg");
508 }
509 
510 std::shared_ptr<array::Array> PS_layer_mass::compute_impl() const {
511 
512  auto result = allocate<array::Scalar>("surface_layer_mass");
513 
514  result->copy_from(model->layer_mass());
515 
516  return result;
517 }
518 
520  m_vars = { { m_sys, "surface_layer_thickness" } };
521  m_vars[0].long_name("thickness of the surface layer (snow and firn)").units("meters");
522 }
523 
524 std::shared_ptr<array::Array> PS_layer_thickness::compute_impl() const {
525 
526  auto result = allocate<array::Scalar>("surface_layer_thickness");
527 
528  result->copy_from(model->layer_thickness());
529 
530  return result;
531 }
532 
534 
535 /*! @brief Report surface melt, averaged over the reporting interval */
536 class SurfaceMelt : public DiagAverageRate<SurfaceModel>
537 {
538 public:
541  kind == AMOUNT
542  ? "surface_melt_flux"
543  : "surface_melt_rate",
544  TOTAL_CHANGE),
545  m_melt_mass(m_grid, "melt_mass"),
546  m_kind(kind)
547  {
548 
549  std::string
550  name = "surface_melt_flux",
551  long_name = "surface melt, averaged over the reporting interval",
552  standard_name = "surface_snow_and_ice_melt_flux",
553  accumulator_units = "kg m-2",
554  internal_units = "kg m-2 second-1",
555  external_units = "kg m-2 year-1";
556  if (kind == MASS) {
557  name = "surface_melt_rate";
558  standard_name = "";
559  accumulator_units = "kg",
560  internal_units = "kg second-1";
561  external_units = "Gt year-1" ;
562  }
563 
564  m_accumulator.metadata()["units"] = accumulator_units;
565 
566  m_vars = { { m_sys, name } };
567  m_vars[0]
568  .long_name(long_name)
569  .standard_name(standard_name)
570  .units(internal_units)
571  .output_units(external_units);
572  m_vars[0]["cell_methods"] = "time: mean";
573  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
574  }
575 
576 protected:
578  const array::Scalar &melt_amount = model->melt();
579 
580  if (m_kind == MASS) {
581  double cell_area = m_grid->cell_area();
582 
583  array::AccessScope list{&m_melt_mass, &melt_amount};
584 
585  for (auto p = m_grid->points(); p; p.next()) {
586  const int i = p.i(), j = p.j();
587  m_melt_mass(i, j) = melt_amount(i, j) * cell_area;
588  }
589  return m_melt_mass;
590  }
591 
592  return melt_amount;
593  }
594 private:
597 };
598 
599 /*! @brief Report surface runoff, averaged over the reporting interval */
600 class SurfaceRunoff : public DiagAverageRate<SurfaceModel>
601 {
602 public:
605  kind == AMOUNT
606  ? "surface_runoff_flux"
607  : "surface_runoff_rate",
608  TOTAL_CHANGE),
609  m_kind(kind),
610  m_runoff_mass(m_grid, "runoff_mass") {
611 
612  std::string
613  name = "surface_runoff_flux",
614  long_name = "surface runoff, averaged over the reporting interval",
615  standard_name = "surface_runoff_flux",
616  accumulator_units = "kg m-2",
617  internal_units = "kg m-2 second-1",
618  external_units = "kg m-2 year-1";
619  if (kind == MASS) {
620  name = "surface_runoff_rate";
621  standard_name = "",
622  accumulator_units = "kg",
623  internal_units = "kg second-1";
624  external_units = "Gt year-1" ;
625  }
626 
627  m_accumulator.metadata()["units"] = accumulator_units;
628 
629  m_vars = { { m_sys, name } };
630  m_vars[0]
631  .long_name(long_name)
632  .standard_name(standard_name)
633  .units(internal_units)
634  .output_units(external_units);
635  m_vars[0]["cell_methods"] = "time: mean";
636  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
637  }
638 
639 protected:
641  const array::Scalar &runoff_amount = model->runoff();
642 
643  if (m_kind == MASS) {
644  double cell_area = m_grid->cell_area();
645 
646  array::AccessScope list{&m_runoff_mass, &runoff_amount};
647 
648  for (auto p = m_grid->points(); p; p.next()) {
649  const int i = p.i(), j = p.j();
650  m_runoff_mass(i, j) = runoff_amount(i, j) * cell_area;
651  }
652  return m_runoff_mass;
653  }
654 
655  return runoff_amount;
656  }
657 private:
660 };
661 
662 /*! @brief Report accumulation (precipitation minus rain), averaged over the reporting interval */
663 class Accumulation : public DiagAverageRate<SurfaceModel>
664 {
665 public:
668  kind == AMOUNT
669  ? "surface_accumulation_flux"
670  : "surface_accumulation_rate",
671  TOTAL_CHANGE),
672  m_kind(kind),
673  m_accumulation_mass(m_grid, "accumulation_mass") {
674 
675  // possible standard name: surface_accumulation_flux
676  std::string
677  name = "surface_accumulation_flux",
678  long_name = "accumulation (precipitation minus rain), averaged over the reporting interval",
679  accumulator_units = "kg m-2",
680  internal_units = "kg m-2 second-1",
681  external_units = "kg m-2 year-1";
682  if (kind == MASS) {
683  name = "surface_accumulation_rate";
684  accumulator_units = "kg",
685  internal_units = "kg second-1";
686  external_units = "Gt year-1" ;
687  }
688 
689  m_accumulator.metadata()["units"] = accumulator_units;
690 
691  m_vars = { { m_sys, name } };
692  m_vars[0]
693  .long_name(long_name)
694  .units(internal_units)
695  .output_units(external_units);
696  m_vars[0]["cell_methods"] = "time: mean";
697  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
698  }
699 
700 protected:
702  const array::Scalar &accumulation_amount = model->accumulation();
703 
704  if (m_kind == MASS) {
705  double cell_area = m_grid->cell_area();
706 
707  array::AccessScope list{&m_accumulation_mass, &accumulation_amount};
708 
709  for (auto p = m_grid->points(); p; p.next()) {
710  const int i = p.i(), j = p.j();
711  m_accumulation_mass(i, j) = accumulation_amount(i, j) * cell_area;
712  }
713  return m_accumulation_mass;
714  }
715 
716  return accumulation_amount;
717  }
718 private:
721 };
722 
723 /*!
724  * Integrate a field over the computational domain.
725  *
726  * If the input has units kg/m^2, the output will be in kg.
727  */
728 static double integrate(const array::Scalar &input) {
729  auto grid = input.grid();
730 
731  double cell_area = grid->cell_area();
732 
733  array::AccessScope list{&input};
734 
735  double result = 0.0;
736 
737  for (auto p = grid->points(); p; p.next()) {
738  const int i = p.i(), j = p.j();
739 
740  result += input(i, j) * cell_area;
741  }
742 
743  return GlobalSum(grid->com, result);
744 }
745 
746 
747 //! \brief Reports the total accumulation rate.
748 class TotalSurfaceAccumulation : public TSDiag<TSFluxDiagnostic, SurfaceModel>
749 {
750 public:
752  : TSDiag<TSFluxDiagnostic, SurfaceModel>(m, "surface_accumulation_rate") {
753 
754  set_units("kg s-1", "kg year-1");
755  m_variable["long_name"] = "surface accumulation rate (PDD model)";
756  }
757 
758  double compute() {
759  return integrate(model->accumulation());
760  }
761 };
762 
763 
764 //! \brief Reports the total melt rate.
765 class TotalSurfaceMelt : public TSDiag<TSFluxDiagnostic, SurfaceModel>
766 {
767 public:
769  : TSDiag<TSFluxDiagnostic, SurfaceModel>(m, "surface_melt_rate") {
770 
771  set_units("kg s-1", "kg year-1");
772  m_variable["long_name"] = "surface melt rate (PDD model)";
773  }
774 
775  double compute() {
776  return integrate(model->melt());
777  }
778 };
779 
780 
781 //! \brief Reports the total top surface ice flux.
782 class TotalSurfaceRunoff : public TSDiag<TSFluxDiagnostic, SurfaceModel>
783 {
784 public:
786  : TSDiag<TSFluxDiagnostic, SurfaceModel>(m, "surface_runoff_rate") {
787 
788  set_units("kg s-1", "kg year-1");
789  m_variable["long_name"] = "surface runoff rate (PDD model)";
790  }
791 
792  double compute() {
793  return integrate(model->runoff());
794  }
795 };
796 
797 } // end of namespace diagnostics
798 
800  using namespace diagnostics;
801 
802  DiagnosticList result = {
803  {"surface_accumulation_flux", Diagnostic::Ptr(new Accumulation(this, AMOUNT))},
804  {"surface_accumulation_rate", Diagnostic::Ptr(new Accumulation(this, MASS))},
805  {"surface_melt_flux", Diagnostic::Ptr(new SurfaceMelt(this, AMOUNT))},
806  {"surface_melt_rate", Diagnostic::Ptr(new SurfaceMelt(this, MASS))},
807  {"surface_runoff_flux", Diagnostic::Ptr(new SurfaceRunoff(this, AMOUNT))},
808  {"surface_runoff_rate", Diagnostic::Ptr(new SurfaceRunoff(this, MASS))},
809  {"climatic_mass_balance", Diagnostic::Ptr(new PS_climatic_mass_balance(this))},
810  {"ice_surface_temp", Diagnostic::Ptr(new PS_ice_surface_temp(this))},
811  {"ice_surface_liquid_water_fraction", Diagnostic::Ptr(new PS_liquid_water_fraction(this))},
812  {"surface_layer_mass", Diagnostic::Ptr(new PS_layer_mass(this))},
813  {"surface_layer_thickness", Diagnostic::Ptr(new PS_layer_thickness(this))}
814  };
815 
816  if (m_config->get_flag("output.ISMIP6")) {
817  result["litemptop"] = Diagnostic::Ptr(new PS_ice_surface_temp(this));
818  }
819 
820  if (m_atmosphere) {
821  result = pism::combine(result, m_atmosphere->diagnostics());
822  }
823 
824  if (m_input_model) {
825  result = pism::combine(result, m_input_model->diagnostics());
826  }
827 
828  return result;
829 }
830 
832  using namespace diagnostics;
833 
834  TSDiagnosticList result = {
835  {"surface_accumulation_rate", TSDiagnostic::Ptr(new TotalSurfaceAccumulation(this))},
836  {"surface_melt_rate", TSDiagnostic::Ptr(new TotalSurfaceMelt(this))},
837  {"surface_runoff_rate", TSDiagnostic::Ptr(new TotalSurfaceRunoff(this))},
838  };
839 
840  if (m_atmosphere) {
841  return pism::combine(result, m_atmosphere->ts_diagnostics());
842  }
843 
844  if (m_input_model) {
845  return pism::combine(result, m_input_model->ts_diagnostics());
846  }
847 
848  return result;
849 }
850 
851 } // end of namespace surface
852 } // end of namespace pism
853 
std::shared_ptr< const Grid > grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
DiagnosticList diagnostics() const
Definition: Component.cc:89
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:118
const SurfaceModel * 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
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:118
High-level PISM I/O class.
Definition: File.hh:56
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
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
std::shared_ptr< TSDiagnostic > Ptr
Definition: Diagnostic.hh:280
Scalar diagnostic reporting a "flux".
Definition: Diagnostic.hh:395
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
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
const array::Scalar & melt() const
Returns melt.
static std::shared_ptr< array::Scalar > allocate_runoff(std::shared_ptr< const Grid > grid)
virtual const array::Scalar & layer_mass_impl() const
const array::Scalar & liquid_water_fraction() const
Returns the liquid water fraction of the ice at the top ice surface.
const array::Scalar & layer_mass() const
Returns mass held in the surface layer.
void update(const Geometry &geometry, double t, double dt)
std::shared_ptr< atmosphere::AtmosphereModel > m_atmosphere
static std::shared_ptr< array::Scalar > allocate_mass_flux(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:73
virtual DiagnosticList diagnostics_impl() const
void dummy_accumulation(const array::Scalar &smb, array::Scalar &result)
virtual void update_impl(const Geometry &geometry, double t, double dt)
std::shared_ptr< array::Scalar > m_melt
void init(const Geometry &geometry)
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
std::shared_ptr< array::Scalar > m_layer_thickness
const array::Scalar & mass_flux() const
static std::shared_ptr< array::Scalar > allocate_temperature(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:92
virtual const array::Scalar & accumulation_impl() const
const array::Scalar & accumulation() const
Returns accumulation.
virtual const array::Scalar & liquid_water_fraction_impl() const
static std::shared_ptr< array::Scalar > allocate_accumulation(std::shared_ptr< const Grid > grid)
virtual MaxTimestep max_timestep_impl(double my_t) const
static std::shared_ptr< array::Scalar > allocate_melt(std::shared_ptr< const Grid > grid)
SurfaceModel(std::shared_ptr< const Grid > g)
std::shared_ptr< array::Scalar > m_runoff
static std::shared_ptr< array::Scalar > allocate_layer_thickness(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:47
void dummy_melt(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< SurfaceModel > m_input_model
static std::shared_ptr< array::Scalar > allocate_layer_mass(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:35
const array::Scalar & temperature() const
std::shared_ptr< array::Scalar > m_layer_mass
virtual const array::Scalar & mass_flux_impl() const
virtual TSDiagnosticList ts_diagnostics_impl() const
virtual const array::Scalar & runoff_impl() const
std::shared_ptr< array::Scalar > m_accumulation
std::shared_ptr< array::Scalar > m_liquid_water_fraction
virtual void init_impl(const Geometry &geometry)
const array::Scalar & layer_thickness() const
Returns thickness of the surface layer. Could be used to compute surface elevation as a sum of elevat...
const array::Scalar & runoff() const
Returns runoff.
virtual const array::Scalar & melt_impl() const
virtual const array::Scalar & layer_thickness_impl() const
void dummy_runoff(const array::Scalar &smb, array::Scalar &result)
virtual const array::Scalar & temperature_impl() const
static std::shared_ptr< array::Scalar > allocate_liquid_water_fraction(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:60
The interface of PISM's surface models.
Definition: SurfaceModel.hh:42
Accumulation(const SurfaceModel *m, AmountKind kind)
Report accumulation (precipitation minus rain), averaged over the reporting interval.
std::shared_ptr< array::Array > compute_impl() const
std::shared_ptr< array::Array > compute_impl() const
std::shared_ptr< array::Array > compute_impl() const
Mass of the surface layer (snow and firn).
std::shared_ptr< array::Array > compute_impl() const
Surface layer (snow and firn) thickness.
std::shared_ptr< array::Array > compute_impl() const
Ice liquid water fraction at the ice surface.
SurfaceMelt(const SurfaceModel *m, AmountKind kind)
Report surface melt, averaged over the reporting interval.
SurfaceRunoff(const SurfaceModel *m, AmountKind kind)
Report surface runoff, averaged over the reporting interval.
Reports the total accumulation rate.
Reports the total top surface ice flux.
#define PISM_ERROR_LOCATION
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
static double integrate(const array::Scalar &input)
static const double g
Definition: exactTestP.cc:36
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:343
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
T combine(const T &a, const T &b)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)