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