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

expr:zonSTAT: wrong result (bug fix).

parent 2c0b6ade
......@@ -5,6 +5,7 @@
2019-03-25 Uwe Schulzweida
* expr:zonSTAT: wrong result (bug fix)
* expr::vertmean: fix wrong warning message about layer bounds
2019-03-25 Uwe Schulzweida
......
......@@ -374,11 +374,11 @@ genZonalID(int vlistID)
{
int zonalID = -1;
int ngrids = vlistNgrids(vlistID);
const int ngrids = vlistNgrids(vlistID);
for (int index = 0; index < ngrids; index++)
{
int gridID = vlistGrid(vlistID, index);
int gridtype = gridInqType(gridID);
const int gridID = vlistGrid(vlistID, index);
const int gridtype = gridInqType(gridID);
if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || gridtype == GRID_GENERIC)
if (gridInqXsize(gridID) > 1 && gridInqYsize(gridID) >= 1)
{
......@@ -409,9 +409,9 @@ Expr(void *process)
cdoOperatorAdd("aexprf", 0, 0, "expr script filename");
// clang-format on
int operatorID = cdoOperatorID();
bool replacesVariables = cdoOperatorF1(operatorID);
bool readsCommandLine = cdoOperatorF2(operatorID);
const int operatorID = cdoOperatorID();
const bool replacesVariables = cdoOperatorF1(operatorID);
const bool readsCommandLine = cdoOperatorF2(operatorID);
operatorInputArg(cdoOperatorEnter(operatorID));
......@@ -420,17 +420,17 @@ Expr(void *process)
char *exprs = readsCommandLine ? exprsFromArgument(exprArgc, exprArgv) : exprsFromFile(exprArgc, exprArgv);
int streamID1 = cdoStreamOpenRead(0);
int vlistID1 = cdoStreamInqVlist(streamID1);
const int streamID1 = cdoStreamOpenRead(0);
const int vlistID1 = cdoStreamInqVlist(streamID1);
exprs = exprsExpand(exprs, vlistID1);
if (Options::cdoVerbose) cdoPrint(exprs);
int nvars1 = vlistNvars(vlistID1);
const int nvars1 = vlistNvars(vlistID1);
int pointID = gridCreate(GRID_GENERIC, 1);
int zonalID = genZonalID(vlistID1);
int surfaceID = getSurfaceID(vlistID1);
const int pointID = gridCreate(GRID_GENERIC, 1);
const int zonalID = genZonalID(vlistID1);
const int surfaceID = getSurfaceID(vlistID1);
paramType *params = params_new(vlistID1);
......@@ -439,7 +439,7 @@ Expr(void *process)
// Set all input variables to 'needed' if replacing is switched off
for (int varID = 0; varID < nvars1; varID++) parse_arg.needed[varID] = !replacesVariables;
int vartsID = params_add_ts(&parse_arg);
const int vartsID = params_add_ts(&parse_arg);
parse_arg.tsID = vartsID;
params_add_coordinates(vlistID1, &parse_arg);
......@@ -535,9 +535,9 @@ Expr(void *process)
for (int varID = parse_arg.nvars1; varID < parse_arg.nparams; varID++)
{
size_t ngp = params[varID].ngp;
size_t nlev = params[varID].nlev;
size_t nItems = cdo::max((size_t) 4, ngp * nlev);
const size_t ngp = params[varID].ngp;
const size_t nlev = params[varID].nlev;
const size_t nItems = cdo::max((size_t) 4, ngp * nlev);
params[varID].data = (double *) Malloc(nItems * sizeof(double));
}
......@@ -545,16 +545,16 @@ Expr(void *process)
{
if (parse_arg.coords[i].needed)
{
int coord = parse_arg.coords[i].coord;
const int coord = parse_arg.coords[i].coord;
if (coord == 'x' || coord == 'y' || coord == 'a' || coord == 'w')
{
int gridID = parse_arg.coords[i].cdiID;
size_t ngp = parse_arg.coords[i].size;
const size_t ngp = parse_arg.coords[i].size;
double *data = (double *) Malloc(ngp * sizeof(double));
parse_arg.coords[i].data = data;
if (coord == 'x' || coord == 'y')
{
int gridtype = gridInqType(gridID);
const int gridtype = gridInqType(gridID);
if (gridtype == GRID_GME) gridID = gridToUnstructured(gridID, 0);
if (gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR) gridID = gridToCurvilinear(gridID, 0);
......@@ -594,8 +594,8 @@ Expr(void *process)
}
else if (coord == 'z')
{
int zaxisID = parse_arg.coords[i].cdiID;
size_t nlev = parse_arg.coords[i].size;
const int zaxisID = parse_arg.coords[i].cdiID;
const size_t nlev = parse_arg.coords[i].size;
double *data = (double *) Malloc(nlev * sizeof(double));
parse_arg.coords[i].data = data;
cdoZaxisInqLevels(zaxisID, data);
......@@ -607,16 +607,16 @@ Expr(void *process)
for (int varID = parse_arg.nvars1; varID < parse_arg.nparams; varID++)
{
int coord = params[varID].coord;
const int coord = params[varID].coord;
if (coord)
{
char *varname = strdup(params[varID].name);
varname[strlen(varname) - 2] = 0;
if (coord == 'x' || coord == 'y' || coord == 'a' || coord == 'w')
{
int coordID = params_get_coordID(&parse_arg, coord, params[varID].gridID);
int gridID = parse_arg.coords[coordID].cdiID;
size_t ngp = parse_arg.coords[coordID].size;
const int coordID = params_get_coordID(&parse_arg, coord, params[varID].gridID);
const int gridID = parse_arg.coords[coordID].cdiID;
const size_t ngp = parse_arg.coords[coordID].size;
double *data = parse_arg.coords[coordID].data;
assert(gridID == params[varID].gridID);
assert(data != nullptr);
......@@ -625,9 +625,9 @@ Expr(void *process)
}
else if (coord == 'z')
{
int coordID = params_get_coordID(&parse_arg, coord, params[varID].zaxisID);
int zaxisID = parse_arg.coords[coordID].cdiID;
size_t nlev = parse_arg.coords[coordID].size;
const int coordID = params_get_coordID(&parse_arg, coord, params[varID].zaxisID);
const int zaxisID = parse_arg.coords[coordID].cdiID;
const size_t nlev = parse_arg.coords[coordID].size;
double *data = parse_arg.coords[coordID].data;
assert(zaxisID == params[varID].zaxisID);
assert(data != nullptr);
......@@ -643,23 +643,23 @@ Expr(void *process)
if (Options::cdoVerbose) vlistPrint(vlistID2);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int streamID2 = cdoStreamOpenWrite(1);
const int streamID2 = cdoStreamOpenWrite(1);
cdoDefVlist(streamID2, vlistID2);
int64_t vdate0 = 0;
int vtime0 = 0;
int calendar = taxisInqCalendar(taxisID1);
const int calendar = taxisInqCalendar(taxisID1);
int nrecs;
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
int64_t vdate = taxisInqVdate(taxisID1);
int vtime = taxisInqVtime(taxisID1);
const int64_t vdate = taxisInqVdate(taxisID1);
const int vtime = taxisInqVtime(taxisID1);
double jdelta = 0;
......@@ -703,7 +703,7 @@ Expr(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
if (parse_arg.needed[varID])
{
size_t offset = params[varID].ngp * levelID;
const size_t offset = params[varID].ngp * levelID;
double *vardata = params[varID].data + offset;
size_t nmiss;
cdoReadRecord(streamID1, vardata, &nmiss);
......@@ -713,10 +713,10 @@ Expr(void *process)
for (int varID = 0; varID < nvars2; varID++)
{
int pidx = varIDmap[varID];
const int pidx = varIDmap[varID];
if (pidx < nvars1) continue;
size_t ngp = params[pidx].ngp;
size_t nlev = params[pidx].nlev;
const size_t ngp = params[pidx].ngp;
const size_t nlev = params[pidx].nlev;
params[pidx].nmiss = 0;
arrayFill(ngp * nlev, params[pidx].data, 0.0);
......@@ -727,19 +727,19 @@ Expr(void *process)
for (int varID = 0; varID < nvars2; varID++)
{
int pidx = varIDmap[varID];
const int pidx = varIDmap[varID];
if (tsID > 0 && params[pidx].steptype == TIME_CONSTANT) continue;
double missval = vlistInqVarMissval(vlistID2, varID);
const double missval = vlistInqVarMissval(vlistID2, varID);
size_t ngp = params[pidx].ngp;
int nlev = (int) params[pidx].nlev;
const size_t ngp = params[pidx].ngp;
const int nlev = (int) params[pidx].nlev;
for (int levelID = 0; levelID < nlev; levelID++)
{
size_t offset = ngp * levelID;
const size_t offset = ngp * levelID;
double *vardata = params[pidx].data + offset;
size_t nmiss = arrayNumMV(ngp, vardata, missval);
const size_t nmiss = arrayNumMV(ngp, vardata, missval);
cdoDefRecord(streamID2, varID, levelID);
cdoWriteRecord(streamID2, vardata, nmiss);
}
......
......@@ -904,13 +904,16 @@ ex_fun_var(const int init, const int funcID, nodeType *p1)
else if (functype == FT_ZON)
{
Field field1, field2;
field1.resize(ngp);
field2.resize(nlat);
void (*exprfunc)(const Field &, Field &) = (void (*)(const Field &, Field &)) fun_sym_tbl[funcID].func;
for (size_t k = 0; k < nlev; k++)
{
fld_field_init(field1, nmiss, missval, ngp, p1data + k * ngp, nullptr);
fld_field_init(field1, nmiss, missval, ngp, &p1data[k * ngp], nullptr);
field1.grid = gridID;
fld_field_init(field2, nmiss, missval, nlat, &pdata[k], nullptr);
fld_field_init(field2, nmiss, missval, nlat, &pdata[k * nlat], nullptr);
exprfunc(field1, field2);
arrayCopy(nlat, field2.vec.data(), &pdata[k * nlat]);
}
}
else if (functype == FT_VERT)
......
......@@ -24,9 +24,8 @@ fld_field_init(Field &field, size_t nmiss, double missval, size_t ngp, double *a
field.gridsize = ngp;
field.nmiss = nmiss;
field.missval = missval;
field.ptr = array;
if (array) for (size_t i = 0; i < ngp; ++i) field.vec[i] = array[i];
if (w) for (size_t i = 0; i < ngp; ++i) field.weightv[i] = w[i];
if (array != nullptr) for (size_t i = 0; i < ngp; ++i) field.vec[i] = array[i];
if (w != nullptr) for (size_t i = 0; i < ngp; ++i) field.weightv[i] = w[i];
}
/*
double *
......
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