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