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

cdilib update

parent 1e7e7da0
......@@ -231,5 +231,7 @@ int gribVersion(unsigned char *buffer, size_t buffersize);
int gribGinfo(long recpos, long recsize, unsigned char *gribbuffer, int *intnum, float *fltnum);
double calculate_pfactor(const double* spectralField, long fieldTruncation, long subsetTruncation);
#endif /* _CGRIBEX_H */
/* Automatically generated by m214003 at 2010-02-08, do not edit */
/* Automatically generated by m214003 at 2010-02-15, do not edit */
/* CGRIBEXLIB_VERSION="1.4.2" */
/* CGRIBEXLIB_VERSION="1.4.3" */
#if defined (HAVE_CONFIG_H)
# include "config.h"
......@@ -138,8 +138,10 @@ void confp3(double pval, int *kexp, int *kmant, int kbits, int kround);
double decfp2(int kexp, int kmant);
void ref2ibm(double *pref, int kbits);
void scaleComplex(double *fpdata, int pcStart, int pcScale, int truncation);
void scatterComplex(double *fpdata, int pcStart, int truncation, int dimSP);
void scaleComplex(double *fpdata, int pcStart, int pcScale, int trunc, int inv);
void scatterComplex(double *fpdata, int pcStart, int trunc, int nsp);
void gatherComplex(double *fpdata, int pcStart, int trunc, int nsp);
void scm0(double *pdl, double *pdr, double *pfl, double *pfr, int klg);
int rowina2(double *p, int ko, int ki, double *pw,
int kcode, double msval, int *kret);
......@@ -330,8 +332,9 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
(lGrib[z++] = (Value) >> 8), \
(lGrib[z++] = (Value)))
#define Put1Int(Value) {ival = Value; if ( ival < 0 ) ival = 128 - ival; Put1Byte(ival);}
#define Put3Int(Value) {ival = Value; if ( ival < 0 ) ival = 8388608 - ival; Put3Byte(ival);}
#define Put1Int(Value) {ival = Value; if ( ival < 0 ) ival = 0x80 - ival; Put1Byte(ival);}
#define Put2Int(Value) {ival = Value; if ( ival < 0 ) ival = 0x8000 - ival; Put2Byte(ival);}
#define Put3Int(Value) {ival = Value; if ( ival < 0 ) ival = 0x800000 - ival; Put3Byte(ival);}
#define Put1Real(Value) \
{ \
......@@ -3196,11 +3199,11 @@ void encodeGDS(GRIBPACK *lGrib, long *gribLen, int *isec2, double *fsec2)
lonIncr = ISEC2_LonIncr;
latIncr = ISEC2_LatIncr;
}
Put2Byte(lonIncr); /* 23-24 i - direction increment */
Put2Byte(lonIncr); /* 23-24 i - direction increment */
if ( ISEC2_GridType == GTYPE_GAUSSIAN )
Put2Byte(ISEC2_NumPar); /* 25-26 Latitudes Pole->Equator */
else
Put2Byte(latIncr); /* 25-26 j - direction increment */
Put2Byte(latIncr); /* 25-26 j - direction increment */
Put1Byte(ISEC2_ScanFlag); /* 27 Scanning mode */
PutnZero(4); /* 28-31 reserved */
......@@ -3493,16 +3496,20 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
/* Uwe Schulzweida, 6/05/2003 : Copy result to fpval to prevent integer overflow */
static char func[] = "encodeBDS";
static int lwarn_cplx = TRUE;
long z = *gribLen;
long i, jloop;
int blockLength, PackStart, Flag;
int ival;
int blockLength, PackStart = 0, Flag = 0;
int binscale = 0;
int nbpv;
int bds_head = 11;
int bds_ext = 0;
/* ibits = BitsPerInt; */
unsigned int max_nbpv_pow2;
int exponent, mantissa;
int unused_bits = 0;
int lspherc = FALSE, lcomplex = FALSE;
int isubset = 0, itemp = 0, itrunc = 0;
double factor = 1, fmin, fmax, zref;
double range, rangec;
double jpepsln = 1.0e-12; /* -----> tolerance used to check equality */
......@@ -3510,10 +3517,57 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
/* on some platforms (eg vpp700, linux) */
extern double _pow2tab[158];
if ( isec4[3] == 64 && lwarn_cplx )
if ( isec2 )
{
lwarn_cplx = FALSE;
Message(func, "complex packing of spectral data unsupported, using simple packing!");
/* If section 2 is present, it says if data is spherical harmonic */
if ( isec2[0] == 50 || isec2[0] == 60 ||
isec2[0] == 70 || isec2[0] == 80 ) lspherc = TRUE;
if ( lspherc )
isec4[2] = 128;
else
isec4[2] = 0;
}
else
{
/* Section 4 says if it's spherical harmonic data.. */
lspherc = ( isec4[2] == 128 );
}
/* Complex packing supported for spherical harmonics. */
lcomplex = ( lspherc && ( isec4[3] == 64 ) ) ||
( lspherc && isec2 && ( isec2[5] == 2 ) );
/* Check input specification is consistent */
if ( lcomplex && isec2 )
{
if ( ( isec4[3] != 64 ) && ( isec2[5] == 2 ) )
{
gprintf(func, " COMPLEX mismatch. isec4[3] = %d\n", isec4[3]);
gprintf(func, " COMPLEX mismatch. isec2[5] = %d\n", isec2[5]);
return (807);
}
else if ( ( isec4[3] == 64 ) && ( isec2[5] != 2 ) )
{
gprintf(func, " COMPLEX mismatch. isec4[3] = %d\n", isec4[3]);
gprintf(func, " COMPLEX mismatch. isec2[5] = %d\n", isec2[5]);
return (807);
}
else if ( lcomplex )
{
/*
Truncation of full spectrum, which is supposed triangular,
has to be diagnosed. Define also sub-set truncation.
*/
isubset = isec4[17];
/* When encoding, use the total number of data. */
itemp = isec4[0];
itrunc = (int) (sqrt(itemp*4 + 1.) - 3) / 2;
}
}
if ( decscale )
......@@ -3522,30 +3576,43 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
for ( i = 0; i < datasize; ++i ) data[i] *= scale;
}
if ( ISEC2_GridType == GTYPE_SPECTRAL )
if ( lspherc )
{
*datstart = 15;
PackStart = 1;
Flag = 128;
blockLength = 11 + 4 + (ISEC4_NumBits*(datasize - 1) + 7)/8;
}
else
{
*datstart = 11;
PackStart = 0;
Flag = 0;
blockLength = 11 + (ISEC4_NumBits*datasize + 7)/8;
if ( lcomplex )
{
int jup, ioff;
jup = isubset;
ioff = (jup+1)*(jup+2);
bds_ext = 4 + 3 + 4*ioff;
PackStart = ioff;
Flag = 192;
}
else
{
bds_ext = 4;
PackStart = 1;
Flag = 128;
}
}
*datstart = bds_head + bds_ext;
blockLength = (*datstart) + (ISEC4_NumBits*(datasize - PackStart) + 7)/8;
if ( (blockLength%2) == 1 ) blockLength++;
if ( ISEC2_GridType == GTYPE_SPECTRAL )
unused_bits = blockLength*8 - 15*8 - ISEC4_NumBits*(datasize - 1);
else
unused_bits = blockLength*8 - 11*8 - ISEC4_NumBits*(datasize);
unused_bits = blockLength*8 - (*datstart)*8 - ISEC4_NumBits*(datasize - PackStart);
Flag += unused_bits;
if ( lspherc && lcomplex )
{
int pcStart, pcScale;
pcStart = isubset;
pcScale = isec4[16];
scaleComplex(data, pcStart, pcScale, itrunc, 0);
gatherComplex(data, pcStart, itrunc, datasize);
}
fmin = fmax = data[PackStart];
#if defined (CRAY)
......@@ -3664,7 +3731,25 @@ int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *ise
Put1Real(zref); /* 6-9 Reference value */
Put1Byte(nbpv); /* 10 Packing size */
if ( PackStart ) Put1Real(data[0]);
if ( lspherc )
{
if ( lcomplex )
{
int jup = isubset;
int ioff = z + bds_ext;
if ( ioff > 0xFFFF ) ioff = 0;
Put2Byte(ioff);
Put2Int(isec4[16]);
Put1Byte(jup);
Put1Byte(jup);
Put1Byte(jup);
for ( i = 0; i < ((jup+1)*(jup+2)); i++ ) Put1Real(data[i]);
}
else
{
Put1Real(data[0]);
}
}
*datsize = ((datasize-PackStart)*ISEC4_NumBits + 7)/8;
......@@ -3729,7 +3814,20 @@ void gribEncode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
encodePDS(lpds, pdsLen, isec1);
gribLen += pdsLen;
/*
if ( ( isec4[3] == 64 ) && ( isec2[5] == 2 ) )
{
static int lwarn_cplx = TRUE;
if ( lwarn_cplx )
Message(func, "Complex packing of spectral data unsupported, using simple packing!");
isec2[5] = 1;
isec4[3] = 0;
lwarn_cplx = FALSE;
}
*/
encodeGDS(lGrib, &gribLen, isec2, fsec2);
/*
----------------------------------------------------------------
......@@ -4418,17 +4516,18 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
{
static char func[] = "decodeBDS";
unsigned char *igrib;
int cplx, jup, kup, mup;
int lspherc = FALSE, lcomplex = FALSE;
int lcompress;
int jup, kup, mup;
int locnd;
int jlend, jlenc;
long i;
int bds_flag, bds_rep, jscale, imiss;
int bds_flag, jscale, imiss;
int bds_ubits;
int ioff = 0;
int iexp, imant;
int pcStart, pcScale;
int zoff;
int bds_aflag;
int bds_head = 11;
double fmin = 0., zscale = 0.;
double *fpdata = fsec4;
int bdsLen;
......@@ -4455,38 +4554,35 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
bds_flag = BDS_Flag;
bds_aflag = (bds_flag >> 4)&1; /* compress */
/* 0------- grid point */
/* 1------- spherical harmonics */
bds_rep = bds_flag >> 7;
lspherc = bds_flag >> 7;
if ( bds_rep == 0 ) isec4[2] = 0;
else isec4[2] = 128;
if ( lspherc ) isec4[2] = 128;
else isec4[2] = 0;
/* -0------ simple packing */
/* -1------ complex packing */
cplx = (bds_flag >> 6)&1;
lcomplex = (bds_flag >> 6)&1;
if ( cplx > 0 ) isec4[3] = 64;
if ( lcomplex ) isec4[3] = 64;
else isec4[3] = 0;
/* ---0---- No additional flags */
/* ---1---- No additional flags */
if ( (bds_flag >> 4)&1 )
lcompress = (bds_flag >> 4)&1; /* compress */
if ( lcompress )
{ isec4[5] = 16; isec4[6] = BDS_Z; zoff = 12; }
else
{ isec4[5] = 0; isec4[6] = 0; zoff = 0; }
if ( isec4[5] ) isec4[6] = BDS_Z; /* compression type: sz = 128 */
else isec4[6] = 0;
{ isec4[5] = 0; isec4[6] = 0; zoff = 0; }
/* ----++++ number of unused bits at end of section) */
bds_ubits = bds_flag & 15;
bds_ubits = bds_flag & 0xF;
/* scale factor (2 bytes) */;
......@@ -4501,15 +4597,16 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
/* convert reference value and scale factor. */
if ( imiss == 0 )
{
fmin = BDS_RefValue;
if ( jscale < 0 )
zscale = 1.0/intpow2(-jscale);
else
zscale = intpow2(jscale);
}
if ( ! (dfunc == 'J') )
if ( imiss == 0 )
{
fmin = BDS_RefValue;
if ( jscale < 0 )
zscale = 1.0/intpow2(-jscale);
else
zscale = intpow2(jscale);
}
/* get number of bits in each data value. */
......@@ -4518,14 +4615,14 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
/* octet number of start of packed data */
/* calculated from start of block 4 - 1 */
locnd = zoff + 11;
locnd = zoff + bds_head;
/* if data is in spherical harmonic form, distinguish */
/* between simple/complex packing (cplx = 0/1) */
/* between simple/complex packing (lcomplex = 0/1) */
if ( bds_rep == 1 )
if ( lspherc )
{
if ( cplx == 0 )
if ( !lcomplex )
{
/* no unpacked binary data present */
......@@ -4534,7 +4631,8 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
/* octet number of start of packed data */
/* calculated from start of block 4 - 1 */
locnd = zoff + 15;
ioff = 1;
locnd += 4*ioff; /* RealCoef */
/* get real (0,0) coefficient in grib format and */
/* convert to floating point. */
......@@ -4544,8 +4642,6 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
if ( imiss ) *fpdata++ = 0.0;
else *fpdata++ = BDS_RealCoef;
}
ioff = bds_rep;
}
else /* complex packed spherical harmonics */
{
......@@ -4565,10 +4661,12 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
/* unpacked binary data */
locnd = zoff + 18;
locnd += 4; /* 2 + power */
locnd += 3; /* j, k, m */
ioff = (jup+1)*(jup+2);
if ( dfunc != 'J' )
for (i = 0; i < ((jup+1)*(jup+2)); i++)
for ( i = 0; i < ioff; i++ )
{
iexp = (bds[locnd+4*i ]);
imant =((bds[locnd+4*i+1]) << 16) +
......@@ -4579,15 +4677,14 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
else *fpdata++ = decfp2(iexp,imant);
}
locnd = zoff + 18 + 4*(jup+1)*(jup+2);
ioff = (jup+1)*(jup+2);
locnd += 4*ioff; /* RealCoef */
}
}
/* Decode data values to floating point and store in fsec4. */
/* First calculate the number of data values. */
/* Take into account that spherical harmonics can be packed */
/* simple (cplx = 0) or complex (cplx = 1) */
/* simple (lcomplex = 0) or complex (lcomplex = 1) */
jlend = bdsLen - locnd;
if ( ISEC4_NumBits == 0 )
......@@ -4601,11 +4698,10 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
jlend = (jlend*8 - bds_ubits) / ISEC4_NumBits;
}
ISEC4_NumValues = jlend + ioff;
ISEC4_NumNonMissValues = 0;
if ( bds_aflag )
if ( lcompress )
{
size_t len;
......@@ -4616,22 +4712,21 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
ISEC4_NumValues = len*8/ISEC4_NumBits;
if ( bds_rep == 1 ) ISEC4_NumValues++;
if ( lspherc ) ISEC4_NumValues++;
}
if ( dfunc == 'J' ) return (bdsLen);
/* check length of output array. */
if ( ! (dfunc == 'J') )
if ( jlend+ioff > fsec4len )
{
*iret = 710;
gprintf(func, " Output array too small. Length = %d", fsec4len);
gprintf(func, " Number of values = %d", jlend+ioff);
gprintf(func, " Return code = %d", *iret);
return (0);
}
if ( jlend+ioff > fsec4len )
{
*iret = 710;
gprintf(func, " Output array too small. Length = %d", fsec4len);
gprintf(func, " Number of values = %d", jlend+ioff);
gprintf(func, " Return code = %d", *iret);
return (0);
}
if ( imiss ) memset((char *)fpdata, 0, jlend*sizeof(double));
else
......@@ -4642,12 +4737,13 @@ int decodeBDS(int decscale, unsigned char *bds, int *isec2, int *isec4,
decode_double_array(igrib, jlend, jlenc, ISEC4_NumBits, fmin, zscale, fpdata);
}
if ( bds_rep == 1 && cplx == 1 )
if ( lspherc && lcomplex )
{
int pcStart, pcScale;
pcStart = isec4[19];
pcScale = isec4[16];
scatterComplex(fsec4, pcStart, ISEC2_PentaJ, ISEC4_NumValues);
scaleComplex(fsec4, pcStart, pcScale, ISEC2_PentaJ);
scaleComplex(fsec4, pcStart, pcScale, ISEC2_PentaJ, 1);
}
......@@ -5489,15 +5585,121 @@ int gribrec_len(int b1, int b2, int b3)
/* ============ */
/* scaleComplex */
/* ============ */
/* calculate_pfactor: source code from grib_api-1.8.0 */
double calculate_pfactor(const double* spectralField, long fieldTruncation, long subsetTruncation)
{
static char func[] = "calculate_pfactor";
/*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
long loop, index, m, n = 0;
double pFactor, zeps = 1.0e-15;
long ismin = (subsetTruncation+1), ismax = (fieldTruncation+1);
double* weights, range, * norms;
double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0, x, y;
double numerator = 0.0, denominator = 0.0, slope;
/*
// Setup the weights
*/
range = (double) (ismax - ismin +1);
weights = (double*) malloc((ismax+1)*sizeof(double));
for( loop = ismin; loop <= ismax; loop++ )
weights[loop] = range / (double) (loop-ismin+1);
/*
// Compute norms
// Handle values 2 at a time (real and imaginary parts).
*/
norms = (double*) malloc((ismax+1)*sizeof(double));
for( loop = 0; loop < ismax+1; loop++ ) norms[loop] = 0.0;
/*
// Form norms for the rows which contain part of the unscaled subset.
*/
index = -2;
for( m = 0; m < subsetTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
index += 2;
if( n >= subsetTruncation ) {
double tval = spectralField[index];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
}
/*
// Form norms for the rows which do not contain part of the unscaled subset.
*/
for( m = subsetTruncation; m <= fieldTruncation; m++ )
for( n = m; n <= fieldTruncation; n++ ) {
double tval = spectralField[index];
index += 2;
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
tval = spectralField[index+1];
tval=tval<0?-tval:tval;
norms[n] = norms[n] > tval ? norms[n] : tval;
}
/*
// Ensure the norms have a value which is not too small in case of
// problems with math functions (e.g. LOG).
*/
void scaleComplex(double *fpdata, int pcStart, int pcScale, int truncation)
for( loop = ismin; loop <= ismax; loop++ ) {
norms[n] = norms[n] > zeps ? norms[n] : zeps;
if( norms[n] == zeps ) weights[n] = 100.0 * zeps;
}
/*
// Do linear fit to find the slope
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double) (loop*(loop+1)) );
y = log( norms[loop] );
weightedSumOverX = weightedSumOverX + x * weights[loop];
weightedSumOverY = weightedSumOverY + y * weights[loop];
sumOfWeights = sumOfWeights + weights[loop];
}
weightedSumOverX = weightedSumOverX / sumOfWeights;
weightedSumOverY = weightedSumOverY / sumOfWeights;
/*
// Perform a least square fit for the equation
*/
for( loop = ismin; loop <= ismax; loop++ ) {
x = log( (double)(loop*(loop+1)) );
y = log( norms[loop] );
numerator =
numerator + weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
denominator =
denominator + weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
}
slope = numerator / denominator;
free(weights);
free(norms);
pFactor = -slope;
if( pFactor < -9999.9 ) pFactor = -9999.9;
if( pFactor > 9999.9 ) pFactor = 9999.9;
return pFactor;
}
void scaleComplex(double *fpdata, int pcStart, int pcScale, int trunc, int inv)
{
static char func[] = "scaleComplex";
double power;
double *scale = (double *) malloc((truncation+1)*sizeof(double));
double *scale = (double *) malloc((trunc+1)*sizeof(double));
int n, m;
int index;
......@@ -5516,20 +5718,23 @@ void scaleComplex(double *fpdata, int pcStart, int pcScale, int truncation)
power = (double) pcScale / 1000.;
scale[0] = 1.0;
for ( n = 1; n <= truncation; n++ )
for ( n = 1; n <= trunc; n++ )
{
if (pcScale != 1000)
scale[n] = 1.0 / pow((double) (n*(n+1)), power);
scale[n] = pow((double) (n*(n+1)), power);
else
scale[n] = 1.0 / (double) (n*(n+1));
scale[n] = (double) (n*(n+1));
}
if ( inv )
for ( n = 1; n <= trunc; n++ ) scale[n] = 1.0 / scale[n];
/* Scale the values */
index = 0;