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

cgribexlib.c update.

parent abd0bd61
/* Automatically generated by m214003 at 2018-11-14, do not edit */
/* Automatically generated by m214003 at 2019-01-29, do not edit */
/* CGRIBEXLIB_VERSION="1.9.2" */
......@@ -1665,14 +1665,14 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
/* Section 1 . Initialise */
/* ----------------------------------------------------------------- */
/* Check conversion type parameter. */
// Check conversion type parameter.
int iround = kround;
if ( iround != 0 && iround != 1 )
{
Error("Invalid conversion type = %d", iround);
/* If not aborting, arbitrarily set rounding to 'up'. */
// If not aborting, arbitrarily set rounding to 'up'.
iround = 1;
}
......@@ -1684,8 +1684,6 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
{
*kexp = 0;
*kmant = 0;
// iexp = 0;
// isign = 0;
goto LABEL900;
}
......@@ -1693,57 +1691,43 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
/* Section 3 . Convert other values. */
/* ----------------------------------------------------------------- */
{
double zeps = kbits != 32 ? 1.0e-12 : 1.0e-8;
const double zeps = (kbits != 32) ? 1.0e-12 : 1.0e-8;
double zref = pval;
/* Sign of value. */
int isign = zref >= 0.0 ? 0 : 128;
// Sign of value.
const int isign = (zref >= 0.0) ? 0 : 128;
zref = fabs(zref);
/* Exponent. */
// Exponent.
int iexp = (int) (log(zref)/log(16.0) + 65.0 + zeps);
/* only ANSI C99 has log2 */
/* iexp = (int) (log2(zref) * 0.25 + 65.0 + zeps); */
// only ANSI C99 has log2
// iexp = (int) (log2(zref) * 0.25 + 65.0 + zeps);
if ( iexp < 0 ) iexp = 0;
if ( iexp > 127 ) iexp = 127;
double rpowref;
/*
rpowref = zref / pow(16.0, (double)(iexp - 70));
*/
rpowref = ldexp(zref, 4 * -(iexp - 70));
/* Mantissa. */
// double rpowref = zref / pow(16.0, (double)(iexp - 70));
double rpowref = ldexp(zref, 4 * -(iexp - 70));
// Mantissa.
if ( iround == 0 )
{
/* Closest number in GRIB format less than original number. */
/* Truncate for positive numbers. */
/* Round up for negative numbers. */
if ( isign == 0 )
*kmant = (int)rpowref;
else
*kmant = (int)lround(rpowref + 0.5);
*kmant = (isign == 0) ? (int)rpowref : (int)lround(rpowref + 0.5);
}
else
{
/* Closest number in GRIB format to the original number */
/* (equal to, greater than or less than original number). */
*kmant = (int)lround(rpowref);
}
/* Check that mantissa value does not exceed 24 bits. */
/* If it does, adjust the exponent upwards and recalculate */
/* the mantissa. */
/* If it does, adjust the exponent upwards and recalculate the mantissa. */
/* 16777215 = 2**24 - 1 */
if ( *kmant > 16777215 )
{
......@@ -1751,24 +1735,19 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
++iexp;
/* Check for exponent overflow during adjustment */
// Check for exponent overflow during adjustment
if ( iexp > 127 )
{
Message("Exponent overflow");
Message("Original number = %30.20f", pval);
Message("Sign = %3d, Exponent = %3d, Mantissa = %12d",
isign, iexp, *kmant);
Message("Sign = %3d, Exponent = %3d, Mantissa = %12d", isign, iexp, *kmant);
Error("Exponent overflow");
/* If not aborting, arbitrarily set value to zero */
// If not aborting, arbitrarily set value to zero
Message("Value arbitrarily set to zero.");
*kexp = 0;
*kmant = 0;
// iexp = 0;
// isign = 0;
goto LABEL900;
}
......@@ -1779,28 +1758,20 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
/* Closest number in GRIB format less than original number. */
/* Truncate for positive numbers. */
/* Round up for negative numbers. */
if ( isign == 0 )
*kmant = (int)rpowref;
else
*kmant = (int)lround(rpowref + 0.5);
*kmant = (isign == 0) ? (int)rpowref : (int)lround(rpowref + 0.5);
}
else
{
/* Closest number in GRIB format to the original number */
/* (equal to, greater or less than original number). */
*kmant = (int)lround(rpowref);
}
/* Repeat calculation (with modified exponent) if still have */
/* mantissa overflow. */
// Repeat calculation (with modified exponent) if still have mantissa overflow.
if ( *kmant > 16777215 ) goto LABEL350;
}
/* Add sign bit to exponent. */
// Add sign bit to exponent.
*kexp = iexp + isign;
}
......@@ -1907,8 +1878,8 @@ double decfp2(int kexp, int kmant)
/* Sign of value. */
int iexp = kexp,
isign = (iexp < 128) * 2 - 1;
int iexp = kexp;
const int isign = (iexp < 128) * 2 - 1;
iexp -= iexp < 128 ? 0 : 128;
......@@ -1918,7 +1889,7 @@ double decfp2(int kexp, int kmant)
iexp -= 64;
double pval = ldexp(1.0, 4 * iexp) * isign * POW_2_M24 * kmant;
const double pval = ldexp(1.0, 4 * iexp) * isign * POW_2_M24 * kmant;
/* ----------------------------------------------------------------- */
/* Section 9. Return to calling routine. */
......@@ -6775,9 +6746,7 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, int inv)
{
double power;
double *scale = (double*) Malloc(((size_t)trunc+1)*sizeof(double));
if ( scale == NULL ) SysError("No Memory!");
if ( pcScale < -10000 || pcScale > 10000 )
......@@ -6790,7 +6759,7 @@ void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, i
if ( pcScale == 0 ) return;
power = (double) pcScale / 1000.;
const double power = (double) pcScale / 1000.;
scale[0] = 1.0;
if (pcScale != 1000)
......@@ -7294,15 +7263,10 @@ C
/* Local variables */
int ilii, ilio, icode;
int iregno, iquano, j210, j220, j230, j240, j225;
T *ztemp = NULL;
T *zline = NULL;
T *zwork = NULL;
ztemp = (T*) Malloc((size_t)klon*(size_t)klat*sizeof(T));
zline = (T*) Malloc(2*(size_t)klon*sizeof(T));
zwork = (T*) Malloc(3*(2*(size_t)klon+3)*sizeof(T));
T *ztemp = (T*) Malloc((size_t)klon*(size_t)klat*sizeof(T));
T *zline = (T*) Malloc(2*(size_t)klon*sizeof(T));
T *zwork = (T*) Malloc(3*(2*(size_t)klon+3)*sizeof(T));
/* Parameter adjustments */
--pfield;
......@@ -7536,9 +7500,7 @@ double TEMPLATE(calculate_pfactor,T)(const T *spectralField, long fieldTruncatio
void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, int inv)
{
double power;
double *scale = (double*) Malloc(((size_t)trunc+1)*sizeof(double));
if ( scale == NULL ) SysError("No Memory!");
if ( pcScale < -10000 || pcScale > 10000 )
......@@ -7551,7 +7513,7 @@ void TEMPLATE(scale_complex,T)(T *fpdata, int pcStart, int pcScale, int trunc, i
if ( pcScale == 0 ) return;
power = (double) pcScale / 1000.;
const double power = (double) pcScale / 1000.;
scale[0] = 1.0;
if (pcScale != 1000)
......@@ -8055,15 +8017,10 @@ C
/* Local variables */
int ilii, ilio, icode;
int iregno, iquano, j210, j220, j230, j240, j225;
T *ztemp = NULL;
T *zline = NULL;
T *zwork = NULL;
ztemp = (T*) Malloc((size_t)klon*(size_t)klat*sizeof(T));
zline = (T*) Malloc(2*(size_t)klon*sizeof(T));
zwork = (T*) Malloc(3*(2*(size_t)klon+3)*sizeof(T));
T *ztemp = (T*) Malloc((size_t)klon*(size_t)klat*sizeof(T));
T *zline = (T*) Malloc(2*(size_t)klon*sizeof(T));
T *zwork = (T*) Malloc(3*(2*(size_t)klon+3)*sizeof(T));
/* Parameter adjustments */
--pfield;
......@@ -8447,13 +8404,12 @@ void TEMPLATE(decode_array_common,T)(const unsigned char *restrict igrib, long j
{
/* code from wgrib routine BDS_unpack */
const unsigned char *bits = igrib;
long i;
unsigned int tbits = 0;
int n_bits = NumBits;
int t_bits = 0;
unsigned jmask = (1U << n_bits) - 1U;
for ( i = 0; i < jlend; i++ )
const unsigned jmask = (1U << n_bits) - 1U;
for (long i = 0; i < jlend; i++)
{
if (n_bits - t_bits > 8)
{
......@@ -8471,7 +8427,7 @@ void TEMPLATE(decode_array_common,T)(const unsigned char *restrict igrib, long j
fpdata[i] = (float)((tbits >> t_bits) & jmask);
}
/* at least this vectorizes :) */
for ( i = 0; i < jlend; i++ )
for (long i = 0; i < jlend; i++)
fpdata[i] = fmin + zscale*fpdata[i];
}
......@@ -8480,18 +8436,16 @@ void TEMPLATE(decode_array_common2,T)(const unsigned char *restrict igrib, long
T fmin, T zscale, T *restrict fpdata)
{
static const unsigned mask[] = {0,1,3,7,15,31,63,127,255};
static const double shift[9]
= {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
static const double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
/* code from wgrib routine BDS_unpack */
const unsigned char *bits = igrib;
long i;
int n_bits = NumBits;
int c_bits, j_bits;
/* older unoptimized code, not often used */
c_bits = 8;
for ( i = 0; i < jlend; i++ )
for (long i = 0; i < jlend; i++)
{
double jj = 0.0;
j_bits = n_bits;
......@@ -8682,13 +8636,12 @@ void TEMPLATE(decode_array_common,T)(const unsigned char *restrict igrib, long j
{
/* code from wgrib routine BDS_unpack */
const unsigned char *bits = igrib;
long i;
unsigned int tbits = 0;
int n_bits = NumBits;
int t_bits = 0;
unsigned jmask = (1U << n_bits) - 1U;
for ( i = 0; i < jlend; i++ )
const unsigned jmask = (1U << n_bits) - 1U;
for (long i = 0; i < jlend; i++)
{
if (n_bits - t_bits > 8)
{
......@@ -8706,7 +8659,7 @@ void TEMPLATE(decode_array_common,T)(const unsigned char *restrict igrib, long j
fpdata[i] = (float)((tbits >> t_bits) & jmask);
}
/* at least this vectorizes :) */
for ( i = 0; i < jlend; i++ )
for (long i = 0; i < jlend; i++)
fpdata[i] = fmin + zscale*fpdata[i];
}
......@@ -8715,18 +8668,16 @@ void TEMPLATE(decode_array_common2,T)(const unsigned char *restrict igrib, long
T fmin, T zscale, T *restrict fpdata)
{
static const unsigned mask[] = {0,1,3,7,15,31,63,127,255};
static const double shift[9]
= {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
static const double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
/* code from wgrib routine BDS_unpack */
const unsigned char *bits = igrib;
long i;
int n_bits = NumBits;
int c_bits, j_bits;
/* older unoptimized code, not often used */
c_bits = 8;
for ( i = 0; i < jlend; i++ )
for (long i = 0; i < jlend; i++)
{
double jj = 0.0;
j_bits = n_bits;
......@@ -10888,17 +10839,17 @@ static
void TEMPLATE(encode_array_common,T)(int numBits, size_t packStart, size_t datasize, GRIBPACK *lGrib,
const T *data, T zref, T factor, size_t *gz)
{
size_t i, z = *gz;
size_t z = *gz;
unsigned int ival;
int cbits, jbits;
unsigned int c;
static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
static const unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
/* code from gribw routine flist2bitstream */
cbits = 8;
c = 0;
for ( i = packStart; i < datasize; i++ )
for (size_t i = packStart; i < datasize; i++)
{
/* note float -> unsigned int .. truncate */
// ival = (unsigned int)(TEMPLATE(round,T)((data[i] - zref) * factor));
......@@ -11451,17 +11402,17 @@ static
void TEMPLATE(encode_array_common,T)(int numBits, size_t packStart, size_t datasize, GRIBPACK *lGrib,
const T *data, T zref, T factor, size_t *gz)
{
size_t i, z = *gz;
size_t z = *gz;
unsigned int ival;
int cbits, jbits;
unsigned int c;
static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
static const unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
/* code from gribw routine flist2bitstream */
cbits = 8;
c = 0;
for ( i = packStart; i < datasize; i++ )
for (size_t i = packStart; i < datasize; i++)
{
/* note float -> unsigned int .. truncate */
// ival = (unsigned int)(TEMPLATE(round,T)((data[i] - zref) * factor));
......
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