tTemperatureIndex.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
---
tTemperatureIndex.cc (26899B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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
37 namespace pism {
38 namespace surface {
39
40 ///// PISM surface model implementing a PDD scheme.
41
42 TemperatureIndex::TemperatureIndex(IceGrid::ConstPtr g,
43 std::shared_ptr<atmosphere::AtmosphereModel> input)
44 : SurfaceModel(g, input),
45 m_mass_flux(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
46 m_firn_depth(m_grid, "firn_depth", WITHOUT_GHOSTS),
47 m_snow_depth(m_grid, "snow_depth", WITHOUT_GHOSTS) {
48
49 m_sd_period = m_config->get_number("surface.pdd.std_dev.period");
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");
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") {
63 m_mbscheme.reset(new PDDrandMassBalance(m_config, m_sys, PDDrandMassBalance::REPEATABLE));
64 } else if (method == "random_process") {
65 m_mbscheme.reset(new PDDrandMassBalance(m_config, m_sys, PDDrandMassBalance::NOT_REPEATABLE));
66 } else {
67 m_mbscheme.reset(new PDDMassBalance(m_config, m_sys));
68 }
69
70 if (use_fausto_params) {
71 m_faustogreve.reset(new FaustoGrevePDDObject(m_grid));
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 int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
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);
82 m_air_temp_sd = IceModelVec2T::ForcingField(m_grid, file,
83 "air_temp_sd", "",
84 max_buffer_size,
85 evaluations_per_year,
86 m_sd_period > 0,
87 LINEAR);
88 m_sd_file_set = true;
89 } else {
90 m_air_temp_sd.reset(new IceModelVec2T(m_grid, "air_temp_sd", 1, 1));
91 m_sd_file_set = false;
92 }
93
94 m_air_temp_sd->set_attrs("climate_forcing",
95 "standard deviation of near-surface air temperature",
96 "Kelvin", "Kelvin", "", 0);
97
98 m_mass_flux.set_attrs("diagnostic",
99 "instantaneous surface mass balance (accumulation/ablation) rate",
100 "kg m-2 s-1", "kg m-2 s-1",
101 "land_ice_surface_specific_mass_balance_flux", 0);
102
103 m_mass_flux.metadata().set_string("comment", "positive values correspond to ice gain");
104
105 m_snow_depth.set_attrs("diagnostic",
106 "snow cover depth (set to zero once a year)",
107 "m", "m", "", 0);
108 m_snow_depth.set(0.0);
109
110 m_firn_depth.set_attrs("diagnostic",
111 "firn cover depth",
112 "m", "m", "", 0);
113 m_firn_depth.metadata().set_number("valid_min", 0.0);
114 m_firn_depth.set(0.0);
115
116 m_temperature = allocate_temperature(g);
117
118 m_accumulation = allocate_accumulation(g);
119 m_melt = allocate_melt(g);
120 m_runoff = allocate_runoff(g);
121 }
122
123 TemperatureIndex::~TemperatureIndex() {
124 // empty
125 }
126
127 void TemperatureIndex::init_impl(const Geometry &geometry) {
128
129 // call the default implementation (not the interface method init())
130 SurfaceModel::init_impl(geometry);
131
132 // report user's modeling choices
133 {
134 m_log->message(2,
135 "* Initializing the default temperature-index, PDD-based surface processes scheme.\n"
136 " Precipitation and 2m air temperature provided by atmosphere are inputs.\n"
137 " Surface mass balance and ice upper surface temperature are outputs.\n"
138 " See PISM User's Manual for control of degree-day factors.\n");
139
140 m_log->message(2,
141 " Computing number of positive degree-days by: %s.\n",
142 m_mbscheme->method().c_str());
143
144 if (m_faustogreve) {
145 m_log->message(2,
146 " Setting PDD parameters from [Faustoetal2009].\n");
147 } else {
148 m_log->message(2,
149 " Using default PDD parameters.\n");
150 }
151 }
152
153 // initialize the spatially-variable air temperature standard deviation
154 {
155 std::string sd_file = m_config->get_string("surface.pdd.std_dev.file");
156 if (sd_file.empty()) {
157 m_log->message(2,
158 " Using constant standard deviation of near-surface temperature.\n");
159 m_air_temp_sd->init_constant(m_base_pddStdDev);
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 auto sd_ref_time = m_config->get_number("surface.pdd.std_dev.reference_year", "seconds");
166
167 m_air_temp_sd->init(sd_file, m_sd_period, sd_ref_time);
168 }
169 }
170
171 // initializing the model state
172 InputOptions input = process_input_options(m_grid->com, m_config);
173
174 std::string firn_file = m_config->get_string("surface.pdd.firn_depth_file");
175
176 if (input.type == INIT_RESTART) {
177 if (not firn_file.empty()) {
178 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
179 "surface.pdd.firn_depth_file is not allowed when"
180 " re-starting from a PISM output file.");
181 }
182
183 m_firn_depth.read(input.filename, input.record);
184 m_snow_depth.read(input.filename, input.record);
185 } else if (input.type == INIT_BOOTSTRAP) {
186
187 m_snow_depth.regrid(input.filename, OPTIONAL, 0.0);
188
189 if (firn_file.empty()) {
190 m_firn_depth.regrid(input.filename, OPTIONAL, 0.0);
191 } else {
192 m_firn_depth.regrid(firn_file, CRITICAL);
193 }
194 } else {
195
196 m_snow_depth.set(0.0);
197
198 if (firn_file.empty()) {
199 m_firn_depth.set(0.0);
200 } else {
201 m_firn_depth.regrid(firn_file, CRITICAL);
202 }
203 }
204
205 {
206 regrid("PDD surface model", m_snow_depth);
207 regrid("PDD surface model", m_firn_depth);
208 }
209
210 // finish up
211 {
212 m_next_balance_year_start = compute_next_balance_year_start(m_grid->ctx()->time()->current());
213
214 m_accumulation->set(0.0);
215 m_melt->set(0.0);
216 m_runoff->set(0.0);
217 }
218 }
219
220 MaxTimestep TemperatureIndex::max_timestep_impl(double my_t) const {
221 return m_atmosphere->max_timestep(my_t);
222 }
223
224 double TemperatureIndex::compute_next_balance_year_start(double time) {
225 // compute the time corresponding to the beginning of the next balance year
226 double
227 balance_year_start_day = m_config->get_number("surface.pdd.balance_year_start_day"),
228 one_day = units::convert(m_sys, 1.0, "days", "seconds"),
229 year_start = m_grid->ctx()->time()->calendar_year_start(time),
230 balance_year_start = year_start + (balance_year_start_day - 1.0) * one_day;
231
232 if (balance_year_start > time) {
233 return balance_year_start;
234 }
235 return m_grid->ctx()->time()->increment_date(balance_year_start, 1);
236 }
237
238 void TemperatureIndex::update_impl(const Geometry &geometry, double t, double dt) {
239
240 // make a copy of the pointer to convince clang static analyzer that its value does not
241 // change during the call
242 FaustoGrevePDDObject *fausto_greve = m_faustogreve.get();
243
244 // update to ensure that temperature and precipitation time series are correct:
245 m_atmosphere->update(geometry, t, dt);
246
247 m_temperature->copy_from(m_atmosphere->mean_annual_temp());
248
249 // set up air temperature and precipitation time series
250 int N = m_mbscheme->get_timeseries_length(dt);
251
252 const double dtseries = dt / N;
253 std::vector<double> ts(N), T(N), S(N), P(N), PDDs(N);
254 for (int k = 0; k < N; ++k) {
255 ts[k] = t + k * dtseries;
256 }
257
258 // update standard deviation time series
259 if (m_sd_file_set) {
260 m_air_temp_sd->update(t, dt);
261 m_air_temp_sd->init_interpolation(ts);
262 }
263
264 const IceModelVec2CellType &mask = geometry.cell_type;
265 const IceModelVec2S &H = geometry.ice_thickness;
266
267 IceModelVec::AccessList list{&mask, &H, m_air_temp_sd.get(), &m_mass_flux,
268 &m_firn_depth, &m_snow_depth,
269 m_accumulation.get(), m_melt.get(), m_runoff.get()};
270
271 const double
272 sigmalapserate = m_config->get_number("surface.pdd.std_dev_lapse_lat_rate"),
273 sigmabaselat = m_config->get_number("surface.pdd.std_dev_lapse_lat_base");
274
275 const IceModelVec2S *latitude = nullptr;
276 if (fausto_greve or sigmalapserate != 0.0) {
277 latitude = &geometry.latitude;
278
279 list.add(*latitude);
280 }
281
282 if (fausto_greve) {
283 const IceModelVec2S
284 *longitude = &geometry.latitude,
285 *surface_altitude = &geometry.ice_surface_elevation;
286
287 fausto_greve->update_temp_mj(*surface_altitude, *latitude, *longitude);
288 }
289
290 LocalMassBalance::DegreeDayFactors ddf = m_base_ddf;
291
292 m_atmosphere->init_timeseries(ts);
293
294 m_atmosphere->begin_pointwise_access();
295
296 const double ice_density = m_config->get_number("constants.ice.density");
297
298 ParallelSection loop(m_grid->com);
299 try {
300 for (Points p(*m_grid); p; p.next()) {
301 const int i = p.i(), j = p.j();
302
303 // the temperature time series from the AtmosphereModel and its modifiers
304 m_atmosphere->temp_time_series(i, j, T);
305
306 if (mask.ice_free_ocean(i, j)) {
307 // ignore precipitation over ice-free ocean
308 for (int k = 0; k < N; ++k) {
309 P[k] = 0.0;
310 }
311 } else {
312 // elsewhere, get precipitation from the atmosphere model
313 m_atmosphere->precip_time_series(i, j, P);
314 }
315
316 // convert precipitation from "kg m-2 second-1" to "m second-1" (PDDMassBalance expects
317 // accumulation in m/second ice equivalent)
318 for (int k = 0; k < N; ++k) {
319 P[k] = P[k] / ice_density;
320 // kg / (m^2 * second) / (kg / m^3) = m / second
321 }
322
323 // interpolate temperature standard deviation time series
324 if (m_sd_file_set) {
325 m_air_temp_sd->interp(i, j, S);
326 } else {
327 double tmp = (*m_air_temp_sd)(i, j);
328 for (int k = 0; k < N; ++k) {
329 S[k] = tmp;
330 }
331 }
332
333 if (fausto_greve) {
334 // we have been asked to set mass balance parameters according to
335 // formula (6) in [\ref Faustoetal2009]; they overwrite ddf set above
336 ddf = fausto_greve->degree_day_factors(i, j, (*latitude)(i, j));
337 }
338
339 // apply standard deviation lapse rate on top of prescribed values
340 if (sigmalapserate != 0.0) {
341 double lat = (*latitude)(i, j);
342 for (int k = 0; k < N; ++k) {
343 S[k] += sigmalapserate * (lat - sigmabaselat);
344 }
345 (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
346 }
347
348 // apply standard deviation param over ice if in use
349 if (m_sd_use_param and mask.icy(i, j)) {
350 for (int k = 0; k < N; ++k) {
351 S[k] = m_sd_param_a * (T[k] - 273.15) + m_sd_param_b;
352 if (S[k] < 0.0) {
353 S[k] = 0.0 ;
354 }
355 }
356 (*m_air_temp_sd)(i, j) = S[0]; // ensure correct SD reporting
357 }
358
359 // Use temperature time series, the "positive" threshhold, and
360 // the standard deviation of the daily variability to get the
361 // number of positive degree days (PDDs)
362 if (mask.ice_free_ocean(i, j)) {
363 for (int k = 0; k < N; ++k) {
364 PDDs[k] = 0.0;
365 }
366 } else {
367 m_mbscheme->get_PDDs(dtseries, S, T, // inputs
368 PDDs); // output
369 }
370
371 // Use temperature time series to remove rainfall from precipitation
372 m_mbscheme->get_snow_accumulation(T, // air temperature (input)
373 P); // precipitation rate (input-output)
374
375 // Use degree-day factors, the number of PDDs, and the snow precipitation to get surface mass
376 // balance (and diagnostics: accumulation, melt, runoff)
377 {
378 double next_snow_depth_reset = m_next_balance_year_start;
379
380 // make copies of firn and snow depth values at this point to avoid accessing 2D
381 // fields in the inner loop
382 double
383 ice = H(i, j),
384 firn = m_firn_depth(i, j),
385 snow = m_snow_depth(i, j);
386
387 // accumulation, melt, runoff over this time-step
388 double
389 A = 0.0,
390 M = 0.0,
391 R = 0.0,
392 SMB = 0.0;
393
394 for (int k = 0; k < N; ++k) {
395 if (ts[k] >= next_snow_depth_reset) {
396 snow = 0.0;
397 while (next_snow_depth_reset <= ts[k]) {
398 next_snow_depth_reset = m_grid->ctx()->time()->increment_date(next_snow_depth_reset, 1);
399 }
400 }
401
402 const double accumulation = P[k] * dtseries;
403
404 LocalMassBalance::Changes changes;
405 changes = m_mbscheme->step(ddf, PDDs[k],
406 ice, firn, snow, accumulation);
407
408 // update ice thickness
409 ice += changes.smb;
410 assert(ice >= 0);
411
412 // update firn depth
413 firn += changes.firn_depth;
414 assert(firn >= 0);
415
416 // update snow depth
417 snow += changes.snow_depth;
418 assert(snow >= 0);
419
420 // update total accumulation, melt, and runoff
421 {
422 A += accumulation;
423 M += changes.melt;
424 R += changes.runoff;
425 SMB += changes.smb;
426 }
427 } // end of the time-stepping loop
428
429 // set firn and snow depths
430 m_firn_depth(i, j) = firn;
431 m_snow_depth(i, j) = snow;
432
433 // set total accumulation, melt, and runoff, and SMB at this point, converting
434 // from "meters, ice equivalent" to "kg / m^2"
435 {
436 (*m_accumulation)(i, j) = A * ice_density;
437 (*m_melt)(i, j) = M * ice_density;
438 (*m_runoff)(i, j) = R * ice_density;
439 // m_mass_flux (unlike m_accumulation, m_melt, and m_runoff), is a
440 // rate. m * (kg / m^3) / second = kg / m^2 / second
441 m_mass_flux(i, j) = SMB * ice_density / dt;
442 }
443 }
444
445 if (mask.ice_free_ocean(i, j)) {
446 m_firn_depth(i, j) = 0.0; // no firn in the ocean
447 m_snow_depth(i, j) = 0.0; // snow over the ocean does not stick
448 }
449 }
450 } catch (...) {
451 loop.failed();
452 }
453 loop.check();
454
455 m_atmosphere->end_pointwise_access();
456
457 m_next_balance_year_start = compute_next_balance_year_start(m_grid->ctx()->time()->current());
458 }
459
460 const IceModelVec2S &TemperatureIndex::mass_flux_impl() const {
461 return m_mass_flux;
462 }
463
464 const IceModelVec2S &TemperatureIndex::temperature_impl() const {
465 return *m_temperature;
466 }
467
468 const IceModelVec2S& TemperatureIndex::accumulation_impl() const {
469 return *m_accumulation;
470 }
471
472 const IceModelVec2S& TemperatureIndex::melt_impl() const {
473 return *m_melt;
474 }
475
476 const IceModelVec2S& TemperatureIndex::runoff_impl() const {
477 return *m_runoff;
478 }
479
480 const IceModelVec2S& TemperatureIndex::firn_depth() const {
481 return m_firn_depth;
482 }
483
484 const IceModelVec2S& TemperatureIndex::snow_depth() const {
485 return m_snow_depth;
486 }
487
488 const IceModelVec2S& TemperatureIndex::air_temp_sd() const {
489 return *m_air_temp_sd;
490 }
491
492 void TemperatureIndex::define_model_state_impl(const File &output) const {
493 SurfaceModel::define_model_state_impl(output);
494 m_firn_depth.define(output, PISM_DOUBLE);
495 m_snow_depth.define(output, PISM_DOUBLE);
496 }
497
498 void TemperatureIndex::write_model_state_impl(const File &output) const {
499 SurfaceModel::write_model_state_impl(output);
500 m_firn_depth.write(output);
501 m_snow_depth.write(output);
502 }
503
504 namespace diagnostics {
505
506 enum AmountKind {AMOUNT, MASS};
507
508 /*! @brief Report surface melt, averaged over the reporting interval */
509 class SurfaceMelt : public DiagAverageRate<TemperatureIndex>
510 {
511 public:
512 SurfaceMelt(const TemperatureIndex *m, AmountKind kind)
513 : DiagAverageRate<TemperatureIndex>(m,
514 kind == AMOUNT
515 ? "surface_melt_flux"
516 : "surface_melt_rate",
517 TOTAL_CHANGE),
518 m_kind(kind) {
519
520 std::string
521 name = "surface_melt_flux",
522 long_name = "surface melt, averaged over the reporting interval",
523 standard_name = "surface_snow_and_ice_melt_flux",
524 accumulator_units = "kg m-2",
525 internal_units = "kg m-2 second-1",
526 external_units = "kg m-2 year-1";
527 if (kind == MASS) {
528 name = "surface_melt_rate";
529 standard_name = "";
530 accumulator_units = "kg",
531 internal_units = "kg second-1";
532 external_units = "Gt year-1" ;
533
534 m_melt_mass.create(m_grid, "melt_mass", WITHOUT_GHOSTS);
535 }
536
537 m_vars = {SpatialVariableMetadata(m_sys, name)};
538 m_accumulator.metadata().set_string("units", accumulator_units);
539
540 set_attrs(long_name, standard_name, internal_units, external_units, 0);
541 m_vars[0].set_string("cell_methods", "time: mean");
542
543 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
544 }
545
546 protected:
547 const IceModelVec2S& model_input() {
548 const IceModelVec2S &melt_amount = model->melt();
549
550 if (m_kind == MASS) {
551 double cell_area = m_grid->cell_area();
552
553 IceModelVec::AccessList list{&m_melt_mass, &melt_amount};
554
555 for (Points p(*m_grid); p; p.next()) {
556 const int i = p.i(), j = p.j();
557 m_melt_mass(i, j) = melt_amount(i, j) * cell_area;
558 }
559 return m_melt_mass;
560 } else {
561 return melt_amount;
562 }
563 }
564 private:
565 IceModelVec2S m_melt_mass;
566 AmountKind m_kind;
567 };
568
569 /*! @brief Report surface runoff, averaged over the reporting interval */
570 class SurfaceRunoff : public DiagAverageRate<TemperatureIndex>
571 {
572 public:
573 SurfaceRunoff(const TemperatureIndex *m, AmountKind kind)
574 : DiagAverageRate<TemperatureIndex>(m,
575 kind == AMOUNT
576 ? "surface_runoff_flux"
577 : "surface_runoff_rate",
578 TOTAL_CHANGE),
579 m_kind(kind) {
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 m_runoff_mass.create(m_grid, "runoff_mass", WITHOUT_GHOSTS);
596 }
597
598 m_vars = {SpatialVariableMetadata(m_sys, name)};
599 m_accumulator.metadata().set_string("units", accumulator_units);
600
601 set_attrs(long_name, standard_name, internal_units, external_units, 0);
602 m_vars[0].set_string("cell_methods", "time: mean");
603
604 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
605 }
606
607 protected:
608 const IceModelVec2S& model_input() {
609 const IceModelVec2S &runoff_amount = model->runoff();
610
611 if (m_kind == MASS) {
612 double cell_area = m_grid->cell_area();
613
614 IceModelVec::AccessList list{&m_runoff_mass, &runoff_amount};
615
616 for (Points p(*m_grid); p; p.next()) {
617 const int i = p.i(), j = p.j();
618 m_runoff_mass(i, j) = runoff_amount(i, j) * cell_area;
619 }
620 return m_runoff_mass;
621 } else {
622 return runoff_amount;
623 }
624 }
625 private:
626 AmountKind m_kind;
627 IceModelVec2S m_runoff_mass;
628 };
629
630 /*! @brief Report accumulation (precipitation minus rain), averaged over the reporting interval */
631 class Accumulation : public DiagAverageRate<TemperatureIndex>
632 {
633 public:
634 Accumulation(const TemperatureIndex *m, AmountKind kind)
635 : DiagAverageRate<TemperatureIndex>(m,
636 kind == AMOUNT
637 ? "surface_accumulation_flux"
638 : "surface_accumulation_rate",
639 TOTAL_CHANGE),
640 m_kind(kind) {
641
642 // possible standard name: surface_accumulation_flux
643 std::string
644 name = "surface_accumulation_flux",
645 long_name = "accumulation (precipitation minus rain), averaged over the reporting interval",
646 accumulator_units = "kg m-2",
647 internal_units = "kg m-2 second-1",
648 external_units = "kg m-2 year-1";
649 if (kind == MASS) {
650 name = "surface_accumulation_rate";
651 accumulator_units = "kg",
652 internal_units = "kg second-1";
653 external_units = "Gt year-1" ;
654
655 m_accumulation_mass.create(m_grid, "accumulation_mass", WITHOUT_GHOSTS);
656 }
657
658
659 m_vars = {SpatialVariableMetadata(m_sys, name)};
660 m_accumulator.metadata().set_string("units", accumulator_units);
661
662 set_attrs(long_name, "", internal_units, external_units, 0);
663 m_vars[0].set_string("cell_methods", "time: mean");
664
665 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
666 }
667
668 protected:
669 const IceModelVec2S& model_input() {
670 const IceModelVec2S &accumulation_amount = model->accumulation();
671
672 if (m_kind == MASS) {
673 double cell_area = m_grid->cell_area();
674
675 IceModelVec::AccessList list{&m_accumulation_mass, &accumulation_amount};
676
677 for (Points p(*m_grid); p; p.next()) {
678 const int i = p.i(), j = p.j();
679 m_accumulation_mass(i, j) = accumulation_amount(i, j) * cell_area;
680 }
681 return m_accumulation_mass;
682 } else {
683 return accumulation_amount;
684 }
685 }
686 private:
687 AmountKind m_kind;
688 IceModelVec2S m_accumulation_mass;
689 };
690
691 /*!
692 * Integrate a field over the computational domain.
693 *
694 * If the input has units kg/m^2, the output will be in kg.
695 */
696 static double integrate(const IceModelVec2S &input) {
697 IceGrid::ConstPtr grid = input.grid();
698
699 double cell_area = grid->cell_area();
700
701 IceModelVec::AccessList list{&input};
702
703 double result = 0.0;
704
705 for (Points p(*grid); p; p.next()) {
706 const int i = p.i(), j = p.j();
707
708 result += input(i, j) * cell_area;
709 }
710
711 return GlobalSum(grid->com, result);
712 }
713
714
715 //! \brief Reports the total accumulation rate.
716 class TotalSurfaceAccumulation : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
717 {
718 public:
719 TotalSurfaceAccumulation(const TemperatureIndex *m)
720 : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_accumulation_rate") {
721
722 set_units("kg s-1", "kg year-1");
723 m_ts.variable().set_string("long_name", "surface accumulation rate (PDD model)");
724 }
725
726 double compute() {
727 return integrate(model->accumulation());
728 }
729 };
730
731
732 //! \brief Reports the total melt rate.
733 class TotalSurfaceMelt : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
734 {
735 public:
736 TotalSurfaceMelt(const TemperatureIndex *m)
737 : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_melt_rate") {
738
739 set_units("kg s-1", "kg year-1");
740 m_ts.variable().set_string("long_name", "surface melt rate (PDD model)");
741 }
742
743 double compute() {
744 return integrate(model->melt());
745 }
746 };
747
748
749 //! \brief Reports the total top surface ice flux.
750 class TotalSurfaceRunoff : public TSDiag<TSFluxDiagnostic, TemperatureIndex>
751 {
752 public:
753 TotalSurfaceRunoff(const TemperatureIndex *m)
754 : TSDiag<TSFluxDiagnostic, TemperatureIndex>(m, "surface_runoff_rate") {
755
756 set_units("kg s-1", "kg year-1");
757 m_ts.variable().set_string("long_name", "surface runoff rate (PDD model)");
758 }
759
760 double compute() {
761 return integrate(model->runoff());
762 }
763 };
764
765 } // end of namespace diagnostics
766
767 DiagnosticList TemperatureIndex::diagnostics_impl() const {
768 using namespace diagnostics;
769
770 DiagnosticList result = {
771 {"surface_accumulation_flux", Diagnostic::Ptr(new Accumulation(this, AMOUNT))},
772 {"surface_accumulation_rate", Diagnostic::Ptr(new Accumulation(this, MASS))},
773 {"surface_melt_flux", Diagnostic::Ptr(new SurfaceMelt(this, AMOUNT))},
774 {"surface_melt_rate", Diagnostic::Ptr(new SurfaceMelt(this, MASS))},
775 {"surface_runoff_flux", Diagnostic::Ptr(new SurfaceRunoff(this, AMOUNT))},
776 {"surface_runoff_rate", Diagnostic::Ptr(new SurfaceRunoff(this, MASS))},
777 {"air_temp_sd", Diagnostic::wrap(*m_air_temp_sd)},
778 {"snow_depth", Diagnostic::wrap(m_snow_depth)},
779 {"firn_depth", Diagnostic::wrap(m_firn_depth)},
780 };
781
782 result = pism::combine(result, SurfaceModel::diagnostics_impl());
783
784 return result;
785 }
786
787 TSDiagnosticList TemperatureIndex::ts_diagnostics_impl() const {
788 using namespace diagnostics;
789
790 TSDiagnosticList result = {
791 {"surface_accumulation_rate", TSDiagnostic::Ptr(new TotalSurfaceAccumulation(this))},
792 {"surface_melt_rate", TSDiagnostic::Ptr(new TotalSurfaceMelt(this))},
793 {"surface_runoff_rate", TSDiagnostic::Ptr(new TotalSurfaceRunoff(this))},
794 };
795
796 result = pism::combine(result, SurfaceModel::ts_diagnostics_impl());
797
798 return result;
799 }
800
801 } // end of namespace surface
802 } // end of namespace pism