From a6dec296fb96b4390f92e71aef4934e6704068a4 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Tue, 29 Jan 2019 11:04:49 +0100
Subject: [PATCH] cgribexlib.c update.
---
src/cgribexlib.c | 143 ++++++++++++++++-------------------------------
1 file changed, 47 insertions(+), 96 deletions(-)
diff --git a/src/cgribexlib.c b/src/cgribexlib.c
index aca6de90d..7f2128f7d 100644
--- a/src/cgribexlib.c
+++ b/src/cgribexlib.c
@@ -1,5 +1,5 @@
-/* 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));
--
GitLab