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

cgribexlib update

parent 335a6935
......@@ -166,8 +166,9 @@
void gribFixZSE(int flag);
void gribSetDebug(int debug);
void gribFixZSE(int flag); /* 1: Fix ZeroShiftError of simple packed spherical harmonics */
void gribSetConst(int flag); /* 1: Don't pack constant fields on regular grids */
void gribSetDebug(int debug); /* 1: Debugging */
void gribSetRound(int round);
void gribSetRefDP(double refval);
void gribSetRefSP(float refval);
......
/* Automatically generated by m214003 at 2010-08-23, do not edit */
/* Automatically generated by m214003 at 2010-09-09, do not edit */
/* CGRIBEXLIB_VERSION="1.4.6" */
......@@ -1276,8 +1276,9 @@ gribExSP(int *isec0, int *isec1, int *isec2, float *fsec2sp, int *isec3,
}
}
int GRB_Debug = 0; /* If set to 1, debugging */
int GRB_Fix_ZSE = 0; /* Fix ZeroShiftError of simple packed spherical harmonics */
int GRB_Fix_ZSE = 0; /* 1: Fix ZeroShiftError of simple packed spherical harmonics */
int GRB_Const = 0; /* 1: Don't pack constant fields on regular grids */
int GRB_Debug = 0; /* 1: Debugging */
void gribSetDebug(int debug)
{
......@@ -1289,6 +1290,7 @@ void gribSetDebug(int debug)
Message(func, "debug level %d", debug);
}
void gribFixZSE(int flag)
{
static const char *func = "gribFixZSE";
......@@ -1299,19 +1301,34 @@ void gribFixZSE(int flag)
Message(func, "Fix ZeroShiftError set to %d", flag);
}
void gribSetConst(int flag)
{
static const char *func = "gribSet Const";
GRB_Const = flag;
if ( GRB_Debug )
Message(func, "Const set to %d", flag);
}
void gribSetRound(int round)
{
}
void gribSetRefDP(double refval)
{
}
void gribSetRefSP(float refval)
{
gribSetRefDP((double) refval);
}
void gribSetValueCheck(int vcheck)
{
}
......@@ -2819,7 +2836,6 @@ void encodeIS(GRIBPACK *lGrib, long *gribLen)
*gribLen = z;
}
/* GRIB block 5 - end block */
static
void encodeES(GRIBPACK *lGrib, long *gribLen, long bdsstart)
......@@ -3118,7 +3134,6 @@ void encodePDS(GRIBPACK *lpds, long pdsLen, int *isec1)
}
}
/* GRIB BLOCK 2 - GRID DESCRIPTION SECTION */
static
void encodeGDS(GRIBPACK *lGrib, long *gribLen, int *isec2, double *fsec2)
......@@ -3251,7 +3266,6 @@ void encodeGDS(GRIBPACK *lGrib, long *gribLen, int *isec2, double *fsec2)
*gribLen = z;
}
/* GRIB BLOCK 3 - BIT MAP SECTION */
static
void encodeBMS(GRIBPACK *lGrib, long *gribLen, double *fsec3, int *isec4, double *data, long *datasize)
......@@ -3345,7 +3359,7 @@ void encodeBMS(GRIBPACK *lGrib, long *gribLen, double *fsec3, int *isec4, double
}
static
void encode_double_array_common(int NumBits, long PackStart, long datasize, GRIBPACK *lGrib,
void encode_double_array_common(int numBits, long packStart, long datasize, GRIBPACK *lGrib,
const double *data, double zref, double factor, long *gz)
{
long i, z = *gz;
......@@ -3358,7 +3372,7 @@ void encode_double_array_common(int NumBits, long PackStart, long datasize, GRIB
cbits = 8;
c = 0;
for ( i = PackStart; i < datasize; i++ )
for ( i = packStart; i < datasize; i++ )
{
/* note float -> unsigned int .. truncate */
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
......@@ -3366,7 +3380,7 @@ void encode_double_array_common(int NumBits, long PackStart, long datasize, GRIB
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
jbits = NumBits;
jbits = numBits;
while ( cbits <= jbits )
{
if ( cbits == 8 )
......@@ -3395,14 +3409,14 @@ void encode_double_array_common(int NumBits, long PackStart, long datasize, GRIB
}
static
void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *lGrib,
void encode_double_array(int numBits, long PackStart, long datasize, GRIBPACK *lGrib,
const double *data, double zref, double factor, long *gz)
{
static const char *func = "encode_double_array";
long i, z = *gz;
unsigned int ipval;
if ( NumBits == 8 )
if ( numBits == 8 )
{
#if defined (CRAY)
#pragma _CRI ivdep
......@@ -3422,7 +3436,7 @@ void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *l
z++;
}
}
else if ( NumBits == 16 )
else if ( numBits == 16 )
{
#if defined (CRAY)
#pragma _CRI ivdep
......@@ -3443,7 +3457,7 @@ void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *l
z += 2;
}
}
else if ( NumBits == 24 )
else if ( numBits == 24 )
{
#if defined (CRAY)
#pragma _CRI ivdep
......@@ -3465,7 +3479,7 @@ void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *l
z += 3;
}
}
else if ( NumBits == 32 )
else if ( numBits == 32 )
{
#if defined (CRAY)
#pragma _CRI ivdep
......@@ -3488,17 +3502,17 @@ void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *l
z += 4;
}
}
else if ( NumBits > 0 && NumBits <= 32 )
else if ( numBits > 0 && numBits <= 32 )
{
encode_double_array_common(NumBits, PackStart, datasize, lGrib,
encode_double_array_common(numBits, PackStart, datasize, lGrib,
data, zref, factor, &z);
}
else if ( NumBits == 0 )
else if ( numBits == 0 )
{
}
else
{
Error(func, "Unimplemented packing factor %d!", NumBits);
Error(func, "Unimplemented packing factor %d!", numBits);
}
*gz = z;
......@@ -3506,7 +3520,7 @@ void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *l
/* GRIB BLOCK 4 - BINARY DATA SECTION */
static
int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *isec4, long datasize, double data[],
int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *isec4, long datasize, double *data,
long *datstart, long *datsize, int code)
{
/* Uwe Schulzweida, 11/04/2003 : Check that number of bits per value is not exceeded */
......@@ -3515,6 +3529,7 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
static const char *func = "encodeBDS";
long z = *gribLen;
long i, jloop;
int numBits;
int ival;
int blockLength, PackStart = 0, Flag = 0;
int binscale = 0;
......@@ -3533,6 +3548,7 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
/* of floating point numbers - needed */
/* on some platforms (eg vpp700, linux) */
extern double _pow2tab[158];
extern int GRB_Const; /* 1: Don't pack constant fields on regular grids */
if ( isec2 )
{
......@@ -3614,12 +3630,7 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
*datstart = bds_head + bds_ext;
blockLength = (*datstart) + (ISEC4_NumBits*(datasize - PackStart) + 7)/8;
if ( (blockLength%2) == 1 ) blockLength++;
unused_bits = blockLength*8 - (*datstart)*8 - ISEC4_NumBits*(datasize - PackStart);
Flag += unused_bits;
numBits = ISEC4_NumBits;
if ( lspherc && lcomplex )
{
......@@ -3651,11 +3662,26 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
zref = fmin;
if ( GRB_Const && !lspherc )
{
if ( IS_EQUAL(fmin, fmax) ) numBits = 0;
}
blockLength = (*datstart) + (numBits*(datasize - PackStart) + 7)/8;
if ( (blockLength%2) == 1 ) blockLength++;
unused_bits = blockLength*8 - (*datstart)*8 - numBits*(datasize - PackStart);
Flag += unused_bits;
/*
Adjust number of bits per value if full integer length to
avoid hitting most significant bit (sign bit).
*/
nbpv = ISEC4_NumBits;
nbpv = numBits;
/* if( nbpv == ibits ) nbpv = nbpv - 1; */
/*
Calculate the binary scaling factor to spread the range of
......@@ -3768,9 +3794,9 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
}
}
*datsize = ((datasize-PackStart)*ISEC4_NumBits + 7)/8;
*datsize = ((datasize-PackStart)*numBits + 7)/8;
encode_double_array(ISEC4_NumBits, PackStart, datasize, lGrib,
encode_double_array(numBits, PackStart, datasize, lGrib,
data, zref, factor, &z);
if ( unused_bits >= 8 )
......@@ -4168,7 +4194,7 @@ int decodePDS(unsigned char *pds, int *isec0, int *isec1)
}
static
int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int dfunc)
int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int *numGridVals, int dfunc)
{
static const char *func = "decodeGDS";
int imisng;
......@@ -4185,6 +4211,8 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
size_t lGribLen = 0;
#endif
*numGridVals = 0;
imisng = 0;
memset(isec2, 0, 22*sizeof(int));
......@@ -4225,9 +4253,13 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
jlenl = (gdsLen - locnl) >> 1;
if ( jlenl == GDS_NumLat )
{
*numGridVals = 0;
ISEC2_Reduced = TRUE;
for ( i = 0; i < jlenl; i++ )
ISEC2_RowLon(i) = GET_UINT2(gds[locnl+2*i], gds[locnl+2*i+1]);
{
ISEC2_RowLon(i) = GET_UINT2(gds[locnl+2*i], gds[locnl+2*i+1]);
*numGridVals += ISEC2_RowLon(i);
}
}
else
{
......@@ -4244,8 +4276,12 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
ISEC2_GridType == GTYPE_GAUSSIAN ||
ISEC2_GridType == GTYPE_LATLON_ROT )
{
if ( ! ReducedGrid ) ISEC2_NumLon = GDS_NumLon;
ISEC2_NumLat = GDS_NumLat;
if ( ! ReducedGrid )
{
ISEC2_NumLon = GDS_NumLon;
*numGridVals = ISEC2_NumLon*ISEC2_NumLat;
}
ISEC2_FirstLat = GDS_FirstLat;
ISEC2_FirstLon = GDS_FirstLon;
ISEC2_ResFlag = GDS_ResFlag;
......@@ -4288,6 +4324,7 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
{
ISEC2_NumLon = GDS_NumLon;
ISEC2_NumLat = GDS_NumLat;
*numGridVals = ISEC2_NumLon*ISEC2_NumLat;
ISEC2_FirstLat = GDS_FirstLat;
ISEC2_FirstLon = GDS_FirstLon;
ISEC2_ResFlag = GDS_ResFlag;
......@@ -4308,6 +4345,7 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
ISEC2_PentaM = GDS_PentaM;
ISEC2_RepType = GDS_RepType;
ISEC2_RepMode = GDS_RepMode;
*numGridVals = (ISEC2_PentaJ+1)*(ISEC2_PentaJ+2);
isec2[ 6] = 0;
isec2[ 7] = 0;
isec2[ 8] = 0;
......@@ -4328,6 +4366,7 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
ISEC2_GME_LonPP = GDS_GME_LonPP;
ISEC2_GME_LonMPL = GDS_GME_LonMPL;
ISEC2_GME_BFlag = GDS_GME_BFlag;
*numGridVals = (ISEC2_GME_NI+1)*(ISEC2_GME_NI+1)*10;
/*
iret = decodeGDS_TR(gds, gdspos, isec0, isec2, imisng);
*/
......@@ -4336,6 +4375,7 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int df
{
ISEC2_NumLon = GDS_NumLon;
ISEC2_NumLat = GDS_NumLat;
*numGridVals = ISEC2_NumLon*ISEC2_NumLat;
Message(func, "Gridtype %d unsupported", ISEC2_GridType);
}
......@@ -4419,7 +4459,7 @@ void decode_double_array_common(unsigned char *igrib, long jlend, int NumBits,
}
static
void decode_double_array(unsigned char *igrib, long jlend, long jlenc, int NumBits,
void decode_double_array(unsigned char *igrib, long jlend, int numBits,
double fmin, double zscale, double *fpdata)
{
long i;
......@@ -4428,9 +4468,10 @@ void decode_double_array(unsigned char *igrib, long jlend, long jlenc, int NumBi
static const char *func = "decode_double_array";
GRIBPACK *lgrib = NULL;
if ( NumBits == 8 || NumBits == 16 ||
NumBits == 24 || NumBits == 32 )
if ( numBits == 8 || numBits == 16 ||
numBits == 24 || numBits == 32 )
{
long jlenc = jlend * numBits / 8;
if ( jlenc > 0 )
{
lgrib = (GRIBPACK *) malloc(jlenc*sizeof(GRIBPACK));
......@@ -4440,88 +4481,88 @@ void decode_double_array(unsigned char *igrib, long jlend, long jlenc, int NumBi
}
}
if ( NumBits == 0 )
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( NumBits == 8 )
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)lgrib[i];
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 16 )
else if ( numBits == 16 )
for ( i = 0; i < jlend; i++ )
{
dval = (((int)lgrib[2*i ] << 8) + (int)lgrib[2*i+1]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 24 )
else if ( numBits == 24 )
for ( i = 0; i < jlend; i++ )
{
dval = (((int)lgrib[3*i ] << 16) + ((int)lgrib[3*i+1] << 8) +
(int)lgrib[3*i+2]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 32 )
else if ( numBits == 32 )
for ( i = 0; i < jlend; i++ )
{
dval = (((unsigned int)lgrib[4*i ] << 24) + ((unsigned int)lgrib[4*i+1] << 16) +
((unsigned int)lgrib[4*i+2] << 8) + (unsigned int)lgrib[4*i+3]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits <= 25 )
else if ( numBits <= 25 )
{
decode_double_array_common(igrib, jlend, NumBits, fmin, zscale, fpdata);
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
fprintf(stderr," Unimplemented packing factor %d\n", NumBits);
fprintf(stderr," Unimplemented packing factor %d\n", numBits);
exit(EXIT_FAILURE);
}
if ( lgrib ) free(lgrib);
#else
if ( NumBits == 0 )
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( NumBits == 8 )
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)igrib[i];
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 16 )
else if ( numBits == 16 )
for ( i = 0; i < jlend; i++ )
{
dval = (((int)igrib[2*i ] << 8) | (int)igrib[2*i+1]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 24 )
else if ( numBits == 24 )
for ( i = 0; i < jlend; i++ )
{
dval = (((int)igrib[3*i ] << 16) + ((int)igrib[3*i+1] << 8) +
(int)igrib[3*i+2]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits == 32 )
else if ( numBits == 32 )
for ( i = 0; i < jlend; i++ )
{
dval = (((unsigned int)igrib[4*i ] << 24) + ((unsigned int)igrib[4*i+1] << 16) +
((unsigned int)igrib[4*i+2] << 8) + (unsigned int)igrib[4*i+3]);
fpdata[i] = fmin + zscale * dval;
}
else if ( NumBits <= 25 )
else if ( numBits <= 25 )
{
decode_double_array_common(igrib, jlend, NumBits, fmin, zscale, fpdata);
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
fprintf(stderr," Unimplemented packing factor %d\n", NumBits);
fprintf(stderr," Unimplemented packing factor %d\n", numBits);
exit(EXIT_FAILURE);
}
#endif
......@@ -4529,7 +4570,7 @@ void decode_double_array(unsigned char *igrib, long jlend, long jlenc, int NumBi
static
int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
double *fsec4, int fsec4len, int dfunc, int bdsLenIn, int llarge, int *iret)
double *fsec4, int fsec4len, int dfunc, int bdsLenIn, int numGridVals, int llarge, int *iret)
{
static const char *func = "decodeBDS";
unsigned char *igrib;
......@@ -4537,7 +4578,7 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
int lcompress;
int jup, kup, mup;
int locnd;
int jlend, jlenc;
int jlend;
long i;
int bds_flag, jscale, imiss;
int bds_ubits;
......@@ -4706,15 +4747,31 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
/* simple (lcomplex = 0) or complex (lcomplex = 1) */
jlend = bdsLen - locnd;
if ( ISEC4_NumBits == 0 )
{
if ( jlend > 1 ) Error(func, "Number of bits per data value = 0");
if ( jlend > 1 )
{
*iret = 2001;
gprintf(func, " Number of bits per data value = 0!");
gprintf(func, " Return code = %d", *iret);
return (0);
}
if ( numGridVals == 0 )
{
*iret = 2002;
gprintf(func, " Constant field unsupported for this grid type!");
gprintf(func, " Return code = %d", *iret);
return (0);
}
jlend = ISEC2_NumLon*ISEC2_NumLat;
jlend = numGridVals;
jlend -= ioff;
}
else
{
jlend = (jlend*8 - bds_ubits) / ISEC4_NumBits;
jlend = (jlend*8 - bds_ubits) / ISEC4_NumBits;
}
ISEC4_NumValues = jlend + ioff;
......@@ -4757,9 +4814,8 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
else
{
igrib += locnd;
jlenc = jlend * ISEC4_NumBits / 8;
decode_double_array(igrib, jlend, jlenc, ISEC4_NumBits, fmin, zscale, fpdata);
decode_double_array(igrib, jlend, ISEC4_NumBits, fmin, zscale, fpdata);
}
if ( lspherc && lcomplex )
......@@ -4807,6 +4863,7 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
int ldebug = FALSE;
int llarge = FALSE, l_iorj = FALSE;
int lsect2 = FALSE, lsect3 = FALSE;
int numGridVals = 0;
static int lmissvalinfo = 1;
*iret = 0;
......@@ -4950,7 +5007,7 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
{
gds = is + isLen + pdsLen;
gdsLen = decodeGDS(gds, isec0, isec2, fsec2, dfunc);
gdsLen = decodeGDS(gds, isec0, isec2, fsec2, &numGridVals, dfunc);
}
/*
......@@ -4983,7 +5040,7 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
bdsLen = ISEC0_GRIB_Len - (isLen + pdsLen + gdsLen + bmsLen);
bdsLen = decodeBDS(ISEC1_DecScaleFactor, bds, isec2, isec4,
fsec4, fsec4len, dfunc, bdsLen, llarge, iret);
fsec4, fsec4len, dfunc, bdsLen, numGridVals, llarge, iret);
if ( *iret != 0 ) return;
......@@ -9061,7 +9118,7 @@ int gribUnzip(unsigned char *dbuf, long dbufsize, unsigned char *sbuf, long sbu
return (gribLen);
}
static const char grb_libvers[] = "1.4.6" " of ""Aug 23 2010"" ""09:40:05";
static const char grb_libvers[] = "1.4.6" " of ""Sep 9 2010"" ""10:26:38";
const char *
cgribexLibraryVersion(void)
{
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment