PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Time.cc
Go to the documentation of this file.
1 // Copyright (C) 2011-2022 Constantine Khroulev
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 <cmath>
20 #include <cassert> // assert
21 #include <cstring> // strlen()
22 #include <limits>
23 
24 #include "pism/util/Time.hh"
25 
26 #include "pism/external/calcalcs/calcalcs.h"
27 
28 #include "pism/util/ConfigInterface.hh"
29 #include "pism/util/VariableMetadata.hh"
30 #include "pism/util/pism_utilities.hh"
31 #include "pism/util/error_handling.hh"
32 #include "pism/util/io/File.hh"
33 #include "pism/util/io/io_helpers.hh"
34 #include "pism/util/Logger.hh"
35 #include "pism/util/io/IO_Flags.hh"
36 
37 namespace pism {
38 
39 //! Get the reference date from a file.
40 static std::string reference_date_from_file(const File &file,
41  const std::string &time_name,
42  const std::string &default_value,
43  bool stop_on_error) {
44 
45  if (file.find_variable(time_name)) {
46  std::string time_units = file.read_text_attribute(time_name, "units");
47 
48  if (not time_units.empty()) {
49  // Check if the time_units includes a reference date.
50  size_t position = time_units.find("since");
51 
52  if (position != std::string::npos) {
53  return string_strip(time_units.substr(position + strlen("since")));
54  }
55 
56  if (stop_on_error) {
57 
59  "%s:units = \"%s\" in '%s' does not contain a reference date",
60  time_name.c_str(),
61  time_units.c_str(),
62  file.filename().c_str());
63  }
64  } else if (stop_on_error) {
66  "the '%s' variable in '%s' has no units",
67  time_name.c_str(), file.filename().c_str());
68  }
69  } else if (stop_on_error) {
71  "'%s' variable is not present in '%s'.",
72  time_name.c_str(), file.filename().c_str());
73  }
74 
75  return default_value;
76 }
77 
78 //! Get the calendar name from a file.
79 static std::string calendar_from_file(const File &file,
80  const std::string &time_name,
81  const std::string &default_value,
82  bool stop_on_error) {
83 
84  if (file.find_variable(time_name)) {
85  std::string calendar_name = file.read_text_attribute(time_name, "calendar");
86 
87  if (not calendar_name.empty()) {
88  return calendar_name;
89  }
90 
91  if (stop_on_error) {
93  "the '%s' variable in '%s' has no calendar attribute",
94  time_name.c_str(), file.filename().c_str());
95  }
96 
97  } else if (stop_on_error) {
99  "'%s' variable is not present in '%s'.",
100  time_name.c_str(), file.filename().c_str());
101  }
102 
103  return default_value;
104 }
105 
106 /*!
107  * Get the reference date from a file or the configuration parameter.
108  */
109 static std::string reference_date(const File *input_file,
110  const Config &config,
111  const Logger &log) {
112 
113  auto default_reference_date = config.get_string("time.reference_date");
114 
115  if (input_file != nullptr) {
116  // input file is not empty
117 
118  auto time = config.get_string("time.dimension_name");
119 
120  if (not config.get_flag("input.bootstrap")) {
121  // restarting from a file: use the reference date in this file
122  bool stop_on_error = true;
123  return reference_date_from_file(*input_file, time, default_reference_date, stop_on_error);
124  }
125 
126  // Bootstrapping: use the configuration parameter and warn about mismatches
127  bool stop_on_error = false;
128  auto ref_date = reference_date_from_file(*input_file, time, default_reference_date, stop_on_error);
129 
130  if (ref_date != default_reference_date) {
131  log.message(2,
132  "WARNING: Using reference date %s\n"
133  " instead of the one present in the input file '%s' (%s)\n",
134  default_reference_date.c_str(), input_file->filename().c_str(), ref_date.c_str());
135  }
136 
137  return ref_date;
138  }
139 
140  return default_reference_date;
141 }
142 
143 /*!
144  * Get the calendar name from a file or a configuration parameter.
145  */
146 static std::string calendar(const File *input_file,
147  const Config &config,
148  const Logger &log) {
149  auto default_calendar = config.get_string("time.calendar");
150 
151  if (input_file != nullptr) {
152  // input file is not empty
153 
154  auto time = config.get_string("time.dimension_name");
155 
156  if (not config.get_flag("input.bootstrap")) {
157  // restarting from a file: use the calendar in this file
158  bool stop_on_error = true;
159  return calendar_from_file(*input_file, time, default_calendar, stop_on_error);
160  }
161 
162  // Bootstrapping: use the configuration parameter and warn about mismatches
163  bool stop_on_error = false;
164  auto calendar = calendar_from_file(*input_file, time, default_calendar, stop_on_error);
165 
166  if (calendar != default_calendar) {
167  log.message(2,
168  "WARNING: Using calendar %s\n"
169  " instead of the one present in the input file '%s' (%s)\n",
170  default_calendar.c_str(), input_file->filename().c_str(), calendar.c_str());
171  }
172 
173  return default_calendar;
174  }
175 
176  return default_calendar;
177 }
178 
179 /*!
180  * Increment the date corresponding to `T` by `years` years.
181  */
182 static double increment_date(const units::Unit &time_units,
183  const std::string &calendar,
184  double T, double years) {
185 
186  double whole_years_double = std::floor(years);
187  if (whole_years_double > std::numeric_limits<int>::max() or
188  whole_years_double < std::numeric_limits<int>::min()) {
190  "time offset of %f years does not fit in an 'int'",
191  whole_years_double);
192  }
193 
194  int whole_years = static_cast<int>(whole_years_double);
195  double year_fraction = years - whole_years;
196  const double day_length = 86400.0;
197 
198  // Get the date corresponding to time T:
199  auto date = time_units.date(T, calendar);
200 
201  // shift the date by the number of whole years requested
202  date.year += whole_years;
203 
204  // check if the resulting year is a leap year:
205  int leap = 0;
206  {
207  calcalcs_cal *cal = ccs_init_calendar(calendar.c_str());
208  assert(cal != NULL);
209  int errcode = ccs_isleap(cal, date.year, &leap);
210  if (errcode != 0) {
211  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "CalCalcs error: %s",
212  ccs_err_str(errcode));
213  }
214  ccs_free_calendar(cal);
215  }
216 
217  double result = 0.0;
218  if (leap == 0 and date.month == 2 and date.day == 29) {
219  // avoid passing an impossible date to UDUNITS (no February 29 in non-leap years):
220  date.day -= 1;
221  result = time_units.time(date, calendar);
222  // add back the day we substracted above
223  result += day_length;
224  } else {
225  result = time_units.time(date, calendar);
226  }
227 
228  int year_length = 360;
229  if (calendar != "360_day") {
230  year_length = (leap == 1) ? 366 : 365;
231  }
232 
233  result += year_fraction * (year_length * day_length);
234 
235  return result;
236 }
237 
238 /*!
239  * Parse the date.
240  *
241  * `input` can be
242  *
243  * - a YYYY-MM-DD date (YYYY can be negative)
244  *
245  * - a number (interpreted as the number of years since the reference date in
246  * `time_units`)
247  *
248  * - a number with units attached ("1 day", etc) interpreted as time since the reference
249  * date in `time_units`
250  */
251 static double parse_date(const std::string &input,
252  const units::Unit &time_units,
253  const std::string &calendar) {
254 
255  std::string spec = string_strip(input);
256 
257  if (spec.empty()) {
259  "got an empty date specification");
260  }
261 
262  // We need to remember if the year was negative in the input string: split() will ignore
263  // empty tokens separated by "-", so "-1000-1-1" will produce ["1000", "1", "1"].
264  bool year_is_negative = (spec[0] == '-');
265 
266  if (year_is_negative and not member(calendar, {"365_day", "360_day", "noleap"})) {
268  "negative dates such as '%s' are not allowed with the '%s' calendar.\n"
269  "Please submit a bug report if this is a problem.",
270  spec.c_str(), calendar.c_str());
271  }
272 
273  auto parts = split(spec, '-');
274 
275  if (parts.size() == 3) {
276  std::vector<int> numbers;
277  for (const auto &p : parts) {
278  try {
279  long int n = parse_integer(p);
280 
284  "%ld does not fit in an 'int'",
285  n);
286  }
287 
288  numbers.push_back(static_cast<int>(n));
289  } catch (RuntimeError &e) {
290  e.add_context("parsing a date specification %s",
291  spec.c_str());
292  throw;
293  }
294  }
295 
296  if (year_is_negative) {
297  numbers[0] *= -1;
298  }
299 
300  // Validate the calendar string and the date in this calendar:
301  {
302  calcalcs_cal *cal = ccs_init_calendar(calendar.c_str());
303  if (cal == NULL) {
305  "calendar string '%s' is invalid",
306  calendar.c_str());
307  }
308 
309  int dummy = 0;
310  int errcode = ccs_date2jday(cal, numbers[0], numbers[1], numbers[2], &dummy);
311  if (errcode != 0) {
312  ccs_free_calendar(cal);
314  "date %s is invalid in the %s calendar",
315  spec.c_str(), calendar.c_str());
316  }
317  ccs_free_calendar(cal);
318  }
319 
320  units::DateTime d{numbers[0], numbers[1], numbers[2], 0, 0, 0.0};
321 
322  return time_units.time(d, calendar);
323  } // end of the block processing dates written as Y-M-D
324 
325  // "spec" must be a number or a number with units attached to it
326  try {
327  // check if strtod() can parse it:
328  char *endptr = NULL;
329  double t = strtod(spec.c_str(), &endptr);
330  if (*endptr == '\0') {
331  // strtod() parsed it successfully: assume that it is in years. This will return
332  // time in seconds.
333  return increment_date(time_units, calendar, 0, t);
334  }
335 
336  // strtod() failed -- assume that this is a number followed by units compatible
337  // with seconds
338  return units::convert(time_units.system(), 1.0, spec, "seconds");
339  } catch (RuntimeError &e) {
340  e.add_context("parsing the date " + spec);
341  throw;
342  }
343 }
344 
345 /*!
346  * Return the start time.
347  */
348 static double start_time(const Config &config,
349  const Logger &log,
350  const File *file,
351  const std::string &reference_date,
352  const std::string &calendar,
353  const units::Unit &time_units) {
354 
355  auto time_start = config.get_string("time.start");
356 
357  if (not time_start.empty()) {
358  return parse_date(time_start, time_units, calendar);
359  }
360 
361  if (file == nullptr) {
362  // 0.0 corresponds to the reference date
363  return 0.0;
364  }
365 
366  // get the calendar in this file
367  auto time_name = config.get_string("time.dimension_name");
368  bool stop_on_error = false;
369  auto file_calendar = calendar_from_file(*file, time_name, calendar, stop_on_error);
370  auto ref_date = reference_date_from_file(*file, time_name, reference_date, stop_on_error);
371 
372  if (file_calendar != calendar) {
374  "calendar in '%s' (%s) does not match the selected calendar (%s)",
375  file->filename().c_str(), file_calendar.c_str(), calendar.c_str());
376  }
377 
378  if (ref_date != reference_date) {
380  "reference date in '%s' (%s) does not match the selected date (%s)",
381  file->filename().c_str(), ref_date.c_str(), reference_date.c_str());
382  }
383 
384  // FIXME: it would make sense to get the length of the time dimension and read the last
385  // number instead.
386  if (file->dimension_length(time_name) > 0) {
387  VariableMetadata time_axis(time_name, time_units.system());
388  time_axis["units"] = time_units.format();
389 
390  std::vector<double> time{};
391  io::read_timeseries(*file, time_axis, log, time);
392 
393  return time.back();
394  }
395 
396  return 0.0;
397 }
398 
399 static double end_time(const Config &config,
400  double time_start,
401  const std::string &calendar,
402  const units::Unit &time_units) {
403  auto time_end = config.get_string("time.end");
404 
405  if (not time_end.empty()) {
406  // parse use time_end and use it
407  return parse_date(time_end, time_units, calendar);
408  }
409 
410  // use time_start and run_length
411  auto run_length = config.get_number("time.run_length", "seconds");
412  return time_start + run_length;
413 }
414 
415 //! Convert model years into seconds using the year length
416 //! corresponding to the current calendar.
417 /*! Do not use this to convert quantities other than time intervals!
418  */
419 double Time::years_to_seconds(double input) const {
420  return input * m_year_length;
421 }
422 
423 //! Convert seconds into model years using the year length
424 //! corresponding to the current calendar.
425 /*! Do not use this to convert quantities other than time intervals!
426  */
427 double Time::seconds_to_years(double input) const {
428  return input / m_year_length;
429 }
430 
431 void Time::init_calendar(const std::string &calendar_string) {
432 
433  if (not pism_is_valid_calendar_name(calendar_string)) {
435  "unsupported calendar: %s", calendar_string.c_str());
436  }
437 
438  m_calendar_string = calendar_string;
439 
440  double seconds_per_day = convert(m_unit_system, 1.0, "day", "seconds");
441  if (calendar_string == "360_day") {
442  m_year_length = 360 * seconds_per_day;
443  } else if (member(calendar_string, {"365_day", "noleap"})) {
444  m_year_length = 365 * seconds_per_day;
445  } else {
446  // use the ~365.2524-day year
447  m_year_length = convert(m_unit_system, 1.0, "year", "seconds");
448  }
449 }
450 
451 //! \brief Sets the current time (in seconds since the reference time).
452 void Time::set(double new_time) {
453  m_time_in_seconds = new_time;
454 }
455 
456 void Time::set_start(double new_start) {
457  m_run_start = new_start;
458 }
459 
460 void Time::set_end(double new_end) {
461  m_run_end = new_end;
462 }
463 
464 //! \brief Current time, in seconds.
465 double Time::current() const {
466  return m_time_in_seconds;
467 }
468 
469 double Time::start() const {
470  return m_run_start;
471 }
472 
473 double Time::end() const {
474  return m_run_end;
475 }
476 
477 std::string Time::units_string() const {
478  return m_time_units.format();
479 }
480 
482  return m_time_units;
483 }
484 
485 //! \brief Returns the calendar string.
486 std::string Time::calendar() const {
487  return m_calendar_string;
488 }
489 
490 void Time::step(double delta_t) {
491  m_time_in_seconds += delta_t;
492 
493  // If we are less than m_t_eps seconds from the end of the run, reset
494  // m_time_in_seconds to avoid taking a very small (and useless) time step.
495  if (m_run_end > m_time_in_seconds and
498  }
499 }
500 
501 std::string Time::run_length() const {
503 }
504 
505 double Time::day_of_the_year_to_year_fraction(unsigned int day) const {
506  const double sperd = 86400.0;
507  return (sperd / m_year_length) * (double) day;
508 }
509 
510 std::vector<double> Time::parse_times(const std::string &spec) const {
511  if (spec.find(',') != std::string::npos) {
512  // a list will always contain a comma because at least two numbers are
513  // needed to specify reporting intervals
514  return parse_list(spec);
515  }
516 
517  // it must be a range specification
518  return parse_range(spec);
519 }
520 
521 std::vector<double> Time::parse_list(const std::string &spec) const {
522  std::vector<double> result;
523 
524  try {
525  for (const auto &s : split(spec, ',')) {
526  result.emplace_back(parse_date(s, m_time_units, m_calendar_string));
527  }
528  } catch (RuntimeError &e) {
529  e.add_context("parsing a list of dates %s", spec.c_str());
530  throw;
531  }
532 
533  return result;
534 }
535 
536 /**
537  * Parses an interval specification string.
538  *
539  * @param[in] spec specification string
540  *
541  */
542 auto Time::parse_interval_length(const std::string &spec) const -> Interval {
543 
544  // check if it is a keyword
545  if (spec == "hourly") {
546  return {3600.0, SIMPLE};
547  }
548 
549  if (spec == "daily") {
550  return {86400.0, SIMPLE};
551  }
552 
553  if (spec == "monthly") {
554  return {1.0, MONTHLY};
555  }
556 
557  if (spec == "yearly") {
558  return {1.0, YEARLY};
559  }
560 
561  if (not m_simple_calendar) {
562  if (spec.find("year") != std::string::npos or spec.find("month") != std::string::npos) {
564  "interval length '%s' with the calendar '%s' is not supported",
565  spec.c_str(), m_calendar_string.c_str());
566  }
567  }
568 
569  try {
570  units::Unit seconds(m_unit_system, "seconds");
571  units::Unit one(m_unit_system, "1");
572  units::Unit tmp(m_unit_system, spec);
573 
574  // Check if these units are compatible with "seconds" or "1". The
575  // latter allows intervals of the form "0.5", which stands for "half
576  // of a model year". This also discards interval specs such as "days
577  // since 1-1-1", even though "days" is compatible with "seconds".
578  if (units::are_convertible(tmp, seconds)) {
579  units::Converter c(tmp, seconds);
580 
581  return {c(1.0), SIMPLE};
582  }
583 
584  if (units::are_convertible(tmp, one)) {
585  units::Converter c(tmp, one);
586 
587  // convert from years to seconds without using UDUNITS-2 (this
588  // way we handle 360-day and 365-day years correctly)
589  return {years_to_seconds(c(1.0)), SIMPLE};
590  }
591  } catch (RuntimeError &e) {
592  e.add_context("processing interval length " + spec);
593  throw;
594  }
595 
597  "interval length '%s' is invalid", spec.c_str());
598 }
599 
600 
601 std::vector<double> Time::parse_range(const std::string &spec) const {
602  double
603  time_start = m_run_start,
604  time_end = m_run_end;
605 
606  Interval I{0.0, SIMPLE};
607 
608  if (spec == "hourly") {
609  I = {3600.0, SIMPLE};
610  } else if (spec == "daily") {
611  I = {86400.0, SIMPLE};
612  } else if (spec == "monthly") {
613  I = {1.0, MONTHLY};
614  } else if (spec == "yearly") {
615  I = {1.0, YEARLY};
616  } else {
617 
618  auto parts = pism::split(spec, ':');
619 
620  if (parts.size() == 1) {
621  I = parse_interval_length(parts[0]);
622 
623  } else if (parts.size() == 3) {
624  time_start = parse_date(parts[0], m_time_units, m_calendar_string);
625  I = parse_interval_length(parts[1]);
626  time_end = parse_date(parts[2], m_time_units, m_calendar_string);
627  } else {
629  "a time range must consist of exactly 3 parts separated by colons (got '%s').",
630  spec.c_str());
631  }
632  }
633 
634  std::vector<double> result;
635  compute_times(time_start, time_end, I, result);
636  return result;
637 }
638 
639 /**
640  * Compute times corresponding to a "simple" time range.
641  *
642  * This is a time range with a fixed step ("hourly" is an example).
643  *
644  * @param[in] time_start beginning of the time interval, in seconds
645  * @param[in] delta step, in seconds
646  * @param[in] time_end end of the interval, in seconds
647  * @param[out] result list of model times
648  *
649  */
650 void Time::compute_times_simple(double time_start, double delta, double time_end,
651  std::vector<double> &result) const {
652  if (time_start >= time_end) {
653  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "a >= b in time range a:dt:b (got a = %s, b = %s)",
654  this->date(time_start).c_str(), this->date(time_end).c_str());
655  }
656 
657  if (delta <= 0) {
658  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "dt <= 0 in time range a:dt:b (got dt = %f seconds)", delta);
659  }
660 
661  int k = 0;
662  double t = time_start;
663 
664  result.clear();
665  do {
666  if (t >= this->start() && t <= this->end()) {
667  result.push_back(t);
668  }
669  k += 1;
670  t = time_start + k * delta;
671  } while (t <= time_end);
672 }
673 
674 double Time::convert_time_interval(double T, const std::string &units) const {
675  if (member(units, {"year", "years", "yr", "a"})) {
676  return this->seconds_to_years(T); // uses year length here
677  }
678  return convert(m_unit_system, T, "seconds", units);
679 }
680 
681 /*!
682 
683  See http://meteora.ucsd.edu/~pierce/calcalcs/index.html and
684  http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#calendar
685 
686  for more details about supported calendars.
687  */
688 Time::Time(MPI_Comm com,
689  Config::ConstPtr config,
690  const Logger &log,
691  units::System::Ptr unit_system)
692  : m_config(config),
693  m_unit_system(unit_system),
694  m_time_units(unit_system, "seconds since 1-1-1") {
695 
696  m_t_eps = config->get_number("time_stepping.resolution", "seconds");
697 
698  auto input_file = config->get_string("input.file");
699 
700  std::unique_ptr<File> file{};
701  if (not input_file.empty()) {
702  file.reset(new File(com, input_file, io::PISM_NETCDF3, io::PISM_READONLY));
703  }
704 
705  // set the reference date
706  auto ref_date = reference_date(file.get(), *config, log);
707 
708  try {
709  // this will validate the reference date
710  m_time_units = units::Unit(m_unit_system, "seconds since " + ref_date);
711  } catch (RuntimeError &e) {
712  e.add_context("setting time units");
713  throw;
714  }
715 
716  m_calendar_string = ::pism::calendar(file.get(), *config, log);
718  m_simple_calendar = member(m_calendar_string, {"360_day", "365_day", "no_leap"});
719 
720  m_run_start = start_time(*config,
721  log,
722  file.get(),
723  ref_date,
725  m_time_units);
726 
727  m_run_end = end_time(*config,
728  m_run_start,
730  m_time_units);
731 
732  if (m_run_start > m_run_end) {
733  auto start = date(m_run_start);
734  auto end = date(m_run_end);
736  "Negative run length. Start: %s, end: %s.",
737  start.c_str(), end.c_str());
738  }
739 
741 
742  auto time_file = config->get_string("time.file");
743  bool continue_run = config->get_flag("time.file.continue");
744 
745  if (not time_file.empty()) {
746  log.message(2,
747  "* Setting time from '%s'...\n",
748  time_file.c_str());
749  init_from_file(com, time_file, log, not continue_run);
750  }
751 }
752 
753 //! \brief Sets the time from a NetCDF file with a time dimension (`-time_file`).
754 /*!
755  * Sets
756  * - calendar
757  * - reference date
758  * - start time
759  * - current time
760  * - end time
761  *
762  * This allows running PISM for the duration of the available forcing.
763  */
764 void Time::init_from_file(MPI_Comm com,
765  const std::string &filename,
766  const Logger &log,
767  bool set_start_time) {
768  try {
769  std::string time_name = m_config->get_string("time.dimension_name");
770 
771  File file(com, filename, io::PISM_NETCDF3, io::PISM_READONLY); // OK to use netcdf3
772 
773  // Set the calendar name from file.
774  std::string new_calendar = file.read_text_attribute(time_name, "calendar");
775  if (not new_calendar.empty()) {
776  init_calendar(new_calendar);
777  }
778 
779  // Set the reference date of internal units.
780  {
781  bool stop_on_error = true;
782  std::string date_string = reference_date_from_file(file,
783  time_name,
784  "irrelevant (not used)",
785  stop_on_error);
786  m_time_units = units::Unit(m_unit_system, "seconds since " + date_string);
787  }
788 
789  // Read time information from the file.
790  std::vector<double> time;
791  std::string time_bounds_name = file.read_text_attribute(time_name, "bounds");
792  if (not time_bounds_name.empty()) {
793  // use the time bounds
794  VariableMetadata bounds(time_bounds_name, m_unit_system);
795  bounds["units"] = m_time_units.format();
796 
797  io::read_time_bounds(file, bounds, log, time);
798  } else {
799  // use the time axis
800  VariableMetadata time_axis(time_name, m_unit_system);
801  time_axis["units"] = m_time_units.format();
802 
803  io::read_timeseries(file, time_axis, log, time);
804  }
805 
806  // Set time.
807  if (set_start_time) {
808  this->set_start(time.front());
809  this->set(time.front());
810  } else {
811  log.message(2, "* Using start time from an -i file to continue an interrupted run.\n");
812  }
813  this->set_end(time.back());
814  } catch (RuntimeError &e) {
815  e.add_context("initializing model time from \"%s\"", filename.c_str());
816  throw;
817  }
818 }
819 
820 double Time::year_fraction(double T) const {
821 
823 
824  units::DateTime D2{D.year, 1, 1, 0, 0, 0.0};
825 
826  auto year_start = m_time_units.time(D2, m_calendar_string);
827 
828  D2.year += 1;
829  auto next_year_start = m_time_units.time(D2, m_calendar_string);
830 
831  return (T - year_start) / (next_year_start - year_start);
832 }
833 
834 std::string Time::date(double T) const {
836 
837  double hour = date.hour + date.minute / 60.0 + date.second / 3600.0;
838 
839  return pism::printf("%04d-%02d-%02d %06.3fh",
840  date.year, date.month, date.day, hour);
841 }
842 
843 double Time::calendar_year_start(double T) const {
845 
846  units::DateTime D2{D.year, 1, 1, 0, 0, 0.0};
847 
849 }
850 
851 
852 double Time::increment_date(double T, double years) const {
854 }
855 
856 void Time::compute_times_monthly(std::vector<double> &result) const {
857 
858  double time = m_run_start;
859 
860  // get the date corresponding to the current time
861  auto date = m_time_units.date(time, m_calendar_string);
862 
863  // beginning of the current month
864  units::DateTime d{date.year, date.month, 1, 0, 0, 0.0};
865 
866  result.clear();
867  while (true) {
868  // find the time corresponding to the beginning of the current month
870 
871  if (time > m_run_end) {
872  break;
873  }
874 
875  if (time >= m_run_start and time <= m_run_end) {
876  result.push_back(time);
877  }
878 
879  if (d.month == 12) {
880  d.year += 1;
881  d.month = 1;
882  } else {
883  d.month += 1;
884  }
885  }
886 }
887 
888 void Time::compute_times_yearly(std::vector<double> &result) const {
889 
890  double time = m_run_start;
891  // get the date corresponding to the current time
892  auto date = m_time_units.date(time, m_calendar_string);
893 
894  units::DateTime d{date.year, 1, 1, 0, 0, 0.0};
895 
896  result.clear();
897  while (true) {
898  // find the time corresponding to the beginning of the current year
900 
901  if (time > m_run_end) {
902  break;
903  }
904 
905  if (time >= m_run_start and time <= m_run_end) {
906  result.push_back(time);
907  }
908 
909  d.year += 1;
910  }
911 }
912 
913 void Time::compute_times(double time_start, double time_end,
914  const Interval &interval,
915  std::vector<double> &result) const {
916  switch (interval.type) {
917  case SIMPLE:
918  compute_times_simple(time_start, interval.dt, time_end, result);
919  break;
920  case MONTHLY:
921  compute_times_monthly(result);
922  break;
923  case YEARLY:
924  compute_times_yearly(result);
925  break;
926  }
927 }
928 
929 /*!
930  * Check if the modeled time interval is a subset of the time interval in a forcing file.
931  *
932  * Returns silently if it is, otherwise throws an exception with an error message.
933  */
934 void check_forcing_duration(const Time &time,
935  double forcing_start,
936  double forcing_end) {
937 
938  double run_start = time.start();
939  double run_end = time.end();
940 
941  if (not (run_start >= forcing_start and
942  run_end <= forcing_end)) {
944  "A time-dependent forcing has to span the whole length of the simulation\n"
945  " Run time: [%s, %s]\n"
946  " Forcing data: [%s, %s]",
947  time.date(run_start).c_str(),
948  time.date(run_end).c_str(),
949  time.date(forcing_start).c_str(),
950  time.date(forcing_end).c_str());
951  }
952 }
953 
954 } // end of namespace pism
calcalcs_cal * ccs_init_calendar(const char *calname)
Definition: calcalcs.c:103
char * ccs_err_str(int error_code)
Definition: calcalcs.c:907
int ccs_isleap(calcalcs_cal *calendar, int year, int *leap)
Definition: calcalcs.c:335
void ccs_free_calendar(calcalcs_cal *cc)
Definition: calcalcs.c:1337
int ccs_date2jday(calcalcs_cal *calendar, int year, int month, int day, int *jday)
Definition: calcalcs.c:441
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
std::shared_ptr< const Config > ConstPtr
std::string get_string(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
bool get_flag(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:361
std::string filename() const
Definition: File.cc:307
unsigned int dimension_length(const std::string &name) const
Get the length of a dimension.
Definition: File.cc:454
std::string read_text_attribute(const std::string &var_name, const std::string &att_name) const
Get a text attribute.
Definition: File.cc:693
High-level PISM I/O class.
Definition: File.hh:56
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
Definition: Logger.cc:49
A basic logging class.
Definition: Logger.hh:40
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
double increment_date(double T, double years) const
Definition: Time.cc:852
std::vector< double > parse_list(const std::string &spec) const
Definition: Time.cc:521
std::string run_length() const
Returns the length of the current run, in years.
Definition: Time.cc:501
const Config::ConstPtr m_config
Definition: Time.hh:151
std::string units_string() const
Internal time units as a string.
Definition: Time.cc:477
void compute_times_monthly(std::vector< double > &result) const
Definition: Time.cc:856
double years_to_seconds(double input) const
Definition: Time.cc:419
void set(double new_time)
Sets the current time (in seconds since the reference time).
Definition: Time.cc:452
double convert_time_interval(double T, const std::string &units) const
Convert time interval from seconds to given units. Handle 'years' using the year length corresponding...
Definition: Time.cc:674
double m_year_length
number of seconds in a year, for "year fraction"
Definition: Time.hh:159
void set_end(double new_end)
Definition: Time.cc:460
void compute_times_simple(double time_start, double delta, double time_end, std::vector< double > &result) const
Definition: Time.cc:650
void init_calendar(const std::string &calendar)
Definition: Time.cc:431
double start() const
Definition: Time.cc:469
double calendar_year_start(double T) const
Returns the model time in seconds corresponding to the beginning of the year T falls into.
Definition: Time.cc:843
void compute_times(double time_start, double time_end, const Interval &interval, std::vector< double > &result) const
Definition: Time.cc:913
@ MONTHLY
Definition: Time.hh:129
@ SIMPLE
Definition: Time.hh:129
@ YEARLY
Definition: Time.hh:129
double end() const
Definition: Time.cc:473
std::vector< double > parse_range(const std::string &spec) const
Definition: Time.cc:601
double m_run_end
run end tim, in seconds since the reference time
Definition: Time.hh:168
double year_fraction(double T) const
Returns the fraction of a year that passed since the last beginning of a year. Only useful in codes w...
Definition: Time.cc:820
double m_time_in_seconds
current time, in seconds since the reference time
Definition: Time.hh:162
Interval parse_interval_length(const std::string &spec) const
Definition: Time.cc:542
std::string m_calendar_string
CF calendar string.
Definition: Time.hh:171
bool m_simple_calendar
Definition: Time.hh:173
double seconds_to_years(double input) const
Definition: Time.cc:427
double current() const
Current time, in seconds.
Definition: Time.cc:465
units::Unit units() const
Definition: Time.cc:481
std::vector< double > parse_times(const std::string &spec) const
Definition: Time.cc:510
Time(MPI_Comm com, Config::ConstPtr config, const Logger &log, units::System::Ptr unit_system)
Definition: Time.cc:688
void compute_times_yearly(std::vector< double > &result) const
Definition: Time.cc:888
double day_of_the_year_to_year_fraction(unsigned int day) const
Convert the day number to the year fraction.
Definition: Time.cc:505
void step(double delta_t)
Advance by delta_t seconds.
Definition: Time.cc:490
void set_start(double new_start)
Definition: Time.cc:456
std::string date(double T) const
Returns the date corresponding to time T.
Definition: Time.cc:834
double m_run_start
run start time, in seconds since the reference time
Definition: Time.hh:165
const units::System::Ptr m_unit_system
Definition: Time.hh:152
units::Unit m_time_units
Definition: Time.hh:153
std::string calendar() const
Returns the calendar string.
Definition: Time.cc:486
void init_from_file(MPI_Comm com, const std::string &filename, const Logger &log, bool set_start_time)
Sets the time from a NetCDF file with a time dimension (-time_file).
Definition: Time.cc:764
double m_t_eps
Time resolution, in seconds.
Definition: Time.hh:156
Time management class.
Definition: Time.hh:55
std::shared_ptr< System > Ptr
Definition: Units.hh:47
DateTime date(double T, const std::string &calendar) const
Definition: Units.cc:153
std::string format() const
Definition: Units.cc:141
System::Ptr system() const
Definition: Units.cc:149
double time(const DateTime &d, const std::string &calendar) const
Definition: Units.cc:168
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
void read_timeseries(const File &file, const VariableMetadata &metadata, const Logger &log, std::vector< double > &data)
Read a time-series variable from a NetCDF file to a vector of doubles.
Definition: io_helpers.cc:860
void read_time_bounds(const File &file, const VariableMetadata &metadata, const Logger &log, std::vector< double > &data)
Definition: io_helpers.cc:974
static double D2(double u_x, double u_y, double u_z, double v_x, double v_y, double v_z)
bool are_convertible(const Unit &u1, const Unit &u2)
Definition: Units.cc:242
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:70
static double increment_date(const units::Unit &time_units, const std::string &calendar, double T, double years)
Definition: Time.cc:182
static double start_time(const Config &config, const Logger &log, const File *file, const std::string &reference_date, const std::string &calendar, const units::Unit &time_units)
Definition: Time.cc:348
std::string printf(const char *format,...)
static const double k
Definition: exactTestP.cc:42
static double end_time(const Config &config, double time_start, const std::string &calendar, const units::Unit &time_units)
Definition: Time.cc:399
std::string string_strip(const std::string &input)
void check_forcing_duration(const Time &time, double forcing_start, double forcing_end)
Definition: Time.cc:934
static std::string calendar_from_file(const File &file, const std::string &time_name, const std::string &default_value, bool stop_on_error)
Get the calendar name from a file.
Definition: Time.cc:79
static std::string reference_date(const File *input_file, const Config &config, const Logger &log)
Definition: Time.cc:109
static std::string calendar(const File *input_file, const Config &config, const Logger &log)
Definition: Time.cc:146
bool member(const std::string &string, const std::set< std::string > &set)
static double parse_date(const std::string &input, const units::Unit &time_units, const std::string &calendar)
Definition: Time.cc:251
long int parse_integer(const std::string &input)
static std::string reference_date_from_file(const File &file, const std::string &time_name, const std::string &default_value, bool stop_on_error)
Get the reference date from a file.
Definition: Time.cc:40
bool pism_is_valid_calendar_name(const std::string &name)
Definition: Time.hh:34
std::vector< std::string > split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a vector of strings.
const int I[]
Definition: ssafd_code.cc:24
IntervalType type
Definition: Time.hh:133