tTime.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.cc (18010B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2019 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 <sstream>
21
22 #include "Time.hh"
23
24 #include "ConfigInterface.hh"
25 #include "Time_Calendar.hh"
26 #include "pism_options.hh"
27 #include "pism_utilities.hh"
28 #include "error_handling.hh"
29 #include "pism/util/io/File.hh"
30 #include "pism/util/Logger.hh"
31
32 namespace pism {
33
34 /**
35 * Select a calendar using the "time.calendar" configuration parameter, the
36 * "-calendar" command-line option, or the "calendar" attribute of the
37 * "time" variable in the file specified using "-time_file".
38 *
39 */
40 std::string calendar_from_options(MPI_Comm com, const Config &config) {
41 // Set the default calendar using the config. parameter or the
42 // "-calendar" option:
43 std::string result = config.get_string("time.calendar");
44
45 // Check if -time_file was set and override the setting above if the
46 // "calendar" attribute is found.
47 options::String time_file("-time_file", "name of the file specifying the run duration");
48 if (time_file.is_set()) {
49 File file(com, time_file, PISM_NETCDF3, PISM_READONLY); // OK to use netcdf3
50
51 std::string time_name = config.get_string("time.dimension_name");
52 if (file.find_variable(time_name)) {
53 std::string tmp = file.read_text_attribute(time_name, "calendar");
54 if (not tmp.empty()) {
55 result = tmp;
56 }
57 }
58 }
59 return result;
60 }
61
62 Time::Ptr time_from_options(MPI_Comm com, Config::ConstPtr config, units::System::Ptr system) {
63 try {
64 std::string calendar = calendar_from_options(com, *config);
65
66 if (calendar == "360_day" or
67 calendar == "365_day" or
68 calendar == "noleap" or
69 calendar == "none") {
70 return Time::Ptr(new Time(config, calendar, system));
71 } else {
72 return Time::Ptr(new Time_Calendar(com, config, calendar, system));
73 }
74
75 } catch (RuntimeError &e) {
76 e.add_context("initializing Time from options");
77 throw;
78 }
79 }
80
81 //! Initialize model time using command-line options and (possibly) files.
82 void initialize_time(MPI_Comm com, const std::string &dimension_name,
83 const Logger &log, Time &time) {
84
85 // Check if we are initializing from a PISM output file:
86 options::String input_file("-i", "Specifies a PISM input file");
87
88 if (input_file.is_set()) {
89 File file(com, input_file, PISM_NETCDF3, PISM_READONLY); // OK to use netcdf3
90 time.init_from_input_file(file, dimension_name, log);
91 }
92
93 time.init(log);
94 }
95
96 //! Get the reference date from a file.
97 std::string reference_date_from_file(const File &file,
98 const std::string &time_name) {
99
100 if (not file.find_variable(time_name)) {
101 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "'%s' variable is not present in '%s'.",
102 time_name.c_str(), file.filename().c_str());
103 }
104 std::string time_units = file.read_text_attribute(time_name, "units");
105
106 // Check if the time_units includes a reference date.
107 size_t position = time_units.find("since");
108 if (position == std::string::npos) {
109 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "time units string '%s' does not contain a reference date",
110 time_units.c_str());
111 }
112
113 return time_units.substr(position);
114 }
115
116
117 //! Convert model years into seconds using the year length
118 //! corresponding to the current calendar.
119 /*! Do not use this to convert quantities other than time intervals!
120 */
121 double Time::years_to_seconds(double input) const {
122 return input * m_year_length;
123 }
124
125 //! Convert seconds into model years using the year length
126 //! corresponding to the current calendar.
127 /*! Do not use this to convert quantities other than time intervals!
128 */
129 double Time::seconds_to_years(double input) const {
130 return input / m_year_length;
131 }
132
133
134 Time::Time(Config::ConstPtr conf,
135 const std::string &calendar_string,
136 units::System::Ptr unit_system)
137 : m_config(conf),
138 m_unit_system(unit_system),
139 m_time_units(m_unit_system, "seconds") {
140
141 init_calendar(calendar_string);
142
143 m_run_start = years_to_seconds(m_config->get_number("time.start_year"));
144 m_run_end = increment_date(m_run_start, (int)m_config->get_number("time.run_length"));
145
146 m_time_in_seconds = m_run_start;
147 }
148
149 Time::~Time() {
150 }
151
152 void Time::init_calendar(const std::string &calendar_string) {
153
154 if (not pism_is_valid_calendar_name(calendar_string)) {
155 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "unsupported calendar: %s", calendar_string.c_str());
156 }
157
158 m_calendar_string = calendar_string;
159
160 double seconds_per_day = convert(m_unit_system, 1.0, "day", "seconds");
161 if (calendar_string == "360_day") {
162 m_year_length = 360 * seconds_per_day;
163 } else if (calendar_string == "365_day" || calendar_string == "noleap") {
164 m_year_length = 365 * seconds_per_day;
165 } else {
166 // use the ~365.2524-day year
167 m_year_length = convert(m_unit_system, 1.0, "year", "seconds");
168 }
169 }
170
171 //! \brief Sets the current time (in seconds since the reference time).
172 void Time::set(double new_time) {
173 m_time_in_seconds = new_time;
174 }
175
176 void Time::set_start(double new_start) {
177 m_run_start = new_start;
178 }
179
180 void Time::set_end(double new_end) {
181 m_run_end = new_end;
182 }
183
184 //! \brief Current time, in seconds.
185 double Time::current() const {
186 return m_time_in_seconds;
187 }
188
189 double Time::start() const {
190 return m_run_start;
191 }
192
193 double Time::end() const {
194 return m_run_end;
195 }
196
197 std::string Time::CF_units_string() const {
198 return "seconds since " + m_config->get_string("time.reference_date");
199 }
200
201 //! \brief Returns the calendar string.
202 std::string Time::calendar() const {
203 return m_calendar_string;
204 }
205
206 void Time::step(double delta_t) {
207 m_time_in_seconds += delta_t;
208
209 // If we are less than 0.001 second from the end of the run, reset
210 // m_time_in_seconds to avoid taking a very small (and useless) time step.
211 if (m_run_end > m_time_in_seconds &&
212 m_run_end - m_time_in_seconds < 1e-3) {
213 m_time_in_seconds = m_run_end;
214 }
215 }
216
217 std::string Time::units_string() const {
218 return "seconds";
219 }
220
221
222 std::string Time::CF_units_to_PISM_units(const std::string &input) const {
223 std::string units = input;
224 size_t n = units.find("since");
225
226 /*!
227 \note This code finds the string "since" in the units_string and
228 terminates it on the first 's' of "since", if this sub-string was found.
229 This is done to ignore the reference date in the time units string (the
230 reference date specification always starts with this word).
231 */
232 if (n != std::string::npos) {
233 units.resize(n);
234 }
235
236 // strip trailing spaces
237 while (ends_with(units, " ") && not units.empty()) {
238 units.resize(units.size() - 1); // this would fail on empty strings
239 }
240
241 return units;
242 }
243
244 bool Time::process_ys(double &result) {
245 options::Real ys("-ys", "Start year", m_config->get_number("time.start_year"));
246 result = years_to_seconds(ys);
247 return ys.is_set();
248 }
249
250 bool Time::process_y(double &result) {
251 options::Real y("-y", "Run length, in years", m_config->get_number("time.run_length"));
252 result = years_to_seconds(y);
253 return y.is_set();
254 }
255
256 bool Time::process_ye(double &result) {
257 options::Real ye("-ye", "End year",
258 m_config->get_number("time.start_year") + m_config->get_number("time.run_length"));
259 result = years_to_seconds(ye);
260 return ye.is_set();
261 }
262
263
264 //! Set start time from a PISM input file.
265 /**
266 * FIXME: This crude implementation does not use reference dates and does not convert units.
267 */
268 void Time::init_from_input_file(const File &file,
269 const std::string &time_name,
270 const Logger &log) {
271 unsigned int time_length = file.dimension_length(time_name);
272
273 bool ys = options::Bool("-ys", "starting time");
274 if (not ys and time_length > 0) {
275 // Set the default starting time to be equal to the last time saved in the input file
276 double T = vector_max(file.read_dimension(time_name));
277 this->set_start(T);
278 this->set(T);
279 log.message(2,
280 "* Time t = %s found in '%s'; setting current time\n",
281 this->date().c_str(), file.filename().c_str());
282 }
283 }
284
285
286 void Time::init(const Logger &log) {
287
288 (void) log;
289
290 double y_seconds, ys_seconds, ye_seconds;
291
292 // At this point the calendar and the year length are set (in the
293 // constructor). The Time_Calendar class will (potentially)
294 // override all this by using settings from -time_file, so that is
295 // fine, too.
296
297 bool y_set = process_y(y_seconds);
298 bool ys_set = process_ys(ys_seconds);
299 bool ye_set = process_ye(ye_seconds);
300
301 if (ys_set and ye_set and y_set) {
302 throw RuntimeError(PISM_ERROR_LOCATION, "all of -y, -ys, -ye are set.");
303 }
304
305 if (y_set and ye_set) {
306 throw RuntimeError(PISM_ERROR_LOCATION, "using -y and -ye together is not allowed.");
307 }
308
309 // Set the start year if -ys is set, use the default otherwise.
310 if (ys_set) {
311 m_run_start = ys_seconds;
312 }
313
314 m_time_in_seconds = m_run_start;
315
316 if (ye_set) {
317 if (ye_seconds < m_time_in_seconds) {
318 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "-ye (%s) is less than -ys (%s) (or input file year or default).\n"
319 "PISM cannot run backward in time.",
320 date(ye_seconds).c_str(), date(m_run_start).c_str());
321 }
322 m_run_end = ye_seconds;
323 } else if (y_set) {
324 m_run_end = m_run_start + y_seconds;
325 } else {
326 m_run_end = increment_date(m_run_start, (int)m_config->get_number("time.run_length"));
327 }
328 }
329
330 std::string Time::date(double T) const {
331 char tmp[256];
332 snprintf(tmp, 256, "%.3f", seconds_to_years(T));
333 return std::string(tmp);
334 }
335
336 std::string Time::date() const {
337 return date(m_time_in_seconds);
338 }
339
340 std::string Time::start_date() const {
341 return date(m_run_start);
342 }
343
344 std::string Time::end_date() const {
345 return date(m_run_end);
346 }
347
348 std::string Time::run_length() const {
349 char tmp[256];
350 snprintf(tmp, 256, "%3.3f", seconds_to_years(m_run_end - m_run_start));
351 return std::string(tmp);
352 }
353
354 double Time::mod(double time, unsigned int period_years) const {
355 if (period_years == 0) {
356 return time;
357 }
358
359 double period_seconds = years_to_seconds(period_years);
360
361 double tmp = time - floor(time / period_seconds) * period_seconds;
362
363 if (fabs(tmp - period_seconds) < 1) {
364 tmp = 0;
365 }
366
367 return tmp;
368 }
369
370 double Time::year_fraction(double T) const {
371 double Y = seconds_to_years(T);
372 return Y - floor(Y);
373 }
374
375 double Time::day_of_the_year_to_day_fraction(unsigned int day) const {
376 const double sperd = 86400.0;
377 return (sperd / m_year_length) * (double) day;
378 }
379
380 double Time::calendar_year_start(double T) const {
381 return T - this->mod(T, 1);
382 }
383
384 double Time::increment_date(double T, int years) const {
385 return T + years_to_seconds(years);
386 }
387
388 std::vector<double> Time::parse_times(const std::string &spec) const {
389 if (spec.find(',') != std::string::npos) {
390 // a list will always contain a comma because at least two numbers are
391 // needed to specify reporting intervals
392 return parse_list(spec);
393 } else {
394 // it must be a range specification
395 return parse_range(spec);
396 }
397 }
398
399 std::vector<double> Time::parse_list(const std::string &spec) const {
400 std::istringstream arg(spec);
401 std::string tmp;
402 std::vector<double> result;
403
404 while(getline(arg, tmp, ',')) {
405 try {
406 double d;
407 parse_date(tmp, &d);
408 result.push_back(d);
409 } catch (RuntimeError &e) {
410 e.add_context("parsing a list of dates %s", spec.c_str());
411 throw;
412 }
413 }
414 return result;
415 }
416
417 /**
418 * Parses an interval specification string.
419 *
420 * @param[in] spec specification string
421 * @param[out] keyword interval type keyword, one of "hourly",
422 * "daily", "monthly", "yearly", "equal"
423 * @param[out] result if `keyword` == "equal", `result` is set to
424 * the interval length in seconds
425 *
426 * @return 0 on success, 1 otherwise
427 */
428 void Time::parse_interval_length(const std::string &spec, std::string &keyword, double *result) const {
429
430 // check if it is a keyword
431 if (spec == "hourly") {
432 keyword = "simple";
433 if (result) {
434 *result = 3600;
435 }
436 return;
437 }
438
439 if (spec == "daily") {
440 keyword = "simple";
441 if (result) {
442 *result = 86400;
443 }
444 return;
445 }
446
447 if (spec == "monthly" || spec == "yearly") {
448 keyword = spec;
449 if (result) {
450 *result = 0;
451 }
452 return;
453 }
454
455 units::Unit seconds(m_time_units.get_system(), "seconds"),
456 one(m_time_units.get_system(), "1"),
457 tmp = one;
458
459 try {
460 tmp = units::Unit(m_time_units.get_system(), spec);
461 } catch (RuntimeError &e) {
462 e.add_context("processing interval length " + spec);
463 throw;
464 }
465
466 // Check if these units are compatible with "seconds" or "1". The
467 // latter allows intervals of the form "0.5", which stands for "half
468 // of a model year". This also discards interval specs such as "days
469 // since 1-1-1", even though "days" is compatible with "seconds".
470 if (units::are_convertible(tmp, seconds)) {
471 units::Converter c(tmp, seconds);
472
473 if (result) {
474 *result = c(1.0);
475 }
476
477 } else if (units::are_convertible(tmp, one)) {
478 units::Converter c(tmp, one);
479
480 if (result) {
481 // this is a rather convoluted way of turning a string into a
482 // floating point number:
483 *result = c(1.0);
484 // convert from years to seconds without using UDUNITS-2 (this
485 // way we handle 360-day and 365-day years correctly)
486 *result = years_to_seconds(*result);
487 }
488
489 } else {
490 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "interval length '%s' is invalid", spec.c_str());
491 }
492 }
493
494
495 std::vector<double> Time::parse_range(const std::string &spec) const {
496 double
497 time_start = m_run_start,
498 time_end = m_run_end,
499 delta = 0;
500 std::string keyword = "simple";
501
502 if (spec == "hourly") {
503 delta = 3600;
504 } else if (spec == "daily") {
505 delta = 86400;
506 } else if (spec == "monthly" || spec == "yearly") {
507 keyword = spec;
508 delta = 0;
509 } else {
510
511 std::vector<std::string> parts = pism::split(spec, ':');
512
513 if (parts.size() == 1) {
514 parse_interval_length(parts[0], keyword, &delta);
515 } else if (parts.size() == 3) {
516 parse_date(parts[0], &time_start);
517 parse_interval_length(parts[1], keyword, &delta);
518 parse_date(parts[2], &time_end);
519 } else {
520 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "a time range must consist of exactly 3 parts separated by colons (got '%s').",
521 spec.c_str());
522 }
523 }
524
525 std::vector<double> result;
526 compute_times(time_start, delta, time_end, keyword, result);
527 return result;
528 }
529
530
531 void Time::parse_date(const std::string &spec, double *result) const {
532
533 if (spec.empty()) {
534 throw RuntimeError(PISM_ERROR_LOCATION, "got an empty date specification");
535 }
536
537 char *endptr = NULL;
538 double d = strtod(spec.c_str(), &endptr);
539 if (*endptr != '\0') {
540 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "date specification '%s' is invalid ('%s' is not an number)",
541 spec.c_str(), spec.c_str());
542 }
543
544 if (result) {
545 *result = years_to_seconds(d);
546 }
547 }
548
549
550 /**
551 * Compute times corresponding to a "simple" time range.
552 *
553 * This is a time range with a fixed step ("hourly" is an example).
554 *
555 * @param[in] time_start beginning of the time interval, in seconds
556 * @param[in] delta step, in seconds
557 * @param[in] time_end end of the interval, in seconds
558 * @param[out] result list of model times
559 *
560 */
561 void Time::compute_times_simple(double time_start, double delta, double time_end,
562 std::vector<double> &result) const {
563 if (time_start >= time_end) {
564 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "a >= b in time range a:dt:b (got a = %s, b = %s)",
565 this->date(time_start).c_str(), this->date(time_end).c_str());
566 }
567
568 if (delta <= 0) {
569 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "dt <= 0 in time range a:dt:b (got dt = %f seconds)", delta);
570 }
571
572 int k = 0;
573 double t = time_start;
574
575 result.clear();
576 do {
577 if (t >= this->start() && t <= this->end()) {
578 result.push_back(t);
579 }
580 k += 1;
581 t = time_start + k * delta;
582 } while (t <= time_end);
583 }
584
585 void Time::compute_times(double time_start, double delta, double time_end,
586 const std::string &keyword,
587 std::vector<double> &result) const {
588 if (keyword == "yearly") {
589 delta = years_to_seconds(1.0);
590 } else if (keyword == "monthly") {
591 delta = years_to_seconds(1.0/12.0);
592 } else if (keyword != "simple") {
593 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "unknown time range keyword: %s",
594 keyword.c_str());
595 }
596
597 compute_times_simple(time_start, delta, time_end, result);
598 }
599
600 double Time::convert_time_interval(double T, const std::string &units) const {
601 if (units == "year" || units == "years" || units == "yr" || units == "a") {
602 return this->seconds_to_years(T); // uses year length here
603 }
604 return convert(m_unit_system, T, "seconds", units);
605 }
606
607 double Time::current_years() const {
608 return seconds_to_years(current());
609 }
610
611 } // end of namespace pism