tTime_Calendar.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
---
tTime_Calendar.cc (18374B)
---
1 // Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 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 2 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 <cassert>
20 #include <sstream>
21 #include <cstdlib>
22 #include <petscsys.h>
23
24 #include "error_handling.hh"
25
26 #include "Time_Calendar.hh"
27 #include "pism_options.hh"
28 #include "pism/util/io/File.hh"
29 #include "pism/external/calcalcs/utCalendar2_cal.h"
30 #include "pism/external/calcalcs/calcalcs.h"
31 #include "ConfigInterface.hh"
32 #include "VariableMetadata.hh"
33 #include "io/io_helpers.hh"
34 #include "pism/util/Logger.hh"
35
36 namespace pism {
37
38 static inline std::string string_strip(std::string input) {
39 if (input.empty() == true) {
40 return input;
41 }
42
43 // strip leading spaces
44 input.erase(0, input.find_first_not_of(" \t"));
45
46 // strip trailing spaces
47 input.substr(input.find_last_not_of(" \t"));
48
49 return input;
50 }
51
52 /*!
53
54 See http://meteora.ucsd.edu/~pierce/calcalcs/index.html and
55 http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#calendar
56
57 for more details about supported calendars.
58 */
59 Time_Calendar::Time_Calendar(MPI_Comm c, Config::ConstPtr conf,
60 const std::string &calendar_string,
61 units::System::Ptr units_system)
62 : Time(conf, calendar_string, units_system),
63 m_com(c) {
64
65 std::string ref_date = m_config->get_string("time.reference_date");
66
67 try {
68 parse_date(ref_date, NULL);
69 } catch (RuntimeError &e) {
70 e.add_context("validating the reference date");
71 throw;
72 }
73
74 try {
75 m_time_units = units::Unit(m_time_units.get_system(), "seconds since " + ref_date);
76 } catch (RuntimeError &e) {
77 e.add_context("setting time units");
78 throw;
79 }
80
81 m_run_start = increment_date(0, (int)m_config->get_number("time.start_year"));
82 m_run_end = increment_date(m_run_start, (int)m_config->get_number("time.run_length"));
83
84 m_time_in_seconds = m_run_start;
85 }
86
87 Time_Calendar::~Time_Calendar() {
88 }
89
90 bool Time_Calendar::process_ys(double &result) {
91
92 options::String ys("-ys", "Start date");
93
94 if (ys.is_set()) {
95 try {
96 parse_date(ys, &result);
97 } catch (RuntimeError &e) {
98 e.add_context("processing the -ys option");
99 throw;
100 }
101 } else {
102 result = m_config->get_number("time.start_year", "seconds");
103 }
104 return ys.is_set();
105 }
106
107 bool Time_Calendar::process_y(double &result) {
108
109 options::Integer y("-y", "Run length, in years (integer)", 0);
110
111 if (y.is_set()) {
112 if (y < 0) {
113 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "-y %d is not allowed (run length can't be negative)",
114 y.value());
115 }
116 result = years_to_seconds(y);
117 } else {
118 result = m_config->get_number("time.run_length", "seconds");
119 }
120 return y.is_set();
121 }
122
123
124 bool Time_Calendar::process_ye(double &result) {
125
126 options::String ye("-ye", "Start date");
127
128 if (ye.is_set()) {
129 try {
130 parse_date(ye, &result);
131 } catch (RuntimeError &e) {
132 e.add_context("processing the -ye option");
133 throw;
134 }
135 } else {
136 result = (m_config->get_number("time.start_year", "seconds") +
137 m_config->get_number("time.run_length", "seconds"));
138 }
139 return ye.is_set();
140 }
141
142 void Time_Calendar::init_from_input_file(const File &nc,
143 const std::string &time_name,
144 const Logger &log) {
145 try {
146 // Set the calendar name from file, unless we are re-starting from a PISM run using the "none"
147 // calendar.
148 std::string new_calendar = nc.read_text_attribute(time_name, "calendar");
149 if (not new_calendar.empty() and
150 not (new_calendar == "none")) {
151 init_calendar(new_calendar);
152 }
153
154 // Set the reference date of internal units.
155 {
156 std::string date_string = reference_date_from_file(nc, time_name);
157 m_time_units = units::Unit(m_unit_system, "seconds " + date_string);
158 }
159
160 // Read time information from the file. (PISM output files don't have time bounds, so we don't
161 // bother checking for them.)
162 std::vector<double> time;
163 {
164 TimeseriesMetadata time_axis(time_name, time_name, m_unit_system);
165 time_axis.set_string("units", m_time_units.format());
166
167 io::read_timeseries(nc, time_axis, *this, log, time);
168 }
169
170 // Set time.
171 this->set_start(time.front());
172 this->set(time.front());
173 log.message(2,
174 "* Time t = %s (calendar: %s) found in '%s'; setting current time\n",
175 this->date().c_str(),
176 this->calendar().c_str(),
177 nc.filename().c_str());
178 } catch (RuntimeError &e) {
179 e.add_context("initializing model time from \"%s\"", nc.filename().c_str());
180 throw;
181 }
182 }
183
184 void Time_Calendar::init(const Logger &log) {
185
186 // process command-line options -y, -ys, -ye
187 Time::init(log);
188
189 options::String time_file("-time_file", "Reads time information from a file");
190
191 if (time_file.is_set()) {
192 log.message(2,
193 "* Setting time from '%s'...\n",
194 time_file->c_str());
195
196 options::ignored(log, "-y");
197 options::ignored(log, "-ys");
198 options::ignored(log, "-ye");
199
200 bool continue_run = options::Bool("-time_file_continue_run",
201 "continue a run using start time in the -i file");
202 init_from_file(time_file, log, not continue_run);
203 }
204 }
205
206 //! \brief Sets the time from a NetCDF file with a time dimension (`-time_file`).
207 /*!
208 * This allows running PISM for the duration of the available forcing.
209 */
210 void Time_Calendar::init_from_file(const std::string &filename, const Logger &log,
211 bool set_start_time) {
212 try {
213 std::string time_name = m_config->get_string("time.dimension_name");
214
215 File file(m_com, filename, PISM_NETCDF3, PISM_READONLY); // OK to use netcdf3
216
217 // Set the calendar name from file.
218 std::string new_calendar = file.read_text_attribute(time_name, "calendar");
219 if (not new_calendar.empty()) {
220 init_calendar(new_calendar);
221 }
222
223 // Set the reference date of internal units.
224 {
225 std::string date_string = reference_date_from_file(file, time_name);
226 m_time_units = units::Unit(m_unit_system, "seconds " + date_string);
227 }
228
229 // Read time information from the file.
230 std::vector<double> time;
231 std::string time_bounds_name = file.read_text_attribute(time_name, "bounds");
232 if (not time_bounds_name.empty()) {
233 // use the time bounds
234 TimeBoundsMetadata bounds(time_bounds_name, time_name, m_unit_system);
235 bounds.set_string("units", m_time_units.format());
236
237 io::read_time_bounds(file, bounds, *this, log, time);
238 } else {
239 // use the time axis
240 TimeseriesMetadata time_axis(time_name, time_name, m_unit_system);
241 time_axis.set_string("units", m_time_units.format());
242
243 io::read_timeseries(file, time_axis, *this, log, time);
244 }
245
246 // Set time.
247 if (set_start_time) {
248 this->set_start(time.front());
249 this->set(time.front());
250 } else {
251 log.message(2, "* Using start time from an -i file to continue an interrupted run.\n");
252 }
253 this->set_end(time.back());
254 } catch (RuntimeError &e) {
255 e.add_context("initializing model time from \"%s\"", filename.c_str());
256 throw;
257 }
258 }
259
260 double Time_Calendar::mod(double time, unsigned int) const {
261 // This class does not support the "mod" operation.
262 return time;
263 }
264
265 double Time_Calendar::year_fraction(double T) const {
266 int year, month, day, hour, minute;
267 double second, year_start, next_year_start;
268
269 utCalendar2_cal(T, m_time_units.get(),
270 &year, &month, &day, &hour, &minute, &second,
271 m_calendar_string.c_str());
272
273 utInvCalendar2_cal(year,
274 1, 1, // month, day
275 0, 0, 0, // hour, minute, second
276 m_time_units.get(),
277 &year_start,
278 m_calendar_string.c_str());
279
280 utInvCalendar2_cal(year + 1,
281 1, 1, // month, day
282 0, 0, 0, // hour, minute, second
283 m_time_units.get(),
284 &next_year_start,
285 m_calendar_string.c_str());
286
287 return (T - year_start) / (next_year_start - year_start);
288 }
289
290 std::string Time_Calendar::date(double T) const {
291 char tmp[256];
292 int year, month, day, hour, minute;
293 double second;
294
295 utCalendar2_cal(T, m_time_units.get(),
296 &year, &month, &day, &hour, &minute, &second,
297 m_calendar_string.c_str());
298
299 snprintf(tmp, 256, "%04d-%02d-%02d", year, month, day);
300
301 return std::string(tmp);
302 }
303
304 std::string Time_Calendar::date() const {
305 return this->date(m_time_in_seconds);
306 }
307
308 std::string Time_Calendar::start_date() const {
309 return this->date(m_run_start);
310 }
311
312 std::string Time_Calendar::end_date() const {
313 return this->date(m_run_end);
314 }
315
316 double Time_Calendar::calendar_year_start(double T) const {
317 int year, month, day, hour, minute;
318 double second, result;
319
320 // Get the date corresponding to time T:
321 utCalendar2_cal(T, m_time_units.get(),
322 &year, &month, &day, &hour, &minute, &second,
323 m_calendar_string.c_str());
324
325 // Get the time in seconds corresponding to the beginning of the
326 // year.
327 utInvCalendar2_cal(year,
328 1, 1, 0, 0, 0.0, // month, day, hour, minute, second
329 m_time_units.get(), &result,
330 m_calendar_string.c_str());
331
332 return result;
333 }
334
335
336 double Time_Calendar::increment_date(double T, int years) const {
337 int year, month, day, hour, minute;
338 double second, result;
339
340 // Get the date corresponding ti time T:
341 utCalendar2_cal(T, m_time_units.get(),
342 &year, &month, &day, &hour, &minute, &second,
343 m_calendar_string.c_str());
344
345 calcalcs_cal *cal = NULL;
346 int errcode, leap = 0;
347 cal = ccs_init_calendar(m_calendar_string.c_str());
348 assert(cal != NULL);
349 errcode = ccs_isleap(cal, year + years, &leap);
350 assert(errcode == 0);
351 ccs_free_calendar(cal);
352
353 if (leap == 0 && month == 2 && day == 29) {
354 PetscErrorCode ierr = PetscPrintf(m_com,
355 "PISM WARNING: date %d year(s) since %d-%d-%d does not exist."
356 " Using %d-%d-%d instead of %d-%d-%d.\n",
357 years,
358 year, month, day,
359 year + years, month, day-1,
360 year + years, month, day);
361 PISM_CHK(ierr, "PetscPrintf");
362 day -= 1;
363 }
364
365 // Get the time in seconds corresponding to the new date.
366 utInvCalendar2_cal(year + years, month, day,
367 hour, minute, second, m_time_units.get(), &result,
368 m_calendar_string.c_str());
369
370 return result;
371 }
372
373 /**
374 * Parses the date string `input`.
375 *
376 * @param[in] input date string of the form "Y-M-D", where year, month,
377 * and day can use as many digits as necessary
378 * @param[out] result
379 *
380 * If `result` is NULL, then it validates the date string, but does
381 * not try to convert to seconds since the reference date. This makes
382 * it possible to validate the reference date itself.
383 *
384 */
385 void Time_Calendar::parse_date(const std::string &input, double *result) const {
386 std::vector<int> numbers;
387 bool year_is_negative = false;
388 std::string tmp, spec = input;
389
390 spec = string_strip(spec);
391
392 std::istringstream arg(spec);
393
394 if (spec.empty() == true) {
395 throw RuntimeError(PISM_ERROR_LOCATION, "got an empty date specification");
396 }
397
398 // If the string starts with "-" then the year is negative. This
399 // would confuse the code below, which treats "-" as a separator, so
400 // we remember that the year is negative and remove "-".
401 if (spec[0] == '-') {
402 year_is_negative = true;
403 spec.substr(1);
404 }
405
406 while(getline(arg, tmp, '-')) {
407
408 // an empty part in the date specification (corresponds to "--")
409 if (tmp.empty() == true) {
410 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "date specification '%s' is invalid (can't have two '-' in a row)",
411 spec.c_str());
412 }
413
414 // check if strtol can parse it:
415 char *endptr = NULL;
416 long int n = strtol(tmp.c_str(), &endptr, 10);
417 if (*endptr != '\0') {
418 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "date specification '%s' is invalid ('%s' is not an integer)",
419 spec.c_str(), tmp.c_str());
420 }
421
422 // FIXME: this may overflow!
423 numbers.push_back((int)n);
424 }
425
426 // wrong number of parts in the YYYY-MM-DD date:
427 if (numbers.size() != 3) {
428 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "date specification '%s' is invalid (should have 3 parts: YYYY-MM-DD, got %d)",
429 spec.c_str(), (int)numbers.size());
430 }
431
432 if (year_is_negative) {
433 numbers[0] *= -1;
434 }
435
436 calcalcs_cal *cal = ccs_init_calendar(m_calendar_string.c_str());
437 if (cal == NULL) {
438 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "calendar string '%s' is invalid", m_calendar_string.c_str());
439 }
440
441 int dummy = 0;
442 int errcode = ccs_date2jday(cal, numbers[0], numbers[1], numbers[2], &dummy);
443 ccs_free_calendar(cal);
444 if (errcode != 0) {
445 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "date %s is invalid in the %s calendar",
446 spec.c_str(), m_calendar_string.c_str());
447 }
448
449 // result is the *output* argument. If it is not NULL, then the user
450 // asked us to convert a date to seconds since the reference time.
451 if (result != NULL) {
452 double time = 0.0;
453 errcode = utInvCalendar2_cal(numbers[0], numbers[1], numbers[2], 0, 0, 0.0,
454 m_time_units.get(), &time, m_calendar_string.c_str());
455 if (errcode != 0) {
456 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "cannot convert '%s' to '%s'",
457 spec.c_str(), m_time_units.format().c_str());
458 }
459
460 *result = time;
461 }
462 }
463
464 void Time_Calendar::parse_interval_length(const std::string &spec,
465 std::string &keyword, double *result) const {
466
467 Time::parse_interval_length(spec, keyword, result);
468
469 // This is called *only* if the 'spec' is *not* one of "monthly",
470 // "yearly", "daily", "hourly", so we don't need to worry about
471 // spec.find("...") finding "year" in "yearly".
472
473 // do not allow intervals specified in terms of "fuzzy" units
474 if (spec.find("year") != std::string::npos || spec.find("month") != std::string::npos) {
475 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "interval length '%s' with the calendar '%s' is not supported",
476 spec.c_str(), m_calendar_string.c_str());
477 }
478 }
479
480
481 void Time_Calendar::compute_times_monthly(std::vector<double> &result) const {
482 int errcode;
483
484 int year, month, day, hour, minute;
485 double second;
486
487 double time = m_run_start;
488 // get the date corresponding to the current time
489 errcode = utCalendar2_cal(time, m_time_units.get(),
490 &year, &month, &day,
491 &hour, &minute, &second,
492 m_calendar_string.c_str());
493 PISM_C_CHK(errcode, 0, "utCalendar2_cal");
494
495 result.clear();
496 while (true) {
497 // find the time corresponding to the beginning of the current
498 // month
499 errcode = utInvCalendar2_cal(year, month, 1, // year, month, day
500 0, 0, 0.0, // hour, minute, second
501 m_time_units.get(), &time,
502 m_calendar_string.c_str());
503 PISM_C_CHK(errcode, 0, "utCalendar2_cal");
504
505 if (time > m_run_end) {
506 break;
507 }
508
509 if (time >= m_run_start && time <= m_run_end) {
510 result.push_back(time);
511 }
512
513 if (month == 12) {
514 year += 1;
515 month = 1;
516 } else {
517 month += 1;
518 }
519
520 }
521 }
522
523 void Time_Calendar::compute_times_yearly(std::vector<double> &result) const {
524 int errcode;
525
526 int year, month, day, hour, minute;
527 double second;
528
529 double time = m_run_start;
530 // get the date corresponding to the current time
531 errcode = utCalendar2_cal(time, m_time_units.get(),
532 &year, &month, &day,
533 &hour, &minute, &second,
534 m_calendar_string.c_str());
535 PISM_C_CHK(errcode, 0, "utCalendar2_cal");
536
537 result.clear();
538 while (true) {
539 // find the time corresponding to the beginning of the current
540 // year
541 errcode = utInvCalendar2_cal(year, 1, 1, // year, month, day
542 0, 0, 0.0, // hour, minute, second
543 m_time_units.get(), &time,
544 m_calendar_string.c_str());
545 PISM_C_CHK(errcode, 0, "utInvCalendar2_cal");
546
547 if (time > m_run_end) {
548 break;
549 }
550
551 if (time >= m_run_start && time <= m_run_end) {
552 result.push_back(time);
553 }
554
555 year += 1;
556 }
557 }
558
559 void Time_Calendar::compute_times(double time_start, double delta, double time_end,
560 const std::string &keyword,
561 std::vector<double> &result) const {
562 if (keyword == "simple") {
563 compute_times_simple(time_start, delta, time_end, result);
564 } else if (keyword == "monthly") {
565 compute_times_monthly(result);
566 } else if (keyword == "yearly") {
567 compute_times_yearly(result);
568 } else {
569 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "'%s' reporting is not implemented", keyword.c_str());
570 }
571 }
572
573 } // end of namespace pism