PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
SurfaceModel.cc
Go to the documentation of this file.
1 // Copyright (C) 2008-2021 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 
23 #include "pism/coupler/SurfaceModel.hh"
24 #include "pism/coupler/AtmosphereModel.hh"
25 #include "pism/util/io/File.hh"
26 #include "pism/util/Vars.hh"
27 #include "pism/util/Time.hh"
28 #include "pism/util/IceGrid.hh"
29 #include "pism/util/pism_options.hh"
30 #include "pism/util/iceModelVec.hh"
31 #include "pism/util/MaxTimestep.hh"
32 #include "pism/util/pism_utilities.hh"
33 #include "pism/util/Context.hh"
34 
35 namespace pism {
36 namespace surface {
37 
39  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "surface_layer_mass", WITHOUT_GHOSTS));
40 
41  result->set_attrs("climate_forcing", "mass held in the surface layer",
42  "kg", "kg", "", 0);
43 
44  result->metadata()["valid_min"] = {0.0};
45 
46  return result;
47 }
48 
50 
51  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "surface_layer_thickness", WITHOUT_GHOSTS));
52 
53  result->set_attrs("climate_forcing",
54  "thickness of the surface process layer at the top surface of the ice",
55  "m", "m", "", 0);
56 
57  result->metadata()["valid_min"] = {0.0};
58 
59  return result;
60 }
61 
63 
65  "ice_surface_liquid_water_fraction", WITHOUT_GHOSTS));
66 
67  result->set_attrs("climate_forcing",
68  "liquid water fraction of the ice at the top surface",
69  "1", "1", "", 0);
70 
71  result->metadata()["valid_range"] = {0.0, 1.0};
72 
73  return result;
74 }
75 
77 
78  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "climatic_mass_balance", WITHOUT_GHOSTS));
79 
80  result->set_attrs("climate_forcing",
81  "surface mass balance (accumulation/ablation) rate",
82  "kg m-2 second-1", "kg m-2 year-1",
83  "land_ice_surface_specific_mass_balance_flux", 0);
84 
85  Config::ConstPtr config = grid->ctx()->config();
86  const double smb_max = config->get_number("surface.given.smb_max", "kg m-2 second-1");
87 
88  result->metadata()["valid_range"] = {-smb_max, smb_max};
89 
90  return result;
91 }
92 
94 
95  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "ice_surface_temp", WITHOUT_GHOSTS));
96 
97  result->set_attrs("climate_forcing",
98  "temperature of the ice at the ice surface but below firn processes",
99  "Kelvin", "Kelvin", "", 0);
100 
101  result->metadata()["valid_range"] = {0.0, 323.15}; // [0C, 50C]
102 
103  return result;
104 }
105 
107 
108  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "surface_accumulation_flux", WITHOUT_GHOSTS));
109 
110  result->set_attrs("diagnostic",
111  "surface accumulation (precipitation minus rain)",
112  "kg m-2", "kg m-2", "", 0);
113 
114  return result;
115 }
116 
118 
119  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "surface_melt_flux", WITHOUT_GHOSTS));
120 
121  result->set_attrs("diagnostic",
122  "surface melt",
123  "kg m-2", "kg m-2", "", 0);
124 
125  return result;
126 }
127 
129 
130  IceModelVec2S::Ptr result(new IceModelVec2S(grid, "surface_runoff_flux", WITHOUT_GHOSTS));
131 
132  result->set_attrs("diagnostic",
133  "surface meltwater runoff",
134  "kg m-2", "kg m-2", "", 0);
135 
136  return result;
137 }
138 
140  : Component(grid) {
141 
148 
149  // default values
150  m_layer_thickness->set(0.0);
151  m_layer_mass->set(0.0);
152  m_liquid_water_fraction->set(0.0);
153  m_accumulation->set(0.0);
154  m_melt->set(0.0);
155  m_runoff->set(0.0);
156 }
157 
158 SurfaceModel::SurfaceModel(IceGrid::ConstPtr g, std::shared_ptr<SurfaceModel> input)
159  : Component(g) {
160  m_input_model = input;
161  // this is a modifier: allocate storage only if necessary (in derived classes)
162 }
163 
165  std::shared_ptr<atmosphere::AtmosphereModel> atmosphere)
166  : SurfaceModel(grid) { // this constructor will allocate storage
167 
168  m_atmosphere = atmosphere;
169 }
170 
171 
172 //! \brief Returns accumulation
173 /*!
174  * Basic surface models currently implemented in PISM do not model accumulation
175  */
177  return accumulation_impl();
178 }
179 
180 //! \brief Returns melt
181 /*!
182  * Basic surface models currently implemented in PISM do not model melt
183  */
185  return melt_impl();
186 }
187 
188 //! \brief Returns runoff
189 /*!
190  * Basic surface models currently implemented in PISM do not model runoff
191  */
193  return runoff_impl();
194 }
195 
197  return mass_flux_impl();
198 }
199 
201  return temperature_impl();
202 }
203 
204 //! \brief Returns the liquid water fraction of the ice at the top ice surface.
205 /*!
206  * Most PISM surface models return 0.
207  */
210 }
211 
212 //! \brief Returns mass held in the surface layer.
213 /*!
214  * Basic surface models currently implemented in PISM do not model the mass of
215  * the surface layer.
216  */
218  return layer_mass_impl();
219 }
220 
221 //! \brief Returns thickness of the surface layer. Could be used to compute surface
222 //! elevation as a sum of elevation of the top surface of the ice and surface layer (firn,
223 //! etc) thickness.
224 /*!
225  * Basic surface models currently implemented in PISM do not model surface
226  * layer thickness.
227  */
229  return layer_thickness_impl();
230 }
231 
233  if (m_input_model) {
234  return m_input_model->accumulation();
235  } else {
236  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
237  }
238 }
239 
241  if (m_input_model) {
242  return m_input_model->melt();
243  } else {
244  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
245  }
246 }
247 
249  if (m_input_model) {
250  return m_input_model->runoff();
251  } else {
252  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
253  }
254 }
255 
257  if (m_input_model) {
258  return m_input_model->mass_flux();
259  } else {
260  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
261  }
262 }
263 
265  if (m_input_model) {
266  return m_input_model->temperature();
267  } else {
268  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
269  }
270 }
271 
273  if (m_input_model) {
274  return m_input_model->liquid_water_fraction();
275  } else {
276  return *m_liquid_water_fraction;
277  }
278 }
279 
281  if (m_input_model) {
282  return m_input_model->layer_mass();
283  } else {
284  return *m_layer_mass;
285  }
286 }
287 
289  if (m_input_model) {
290  return m_input_model->layer_thickness();
291  } else {
292  return *m_layer_thickness;
293  }
294 }
295 
297  if (m_atmosphere) {
298  return m_atmosphere->ts_diagnostics();
299  }
300 
301  if (m_input_model) {
302  return m_input_model->ts_diagnostics();
303  }
304 
305  return {};
306 }
307 
308 void SurfaceModel::init(const Geometry &geometry) {
309  this->init_impl(geometry);
310 }
311 
312 void SurfaceModel::init_impl(const Geometry &geometry) {
313  if (m_atmosphere) {
314  m_atmosphere->init(geometry);
315  }
316 
317  if (m_input_model) {
318  m_input_model->init(geometry);
319  }
320 }
321 
322 void SurfaceModel::update(const Geometry &geometry, double t, double dt) {
323  this->update_impl(geometry, t, dt);
324 }
325 
326 void SurfaceModel::update_impl(const Geometry &geometry, double t, double dt) {
327  if (m_atmosphere) {
328  m_atmosphere->update(geometry, t, dt);
329  }
330 
331  if (m_input_model) {
332  m_input_model->update(geometry, t, dt);
333  }
334 }
335 
336 void SurfaceModel::define_model_state_impl(const File &output) const {
337  if (m_atmosphere) {
338  m_atmosphere->define_model_state(output);
339  }
340 
341  if (m_input_model) {
342  m_input_model->define_model_state(output);
343  }
344 }
345 
346 void SurfaceModel::write_model_state_impl(const File &output) const {
347  if (m_atmosphere) {
348  m_atmosphere->write_model_state(output);
349  }
350 
351  if (m_input_model) {
352  m_input_model->write_model_state(output);
353  }
354 }
355 
357  if (m_atmosphere) {
358  return m_atmosphere->max_timestep(t);
359  }
360 
361  if (m_input_model) {
362  return m_input_model->max_timestep(t);
363  }
364 
365  return MaxTimestep("surface model");
366 }
367 
368 /*!
369  * Use the surface mass balance to compute dummy accumulation.
370  *
371  * This is used by surface models that compute the SMB but do not provide accumulation,
372  * melt, and runoff.
373  *
374  * We assume that the positive part of the SMB is accumulation and the negative part is
375  * runoff. This ensures that outputs of PISM's surface models satisfy "SMB = accumulation
376  * - runoff".
377  */
379 
380  IceModelVec::AccessList list{&result, &smb};
381 
382  for (Points p(*m_grid); p; p.next()) {
383  const int i = p.i(), j = p.j();
384  result(i,j) = std::max(smb(i,j), 0.0);
385  }
386 }
387 
388 /*!
389  * Use the surface mass balance to compute dummy runoff.
390  *
391  * This is used by surface models that compute the SMB but do not provide accumulation,
392  * melt, and runoff.
393  *
394  * We assume that the positive part of the SMB is accumulation and the negative part is
395  * runoff. This ensures that outputs of PISM's surface models satisfy "SMB = accumulation
396  * - runoff".
397  */
399 
400  IceModelVec::AccessList list{&result, &smb};
401 
402  for (Points p(*m_grid); p; p.next()) {
403  const int i = p.i(), j = p.j();
404  result(i,j) = std::max(-smb(i,j), 0.0);
405  }
406 }
407 
408 /*!
409  * Use the surface mass balance to compute dummy runoff.
410  *
411  * This is used by surface models that compute the SMB but do not provide accumulation,
412  * melt, and runoff.
413  *
414  * We assume that all melt runs off, i.e. runoff = melt, but treat melt as a "derived"
415  * quantity.
416  */
418  dummy_runoff(smb, result);
419 }
420 
421 namespace diagnostics {
422 
423 // SurfaceModel diagnostics (these don't need to be in the header)
424 
425 /*! @brief Climatic mass balance */
426 class PS_climatic_mass_balance : public Diag<SurfaceModel>
427 {
428 public:
430 protected:
432 };
433 
434 /*! @brief Ice surface temperature. */
435 class PS_ice_surface_temp : public Diag<SurfaceModel>
436 {
437 public:
439 protected:
441 };
442 
443 /*! @brief Ice liquid water fraction at the ice surface. */
444 class PS_liquid_water_fraction : public Diag<SurfaceModel>
445 {
446 public:
448 protected:
450 };
451 
452 /*! @brief Mass of the surface layer (snow and firn). */
453 class PS_layer_mass : public Diag<SurfaceModel>
454 {
455 public:
456  PS_layer_mass(const SurfaceModel *m);
457 protected:
459 };
460 
461 /*! @brief Surface layer (snow and firn) thickness. */
462 class PS_layer_thickness : public Diag<SurfaceModel>
463 {
464 public:
466 protected:
468 };
469 
471  : Diag<SurfaceModel>(m) {
472 
473  /* set metadata: */
474  m_vars = {SpatialVariableMetadata(m_sys, "climatic_mass_balance")};
475 
476  set_attrs("surface mass balance (accumulation/ablation) rate",
477  "land_ice_surface_specific_mass_balance_flux",
478  "kg m-2 second-1", "kg m-2 year-1", 0);
479 }
480 
482 
483  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS));
484  result->metadata(0) = m_vars[0];
485 
486  result->copy_from(model->mass_flux());
487 
488  return result;
489 }
490 
492  : Diag<SurfaceModel>(m) {
493 
494 
495  auto ismip6 = m_config->get_flag("output.ISMIP6");
496 
497  /* set metadata: */
499  ismip6 ? "litemptop" : "ice_surface_temp")};
500 
501  set_attrs("ice temperature at the top ice surface",
502  "temperature_at_top_of_ice_sheet_model",
503  "K", "K", 0);
504 }
505 
507 
508  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "ice_surface_temp", WITHOUT_GHOSTS));
509  result->metadata(0) = m_vars[0];
510 
511  result->copy_from(model->temperature());
512 
513  return result;
514 }
515 
517  : Diag<SurfaceModel>(m) {
518 
519  /* set metadata: */
520  m_vars = {SpatialVariableMetadata(m_sys, "ice_surface_liquid_water_fraction")};
521 
522  set_attrs("ice liquid water fraction at the ice surface", "",
523  "1", "1", 0);
524 }
525 
527 
528  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "ice_surface_liquid_water_fraction", WITHOUT_GHOSTS));
529  result->metadata(0) = m_vars[0];
530 
531  result->copy_from(model->liquid_water_fraction());
532 
533  return result;
534 }
535 
537  : Diag<SurfaceModel>(m) {
538 
539  /* set metadata: */
540  m_vars = {SpatialVariableMetadata(m_sys, "surface_layer_mass")};
541 
542  set_attrs("mass of the surface layer (snow and firn)", "",
543  "kg", "kg", 0);
544 }
545 
547 
548  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "surface_layer_mass", WITHOUT_GHOSTS));
549  result->metadata(0) = m_vars[0];
550 
551  result->copy_from(model->layer_mass());
552 
553  return result;
554 }
555 
557  : Diag<SurfaceModel>(m) {
558 
559  /* set metadata: */
560  m_vars = {SpatialVariableMetadata(m_sys, "surface_layer_thickness")};
561 
562  set_attrs("thickness of the surface layer (snow and firn)", "",
563  "meters", "meters", 0);
564 }
565 
567 
568  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "surface_layer_thickness", WITHOUT_GHOSTS));
569  result->metadata(0) = m_vars[0];
570 
571  result->copy_from(model->layer_thickness());
572 
573  return result;
574 }
575 } // end of namespace diagnostics
576 
578  using namespace diagnostics;
579 
580  DiagnosticList result = {
581  {"climatic_mass_balance", Diagnostic::Ptr(new PS_climatic_mass_balance(this))},
582  {"ice_surface_temp", Diagnostic::Ptr(new PS_ice_surface_temp(this))},
583  {"ice_surface_liquid_water_fraction", Diagnostic::Ptr(new PS_liquid_water_fraction(this))},
584  {"surface_layer_mass", Diagnostic::Ptr(new PS_layer_mass(this))},
585  {"surface_layer_thickness", Diagnostic::Ptr(new PS_layer_thickness(this))}
586  };
587 
588  if (m_config->get_flag("output.ISMIP6")) {
589  result["litemptop"] = Diagnostic::Ptr(new PS_ice_surface_temp(this));
590  }
591 
592  if (m_atmosphere) {
593  result = pism::combine(result, m_atmosphere->diagnostics());
594  }
595 
596  if (m_input_model) {
597  result = pism::combine(result, m_input_model->diagnostics());
598  }
599 
600  return result;
601 }
602 
603 } // end of namespace surface
604 } // end of namespace pism
605 
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
DiagnosticList diagnostics() const
Definition: Component.cc:89
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:101
std::shared_ptr< const Config > ConstPtr
const SurfaceModel * model
Definition: Diagnostic.hh:166
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:161
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:112
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:64
IceGrid::ConstPtr m_grid
the grid
Definition: Diagnostic.hh:106
void set_attrs(const std::string &long_name, const std::string &standard_name, const std::string &units, const std::string &glaciological_units, unsigned int N=0)
A method for setting common variable attributes.
Definition: Diagnostic.cc:132
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:110
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
std::shared_ptr< IceModelVec2S > Ptr
Definition: iceModelVec.hh:341
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
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
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
static IceModelVec2S::Ptr allocate_mass_flux(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:76
const IceModelVec2S & liquid_water_fraction() const
Returns the liquid water fraction of the ice at the top ice surface.
IceModelVec2S::Ptr m_liquid_water_fraction
void update(const Geometry &geometry, double t, double dt)
std::shared_ptr< atmosphere::AtmosphereModel > m_atmosphere
const IceModelVec2S & temperature() const
virtual DiagnosticList diagnostics_impl() const
virtual const IceModelVec2S & temperature_impl() const
virtual void update_impl(const Geometry &geometry, double t, double dt)
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).
static IceModelVec2S::Ptr allocate_liquid_water_fraction(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:62
IceModelVec2S::Ptr m_melt
static IceModelVec2S::Ptr allocate_melt(IceGrid::ConstPtr grid)
const IceModelVec2S & layer_thickness() const
Returns thickness of the surface layer. Could be used to compute surface elevation as a sum of elevat...
virtual MaxTimestep max_timestep_impl(double my_t) const
SurfaceModel(IceGrid::ConstPtr g)
IceModelVec2S::Ptr m_layer_mass
virtual const IceModelVec2S & melt_impl() const
IceModelVec2S::Ptr m_accumulation
const IceModelVec2S & melt() const
Returns melt.
static IceModelVec2S::Ptr allocate_temperature(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:93
std::shared_ptr< SurfaceModel > m_input_model
void dummy_accumulation(const IceModelVec2S &smb, IceModelVec2S &result)
const IceModelVec2S & mass_flux() const
static IceModelVec2S::Ptr allocate_accumulation(IceGrid::ConstPtr grid)
virtual const IceModelVec2S & mass_flux_impl() const
static IceModelVec2S::Ptr allocate_layer_mass(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:38
const IceModelVec2S & layer_mass() const
Returns mass held in the surface layer.
virtual TSDiagnosticList ts_diagnostics_impl() const
static IceModelVec2S::Ptr allocate_runoff(IceGrid::ConstPtr grid)
const IceModelVec2S & runoff() const
Returns runoff.
virtual const IceModelVec2S & layer_mass_impl() const
virtual const IceModelVec2S & liquid_water_fraction_impl() const
const IceModelVec2S & accumulation() const
Returns accumulation.
IceModelVec2S::Ptr m_runoff
virtual void init_impl(const Geometry &geometry)
void dummy_runoff(const IceModelVec2S &smb, IceModelVec2S &result)
virtual const IceModelVec2S & runoff_impl() const
IceModelVec2S::Ptr m_layer_thickness
virtual const IceModelVec2S & accumulation_impl() const
static IceModelVec2S::Ptr allocate_layer_thickness(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:49
virtual const IceModelVec2S & layer_thickness_impl() const
void dummy_melt(const IceModelVec2S &smb, IceModelVec2S &result)
The interface of PISM's surface models.
Definition: SurfaceModel.hh:45
Mass of the surface layer (snow and firn).
Surface layer (snow and firn) thickness.
Ice liquid water fraction at the ice surface.
#define PISM_ERROR_LOCATION
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double g
Definition: exactTestP.cc:39
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:346
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
T combine(const T &a, const T &b)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49