mirror of https://github.com/ecmwf/eccodes.git
823 lines
21 KiB
C
823 lines
21 KiB
C
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <unistd.h>
|
|
#include <string.h>
|
|
#include "grib_api.h"
|
|
#include <sys/param.h>
|
|
|
|
#include <dirent.h>
|
|
#include <sys/types.h>
|
|
#include <time.h>
|
|
#include <sys/time.h>
|
|
#include <sys/stat.h>
|
|
#include <sys/resource.h>
|
|
|
|
#include <math.h>
|
|
#include "test.h"
|
|
|
|
#if 0
|
|
double maxRelativeError = 1e-4; /*1e-7;*/
|
|
double maxAbsoluteError = 1e-50;
|
|
#else
|
|
const double maxRelativeError = 1e-7; /*1e-7;*/
|
|
const double maxAbsoluteError = 0;
|
|
#endif
|
|
|
|
|
|
#define ISECTION_0 2
|
|
#define ISECTION_1 1024 /* beware of for ocean data */
|
|
#define ISECTION_2 3000
|
|
#define ISECTION_3 2
|
|
#define ISECTION_4 512
|
|
|
|
#define RSECTION_2 512
|
|
#define RSECTION_3 2
|
|
#define RSECTION_4 1
|
|
|
|
|
|
struct { int code; char *message; } gribex_errors[] = {
|
|
{201 ," An invalid coding/decoding option was requested.",},
|
|
{202 ,"The number of bits per data value exceeds a word length.",},
|
|
{203 ,"Missing data indicated and data field contains non-zero values.",},
|
|
{204 ,"Number of bits per data value is invalid.",},
|
|
{205 ,"Number of bits per data value is same as the number of bits per computer word.",},
|
|
{301 ,"Error in encoding or decoding the letters GRIB.",},
|
|
{302 ,"Error decoding the length of the GRIB message.",},
|
|
{303 ,"Error encoding or decoding the GRIB edition number.",},
|
|
{304 ,"Error decoding octets 22 and 23, experimental edition check.",},
|
|
{305 ,"Input data is not GRIB or pseudo-GRIB.",},
|
|
{401 ,"Error encoding or decoding the length of section 1.",},
|
|
{402 ,"Error encoding or decoding the parameter version number.",},
|
|
{403 ,"Error encoding or decoding six header fields in section 1, from 'identification of centre' to 'indicator of type of level'.",},
|
|
{404 ,"Error encoding or decoding height, pressure, etc of levels.",},
|
|
{405 ,"Error encoding or decoding six fields in section 1, from 'year of century' to 'indicator of unit of time range'.",},
|
|
{406 ,"Error encoding or decoding period of time.",},
|
|
{407 ,"Error encoding or decoding time range indicator.",},
|
|
{408 ,"Error encoding or decoding number averaged.",},
|
|
{409 ,"Error encoding or decoding number missing from averages etc.",},
|
|
{410 ,"Error encoding or decoding century of data or reserved field.",},
|
|
{411 ,"Error encoding or decoding units decimal scale factor.",},
|
|
{412 ,"Error encoding or decoding ECMWF local data.",},
|
|
{413 ,"Grib edition not handled.",},
|
|
{499 ,"Error found when checking values for section 1 against expected valid GRIB values.",},
|
|
{501 ,"Error encoding or decoding length of section 2.",},
|
|
{502 ,"Error encoding or decoding number of vertical coordinate parameters.",},
|
|
{503 ,"Error encoding or decoding location of list of vertical coordinate parameters or list of numbers of points.",},
|
|
{504 ,"Error encoding or decoding data representation type.",},
|
|
{505 ,"Error encoding or decoding number of points along a parallel or meridian.",},
|
|
{506 ,"Error encoding or decoding latitude or longitude of the first grid point.",},
|
|
{507 ,"Error encoding or decoding the components flag.",},
|
|
{508 ,"Error encoding or decoding latitude or longitude of the last grid point.",},
|
|
{509 ,"Error encoding or decoding the i direction increment.",},
|
|
{510 ,"Error encoding or decoding the number of parallels between the pole and the Equator.",},
|
|
{511 ,"Error encoding or decoding scanning mode flags.",},
|
|
{513 ,"Error encoding or decoding the j direction increment.",},
|
|
{514 ,"Error encoding or decoding J,K,M pentagonal resolution parameters.",},
|
|
{515 ,"Error encoding or decoding representation type or mode.",},
|
|
{517 ,"Error encoding or decoding latitude or longitude of southern pole.",},
|
|
{518 ,"Error encoding or decoding angle of rotation.",},
|
|
{519 ,"Error encoding or decoding latitude or longitude of pole of stretching.",},
|
|
{520 ,"Error encoding or decoding stretching factor.",},
|
|
{521 ,"Error encoding or decoding vertical coordinate parameters.",},
|
|
{522 ,"Error encoding or decoding list of numbers of points.",},
|
|
{523 ,"Error encoding or decoding number of points along X or Y axis.",},
|
|
{524 ,"Error encoding or decoding X or Y axis grid lengths.",},
|
|
{525 ,"Error encoding or decoding the projection centre flag.",},
|
|
{526 ,"Error encoding or decoding the latitude or longitude of sub-satellite point.",},
|
|
{527 ,"Error encoding or decoding the diameter of the Earth in X or Y direction.",},
|
|
{528 ,"Error encoding or decoding the X or Y coordinate of the sub-satellite point.",},
|
|
{529 ,"Error encoding or decoding the orientation of the grid or camera angle.",},
|
|
{530 ,"Error encoding or decoding the X or Y coordinate of the origin of sector.",},
|
|
{531 ,"Error inserting/extracting Latin1 or Latin2 of secants points.",},
|
|
{532 ,"Error inserting/extracting latitude/longitude of southern pole.",},
|
|
{540 ,"Error encoding gaussian grid section 2.",},
|
|
{541 ,"Error decoding gaussian grid section 2.",},
|
|
{542 ,"Error encoding lat/long grid section 2.",},
|
|
{543 ,"Error decoding lat/long grid section .",},
|
|
{544 ,"Error encoding spectral section 2.",},
|
|
{545 ,"Error decoding spectral section 2.",},
|
|
{546 ,"Error encoding space view section 2, or error encoding ECMWF ocean section 2.",},
|
|
{547 ,"Error decoding space view section 2, or error decoding ECMWF ocean section 2.",},
|
|
{598 ,"Representation type not handled.",},
|
|
{599 ,"Error found when checking values for section 2 against expected valid GRIB values.",},
|
|
{601 ,"Error encoding or decoding length of section 3.",},
|
|
{602 ,"Error encoding or decoding the number of unused bits at the end of section 3.",},
|
|
{603 ,"Error encoding or decoding the bitmap reference table.",},
|
|
{604 ,"Error encoding or decoding the primary bitmap.",},
|
|
{605 ,"Unable to convert quasi-regular (reduced) gaussian grid with a bitmap.",},
|
|
{699 ,"Error found when checking values for section 3 against expected valid GRIB values.",},
|
|
{701 ,"Error encoding or decoding the length of section 4.",},
|
|
{703 ,"Second-order packing implies additional flags.",},
|
|
{704 ,"Function A or B invalid for second-order packed field (grid-point).",},
|
|
{705 ,"Only simple packing handled for grid data.",},
|
|
{706 ,"Error in decoding the section 4 flag field.",},
|
|
{707 ,"Error encoding or decoding the field scale factor.",},
|
|
{708 ,"Error encoding or decoding the field reference value.",},
|
|
{709 ,"Error encoding or decoding the number of bits per data value.",},
|
|
{710 ,"Output array too small.",},
|
|
{711 ,"Error encoding or decoding a REAL coefficient.",},
|
|
{712 ,"Error encoding or decoding data values.",},
|
|
{713 ,"Error encoding or decoding the 'flag and unused bits' field.",},
|
|
{714 ,"Decoding option is 'X' and number of values is illegal.",},
|
|
{715 ,"Decoding option is 'X' and scanning mode is not north to south and west to east.",},
|
|
{716 ,"Decoding option is 'X' and field is not a gaussian or latitude/longitude grid.",},
|
|
{717 ,"Decoding option is 'X' and a bitmap is encountered.",},
|
|
{720 ,"Error encoding or decoding the octet number at which packed data begins.",},
|
|
{721 ,"Error encoding or decoding the 'extended flag' field.",},
|
|
{722 ,"Error encoding or decoding the first or second dimension of a matrix.",},
|
|
{723 ,"Error encoding or decoding six fields, from first dimension coordinate value onwards.",},
|
|
{724 ,"Error encoding or decoding first or second dimension coefficients.",},
|
|
{725 ,"Error encoding a secondary bitmap.",},
|
|
{726 ,"Output array is not big enough to allow expansion using bitmap(s), or constant (0-bit) field without section 2 not supported.",},
|
|
{727 ,"Invalid predefined bitmap number given, or constant (0-bit) spectral field not supported.",},
|
|
{728 ,"Open of predefined bitmap failed, or constant (0-bit) field function G or B not supported.",},
|
|
{729 ,"Problem reading a predefined bitmap size, or output array is not big enough to allow expansion using bitmaps.",},
|
|
{730 ,"Problem allocating memory for a predefined bitmap.",},
|
|
{731 ,"Problem reading the number of non-missing points in a predefined bitmap.",},
|
|
{732 ,"Problem reading a predefined bitmap.",},
|
|
{733 ,"Problem closing a predefined bitmap file.",},
|
|
{734 ,"Error inserting padding zeroes at end of section 4.",},
|
|
{796 ,"Error extracting pointer to packed data.",},
|
|
{797 ,"Error extracting sub-set truncation or constant (0-bit) field function G or B not supported.",},
|
|
{798 ,"Function is X but no section 2 included.",},
|
|
{799 ,"Error found when checking values for section 4 against expected valid GRIB values.",},
|
|
{800 ,"Error encoding or decoding the 7777 group.",},
|
|
{801 ,"Error inserting or extracting the 7777 group.",},
|
|
{802 ,"Error encoding or decoding the length of the GRIB message.",},
|
|
{805 ,"End of message 7777 group not found.",},
|
|
{806 ,"Error in decoding a primary or secondary bitmap.",},
|
|
{807 ,"Inconsistent values specified for complex packing in KSEC2(6) and KSEC4(4).",},
|
|
{808 ,"Error converting quasi-regular (reduced) gaussian grid to regular.",},
|
|
{809 ,"Error padding unused part of GRIB to zero.",},
|
|
{810 ,"Error inserting dummy zero, or error padding GRIB to multiple of 120.",},
|
|
{811 ,"Unable to handle secondary bitmaps for J option.",},
|
|
};
|
|
|
|
extern double grib_power(long,long);
|
|
|
|
|
|
double bucket_size(grib_handle* g)
|
|
{
|
|
long bits;long binary_scale_factor;long decimal_scale_factor;
|
|
unsigned long zero = 0;
|
|
unsigned long one = 1;
|
|
double z,o,s,d,reference_value;
|
|
int err=0;
|
|
|
|
err=grib_get_double(g,"referenceValue",&reference_value);
|
|
if (err!=GRIB_SUCCESS) return 0.0;
|
|
|
|
grib_get_long(g,"binaryScaleFactor",&binary_scale_factor);
|
|
if (err!=GRIB_SUCCESS) return 0.0;
|
|
|
|
grib_get_long(g,"numberOfBitsContainingEachPackedValue",&bits);
|
|
if (err!=GRIB_SUCCESS) return 0.0;
|
|
|
|
grib_get_long(g,"decimalScaleFactor",&decimal_scale_factor);
|
|
if (err!=GRIB_SUCCESS) return 0.0;
|
|
|
|
s = grib_power(binary_scale_factor,2);
|
|
d = grib_power(-decimal_scale_factor,10) ;
|
|
|
|
z = (double) (((zero*s)+reference_value)*d);
|
|
o = (double) (((one *s)+reference_value)*d);
|
|
|
|
return o-z;
|
|
}
|
|
|
|
|
|
static double err(double A, double B)
|
|
{
|
|
double relativeError;
|
|
|
|
if(fabs(A) <= maxAbsoluteError || fabs(B) <= maxAbsoluteError)
|
|
relativeError = fabs(A-B);
|
|
else if (fabs(B) > fabs(A))
|
|
relativeError = fabs((A - B) / B);
|
|
else
|
|
relativeError = fabs((A - B) / A);
|
|
|
|
return relativeError;
|
|
}
|
|
|
|
|
|
|
|
|
|
#define NUMBER(a) (sizeof(a)/sizeof(a[0]))
|
|
|
|
extern void grsrnd_(long*);
|
|
extern void gribex_(long*, long*, long*, double*, long*, double*, long*, double*,
|
|
long*, char *, long*, long*, const char*, long*,long );
|
|
|
|
static int cgribex(
|
|
long miss,
|
|
|
|
long ksec0[],long ksec1[],
|
|
long ksec2[],double rsec2[],
|
|
long ksec3[],double rsec3[],
|
|
long ksec4[],double rsec4[],
|
|
long sec4len,
|
|
|
|
char *buffer,
|
|
long *buflen,
|
|
char *op)
|
|
{
|
|
long ret = 1;
|
|
long in_bufsize = (*buflen) / sizeof(long);
|
|
long out_bufsize = in_bufsize;
|
|
long round = 0;
|
|
|
|
rsec3[1] = miss,
|
|
ksec3[1] = miss,
|
|
|
|
grsrnd_(&round); /* no rounding !!! */
|
|
|
|
gribex_(ksec0, ksec1, ksec2, rsec2, ksec3, rsec3, ksec4, rsec4,
|
|
&sec4len, buffer, &in_bufsize, &out_bufsize, op, &ret,1 );
|
|
|
|
*buflen = out_bufsize * sizeof(long);
|
|
|
|
return ret;
|
|
}
|
|
|
|
static double proc_cpu()
|
|
{
|
|
struct rusage rup;
|
|
|
|
if(getrusage(RUSAGE_SELF,&rup) != -1)
|
|
{
|
|
return (rup.ru_utime.tv_sec + rup.ru_utime.tv_usec / 1000000.0 +
|
|
rup.ru_stime.tv_sec + rup.ru_stime.tv_usec / 1000000.0);
|
|
}
|
|
return clock()/(double)CLOCKS_PER_SEC;
|
|
}
|
|
|
|
int compare_values(double* v1,double* v2,int count,const char *what)
|
|
{
|
|
int i,j=0;
|
|
double maxError = 0;
|
|
for(i = 0; i < count ;i++)
|
|
{
|
|
double e = err(v1[i],v2[i]);
|
|
if(e > maxError)
|
|
{
|
|
j = i;
|
|
maxError = e;
|
|
}
|
|
}
|
|
fprintf(stdout,"value no %d out of %d %s v1=%g v2=%g abs=%.30e rel=%g\n",j,count,what,v1[j], v2[j],v1[j]-v2[j],maxError);
|
|
if(maxError > maxRelativeError)
|
|
{
|
|
fprintf(stderr,"value no %d out of %d %s v1=%g v2=%g abs=%g rel=%g\n",j,count,what,v1[j], v2[j],v1[j]-v2[j],maxError);
|
|
return 1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
void get_max_relative_error(double* v1,double* v2,int count,double* maxError,int* id) {
|
|
int i;
|
|
*id=0;
|
|
*maxError = 0;
|
|
for(i = 0; i < count ;i++) {
|
|
double e = err(v1[i],v2[i]);
|
|
if(e > *maxError)
|
|
{
|
|
*id = i;
|
|
*maxError = e;
|
|
}
|
|
}
|
|
}
|
|
|
|
int compare_with_gribex(grib_handle* g)
|
|
{
|
|
|
|
long ed = 0;
|
|
int err = 0;
|
|
long missing = 0;
|
|
long is_int;
|
|
long miss;
|
|
int i;
|
|
size_t count;
|
|
long sec4len;
|
|
long len;
|
|
double now,gtime,atime;
|
|
size_t size;
|
|
const void *buffer;
|
|
long ksec0[ISECTION_0];
|
|
long ksec1[ISECTION_1];
|
|
long ksec2[ISECTION_2];
|
|
long ksec3[ISECTION_3];
|
|
long ksec4[ISECTION_4];
|
|
|
|
double rsec2[RSECTION_2];
|
|
double rsec3[RSECTION_3];
|
|
double *rsec4;
|
|
double *values;
|
|
|
|
|
|
if((err = grib_get_long(g,"editionNumber",&ed)))
|
|
{
|
|
fprintf(stderr,"Cannot get editionNumber %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if(ed > 1)
|
|
return 0;
|
|
|
|
|
|
if((err = grib_get_long(g,"missingValue",&miss)))
|
|
{
|
|
fprintf(stderr,"Cannot get missingValue %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if((err = grib_get_long(g,"integerPointValues",&is_int)))
|
|
{
|
|
fprintf(stderr,"Cannot get integerPointValues %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if((err = grib_get_size(g,"values",&count)))
|
|
{
|
|
fprintf(stderr,"Cannot get number of values %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
printf("Number of values %ld\n",(long)count);
|
|
|
|
/* printf("number of values is %d\n",count); */
|
|
|
|
rsec4 = (double*)malloc(sizeof(double)*(count+100000));
|
|
values = (double*)malloc(sizeof(double)*(count+100000));
|
|
|
|
sec4len = count +100000;
|
|
if((err = grib_get_message(g,&buffer,&size)))
|
|
{
|
|
fprintf(stderr,"Cannot get message %d\n",err);
|
|
goto foo;
|
|
}
|
|
len = size;
|
|
|
|
now = proc_cpu();
|
|
err = cgribex(
|
|
miss,
|
|
ksec0,ksec1,ksec2,rsec2,ksec3,rsec3,
|
|
ksec4,rsec4,sec4len,
|
|
(char*)buffer,&len,"D");
|
|
|
|
gtime = proc_cpu() - now;
|
|
|
|
if(err > 0)
|
|
{
|
|
char *e = "unknown";
|
|
for(i = 0; i < NUMBER(gribex_errors); i++)
|
|
if(gribex_errors[i].code == err)
|
|
e = gribex_errors[i].message;
|
|
|
|
fprintf(stderr,"GRIBEX error %d: %s\n",err,e);
|
|
fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,(long)ksec4[0]);
|
|
|
|
goto foo;
|
|
}
|
|
|
|
if(ksec4[0] == 0){
|
|
fprintf(stderr,"No GRIB processed - no values in KSEC4\n");
|
|
goto foo;
|
|
}
|
|
|
|
|
|
now = proc_cpu();
|
|
if((err = grib_get_double_array(g,"values",values,&count)))
|
|
{
|
|
fprintf(stderr,"Cannot get values %d (%s) count=%ld\n",err,grib_get_error_message(err),(long)count);
|
|
return err;
|
|
}
|
|
atime = proc_cpu() - now;
|
|
|
|
if( ksec4[0] < 0)
|
|
{
|
|
missing = 1;
|
|
ksec4[0] = -ksec4[0];
|
|
}
|
|
|
|
if(count != ksec4[0])
|
|
{
|
|
fprintf(stderr,"Count mismatch api=%ld gribex=%ld\n",
|
|
(long)count,ksec4[0]);
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
#if 0
|
|
printf("Values: %d gribex: %g api:%g\n",count,gtime,atime);
|
|
|
|
if(atime > gtime*1.1)
|
|
{
|
|
fprintf(stderr,"Too slow: gribex: %g api:%g\n",gtime,atime);
|
|
exit(1);
|
|
}
|
|
#endif
|
|
|
|
if(missing)
|
|
{
|
|
for(i = 0; i < count ;i++)
|
|
{
|
|
if(values[i] != miss)
|
|
{
|
|
fprintf(stderr,"Missing field: value %g not set to %ld\n",values[i],miss);
|
|
err = 1;
|
|
goto foo;
|
|
|
|
}
|
|
}
|
|
}
|
|
else if(is_int)
|
|
{
|
|
long long* lval = (long long*)rsec4;
|
|
for(i = 0; i < count ;i++)
|
|
rsec4[i] = lval[i];
|
|
|
|
if(compare_values(rsec4,values,count,"gribex"))
|
|
{
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(compare_values(rsec4,values,count,"gribex"))
|
|
{
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
}
|
|
|
|
foo:
|
|
free(rsec4);
|
|
free(values);
|
|
|
|
return err;
|
|
}
|
|
|
|
int compare_two_gribex(grib_handle* g)
|
|
{
|
|
|
|
long ed = 0;
|
|
int err = 0;
|
|
long is_int;
|
|
long miss;
|
|
int i;
|
|
size_t count;
|
|
long sec4len;
|
|
long len;
|
|
size_t size;
|
|
const void *buffer;
|
|
long ksec0[ISECTION_0];
|
|
long ksec1[ISECTION_1];
|
|
long ksec2[ISECTION_2];
|
|
long ksec3[ISECTION_3];
|
|
long ksec4[ISECTION_4];
|
|
|
|
double rsec2[RSECTION_2];
|
|
double rsec3[RSECTION_3];
|
|
double *rsec4_1 =0 ;
|
|
double *rsec4_2 =0 ;
|
|
double *rsec4_3 =0 ;
|
|
char *b;
|
|
|
|
|
|
if((err = grib_get_long(g,"editionNumber",&ed)))
|
|
{
|
|
fprintf(stderr,"Cannot get editionNumber %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if(ed > 1)
|
|
return 0;
|
|
|
|
|
|
if((err = grib_get_long(g,"missingValue",&miss)))
|
|
{
|
|
fprintf(stderr,"Cannot get missingValue %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if((err = grib_get_long(g,"integerPointValues",&is_int)))
|
|
{
|
|
fprintf(stderr,"Cannot get integerPointValues %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
if(is_int)
|
|
return 0;
|
|
|
|
if((err = grib_get_size(g,"values",&count)))
|
|
{
|
|
fprintf(stderr,"Cannot get number of values %d\n",err);
|
|
return err;
|
|
}
|
|
|
|
printf("Number of values %ld\n",(long)count);
|
|
|
|
/* printf("number of values is %d\n",count); */
|
|
|
|
rsec4_1 = (double*)malloc(sizeof(double)*(count+100000));
|
|
rsec4_2 = (double*)malloc(sizeof(double)*(count+100000));
|
|
rsec4_3 = (double*)malloc(sizeof(double)*(count+100000));
|
|
b = (char*)malloc(sizeof(char)*(count+100000));
|
|
|
|
sec4len = count +100000;
|
|
if((err = grib_get_message(g,&buffer,&size)))
|
|
{
|
|
fprintf(stderr,"Cannot get message %d\n",err);
|
|
goto foo;
|
|
}
|
|
len = size;
|
|
|
|
err = cgribex(
|
|
miss,
|
|
ksec0,ksec1,ksec2,rsec2,ksec3,rsec3,
|
|
ksec4,rsec4_1,sec4len,
|
|
(char*)buffer,&len,"D");
|
|
|
|
|
|
if(err > 0)
|
|
{
|
|
char *e = "unknown";
|
|
for(i = 0; i < NUMBER(gribex_errors); i++)
|
|
if(gribex_errors[i].code == err)
|
|
e = gribex_errors[i].message;
|
|
|
|
fprintf(stderr,"GRIBEX error %d: %s\n",err,e);
|
|
fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]);
|
|
|
|
goto foo;
|
|
}
|
|
|
|
memcpy(rsec4_3,rsec4_1,sizeof(double)*(count+100000));
|
|
|
|
if(ksec4[0] == 0){
|
|
fprintf(stderr,"No GRIB processed - no values in KSEC4\n");
|
|
goto foo;
|
|
}
|
|
|
|
if( ksec4[0] < 0)
|
|
goto foo;
|
|
|
|
if(count != ksec4[0])
|
|
{
|
|
fprintf(stderr,"Count mismatch api=%ld gribex=%ld\n",
|
|
(long)count,ksec4[0]);
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
len = count +100000;
|
|
err = cgribex(
|
|
miss,
|
|
ksec0,ksec1,ksec2,rsec2,ksec3,rsec3,
|
|
ksec4,rsec4_1,sec4len,
|
|
(char*)b,&len,"C");
|
|
|
|
|
|
if(err > 0)
|
|
{
|
|
char *e = "unknown";
|
|
for(i = 0; i < NUMBER(gribex_errors); i++)
|
|
if(gribex_errors[i].code == err)
|
|
e = gribex_errors[i].message;
|
|
|
|
fprintf(stderr,"GRIBEX error %d: %s\n",err,e);
|
|
fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]);
|
|
|
|
goto foo;
|
|
}
|
|
|
|
sec4len = count +100000;
|
|
err = cgribex(
|
|
miss,
|
|
ksec0,ksec1,ksec2,rsec2,ksec3,rsec3,
|
|
ksec4,rsec4_2,sec4len,
|
|
(char*)buffer,&len,"D");
|
|
|
|
|
|
if(err > 0)
|
|
{
|
|
char *e = "unknown";
|
|
for(i = 0; i < NUMBER(gribex_errors); i++)
|
|
if(gribex_errors[i].code == err)
|
|
e = gribex_errors[i].message;
|
|
|
|
fprintf(stderr,"GRIBEX error %d: %s\n",err,e);
|
|
fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]);
|
|
|
|
goto foo;
|
|
}
|
|
|
|
if(compare_values(rsec4_3,rsec4_2,count,"gribex2"))
|
|
{
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
foo:
|
|
free(rsec4_1);
|
|
free(rsec4_2);
|
|
free(rsec4_3);
|
|
free(b);
|
|
|
|
return err;
|
|
}
|
|
|
|
int compare_lists(grib_handle* g,grib_handle *h,const char* what)
|
|
{
|
|
int err;
|
|
size_t count = 0;
|
|
size_t ncount = 0;
|
|
double *values1;
|
|
double *values2;
|
|
|
|
int e1 = grib_get_size(g,what,&count);
|
|
int e2 = grib_get_size(h,what,&ncount);
|
|
|
|
|
|
if(e1 == GRIB_NOT_FOUND && e2 == GRIB_NOT_FOUND)
|
|
return 0;
|
|
|
|
if(e1)
|
|
{
|
|
fprintf(stderr,"Cannot get number of '%s' %d\n",what,e1);
|
|
return 1;
|
|
}
|
|
|
|
if(e2)
|
|
{
|
|
fprintf(stderr,"Cannot get number of '%s' %d\n",what,e2);
|
|
return 1;
|
|
}
|
|
|
|
if(ncount != count)
|
|
{
|
|
fprintf(stderr,"'%s' count mismatch %ld -> %ld\n",what,(long)count,(long)ncount);
|
|
return 1;
|
|
}
|
|
|
|
values1 = (double*)malloc(sizeof(double)*(count));
|
|
values2 = (double*)malloc(sizeof(double)*(count));
|
|
|
|
|
|
ncount = count;
|
|
if((err = grib_get_double_array(g,what,values1,&count)))
|
|
{
|
|
fprintf(stderr,"Cannot get %s1 %d (%s) count=%ld\n",what,err,grib_get_error_message(err),(long)count);
|
|
goto foo;
|
|
}
|
|
|
|
if(ncount != count)
|
|
{
|
|
fprintf(stderr,"%s count changes %ld -> %ld\n",what,(long)count,(long)ncount);
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
|
|
ncount = count;
|
|
if((err = grib_get_double_array(h,what,values2,&count)))
|
|
{
|
|
fprintf(stderr,"Cannot get %s2 %d (%s) count=%ld\n",what,err,grib_get_error_message(err),(long)ncount);
|
|
goto foo;
|
|
}
|
|
|
|
if(ncount != count)
|
|
{
|
|
fprintf(stderr,"%s count changes %ld -> %ld\n",what,(long)count,(long)ncount);
|
|
err = 1;
|
|
goto foo;
|
|
}
|
|
|
|
err = compare_values(values1,values2,count,what);
|
|
|
|
foo:
|
|
free(values2);
|
|
free(values1);
|
|
|
|
return err;
|
|
}
|
|
|
|
int compare_with_handle(grib_handle* g,grib_handle *h)
|
|
{
|
|
|
|
double s1,d1,r1,l1;
|
|
double s2,d2,r2,l2;
|
|
|
|
|
|
l1 = s1 = r1 = d1 = 1e-100;
|
|
grib_get_double(g,"binaryScaleFactor",&s1);
|
|
grib_get_double(g,"referenceValue",&r1);
|
|
grib_get_double(g,"decimalScaleFactor",&d1);
|
|
grib_get_double(g,"laplacianOperator",&l1);
|
|
|
|
|
|
l2 = s2 = r2 = d2 = 1e-100;
|
|
grib_get_double(h,"binaryScaleFactor",&s2);
|
|
grib_get_double(h,"referenceValue",&r2);
|
|
grib_get_double(h,"decimalScaleFactor",&d2);
|
|
grib_get_double(g,"laplacianOperator",&l2);
|
|
|
|
if(s1 != s2)
|
|
{
|
|
fprintf(stderr,
|
|
"Scale factor 1 = %g "
|
|
"Scale factor 2 = %g diff=%g"
|
|
"\n" ,s1,s2,s1-s2);
|
|
}
|
|
|
|
if(d1 != d2)
|
|
{
|
|
fprintf(stderr,
|
|
"Decimal Scale factor 1 = %g "
|
|
"Decimal Scale factor 2 = %g diff=%g"
|
|
"\n" ,d1,d2,d1-d2);
|
|
}
|
|
|
|
if(r1 != r2)
|
|
{
|
|
fprintf(stderr,
|
|
"Ref. value 1 = %g "
|
|
"Ref. value 2 = %g diff=%g"
|
|
"\n" ,r1,r2,r1-r2);
|
|
}
|
|
|
|
if(l1 != l2)
|
|
{
|
|
fprintf(stderr,
|
|
"Laplacian 1 = %g "
|
|
"Laplacian 2 = %g diff=%g"
|
|
"\n" ,l1,l2,l1-l2);
|
|
}
|
|
|
|
|
|
|
|
if(compare_lists(g,h,"values"))
|
|
return 1;
|
|
|
|
if(compare_lists(g,h,"pv"))
|
|
return 1;
|
|
|
|
if(compare_lists(g,h,"pl"))
|
|
return 1;
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
int save_message(grib_handle* h,const char* path)
|
|
{
|
|
FILE* out = fopen(path,"w");
|
|
const void *p;
|
|
size_t len;
|
|
|
|
if(!out) {
|
|
perror(path);
|
|
return -1;
|
|
}
|
|
|
|
printf("save to %s\n",path);
|
|
|
|
grib_get_message(h,&p,&len);
|
|
|
|
if(fwrite(p,1,len,out) != len)
|
|
{
|
|
perror(path);
|
|
fclose(out);
|
|
return -1;
|
|
}
|
|
|
|
if(fclose(out))
|
|
{
|
|
perror(path);
|
|
return -1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
const char* problem(const char* path,const char* name,int field)
|
|
{
|
|
static char buf[1024];
|
|
const char *p = path;
|
|
const char *q = path;
|
|
|
|
while(*p)
|
|
{
|
|
if(*p == '/') q=p+1;
|
|
p++;
|
|
}
|
|
|
|
sprintf(buf,"/vol/dataserv/grib_samples/problems/%s_%s_%d.grib",q,name,field);
|
|
return buf;
|
|
}
|