Skip to content

GitLab

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

grib2: Added support for variables with different perturbationNumber.

parent 59c4fefe
Hide whitespace changes
Inline Side-by-side
Showing with 69 additions and 49 deletions
ChangeLog
View file @ ac1e7627
......@@ -4,6 +4,10 @@
* using EXSE library version 1.4.1
* Version 1.9.6 released
2018-11-16 Uwe Schulzweida
* grib2: Added support for variables with different perturbationNumber
2018-11-09 Uwe Schulzweida
* Added check for unsupported NetCDF4/HDF5 library combination (NetCDF <= 4.4.0 with libhdf5 >= 1.10.0)
......
src/cdi_int.h
View file @ ac1e7627
......@@ -153,15 +153,8 @@ typedef struct
int ilevel;
int ilevel2;
int ltype;
short perturbationNumber;
short tsteptype;
#ifdef HIRLAM_EXTENSIONS
// NOTE: tsteptype MUST be part of attributes used to compare variables!
// Modern NWP models (HARMONIE, HIRLAM) use timeRangeIndicator to specify
// if the field is instantanous or accumulated.
// Both types are typically in the same GRIB-file.
// (181; 105, 0, timeRangeIndicator=0) .. instantanous rain
// (181; 105, 0, timeRangeIndicator=4) .. accumulated rain .. both can be in the same grib file
#endif // HIRLAM_EXTENSIONS
short varID;
short levelID;
short used;
......
src/stream_cgribex.c
View file @ ac1e7627
......@@ -564,7 +564,7 @@ void cgribexAddRecord(stream_t *streamptr, cgribexrec_t *cgribexp, int param, si
varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, 0, 0,
datatype, &varID, &levelID, tsteptype, numavg, leveltype, -1,
NULL, NULL, NULL, NULL, NULL, NULL);
NULL, NULL, NULL, NULL, NULL, NULL, 0);
record->varID = (short)varID;
record->levelID = (short)levelID;
......
src/stream_ext.c
View file @ ac1e7627
......@@ -198,7 +198,7 @@ void extAddRecord(stream_t *streamptr, int param, int level, size_t xysize,
int varID, levelID = 0;
varAddRecord(recID, param, gridID, leveltype, 0, level, 0, 0, 0,
extInqDatatype(prec, number), &varID, &levelID, TSTEP_INSTANT, 0, 0, -1,
NULL, NULL, NULL, NULL, NULL, NULL);
NULL, NULL, NULL, NULL, NULL, NULL, 0);
record->varID = (short)varID;
record->levelID = (short)levelID;
......
src/stream_gribapi.c
View file @ ac1e7627
......@@ -26,25 +26,18 @@ extern int CDI_inventory_mode;
static const var_tile_t dummy_tiles = { 0, -1, -1, -1, -1, -1 };
typedef struct {
int param;
int level1;
int level2;
int ltype;
int tsteptype;
int perturbationNumber;
size_t gridsize;
#ifdef HIRLAM_EXTENSIONS
// NOTE: tsteptype MUST be part of attributes used to compare variables!
// Modern NWP models (HARMONIE, HIRLAM) use timeRangeIndicator to specify
// if the field is instantanous or accumulated.
// Both types are typically in the same GRIB-file.
// (181; 105, 0, timeRangeIndicator=0) .. instantanous rain
// (181; 105, 0, timeRangeIndicator=4) .. accumulated rain .. both can be in the same grib file
#endif // HIRLAM_EXTENSIONS
char name[32];
var_tile_t tiles;
} compvar2_t;
......@@ -411,6 +404,27 @@ void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t len
else if ( length == 2 && memcmp(string, "~", length) == 0 ) string[0] = 0;
}
static
int gribapiGetEnsembleInfo(grib_handle *gh, long *typeOfEnsembleForecast, long *numberOfForecastsInEnsemble, long *perturbationNumber)
{
int status = 0;
if (grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast) == 0)
{
GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble), 0);
GRIB_CHECK(grib_get_long(gh, "perturbationNumber", perturbationNumber), 0);
if (*perturbationNumber > 0) status = 1;
}
if (status == 0)
{
*typeOfEnsembleForecast = 0;
*perturbationNumber = 0;
*numberOfForecastsInEnsemble = 0;
}
return status;
}
static
void gribapiGetKeys(grib_handle *gh, int varID)
{
......@@ -446,18 +460,13 @@ void gribapiGetKeys(grib_handle *gh, int varID)
Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
*/
long typeOfEnsembleForecast = 0;
if ( grib_get_long(gh, "typeOfEnsembleForecast", &typeOfEnsembleForecast) == 0 )
long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
if ( perturbationNumber > 0 )
{
long perturbationNumber = 0, numberOfForecastsInEnsemble = 0;
GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", &numberOfForecastsInEnsemble), 0);
GRIB_CHECK(grib_get_long(gh, "perturbationNumber", &perturbationNumber), 0);
if ( perturbationNumber > 0 )
{
varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
}
varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
}
long section2Length;
......@@ -584,7 +593,7 @@ static
void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
size_t recsize, off_t position, int datatype, int comptype, const char *varname,
int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
const var_tile_t *tiles, int lread_additional_keys)
const var_tile_t *tiles, int perturbationNumber, int lread_additional_keys)
{
char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
......@@ -609,6 +618,8 @@ void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
record->gridsize = gribapiGetGridsize(gh);
record->tiles = tiles ? *tiles : dummy_tiles;
record->perturbationNumber = perturbationNumber;
strncpy(record->varname, varname, sizeof(record->varname)-1);
record->varname[sizeof(record->varname) - 1] = 0;
......@@ -686,7 +697,7 @@ void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
int tile_index = 0, varID;
varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
varname, stdname, longname, units, tiles, &tile_index);
varname, stdname, longname, units, tiles, &tile_index, perturbationNumber);
record->varID = (short)varID;
record->levelID = (short)levelID;
......@@ -757,7 +768,7 @@ void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
}
static
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype,
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype, long perturbationNumber,
size_t gridsize, char *name, var_tile_t tiles_data)
{
compvar2_t compVar;
......@@ -771,6 +782,7 @@ compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int t
compVar.level2 = level2;
compVar.ltype = leveltype;
compVar.tsteptype = tsteptype;
compVar.perturbationNumber = (int)perturbationNumber;
compVar.gridsize = gridsize;
//memset(compVar.name, 0, maxlen);
memcpy(compVar.name, name, len);
......@@ -789,6 +801,7 @@ int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
compVar0.level2 = record.ilevel2;
compVar0.ltype = record.ltype;
compVar0.tsteptype = record.tsteptype;
compVar0.perturbationNumber = record.perturbationNumber;
compVar0.gridsize = record.gridsize;
memcpy(compVar0.name, record.varname, sizeof(compVar.name));
......@@ -962,8 +975,6 @@ int gribapiScanTimestep1(stream_t * streamptr)
varname[0] = 0;
gribapiGetString(gh, "shortName", varname, sizeof(varname));
int tsteptype = gribapiGetTsteptype(gh);
int64_t vdate = 0;
int vtime = 0;
gribapiGetValidityDateTime(gh, &vdate, &vtime);
......@@ -982,10 +993,15 @@ int gribapiScanTimestep1(stream_t * streamptr)
}
}
long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
if ( nrecs )
{
int tsteptype = gribapiGetTsteptype(gh);
size_t gridsize = gribapiGetGridsize(gh);
checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, tiles),
checkTimeResult result = checkTime(streamptr,
gribapiVarSet(param, level1, level2, leveltype1, tsteptype, perturbationNumber, gridsize, varname, tiles),
&datetime, &datetime0);
if ( result == CHECKTIME_STOP )
{
......@@ -1031,7 +1047,7 @@ int gribapiScanTimestep1(stream_t * streamptr)
var_tile_t *ptiles = NULL;
if ( memcmp(&tiles, &dummy_tiles, sizeof(var_tile_t)) != 0 ) ptiles = &tiles;
gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, ptiles, 1);
leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, ptiles, perturbationNumber, 1);
grib_handle_delete(gh);
gh = NULL;
......@@ -1200,7 +1216,6 @@ int gribapiScanTimestep2(stream_t * streamptr)
datetime0.time = vtime;
}
int tsteptype = gribapiGetTsteptype(gh);
/*
if ( ISEC1_AvgNum )
{
......@@ -1220,8 +1235,11 @@ int gribapiScanTimestep2(stream_t * streamptr)
.time = vtime
};
long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
int tsteptype = gribapiGetTsteptype(gh);
size_t gridsize = gribapiGetGridsize(gh);
compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, tiles);
compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, perturbationNumber, gridsize, varname, tiles);
for ( recID = 0; recID < nrecords; recID++ )
if ( gribapiVarCompare(compVar, streamptr->tsteps[tsID].records[recID], 0) == 0 ) break;
......@@ -1450,10 +1468,11 @@ int gribapiScanTimestep(stream_t * streamptr)
.time = vtime
};
long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
int tsteptype = gribapiGetTsteptype(gh);
size_t gridsize = gribapiGetGridsize(gh);
compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, tiles);
compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, perturbationNumber, gridsize, varname, tiles);
for ( vrecID = 0; vrecID < nrecs; vrecID++ )
{
......
src/stream_ieg.c
View file @ ac1e7627
......@@ -663,7 +663,7 @@ void iegAddRecord(stream_t *streamptr, int param, int *pdb, int *gdb, double *vc
int levelID = 0;
varAddRecord(recID, param, gridID, leveltype, lbounds, level1, level2, 0, 0,
datatype, &varID, &levelID, TSTEP_INSTANT, 0, 0, -1,
NULL, NULL, NULL, NULL, NULL, NULL);
NULL, NULL, NULL, NULL, NULL, NULL, 0);
record->varID = (short)varID;
record->levelID = (short)levelID;
......
src/stream_record.c
View file @ ac1e7627
......@@ -26,6 +26,7 @@ void recordInitEntry(record_t *record)
record->ilevel = CDI_UNDEFID;
record->used = false;
record->tsteptype = CDI_UNDEFID;
record->perturbationNumber = 0;
record->varID = CDI_UNDEFID;
record->levelID = CDI_UNDEFID;
memset(record->varname, 0, sizeof(record->varname));
......
src/stream_srv.c
View file @ ac1e7627
......@@ -204,7 +204,7 @@ void srv_add_record(stream_t *streamptr, int param, int level, size_t xsize, siz
int varID;
varAddRecord(recID, param, gridID, leveltype, 0, level, 0, 0, 0,
datatype, &varID, &levelID, TSTEP_INSTANT, 0, 0, -1,
NULL, NULL, NULL, NULL, NULL, NULL);
NULL, NULL, NULL, NULL, NULL, NULL, 0);
xassert(varID <= SHRT_MAX && levelID <= SHRT_MAX);
record->varID = (short)varID;
......
src/varscan.c
View file @ ac1e7627
......@@ -49,6 +49,7 @@ typedef struct
int param;
int prec;
int tsteptype;
int perturbationNumber;
int gridID;
int zaxistype;
int ltype1; /* GRIB first level type */
......@@ -98,6 +99,7 @@ void paramInitEntry(unsigned varID, int param)
vartable[varID].param = param;
vartable[varID].prec = 0;
vartable[varID].tsteptype = TSTEP_INSTANT;
vartable[varID].perturbationNumber = 0;
vartable[varID].timave = 0;
vartable[varID].gridID = CDI_UNDEFID;
vartable[varID].zaxistype = 0;
......@@ -125,10 +127,9 @@ void paramInitEntry(unsigned varID, int param)
vartable[varID].tiles = NULL;
}
/* Test if a variable specified by the given meta-data has already
* been registered in "vartable". */
/* Test if a variable specified by the given meta-data has already been registered in "vartable". */
static unsigned
varGetEntry(int param, int gridID, int zaxistype, int ltype1, int tsteptype, const char *name, const var_tile_t *tiles)
varGetEntry(int param, int gridID, int zaxistype, int ltype1, int tsteptype, const char *name, const var_tile_t *tiles, int perturbationNumber)
{
for ( unsigned varID = 0; varID < varTablesize; varID++ )
{
......@@ -145,6 +146,7 @@ varGetEntry(int param, int gridID, int zaxistype, int ltype1, int tsteptype, con
if ( (vartable[varID].zaxistype == zaxistype) &&
(vartable[varID].ltype1 == ltype1 ) &&
(vartable[varID].tsteptype == tsteptype) &&
(vartable[varID].perturbationNumber == perturbationNumber) &&
(vartable[varID].gridID == gridID ) &&
(vt_no_of_tiles == no_of_tiles) )
{
......@@ -425,10 +427,10 @@ void varAddRecord(int recID, int param, int gridID, int zaxistype, int lbounds,
int level1, int level2, int level_sf, int level_unit, int prec,
int *pvarID, int *plevelID, int tsteptype, int numavg, int ltype1, int ltype2,
const char *name, const char *stdname, const char *longname, const char *units,
const var_tile_t *tiles, int *tile_index)
const var_tile_t *tiles, int *tile_index, int perturbationNumber)
{
unsigned varID = (cdiSplitLtype105 != 1 || zaxistype != ZAXIS_HEIGHT) ?
varGetEntry(param, gridID, zaxistype, ltype1, tsteptype, name, tiles) : (unsigned) CDI_UNDEFID;
varGetEntry(param, gridID, zaxistype, ltype1, tsteptype, name, tiles, perturbationNumber) : (unsigned) CDI_UNDEFID;
if ( varID == (unsigned) CDI_UNDEFID )
{
......@@ -442,6 +444,7 @@ void varAddRecord(int recID, int param, int gridID, int zaxistype, int lbounds,
vartable[varID].level_sf = level_sf;
vartable[varID].level_unit = level_unit;
vartable[varID].tsteptype = tsteptype;
vartable[varID].perturbationNumber = perturbationNumber;
if ( numavg ) vartable[varID].timave = 1;
......
src/varscan.h
View file @ ac1e7627
......@@ -10,7 +10,7 @@ void varAddRecord(int recID, int param, int gridID, int zaxistype, int lbounds,
int level1, int level2, int level_sf, int level_unit, int prec,
int *pvarID, int *plevelID, int tsteptype, int numavg, int ltype1, int ltype2,
const char *name, const char *stdname, const char *longname, const char *units,
const var_tile_t *tiles, int *tile_index);
const var_tile_t *tiles, int *tile_index, int perturbationNumber);
void varDefVCT(size_t vctsize, double *vctptr);
void varDefZAxisReference(int nlev, int nvgrid, unsigned char uuid[CDI_UUID_SIZE]);
......
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