Commit fb7bb4c1 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

replace calendar module by a new one

parent e56472f8
2007-09-?? Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* using GRIB library version 1.0.6
* replace calendar module by a new one because of date/time problems
with years < 0 [report: Veronika Gayler]
* add module timebase
* use always decode_date and encode_date to decode/encode date
* Version 1.0.8 released
......
......@@ -585,7 +585,7 @@ static void printShortinfo(int streamID, int vlistID, int vardis)
decode_date(vdate, &year, &month, &day);
decode_time(vtime, &hour, &minute);
fprintf(stdout, " %4.4d-%2.2d-%2.2d %2.2d:%2.2d", year, month, day, hour, minute);
fprintf(stdout, " %6.4d-%2.2d-%2.2d %2.2d:%2.2d", year, month, day, hour, minute);
ntimeout++;
tsID++;
}
......
......@@ -17,6 +17,7 @@ libcdi_a_SOURCES = \
util.c \
varscan.c \
calendar.c \
timebase.c \
vlist.c \
vlist_var.c \
basetime.c \
......@@ -58,6 +59,7 @@ libcdi_a_SOURCES = \
cdi.h \
cdi_limits.h \
calendar.h \
timebase.h \
vlist.h \
taxis.h \
stream_grb.h \
......
......@@ -99,6 +99,7 @@ libcdi_a_SOURCES = \
util.c \
varscan.c \
calendar.c \
timebase.c \
vlist.c \
vlist_var.c \
basetime.c \
......@@ -140,6 +141,7 @@ libcdi_a_SOURCES = \
cdi.h \
cdi_limits.h \
calendar.h \
timebase.h \
vlist.h \
taxis.h \
stream_grb.h \
......@@ -172,8 +174,8 @@ am_libcdi_a_OBJECTS = cdiFortran.$(OBJEXT) cdi_error.$(OBJEXT) \
taxis.$(OBJEXT) error.$(OBJEXT) version.$(OBJEXT) \
dmemory.$(OBJEXT) model.$(OBJEXT) institution.$(OBJEXT) \
table.$(OBJEXT) util.$(OBJEXT) varscan.$(OBJEXT) \
calendar.$(OBJEXT) vlist.$(OBJEXT) vlist_var.$(OBJEXT) \
basetime.$(OBJEXT) stream_history.$(OBJEXT) \
calendar.$(OBJEXT) timebase.$(OBJEXT) vlist.$(OBJEXT) \
vlist_var.$(OBJEXT) basetime.$(OBJEXT) stream_history.$(OBJEXT) \
stream_grb.$(OBJEXT) stream_cdf.$(OBJEXT) stream_srv.$(OBJEXT) \
stream_ext.$(OBJEXT) stream_ieg.$(OBJEXT) grid.$(OBJEXT) \
grid_rot.$(OBJEXT) grid_gme.$(OBJEXT) zaxis.$(OBJEXT) \
......@@ -207,10 +209,11 @@ am__depfiles_maybe = depfiles
@AMDEP_TRUE@ ./$(DEPDIR)/stream_record.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_srv.Po ./$(DEPDIR)/stream_var.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/swap.Po ./$(DEPDIR)/table.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/taxis.Po ./$(DEPDIR)/tsteps.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/util.Po ./$(DEPDIR)/varscan.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/version.Po ./$(DEPDIR)/vlist.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/vlist_var.Po ./$(DEPDIR)/zaxis.Po
@AMDEP_TRUE@ ./$(DEPDIR)/taxis.Po ./$(DEPDIR)/timebase.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/tsteps.Po ./$(DEPDIR)/util.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/varscan.Po ./$(DEPDIR)/version.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/vlist.Po ./$(DEPDIR)/vlist_var.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/zaxis.Po
COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
CCLD = $(CC)
......@@ -322,6 +325,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/swap.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/table.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/taxis.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/timebase.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/tsteps.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/util.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/varscan.Po@am__quote@
......
#ifndef _BASETIME_H
#define _BASETIME_H
#ifndef _CALENDAR_H
#include "calendar.h"
#endif
typedef struct {
int ncvarid;
......
#if defined (HAVE_CONFIG_H)
# include "config.h"
#endif
#include <stdio.h> /* for NULL */
#include <math.h> /* for modf(), floor(), log10(), ceil() */
#include <float.h> /* for DBL_MANT_DIG */
#include <stdio.h>
#include "calendar.h"
#include "cdi.h" /* CALENDAR_ */
#include "timebase.h"
#undef ABS
#define ABS(a) ((a) < 0 ? -(a) : (a))
static int month_360[12] = {30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30};
static int month_365[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
static int month_366[12] = {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
/*
* The following two functions convert between Julian day number and
* Gregorian/Julian dates (Julian dates are used prior to October 15,
* 1582; Gregorian dates are used after that). Julian day number 0 is
* midday, January 1, 4713 BCE. The Gregorian calendar was adopted
* midday, October 15, 1582.
*
* Author: Robert Iles, March 1994
*
* C Porter: Steve Emmerson, October 1995
*
* Original: http://www.nag.co.uk:70/nagware/Examples/calendar.f90
*
* There is no warranty on this code.
*/
/*
* Convert a Julian day number to a Gregorian/Julian date.
*/
static void julday_to_gregdate(
unsigned long julday, /* Julian day number to convert */
int *year, /* Gregorian year (out) */
int *month, /* Gregorian month (1-12) (out) */
int *day /* Gregorian day (1-31) (out) */
)
{
long ja, jb;
int jc;
long jd;
int je, iday, imonth, iyear;
double xc;
if (julday < 2299161)
ja = julday;
else
{
int ia = (int) (((julday - 1867216) - 0.25) / 36524.25);
ja = julday + 1 + ia - (int)(0.25 * ia);
}
jb = ja + 1524;
xc = ((jb - 2439870) - 122.1) / 365.25;
jc = (long) (6680.0 + xc);
jd = 365 * jc + (int)(0.25 * jc);
je = (int) ((jb - jd) / 30.6001);
iday = (int) (jb - jd - (int)(30.6001 * je));
imonth = je - 1;
if (imonth > 12)
imonth -= 12;
int calendar_dpy(int calendar)
{
int dpy = 0;
iyear = jc - 4715;
if (imonth > 2)
iyear -= 1;
if (iyear <= 0)
iyear -= 1;
if ( calendar == CALENDAR_360DAYS ) dpy = 360;
else if ( calendar == CALENDAR_365DAYS ) dpy = 365;
else if ( calendar == CALENDAR_366DAYS ) dpy = 366;
*year = iyear;
*month = imonth;
*day = iday;
return (dpy);
}
/*
* Convert a Gregorian/Julian date to a Julian day number.
*
* The Gregorian calendar was adopted midday, October 15, 1582.
*/
static unsigned long gregdate_to_julday(
int year, /* Gregorian year */
int month, /* Gregorian month (1-12) */
int day /* Gregorian day (1-31) */
)
int days_per_month(int calendar, int year, int month)
{
long igreg = 15 + 31 * (10 + (12 * 1582));
int iy; /* signed, origin 0 year */
int ja; /* Julian century */
int jm; /* Julian month */
int jy; /* Julian year */
unsigned long julday; /* returned Julian day number */
/*
* Because there is no 0 BC or 0 AD, assume the user wants the start of
* the common era if they specify year 0.
*/
if (year == 0)
year = 1;
iy = year;
if (year < 0)
iy++;
if (month > 2)
{
jy = iy;
jm = month + 1;
}
else
{
jy = iy - 1;
jm = month + 13;
}
int dayspermonth = 0;
int *dpm = NULL;
int dpy;
/*
* Note: SLIGHTLY STRANGE CONSTRUCTIONS REQUIRED TO AVOID PROBLEMS WITH
* OPTIMISATION OR GENERAL ERRORS UNDER VMS!
*/
julday = day + (int)(30.6001 * jm);
if (jy >= 0)
{
julday += 365 * jy;
julday += jy / 4;
}
else
{
double xi = 365.25 * jy;
dpy = calendar_dpy(calendar);
if ((int)xi != xi)
xi -= 1;
julday += (int)xi;
}
julday += 1720995;
if ( dpy == 360 ) dpm = month_360;
else if ( dpy == 365 ) dpm = month_365;
else dpm = month_366;
if (day + (31* (month + (12 * iy))) >= igreg)
if ( month >= 1 && month <= 12 )
dayspermonth = dpm[month-1];
else
fprintf(stderr, "days_per_month: month %d out of range\n", month);
if ( dpy == 0 && month == 2 )
{
ja = jy/100;
julday -= ja;
julday += 2;
julday += ja/4;
if ( (year%4 == 0 && year%100 != 0) || year%400 == 0 )
dayspermonth = 29;
else
dayspermonth = 28;
}
return julday;
}
/*
* Encode a date as a double-precision value.
*/
static double utencdate(int year, int month, int day)
{
return ((long)gregdate_to_julday(year, month, day) -
(long)gregdate_to_julday(2001, 1, 1)) * 86400.0;
return (dayspermonth);
}
/*
* Encode a time as a double-precision value.
*/
static double utencclock(int hours, int minutes, double seconds)
int days_per_year(int calendar, int year)
{
return (hours*60 + minutes)*60 + seconds;
}
int daysperyear;
int dpy;
dpy = calendar_dpy(calendar);
/*
* Decompose a value into a set of values accounting for uncertainty.
*/
static void decomp(
double value,
double uncer, /* >= 0 */
int nbasis,
double *basis, /* all values > 0 */
double *count
)
{
int i;
for (i = 0; i < nbasis; i++)
if ( dpy == 0 )
{
double r = fmod(value, basis[i]); /* remainder */
/* Adjust remainder to minimum magnitude. */
if (ABS(2*r) > basis[i])
r += r > 0
? -basis[i]
: basis[i];
if (ABS(r) <= uncer)
{
/* The value equals a basis multiple within the uncertainty. */
double half = value < 0 ? -basis[i]/2 : basis[i]/2;
modf((value+half)/basis[i], count+i);
break;
}
value = basis[i] * modf(value/basis[i], count+i);
if ( year == 1582 )
dpy = 355;
else if ( (year%4 == 0 && year%100 != 0) || year%400 == 0 )
dpy = 366;
else
dpy = 365;
}
for (i++; i < nbasis; i++)
count[i] = 0;
daysperyear = dpy;
return (daysperyear);
}
/*
* Decode a time from a double-precision value.
*/
static void dectime(
double value,
int *year,
int *month,
int *day,
int *hour,
int *minute,
float *second)
static void decode_day(int dpy, int days, int *year, int *month, int *day)
{
long days;
int hours;
int minutes;
double seconds;
double uncer; /* uncertainty of input value */
typedef union
{
double vec[7];
struct
{
double days;
double hours12;
double hours;
double minutes10;
double minutes;
double seconds10;
double seconds;
} ind;
} Basis;
Basis counts;
static Basis basis;
basis.ind.days = 86400;
basis.ind.hours12 = 43200;
basis.ind.hours = 3600;
basis.ind.minutes10 = 600;
basis.ind.minutes = 60;
basis.ind.seconds10 = 10;
basis.ind.seconds = 1;
uncer = ldexp(value < 0 ? -value : value, -DBL_MANT_DIG);
days = (long) floor(value/86400.0);
value -= days * 86400.0; /* make positive excess */
decomp(value, uncer, sizeof(basis.vec)/sizeof(basis.vec[0]),
basis.vec, counts.vec);
days += (long) counts.ind.days;
hours = (int)counts.ind.hours12 * 12 + (int)counts.ind.hours;
minutes = (int)counts.ind.minutes10 * 10 + (int)counts.ind.minutes;
seconds = (int)counts.ind.seconds10 * 10 + counts.ind.seconds;
*second = seconds;
*minute = minutes;
*hour = hours;
julday_to_gregdate(gregdate_to_julday(2001, 1, 1) + days, year, month, day);
}
int i = 0;
int *dpm = NULL;
*year = (days-1) / dpy;
days -= (*year*dpy);
/*
* Convert a Gregorian/Julian date and time into a temporal value.
*
*/
void InvCalendar(int year, int month, int day, int hour,
int minute, double second, sUnit unit, double *value)
{
*value = (utencdate(year, month, day) +
utencclock(hour, minute, second) - unit.origin) /
unit.factor;
}
if ( dpy == 360 ) dpm = month_360;
else if ( dpy == 365 ) dpm = month_365;
else if ( dpy == 366 ) dpm = month_366;
/*
* Convert a temporal value into UTC Gregorian/Julian date and time.
*
*/
void Calendar(double value, sUnit unit, int *year, int *month,
int *day, int *hour, int *minute, double *second)
{
float sec;
if ( dpm )
for ( i = 0; i < 12; i++ )
{
if ( days > dpm[i] ) days -= dpm[i];
else break;
}
dectime(unit.origin + value*unit.factor, year, month, day, hour, minute, &sec);
*second = sec;
*month = i + 1;
*day = days;
}
/*
* Encode a date as a double-precision value.
*/
static double encdateXXX(int dpy, int year, int month, int day)
static int encode_day(int dpy, int year, int month, int day)
{
int i;
int *dpm = NULL;
......@@ -315,170 +111,204 @@ static double encdateXXX(int dpy, int year, int month, int day)
if ( dpm ) for ( i = 0; i < month-1; i++ ) rval += dpm[i];
return (rval * 86400.0);
return (rval);
}
/*
* Convert a XXX days/year date and time into a temporal value.
*
*/
static void InvCalendarXXX(int dpy, int year, int month, int day, int hour,
int minute, double second, sUnit unit, double *value)
int date_to_calday(int calendar, int date)
{
*value = (encdateXXX(dpy, year, month, day) +
utencclock(hour, minute, second) - unit.origin) /
unit.factor;
}
int calday;
int dpy;
int year, month, day;
dpy = calendar_dpy(calendar);
/*
* Convert a temporal value into XXX days/year date and time.
*
*/
static void dectimeXXX(int dpy, double value, int *year, int *month,
int *day, int *hour, int *minute, float *second)
{
int i = 0;
int *dpm = NULL;
long days;
int hours;
int minutes;
double seconds;
double uncer; /* uncertainty of input value */
typedef union
{
double vec[7];
struct
{
double days;
double hours12;
double hours;
double minutes10;
double minutes;
double seconds10;
double seconds;
} ind;
} Basis;
Basis counts;
static Basis basis;
basis.ind.days = 86400;
basis.ind.hours12 = 43200;
basis.ind.hours = 3600;
basis.ind.minutes10 = 600;
basis.ind.minutes = 60;
basis.ind.seconds10 = 10;
basis.ind.seconds = 1;
uncer = ldexp(value < 0 ? -value : value, -DBL_MANT_DIG);
days = (long) floor(value/86400.0);
value -= days * 86400.0; /* make positive excess */
decomp(value, uncer, sizeof(basis.vec)/sizeof(basis.vec[0]),
basis.vec, counts.vec);
days += (long) counts.ind.days;
hours = (int)counts.ind.hours12 * 12 + (int)counts.ind.hours;
minutes = (int)counts.ind.minutes10 * 10 + (int)counts.ind.minutes;
seconds = (int)counts.ind.seconds10 * 10 + counts.ind.seconds;
*second = seconds;
*minute = minutes;
*hour = hours;
*year = (days-1) / dpy;
days -= (*year*dpy);
if ( dpy == 360 ) dpm = month_360;
else if ( dpy == 365 ) dpm = month_365;
else if ( dpy == 366 ) dpm = month_366;
if ( dpm )
for ( i = 0; i < 12; i++ )
{
if ( days > dpm[i] ) days -= dpm[i];
else break;
}
decode_date(date, &year, &month, &day);
*month = i + 1;
if ( dpy == 360 || dpy == 365 || dpy == 366 )
calday = encode_day(dpy, year, month, day);
else
calday = encode_julday(year, month, day);
*day = days;
return (calday);
}
/*
* Convert a temporal value into XXX days/year date and time.
*
*/
static void CalendarXXX(int dpy, double value, sUnit unit, int *year, int *month,
int *day, int *hour, int *minute, double *second)
int calday_to_date(int calendar, int calday)
{
float sec;
int date;
int dpy;
int year, month, day;
dpy = calendar_dpy(calendar);
if ( dpy == 360 || dpy == 365 || dpy == 366 )
decode_day(dpy, calday, &year, &month, &day);
else
decode_julday(calday, &year, &month, &day);
date = encode_date(year, month, day);
dectimeXXX(dpy, unit.origin + value*unit.factor, year, month, day, hour, minute, &sec);
*second = sec;
return (date);
}
/*
* Convert a cdi date and time into a temporal value.
*
*/
void cdiInvCalendar(int dpy, int year, int month, int day, int hour,
int minute, double second, sUnit unit, double *value)
void encode_caldaysec(int calendar, int year, int month, int day, int hour, int minute, int *julday, int *secofday)
{