Skip to content

GitLab

Commit 295be48d authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

cgribexlib update

parent d4817419
Hide whitespace changes
Inline Side-by-side
Showing with 220 additions and 203 deletions
ChangeLog
View file @ 295be48d
2013-03-?? Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* Version 1.6.0 released
* using CGRIBEX library version 1.6.0
2013-01-08 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* global netCDF attribute "source" missing (bug fix)
......
src/cgribexlib.c
View file @ 295be48d
/* Automatically generated by m214003 at 2012-12-17, do not edit */
/* Automatically generated by m214003 at 2013-01-14, do not edit */
/* CGRIBEXLIB_VERSION="1.5.6" */
/* CGRIBEXLIB_VERSION="1.6.0" */
#ifdef _ARCH_PWR6
#pragma options nostrict
......@@ -63,6 +63,7 @@
# define VECTORCODE
#endif
#if defined (VECTORCODE)
#if defined (INT32)
# define GRIBPACK unsigned INT32
......@@ -1698,6 +1699,208 @@ void encode_double_array_unrolled(int numBits, size_t packStart, size_t datasize
#undef __UNROLL_DEPTH_2
}
static
void decode_double_array_common(unsigned char *igrib, long jlend, int NumBits,
double fmin, double zscale, double *fpdata)
{
/* code from wgrib routine BDS_unpack */
unsigned char *bits = igrib;
unsigned int jmask;
long i;
unsigned int tbits = 0;
int n_bits = NumBits;
int t_bits = 0;
jmask = (1 << n_bits) - 1;
for ( i = 0; i < jlend; i++ )
{
if (n_bits - t_bits > 8)
{
tbits = (tbits << 16) | (bits[0] << 8) | (bits[1]);
bits += 2;
t_bits += 16;
}
while ( t_bits < n_bits )
{
tbits = (tbits * 256) + *bits++;
t_bits += 8;
}
t_bits -= n_bits;
fpdata[i] = (tbits >> t_bits) & jmask;
}
/* at least this vectorizes :) */
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin + zscale*fpdata[i];
}
static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
static double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
static
void decode_double_array_common2(unsigned char *igrib, long jlend, int NumBits,
double fmin, double zscale, double *fpdata)
{
/* code from wgrib routine BDS_unpack */
unsigned char *bits = igrib;
long i;
int n_bits = NumBits;
int c_bits, j_bits;
double jj;
/* older unoptimized code, not often used */
c_bits = 8;
for ( i = 0; i < jlend; i++ )
{
jj = 0.0;
j_bits = n_bits;
while (c_bits <= j_bits)
{
if (c_bits == 8)
{
jj = jj * 256.0 + (double) (*bits++);
j_bits -= 8;
}
else
{
jj = (jj * shift[c_bits]) + (double) (*bits & mask[c_bits]);
bits++;
j_bits -= c_bits;
c_bits = 8;
}
}
if (j_bits)
{
c_bits -= j_bits;
jj = (jj * shift[j_bits]) + (double) ((*bits >> c_bits) & mask[j_bits]);
}
fpdata[i] = fmin + zscale*jj;
}
}
static
void decode_double_array_byte(unsigned char *igrib, long jlend, int numBits,
double fmin, double zscale, double *fpdata)
{
long i;
double dval;
#if defined (VECTORCODE)
GRIBPACK *lgrib = NULL;
if ( numBits%8 == 0 )
{
long jlenc = jlend * numBits / 8;
if ( jlenc > 0 )
{
lgrib = (GRIBPACK *) malloc(jlenc*sizeof(GRIBPACK));
if ( lgrib == NULL ) SysError("No Memory!");
(void) UNPACK_GRIB(igrib, lgrib, jlenc, -1L);
}
}
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)lgrib[i];
fpdata[i] = fmin + zscale * dval;
}
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 )
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 )
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 )
{
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else if ( numBits > 25 && numBits < 32 )
{
decode_double_array_common2(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
Error("Unimplemented packing factor %d!", numBits);
}
if ( lgrib ) free(lgrib);
#else
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)igrib[i];
fpdata[i] = fmin + zscale * dval;
}
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 )
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 )
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 )
{
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else if ( numBits > 25 && numBits < 32 )
{
decode_double_array_common2(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
Error("Unimplemented packing factor %d!", numBits);
}
#endif
}
static
void decode_double_array_unrolled(unsigned char *igrib, long jlend, int numBits,
double fmin, double zscale, double *fpdata)
{
decode_double_array_byte(igrib, jlend, numBits, fmin, zscale, fpdata);
}
#define NINT(x) ((x) < 0 ? (int)((x)-.5) : (int)((x)+.5))
......@@ -4856,6 +5059,8 @@ void gribEncode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
int gribVersion(unsigned char *is, size_t buffersize)
{
if ( buffersize < 8 )
......@@ -5370,204 +5575,6 @@ int decodeGDS(unsigned char *gds, int *isec0, int *isec2, double *fsec2, int *n
return (gdsLen);
}
static
void decode_double_array_common(unsigned char *igrib, long jlend, int NumBits,
double fmin, double zscale, double *fpdata)
{
/* code from wgrib routine BDS_unpack */
unsigned char *bits = igrib;
unsigned int jmask;
long i;
unsigned int tbits = 0;
int n_bits = NumBits;
int t_bits = 0;
jmask = (1 << n_bits) - 1;
for ( i = 0; i < jlend; i++ )
{
if (n_bits - t_bits > 8)
{
tbits = (tbits << 16) | (bits[0] << 8) | (bits[1]);
bits += 2;
t_bits += 16;
}
while ( t_bits < n_bits )
{
tbits = (tbits * 256) + *bits++;
t_bits += 8;
}
t_bits -= n_bits;
fpdata[i] = (tbits >> t_bits) & jmask;
}
/* at least this vectorizes :) */
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin + zscale*fpdata[i];
}
static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
static double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
static
void decode_double_array_common2(unsigned char *igrib, long jlend, int NumBits,
double fmin, double zscale, double *fpdata)
{
/* code from wgrib routine BDS_unpack */
unsigned char *bits = igrib;
long i;
int n_bits = NumBits;
int c_bits, j_bits;
double jj;
/* older unoptimized code, not often used */
c_bits = 8;
for ( i = 0; i < jlend; i++ )
{
jj = 0.0;
j_bits = n_bits;
while (c_bits <= j_bits)
{
if (c_bits == 8)
{
jj = jj * 256.0 + (double) (*bits++);
j_bits -= 8;
}
else
{
jj = (jj * shift[c_bits]) + (double) (*bits & mask[c_bits]);
bits++;
j_bits -= c_bits;
c_bits = 8;
}
}
if (j_bits)
{
c_bits -= j_bits;
jj = (jj * shift[j_bits]) + (double) ((*bits >> c_bits) & mask[j_bits]);
}
fpdata[i] = fmin + zscale*jj;
}
}
static
void decode_double_array(unsigned char *igrib, long jlend, int numBits,
double fmin, double zscale, double *fpdata)
{
long i;
double dval;
#if defined (VECTORCODE)
GRIBPACK *lgrib = NULL;
if ( numBits == 8 || numBits == 16 ||
numBits == 24 || numBits == 32 )
{
long jlenc = jlend * numBits / 8;
if ( jlenc > 0 )
{
lgrib = (GRIBPACK *) malloc(jlenc*sizeof(GRIBPACK));
if ( lgrib == NULL ) SysError("No Memory!");
(void) UNPACK_GRIB(igrib, lgrib, jlenc, -1L);
}
}
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)lgrib[i];
fpdata[i] = fmin + zscale * dval;
}
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 )
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 )
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 )
{
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else if ( numBits > 25 && numBits < 32 )
{
decode_double_array_common2(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
fprintf(stderr," Unimplemented packing factor %d!\n", numBits);
exit(EXIT_FAILURE);
}
if ( lgrib ) free(lgrib);
#else
if ( numBits == 0 )
{
for ( i = 0; i < jlend; i++ )
fpdata[i] = fmin;
}
else if ( numBits == 8 )
for ( i = 0; i < jlend; i++ )
{
dval = (int)igrib[i];
fpdata[i] = fmin + zscale * dval;
}
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 )
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 )
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 )
{
decode_double_array_common(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else if ( numBits > 25 && numBits < 32 )
{
decode_double_array_common2(igrib, jlend, numBits, fmin, zscale, fpdata);
}
else
{
fprintf(stderr, "Unimplemented packing factor %d!\n", numBits);
exit(EXIT_FAILURE);
}
#endif
}
static
int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
double *fsec4, int fsec4len, int dfunc, int bdsLenIn, int numGridVals, int llarge, int *iret)
......@@ -5814,7 +5821,11 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
{
igrib += locnd;
decode_double_array(igrib, jlend, ISEC4_NumBits, fmin, zscale, fpdata);
#if defined (_ARCH_PWR6)
decode_double_array_unrolled(igrib, jlend, ISEC4_NumBits, fmin, zscale, fpdata);
#else
decode_double_array_byte (igrib, jlend, ISEC4_NumBits, fmin, zscale, fpdata);
#endif
}
if ( lspherc && lcomplex )
......@@ -6155,7 +6166,8 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
}
if ( dfunc == 'R' && *iret == -801 )
gprintf(__func__, "Number of values (%d) and sum of lons per row (%d) differ, abort conversion to regular grid!", ISEC4_NumValues, nvalues);
gprintf(__func__, "Number of values (%d) and sum of lons per row (%d) differ, abort conversion to regular grid!",
ISEC4_NumValues, nvalues);
if ( dfunc == 'R' && *iret != -801 )
{
......@@ -10176,7 +10188,7 @@ int gribUnzip(unsigned char *dbuf, long dbufsize, unsigned char *sbuf, long sbu
return (gribLen);
}
static const char grb_libvers[] = "1.5.6" " of ""Dec 17 2012"" ""13:44:05";
static const char grb_libvers[] = "1.6.0" " of ""Jan 14 2013"" ""11:39:57";
const char *
cgribexLibraryVersion(void)
{
......
Markdown is supported
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