tio_helpers.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
---
tio_helpers.cc (55138B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019, 2020 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
20 #include <memory>
21 #include <cassert>
22
23 #include "io_helpers.hh"
24 #include "File.hh"
25 #include "pism/util/IceGrid.hh"
26 #include "pism/util/VariableMetadata.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/ConfigInterface.hh"
30 #include "pism/util/io/LocalInterpCtx.hh"
31 #include "pism/util/Time.hh"
32 #include "pism/util/Logger.hh"
33 #include "pism/util/Context.hh"
34 #include "pism/util/projection.hh"
35 #include "pism/util/interpolation.hh"
36 #include "pism/util/Profiling.hh"
37
38 namespace pism {
39 namespace io {
40
41 //! \brief Bi-(or tri-)linear interpolation.
42 /*!
43 * This is the interpolation code itself.
44 *
45 * Note that its inputs are (essentially)
46 * - the definition of the input grid
47 * - the definition of the output grid
48 * - input array (lic->buffer)
49 * - output array (double *output_array)
50 *
51 * The `output_array` is expected to be big enough to contain
52 * `grid.xm()*`grid.ym()*length(zlevels_out)` numbers.
53 *
54 * We should be able to switch to using an external interpolation library
55 * fairly easily...
56 */
57 static void regrid(const IceGrid& grid, const std::vector<double> &zlevels_out,
58 LocalInterpCtx *lic, double *output_array) {
59 // We'll work with the raw storage here so that the array we are filling is
60 // indexed the same way as the buffer we are pulling from (input_array)
61
62 const int X = 1, Z = 3; // indices, just for clarity
63
64 unsigned int nlevels = zlevels_out.size();
65 double *input_array = &(lic->buffer[0]);
66
67 // array sizes for mapping from logical to "flat" indices
68 int
69 x_count = lic->count[X],
70 z_count = lic->count[Z];
71
72 for (Points p(grid); p; p.next()) {
73 const int i_global = p.i(), j_global = p.j();
74
75 const int i = i_global - grid.xs(), j = j_global - grid.ys();
76
77 // Indices of neighboring points.
78 const int
79 X_m = lic->x->left(i),
80 X_p = lic->x->right(i),
81 Y_m = lic->y->left(j),
82 Y_p = lic->y->right(j);
83
84 for (unsigned int k = 0; k < nlevels; k++) {
85
86 double
87 a_mm = 0.0,
88 a_mp = 0.0,
89 a_pm = 0.0,
90 a_pp = 0.0;
91
92 if (nlevels > 1) {
93 const int
94 Z_m = lic->z->left(k),
95 Z_p = lic->z->right(k);
96
97 const double alpha_z = lic->z->alpha(k);
98
99 // We pretend that there are always 8 neighbors (4 in the map plane,
100 // 2 vertical levels). And compute the indices into the input_array for
101 // those neighbors.
102 const int
103 mmm = (Y_m * x_count + X_m) * z_count + Z_m,
104 mmp = (Y_m * x_count + X_m) * z_count + Z_p,
105 mpm = (Y_m * x_count + X_p) * z_count + Z_m,
106 mpp = (Y_m * x_count + X_p) * z_count + Z_p,
107 pmm = (Y_p * x_count + X_m) * z_count + Z_m,
108 pmp = (Y_p * x_count + X_m) * z_count + Z_p,
109 ppm = (Y_p * x_count + X_p) * z_count + Z_m,
110 ppp = (Y_p * x_count + X_p) * z_count + Z_p;
111
112 // linear interpolation in the z-direction
113 a_mm = input_array[mmm] * (1.0 - alpha_z) + input_array[mmp] * alpha_z;
114 a_mp = input_array[mpm] * (1.0 - alpha_z) + input_array[mpp] * alpha_z;
115 a_pm = input_array[pmm] * (1.0 - alpha_z) + input_array[pmp] * alpha_z;
116 a_pp = input_array[ppm] * (1.0 - alpha_z) + input_array[ppp] * alpha_z;
117 } else {
118 // we don't need to interpolate vertically for the 2-D case
119 a_mm = input_array[Y_m * x_count + X_m];
120 a_mp = input_array[Y_m * x_count + X_p];
121 a_pm = input_array[Y_p * x_count + X_m];
122 a_pp = input_array[Y_p * x_count + X_p];
123 }
124
125 // interpolation coefficient in the x direction
126 const double x_alpha = lic->x->alpha(i);
127 // interpolation coefficient in the y direction
128 const double y_alpha = lic->y->alpha(j);
129
130 // interpolate in x direction
131 const double a_m = a_mm * (1.0 - x_alpha) + a_mp * x_alpha;
132 const double a_p = a_pm * (1.0 - x_alpha) + a_pp * x_alpha;
133
134 int index = (j * grid.xm() + i) * nlevels + k;
135
136 // index into the new array and interpolate in x direction
137 output_array[index] = a_m * (1.0 - y_alpha) + a_p * y_alpha;
138 // done with the point at (x,y,z)
139 }
140 }
141 }
142
143 static void compute_start_and_count(const File& file,
144 units::System::Ptr unit_system,
145 const std::string &short_name,
146 unsigned int t_start, unsigned int t_count,
147 unsigned int x_start, unsigned int x_count,
148 unsigned int y_start, unsigned int y_count,
149 unsigned int z_start, unsigned int z_count,
150 std::vector<unsigned int> &start,
151 std::vector<unsigned int> &count,
152 std::vector<unsigned int> &imap) {
153 std::vector<std::string> dims = file.dimensions(short_name);
154 unsigned int ndims = dims.size();
155
156 assert(ndims > 0);
157 if (ndims == 0) {
158 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
159 "Cannot compute start and count: variable %s is a scalar.",
160 short_name.c_str());
161 }
162
163 // Resize output vectors:
164 start.resize(ndims);
165 count.resize(ndims);
166 imap.resize(ndims);
167
168 // Assemble start, count and imap:
169 for (unsigned int j = 0; j < ndims; j++) {
170 std::string dimname = dims[j];
171
172 AxisType dimtype = file.dimension_type(dimname, unit_system);
173
174 switch (dimtype) {
175 case T_AXIS:
176 start[j] = t_start;
177 count[j] = t_count;
178 imap[j] = x_count * y_count * z_count;
179 break;
180 case Y_AXIS:
181 start[j] = y_start;
182 count[j] = y_count;
183 imap[j] = x_count * z_count;
184 break;
185 case X_AXIS:
186 start[j] = x_start;
187 count[j] = x_count;
188 imap[j] = z_count;
189 break;
190 default:
191 case Z_AXIS:
192 start[j] = z_start;
193 count[j] = z_count;
194 imap[j] = 1;
195 break;
196 }
197 }
198 }
199
200 //! \brief Define a dimension \b and the associated coordinate variable. Set attributes.
201 void define_dimension(const File &file, unsigned long int length,
202 const VariableMetadata &metadata) {
203 std::string name = metadata.get_name();
204 try {
205 file.define_dimension(name, length);
206
207 std::vector<std::string> dims(1, name);
208 file.define_variable(name, PISM_DOUBLE, dims);
209
210 write_attributes(file, metadata, PISM_DOUBLE);
211
212 } catch (RuntimeError &e) {
213 e.add_context("defining dimension '%s' in '%s'", name.c_str(), file.filename().c_str());
214 throw;
215 }
216 }
217
218
219 //! Prepare a file for output.
220 void define_time(const File &file, const Context &ctx) {
221 const Time &time = *ctx.time();
222 const Config &config = *ctx.config();
223
224 define_time(file,
225 config.get_string("time.dimension_name"),
226 time.calendar(),
227 time.CF_units_string(),
228 ctx.unit_system());
229 }
230
231 /*!
232 * Define a time dimension and the corresponding coordinate variable. Does nothing if the time
233 * variable is already present.
234 */
235 void define_time(const File &file, const std::string &name, const std::string &calendar,
236 const std::string &units, units::System::Ptr unit_system) {
237 try {
238 if (file.find_variable(name)) {
239 return;
240 }
241
242 // time
243 VariableMetadata time(name, unit_system);
244 time.set_string("long_name", "time");
245 time.set_string("calendar", calendar);
246 time.set_string("units", units);
247 time.set_string("axis", "T");
248
249 define_dimension(file, PISM_UNLIMITED, time);
250 } catch (RuntimeError &e) {
251 e.add_context("defining the time dimension in \"" + file.filename() + "\"");
252 throw;
253 }
254 }
255
256 //! Prepare a file for output.
257 void append_time(const File &file, const Config &config, double time_seconds) {
258 append_time(file, config.get_string("time.dimension_name"),
259 time_seconds);
260 }
261
262 //! \brief Append to the time dimension.
263 void append_time(const File &file, const std::string &name, double value) {
264 try {
265 unsigned int start = file.dimension_length(name);
266
267 file.write_variable(name, {start}, {1}, &value);
268
269 // PIO's I/O type PnetCDF requires this
270 file.sync();
271 } catch (RuntimeError &e) {
272 e.add_context("appending to the time dimension in \"" + file.filename() + "\"");
273 throw;
274 }
275 }
276
277 //! \brief Define dimensions a variable depends on.
278 static void define_dimensions(const SpatialVariableMetadata& var,
279 const IceGrid& grid, const File &file) {
280
281 // x
282 std::string x_name = var.get_x().get_name();
283 if (not file.find_dimension(x_name)) {
284 define_dimension(file, grid.Mx(), var.get_x());
285 file.write_attribute(x_name, "spacing_meters", PISM_DOUBLE,
286 {grid.x(1) - grid.x(0)});
287 file.write_attribute(x_name, "not_written", PISM_INT, {1.0});
288 }
289
290 // y
291 std::string y_name = var.get_y().get_name();
292 if (not file.find_dimension(y_name)) {
293 define_dimension(file, grid.My(), var.get_y());
294 file.write_attribute(y_name, "spacing_meters", PISM_DOUBLE,
295 {grid.y(1) - grid.y(0)});
296 file.write_attribute(y_name, "not_written", PISM_INT, {1.0});
297 }
298
299 // z
300 std::string z_name = var.get_z().get_name();
301 if (not z_name.empty()) {
302 if (not file.find_dimension(z_name)) {
303 const std::vector<double>& levels = var.get_levels();
304 // make sure we have at least one level
305 unsigned int nlevels = std::max(levels.size(), (size_t)1);
306 define_dimension(file, nlevels, var.get_z());
307 file.write_attribute(z_name, "not_written", PISM_INT, {1.0});
308
309 bool spatial_dim = not var.get_z().get_string("axis").empty();
310
311 if (nlevels > 1 and spatial_dim) {
312 double dz_max = levels[1] - levels[0];
313 double dz_min = levels.back() - levels.front();
314
315 for (unsigned int k = 0; k < nlevels - 1; ++k) {
316 double dz = levels[k+1] - levels[k];
317 dz_max = std::max(dz_max, dz);
318 dz_min = std::min(dz_min, dz);
319 }
320
321 file.write_attribute(z_name, "spacing_min_meters", PISM_DOUBLE,
322 {dz_min});
323 file.write_attribute(z_name, "spacing_max_meters", PISM_DOUBLE,
324 {dz_max});
325 }
326 }
327 }
328 }
329
330 static void write_dimension_data(const File &file, const std::string &name,
331 const std::vector<double> &data) {
332 bool written = file.attribute_type(name, "not_written") == PISM_NAT;
333 if (not written) {
334 file.write_variable(name, {0}, {(unsigned int)data.size()}, data.data());
335 file.redef();
336 file.remove_attribute(name, "not_written");
337 }
338 }
339
340 void write_dimensions(const SpatialVariableMetadata& var,
341 const IceGrid& grid, const File &file) {
342 // x
343 std::string x_name = var.get_x().get_name();
344 if (file.find_dimension(x_name)) {
345 write_dimension_data(file, x_name, grid.x());
346 }
347
348 // y
349 std::string y_name = var.get_y().get_name();
350 if (file.find_dimension(y_name)) {
351 write_dimension_data(file, y_name, grid.y());
352 }
353
354 // z
355 std::string z_name = var.get_z().get_name();
356 if (file.find_dimension(z_name)) {
357 write_dimension_data(file, z_name, var.get_levels());
358 }
359 }
360
361 /**
362 * Check if the storage order of a variable in the current file
363 * matches the memory storage order used by PISM.
364 *
365 * @param[in] file input file
366 * @param var_name name of the variable to check
367 * @returns false if storage orders match, true otherwise
368 */
369 static bool use_transposed_io(const File &file,
370 units::System::Ptr unit_system,
371 const std::string &var_name) {
372
373 std::vector<std::string> dimnames = file.dimensions(var_name);
374
375 std::vector<AxisType> storage, memory = {Y_AXIS, X_AXIS};
376
377 for (unsigned int j = 0; j < dimnames.size(); ++j) {
378 AxisType dimtype = file.dimension_type(dimnames[j], unit_system);
379
380 if (j == 0 && dimtype == T_AXIS) {
381 // ignore the time dimension, but only if it is the first
382 // dimension in the list
383 continue;
384 }
385
386 if (dimtype == X_AXIS || dimtype == Y_AXIS) {
387 storage.push_back(dimtype);
388 } else if (dimtype == Z_AXIS) {
389 memory.push_back(dimtype); // now memory = {Y_AXIS, X_AXIS, Z_AXIS}
390 // assume that this variable has only one Z_AXIS in the file
391 storage.push_back(dimtype);
392 } else {
393 // an UNKNOWN_AXIS or T_AXIS at index != 0 was found, use mapped I/O
394 return true;
395 }
396 }
397
398 // we support 2D and 3D in-memory arrays, but not 4D
399 assert(memory.size() <= 3);
400
401 if (storage == memory) {
402 // same storage order, do not use mapped I/O
403 return false;
404 } else {
405 // different storage orders, use mapped I/O
406 return true;
407 }
408 }
409
410 //! \brief Read an array distributed according to the grid.
411 static void read_distributed_array(const File &file, const IceGrid &grid,
412 const std::string &var_name,
413 unsigned int z_count, unsigned int t_start,
414 double *output) {
415 try {
416 std::vector<unsigned int> start, count, imap;
417 const unsigned int t_count = 1;
418 compute_start_and_count(file,
419 grid.ctx()->unit_system(),
420 var_name,
421 t_start, t_count,
422 grid.xs(), grid.xm(),
423 grid.ys(), grid.ym(),
424 0, z_count,
425 start, count, imap);
426
427 bool transposed_io = use_transposed_io(file, grid.ctx()->unit_system(), var_name);
428 if (transposed_io) {
429 file.read_variable_transposed(var_name, start, count, imap, output);
430 } else {
431 file.read_variable(var_name, start, count, output);
432 }
433
434 } catch (RuntimeError &e) {
435 e.add_context("reading variable '%s' from '%s'", var_name.c_str(),
436 file.filename().c_str());
437 throw;
438 }
439 }
440
441 static void regrid_vec_generic(const File &file, const IceGrid &grid,
442 const std::string &variable_name,
443 const std::vector<double> &zlevels_out,
444 unsigned int t_start,
445 bool fill_missing,
446 double default_value,
447 InterpolationType interpolation_type,
448 double *output) {
449 const int X = 1, Y = 2, Z = 3; // indices, just for clarity
450
451 const Profiling& profiling = grid.ctx()->profiling();
452
453 try {
454 grid_info gi(file, variable_name, grid.ctx()->unit_system(), grid.registration());
455 LocalInterpCtx lic(gi, grid, zlevels_out, interpolation_type);
456
457 std::vector<double> &buffer = lic.buffer;
458
459 const unsigned int t_count = 1;
460 std::vector<unsigned int> start, count, imap;
461 compute_start_and_count(file,
462 grid.ctx()->unit_system(),
463 variable_name,
464 t_start, t_count,
465 lic.start[X], lic.count[X],
466 lic.start[Y], lic.count[Y],
467 lic.start[Z], lic.count[Z],
468 start, count, imap);
469
470 bool transposed_io = use_transposed_io(file, grid.ctx()->unit_system(), variable_name);
471 profiling.begin("io.regridding.read");
472 if (transposed_io) {
473 file.read_variable_transposed(variable_name, start, count, imap, &buffer[0]);
474 } else {
475 file.read_variable(variable_name, start, count, &buffer[0]);
476 }
477 profiling.end("io.regridding.read");
478
479 // Replace missing values if the _FillValue attribute is present,
480 // and if we have missing values to replace.
481 if (fill_missing) {
482 std::vector<double> attribute = file.read_double_attribute(variable_name, "_FillValue");
483 if (attribute.size() == 1) {
484 const double fill_value = attribute[0],
485 epsilon = 1e-12;
486 for (unsigned int i = 0; i < buffer.size(); ++i) {
487 if (fabs(buffer[i] - fill_value) < epsilon) {
488 buffer[i] = default_value;
489 }
490 }
491 }
492 }
493
494 // interpolate
495 profiling.begin("io.regridding.interpolate");
496 regrid(grid, zlevels_out, &lic, output);
497 profiling.end("io.regridding.interpolate");
498 } catch (RuntimeError &e) {
499 e.add_context("reading variable '%s' (using linear interpolation) from '%s'",
500 variable_name.c_str(), file.filename().c_str());
501 throw;
502 }
503 }
504
505 //! \brief Read a PETSc Vec from a file, using bilinear (or trilinear)
506 //! interpolation to put it on the grid defined by "grid" and zlevels_out.
507 static void regrid_vec(const File &file, const IceGrid &grid, const std::string &var_name,
508 const std::vector<double> &zlevels_out,
509 unsigned int t_start,
510 InterpolationType interpolation_type,
511 double *output) {
512 regrid_vec_generic(file, grid,
513 var_name,
514 zlevels_out,
515 t_start,
516 false, 0.0,
517 interpolation_type,
518 output);
519 }
520
521 /** Regrid `var_name` from a file, replacing missing values with `default_value`.
522 *
523 * @param[in] file input file
524 * @param grid computational grid; used to initialize interpolation
525 * @param var_name variable to regrid
526 * @param zlevels_out vertical levels of the resulting grid
527 * @param t_start time index of the record to regrid
528 * @param default_value default value to replace `_FillValue` with
529 * @param[out] output resulting interpolated field
530 */
531 static void regrid_vec_fill_missing(const File &file, const IceGrid &grid,
532 const std::string &var_name,
533 const std::vector<double> &zlevels_out,
534 unsigned int t_start,
535 double default_value,
536 InterpolationType interpolation_type,
537 double *output) {
538 regrid_vec_generic(file, grid,
539 var_name,
540 zlevels_out,
541 t_start,
542 true, default_value,
543 interpolation_type,
544 output);
545 }
546
547 //! Define a NetCDF variable corresponding to a VariableMetadata object.
548 void define_spatial_variable(const SpatialVariableMetadata &var,
549 const IceGrid &grid, const File &file,
550 IO_Type default_type) {
551 std::vector<std::string> dims;
552 std::string name = var.get_name();
553
554 if (file.find_variable(name)) {
555 return;
556 }
557
558 define_dimensions(var, grid, file);
559
560 std::string
561 x = var.get_x().get_name(),
562 y = var.get_y().get_name(),
563 z = var.get_z().get_name();
564
565 if (not var.get_time_independent()) {
566 dims.push_back(grid.ctx()->config()->get_string("time.dimension_name"));
567 }
568
569 dims.push_back(y);
570 dims.push_back(x);
571
572 if (not z.empty()) {
573 dims.push_back(z);
574 }
575
576 assert(dims.size() > 1);
577
578 IO_Type type = var.get_output_type();
579 if (type == PISM_NAT) {
580 type = default_type;
581 }
582 file.define_variable(name, type, dims);
583
584 write_attributes(file, var, type);
585
586 // add the "grid_mapping" attribute if the grid has an associated mapping. Variables lat, lon,
587 // lat_bnds, and lon_bnds should not have the grid_mapping attribute to support CDO (see issue
588 // #384).
589 const VariableMetadata &mapping = grid.get_mapping_info().mapping;
590 if (mapping.has_attributes() and not
591 (name == "lat_bnds" or name == "lon_bnds" or name == "lat" or name == "lon")) {
592 file.write_attribute(var.get_name(), "grid_mapping",
593 mapping.get_name());
594 }
595
596 if (var.get_time_independent()) {
597 // mark this variable as "not written" so that write_spatial_variable can avoid
598 // writing it more than once.
599 file.write_attribute(var.get_name(), "not_written", PISM_INT, {1.0});
600 }
601 }
602
603 //! Read a variable from a file into an array `output`.
604 /*! This also converts data from input units to internal units if needed.
605 */
606 void read_spatial_variable(const SpatialVariableMetadata &variable,
607 const IceGrid& grid, const File &file,
608 unsigned int time, double *output) {
609
610 const Logger &log = *grid.ctx()->log();
611
612 // Find the variable:
613 auto var = file.find_variable(variable.get_name(), variable.get_string("standard_name"));
614
615 if (not var.exists) {
616 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't find '%s' (%s) in '%s'.",
617 variable.get_name().c_str(),
618 variable.get_string("standard_name").c_str(),
619 file.filename().c_str());
620 }
621
622 // Sanity check: the variable in an input file should have the expected
623 // number of spatial dimensions.
624 {
625 // Set of spatial dimensions this field has.
626 std::set<int> axes;
627 axes.insert(X_AXIS);
628 axes.insert(Y_AXIS);
629 if (not variable.get_z().get_name().empty()) {
630 axes.insert(Z_AXIS);
631 }
632
633 std::vector<std::string> input_dims;
634 int input_ndims = 0; // number of spatial dimensions (input file)
635 size_t matching_dim_count = 0; // number of matching dimensions
636
637 input_dims = file.dimensions(var.name);
638 for (auto d : input_dims) {
639 AxisType tmp = file.dimension_type(d, variable.unit_system());
640
641 if (tmp != T_AXIS) {
642 ++input_ndims;
643 }
644
645 if (axes.find(tmp) != axes.end()) {
646 ++matching_dim_count;
647 }
648 }
649
650 if (axes.size() != matching_dim_count) {
651
652 // Print the error message and stop:
653 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
654 "found the %dD variable %s (%s) in '%s' while trying to read\n"
655 "'%s' ('%s'), which is %d-dimensional.",
656 input_ndims,
657 var.name.c_str(),
658 join(input_dims, ",").c_str(),
659 file.filename().c_str(),
660 variable.get_name().c_str(),
661 variable.get_string("long_name").c_str(),
662 static_cast<int>(axes.size()));
663 }
664 }
665
666 // make sure we have at least one level
667 const std::vector<double>& zlevels = variable.get_levels();
668 unsigned int nlevels = std::max(zlevels.size(), (size_t)1);
669
670 read_distributed_array(file, grid, var.name, nlevels, time, output);
671
672 std::string input_units = file.read_text_attribute(var.name, "units");
673 const std::string &internal_units = variable.get_string("units");
674
675 if (input_units.empty() and not internal_units.empty()) {
676 const std::string &long_name = variable.get_string("long_name");
677 log.message(2,
678 "PISM WARNING: Variable '%s' ('%s') does not have the units attribute.\n"
679 " Assuming that it is in '%s'.\n",
680 variable.get_name().c_str(), long_name.c_str(),
681 internal_units.c_str());
682 input_units = internal_units;
683 }
684
685 // Convert data:
686 size_t size = grid.xm() * grid.ym() * nlevels;
687
688 units::Converter(variable.unit_system(),
689 input_units, internal_units).convert_doubles(output, size);
690 }
691
692 //! \brief Write a double array to a file.
693 /*!
694 Converts units if internal and "glaciological" units are different.
695 */
696 void write_spatial_variable(const SpatialVariableMetadata &var,
697 const IceGrid& grid,
698 const File &file,
699 const double *input) {
700
701 auto name = var.get_name();
702
703 if (not file.find_variable(name)) {
704 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't find '%s' in '%s'.",
705 name.c_str(),
706 file.filename().c_str());
707 }
708
709 write_dimensions(var, grid, file);
710
711 // avoid writing time-independent variables more than once (saves time when writing to
712 // extra_files)
713 if (var.get_time_independent()) {
714 bool written = file.attribute_type(var.get_name(), "not_written") == PISM_NAT;
715 if (written) {
716 return;
717 } else {
718 file.redef();
719 file.remove_attribute(var.get_name(), "not_written");
720 }
721 }
722
723 // make sure we have at least one level
724 unsigned int nlevels = std::max(var.get_levels().size(), (size_t)1);
725
726 std::string
727 units = var.get_string("units"),
728 glaciological_units = var.get_string("glaciological_units");
729
730 if (units != glaciological_units) {
731 size_t data_size = grid.xm() * grid.ym() * nlevels;
732
733 // create a temporary array, convert to glaciological units, and
734 // save
735 std::vector<double> tmp(data_size);
736 for (size_t k = 0; k < data_size; ++k) {
737 tmp[k] = input[k];
738 }
739
740 units::Converter(var.unit_system(),
741 units,
742 glaciological_units).convert_doubles(&tmp[0], tmp.size());
743
744 file.write_distributed_array(name, grid, nlevels, &tmp[0]);
745 } else {
746 file.write_distributed_array(name, grid, nlevels, input);
747 }
748 }
749
750 //! \brief Regrid from a NetCDF file into a distributed array `output`.
751 /*!
752 - if `flag` is `CRITICAL` or `CRITICAL_FILL_MISSING`, stops if the
753 variable was not found in the input file
754 - if `flag` is one of `CRITICAL_FILL_MISSING` and
755 `OPTIONAL_FILL_MISSING`, replace _FillValue with `default_value`.
756 - sets `v` to `default_value` if `flag` is `OPTIONAL` and the
757 variable was not found in the input file
758 - uses the last record in the file
759 */
760 void regrid_spatial_variable(SpatialVariableMetadata &var,
761 const IceGrid& grid, const File &file,
762 RegriddingFlag flag, bool report_range,
763 bool allow_extrapolation,
764 double default_value,
765 InterpolationType interpolation_type,
766 double *output) {
767 unsigned int t_length = file.nrecords(var.get_name(),
768 var.get_string("standard_name"),
769 var.unit_system());
770
771 regrid_spatial_variable(var, grid, file, t_length - 1, flag,
772 report_range, allow_extrapolation,
773 default_value, interpolation_type, output);
774 }
775
776 static void compute_range(MPI_Comm com, double *data, size_t data_size, double *min, double *max) {
777 double
778 min_result = data[0],
779 max_result = data[0];
780 for (size_t k = 0; k < data_size; ++k) {
781 min_result = std::min(min_result, data[k]);
782 max_result = std::max(max_result, data[k]);
783 }
784
785 if (min) {
786 *min = GlobalMin(com, min_result);
787 }
788
789 if (max) {
790 *max = GlobalMax(com, max_result);
791 }
792 }
793
794 /*! @brief Check that x, y, and z coordinates of the input grid are strictly increasing. */
795 void check_input_grid(const grid_info &input) {
796 if (not is_increasing(input.x)) {
797 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
798 "input x coordinate has to be strictly increasing");
799 }
800
801 if (not is_increasing(input.y)) {
802 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
803 "input y coordinate has to be strictly increasing");
804 }
805
806 if (not is_increasing(input.z)) {
807 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
808 "input vertical grid has to be strictly increasing");
809 }
810 }
811
812 /*!
813 * Check the overlap of the input grid and the internal grid.
814 *
815 * Set `allow_extrapolation` to `true` to "extend" the vertical grid during "bootstrapping".
816 */
817 static void check_grid_overlap(const grid_info &input, const IceGrid &internal,
818 const std::vector<double> &z_internal) {
819
820 // Grid spacing (assume that the grid is equally-spaced) and the
821 // extent of the domain. To compute the extent of the domain, assume
822 // that the grid is cell-centered, so edge of the domain is half of
823 // the grid spacing away from grid points at the edge.
824 //
825 // Note that x_min is not the same as internal.x(0).
826 const double
827 x_min = internal.x0() - internal.Lx(),
828 x_max = internal.x0() + internal.Lx(),
829 y_min = internal.y0() - internal.Ly(),
830 y_max = internal.y0() + internal.Ly(),
831 input_x_min = input.x0 - input.Lx,
832 input_x_max = input.x0 + input.Lx,
833 input_y_min = input.y0 - input.Ly,
834 input_y_max = input.y0 + input.Ly;
835
836 // tolerance (one micron)
837 double eps = 1e-6;
838
839 // horizontal grid extent
840 if (not (x_min >= input_x_min - eps and x_max <= input_x_max + eps and
841 y_min >= input_y_min - eps and y_max <= input_y_max + eps)) {
842 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
843 "PISM's computational domain is not a subset of the domain in '%s'\n"
844 "PISM grid: x: [%3.3f, %3.3f] y: [%3.3f, %3.3f] meters\n"
845 "input file grid: x: [%3.3f, %3.3f] y: [%3.3f, %3.3f] meters",
846 input.filename.c_str(),
847 x_min, x_max, y_min, y_max,
848 input_x_min, input_x_max, input_y_min, input_y_max);
849 }
850
851 if (z_internal.size() == 0) {
852 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
853 "Interval vertical grid has 0 levels. This should never happen.");
854 }
855
856 if (z_internal.size() == 1 and input.z.size() > 1) {
857 // internal field is 2D or 3D with one level, input variable is 3D with more than one
858 // vertical level
859 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
860 "trying to read in a 2D field but the input file %s contains\n"
861 "a 3D field with %d levels",
862 input.filename.c_str(),
863 static_cast<int>(input.z.size()));
864 }
865
866 if (z_internal.size() > 1 and input.z.size() <= 1) {
867 // internal field is 3D with more than one vertical level, input variable is 2D or 3D
868 // with 1 level
869 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
870 "trying to read in a 3D field but the input file %s contains\n"
871 "a 2D field", input.filename.c_str());
872 }
873
874 if (z_internal.size() > 1 and input.z.size() > 0) {
875 // both internal field and input variable are 3D: check vertical grid extent
876 // Note: in PISM 2D fields have one vertical level (z = 0).
877 const double
878 input_z_min = input.z.front(),
879 input_z_max = input.z.back(),
880 z_min = z_internal.front(),
881 z_max = z_internal.back();
882
883 if (not (z_min >= input.z_min - eps and z_max <= input.z_max + eps)) {
884 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
885 "PISM's computational domain is not a subset of the domain in '%s'\n"
886 "PISM grid: z: [%3.3f, %3.3f] meters\n"
887 "input file grid: z: [%3.3f, %3.3f] meters",
888 input.filename.c_str(),
889 z_min, z_max,
890 input_z_min, input_z_max);
891 }
892 }
893 }
894
895 void regrid_spatial_variable(SpatialVariableMetadata &variable,
896 const IceGrid& grid, const File &file,
897 unsigned int t_start, RegriddingFlag flag,
898 bool report_range,
899 bool allow_extrapolation,
900 double default_value,
901 InterpolationType interpolation_type,
902 double *output) {
903 const Logger &log = *grid.ctx()->log();
904
905 units::System::Ptr sys = variable.unit_system();
906 const std::vector<double>& levels = variable.get_levels();
907 const size_t data_size = grid.xm() * grid.ym() * levels.size();
908
909 // Find the variable
910 auto var = file.find_variable(variable.get_name(), variable.get_string("standard_name"));
911
912 if (var.exists) { // the variable was found successfully
913
914 {
915 grid_info input_grid(file, var.name, sys, grid.registration());
916
917 check_input_grid(input_grid);
918
919 if (not allow_extrapolation) {
920 check_grid_overlap(input_grid, grid, levels);
921 }
922 }
923
924 if (flag == OPTIONAL_FILL_MISSING or flag == CRITICAL_FILL_MISSING) {
925 log.message(2,
926 "PISM WARNING: Replacing missing values with %f [%s] in variable '%s' read from '%s'.\n",
927 default_value, variable.get_string("units").c_str(), variable.get_name().c_str(),
928 file.filename().c_str());
929
930 regrid_vec_fill_missing(file, grid, var.name, levels,
931 t_start, default_value, interpolation_type, output);
932 } else {
933 regrid_vec(file, grid, var.name, levels, t_start, interpolation_type, output);
934 }
935
936 // Now we need to get the units string from the file and convert
937 // the units, because check_range and report_range expect data to
938 // be in PISM (MKS) units.
939
940 std::string input_units = file.read_text_attribute(var.name, "units");
941 std::string internal_units = variable.get_string("units");
942
943 if (input_units.empty() and not internal_units.empty()) {
944 log.message(2,
945 "PISM WARNING: Variable '%s' ('%s') does not have the units attribute.\n"
946 " Assuming that it is in '%s'.\n",
947 variable.get_name().c_str(),
948 variable.get_string("long_name").c_str(),
949 internal_units.c_str());
950 input_units = internal_units;
951 }
952
953 // Convert data:
954 units::Converter(sys, input_units, internal_units).convert_doubles(output, data_size);
955
956 // Check the range and report it if necessary.
957 {
958 double min = 0.0, max = 0.0;
959 read_valid_range(file, var.name, variable);
960
961 compute_range(grid.com, output, data_size, &min, &max);
962
963 // Check the range and warn the user if needed:
964 variable.check_range(file.filename(), min, max);
965 if (report_range) {
966 // We can report the success, and the range now:
967 log.message(2, " FOUND ");
968
969 variable.report_range(log, min, max, var.found_using_standard_name);
970 }
971 }
972 } else { // couldn't find the variable
973 if (flag == CRITICAL or flag == CRITICAL_FILL_MISSING) {
974 // if it's critical, print an error message and stop
975 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't find '%s' in the regridding file '%s'.",
976 variable.get_name().c_str(), file.filename().c_str());
977 }
978
979 // If it is optional, fill with the provided default value.
980 // units::Converter constructor will make sure that units are compatible.
981 units::Converter c(sys,
982 variable.get_string("units"),
983 variable.get_string("glaciological_units"));
984
985 std::string spacer(variable.get_name().size(), ' ');
986 log.message(2,
987 " absent %s / %-10s\n"
988 " %s \\ not found; using default constant %7.2f (%s)\n",
989 variable.get_name().c_str(),
990 variable.get_string("long_name").c_str(),
991 spacer.c_str(), c(default_value),
992 variable.get_string("glaciological_units").c_str());
993
994 for (size_t k = 0; k < data_size; ++k) {
995 output[k] = default_value;
996 }
997 } // end of if (exists)
998 }
999
1000
1001
1002 //! Define a NetCDF variable corresponding to a time-series.
1003 void define_timeseries(const TimeseriesMetadata& var,
1004 const File &file, IO_Type nctype) {
1005
1006 std::string name = var.get_name();
1007 std::string dimension_name = var.get_dimension_name();
1008
1009 if (file.find_variable(name)) {
1010 return;
1011 }
1012
1013 if (not file.find_dimension(dimension_name)) {
1014 define_dimension(file, PISM_UNLIMITED,
1015 VariableMetadata(dimension_name, var.unit_system()));
1016 }
1017
1018 if (not file.find_variable(name)) {
1019 file.define_variable(name, nctype, {dimension_name});
1020 }
1021
1022 write_attributes(file, var, nctype);
1023 }
1024
1025 //! Read a time-series variable from a NetCDF file to a vector of doubles.
1026 void read_timeseries(const File &file, const TimeseriesMetadata &metadata,
1027 const Time &time, const Logger &log, std::vector<double> &data) {
1028
1029 std::string name = metadata.get_name();
1030
1031 try {
1032 // Find the variable:
1033 std::string
1034 long_name = metadata.get_string("long_name"),
1035 standard_name = metadata.get_string("standard_name"),
1036 dimension_name = metadata.get_dimension_name();
1037
1038 auto var = file.find_variable(name, standard_name);
1039
1040 if (not var.exists) {
1041 throw RuntimeError(PISM_ERROR_LOCATION, "variable " + name + " is missing");
1042 }
1043
1044 std::vector<std::string> dims = file.dimensions(var.name);
1045 if (dims.size() != 1) {
1046 throw RuntimeError(PISM_ERROR_LOCATION, "a time-series variable has to be one-dimensional");
1047 }
1048
1049 unsigned int length = file.dimension_length(dimension_name);
1050 if (length <= 0) {
1051 throw RuntimeError(PISM_ERROR_LOCATION, "dimension " + dimension_name + " has length zero");
1052 }
1053
1054 data.resize(length); // memory allocation happens here
1055
1056 file.read_variable(var.name, {0}, {length}, data.data());
1057
1058 units::System::Ptr system = metadata.unit_system();
1059 bool input_has_units = false;
1060 units::Unit internal_units(system, metadata.get_string("units")),
1061 input_units(system, "1");
1062
1063 std::string input_units_string = file.read_text_attribute(var.name, "units");
1064
1065 if (input_units_string.empty()) {
1066 input_has_units = false;
1067 } else {
1068 input_units_string = time.CF_units_to_PISM_units(input_units_string);
1069
1070 input_units = units::Unit(system, input_units_string);
1071 input_has_units = true;
1072 }
1073
1074 if (metadata.has_attribute("units") && not input_has_units) {
1075 std::string units_string = internal_units.format();
1076 log.message(2,
1077 "PISM WARNING: Variable '%s' ('%s') does not have the units attribute.\n"
1078 " Assuming that it is in '%s'.\n",
1079 name.c_str(), long_name.c_str(),
1080 units_string.c_str());
1081 input_units = internal_units;
1082 }
1083
1084 units::Converter(input_units, internal_units).convert_doubles(&data[0], data.size());
1085
1086 } catch (RuntimeError &e) {
1087 e.add_context("reading time-series variable '%s' from '%s'", name.c_str(),
1088 file.filename().c_str());
1089 throw;
1090 }
1091 }
1092
1093 void write_timeseries(const File &file, const TimeseriesMetadata &metadata, size_t t_start,
1094 double data, IO_Type nctype) {
1095 std::vector<double> vector_data(1, data);
1096
1097 // this call will handle errors
1098 write_timeseries(file, metadata, t_start, vector_data, nctype);
1099 }
1100
1101 /** @brief Write a time-series `data` to a file.
1102 *
1103 * Always use glaciological units when saving time-series.
1104 */
1105 void write_timeseries(const File &file, const TimeseriesMetadata &metadata, size_t t_start,
1106 const std::vector<double> &data,
1107 IO_Type nctype) {
1108
1109 std::string name = metadata.get_name();
1110 try {
1111 if (not file.find_variable(name)) {
1112 define_timeseries(metadata, file, nctype);
1113 }
1114
1115 // create a copy of "data":
1116 std::vector<double> tmp = data;
1117
1118 units::System::Ptr system = metadata.unit_system();
1119 // convert to glaciological units:
1120 units::Converter(system,
1121 metadata.get_string("units"),
1122 metadata.get_string("glaciological_units")).convert_doubles(&tmp[0], tmp.size());
1123
1124 file.write_variable(name, {(unsigned int)t_start}, {(unsigned int)tmp.size()}, tmp.data());
1125
1126 } catch (RuntimeError &e) {
1127 e.add_context("writing time-series variable '%s' to '%s'", name.c_str(),
1128 file.filename().c_str());
1129 throw;
1130 }
1131 }
1132
1133 void define_time_bounds(const TimeBoundsMetadata& var,
1134 const File &file, IO_Type nctype) {
1135 std::string name = var.get_name();
1136 std::string dimension_name = var.get_dimension_name();
1137 std::string bounds_name = var.get_bounds_name();
1138
1139 if (file.find_variable(name)) {
1140 return;
1141 }
1142
1143 if (not file.find_dimension(dimension_name)) {
1144 file.define_dimension(dimension_name, PISM_UNLIMITED);
1145 }
1146
1147 if (not file.find_dimension(bounds_name)) {
1148 file.define_dimension(bounds_name, 2);
1149 }
1150
1151 file.define_variable(name, nctype, {dimension_name, bounds_name});
1152
1153 write_attributes(file, var, nctype);
1154 }
1155
1156 void read_time_bounds(const File &file,
1157 const TimeBoundsMetadata &metadata,
1158 const Time &time, const Logger &log,
1159 std::vector<double> &data) {
1160
1161 std::string name = metadata.get_name();
1162
1163 try {
1164 units::System::Ptr system = metadata.unit_system();
1165 units::Unit internal_units(system, metadata.get_string("units"));
1166
1167 // Find the variable:
1168 if (not file.find_variable(name)) {
1169 throw RuntimeError(PISM_ERROR_LOCATION, "variable " + name + " is missing");
1170 }
1171
1172 std::vector<std::string> dims = file.dimensions(name);
1173
1174 if (dims.size() != 2) {
1175 throw RuntimeError(PISM_ERROR_LOCATION, "variable " + name + " has to has two dimensions");
1176 }
1177
1178 std::string
1179 &dimension_name = dims[0],
1180 &bounds_name = dims[1];
1181
1182 // Check that we have 2 vertices (interval end-points) per time record.
1183 unsigned int length = file.dimension_length(bounds_name);
1184 if (length != 2) {
1185 throw RuntimeError(PISM_ERROR_LOCATION, "time-bounds variable " + name + " has to have exactly 2 bounds per time record");
1186 }
1187
1188 // Get the number of time records.
1189 length = file.dimension_length(dimension_name);
1190 if (length <= 0) {
1191 throw RuntimeError(PISM_ERROR_LOCATION, "dimension " + dimension_name + " has length zero");
1192 }
1193
1194 data.resize(2*length); // memory allocation happens here
1195
1196 file.read_variable(name, {0, 0}, {length, 2}, data.data());
1197
1198 // Find the corresponding 'time' variable. (We get units from the 'time'
1199 // variable, because according to CF-1.5 section 7.1 a "boundary variable"
1200 // may not have metadata set.)
1201 if (not file.find_variable(dimension_name)) {
1202 throw RuntimeError(PISM_ERROR_LOCATION, "time coordinate variable " + dimension_name + " is missing");
1203 }
1204
1205 bool input_has_units = false;
1206 units::Unit input_units(internal_units.get_system(), "1");
1207
1208 std::string input_units_string = file.read_text_attribute(dimension_name, "units");
1209 input_units_string = time.CF_units_to_PISM_units(input_units_string);
1210
1211 if (input_units_string.empty() == true) {
1212 input_has_units = false;
1213 } else {
1214 input_units = units::Unit(internal_units.get_system(), input_units_string);
1215 input_has_units = true;
1216 }
1217
1218 if (metadata.has_attribute("units") && not input_has_units) {
1219 std::string units_string = internal_units.format();
1220 log.message(2,
1221 "PISM WARNING: Variable '%s' does not have the units attribute.\n"
1222 " Assuming that it is in '%s'.\n",
1223 dimension_name.c_str(),
1224 units_string.c_str());
1225 input_units = internal_units;
1226 }
1227
1228 units::Converter(input_units, internal_units).convert_doubles(&data[0], data.size());
1229
1230 // FIXME: check that time intervals described by the time bounds
1231 // variable are contiguous (without gaps) and stop if they are not.
1232 } catch (RuntimeError &e) {
1233 e.add_context("reading time bounds variable '%s' from '%s'", name.c_str(),
1234 file.filename().c_str());
1235 throw;
1236 }
1237 }
1238
1239 void write_time_bounds(const File &file, const TimeBoundsMetadata &metadata,
1240 size_t t_start, const std::vector<double> &data, IO_Type nctype) {
1241 std::string name = metadata.get_name();
1242 try {
1243 bool variable_exists = file.find_variable(name);
1244 if (not variable_exists) {
1245 define_time_bounds(metadata, file, nctype);
1246 }
1247
1248 // make a copy of "data"
1249 std::vector<double> tmp = data;
1250
1251 // convert to glaciological units:
1252 units::System::Ptr system = metadata.unit_system();
1253 units::Converter(system,
1254 metadata.get_string("units"),
1255 metadata.get_string("glaciological_units")).convert_doubles(&tmp[0], tmp.size());
1256
1257 std::vector<unsigned int>
1258 start{static_cast<unsigned int>(t_start), 0},
1259 count{static_cast<unsigned int>(tmp.size()) / 2, 2};
1260
1261 file.write_variable(name, start, count, &tmp[0]);
1262
1263 } catch (RuntimeError &e) {
1264 e.add_context("writing time-bounds variable '%s' to '%s'", name.c_str(),
1265 file.filename().c_str());
1266 throw;
1267 }
1268 }
1269
1270 bool file_exists(MPI_Comm com, const std::string &filename) {
1271 int file_exists_flag = 0, rank = 0;
1272 MPI_Comm_rank(com, &rank);
1273
1274 if (rank == 0) {
1275 // Check if the file exists:
1276 if (FILE *f = fopen(filename.c_str(), "r")) {
1277 file_exists_flag = 1;
1278 fclose(f);
1279 } else {
1280 file_exists_flag = 0;
1281 }
1282 }
1283 MPI_Bcast(&file_exists_flag, 1, MPI_INT, 0, com);
1284
1285 if (file_exists_flag == 1) {
1286 return true;
1287 } else {
1288 return false;
1289 }
1290 }
1291
1292 void read_attributes(const File &file,
1293 const std::string &variable_name,
1294 VariableMetadata &variable) {
1295 try {
1296 bool variable_exists = file.find_variable(variable_name);
1297
1298 if (not variable_exists) {
1299 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "variable \"%s\" is missing", variable_name.c_str());
1300 }
1301
1302 variable.clear_all_strings();
1303 variable.clear_all_doubles();
1304
1305 unsigned int nattrs = file.nattributes(variable_name);
1306
1307 for (unsigned int j = 0; j < nattrs; ++j) {
1308 std::string attribute_name = file.attribute_name(variable_name, j);
1309 IO_Type nctype = file.attribute_type(variable_name, attribute_name);
1310
1311 if (nctype == PISM_CHAR) {
1312 variable.set_string(attribute_name,
1313 file.read_text_attribute(variable_name, attribute_name));
1314 } else {
1315 variable.set_numbers(attribute_name,
1316 file.read_double_attribute(variable_name, attribute_name));
1317 }
1318 } // end of for (int j = 0; j < nattrs; ++j)
1319 } catch (RuntimeError &e) {
1320 e.add_context("reading attributes of variable '%s' from '%s'",
1321 variable_name.c_str(), file.filename().c_str());
1322 throw;
1323 }
1324 }
1325
1326 //! Write variable attributes to a NetCDF file.
1327 /*!
1328 - If both valid_min and valid_max are set, then valid_range is written
1329 instead of the valid_min, valid_max pair.
1330
1331 - Skips empty text attributes.
1332 */
1333 void write_attributes(const File &file, const VariableMetadata &variable, IO_Type nctype) {
1334 std::string var_name = variable.get_name();
1335
1336 try {
1337 std::string
1338 units = variable.get_string("units"),
1339 glaciological_units = variable.get_string("glaciological_units");
1340
1341 bool use_glaciological_units = units != glaciological_units;
1342
1343 // units, valid_min, valid_max and valid_range need special treatment:
1344 if (variable.has_attribute("units")) {
1345 std::string output_units = use_glaciological_units ? glaciological_units : units;
1346
1347 file.write_attribute(var_name, "units", output_units);
1348 }
1349
1350 std::vector<double> bounds(2);
1351 if (variable.has_attribute("valid_range")) {
1352 bounds = variable.get_numbers("valid_range");
1353 } else {
1354 if (variable.has_attribute("valid_min")) {
1355 bounds[0] = variable.get_number("valid_min");
1356 }
1357 if (variable.has_attribute("valid_max")) {
1358 bounds[1] = variable.get_number("valid_max");
1359 }
1360 }
1361
1362 double fill_value = 0.0;
1363 if (variable.has_attribute("_FillValue")) {
1364 fill_value = variable.get_number("_FillValue");
1365 }
1366
1367 // We need to save valid_min, valid_max and valid_range in the units
1368 // matching the ones in the output.
1369 if (use_glaciological_units) {
1370
1371 units::Converter c(variable.unit_system(), units, glaciological_units);
1372
1373 bounds[0] = c(bounds[0]);
1374 bounds[1] = c(bounds[1]);
1375 fill_value = c(fill_value);
1376 }
1377
1378 if (variable.has_attribute("_FillValue")) {
1379 file.write_attribute(var_name, "_FillValue", nctype, {fill_value});
1380 }
1381
1382 if (variable.has_attribute("valid_range")) {
1383 file.write_attribute(var_name, "valid_range", nctype, bounds);
1384 } else if (variable.has_attribute("valid_min") and
1385 variable.has_attribute("valid_max")) {
1386 file.write_attribute(var_name, "valid_range", nctype, bounds);
1387 } else if (variable.has_attribute("valid_min")) {
1388 file.write_attribute(var_name, "valid_min", nctype, {bounds[0]});
1389 } else if (variable.has_attribute("valid_max")) {
1390 file.write_attribute(var_name, "valid_max", nctype, {bounds[1]});
1391 }
1392
1393 // Write text attributes:
1394 for (auto s : variable.get_all_strings()) {
1395 std::string
1396 name = s.first,
1397 value = s.second;
1398
1399 if (name == "units" or
1400 name == "glaciological_units" or
1401 value.empty()) {
1402 continue;
1403 }
1404
1405 file.write_attribute(var_name, name, value);
1406 }
1407
1408 // Write double attributes:
1409 for (auto d : variable.get_all_doubles()) {
1410 std::string name = d.first;
1411 std::vector<double> values = d.second;
1412
1413 if (name == "valid_min" or
1414 name == "valid_max" or
1415 name == "valid_range" or
1416 name == "_FillValue" or
1417 values.empty()) {
1418 continue;
1419 }
1420
1421 file.write_attribute(var_name, name, nctype, values);
1422 }
1423
1424 } catch (RuntimeError &e) {
1425 e.add_context("writing attributes of variable '%s' to '%s'",
1426 var_name.c_str(), file.filename().c_str());
1427 throw;
1428 }
1429 }
1430
1431 //! Read the valid range information from a file.
1432 /*! Reads `valid_min`, `valid_max` and `valid_range` attributes; if \c
1433 valid_range is found, sets the pair `valid_min` and `valid_max` instead.
1434 */
1435 void read_valid_range(const File &file, const std::string &name, VariableMetadata &variable) {
1436 try {
1437 // Never reset valid_min/max if they were set internally
1438 if (variable.has_attribute("valid_min") or
1439 variable.has_attribute("valid_max")) {
1440 return;
1441 }
1442
1443 // Read the units.
1444 std::string file_units = file.read_text_attribute(name, "units");
1445
1446 if (file_units.empty()) {
1447 // If the variable in the file does not have the units attribute we assume that
1448 // units in the file match internal (PISM) units.
1449 file_units = variable.get_string("units");
1450 }
1451
1452 units::Converter c(variable.unit_system(),
1453 file_units,
1454 variable.get_string("units"));
1455
1456 std::vector<double> bounds = file.read_double_attribute(name, "valid_range");
1457 if (bounds.size() == 2) { // valid_range is present
1458 variable.set_number("valid_min", c(bounds[0]));
1459 variable.set_number("valid_max", c(bounds[1]));
1460 } else { // valid_range has the wrong length or is missing
1461 bounds = file.read_double_attribute(name, "valid_min");
1462 if (bounds.size() == 1) { // valid_min is present
1463 variable.set_number("valid_min", c(bounds[0]));
1464 }
1465
1466 bounds = file.read_double_attribute(name, "valid_max");
1467 if (bounds.size() == 1) { // valid_max is present
1468 variable.set_number("valid_max", c(bounds[0]));
1469 }
1470 }
1471 } catch (RuntimeError &e) {
1472 e.add_context("reading valid range of variable '%s' from '%s'", name.c_str(),
1473 file.filename().c_str());
1474 throw;
1475 }
1476 }
1477
1478 /*!
1479 * Return the name of the time dimension corresponding to a NetCDF variable.
1480 *
1481 * Returns an empty string if this variable is time-independent.
1482 */
1483 std::string time_dimension(units::System::Ptr unit_system,
1484 const File &file,
1485 const std::string &variable_name) {
1486
1487 auto dims = file.dimensions(variable_name);
1488
1489 for (auto d : dims) {
1490 if (file.dimension_type(d, unit_system) == T_AXIS) {
1491 return d;
1492 }
1493 }
1494
1495 return "";
1496 }
1497
1498 //! \brief Moves the file aside (file.nc -> file.nc~).
1499 /*!
1500 * Note: only one processor does the renaming.
1501 */
1502 void move_if_exists(MPI_Comm com, const std::string &file_to_move, int rank_to_use) {
1503 int stat = 0, rank = 0;
1504 MPI_Comm_rank(com, &rank);
1505 std::string backup_filename = file_to_move + "~";
1506
1507 if (rank == rank_to_use) {
1508 bool exists = false;
1509
1510 // Check if the file exists:
1511 if (FILE *f = fopen(file_to_move.c_str(), "r")) {
1512 fclose(f);
1513 exists = true;
1514 } else {
1515 exists = false;
1516 }
1517
1518 if (exists) {
1519 stat = rename(file_to_move.c_str(), backup_filename.c_str());
1520 }
1521 } // end of "if (rank == rank_to_use)"
1522
1523 int global_stat = 0;
1524 MPI_Allreduce(&stat, &global_stat, 1, MPI_INT, MPI_SUM, com);
1525
1526 if (global_stat != 0) {
1527 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1528 "PISM ERROR: can't move '%s' to '%s'",
1529 file_to_move.c_str(), backup_filename.c_str());
1530 }
1531 }
1532
1533 //! \brief Check if a file is present are remove it.
1534 /*!
1535 * Note: only processor 0 does the job.
1536 */
1537 void remove_if_exists(MPI_Comm com, const std::string &file_to_remove, int rank_to_use) {
1538 int stat = 0, rank = 0;
1539 MPI_Comm_rank(com, &rank);
1540
1541 if (rank == rank_to_use) {
1542 bool exists = false;
1543
1544 // Check if the file exists:
1545 if (FILE *f = fopen(file_to_remove.c_str(), "r")) {
1546 fclose(f);
1547 exists = true;
1548 } else {
1549 exists = false;
1550 }
1551
1552 if (exists) {
1553 stat = remove(file_to_remove.c_str());
1554 }
1555 } // end of "if (rank == rank_to_use)"
1556
1557 int global_stat = 0;
1558 MPI_Allreduce(&stat, &global_stat, 1, MPI_INT, MPI_SUM, com);
1559
1560 if (global_stat != 0) {
1561 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1562 "PISM ERROR: can't remove '%s'", file_to_remove.c_str());
1563 }
1564 }
1565
1566 } // end of namespace io
1567 } // end of namespace pism