1 /*********************************************************************
2  *   Copyright 2008, University Corporation for Atmospheric Research
3  *   See netcdf/COPYRIGHT file for copying and redistribution conditions.
4  *   $Id: nctime.c,v 1.9 2010/05/05 22:15:39 dmh Exp $
5  *********************************************************************/
6 
7 /*
8  * This code was extracted with permission from the CDMS time
9  * conversion and arithmetic routines developed by Bob Drach, Lawrence
10  * Livermore National Laboratory as part of the cdtime library.  Russ
11  * Rew of the UCAR Unidata Program made changes and additions to
12  * support the "-t" option of the netCDF ncdump utility, including a
13  * 366-day climate calendar.
14  *
15  * For the complete time conversion and climate calendar facilities of
16  * the CDMS library, get the original sources from LLNL.
17  */
18 
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include <ctype.h>
22 #include <math.h>
23 #include <string.h>
24 #include <stdarg.h>
25 #include <assert.h>
26 #include "nctime.h"
27 
28 static int cuErrOpts;			     /* Error options */
29 static int cuErrorOccurred = 0;		     /* True iff cdError was called */
30 
31 #define CD_DEFAULT_BASEYEAR "1979"	     /* Default base year for relative time (no 'since' clause) */
32 #define VALCMP(a,b) ((a)<(b)?-1:(b)<(a)?1:0)
33 
34 /* forward declarations */
35 static void CdMonthDay(int *doy, CdTime *date);
36 static void CdDayOfYear(CdTime *date, int *doy);
37 static void cdComp2Rel(cdCalenType timetype, cdCompTime comptime, char* relunits, double* reltime);
38 static void cdRel2CompMixed(double reltime, cdUnitTime unit, cdCompTime basetime, cdCompTime *comptime);
39 static void cdRel2Comp(cdCalenType timetype, char* relunits, double reltime, cdCompTime* comptime);
40 
41 /* Trim trailing whitespace, up to n characters. */
42 /* If no whitespace up to the last character, set */
43 /* the last character to null, else set the first */
44 /* whitespace character to null. */
45 static void
cdTrim(char * s,int n)46 cdTrim(char* s, int n)
47 {
48 	char* c;
49 
50 	if(s==NULL)
51 		return;
52 	for(c=s; *c && c<s+n-1 && !isspace((int)*c); c++);
53 	*c='\0';
54 	return;
55 }
56 
57 static void
cdError(char * fmt,...)58 cdError(char *fmt, ...)
59 {
60 	va_list args;
61 
62 	cuErrorOccurred = 1;
63 	if(cuErrOpts & CU_VERBOSE){
64 		va_start(args,fmt);
65 		fprintf(stderr, "CDMS error: ");
66 		vfprintf(stderr, fmt, args);
67 		fprintf(stderr, "\n");
68 		va_end(args);
69 	}
70 	if(cuErrOpts & CU_FATAL)
71 		exit(1);
72 	return;
73 }
74 
75 #define ISLEAP(year,timeType)	((timeType & Cd366) || (((timeType) & CdHasLeap) && (!((year) % 4) && (((timeType) & CdJulianType) || (((year) % 100) || !((year) % 400))))))
76 
77 static int mon_day_cnt[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
78 static int days_sum[12] = {0,31,59,90,120,151,181,212,243,273,304,334};
79 
80 /* Compute month and day from year and day-of-year.
81  *
82  *	Input:
83  *		doy	     (int)  (day-of-year)
84  *		date->year   (long)  (year since 0 BC)
85  *              date->timeType (CdTimetype) (time type)
86  *              date->baseYear   base year for relative times
87  *	Output:
88  *		date->month  (short)  (month in year)
89  *		date->day    (short)  (day in month)
90  *
91  *
92  * Derived from NRL NEONS V3.6.
93  */
94 
95 static void
CdMonthDay(int * doy,CdTime * date)96 CdMonthDay(int *doy, CdTime *date)
97 {
98 	int i;				/* month counter */
99 	int idoy;			/* day of year counter */
100 	long year;
101 
102 	if ((idoy = *doy) < 1) {
103 		date->month = 0;
104 		date->day   = 0;
105 		return;
106 	}
107 
108 	if(!(date->timeType & CdChronCal))   /* Ignore year for Clim calendar */
109 		year = 0;
110 	else if(!(date->timeType & CdBase1970))	/* year is offset from base for relative time */
111 		year = date->baseYear + date->year;
112 	else
113 		year = date->year;
114 
115 	if (ISLEAP(year,date->timeType)) {
116 		mon_day_cnt[1] = 29;
117 	} else {
118 		mon_day_cnt[1] = 28;
119 	}
120 	date->month	= 0;
121 	for (i = 0; i < 12; i++) {
122 		int delta;
123 		(date->month)++;
124 		date->day = (short)idoy;
125 		delta = ((date->timeType & Cd365) || (date->timeType & Cd366) ? (mon_day_cnt[date->month-1]) : 30);
126 	        idoy -= delta;
127 		if(idoy <= 0)
128 		    return;
129 	}
130 	return;
131 }
132 
133 /* Compute day-of-year from year, month and day
134  *
135  *	Input:
136  *		date->year  (long)  (year since 0 BC)
137  *		date->month (short)  (month in year)
138  *		date->day   (short)  (day in month)
139  *              date->baseYear   base year for relative times
140  *	Output: doy         (int)  (day-of-year)
141  *
142  * Derived from NRL NEONS V3.6
143  */
144 
145 static void
CdDayOfYear(CdTime * date,int * doy)146 CdDayOfYear(CdTime *date, int *doy)
147 {
148 	int leap_add = 0;		/* add 1 day if leap year */
149 	int month;			/* month */
150 	long year;
151 
152    	month	= date->month;
153 	if (month < 1 || month > 12) {
154 		cdError( "Day-of-year error; month: %d\n", month);
155 		month = 1;
156 	}
157 
158 	if(!(date->timeType & CdChronCal))   /* Ignore year for Clim calendar */
159 		year = 0;
160 	else if(!(date->timeType & CdBase1970))	/* year is offset from base for relative time */
161 		year = date->baseYear + date->year;
162 	else
163 		year = date->year;
164 
165 	if (ISLEAP(year,date->timeType) && month > 2) leap_add = 1;
166 	if( ((date->timeType) & Cd365) || ((date->timeType) & Cd366) ) {
167 	    *doy = days_sum[month-1] + date->day + leap_add ;
168 	} else {		/* date->timeType & Cd360 */
169 	    *doy = 30*(month-1) + date->day + leap_add ;
170 	}
171 	return;
172 }
173 
174 /* Convert epochal time (hours since 00 jan 1, 1970)
175  *   to human time (structured)
176  *
177  * Input:
178  *   etime = epochal time representation
179  *   timeType = time type (e.g., CdChron, CdClim, etc.) as defined in cdms.h
180  *   baseYear = base real, used for relative time types only
181  *
182  * Output: htime = human (structured) time representation
183  *
184  * Derived from NRL Neons V3.6
185  */
186 void
Cde2h(double etime,CdTimeType timeType,long baseYear,CdTime * htime)187 Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime)
188 {
189 	long 	ytemp;			/* temporary year holder */
190 	int 	yr_day_cnt;		/* count of days in year */
191 	int 	doy;			/* day of year */
192 	int     daysInLeapYear;		     /* number of days in a leap year */
193 	int     daysInYear;		     /* days in non-leap year */
194 
195 	doy	= (int) floor(etime / 24.) + 1;
196 	htime->hour	= etime - (double) (doy - 1) * 24.;
197 
198 					     /* Correct for goofy floor func on J90 */
199 	if(htime->hour >= 24.){
200 		doy += 1;
201 		htime->hour -= 24.;
202 	}
203 
204 	htime->baseYear = (timeType & CdBase1970) ? 1970 : baseYear;
205 	if(!(timeType & CdChronCal)) htime->baseYear = 0; /* Set base year to 0 for Clim */
206 	if(timeType & Cd366) {
207 	    daysInLeapYear = 366;
208 	    daysInYear = 366;
209 	} else {
210 	    daysInLeapYear = (timeType & Cd365) ? 366 : 360;
211 	    daysInYear = (timeType & Cd365) ? 365 : 360;
212 	}
213 
214 	if (doy > 0) {
215 		for (ytemp = htime->baseYear; ; ytemp++) {
216 			yr_day_cnt = ISLEAP(ytemp,timeType) ? daysInLeapYear : daysInYear;
217 			if (doy <= yr_day_cnt) break;
218 			doy -= yr_day_cnt;
219 		}
220 	} else {
221 		for (ytemp = htime->baseYear-1; ; ytemp--) {
222 			yr_day_cnt = ISLEAP(ytemp,timeType) ? daysInLeapYear : daysInYear;
223 			doy += yr_day_cnt;
224 			if (doy > 0) break;
225 		}
226 	}
227         htime->year = (timeType & CdBase1970) ? ytemp : (ytemp - htime->baseYear);
228 	if(!(timeType & CdChronCal)) htime->year = 0; /* Set year to 0 for Clim */
229 	htime->timeType = timeType;
230 	CdMonthDay(&doy,htime);
231 
232         return;
233 }
234 
235 /* Add 'nDel' times 'delTime' to epochal time 'begEtm',
236  * return the result in epochal time 'endEtm'.
237  */
238 static void
CdAddDelTime(double begEtm,long nDel,CdDeltaTime delTime,CdTimeType timeType,long baseYear,double * endEtm)239 CdAddDelTime(double begEtm, long nDel, CdDeltaTime delTime, CdTimeType timeType,
240 	     long baseYear, double *endEtm)
241 {
242 	double delHours;
243 	long delMonths, delYears;
244 	CdTime bhtime, ehtime;
245 
246 	extern void Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime);
247 	extern void Cdh2e(CdTime *htime, double *etime);
248 
249 	switch(delTime.units){
250 	  case CdYear:
251 		delMonths = 12;
252 		break;
253 	  case CdSeason:
254 		delMonths = 3;
255 		break;
256 	  case CdMonth:
257 		delMonths = 1;
258 		break;
259 	  case CdWeek:
260 		delHours = 168.0;
261 		break;
262 	  case CdDay:
263 		delHours = 24.0;
264 		break;
265 	  case CdHour:
266 		delHours = 1.0;
267 		break;
268 	  case CdMinute:
269 		delHours = 1./60.;
270 		break;
271 	  case CdSecond:
272 		delHours = 1./3600.;
273 		break;
274 	  default:
275 		cdError("Invalid delta time units: %d\n",delTime.units);
276 		return;
277 	}
278 
279 	switch(delTime.units){
280 	  case CdYear: case CdSeason: case CdMonth:
281 		Cde2h(begEtm,timeType,baseYear,&bhtime);
282 		delMonths = delMonths * nDel * delTime.count + bhtime.month - 1;
283 		delYears = (delMonths >= 0 ? (delMonths/12) : (delMonths+1)/12 - 1);
284 		ehtime.year = bhtime.year + delYears;
285 		ehtime.month = (short)(delMonths - (12 * delYears) + 1);
286 		ehtime.day = 1;
287 		ehtime.hour = 0.0;
288 		ehtime.timeType = timeType;
289 		ehtime.baseYear = !(timeType & CdChronCal) ? 0 :
290 			(timeType & CdBase1970) ? 1970 : baseYear; /* base year is 0 for Clim, */
291 								   /* 1970 for Chron, */
292 								   /* or input base year for Rel */
293 		Cdh2e(&ehtime,endEtm);
294 		break;
295 	  case CdWeek: case CdDay: case CdHour: case CdMinute: case CdSecond:
296 		delHours = delHours * (double)(nDel * delTime.count);
297 		*endEtm = begEtm + delHours;
298 		break;
299 	  default: break;
300 	}
301 	return;
302 }
303 
304 /* Parse relative units, returning the unit and base component time. */
305 /* Function returns 1 if error, 0 on success */
306 int
cdParseRelunits(cdCalenType timetype,char * relunits,cdUnitTime * unit,cdCompTime * base_comptime)307 cdParseRelunits(cdCalenType timetype, char* relunits, cdUnitTime* unit, cdCompTime* base_comptime)
308 {
309 	char charunits[CD_MAX_RELUNITS];
310 	char basetime_1[CD_MAX_CHARTIME];
311 	char basetime_2[CD_MAX_CHARTIME];
312 	char basetime[CD_MAX_CHARTIME];
313 	int nconv1, nconv2, nconv;
314 
315 					     /* Parse the relunits */
316 	/* Allow ISO-8601 "T" date-time separator as well as blank separator */
317 	nconv1 = sscanf(relunits,"%s since %[^T]T%s",charunits,basetime_1,basetime_2);
318 	if(nconv1==EOF || nconv1==0){
319 		cdError("Error on relative units conversion, string = %s\n",relunits);
320 		return 1;
321 	}
322 	nconv2 = sscanf(relunits,"%s since %s %s",charunits,basetime_1,basetime_2);
323 	if(nconv2==EOF || nconv2==0){
324 		cdError("Error on relative units conversion, string = %s\n",relunits);
325 		return 1;
326 	}
327 	if(nconv1 < nconv2) {
328 	    nconv = nconv2;
329 	} else {
330 	    nconv = sscanf(relunits,"%s since %[^T]T%s",charunits,basetime_1,basetime_2);
331 	}
332 					     /* Get the units */
333 	cdTrim(charunits,CD_MAX_RELUNITS);
334 	if(!strncmp(charunits,"sec",3) || !strcmp(charunits,"s")){
335 		*unit = cdSecond;
336 	}
337 	else if(!strncmp(charunits,"min",3) || !strcmp(charunits,"mn")){
338 		*unit = cdMinute;
339 	}
340 	else if(!strncmp(charunits,"hour",4) || !strcmp(charunits,"hr")){
341 		*unit = cdHour;
342 	}
343 	else if(!strncmp(charunits,"day",3) || !strcmp(charunits,"dy")){
344 		*unit = cdDay;
345 	}
346 	else if(!strncmp(charunits,"week",4) || !strcmp(charunits,"wk")){
347 		*unit = cdWeek;
348 	}
349 	else if(!strncmp(charunits,"month",5) || !strcmp(charunits,"mo")){
350 		*unit = cdMonth;
351 	}
352 	else if(!strncmp(charunits,"season",6)){
353 		*unit = cdSeason;
354 	}
355 	else if(!strncmp(charunits,"year",4) || !strcmp(charunits,"yr")){
356 		if(!(timetype & cdStandardCal)){
357 			cdError("Error on relative units conversion: climatological units cannot be 'years'.\n");
358 			return 1;
359 		}
360 		*unit = cdYear;
361 	}
362 	else {
363 		cdError("Error on relative units conversion: invalid units = %s\n",charunits);
364 		return 1;
365 	}
366 
367 					     /* Build the basetime, if any (default is 1979), */
368 					     /* or month 1 for climatological time. */
369 	if(nconv == 1){
370 		if(timetype & cdStandardCal)
371 			strcpy(basetime,CD_DEFAULT_BASEYEAR);
372 		else
373 			strcpy(basetime,"1");
374 	}
375 					     /* Convert the basetime to component, then epochal (hours since 1970) */
376 	else{
377 		if(nconv == 2){
378 			cdTrim(basetime_1,CD_MAX_CHARTIME);
379 			strcpy(basetime,basetime_1);
380 		}
381 		else{
382 			cdTrim(basetime_1,CD_MAX_CHARTIME);
383 			cdTrim(basetime_2,CD_MAX_CHARTIME);
384 			sprintf(basetime,"%s %s",basetime_1,basetime_2);
385 		}
386 	}
387 
388 	cdChar2Comp(timetype, basetime, base_comptime);
389 
390 	return 0;
391 }
392 
393 /* ca - cb in Gregorian calendar */
394 /* Result is in hours. */
395 static double
cdDiffGregorian(cdCompTime ca,cdCompTime cb)396 cdDiffGregorian(cdCompTime ca, cdCompTime cb){
397 
398 	double rela, relb;
399 
400 	cdComp2Rel(cdStandard, ca, "hours", &rela);
401 	cdComp2Rel(cdStandard, cb, "hours", &relb);
402 	return (rela - relb);
403 }
404 
405 /* Return -1, 0, 1 as ca is less than, equal to, */
406 /* or greater than cb, respectively. */
407 static int
cdCompCompare(cdCompTime ca,cdCompTime cb)408 cdCompCompare(cdCompTime ca, cdCompTime cb){
409 
410 	int test;
411 
412 	if ((test = VALCMP(ca.year, cb.year)))
413 		return test;
414 	else if ((test = VALCMP(ca.month, cb.month)))
415 		return test;
416 	else if ((test = VALCMP(ca.day, cb.day)))
417 		return test;
418 	else
419 		return (VALCMP(ca.hour, cb.hour));
420 }
421 
422 /* ca - cb in Julian calendar.  Result is in hours. */
423 static double
cdDiffJulian(cdCompTime ca,cdCompTime cb)424 cdDiffJulian(cdCompTime ca, cdCompTime cb){
425 
426 	double rela, relb;
427 
428 	cdComp2Rel(cdJulian, ca, "hours", &rela);
429 	cdComp2Rel(cdJulian, cb, "hours", &relb);
430 	return (rela - relb);
431 }
432 
433 /* ca - cb in mixed Julian/Gregorian calendar. */
434 /* Result is in hours. */
435 static double
cdDiffMixed(cdCompTime ca,cdCompTime cb)436 cdDiffMixed(cdCompTime ca, cdCompTime cb){
437 
438 	static cdCompTime ZA = {1582, 10, 5, 0.0};
439 	static cdCompTime ZB = {1582, 10, 15, 0.0};
440 	double result;
441 
442 	if (cdCompCompare(cb, ZB) == -1){
443 		if (cdCompCompare(ca, ZB) == -1) {
444 			result = cdDiffJulian(ca, cb);
445 		}
446 		else {
447 			result = cdDiffGregorian(ca, ZB) + cdDiffJulian(ZA, cb);
448 		}
449 	}
450 	else {
451 		if (cdCompCompare(ca, ZB) == -1){
452 			result = cdDiffJulian(ca, ZA) + cdDiffGregorian(ZB, cb);
453 		}
454 		else {
455 			result = cdDiffGregorian(ca, cb);
456 		}
457 	}
458 	return result;
459 }
460 
461 /* Divide ('endEtm' - 'begEtm') by 'delTime',
462  * return the integer portion of the result in 'nDel'.
463  */
464 static void
CdDivDelTime(double begEtm,double endEtm,CdDeltaTime delTime,CdTimeType timeType,long baseYear,long * nDel)465 CdDivDelTime(double begEtm, double endEtm, CdDeltaTime delTime, CdTimeType timeType,
466 	     long baseYear, long *nDel)
467 {
468 	double delHours, frange;
469 	long delMonths, range;
470 	CdTime bhtime, ehtime;
471 	int hoursInYear;
472 
473 	extern void Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime);
474 
475 	switch(delTime.units){
476 	  case CdYear:
477 		delMonths = 12;
478 		break;
479 	  case CdSeason:
480 		delMonths = 3;
481 		break;
482 	  case CdMonth:
483 		delMonths = 1;
484 		break;
485 	  case CdWeek:
486 		delHours = 168.0;
487 		break;
488 	  case CdDay:
489 		delHours = 24.0;
490 		break;
491 	  case CdHour:
492 		delHours = 1.0;
493 		break;
494 	  case CdMinute:
495 		delHours = 1./60.;
496 		break;
497 	  case CdSecond:
498 		delHours = 1./3600.;
499 		break;
500 	  default:
501 		cdError("Invalid delta time units: %d\n",delTime.units);
502 		return;
503 	}
504 
505 	switch(delTime.units){
506 	  case CdYear: case CdSeason: case CdMonth:
507 		delMonths *= delTime.count;
508 		Cde2h(begEtm,timeType,baseYear,&bhtime);
509 		Cde2h(endEtm,timeType,baseYear,&ehtime);
510 		if(timeType & CdChronCal){   /* Chron and Rel time */
511 			range = 12*(ehtime.year - bhtime.year)
512 				+ (ehtime.month - bhtime.month);
513 		}
514 		else{			     /* Clim time, ignore year */
515 			range = (ehtime.month - bhtime.month);
516 			if(range < 0) range += 12;
517 		}
518 		*nDel = abs((int)range)/delMonths;
519 		break;
520 	  case CdWeek: case CdDay: case CdHour: case CdMinute: case CdSecond:
521 		delHours *= (double)delTime.count;
522 		if(timeType & CdChronCal){   /* Chron and Rel time */
523 			frange = fabs(endEtm - begEtm);
524 		}
525 		else{			     /* Clim time, ignore year, but */
526 					     /* wraparound relative to hours-in-year*/
527 			frange = endEtm - begEtm;
528 			if(timeType & Cd366) {
529 			    hoursInYear = 8784;
530 			} else {
531 			    hoursInYear = (timeType & Cd365) ? 8760. : 8640.;
532 			}
533 					     /* Normalize frange to interval [0,hoursInYear) */
534 			if(frange < 0.0 || frange >= hoursInYear)
535 				frange -= hoursInYear * floor(frange/hoursInYear);
536 		}
537 		*nDel = (long)((frange + 1.e-10*delHours)/delHours);
538 		break;
539 	    default: break;
540 	}
541 	return;
542 }
543 
544 /* Value is in hours. Translate to units. */
545 static double
cdFromHours(double value,cdUnitTime unit)546 cdFromHours(double value, cdUnitTime unit){
547 	double result;
548 
549 	switch(unit){
550 	case cdSecond:
551 		result = value * 3600.0;
552 		break;
553 	case cdMinute:
554 		result = value * 60.0;
555 		break;
556 	case cdHour:
557 		result = value;
558 		break;
559 	case cdDay:
560 		result = value/24.0;
561 		break;
562 	case cdWeek:
563 		result = value/168.0;
564 		break;
565 	case cdMonth:
566 	case cdSeason:
567 	case cdYear:
568 	case cdFraction:
569 	default:
570     		cdError("Error on conversion from hours to vague unit");
571 		result = 0;
572 		break;
573 	}
574 	return result;
575 }
576 					     /* Map to old timetypes */
577 static int
cdToOldTimetype(cdCalenType newtype,CdTimeType * oldtype)578 cdToOldTimetype(cdCalenType newtype, CdTimeType* oldtype)
579 {
580 	switch(newtype){
581 	  case cdStandard:
582 		*oldtype = CdChron;
583 		break;
584 	  case cdJulian:
585 		*oldtype = CdJulianCal;
586 		break;
587 	  case cdNoLeap:
588 		*oldtype = CdChronNoLeap;
589 		break;
590 	  case cd360:
591 		*oldtype = CdChron360;
592 		break;
593 	  case cd366:
594 		*oldtype = CdChron366;
595 		break;
596 	  case cdClim:
597 		*oldtype = CdClim;
598 		break;
599 	  case cdClimLeap:
600 		*oldtype = CdClimLeap;
601 		break;
602 	  case cdClim360:
603 		*oldtype = CdClim360;
604 		break;
605 	  default:
606 		cdError("Error on relative units conversion, invalid timetype = %d",newtype);
607 		return 1;
608 	}
609 	return 0;
610 }
611 
612 /* Convert human time to epochal time (hours since 00 jan 1, 1970)
613  *
614  * Input: htime = human time representation
615  *
616  * Output: etime = epochal time representation
617  *
618  * Derived from NRL Neons V3.6
619  */
620 void
Cdh2e(CdTime * htime,double * etime)621 Cdh2e(CdTime *htime, double *etime)
622 {
623 	long 	ytemp, year;			/* temporary year holder */
624 	int	day_cnt;		/* count of days */
625 	int 	doy;			/* day of year */
626 	long    baseYear;		     /* base year for epochal time */
627 	int     daysInLeapYear;		     /* number of days in a leap year */
628 	int     daysInYear;		     /* days in non-leap year */
629 
630 	CdDayOfYear(htime,&doy);
631 
632 	day_cnt	= 0;
633 
634 	baseYear = ((htime->timeType) & CdBase1970) ? 1970 : htime->baseYear;
635 	year = ((htime->timeType) & CdBase1970) ? htime->year : (htime->year + htime->baseYear);
636 	if(!((htime->timeType) & CdChronCal)) baseYear = year = 0;	/* set year and baseYear to 0 for Clim */
637 	if((htime->timeType) & Cd366) {
638 	    daysInLeapYear = 366;
639 	    daysInYear = 366;
640 	} else {
641 	    daysInLeapYear = ((htime->timeType) & Cd365) ? 366 : 360;
642 	    daysInYear = ((htime->timeType) & Cd365) ? 365 : 360;
643 	}
644 
645 	if (year > baseYear) {
646 		for (ytemp = year - 1; ytemp >= baseYear; ytemp--) {
647 			day_cnt += ISLEAP(ytemp,htime->timeType) ? daysInLeapYear : daysInYear;
648 		}
649 	} else if (year < baseYear) {
650 		for (ytemp = year; ytemp < baseYear; ytemp++) {
651 			day_cnt -= ISLEAP(ytemp,htime->timeType) ? daysInLeapYear : daysInYear;
652 		}
653 	}
654 	*etime	= (double) (day_cnt + doy - 1) * 24. + htime->hour;
655         return;
656 }
657 
658 /* Validate the component time, return 0 if valid, 1 if not */
659 static int
cdValidateTime(cdCalenType timetype,cdCompTime comptime)660 cdValidateTime(cdCalenType timetype, cdCompTime comptime)
661 {
662 	if(comptime.month<1 || comptime.month>12){
663 		cdError("Error on time conversion: invalid month = %hd\n",comptime.month);
664 		return 1;
665 	}
666 	if(comptime.day<1 || comptime.day>31){
667 		cdError("Error on time conversion: invalid day = %hd\n",comptime.day);
668 		return 1;
669 	}
670 	if(comptime.hour<0.0 || comptime.hour>24.0){
671 		cdError("Error on time conversion: invalid hour = %lf\n",comptime.hour);
672 		return 1;
673 	}
674 	return 0;
675 }
676 
677 void
cdChar2Comp(cdCalenType timetype,char * chartime,cdCompTime * comptime)678 cdChar2Comp(cdCalenType timetype, char* chartime, cdCompTime* comptime)
679 {
680 	double sec;
681 	int ihr, imin, nconv;
682 	long year;
683 	short day;
684 	short month;
685 
686 	comptime->year = CD_NULL_YEAR;
687 	comptime->month = CD_NULL_MONTH;
688 	comptime->day = CD_NULL_DAY;
689 	comptime->hour = CD_NULL_HOUR;
690 
691 	if(timetype & cdStandardCal){
692 		nconv = sscanf(chartime,"%ld-%hd-%hd %d:%d:%lf",&year,&month,&day,&ihr,&imin,&sec);
693 		if(nconv==EOF || nconv==0){
694 			cdError("Error on character time conversion, string = %s\n",chartime);
695 			return;
696 		}
697 		if(nconv >= 1){
698 			comptime->year = year;
699 		}
700 		if(nconv >= 2){
701 			comptime->month = month;
702 		}
703 		if(nconv >= 3){
704 			comptime->day = day;
705 		}
706 		if(nconv >= 4){
707 			if(ihr<0 || ihr>23){
708 				cdError("Error on character time conversion: invalid hour = %d\n",ihr);
709 				return;
710 			}
711 			comptime->hour = (double)ihr;
712 		}
713 		if(nconv >= 5){
714 			if(imin<0 || imin>59){
715 				cdError("Error on character time conversion: invalid minute = %d\n",imin);
716 				return;
717 			}
718 			comptime->hour += (double)imin/60.;
719 		}
720 		if(nconv >= 6){
721 			if(sec<0.0 || sec>60.0){
722 				cdError("Error on character time conversion: invalid second = %lf\n",sec);
723 				return;
724 			}
725 			comptime->hour += sec/3600.;
726 		}
727 	}
728 	else{				     /* Climatological */
729 		nconv = sscanf(chartime,"%hd-%hd %d:%d:%lf",&month,&day,&ihr,&imin,&sec);
730 		if(nconv==EOF || nconv==0){
731 			cdError("Error on character time conversion, string = %s",chartime);
732 			return;
733 		}
734 		if(nconv >= 1){
735 			comptime->month = month;
736 		}
737 		if(nconv >= 2){
738 			comptime->day = day;
739 		}
740 		if(nconv >= 3){
741 			if(ihr<0 || ihr>23){
742 				cdError("Error on character time conversion: invalid hour = %d\n",ihr);
743 				return;
744 			}
745 			comptime->hour = (double)ihr;
746 		}
747 		if(nconv >= 4){
748 			if(imin<0 || imin>59){
749 				cdError("Error on character time conversion: invalid minute = %d\n",imin);
750 				return;
751 			}
752 			comptime->hour += (double)imin/60.;
753 		}
754 		if(nconv >= 5){
755 			if(sec<0.0 || sec>60.0){
756 				cdError("Error on character time conversion: invalid second = %lf\n",sec);
757 				return;
758 			}
759 			comptime->hour += sec/3600.;
760 		}
761 	}
762 	(void)cdValidateTime(timetype,*comptime);
763 	return;
764 }
765 
766 /* Convert ct to relunits (unit, basetime) */
767 /* in the mixed Julian/Gregorian calendar. */
768 /* unit is anything but year, season, month. unit and basetime are */
769 /* from the parsed relunits. Return result in reltime. */
770 static void
cdComp2RelMixed(cdCompTime ct,cdUnitTime unit,cdCompTime basetime,double * reltime)771 cdComp2RelMixed(cdCompTime ct, cdUnitTime unit, cdCompTime basetime, double *reltime){
772 
773 	double hourdiff;
774 
775 	hourdiff = cdDiffMixed(ct, basetime);
776 	*reltime = cdFromHours(hourdiff, unit);
777 	return;
778 }
779 
780 static void
cdComp2Rel(cdCalenType timetype,cdCompTime comptime,char * relunits,double * reltime)781 cdComp2Rel(cdCalenType timetype, cdCompTime comptime, char* relunits, double* reltime)
782 {
783 	cdCompTime base_comptime;
784 	CdDeltaTime deltime;
785 	CdTime humantime;
786 	CdTimeType old_timetype;
787 	cdUnitTime unit;
788 	double base_etm, etm, delta;
789 	long ndel, hoursInYear;
790 
791 					     /* Parse the relunits */
792 	if(cdParseRelunits(timetype, relunits, &unit, &base_comptime))
793 		return;
794 
795 					     /* Handle mixed Julian/Gregorian calendar */
796 	if (timetype == cdMixed){
797 		switch(unit){
798 		case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
799 			cdComp2RelMixed(comptime, unit, base_comptime, reltime);
800 			return;
801 		case cdYear: case cdSeason: case cdMonth:
802 			timetype = cdStandard;
803 			break;
804 		case cdFraction:
805 		        cdError("invalid unit in conversion");
806 		        break;
807 		default: break;
808 		}
809 	}
810 
811 					     /* Convert basetime to epochal */
812 	humantime.year = base_comptime.year;
813 	humantime.month = base_comptime.month;
814 	humantime.day = base_comptime.day;
815 	humantime.hour = base_comptime.hour;
816 	humantime.baseYear = 1970;
817 					     /* Map to old-style timetype */
818 	if(cdToOldTimetype(timetype,&old_timetype))
819 		return;
820 	humantime.timeType = old_timetype;
821 	Cdh2e(&humantime,&base_etm);
822 
823 					     /* Map end time to epochal */
824 	humantime.year = comptime.year;
825 	humantime.month = comptime.month;
826 	humantime.day = comptime.day;
827 	humantime.hour = comptime.hour;
828 	Cdh2e(&humantime,&etm);
829 					     /* Calculate relative time value for months or hours */
830 	deltime.count = 1;
831 	/* Coverity[MIXED_ENUMS] */
832 	deltime.units = (CdTimeUnit)unit;
833 	switch(unit){
834 	  case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
835 		delta = etm - base_etm;
836 		if(!(timetype & cdStandardCal)){	/* Climatological time */
837 			hoursInYear = (timetype & cd365Days) ? 8760. : (timetype & cdHasLeap) ? 8784. : 8640.;
838 					     /* Normalize delta to interval [0,hoursInYear) */
839 			if(delta < 0.0 || delta >= hoursInYear) {
840 				double down = ((double)delta)/((double)hoursInYear);
841 				down = floor(down);
842 				down = down * (double)hoursInYear;
843 				delta = delta - down;
844 			}
845 		}
846 		break;
847 	  case cdYear: case cdSeason: case cdMonth:
848 		CdDivDelTime(base_etm, etm, deltime, old_timetype, 1970, &ndel);
849 		break;
850 	  case cdFraction:
851 	        cdError("invalid unit in conversion");
852 		break;
853 	  default: break;
854 	}
855 
856 					     /* Convert to output units */
857 	switch(unit){
858 	  case cdSecond:
859 		*reltime = 3600.0 * delta;
860 		break;
861 	  case cdMinute:
862 		*reltime = 60.0 * delta;
863 		break;
864 	  case cdHour:
865 		*reltime = delta;
866 		break;
867 	  case cdDay:
868 		*reltime = delta/24.0;
869 		break;
870 	  case cdWeek:
871 		*reltime = delta/168.0;
872 		break;
873 	  case cdMonth: case cdSeason: case cdYear: /* Already in correct units */
874 		if(timetype & cdStandardCal)
875 			*reltime = (base_etm <= etm) ? (double)ndel : (double)(-ndel);
876 		else			     /* Climatological time is already normalized*/
877 			*reltime = (double)ndel;
878 		break;
879 	  default:
880 		cdError("invalid unit in conversion");
881 		break;
882 	}
883 
884 	return;
885 }
886 
887 /* Add (value,unit) to comptime. */
888 /* value is in hours. */
889 /* calendar is anything but cdMixed. */
890 static void
cdCompAdd(cdCompTime comptime,double value,cdCalenType calendar,cdCompTime * result)891 cdCompAdd(cdCompTime comptime, double value, cdCalenType calendar, cdCompTime *result){
892 
893 	double reltime;
894 
895 	cdComp2Rel(calendar, comptime, "hours", &reltime);
896 	reltime += value;
897 	cdRel2Comp(calendar, "hours", reltime, result);
898 	return;
899 }
900 
901 /* Add value in hours to ct, in the mixed Julian/Gregorian
902  * calendar. */
903 static void
cdCompAddMixed(cdCompTime ct,double value,cdCompTime * result)904 cdCompAddMixed(cdCompTime ct, double value, cdCompTime *result){
905 
906 	static cdCompTime ZA = {1582, 10, 5, 0.0};
907 	static cdCompTime ZB = {1582, 10, 15, 0.0};
908 	double xj, xg;
909 
910 	if (cdCompCompare(ct, ZB) == -1){
911 		xj = cdDiffJulian(ZA, ct);
912 		if (value <= xj){
913 			cdCompAdd(ct, value, cdJulian, result);
914 		}
915 		else {
916 			cdCompAdd(ZB, value-xj, cdStandard, result);
917 		}
918 	}
919 	else {
920 		xg = cdDiffGregorian(ZB, ct);
921 		if (value > xg){
922 			cdCompAdd(ct, value, cdStandard, result);
923 		}
924 		else {
925 			cdCompAdd(ZA, value-xg, cdJulian, result);
926 		}
927 	}
928 	return;
929 }
930 
931 /* Return value expressed in hours. */
932 static double
cdToHours(double value,cdUnitTime unit)933 cdToHours(double value, cdUnitTime unit){
934 
935 	double result = 0;
936 
937 	switch(unit){
938 	case cdSecond:
939 		result = value/3600.0;
940 		break;
941 	case cdMinute:
942 		result = value/60.0;
943 		break;
944 	case cdHour:
945 		result = value;
946 		break;
947 	case cdDay:
948 		result = 24.0 * value;
949 		break;
950 	case cdWeek:
951 		result = 168.0 * value;
952 		break;
953 	default:
954 	        cdError("invalid unit in conversion");
955 		break;
956 
957 	}
958 	return result;
959 }
960 
961 /* Convert relative time (reltime, unit, basetime) to comptime in the
962  * mixed Julian/Gregorian calendar. unit is anything but year, season,
963  * month. unit and basetime are from the parsed relunits. Return
964  * result in comptime. */
965 static void
cdRel2CompMixed(double reltime,cdUnitTime unit,cdCompTime basetime,cdCompTime * comptime)966 cdRel2CompMixed(double reltime, cdUnitTime unit, cdCompTime basetime, cdCompTime *comptime){
967 
968 	reltime = cdToHours(reltime, unit);
969 	cdCompAddMixed(basetime, reltime, comptime);
970 	return;
971 }
972 
973 
974 static void
cdRel2Comp(cdCalenType timetype,char * relunits,double reltime,cdCompTime * comptime)975 cdRel2Comp(cdCalenType timetype, char* relunits, double reltime, cdCompTime* comptime)
976 {
977 	CdDeltaTime deltime;
978 	CdTime humantime;
979 	CdTimeType old_timetype;
980 	cdCompTime base_comptime;
981 	cdUnitTime unit, baseunits;
982 	double base_etm, result_etm;
983 	double delta;
984 	long idelta;
985 
986 					     /* Parse the relunits */
987 	if(cdParseRelunits(timetype, relunits, &unit, &base_comptime))
988 		return;
989 
990 	if (timetype == cdMixed){
991 		switch(unit){
992 		case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
993 			cdRel2CompMixed(reltime, unit, base_comptime, comptime);
994 			return;
995 		case cdYear: case cdSeason: case cdMonth:
996 			timetype = cdStandard;
997 			break;
998 		case cdFraction:
999 		        cdError("invalid unit in conversion");
1000 		        break;
1001 		default: break;
1002 		}
1003 	}
1004 
1005 	baseunits =cdBadUnit;
1006 	switch(unit){
1007 	  case cdSecond:
1008 		delta = reltime/3600.0;
1009 		baseunits = cdHour;
1010 		break;
1011 	  case cdMinute:
1012 		delta = reltime/60.0;
1013 		baseunits = cdHour;
1014 		break;
1015 	  case cdHour:
1016 		delta = reltime;
1017 		baseunits = cdHour;
1018 		break;
1019 	  case cdDay:
1020 		delta = 24.0 * reltime;
1021 		baseunits = cdHour;
1022 		break;
1023 	  case cdWeek:
1024 		delta = 168.0 * reltime;
1025 		baseunits = cdHour;
1026 		break;
1027 	  case cdMonth:
1028 		idelta = (long)(reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1029 		baseunits = cdMonth;
1030 		break;
1031 	  case cdSeason:
1032 		idelta = (long)(3.0 * reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1033 		baseunits = cdMonth;
1034 		break;
1035 	  case cdYear:
1036 		idelta = (long)(12 * reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1037 		baseunits = cdMonth;
1038 		break;
1039 	  default:
1040 	        cdError("invalid unit in conversion");
1041 		break;
1042 	}
1043 
1044 	deltime.count = 1;
1045 	/* Coverity[MIXED_ENUMS] */
1046 	deltime.units = (CdTimeUnit)baseunits;
1047 
1048 	humantime.year = base_comptime.year;
1049 	humantime.month = base_comptime.month;
1050 	humantime.day = base_comptime.day;
1051 	humantime.hour = base_comptime.hour;
1052 	humantime.baseYear = 1970;
1053 					     /* Map to old-style timetype */
1054 	if(cdToOldTimetype(timetype,&old_timetype))
1055 		return;
1056 	humantime.timeType = old_timetype;
1057 
1058 	Cdh2e(&humantime,&base_etm);
1059 					     /* If months, seasons, or years, */
1060 	if(baseunits == cdMonth){
1061 
1062 					     /* Calculate new epochal time from integer months. */
1063 					     /* Convert back to human, then comptime. */
1064 					     /* For zero reltime, just return the basetime*/
1065 		if(reltime != 0.0){
1066 			CdAddDelTime(base_etm,idelta,deltime,old_timetype,1970,&result_etm);
1067 			Cde2h(result_etm, old_timetype, 1970, &humantime);
1068 		}
1069 	}
1070 					     /* Calculate new epochal time. */
1071 					     /* Convert back to human, then comptime. */
1072 	else if(baseunits == cdHour){
1073 		Cde2h(base_etm+delta, old_timetype, 1970, &humantime);
1074 
1075 	}
1076 	comptime->year = humantime.year;
1077 	comptime->month = humantime.month;
1078 	comptime->day = humantime.day;
1079 	comptime->hour = humantime.hour;
1080 
1081 	return;
1082 }
1083 
1084 /* rkr: output as ISO 8601 strings */
1085 static void
cdComp2Iso(cdCalenType timetype,int separator,cdCompTime comptime,char * time)1086 cdComp2Iso(cdCalenType timetype, int separator, cdCompTime comptime, char* time)
1087 {
1088 	double dtmp, sec;
1089 	int ihr, imin, isec;
1090 	int nskip;
1091 
1092 	if(cdValidateTime(timetype,comptime))
1093 		return;
1094 
1095 	ihr = (int)comptime.hour;
1096 	dtmp = 60.0 * (comptime.hour - (double)ihr);
1097 	imin = (int)dtmp;
1098 	sec = 60.0 * (dtmp - (double)imin);
1099 	isec = (int)sec;
1100 
1101 	if(sec == isec)
1102 	    if(isec == 0)
1103 		if(imin == 0)
1104 		    if(ihr == 0)
1105 			nskip = 4;
1106 		    else
1107 			nskip = 3;
1108 		else
1109 		    nskip = 2;
1110 	    else
1111 		nskip = 1;
1112 	else
1113 	    nskip = 0;
1114 
1115 	if(timetype & cdStandardCal){
1116 	    switch (nskip) {
1117 	    case 0:		/* sec != 0 && (int)sec != sec */
1118 		sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d:%lf",
1119 			comptime.year,comptime.month,comptime.day,separator,ihr,imin,sec);
1120 		break;
1121 	    case 1:
1122 		sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d:%2.2d",
1123 			comptime.year,comptime.month,comptime.day,separator,ihr,imin,isec);
1124 		break;
1125 	    case 2:
1126 		sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d",
1127 			comptime.year,comptime.month,comptime.day,separator,ihr,imin);
1128 		break;
1129 	    case 3:
1130 		sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d",
1131 			comptime.year,comptime.month,comptime.day,separator,ihr);
1132 		break;
1133 	    case 4:
1134 		sprintf(time,"%4.4ld-%2.2hd-%2.2hd",
1135 			comptime.year,comptime.month,comptime.day);
1136 		break;
1137 	    }
1138 	}
1139 	else {				     /* Climatological */
1140 	    switch (nskip) {
1141 	    case 0:		/* sec != 0 && (int)sec != sec */
1142 		sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d:%lf",
1143 			comptime.month,comptime.day,separator,ihr,imin,sec);
1144 		break;
1145 	    case 1:
1146 		sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d:%2.2d",
1147 			comptime.month,comptime.day,separator,ihr,imin,isec);
1148 		break;
1149 	    case 2:
1150 		sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d",
1151 			comptime.month,comptime.day,separator,ihr,imin);
1152 		break;
1153 	    case 3:
1154 		sprintf(time,"%2.2hd-%2.2hd%c%2.2d",
1155 			comptime.month,comptime.day,separator,ihr);
1156 		break;
1157 	    case 4:
1158 		sprintf(time,"%2.2hd-%2.2hd",
1159 			comptime.month,comptime.day);
1160 		break;
1161 	    }
1162 	}
1163 	return;
1164 }
1165 
1166 /* rkr: added for output closer to ISO 8601 */
1167 void
cdRel2Iso(cdCalenType timetype,char * relunits,int separator,double reltime,char * chartime)1168 cdRel2Iso(cdCalenType timetype, char* relunits, int separator, double reltime, char* chartime)
1169 {
1170 	cdCompTime comptime;
1171 
1172 	cdRel2Comp(timetype, relunits, reltime, &comptime);
1173 	cdComp2Iso(timetype, separator, comptime, chartime);
1174 
1175 	return;
1176 }
1177 
1178 int
cdSetErrOpts(int opts)1179 cdSetErrOpts(int opts)
1180 {
1181     int old = cuErrOpts;
1182     cuErrOpts = opts;
1183     return old;
1184 }
1185