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

cgribexlib update.

parent 6562b5b2
......@@ -196,6 +196,7 @@ const char *cgribexLibraryVersion(void);
void gribDebug(int debug);
void gribSetCalendar(int calendar);
void gribDateTimeX(int *isec1, int *date, int *time, int *startDate, int *startTime);
void gribDateTime(int *isec1, int *date, int *time);
int gribRefDate(int *isec1);
int gribRefTime(int *isec1);
......
/* Automatically generated by m214003 at 2019-07-12, do not edit */
/* Automatically generated by m214003 at 2019-09-09, do not edit */
/* CGRIBEXLIB_VERSION="1.9.4" */
......@@ -47,6 +47,7 @@
#include "config.h"
#endif
#include <inttypes.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
......@@ -89,8 +90,17 @@
# define GRIBPACK unsigned char
#endif
#define U_BYTEORDER static union {unsigned long l; unsigned char c[sizeof(long)];} u_byteorder = {1}
#define IS_BIGENDIAN() (u_byteorder.c[sizeof(long) - 1])
#ifndef HOST_ENDIANNESS
#ifdef __cplusplus
static const uint32_t HOST_ENDIANNESS_temp[1] = { UINT32_C(0x00030201) };
#define HOST_ENDIANNESS (((const unsigned char *)HOST_ENDIANNESS_temp)[0])
#else
#define HOST_ENDIANNESS (((const unsigned char *)&(const uint32_t[1]){UINT32_C(0x00030201)})[0])
#endif
#endif
#define U_BYTEORDER
#define IS_BIGENDIAN() (HOST_ENDIANNESS == 0)
#if defined (__xlC__) /* performance problems on IBM */
#ifndef DBL_IS_NAN
......@@ -129,6 +139,25 @@ extern "C" {
#define intpow2(x) (ldexp(1.0, (x)))
static inline int
gribrec_len(unsigned b1, unsigned b2, unsigned b3)
{
/*
If bit 7 of b1 is set, we have to rescale by factor of 120.
This is a fixup to get round the restriction on product lengths
due to the count being only 24 bits. It is only possible because
the (default) rounding for GRIB products is 120 bytes.
*/
int needRescaling = b1 & (1 << 7);
int gribsize = (int)((((b1&127) << 16)+(b2<<8) + b3));
if ( needRescaling ) gribsize *= 120;
return gribsize;
}
unsigned correct_bdslen(unsigned bdslen, long recsize, long gribpos);
/* CDI converter routines */
......@@ -159,7 +188,8 @@ int cdiEncodeTime(int hour, int minute, int second);
extern FILE *grprsm;
extern int CGRIBEX_Debug;
extern int CGRIBEX_Debug, CGRIBEX_Fix_ZSE, CGRIBEX_Const;
extern int CGRIBEX_grib_calendar;
void gprintf(const char *caller, const char *fmt, ...);
......@@ -177,7 +207,6 @@ void scatter_complex_float(float *fpdata, int pcStart, int trunc, int nsp);
void gather_complex_double(double *fpdata, size_t pcStart, size_t trunc, size_t nsp);
void gather_complex_float(float *fpdata, size_t pcStart, size_t trunc, size_t nsp);
void scm0_double(double *pdl, double *pdr, double *pfl, double *pfr, int klg);
int qu2reg2(double *pfield, int *kpoint, int klat, int klon,
double *ztemp, double msval, int *kret);
int qu2reg3_double(double *pfield, int *kpoint, int klat, int klon,
......@@ -337,7 +366,6 @@ int grib2Sections(unsigned char *gribbuffer, long gribbufsize, unsigned char **i
#define BMS_Len ((bms) == NULL ? 0 : GET_UINT3(bms[0], bms[1], bms[2]))
#define BMS_UnusedBits (bms[3])
#define BMS_Numeric
#define BMS_Bitmap ((bms) == NULL ? NULL : (bms)+6)
#define BMS_BitmapSize (((((bms[0]<<16)+(bms[1]<<8)+bms[2]) - 6)<<3) - bms[3])
......@@ -368,7 +396,7 @@ int grib2Sections(unsigned char *gribbuffer, long gribbufsize, unsigned char **i
#define PutnZero(n) \
{ \
for ( size_t i = z >= 0 ? (size_t)z : 0; i < (size_t)(z+n); i++ ) lGrib[i] = 0; \
for ( size_t i___ = z >= 0 ? (size_t)z : 0; i___ < (size_t)(z+n); i___++ ) lGrib[i___] = 0; \
z += n; \
}
......@@ -1186,8 +1214,18 @@ xlc_r -g -O3 -qhot -q64 -qarch=auto -qtune=auto -qreport -DTEST_ENCODE encode_ar
#include <stdio.h>
#include <stdlib.h>
#define GRIBPACK unsigned char
#define IS_BIGENDIAN() (u_byteorder.c[sizeof(long) - 1])
#define U_BYTEORDER static union {unsigned long l; unsigned char c[sizeof(long)];} u_byteorder = {1}
#ifndef HOST_ENDIANNESS
#ifdef __cplusplus
static const uint32_t HOST_ENDIANNESS_temp[1] = { UINT32_C(0x00030201) };
#define HOST_ENDIANNESS (((const unsigned char *)HOST_ENDIANNESS_temp)[0])
#else
#define HOST_ENDIANNESS (((const unsigned char *)&(const uint32_t[1]){UINT32_C(0x00030201)})[0])
#endif
#endif
#define U_BYTEORDER
#define IS_BIGENDIAN() (HOST_ENDIANNESS == 0)
#define Error(x,y)
#endif
......@@ -1515,7 +1553,6 @@ double dtime()
return (tseconds);
}
#define NRUN 10000
static
void pout(char *name, int s, unsigned char *lgrib, long datasize, double tt)
......@@ -1526,7 +1563,11 @@ void pout(char *name, int s, unsigned char *lgrib, long datasize, double tt)
int main(void)
{
long datasize = 1000000;
enum {
datasize = 1000000,
NRUN = 10000,
};
double t_begin, t_end;
float *dataf = (float*) malloc(datasize*sizeof(float));
......@@ -1937,49 +1978,67 @@ double decfp2(int kexp, int kmant)
#include <stdarg.h>
int gribRefDate(int *isec1)
static
void gribDecodeRefDate(int *isec1, int *year, int *month, int *day)
{
int century = ISEC1_Century;
int ryear = ISEC1_Year;
if ( century < 0 ) century = -century;
century -= 1;
int ryear = ISEC1_Year;
/* if ( century != 0 ) */
if ( century == -255 && ryear == 127 )
{
if ( ryear == 100 )
{
ryear = 0;
century += 1;
}
century = 0;
ryear = 0;
}
else
{
/* if ( century != 0 ) */
{
if ( ryear == 100 )
{
ryear = 0;
century += 1;
}
if ( ryear != 255 )
{
ryear = century*100 + ryear;
if ( ISEC1_Century < 0 ) ryear = -ryear;
}
else
ryear = 1;
if ( ryear != 255 )
{
ryear = century*100 + ryear;
if ( ISEC1_Century < 0 ) ryear = -ryear;
}
else
ryear = 1;
}
}
int rmonth = ISEC1_Month;
int rday = ISEC1_Day;
*year = ryear;
*month = ISEC1_Month;
*day = ISEC1_Day;
}
int date = cdiEncodeDate(ryear, rmonth, rday);
return date ;
int gribRefDate(int *isec1)
{
int ryear, rmonth, rday;
gribDecodeRefDate(isec1, &ryear, &rmonth, &rday);
return cdiEncodeDate(ryear, rmonth, rday);
}
int gribRefTime(int *isec1)
static
void gribDecodeRefTime(int *isec1, int *hour, int *minute, int *second)
{
int rhour = ISEC1_Hour;
int rminute = ISEC1_Minute;
*hour = ISEC1_Hour;
*minute = ISEC1_Minute;
*second = 0;
}
int time = cdiEncodeTime(rhour, rminute, 0);
return time;
int gribRefTime(int *isec1)
{
int rhour, rminute, rsecond;
gribDecodeRefTime(isec1, &rhour, &rminute, &rsecond);
return cdiEncodeTime(rhour, rminute, rsecond);
}
......@@ -1987,7 +2046,7 @@ bool gribTimeIsFC(int *isec1)
{
bool isFC = false;
int time_period = (ISEC1_TimeRange == 10) ? (ISEC1_TimePeriod1<<8) + ISEC1_TimePeriod2 : ISEC1_TimePeriod1;
int time_period = (ISEC1_TimeRange == 10) ? (ISEC1_TimePeriod1<<8) + ISEC1_TimePeriod2 : ISEC1_TimePeriod1;
if ( time_period > 0 && ISEC1_Day > 0 )
{
......@@ -1997,94 +2056,80 @@ bool gribTimeIsFC(int *isec1)
return isFC;
}
void gribDateTime(int *isec1, int *date, int *time)
static
int getTimeUnitFactor(int timeUnit)
{
static bool lprint = true;
int64_t julday;
int secofday;
int64_t addsec = 0;
int64_t time_period = 0;
extern int grib_calendar;
int century = ISEC1_Century;
int ryear = ISEC1_Year;
if ( century == -255 && ryear == 127 )
{
century = 0;
ryear = 0;
switch ( timeUnit )
{
case ISEC1_TABLE4_MINUTE: return 60; break;
case ISEC1_TABLE4_QUARTER: return 900; break;
case ISEC1_TABLE4_30MINUTES: return 1800; break;
case ISEC1_TABLE4_HOUR: return 3600; break;
case ISEC1_TABLE4_3HOURS: return 10800; break;
case ISEC1_TABLE4_6HOURS: return 21600; break;
case ISEC1_TABLE4_12HOURS: return 43200; break;
case ISEC1_TABLE4_DAY: return 86400; break;
default:
if ( lprint )
{
gprintf(__func__, "Time unit %d unsupported", timeUnit);
lprint = false;
}
break;
}
else
{
if ( century < 0 ) century = -century;
century -= 1;
/* if ( century != 0 ) */
{
if ( ryear == 100 )
{
ryear = 0;
century += 1;
}
return 0;
}
if ( ryear != 255 )
{
ryear = century*100 + ryear;
if ( ISEC1_Century < 0 ) ryear = -ryear;
}
else
ryear = 1;
}
}
int rmonth = ISEC1_Month;
int rday = ISEC1_Day;
void gribDateTimeX(int *isec1, int *date, int *time, int *startDate, int *startTime)
{
*startDate = 0;
*startTime = 0;
int ryear, rmonth, rday;
gribDecodeRefDate(isec1, &ryear, &rmonth, &rday);
int rhour = ISEC1_Hour;
int rminute = ISEC1_Minute;
int second = 0;
int rhour, rminute, rsecond;
gribDecodeRefTime(isec1, &rhour, &rminute, &rsecond);
/* printf("ref %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute); */
// printf("ref %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute);
int64_t time_period = 0, time_period_x = 0;
if ( ISEC1_TimeRange == 10 )
time_period = (ISEC1_TimePeriod1<<8) + ISEC1_TimePeriod2;
else if ( ISEC1_TimeRange >=2 && ISEC1_TimeRange <= 5 )
time_period = ISEC1_TimePeriod2;
{
time_period_x = ISEC1_TimePeriod1;
time_period = ISEC1_TimePeriod2;
}
else if ( ISEC1_TimeRange == 0 )
time_period = ISEC1_TimePeriod1;
if ( time_period > 0 && rday > 0 )
{
encode_caldaysec(grib_calendar, ryear, rmonth, rday, rhour, rminute, second, &julday, &secofday);
addsec = 0;
switch ( ISEC1_TimeUnit )
{
case ISEC1_TABLE4_MINUTE: addsec = 60 * time_period; break;
case ISEC1_TABLE4_QUARTER: addsec = 900 * time_period; break;
case ISEC1_TABLE4_30MINUTES: addsec = 1800 * time_period; break;
case ISEC1_TABLE4_HOUR: addsec = 3600 * time_period; break;
case ISEC1_TABLE4_3HOURS: addsec = 10800 * time_period; break;
case ISEC1_TABLE4_6HOURS: addsec = 21600 * time_period; break;
case ISEC1_TABLE4_12HOURS: addsec = 43200 * time_period; break;
case ISEC1_TABLE4_DAY: addsec = 86400 * time_period; break;
default:
if ( lprint )
{
gprintf(__func__, "Time unit %d unsupported", ISEC1_TimeUnit);
lprint = false;
}
break;
}
int64_t julday;
int secofday;
encode_caldaysec(CGRIBEX_grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);
julday_add_seconds(addsec, &julday, &secofday);
int time_unit_factor = getTimeUnitFactor(ISEC1_TimeUnit);
decode_caldaysec(grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &second);
if (time_period_x > 0)
{
int64_t julday_x = julday;
int secofday_x = secofday;
julday_add_seconds(time_unit_factor * time_period_x, &julday_x, &secofday_x);
decode_caldaysec(CGRIBEX_grib_calendar, julday_x, secofday_x, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
*startDate = cdiEncodeDate(ryear, rmonth, rday);
*startTime = cdiEncodeTime(rhour, rminute, 0);
}
julday_add_seconds(time_unit_factor * time_period, &julday, &secofday);
decode_caldaysec(CGRIBEX_grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
}
/*
printf("new %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute);
*/
// printf("new %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute);
*date = cdiEncodeDate(ryear, rmonth, rday);
*time = cdiEncodeTime(rhour, rminute, 0);
......@@ -2092,6 +2137,13 @@ void gribDateTime(int *isec1, int *date, int *time)
}
void gribDateTime(int *isec1, int *date, int *time)
{
int sdate, stime;
gribDateTimeX(isec1, date, time, &sdate, &stime);
}
void gprintf(const char *caller, const char *fmt, ...)
{
va_list args;
......@@ -2102,7 +2154,7 @@ void gprintf(const char *caller, const char *fmt, ...)
fprintf(grprsm, "%-18s : ", caller);
vfprintf(grprsm, fmt, args);
fprintf(grprsm, "\n");
fputs("\n", grprsm);
va_end(args);
}
......@@ -3003,7 +3055,7 @@ void gribPrintSec1(int *isec0, int *isec1)
}
}
void printQuasi(int *isec2)
static void printQuasi(int *isec2)
{
/*
......@@ -3861,7 +3913,7 @@ size_t gribReadSize(int fileID)
unsigned pdssize = read3ByteMSBFirst(fileptr);
if ( CGRIBEX_Debug ) Message("pdssize = %u", pdssize);
int flag;
int flag = 0;
for ( int i = 0; i < 5; ++i ) flag = filePtrGetc(fileptr);
if ( CGRIBEX_Debug ) Message("flag = %d", flag);
......@@ -4002,25 +4054,12 @@ int gribWrite(int fileID, unsigned char *buffer, size_t buffersize)
FILE *grprsm = NULL;
double fref;
double fmaxval;
int nfref;
int nfmaxval;
int nrnd;
int ndbg;
int nvck;
int nonoff;
int noabort;
int num2ok;
int next2o;
int nloc2o;
int nsubce;
int grib_calendar = -1;
int CGRIBEX_grib_calendar = -1;
void gribSetCalendar(int calendar)
{
grib_calendar = calendar;
CGRIBEX_grib_calendar = calendar;
}
......@@ -4107,67 +4146,31 @@ C ----------------------------------------------------------------
*/
/*
Common area variables have not been set. Set them.
User supplied reference value.
*/
fref = 0.0;
/*
Reference value supplied by user flag. Set to off.
*/
nfref = 0;
/*
User supplied maximum value.
*/
fmaxval = 0.0;
/*
Maximum value supplied by user flag. Set to off.
*/
nfmaxval = 0;
/*
Set rounding to 120 bytes on.
*/
nrnd = 1;
/*
Set GRIB calendar.
*/
if ( grib_calendar == -1 )
if ( CGRIBEX_grib_calendar == -1 )
{
grib_calendar = CALENDAR_PROLEPTIC;
CGRIBEX_grib_calendar = CALENDAR_PROLEPTIC;
envString = getenv("GRIB_CALENDAR");
if ( envString )
{
if ( strncmp(envString, "standard", 8) == 0 )
grib_calendar = CALENDAR_STANDARD;
CGRIBEX_grib_calendar = CALENDAR_STANDARD;
else if ( strncmp(envString, "proleptic", 9) == 0 )
grib_calendar = CALENDAR_PROLEPTIC;
CGRIBEX_grib_calendar = CALENDAR_PROLEPTIC;
else if ( strncmp(envString, "360days", 7) == 0 )
grib_calendar = CALENDAR_360DAYS;
CGRIBEX_grib_calendar = CALENDAR_360DAYS;
else if ( strncmp(envString, "365days", 7) == 0 )
grib_calendar = CALENDAR_365DAYS;
CGRIBEX_grib_calendar = CALENDAR_365DAYS;
else if ( strncmp(envString, "366days", 7) == 0 )
grib_calendar = CALENDAR_366DAYS;
CGRIBEX_grib_calendar = CALENDAR_366DAYS;
else if ( strncmp(envString, "none", 4) == 0 )
grib_calendar = CALENDAR_NONE;
CGRIBEX_grib_calendar = CALENDAR_NONE;
}
}
/*
Set debug print off.
*/
ndbg = 0;
envString = getenv("GRIBEX_DEBUG");
if ( envString != NULL )
{
if ( !strncmp(envString, "ON", 2) )
ndbg = 1;
else if( *envString == '1')
ndbg = 1;
else if( *envString == '2')
ndbg = 2;
else
ndbg = 0;
}
/*
Set GRIBEX compatibility mode.
*/
......@@ -4177,19 +4180,6 @@ C ----------------------------------------------------------------
if ( atoi(envString) == 1 ) CGRIBEX_Const = 0;
}
/*
Set GRIB value checking on.
*/
nvck = 1;
envString = getenv("GRIBEX_CHECK");
if ( envString )
{
if ( !strncmp(envString, "OFF", 3) )
nvck = 0;
else
nvck = 1;
}
/*
See if output stream needs changing
*/
......@@ -4226,38 +4216,10 @@ C ----------------------------------------------------------------
}
}
}
/*
Set P factor switch to default, user supplies the P factor.
*/
nonoff = 0;
/*
Set abort flag to NO abort
*/
noabort = 1;
/*
Mark common area values set by user.
*/
lfirst = false;
/*
Exhaustive use of all possible second-order packing methods
for HOPER='K'. Set to off.
*/
num2ok = 0;
/*
Use of extended second-order packing methods for grid-point
encoding (HOPER='C' and 'K'). Set to on.
*/
next2o = 1;
/*
Use of non-local second-order packing methods for grid-point
encoding (HOPER='C' and 'K'). Set to on.
*/
nloc2o = 1;
/*
Use of (all valid) sub-centre values for ECMWF fields encoding .
encoding. Set to off.
*/
nsubce = 0;
}
/* pack 8-bit bytes from 64-bit words to a packed buffer */
......@@ -4299,7 +4261,7 @@ long packInt64(unsigned INT64 *up, unsigned char *cp, long bc, long tc)
ip6 = ip0 + 6;
ip7 = ip0 + 7;
up0 = (unsigned INT64 *) (cp + head);
up0 = (unsigned INT64 *)(void *)(cp + head);
/* Here we should process any bytes until the first word boundary
* of our destination buffer
......@@ -4376,7 +4338,6 @@ long unpackInt64(const unsigned char *cp, unsigned INT64 *up, long bc, long tc)
{
U_BYTEORDER;
const unsigned char *cp0;
unsigned INT64 *up0;
unsigned INT64 *ip0, *ip1, *ip2, *ip3, *ip4, *ip5, *ip6, *ip7;
long head, trail, inner, i, j;
long offset;
......@@ -4410,7 +4371,7 @@ long unpackInt64(const unsigned char *cp, unsigned INT64 *up, long bc, long tc)
ip6 = ip0 + 6;
ip7 = ip0 + 7;
up0 = (unsigned INT64 *) (cp + head);
const unsigned INT64 *up0 = (const unsigned INT64 *)(const void *)(cp + head);
/* Process any bytes until the first word boundary
* of our source buffer
......@@ -4509,7 +4470,7 @@ long packInt32(unsigned INT32 *up, unsigned char *cp, long bc, long tc)
ip2 = ip0 + 2;
ip3 = ip0 + 3;
up0 = (unsigned INT32 *) (cp + head);
up0 = (unsigned INT32 *)(void *)(cp + head);
/* Here we should process any bytes until the first word boundary
* of our destination buffer
......@@ -4580,7 +4541,6 @@ long unpackInt32(const unsigned char *cp, unsigned INT32 *up, long bc, long tc)
{
U_BYTEORDER;
const unsigned char *cp0;
unsigned INT32 *up0;