Skip to content
GitLab
Commit 104cf06d authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

No commit message 

No commit message
parent 42e313da
Hide whitespace changes
Inline Side-by-side
src/cgribex.h
View file @ 104cf06d
#ifndef _GRIB_H
#define _GRIB_H
#ifndef _CGRIBEX_H
#define _CGRIBEX_H
#include <stdio.h>
#include <sys/types.h>
......@@ -72,7 +72,7 @@
#define ISEC1_TimePeriod2 (isec1[16]) /* P2 Time period */
#define ISEC1_TimeRange (isec1[17]) /* Time range indicator */
#define ISEC1_AvgNum (isec1[18]) /* Number of products included in an average */
#define ISEC1_AvgMiss (isec1[19]) /* Number of products missing form an average */
#define ISEC1_AvgMiss (isec1[19]) /* Number of products missing from an average */
#define ISEC1_Century (isec1[20]) /* Century */
#define ISEC1_SubCenterID (isec1[21]) /* Subcenter identifier */
#define ISEC1_DecScaleFactor (isec1[22]) /* Decimal scale factor */
......@@ -179,7 +179,7 @@ void gribExDP(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
int kleng, int *kword, char *hoper, int *kret);
const char *gribLibraryVersion(void);
const char *cgribexLibraryVersion(void);
void gribDebug(int debug);
......@@ -205,8 +205,11 @@ void gribPrintGDS(int nrec, long recpos, long recsize, unsigned char *gribbuffe
void gribPrintBMS(int nrec, long recpos, long recsize, unsigned char *gribbuffer);
void gribPrintBDS(int nrec, long recpos, long recsize, unsigned char *gribbuffer);
int gribSections(unsigned char *gribbuffer, long recsize, unsigned char **pdsp,
unsigned char **gdsp, unsigned char **bmsp, unsigned char **bdsp);
int grib1Sections(unsigned char *gribbuffer, long recsize, unsigned char **pdsp,
unsigned char **gdsp, unsigned char **bmsp, unsigned char **bdsp);
int grib2Sections(unsigned char *gribbuffer, long recsize, unsigned char **idsp,
unsigned char **lusp, unsigned char **gdsp, unsigned char **pdsp,
unsigned char **drsp, unsigned char **bmsp, unsigned char **bdsp);
int gribGetZip(long recsize, unsigned char *gribbuffer, long *urecsize);
......@@ -226,6 +229,7 @@ int gribFileSeek(int fileID, long *offset);
int gribReadSize(int fileID);
int gribVersion(unsigned char *buffer, size_t buffersize);
int gribGinfo(long recpos, long recsize, unsigned char *gribbuffer, int *intnum, float *fltnum);
int gribGinfo(long recpos, long recsize, unsigned char *gribbuffer, int *intnum, float *fltnum);
#endif /* _CGRIBEX_H */
#endif /* _GRIB_H */
src/cgribexlib.c
View file @ 104cf06d
/* Automatically generated by m214003 at 2009-07-03, do not edit */
/* Automatically generated by m214003 at 2009-09-03, do not edit */
/* GRIBLIB_VERSION="1.4.0" */
/* CGRIBEXLIB_VERSION="1.4.0" */
#if defined (HAVE_CONFIG_H)
# include "config.h"
......@@ -30,8 +30,8 @@
#include <float.h>
#if ! defined (_GRIB_H)
# include "grib.h"
#if ! defined (_CGRIBEX_H)
# include "cgribex.h"
#endif
#if ! defined (_ERROR_H)
# include "error.h"
......@@ -52,6 +52,12 @@
# define UCHAR unsigned char
#endif
#if defined (CRAY) || defined (SX) || defined (__uxpch__) || defined (__xlC__)
# define VECTORCODE
#endif
#if defined (VECTORCODE)
#if defined (INT32)
# define GRIBPACK unsigned INT32
# define PACK_GRIB packInt32
......@@ -61,6 +67,9 @@
# define PACK_GRIB packInt64
# define UNPACK_GRIB unpackInt64
#endif
#else
# define GRIBPACK unsigned char
#endif
#ifndef DBL_IS_NAN
......@@ -115,9 +124,9 @@ int qu2reg3(double *pfield, int *kpoint, int klat, int klon,
double msval, int *kret, int omisng, int operio, int oveggy);
#if defined (INT32)
long packInt32(unsigned INT32 *up, char *cp, long bc, long tc);
long packInt32(unsigned INT32 *up, unsigned char *cp, long bc, long tc);
#endif
long packInt64(unsigned INT64 *up, char *cp, long bc, long tc);
long packInt64(unsigned INT64 *up, unsigned char *cp, long bc, long tc);
#if defined (INT32)
long unpackInt32(unsigned char *cp, unsigned INT32 *up, long bc, long tc);
#endif
......@@ -151,7 +160,7 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
#define BUDG_START(s) (s[0]=='B' && s[1]=='U' && s[2]=='D' && s[3]=='G')
#define TIDE_START(s) (s[0]=='T' && s[1]=='I' && s[2]=='D' && s[3]=='E')
#define GRIB_START(s) (s[0]=='G' && s[1]=='R' && s[2]=='I' && s[3]=='B')
#define GRIB_END(s) (s[0]=='7' && s[1]=='7' && s[2]=='7' && s[3]=='7')
#define GRIB_FIN(s) (s[0]=='7' && s[1]=='7' && s[2]=='7' && s[3]=='7')
/* GRIB1 Section 0: Indicator Section (IS) */
......@@ -1140,7 +1149,7 @@ gribExDP(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
fsec3, isec4, fsec4, klenp, kgrib,
kleng, kword, yfunc, kret);
else if ( yfunc == 'V' )
fprintf(stderr, " c-gribex: Version is %s\n", gribLibraryVersion());
fprintf(stderr, " c_gribex: Version is %s\n", cgribexLibraryVersion());
else
{
Error(func, "oper %c unsupported\n", yfunc);
......@@ -1215,7 +1224,7 @@ gribExSP(int *isec0, int *isec1, int *isec2, float *fsec2sp, int *isec3,
free(fsec4dp);
}
else if ( yfunc == 'V' )
fprintf(stderr, " c-gribex: Version is %s\n", gribLibraryVersion());
fprintf(stderr, " c-gribex: Version is %s\n", cgribexLibraryVersion());
else
{
Error(func, "oper %c unsupported\n", yfunc);
......@@ -2546,7 +2555,6 @@ void gribPrintSec4DP(int *isec0, int *isec4, double *fsec4)
Uwe Schulzweida MPIfM 01/04/2001
*/
int inum;
int j;
......@@ -2742,10 +2750,10 @@ int BitsPerInt = (int) (sizeof(int) * 8);
/* GRIB block 0 - indicator block */
void encodeIS(GRIBPACK *lGrib, int *gribLen)
static
void encodeIS(GRIBPACK *lGrib, long *gribLen)
{
int z = *gribLen;
long z = *gribLen;
lGrib[0] = 'G';
lGrib[1] = 'R';
......@@ -2766,10 +2774,10 @@ void encodeIS(GRIBPACK *lGrib, int *gribLen)
/* GRIB block 5 - end block */
void encodeES(GRIBPACK *lGrib, int *gribLen)
static
void encodeES(GRIBPACK *lGrib, long *gribLen)
{
int z = *gribLen;
long z = *gribLen;
lGrib[z++] = '7';
lGrib[z++] = '7';
......@@ -2797,9 +2805,10 @@ void encodeES(GRIBPACK *lGrib, int *gribLen)
#define DWD_extension_254_len 26
#define ECMWF_extension_1_len 24
int getLocalExtLen(int *isec1)
static
long getLocalExtLen(int *isec1)
{
int extlen = 0;
long extlen = 0;
if ( ISEC1_LocalFLag )
{
......@@ -2826,21 +2835,22 @@ int getLocalExtLen(int *isec1)
return (extlen);
}
int getPdsLen(int *isec1)
static
long getPdsLen(int *isec1)
{
int pdslen = 28;
long pdslen = 28;
pdslen += getLocalExtLen(isec1);
return (pdslen);
}
void encodePDS_DWD_local_Extension_254(GRIBPACK *lGrib, int *zs, int *isec1)
static
void encodePDS_DWD_local_Extension_254(GRIBPACK *lGrib, long *zs, int *isec1)
{
int i, localextlen, isvn;
int z = *zs;
int isvn;
long localextlen, i;
long z = *zs;
localextlen = getLocalExtLen(isec1);
for ( i = 0; i < localextlen-2; i++ )
......@@ -2854,11 +2864,12 @@ void encodePDS_DWD_local_Extension_254(GRIBPACK *lGrib, int *zs, int *isec1)
*zs = z;
}
void encodePDS_DWD_local_Extension_253(GRIBPACK *lGrib, int *zs, int *isec1)
static
void encodePDS_DWD_local_Extension_253(GRIBPACK *lGrib, long *zs, int *isec1)
{
int i, localextlen, isvn;
int z = *zs;
int isvn;
long localextlen, i;
long z = *zs;
localextlen = DWD_extension_254_len;
for ( i = 0; i < localextlen-2; i++ )
......@@ -2879,12 +2890,13 @@ void encodePDS_DWD_local_Extension_253(GRIBPACK *lGrib, int *zs, int *isec1)
*zs = z;
}
void encodePDS_ECMWF_local_Extension_1 (GRIBPACK *lGrib, int *zs, int *isec1)
static
void encodePDS_ECMWF_local_Extension_1(GRIBPACK *lGrib, long *zs, int *isec1)
{
int i, localextlen, isvn;
int z = *zs;
static char func[] = "encodePDS_ECMWF_local_Extension_1";
// static char func[] = "encodePDS_ECMWF_local_Extension_1";
// int isvn;
long localextlen, i;
long z = *zs;
localextlen = getLocalExtLen(isec1);
for ( i = 0; i < localextlen-12; i++ )
......@@ -2913,11 +2925,11 @@ void encodePDS_ECMWF_local_Extension_1 (GRIBPACK *lGrib, int *zs, int *isec1)
/* GRIB BLOCK 1 - PRODUCT DESCRIPTION SECTION */
void encodePDS(GRIBPACK *lpds, int pdsLen, int *isec1)
static
void encodePDS(GRIBPACK *lpds, long pdsLen, int *isec1)
{
GRIBPACK *lGrib = lpds;
int z = 0;
long z = 0;
int ival, century, year;
century = ISEC1_Century;
......@@ -3028,7 +3040,7 @@ void encodePDS(GRIBPACK *lpds, int pdsLen, int *isec1)
}
else
{
int i, localextlen;
long i, localextlen;
localextlen = getLocalExtLen(isec1);
for ( i = 0; i < localextlen; i++ )
{
......@@ -3040,13 +3052,13 @@ void encodePDS(GRIBPACK *lpds, int pdsLen, int *isec1)
/* GRIB BLOCK 2 - GRID DESCRIPTION SECTION */
void encodeGDS(GRIBPACK *lGrib, int *gribLen, int *isec2, double *fsec2)
static
void encodeGDS(GRIBPACK *lGrib, long *gribLen, int *isec2, double *fsec2)
{
static char func[] = "encodeGDS";
int z = *gribLen;
long z = *gribLen;
int exponent, mantissa;
int i;
long i;
int ival;
int pvoffset = 255;
int gdslen = 32;
......@@ -3173,17 +3185,17 @@ void encodeGDS(GRIBPACK *lGrib, int *gribLen, int *isec2, double *fsec2)
/* GRIB BLOCK 3 - BIT MAP SECTION */
void encodeBMS(GRIBPACK *lGrib, int *gribLen, double *fsec3, int *isec4, double *data, int *datasize)
static
void encodeBMS(GRIBPACK *lGrib, long *gribLen, double *fsec3, int *isec4, double *data, long *datasize)
{
static char func[] = "encodeBMS";
GRIBPACK *bitmap;
int bitmapSize;
int imaskSize;
int i;
int bmsLen, bmsUnusedBits;
int fsec4size;
int z = *gribLen;
long bitmapSize;
long imaskSize;
long i;
long bmsLen, bmsUnusedBits;
long fsec4size;
long z = *gribLen;
unsigned int *imask;
/* unsigned int c, imask; */
......@@ -3256,28 +3268,185 @@ void encodeBMS(GRIBPACK *lGrib, int *gribLen, double *fsec3, int *isec4, double
*datasize = fsec4size;
}
static
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;
unsigned int ipval;
int cbits, jbits;
unsigned int c;
static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
/* code from gribw routine flist2bitstream */
/* GRIB BLOCK 4 - BINARY DATA SECTION */
cbits = 8;
c = 0;
for ( i = PackStart; i < datasize; i++ )
{
/* note float -> unsigned int .. truncate */
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
/*
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
jbits = NumBits;
while ( cbits <= jbits )
{
if ( cbits == 8 )
{
jbits -= 8;
lGrib[z++] = (ipval >> jbits) & 255;
}
else
{
jbits -= cbits;
lGrib[z++] = (c << cbits) + ((ipval >> jbits) & mask[cbits]);
cbits = 8;
c = 0;
}
}
/* now jbits < cbits */
if ( jbits )
{
c = (c << jbits) + (ipval & mask[jbits]);
cbits -= jbits;
}
}
if (cbits != 8) lGrib[z++] = c << cbits;
*gz = z;
}
static
void encode_double_array(int NumBits, long PackStart, long datasize, GRIBPACK *lGrib,
const double *data, double zref, double factor, long *gz)
{
static char func[] = "encode_double_array";
long i, z = *gz;
unsigned int ipval;
if ( NumBits == 8 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
/*
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
lGrib[z ] = ipval;
z++;
}
}
else if ( NumBits == 16 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
/*
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
lGrib[z ] = ipval >> 8;
lGrib[z+1] = ipval;
z += 2;
}
}
else if ( NumBits == 24 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
/*
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
lGrib[z ] = ipval >> 16;
lGrib[z+1] = ipval >> 8;
lGrib[z+2] = ipval;
z += 3;
}
}
else if ( NumBits == 32 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
ipval = (unsigned int) ((data[i] - zref) * factor + 0.5);
/*
if ( ipval > max_nbpv_pow2 ) ipval = max_nbpv_pow2;
if ( ipval < 0 ) ipval = 0;
*/
lGrib[z ] = ipval >> 24;
lGrib[z+1] = ipval >> 16;
lGrib[z+2] = ipval >> 8;
lGrib[z+3] = ipval;
z += 4;
}
}
else if ( NumBits > 0 && NumBits <= 32 )
{
encode_double_array_common(NumBits, PackStart, datasize, lGrib,
data, zref, factor, &z);
}
else if ( NumBits == 0 )
{
}
else
{
Error(func, "Unimplemented packing factor %d", NumBits);
}
*gz = z;
}
int encodeBDS(GRIBPACK *lGrib, int *gribLen, int decscale, int *isec2, int *isec4, int datasize, double data[],
int *datstart, int *datsize)
/* GRIB BLOCK 4 - BINARY DATA SECTION */
static
int encodeBDS(GRIBPACK *lGrib, long *gribLen, int decscale, int *isec2, int *isec4, long datasize, double data[],
long *datstart, long *datsize)
{
/* Uwe Schulzweida, 11/04/2003 : Check that number of bits per value is not exceeded */
/* Uwe Schulzweida, 6/05/2003 : Copy result to fpval to prevent integer overflow */
static char func[] = "encodeBDS";
static int lwarn_cplx = TRUE;
int z = *gribLen;
int i, jloop;
long z = *gribLen;
long i, jloop;
int blockLength, PackStart, Flag;
int binscale = 0;
unsigned int pval;
int nbpv;
/* ibits = BitsPerInt; */
double max_nbpv_pow2;
/* unsigned int max_nbpv_pow2; */
int exponent, mantissa;
int unused_bits = 0;
double fpval;
double factor = 1, fmin, fmax, zref;
double range, rangec;
double jpepsln = 1.0e-12; /* -----> tolerance used to check equality */
......@@ -3334,6 +3503,10 @@ int encodeBDS(GRIBPACK *lGrib, int *gribLen, int decscale, int *isec2, int *isec
{
if ( fmin > data[i] ) fmin = data[i];
if ( fmax < data[i] ) fmax = data[i];
/*
fmin = fmin < data[i] ? fmin : data[i];
fmax = fmax > data[i] ? fmax : data[i];
*/
}
zref = fmin;
......@@ -3413,7 +3586,7 @@ int encodeBDS(GRIBPACK *lGrib, int *gribLen, int decscale, int *isec2, int *isec
else factor = 1.0/intpow2( binscale);
}
max_nbpv_pow2 = (double) (intpow2(nbpv) - 1);
/* max_nbpv_pow2 = (unsigned int) (intpow2(nbpv) - 1); */
ref2ibm(&zref, BitsPerInt);
......@@ -3428,138 +3601,8 @@ int encodeBDS(GRIBPACK *lGrib, int *gribLen, int decscale, int *isec2, int *isec
*datsize = ((datasize-PackStart)*ISEC4_NumBits + 7)/8;
if ( ISEC4_NumBits == 8 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
fpval = ((data[i] - zref) * factor + 0.5);
if ( fpval > max_nbpv_pow2 ) fpval = max_nbpv_pow2;
if ( fpval < 0 ) fpval = 0;
pval = (unsigned int) fpval;
lGrib[z ] = pval;
z++;
}
}
else if ( ISEC4_NumBits == 16 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
fpval = ((data[i] - zref) * factor + 0.5);
if ( fpval > max_nbpv_pow2 ) fpval = max_nbpv_pow2;
if ( fpval < 0 ) fpval = 0;
pval = (unsigned int) fpval;
lGrib[z ] = pval >> 8;
lGrib[z+1] = pval;
z += 2;
}
}
else if ( ISEC4_NumBits == 24 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
fpval = ((data[i] - zref) * factor + 0.5);
if ( fpval > max_nbpv_pow2 ) fpval = max_nbpv_pow2;
if ( fpval < 0 ) fpval = 0;
pval = (unsigned int) fpval;
lGrib[z ] = pval >> 16;
lGrib[z+1] = pval >> 8;
lGrib[z+2] = pval;
z += 3;
}
}
else if ( ISEC4_NumBits == 32 )
{
#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
for ( i = PackStart; i < datasize; i++ )
{
fpval = ((data[i] - zref) * factor + 0.5);
if ( fpval > max_nbpv_pow2 ) fpval = max_nbpv_pow2;
if ( fpval < 0 ) fpval = 0;
pval = (unsigned int) fpval;
lGrib[z ] = pval >> 24;
lGrib[z+1] = pval >> 16;
lGrib[z+2] = pval >> 8;
lGrib[z+3] = pval;
z += 4;
}
}
else if ( ISEC4_NumBits > 0 && ISEC4_NumBits <= 32 )
{
int cbits, jbits;
unsigned int c;
static 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++ )
{
/* note float -> unsigned int .. truncate */
fpval = ((data[i] - zref) * factor + 0.5);
if ( fpval > max_nbpv_pow2 ) fpval = max_nbpv_pow2;
if ( fpval < 0 ) fpval = 0;
pval = (unsigned int) fpval;
jbits = ISEC4_NumBits;
while ( cbits <= jbits )
{
if ( cbits == 8 )
{