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

replace vlistDefVarCode() by vlistDefVarParam()

parent 0f8da6c7
......@@ -3,6 +3,10 @@
* using CDI library version 1.7.0
* Version 1.7.0 released
2015-06-03 Uwe Schulzweida
* setpartabn: replace vlistDefVarCode() by vlistDefVarParam()
2015-06-01 Uwe Schulzweida
* vertvar, vertstd: changed to weighted var/std if layer bounds are available
......
......@@ -105,7 +105,6 @@ void *CDIwrite(void *argument)
unsigned int seed = 1;
const char *gridfile;
char sinfo[64];
char *envstr;
off_t nvalues = 0;
double file_size = 0, data_size = 0;
double tw, tw0, t0, twsum = 0;
......@@ -117,7 +116,7 @@ void *CDIwrite(void *argument)
cdoInitialize(argument);
envstr = getenv("MEMTYPE");
char *envstr = getenv("MEMTYPE");
if ( envstr )
{
if ( strcmp(envstr, "float") == 0 ) memtype = MEMTYPE_FLOAT;
......@@ -212,7 +211,7 @@ void *CDIwrite(void *argument)
for ( i = 0; i < nvars; ++i )
{
varID = vlistDefVar(vlistID, gridID, zaxisID, TSTEP_INSTANT);
vlistDefVarCode(vlistID, varID, varID+1);
vlistDefVarParam(vlistID, varID, cdiEncodeParam(varID+1, 255, 255));
// vlistDefVarName(vlistID, varID, );
}
......
......@@ -122,19 +122,13 @@ void pl_index(int *kmax, int *kmin, double pmax, double pmin, long nlevs, double
void *Cloudlayer(void *argument)
{
int streamID1, streamID2;
int vlistID1, vlistID2;
int taxisID1, taxisID2;
int gridID, zaxisID, tsID;
int surfaceID;
int nlevel, nhlev, nlevs, nrecs, recID, code;
int gridID, zaxisID;
int nlevel, nlevs, nrecs, recID, code;
int varID, levelID;
int zrev = FALSE;
int nvars;
int i;
int offset;
int nmiss;
int aclcac_code;
int aclcacID = -1;
int nvars2 = 0;
int aclcac_code_found = 0;
......@@ -142,9 +136,7 @@ void *Cloudlayer(void *argument)
char varname[CDI_MAX_NAME];
double sfclevel = 0;
double *plevs = NULL;
double *aclcac = NULL;
double *cloud[NVARS];
double missval;
double pmin = 0, pmax = 0;
cdoInitialize(argument);
......@@ -161,15 +153,15 @@ void *Cloudlayer(void *argument)
nvars2 = NVARS;
}
streamID1 = streamOpenRead(cdoStreamName(0));
int streamID1 = streamOpenRead(cdoStreamName(0));
vlistID1 = streamInqVlist(streamID1);
int vlistID1 = streamInqVlist(streamID1);
int gridsize = vlist_check_gridsize(vlistID1);
aclcac_code = 223;
int aclcac_code = 223;
nvars = vlistNvars(vlistID1);
int nvars = vlistNvars(vlistID1);
for ( varID = 0; varID < nvars; ++varID )
{
gridID = vlistInqVarGrid(vlistID1, varID);
......@@ -204,14 +196,14 @@ void *Cloudlayer(void *argument)
cdoAbort("Cloud cover (parameter 223) not found!");
}
missval = vlistInqVarMissval(vlistID1, aclcacID);
double missval = vlistInqVarMissval(vlistID1, aclcacID);
gridID = vlistInqVarGrid(vlistID1, aclcacID);
zaxisID = vlistInqVarZaxis(vlistID1, aclcacID);
nlevel = zaxisInqSize(zaxisID);
nhlev = nlevel+1;
int nhlev = nlevel+1;
aclcac = (double*) malloc(gridsize*nlevel*sizeof(double));
double *aclcac = (double*) malloc(gridsize*nlevel*sizeof(double));
for ( varID = 0; varID < nvars2; ++varID )
cloud[varID] = (double*) malloc(gridsize*sizeof(double));
......@@ -286,15 +278,15 @@ void *Cloudlayer(void *argument)
cdoAbort("Unsupported Z-Axis type!");
surfaceID = zaxisCreate(ZAXIS_SURFACE, 1);
int surfaceID = zaxisCreate(ZAXIS_SURFACE, 1);
zaxisDefLevels(surfaceID, &sfclevel);
vlistID2 = vlistCreate();
int vlistID2 = vlistCreate();
if ( nvars2 == 1 )
{
varID = vlistDefVar(vlistID2, gridID, surfaceID, TSTEP_INSTANT);
vlistDefVarCode(vlistID2, varID, 33);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(33, 128, 255));
vlistDefVarName(vlistID2, varID, "cld_lay");
vlistDefVarLongname(vlistID2, varID, "cloud layer");
vlistDefVarMissval(vlistID2, varID, missval);
......@@ -302,33 +294,33 @@ void *Cloudlayer(void *argument)
else
{
varID = vlistDefVar(vlistID2, gridID, surfaceID, TSTEP_INSTANT);
vlistDefVarCode(vlistID2, varID, 34);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(34, 128, 255));
vlistDefVarName(vlistID2, varID, "low_cld");
vlistDefVarLongname(vlistID2, varID, "low cloud");
vlistDefVarMissval(vlistID2, varID, missval);
varID = vlistDefVar(vlistID2, gridID, surfaceID, TSTEP_INSTANT);
vlistDefVarCode(vlistID2, varID, 35);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(35, 128, 255));
vlistDefVarName(vlistID2, varID, "mid_cld");
vlistDefVarLongname(vlistID2, varID, "mid cloud");
vlistDefVarMissval(vlistID2, varID, missval);
varID = vlistDefVar(vlistID2, gridID, surfaceID, TSTEP_INSTANT);
vlistDefVarCode(vlistID2, varID, 36);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(36, 128, 255));
vlistDefVarName(vlistID2, varID, "hih_cld");
vlistDefVarLongname(vlistID2, varID, "high cloud");
vlistDefVarMissval(vlistID2, varID, missval);
}
taxisID1 = vlistInqTaxis(vlistID1);
taxisID2 = taxisDuplicate(taxisID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
streamDefVlist(streamID2, vlistID2);
tsID = 0;
int tsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
taxisCopyTimestep(taxisID2, taxisID1);
......
......@@ -33,12 +33,12 @@ void init_vars(int vlistID, int gridID, int zaxisID, int nvars)
int code[] = {11, 17, 33, 34, 1, 2/*, 3*/};
char *name[] = {"temp", "depoint", "u", "v", "height", "pressure" /*, "station"*/};
char *units[] = {"Celsius", "", "m/s", "m/s", "m", "hPa" /*, ""*/};
int i, varID;
int varID;
for ( i = 0; i < nvars; ++i )
for ( int i = 0; i < nvars; ++i )
{
varID = vlistDefVar(vlistID, gridID, zaxisID, TSTEP_INSTANT);
vlistDefVarCode(vlistID, varID, code[i]);
vlistDefVarParam(vlistID, varID, cdiEncodeParam(code[i], 255, 255));
vlistDefVarName(vlistID, varID, name[i]);
vlistDefVarUnits(vlistID, varID, units[i]);
vlistDefVarDatatype(vlistID, varID, DATATYPE_FLT32);
......@@ -48,31 +48,26 @@ void init_vars(int vlistID, int gridID, int zaxisID, int nvars)
static
void init_data(int vlistID, int nvars, double *data[])
{
int varID, i, gridsize;
int gridsize;
double missval;
for ( varID = 0; varID < nvars; ++varID )
for ( int varID = 0; varID < nvars; ++varID )
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
missval = vlistInqVarMissval(vlistID, varID);
for ( i = 0; i < gridsize; ++i )
{
data[varID][i] = missval;
}
for ( int i = 0; i < gridsize; ++i ) data[varID][i] = missval;
}
}
static
void write_data(int streamID, int vlistID, int nvars, double *data[])
{
int i;
int varID;
int nmiss;
int gridsize;
double missval;
for ( varID = 0; varID < nvars; ++varID )
for ( int varID = 0; varID < nvars; ++varID )
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
missval = vlistInqVarMissval(vlistID, varID);
......@@ -80,7 +75,7 @@ void write_data(int streamID, int vlistID, int nvars, double *data[])
streamDefRecord(streamID, varID, 0);
nmiss = 0;
for ( i = 0; i < gridsize; ++i )
for ( int i = 0; i < gridsize; ++i )
if ( DBL_IS_EQUAL(data[varID][i], missval) ) nmiss++;
streamWriteRecord(streamID, data[varID], nmiss);
......@@ -90,11 +85,9 @@ void write_data(int streamID, int vlistID, int nvars, double *data[])
static
int getDate(const char *name)
{
int date = 0;
char *pname;
pname = strchr(name, '_');
char *pname = strchr(name, '_');
int date = 0;
if ( pname ) date = atoi(pname+1);
return(date);
......@@ -104,19 +97,11 @@ int getDate(const char *name)
void *Importobs(void *argument)
{
int operatorID;
char line[MAX_LINE_LEN];
int streamID;
int tsID;
int gridID, zaxisID, taxisID, vlistID;
int i, j;
int nvars = MAX_VARS;
int vdate = 0, vtime = 0;
int vdate0 = 0, vtime0 = 0;
int gridsize, xsize, ysize;
double *xvals = NULL, *yvals = NULL;
int vtime = 0;
double *data[MAX_VARS];
FILE *fp;
char dummy[32], station[32], datetime[32];
float lat, lon, height1, pressure, height2, value;
double latmin = 90, latmax = -90, lonmin = 360, lonmax = -360;
......@@ -129,23 +114,23 @@ void *Importobs(void *argument)
cdoOperatorAdd("import_obs", 0, 0, "grid description file or name");
operatorID = cdoOperatorID();
int operatorID = cdoOperatorID();
operatorInputArg(cdoOperatorEnter(operatorID));
gridID = cdoDefineGrid(operatorArgv()[0]);
int gridID = cdoDefineGrid(operatorArgv()[0]);
if ( gridInqType(gridID) != GRID_LONLAT )
cdoAbort("Unsupported grid type: %s", gridNamePtr(gridInqType(gridID)));
gridsize = gridInqSize(gridID);
xsize = gridInqXsize(gridID);
ysize = gridInqYsize(gridID);
int gridsize = gridInqSize(gridID);
int xsize = gridInqXsize(gridID);
int ysize = gridInqYsize(gridID);
// printf("gridsize=%d, xsize=%d, ysize=%d\n", gridsize, xsize, ysize);
xvals = (double*) malloc(gridsize*sizeof(double));
yvals = (double*) malloc(gridsize*sizeof(double));
double *xvals = (double*) malloc(gridsize*sizeof(double));
double *yvals = (double*) malloc(gridsize*sizeof(double));
gridInqXvals(gridID, xvals);
gridInqYvals(gridID, yvals);
......@@ -159,100 +144,100 @@ void *Importobs(void *argument)
grid_to_degree(units, gridsize, yvals, "grid center lat");
}
fp = fopen(cdoStreamName(0)->args, "r");
FILE *fp = fopen(cdoStreamName(0)->args, "r");
if ( fp == NULL ) { perror(cdoStreamName(0)->args); exit(EXIT_FAILURE); }
vdate = getDate(cdoStreamName(0)->args);
int vdate = getDate(cdoStreamName(0)->args);
if ( vdate <= 999999 ) vdate = vdate*100 + 1;
streamID = streamOpenWrite(cdoStreamName(1), cdoFiletype());
int streamID = streamOpenWrite(cdoStreamName(1), cdoFiletype());
zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
int zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
taxisID = taxisCreate(TAXIS_ABSOLUTE);
int taxisID = taxisCreate(TAXIS_ABSOLUTE);
vlistID = vlistCreate();
int vlistID = vlistCreate();
vlistDefTaxis(vlistID, taxisID);
for ( i = 0; i < nvars; ++i ) data[i] = (double*) malloc(gridsize*sizeof(double));
init_vars(vlistID, gridID, zaxisID, nvars);
streamDefVlist(streamID, vlistID);
int vdate0 = 0;
int vtime0 = 0;
//ntime = 0;
int tsID = 0;
while ( readline(fp, line, MAX_LINE_LEN) )
{
for ( i = 0; i < nvars; ++i ) data[i] = (double*) malloc(gridsize*sizeof(double));
init_vars(vlistID, gridID, zaxisID, nvars);
streamDefVlist(streamID, vlistID);
vdate0 = 0;
vtime0 = 0;
//ntime = 0;
tsID = 0;
while ( readline(fp, line, MAX_LINE_LEN) )
{
sscanf(line, "%s %s %s %g %g %g %d %g %g %g",
dummy, station, datetime, &lat, &lon, &height1, &code, &pressure, &height2, &value);
sscanf(datetime, "%d_%d", &vdate, &vtime);
if ( vdate != vdate0 || vtime != vtime0 )
{
if ( tsID > 0 ) write_data(streamID, vlistID, nvars, data);
vdate0 = vdate;
vtime0 = vtime;
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
taxisDefVdate(taxisID, vdate);
taxisDefVtime(taxisID, vtime);
streamDefTimestep(streamID, tsID);
sscanf(line, "%s %s %s %g %g %g %d %g %g %g",
dummy, station, datetime, &lat, &lon, &height1, &code, &pressure, &height2, &value);
sscanf(datetime, "%d_%d", &vdate, &vtime);
if ( vdate != vdate0 || vtime != vtime0 )
{
if ( tsID > 0 ) write_data(streamID, vlistID, nvars, data);
vdate0 = vdate;
vtime0 = vtime;
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
taxisDefVdate(taxisID, vdate);
taxisDefVtime(taxisID, vtime);
streamDefTimestep(streamID, tsID);
init_data(vlistID, nvars, data);
init_data(vlistID, nvars, data);
tsID++;
}
tsID++;
}
if ( lon < lonmin ) lonmin = lon;
if ( lon > lonmax ) lonmax = lon;
if ( lat < latmin ) latmin = lat;
if ( lat > latmax ) latmax = lat;
if ( lon < lonmin ) lonmin = lon;
if ( lon > lonmax ) lonmax = lon;
if ( lat < latmin ) latmin = lat;
if ( lat > latmax ) latmax = lat;
dy = yvals[1] - yvals[0];
for ( j = 0; j < ysize; ++j )
if ( lat >= (yvals[j]-dy/2) && lat < (yvals[j]+dy/2) ) break;
dy = yvals[1] - yvals[0];
for ( j = 0; j < ysize; ++j )
if ( lat >= (yvals[j]-dy/2) && lat < (yvals[j]+dy/2) ) break;
dx = xvals[1] - xvals[0];
if ( lon < (xvals[0] - dx/2) && lon < 0 ) lon+=360;
for ( i = 0; i < xsize; ++i )
if ( lon >= (xvals[i]-dx/2) && lon < (xvals[i]+dx/2) ) break;
dx = xvals[1] - xvals[0];
if ( lon < (xvals[0] - dx/2) && lon < 0 ) lon+=360;
for ( i = 0; i < xsize; ++i )
if ( lon >= (xvals[i]-dx/2) && lon < (xvals[i]+dx/2) ) break;
index = -1;
if ( code == 11 ) index = 0;
if ( code == 17 ) index = 1;
if ( code == 33 ) index = 2;
if ( code == 34 ) index = 3;
//printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
if ( i < xsize && j < ysize && index >= 0 )
{
pstation = station;
while (isalpha(*pstation)) pstation++;
// printf("station %s %d\n", pstation, atoi(pstation));
data[index][j*xsize+i] = value;
data[ 4][j*xsize+i] = height1;
data[ 5][j*xsize+i] = pressure;
// data[ 6][j*xsize+i] = atoi(pstation);
}
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
}
write_data(streamID, vlistID, nvars, data);
for ( i = 0; i < nvars; ++i ) free(data[i]);
index = -1;
if ( code == 11 ) index = 0;
if ( code == 17 ) index = 1;
if ( code == 33 ) index = 2;
if ( code == 34 ) index = 3;
//printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
if ( i < xsize && j < ysize && index >= 0 )
{
pstation = station;
while (isalpha(*pstation)) pstation++;
// printf("station %s %d\n", pstation, atoi(pstation));
data[index][j*xsize+i] = value;
data[ 4][j*xsize+i] = height1;
data[ 5][j*xsize+i] = pressure;
// data[ 6][j*xsize+i] = atoi(pstation);
}
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
}
printf("lonmin=%g, lonmax=%g, latmin=%g, latmax=%g\n", lonmin, lonmax, latmin, latmax);
write_data(streamID, vlistID, nvars, data);
for ( i = 0; i < nvars; ++i ) free(data[i]);
if ( cdoVerbose )
printf("lonmin=%g, lonmax=%g, latmin=%g, latmax=%g\n", lonmin, lonmax, latmin, latmax);
processDefVarNum(vlistNvars(vlistID), streamID);
......
......@@ -68,15 +68,12 @@ static int input_darray(int nval, double *array)
void *Input(void *argument)
{
int INPUT, INPUTSRV, INPUTEXT;
int operatorID;
int varID = 0;
int nlevs = 1;
int gridsize0 = 0, gridsize = 0;
int gridID = -1, zaxisID, vdate = 0, vtime = 0;
int nrecs;
int levelID;
int tsID, taxisID = 0;
int taxisID = 0;
int streamID = -1;
int vlistID = -1;
int nmiss = 0;
......@@ -91,11 +88,11 @@ void *Input(void *argument)
cdoInitialize(argument);
INPUT = cdoOperatorAdd("input", 0, 0, NULL);
INPUTSRV = cdoOperatorAdd("inputsrv", 0, 0, NULL);
INPUTEXT = cdoOperatorAdd("inputext", 0, 0, NULL);
int INPUT = cdoOperatorAdd("input", 0, 0, NULL);
int INPUTSRV = cdoOperatorAdd("inputsrv", 0, 0, NULL);
int INPUTEXT = cdoOperatorAdd("inputext", 0, 0, NULL);
operatorID = cdoOperatorID();
int operatorID = cdoOperatorID();
if ( operatorID == INPUT )
{
......@@ -105,9 +102,9 @@ void *Input(void *argument)
}
levels[0] = 0;
nrecs = 0;
int nrecs = 0;
tsID = 0;
int tsID = 0;
while ( TRUE )
{
if ( operatorID == INPUT )
......@@ -234,7 +231,7 @@ void *Input(void *argument)
vlistID = vlistCreate();
varID = vlistDefVar(vlistID, gridID, zaxisID, TSTEP_INSTANT);
vlistDefVarCode(vlistID, varID, code);
vlistDefVarParam(vlistID, varID, cdiEncodeParam(code, 255, 255));
missval = vlistInqVarMissval(vlistID, varID);
......
......@@ -257,7 +257,6 @@ void uv_to_p_grid(int nlon, int nlat, double *grid1x, double *grid1y,
}
void *Mrotuvb(void *argument)
{
int streamID1, streamID2, streamID3;
......
......@@ -38,7 +38,7 @@ void *Pressure(void *argument)
{
int mode;
enum {ECHAM_MODE, WMO_MODE};
int geop_code = 0, temp_code = 0, ps_code = 0, lsp_code = 0;
int ps_code = 0, lsp_code = 0;
int streamID2;
int vlistID2;
int recID, nrecs;
......@@ -46,10 +46,10 @@ void *Pressure(void *argument)
int tsID, varID, levelID;
int nvars;
int zaxisIDp, zaxisIDh = -1, nzaxis;
int ngrids, gridID, zaxisID;
int nhlev = 0, nhlevf = 0, nhlevh = 0, nlevel;
int gridID, zaxisID;
int nhlev = 0, nhlevf = 0, nhlevh = 0, nlevel = 0;
int nvct;
int geopID = -1, tempID = -1, psID = -1, lnpsID = -1, pvarID = -1;
int psID = -1, lnpsID = -1, pvarID = -1;
int code, param;
char paramstr[32];
char varname[CDI_MAX_NAME];
......@@ -246,15 +246,11 @@ void *Pressure(void *argument)
if ( tableNum == 2 )
{
mode = WMO_MODE;
geop_code = 6;
temp_code = 11;
ps_code = 1;
}
else if ( tableNum == 128 )
{
mode = ECHAM_MODE;
geop_code = 129;
temp_code = 130;
ps_code = 134;
lsp_code = 152;
}
......@@ -264,8 +260,6 @@ void *Pressure(void *argument)
else
{
mode = ECHAM_MODE;
geop_code = 129;
temp_code = 130;
ps_code = 134;
lsp_code = 152;
}
......@@ -290,17 +284,13 @@ void *Pressure(void *argument)
if ( mode == ECHAM_MODE )
{
if ( code == geop_code && nlevel == 1 ) geopID = varID;
else if ( code == temp_code && nlevel == nhlev ) tempID = varID;
else if ( code == ps_code && nlevel == 1 ) psID = varID;
if ( code == ps_code && nlevel == 1 ) psID = varID;
else if ( code == lsp_code && nlevel == 1 ) lnpsID = varID;
/* else if ( code == 156 ) gheightID = varID; */
}
else if ( mode == WMO_MODE )
{
if ( code == geop_code && nlevel == 1 ) geopID = varID;
else if ( code == temp_code && nlevel == nhlev ) tempID = varID;
else if ( code == ps_code && nlevel == 1 ) psID = varID;
if ( code == ps_code && nlevel == 1 ) psID = varID;
}
}
......@@ -331,7 +321,7 @@ void *Pressure(void *argument)
vlistID2 = vlistCreate();
varID = vlistDefVar(vlistID2, gridID, zaxisIDp, TSTEP_INSTANT);
vlistDefVarCode(vlistID2, varID, 1);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(1, 255, 255));
vlistDefVarName(vlistID2, varID, "pressure");
vlistDefVarStdname(vlistID2, varID, "air_pressure");
vlistDefVarUnits(vlistID2, varID, "Pa");
......@@ -394,7 +384,7 @@ void *Pressure(void *argument)
pout = deltap;
}
else
else if ( operatorID == PRESSURE_HL )
{
nlevel = nhlevh;
pout = half_press;
......