PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
TemperatureIndex.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 PISM Authors
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> // std::min
20 
21 #include "TemperatureIndex.hh"
22 #include "localMassBalance.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/pism_options.hh"
25 #include "pism/util/Vars.hh"
26 #include "pism/util/Time.hh"
27 #include "pism/coupler/AtmosphereModel.hh"
28 #include "pism/util/Mask.hh"
29 #include "pism/util/io/File.hh"
30 
31 #include "pism/util/error_handling.hh"
32 #include "pism/util/io/io_helpers.hh"
33 #include "pism/util/pism_utilities.hh"
34 #include "pism/util/IceModelVec2CellType.hh"
35 #include "pism/geometry/Geometry.hh"
36 #include "pism/util/Context.hh"
37 
38 namespace pism {
39 namespace surface {
40 
41 ///// PISM surface model implementing a PDD scheme.
42 
44  std::shared_ptr<atmosphere::AtmosphereModel> input)
45  : SurfaceModel(g, input),
46  m_mass_flux(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
47  m_firn_depth(m_grid, "firn_depth", WITHOUT_GHOSTS),
48  m_snow_depth(m_grid, "snow_depth", WITHOUT_GHOSTS) {
49 
50  m_base_ddf.snow = m_config->get_number("surface.pdd.factor_snow");
51  m_base_ddf.ice = m_config->get_number("surface.pdd.factor_ice");
52  m_base_ddf.refreeze_fraction = m_config->get_number("surface.pdd.refreeze");
53  m_base_pddStdDev = m_config->get_number("surface.pdd.std_dev.value");
54  m_sd_use_param = m_config->get_flag("surface.pdd.std_dev.use_param");
55  m_sd_param_a = m_config->get_number("surface.pdd.std_dev.param_a");
56  m_sd_param_b = m_config->get_number("surface.pdd.std_dev.param_b");
57 
58  bool use_fausto_params = m_config->get_flag("surface.pdd.fausto.enabled");
59 
60  std::string method = m_config->get_string("surface.pdd.method");
61 
62  if (method == "repeatable_random_process") {
64  } else if (method == "random_process") {
66  } else {
68  }
69 
70  if (use_fausto_params) {
72  m_base_pddStdDev = 2.53;
73  }
74 
75  std::string sd_file = m_config->get_string("surface.pdd.std_dev.file");
76 
77  if (not sd_file.empty()) {
78  bool sd_periodic = m_config->get_flag("surface.pdd.std_dev.periodic");
79  int max_buffer_size = (unsigned int) m_config->get_number("input.forcing.buffer_size");
80 
81  File file(m_grid->com, sd_file, PISM_NETCDF3, PISM_READONLY);
83  "air_temp_sd", "",
84  max_buffer_size,
85  sd_periodic,
86  LINEAR);
87  m_sd_file_set = true;
88  } else {
89  m_air_temp_sd.reset(new IceModelVec2T(m_grid, "air_temp_sd", 1));
90  m_sd_file_set = false;
91  }
92 
93  m_air_temp_sd->set_attrs("climate_forcing",
94  "standard deviation of near-surface air temperature",
95  "Kelvin", "Kelvin", "", 0);
96 
97  m_mass_flux.set_attrs("diagnostic",
98  "instantaneous surface mass balance (accumulation/ablation) rate",
99  "kg m-2 s-1", "kg m-2 s-1",
100  "land_ice_surface_specific_mass_balance_flux", 0);
101 
102  m_mass_flux.metadata()["comment"] = "positive values correspond to ice gain";
103 
104  m_snow_depth.set_attrs("diagnostic",
105  "snow cover depth (set to zero once a year)",
106  "m", "m", "", 0);
107  m_snow_depth.set(0.0);
108 
109  m_firn_depth.set_attrs("diagnostic",
110  "firn cover depth",
111  "m", "m", "", 0);
112  m_firn_depth.metadata()["valid_min"] = {0.0};
113  m_firn_depth.set(0.0);
114 
116 
120 }
121 
122 void TemperatureIndex::init_impl(const Geometry &geometry) {
123 
124  // call the default implementation (not the interface method init())
125  SurfaceModel::init_impl(geometry);
126 
127  // report user's modeling choices
128  {
129  m_log->message(2,
130  "* Initializing the default temperature-index, PDD-based surface processes scheme.\n"
131  " Precipitation and 2m air temperature provided by atmosphere are inputs.\n"
132  " Surface mass balance and ice upper surface temperature are outputs.\n"
133  " See PISM User's Manual for control of degree-day factors.\n");
134 
135  m_log->message(2,
136  " Computing number of positive degree-days by: %s.\n",
137  m_mbscheme->method().c_str());
138 
139  if (m_faustogreve) {
140  m_log->message(2,
141  " Setting PDD parameters from [Faustoetal2009].\n");
142  } else {
143  m_log->message(2,
144  " Using default PDD parameters.\n");
145  }
146  }
147 
148  // initialize the spatially-variable air temperature standard deviation
149  {
150  std::string sd_file = m_config->get_string("surface.pdd.std_dev.file");
151  if (sd_file.empty()) {
152  m_log->message(2,
153  " Using constant standard deviation of near-surface temperature.\n");
154 
155  SpatialVariableMetadata attributes = m_air_temp_sd->metadata();
156  // replace with a constant IceModelVec2T:
158  // restore metadata:
159  m_air_temp_sd->metadata() = attributes;
160  } else {
161  m_log->message(2,
162  " Reading standard deviation of near-surface temperature from '%s'...\n",
163  sd_file.c_str());
164 
165  bool sd_periodic = m_config->get_flag("surface.pdd.std_dev.periodic");
166  m_air_temp_sd->init(sd_file, sd_periodic);
167  }
168  }
169 
170  // initializing the model state
172 
173  std::string firn_file = m_config->get_string("surface.pdd.firn_depth_file");
174 
175  if (input.type == INIT_RESTART) {
176  if (not firn_file.empty()) {
178  "surface.pdd.firn_depth_file is not allowed when"
179  " re-starting from a PISM output file.");
180  }
181 
182  m_firn_depth.read(input.filename, input.record);
183  m_snow_depth.read(input.filename, input.record);
184  } else if (input.type == INIT_BOOTSTRAP) {
185 
186  m_snow_depth.regrid(input.filename, OPTIONAL, 0.0);
187 
188  if (firn_file.empty()) {
189  m_firn_depth.regrid(input.filename, OPTIONAL, 0.0);
190  } else {
191  m_firn_depth.regrid(firn_file, CRITICAL);
192  }
193  } else {
194 
195  m_snow_depth.set(0.0);
196 
197  if (firn_file.empty()) {
198  m_firn_depth.set(0.0);
199  } else {
200  m_firn_depth.regrid(firn_file, CRITICAL);
201  }
202  }
203 
204  {
205  regrid("PDD surface model", m_snow_depth);
206  regrid("PDD surface model", m_firn_depth);
207  }
208 
209  // finish up
210  {
212 
213  m_accumulation->set(0.0);
214  m_melt->set(0.0);
215  m_runoff->set(0.0);
216  }
217 }
218 
220  return m_atmosphere->max_timestep(my_t);
221 }
222 
224  // compute the time corresponding to the beginning of the next balance year
225  double
226  balance_year_start_day = m_config->get_number("surface.pdd.balance_year_start_day"),
227  one_day = units::convert(m_sys, 1.0, "days", "seconds"),
228  year_start = m_grid->ctx()->time()->calendar_year_start(time),
229  balance_year_start = year_start + (balance_year_start_day - 1.0) * one_day;
230 
231  if (balance_year_start > time) {
232  return balance_year_start;
233  }
234  return m_grid->ctx()->time()->increment_date(balance_year_start, 1);
235 }
236 
237 void TemperatureIndex::update_impl(const Geometry &geometry, double t, double dt) {
238 
239  // make a copy of the pointer to convince clang static analyzer that its value does not
240  // change during the call
241  FaustoGrevePDDObject *fausto_greve = m_faustogreve.get();
242 
243  // update to ensure that temperature and precipitation time series are correct:
244  m_atmosphere->update(geometry, t, dt);
245 
246  m_temperature->copy_from(m_atmosphere->mean_annual_temp());
247 
248  // set up air temperature and precipitation time series
249  int N = m_mbscheme->get_timeseries_length(dt);
250 
251  const double dtseries = dt / N;
252  std::vector<double> ts(N), T(N), S(N), P(N), PDDs(N);
253  for (int k = 0; k < N; ++k) {
254  ts[k] = t + k * dtseries;
255  }
256 
257  // update standard deviation time series
258  if (m_sd_file_set) {
259  m_air_temp_sd->update(t, dt);
260  m_air_temp_sd->init_interpolation(ts);
261  }
262 
263  const IceModelVec2CellType &mask = geometry.cell_type;
264  const IceModelVec2S &H = geometry.ice_thickness;
265 
266  IceModelVec::AccessList list{&mask, &H, m_air_temp_sd.get(), &m_mass_flux,
268  m_accumulation.get(), m_melt.get(), m_runoff.get()};
269 
270  const double
271  sigmalapserate = m_config->get_number("surface.pdd.std_dev.lapse_lat_rate"),
272  sigmabaselat = m_config->get_number("surface.pdd.std_dev.lapse_lat_base");
273 
274  const IceModelVec2S *latitude = nullptr;
275  if (fausto_greve or sigmalapserate != 0.0) {
276  latitude = &geometry.latitude;
277 
278  list.add(*latitude);
279  }
280 
281  if (fausto_greve) {
282  const IceModelVec2S
283  *longitude = &geometry.latitude,
284  *surface_altitude = &geometry.ice_surface_elevation;
285 
286  fausto_greve->update_temp_mj(*surface_altitude, *latitude, *longitude);
287  }
288 
290 
291  m_atmosphere->init_timeseries(ts);
292 
293  m_atmosphere->begin_pointwise_access();
294 
295  const double ice_density = m_config->get_number("constants.ice.density");
296 
297  ParallelSection loop(m_grid->com);
298  try {
299  for (Points p(*m_grid); p; p.next()) {
300  const int i = p.i(), j = p.j();
301 
302  // the temperature time series from the AtmosphereModel and its modifiers
303  m_atmosphere->temp_time_series(i, j, T);
304 
305  if (mask.ice_free_ocean(i, j)) {
306  // ignore precipitation over ice-free ocean
307  for (int k = 0; k < N; ++k) {
308  P[k] = 0.0;
309  }
310  } else {
311  // elsewhere, get precipitation from the atmosphere model
312  m_atmosphere->precip_time_series(i, j, P);
313  }
314 
315  // convert precipitation from "kg m-2 second-1" to "m second-1" (PDDMassBalance expects
316  // accumulation in m/second ice equivalent)
317  for (int k = 0; k < N; ++k) {
318  P[k] = P[k] / ice_density;
319  // kg / (m^2 * second) / (kg / m^3) = m / second
320  }
321 
322  // interpolate temperature standard deviation time series
323  if (m_sd_file_set) {
324  m_air_temp_sd->interp(i, j, S);
325  } else {
326  double tmp = (*m_air_temp_sd)(i, j);
327  for (int k = 0; k < N; ++k) {
328  S[k] = tmp;
329  }
330  }
331 
332  if (fausto_greve) {
333  // we have been asked to set mass balance parameters according to
334  // formula (6) in [\ref Faustoetal2009]; they overwrite ddf set above
335  ddf = fausto_greve->degree_day_factors(i, j, (*latitude)(i, j));
336  }
337 
338  // apply standard deviation lapse rate on top of prescribed values
339  if (sigmalapserate != 0.0) {
340  double lat = (*latitude)(i, j);
341  for (int k = 0; k < N; ++k) {
342  S[k] += sigmalapserate * (lat - sigmabaselat);
343  }
344  (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
345  }
346 
347  // apply standard deviation param over ice if in use
348  if (m_sd_use_param and mask.icy(i, j)) {
349  for (int k = 0; k < N; ++k) {
350  S[k] = m_sd_param_a * (T[k] - 273.15) + m_sd_param_b;
351  if (S[k] < 0.0) {
352  S[k] = 0.0 ;
353  }
354  }
355  (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
356  }
357 
358  // Use temperature time series, the "positive" threshhold, and
359  // the standard deviation of the daily variability to get the
360  // number of positive degree days (PDDs)
361  if (mask.ice_free_ocean(i, j)) {
362  for (int k = 0; k < N; ++k) {
363  PDDs[k] = 0.0;
364  }
365  } else {
366  m_mbscheme->get_PDDs(dtseries, S, T, // inputs
367  PDDs); // output
368  }
369 
370  // Use temperature time series to remove rainfall from precipitation
371  m_mbscheme->get_snow_accumulation(T, // air temperature (input)
372  P); // precipitation rate (input-output)
373 
374  // Use degree-day factors, the number of PDDs, and the snow precipitation to get surface mass
375  // balance (and diagnostics: accumulation, melt, runoff)
376  {
377  double next_snow_depth_reset = m_next_balance_year_start;
378 
379  // make copies of firn and snow depth values at this point to avoid accessing 2D
380  // fields in the inner loop
381  double
382  ice = H(i, j),
383  firn = m_firn_depth(i, j),
384  snow = m_snow_depth(i, j);
385 
386  // accumulation, melt, runoff over this time-step
387  double
388  A = 0.0,
389  M = 0.0,
390  R = 0.0,
391  SMB = 0.0;
392 
393  for (int k = 0; k < N; ++k) {
394  if (ts[k] >= next_snow_depth_reset) {
395  snow = 0.0;
396  while (next_snow_depth_reset <= ts[k]) {
397  next_snow_depth_reset = m_grid->ctx()->time()->increment_date(next_snow_depth_reset, 1);
398  }
399  }
400 
401  const double accumulation = P[k] * dtseries;
402 
404  changes = m_mbscheme->step(ddf, PDDs[k],
405  ice, firn, snow, accumulation);
406 
407  // update ice thickness
408  ice += changes.smb;
409  assert(ice >= 0);
410 
411  // update firn depth
412  firn += changes.firn_depth;
413  assert(firn >= 0);
414 
415  // update snow depth
416  snow += changes.snow_depth;
417  assert(snow >= 0);
418 
419  // update total accumulation, melt, and runoff
420  {
421  A += accumulation;
422  M += changes.melt;
423  R += changes.runoff;
424  SMB += changes.smb;
425  }
426  } // end of the time-stepping loop
427 
428  // set firn and snow depths
429  m_firn_depth(i, j) = firn;
430  m_snow_depth(i, j) = snow;
431 
432  // set total accumulation, melt, and runoff, and SMB at this point, converting
433  // from "meters, ice equivalent" to "kg / m^2"
434  {
435  (*m_accumulation)(i, j) = A * ice_density;
436  (*m_melt)(i, j) = M * ice_density;
437  (*m_runoff)(i, j) = R * ice_density;
438  // m_mass_flux (unlike m_accumulation, m_melt, and m_runoff), is a
439  // rate. m * (kg / m^3) / second = kg / m^2 / second
440  m_mass_flux(i, j) = SMB * ice_density / dt;
441  }
442  }
443 
444  if (mask.ice_free_ocean(i, j)) {
445  m_firn_depth(i, j) = 0.0; // no firn in the ocean
446  m_snow_depth(i, j) = 0.0; // snow over the ocean does not stick
447  }
448  }
449  } catch (...) {
450  loop.failed();
451  }
452  loop.check();
453 
454  m_atmosphere->end_pointwise_access();
455 
457 }
458 
460  return m_mass_flux;
461 }
462 
464  return *m_temperature;
465 }
466 
468  return *m_accumulation;
469 }
470 
472  return *m_melt;
473 }
474 
476  return *m_runoff;
477 }
478 
480  return m_firn_depth;
481 }
482 
484  return m_snow_depth;
485 }
486 
488  return *m_air_temp_sd;
489 }
490 
495 }
496 
499  m_firn_depth.write(output);
500  m_snow_depth.write(output);
501 }
502 
503 namespace diagnostics {
504 
506 
507 /*! @brief Report surface melt, averaged over the reporting interval */
508 class SurfaceMelt : public DiagAverageRate<TemperatureIndex>
509 {
510 public:
513  kind == AMOUNT
514  ? "surface_melt_flux"
515  : "surface_melt_rate",
516  TOTAL_CHANGE),
517  m_melt_mass(m_grid, "melt_mass", WITHOUT_GHOSTS),
518  m_kind(kind)
519  {
520 
521  std::string
522  name = "surface_melt_flux",
523  long_name = "surface melt, averaged over the reporting interval",
524  standard_name = "surface_snow_and_ice_melt_flux",
525  accumulator_units = "kg m-2",
526  internal_units = "kg m-2 second-1",
527  external_units = "kg m-2 year-1";
528  if (kind == MASS) {
529  name = "surface_melt_rate";
530  standard_name = "";
531  accumulator_units = "kg",
532  internal_units = "kg second-1";
533  external_units = "Gt year-1" ;
534  }
535 
537  m_accumulator.metadata()["units"] = accumulator_units;
538 
539  set_attrs(long_name, standard_name, internal_units, external_units, 0);
540  m_vars[0]["cell_methods"] = "time: mean";
541 
542  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
543  }
544 
545 protected:
547  const IceModelVec2S &melt_amount = model->melt();
548 
549  if (m_kind == MASS) {
550  double cell_area = m_grid->cell_area();
551 
552  IceModelVec::AccessList list{&m_melt_mass, &melt_amount};
553 
554  for (Points p(*m_grid); p; p.next()) {
555  const int i = p.i(), j = p.j();
556  m_melt_mass(i, j) = melt_amount(i, j) * cell_area;
557  }
558  return m_melt_mass;
559  } else {
560  return melt_amount;
561  }
562  }
563 private:
566 };
567 
568 /*! @brief Report surface runoff, averaged over the reporting interval */
569 class SurfaceRunoff : public DiagAverageRate<TemperatureIndex>
570 {
571 public:
574  kind == AMOUNT
575  ? "surface_runoff_flux"
576  : "surface_runoff_rate",
577  TOTAL_CHANGE),
578  m_kind(kind),
579  m_runoff_mass(m_grid, "runoff_mass", WITHOUT_GHOSTS) {
580 
581  std::string
582  name = "surface_runoff_flux",
583  long_name = "surface runoff, averaged over the reporting interval",
584  standard_name = "surface_runoff_flux",
585  accumulator_units = "kg m-2",
586  internal_units = "kg m-2 second-1",
587  external_units = "kg m-2 year-1";
588  if (kind == MASS) {
589  name = "surface_runoff_rate";
590  standard_name = "",
591  accumulator_units = "kg",
592  internal_units = "kg second-1";
593  external_units = "Gt year-1" ;
594  }
595 
597  m_accumulator.metadata()["units"] = accumulator_units;
598 
599  set_attrs(long_name, standard_name, internal_units, external_units, 0);
600  m_vars[0]["cell_methods"] = "time: mean";
601 
602  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
603  }
604 
605 protected:
607  const IceModelVec2S &runoff_amount = model->runoff();
608 
609  if (m_kind == MASS) {
610  double cell_area = m_grid->cell_area();
611 
612  IceModelVec::AccessList list{&m_runoff_mass, &runoff_amount};
613 
614  for (Points p(*m_grid); p; p.next()) {
615  const int i = p.i(), j = p.j();
616  m_runoff_mass(i, j) = runoff_amount(i, j) * cell_area;
617  }
618  return m_runoff_mass;
619  } else {
620  return runoff_amount;
621  }
622  }
623 private:
626 };
627 
628 /*! @brief Report accumulation (precipitation minus rain), averaged over the reporting interval */
629 class Accumulation : public DiagAverageRate<TemperatureIndex>
630 {
631 public:
634  kind == AMOUNT
635  ? "surface_accumulation_flux"
636  : "surface_accumulation_rate",
637  TOTAL_CHANGE),
638  m_kind(kind),
639  m_accumulation_mass(m_grid, "accumulation_mass", WITHOUT_GHOSTS) {
640 
641  // possible standard name: surface_accumulation_flux
642  std::string
643  name = "surface_accumulation_flux",
644  long_name = "accumulation (precipitation minus rain), averaged over the reporting interval",
645  accumulator_units = "kg m-2",
646  internal_units = "kg m-2 second-1",
647  external_units = "kg m-2 year-1";
648  if (kind == MASS) {
649  name = "surface_accumulation_rate";
650  accumulator_units = "kg",
651  internal_units = "kg second-1";
652  external_units = "Gt year-1" ;
653  }
654 
655 
657  m_accumulator.metadata()["units"] = accumulator_units;
658 
659  set_attrs(long_name, "", internal_units, external_units, 0);
660  m_vars[0]["cell_methods"] = "time: mean";
661 
662  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
663  }
664 
665 protected:
667  const IceModelVec2S &accumulation_amount = model->accumulation();
668 
669  if (m_kind == MASS) {
670  double cell_area = m_grid->cell_area();
671 
672  IceModelVec::AccessList list{&m_accumulation_mass, &accumulation_amount};
673 
674  for (Points p(*m_grid); p; p.next()) {
675  const int i = p.i(), j = p.j();
676  m_accumulation_mass(i, j) = accumulation_amount(i, j) * cell_area;
677  }
678  return m_accumulation_mass;
679  } else {
680  return accumulation_amount;
681  }
682  }
683 private:
686 };
687 
688 /*!
689  * Integrate a field over the computational domain.
690  *
691  * If the input has units kg/m^2, the output will be in kg.
692  */
693 static double integrate(const IceModelVec2S &input) {
694  IceGrid::ConstPtr grid = input.grid();
695 
696  double cell_area = grid->cell_area();
697 
698  IceModelVec::AccessList list{&input};
699 
700  double result = 0.0;
701 
702  for (Points p(*grid); p; p.next()) {
703  const int i = p.i(), j = p.j();
704 
705  result += input(i, j) * cell_area;
706  }
707 
708  return GlobalSum(grid->com, result);
709 }
710 
711 
712 //! \brief Reports the total accumulation rate.
713 class TotalSurfaceAccumulation : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
714 {
715 public:
717  : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_accumulation_rate") {
718 
719  set_units("kg s-1", "kg year-1");
720  m_variable["long_name"] = "surface accumulation rate (PDD model)";
721  }
722 
723  double compute() {
724  return integrate(model->accumulation());
725  }
726 };
727 
728 
729 //! \brief Reports the total melt rate.
730 class TotalSurfaceMelt : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
731 {
732 public:
734  : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_melt_rate") {
735 
736  set_units("kg s-1", "kg year-1");
737  m_variable["long_name"] = "surface melt rate (PDD model)";
738  }
739 
740  double compute() {
741  return integrate(model->melt());
742  }
743 };
744 
745 
746 //! \brief Reports the total top surface ice flux.
747 class TotalSurfaceRunoff : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
748 {
749 public:
751  : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_runoff_rate") {
752 
753  set_units("kg s-1", "kg year-1");
754  m_variable["long_name"] = "surface runoff rate (PDD model)";
755  }
756 
757  double compute() {
758  return integrate(model->runoff());
759  }
760 };
761 
762 } // end of namespace diagnostics
763 
765  using namespace diagnostics;
766 
767  DiagnosticList result = {
768  {"surface_accumulation_flux", Diagnostic::Ptr(new Accumulation(this, AMOUNT))},
769  {"surface_accumulation_rate", Diagnostic::Ptr(new Accumulation(this, MASS))},
770  {"surface_melt_flux", Diagnostic::Ptr(new SurfaceMelt(this, AMOUNT))},
771  {"surface_melt_rate", Diagnostic::Ptr(new SurfaceMelt(this, MASS))},
772  {"surface_runoff_flux", Diagnostic::Ptr(new SurfaceRunoff(this, AMOUNT))},
773  {"surface_runoff_rate", Diagnostic::Ptr(new SurfaceRunoff(this, MASS))},
774  {"air_temp_sd", Diagnostic::wrap(*m_air_temp_sd)},
775  {"snow_depth", Diagnostic::wrap(m_snow_depth)},
776  {"firn_depth", Diagnostic::wrap(m_firn_depth)},
777  };
778 
779  result = pism::combine(result, SurfaceModel::diagnostics_impl());
780 
781  return result;
782 }
783 
785  using namespace diagnostics;
786 
787  TSDiagnosticList result = {
788  {"surface_accumulation_rate", TSDiagnostic::Ptr(new TotalSurfaceAccumulation(this))},
789  {"surface_melt_rate", TSDiagnostic::Ptr(new TotalSurfaceMelt(this))},
790  {"surface_runoff_rate", TSDiagnostic::Ptr(new TotalSurfaceRunoff(this))},
791  };
792 
794 
795  return result;
796 }
797 
798 } // end of namespace surface
799 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:140
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
virtual void regrid(const std::string &module_name, IceModelVec &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:151
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
DiagnosticList diagnostics() const
Definition: Component.cc:89
const Model * model
Definition: Diagnostic.hh:166
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
double m_fill_value
fill value (used often enough to justify storing it)
Definition: Diagnostic.hh:114
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
double to_internal(double x) const
Definition: Diagnostic.cc:58
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
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2S latitude
Definition: Geometry.hh:43
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
void add(double alpha, const IceModelVec2S &x)
static std::shared_ptr< IceModelVec2T > ForcingField(IceGrid::ConstPtr grid, const File &file, const std::string &short_name, const std::string &standard_name, int max_buffer_size, bool periodic, InterpolationType interpolation_type=PIECEWISE_CONSTANT)
static std::shared_ptr< IceModelVec2T > Constant(IceGrid::ConstPtr grid, const std::string &short_name, double value)
A class for storing and accessing 2D time-series (for climate forcing)
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void read(const std::string &filename, unsigned int time)
Definition: iceModelVec.cc:833
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void failed()
Indicates a failure of a parallel section.
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).
VariableMetadata m_variable
Definition: Diagnostic.hh:323
void set_units(const std::string &units, const std::string &glaciological_units)
Definition: Diagnostic.cc:207
std::shared_ptr< TSDiagnostic > Ptr
Definition: Diagnostic.hh:280
Scalar diagnostic reporting a "flux".
Definition: Diagnostic.hh:395
void update_temp_mj(const IceModelVec2S &surfelev, const IceModelVec2S &lat, const IceModelVec2S &lon)
Updates mean July near-surface air temperature.
LocalMassBalance::DegreeDayFactors degree_day_factors(int i, int j, double latitude)
A PDD implementation which computes the local mass balance based on an expectation integral.
An alternative PDD implementation which simulates a random process to get the number of PDDs.
std::shared_ptr< atmosphere::AtmosphereModel > m_atmosphere
virtual DiagnosticList diagnostics_impl() const
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_melt(IceGrid::ConstPtr grid)
static IceModelVec2S::Ptr allocate_temperature(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:93
static IceModelVec2S::Ptr allocate_accumulation(IceGrid::ConstPtr grid)
virtual TSDiagnosticList ts_diagnostics_impl() const
static IceModelVec2S::Ptr allocate_runoff(IceGrid::ConstPtr grid)
const IceModelVec2S & accumulation() const
Returns accumulation.
virtual void init_impl(const Geometry &geometry)
The interface of PISM's surface models.
Definition: SurfaceModel.hh:45
virtual const IceModelVec2S & runoff_impl() const
virtual const IceModelVec2S & accumulation_impl() const
IceModelVec2S m_firn_depth
firn depth
virtual TSDiagnosticList ts_diagnostics_impl() const
virtual void init_impl(const Geometry &geometry)
IceModelVec2S m_mass_flux
cached surface mass balance rate
virtual void update_impl(const Geometry &geometry, double t, double dt)
virtual const IceModelVec2S & mass_flux_impl() const
virtual const IceModelVec2S & melt_impl() const
TemperatureIndex(IceGrid::ConstPtr g, std::shared_ptr< atmosphere::AtmosphereModel > input)
std::shared_ptr< IceModelVec2T > m_air_temp_sd
standard deviation of the daily variability of the air temperature
std::unique_ptr< FaustoGrevePDDObject > m_faustogreve
if not NULL then user wanted fausto PDD stuff
virtual const IceModelVec2S & temperature_impl() const
double compute_next_balance_year_start(double time)
virtual MaxTimestep max_timestep_impl(double t) const
IceModelVec2S m_snow_depth
snow depth (reset once a year)
double m_base_pddStdDev
K; daily amount of randomness.
const IceModelVec2S & firn_depth() const
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
IceModelVec2S::Ptr m_runoff
total runoff during the last time step
const IceModelVec2S & air_temp_sd() const
const IceModelVec2S & snow_depth() const
IceModelVec2S::Ptr m_accumulation
total accumulation during the last time step
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual DiagnosticList diagnostics_impl() const
std::unique_ptr< LocalMassBalance > m_mbscheme
mass balance scheme to use
IceModelVec2S::Ptr m_melt
total melt during the last time step
LocalMassBalance::DegreeDayFactors m_base_ddf
holds degree-day factors in location-independent case
A class implementing a temperature-index (positive degree-day) scheme to compute melt and runoff,...
Accumulation(const TemperatureIndex *m, AmountKind kind)
Report accumulation (precipitation minus rain), averaged over the reporting interval.
SurfaceMelt(const TemperatureIndex *m, AmountKind kind)
Report surface melt, averaged over the reporting interval.
SurfaceRunoff(const TemperatureIndex *m, AmountKind kind)
Report surface runoff, averaged over the reporting interval.
Reports the total top surface ice flux.
#define PISM_ERROR_LOCATION
static double integrate(const IceModelVec2S &input)
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:72
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition: Component.cc:45
static const double g
Definition: exactTestP.cc:39
@ INIT_BOOTSTRAP
Definition: Component.hh:39
@ INIT_RESTART
Definition: Component.hh:39
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:346
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
@ OPTIONAL
Definition: IO_Flags.hh:70
@ CRITICAL
Definition: IO_Flags.hh:70
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
T combine(const T &a, const T &b)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
InitializationType type
initialization type
Definition: Component.hh:44
std::string filename
name of the input file (if applicable)
Definition: Component.hh:46
unsigned int record
index of the record to re-start from
Definition: Component.hh:48
double ice
m day^-1 K^-1; ice-equivalent amount of ice melted, per PDD
double refreeze_fraction
fraction of melted snow which refreezes as ice
double snow
m day^-1 K^-1; ice-equivalent amount of snow melted, per PDD
A struct which holds degree day factors.
static double S(unsigned n)
Definition: test_cube.c:58