Commit 1356c277 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Exprf: added support for vertstat

parent a6145424
......@@ -96,7 +96,7 @@ void genLayerBounds(int nlev, double *levels, double *lbounds, double *ubounds)
}
}
static
int getLayerThickness(int genbounds, int index, int zaxisID, int nlev, double *thickness, double *weights)
{
int status = 0;
......
......@@ -35,6 +35,7 @@
static const char *ExIn[] = {"expr", "init"};
static const char *tmpvnm = "_tmp_";
int pointID = -1;
int surfaceID = -1;
#define COMPLT(x,y) ((x) < (y) ? 1 : 0)
#define COMPGT(x,y) ((x) > (y) ? 1 : 0)
......@@ -95,7 +96,7 @@ static func_t fun_sym_tbl[] =
{0, 0, "atanh", (void (*)()) atanh},
{0, 0, "gamma", (void (*)()) tgamma},
// cdo functions (Reduce grid to point)
// cdo field functions (Reduce grid to point)
{1, 0, "fldmin", (void (*)()) fldmin},
{1, 0, "fldmax", (void (*)()) fldmax},
{1, 0, "fldsum", (void (*)()) fldsum},
......@@ -105,6 +106,17 @@ static func_t fun_sym_tbl[] =
{1, 1, "fldstd1", (void (*)()) fldstd1},
{1, 1, "fldvar", (void (*)()) fldvar},
{1, 1, "fldvar1", (void (*)()) fldvar1},
// cdo field functions (Reduce grid to point)
{2, 0, "vertmin", (void (*)()) fldmin},
{2, 0, "vertmax", (void (*)()) fldmax},
{2, 0, "vertsum", (void (*)()) fldsum},
{2, 1, "vertmean", (void (*)()) fldmean},
{2, 1, "vertavg", (void (*)()) fldavg},
{2, 1, "vertstd", (void (*)()) fldstd},
{2, 1, "vertstd1", (void (*)()) fldstd1},
{2, 1, "vertvar", (void (*)()) fldvar},
{2, 1, "vertvar1", (void (*)()) fldvar1},
};
static int NumFunc = sizeof(fun_sym_tbl) / sizeof(fun_sym_tbl[0]);
......@@ -699,6 +711,9 @@ nodeType *expr(int init, int oper, nodeType *p1, nodeType *p2)
static
nodeType *ex_fun_con(int funcID, nodeType *p1)
{
int functype = fun_sym_tbl[funcID].type;
if ( functype != 0 ) cdoAbort("Function %s not available for constant values!", fun_sym_tbl[funcID].name);
nodeType *p = (nodeType*) Calloc(1, sizeof(nodeType));
p->type = typeCon;
......@@ -734,11 +749,16 @@ nodeType *ex_fun_var(int init, int funcID, nodeType *p1)
p->param.gridID = pointID;
p->param.ngp = 1;
}
else if ( functype == 2 )
{
p->param.zaxisID = surfaceID;
p->param.nlev = 1;
}
if ( ! init )
{
p->param.data = (double*) Malloc(p->param.ngp*nlev*sizeof(double));
double *restrict pdata = p->param.data;
p->param.data = (double*) Malloc(p->param.ngp*p->param.nlev*sizeof(double));
double *restrict pdata = p->param.data;
double *restrict p1data = p1->param.data;
if ( functype == 0 )
......@@ -779,9 +799,25 @@ nodeType *ex_fun_var(int init, int funcID, nodeType *p1)
}
if ( weights ) Free(weights);
}
else if ( functype == 2 )
{
field_t field;
double *weights = NULL;
if ( funcflag == 1 ) weights = vert_weights(p1->param.zaxisID, nlev);
double *array = (double*) Malloc(nlev*sizeof(double));
double (*exprfunc)(field_t) = (double (*)(field_t)) fun_sym_tbl[funcID].func;
for ( size_t i = 0; i < ngp; i++ )
{
for ( size_t k = 0; k < nlev; k++ ) array[k] = p1data[k*ngp+i];
fld_field_init(&field, nmiss, missval, nlev, array, weights);
pdata[i] = exprfunc(field);
}
if ( array ) Free(array);
if ( weights ) Free(weights);
}
nmiss = 0;
for ( size_t i = 0; i < p->param.ngp*nlev; i++ )
for ( size_t i = 0; i < p->param.ngp*p->param.nlev; i++ )
if ( DBL_IS_EQUAL(pdata[i], missval) ) nmiss++;
p->param.nmiss = nmiss;
......@@ -1075,6 +1111,7 @@ int param_search_name(int nparam, paramType *params, const char *name)
nodeType *expr_run(nodeType *p, parse_param_t *parse_arg)
{
pointID = parse_arg->pointID;
surfaceID = parse_arg->surfaceID;
int init = parse_arg->init;
paramType *params = parse_arg->params;
//paramType *param2 = &parse_arg->param2;
......
......@@ -50,48 +50,27 @@ double *fld_weights(int gridID, size_t ngp)
return weights;
}
/*
double fun_fldmin(size_t nmiss, double missval, size_t ngp, double *array, double *weights)
{
field_t field;
fld_field_init(&field, nmiss, missval, ngp, array, weights);
return fldmin(field);
}
int getLayerThickness(int genbounds, int index, int zaxisID, int nlev, double *thickness, double *weights);
double fun_fldmax(size_t nmiss, double missval, size_t ngp, double *array, double *weights)
double *vert_weights(int zaxisID, size_t nlev)
{
field_t field;
fld_field_init(&field, nmiss, missval, ngp, array, weights);
return fldmax(field);
}
double fun_fldsum(size_t nmiss, double missval, size_t ngp, double *array, double *weights)
{
field_t field;
fld_field_init(&field, nmiss, missval, ngp, array, weights);
return fldsum(field);
}
double fun_fldmean(size_t nmiss, double missval, size_t ngp, double *array, double *weights)
{
field_t field;
fld_field_init(&field, nmiss, missval, ngp, array, weights);
return fldmean(field);
}
static bool lwarn = true;
double *weights = (double*) Malloc(nlev*sizeof(double));
double *thickness = (double*) Malloc(nlev*sizeof(double));
for ( size_t i = 0; i < nlev; ++i ) weights[i] = 1;
if ( nlev > 1 )
{
int wstatus = getLayerThickness(0, 0, zaxisID, nlev, thickness, weights);
if ( wstatus != 0 && lwarn && nlev > 1 )
{
lwarn = false;
cdoWarning("Layer bounds not available, using constant vertical weights!");
}
}
double fun_fldavg(size_t nmiss, double missval, size_t ngp, double *array, double *weights)
{
field_t field;
fld_field_init(&field, nmiss, missval, ngp, array, weights);
Free(thickness);
return fldavg(field);
return weights;
}
*/
......@@ -18,13 +18,8 @@
#ifndef _EXPR_FUN_H
#define _EXPR_FUN_H
double *fld_weights(int gridID, size_t ngp);
void fld_field_init(field_t *field, size_t nmiss, double missval, size_t ngp, double *array, double *w);
/*
double fun_fldmin(size_t nmiss, double missval, size_t ngp, double *array, double *weights);
double fun_fldmax(size_t nmiss, double missval, size_t ngp, double *array, double *weights);
double fun_fldsum(size_t nmiss, double missval, size_t ngp, double *array, double *weights);
double fun_fldmean(size_t nmiss, double missval, size_t ngp, double *array, double *weights);
double fun_fldavg(size_t nmiss, double missval, size_t ngp, double *array, double *weights);
*/
double *fld_weights(int gridID, size_t ngp);
double *vert_weights(int zaxisID, size_t nlev);
#endif
......@@ -305,6 +305,9 @@ static
void prevarsum(const double *restrict array, const double *restrict w, size_t len, int nmiss,
double missval, double *restrict rsum, double *restrict rsumw, double *restrict rsumq, double *restrict rsumwq)
{
assert(array!=NULL);
assert(w!=NULL);
*rsum = *rsumw = 0;
*rsumq = *rsumwq = 0;
......
Supports Markdown
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