Commit 35eedf40 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

expr: added support for calculations with var1[n] and var2[1]

parent 1f102f51
......@@ -3,6 +3,10 @@
* using CDI library version 1.7.0
* Version 1.7.0 released
2015-08-19 Uwe Schulzweida
* expr: added support for calculations with var1[n] and var2[1]
2015-08-18 Uwe Schulzweida
* ap2pl: use upper level of air_pressure if surface pressure not found
......
......@@ -109,33 +109,11 @@ nodeType *expr_con_con(int oper, nodeType *p1, nodeType *p2)
}
static
nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
void oper_expr_con_var(int oper, int nmiss, long n, double missval1, double missval2,
double *restrict odat, double cval, const double *restrict idat)
{
int gridID = p2->gridID;
int zaxisID = p2->zaxisID;
int nmiss = p2->nmiss;
double missval1 = p2->missval;
double missval2 = p2->missval;
int ngp = gridInqSize(gridID);
int nlev = zaxisInqSize(zaxisID);
long n = ngp*nlev;
long i;
nodeType *p = (nodeType*) malloc(sizeof(nodeType));
p->type = typeVar;
p->tmpvar = 1;
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval1;
p->data = (double*) malloc(n*sizeof(double));
double *restrict odat = p->data;
const double *restrict idat = p2->data;
double cval = p1->u.con.value;
switch ( oper )
{
case '+':
......@@ -196,46 +174,14 @@ nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
cdoAbort("%s: operator %c unsupported!", __func__, oper);
break;
}
nmiss = 0;
for ( i = 0; i < n; i++ )
if ( DBL_IS_EQUAL(p->data[i], missval1) ) nmiss++;
p->nmiss = nmiss;
if ( p2->tmpvar ) free(p2->data);
return (p);
}
static
nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
void oper_expr_var_con(int oper, int nmiss, long n, double missval1, double missval2,
double *restrict odat, const double *restrict idat, double cval)
{
int gridID = p1->gridID;
int zaxisID = p1->zaxisID;
int nmiss = p1->nmiss;
double missval1 = p1->missval;
double missval2 = p1->missval;
int ngp = gridInqSize(gridID);
int nlev = zaxisInqSize(zaxisID);
long n = ngp*nlev;
long i;
nodeType *p = (nodeType*) malloc(sizeof(nodeType));
p->type = typeVar;
p->tmpvar = 1;
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval1;
p->data = (double*) malloc(n*sizeof(double));
double *restrict odat = p->data;
const double *restrict idat = p1->data;
double cval = p2->u.con.value;
switch ( oper )
{
case '+':
......@@ -298,6 +244,155 @@ nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
cdoAbort("%s: operator %c unsupported!", __func__, oper);
break;
}
}
static
void oper_expr_var_var(int oper, int nmiss, long ngp, double missval1, double missval2,
double *restrict odat, const double *restrict idat1, const double *restrict idat2)
{
long i;
switch ( oper )
{
case '+':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = ADD(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] + idat2[i];
break;
case '-':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = SUB(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] - idat2[i];
break;
case '*':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MUL(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] * idat2[i];
break;
case '/':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = DIV(idat1[i], idat2[i]);
else
{
for ( i = 0; i < ngp; ++i )
{
if ( IS_EQUAL(idat2[i], 0.) ) odat[i] = missval1;
else odat[i] = idat1[i] / idat2[i];
}
}
break;
case '^':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = POW(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = pow(idat1[i], idat2[i]);
break;
case '<':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLT(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLT(idat1[i], idat2[i]);
break;
case '>':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPGT(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPGT(idat1[i], idat2[i]);
break;
case LE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLE(idat1[i], idat2[i]);
break;
case GE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPGE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPGE(idat1[i], idat2[i]);
break;
case NE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPNE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPNE(idat1[i], idat2[i]);
break;
case EQ:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPEQ(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPEQ(idat1[i], idat2[i]);
break;
case LEG:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLEG(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLEG(idat1[i], idat2[i]);
break;
case AND:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPAND(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPAND(idat1[i], idat2[i]);
break;
case OR:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPOR(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPOR(idat1[i], idat2[i]);
break;
default:
cdoAbort("%s: operator %d (%c) unsupported!", __func__, (int)oper, oper);
break;
}
}
static
nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
{
int gridID = p2->gridID;
int zaxisID = p2->zaxisID;
int nmiss = p2->nmiss;
double missval1 = p2->missval;
double missval2 = p2->missval;
int ngp = gridInqSize(gridID);
int nlev = zaxisInqSize(zaxisID);
long n = ngp*nlev;
long i;
nodeType *p = (nodeType*) malloc(sizeof(nodeType));
p->type = typeVar;
p->tmpvar = 1;
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval1;
p->data = (double*) malloc(n*sizeof(double));
double *restrict odat = p->data;
const double *restrict idat = p2->data;
double cval = p1->u.con.value;
oper_expr_con_var(oper, nmiss, n, missval1, missval2, odat, cval, idat);
nmiss = 0;
for ( i = 0; i < n; i++ )
if ( DBL_IS_EQUAL(p->data[i], missval1) ) nmiss++;
p->nmiss = nmiss;
if ( p2->tmpvar ) free(p2->data);
return (p);
}
static
nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
{
int gridID = p1->gridID;
int zaxisID = p1->zaxisID;
int nmiss = p1->nmiss;
double missval1 = p1->missval;
double missval2 = p1->missval;
int ngp = gridInqSize(gridID);
int nlev = zaxisInqSize(zaxisID);
long n = ngp*nlev;
long i;
nodeType *p = (nodeType*) malloc(sizeof(nodeType));
p->type = typeVar;
p->tmpvar = 1;
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval1;
p->data = (double*) malloc(n*sizeof(double));
double *restrict odat = p->data;
const double *restrict idat = p1->data;
double cval = p2->u.con.value;
oper_expr_var_con(oper, nmiss, n, missval1, missval2, odat, idat, cval);
nmiss = 0;
for ( i = 0; i < n; i++ )
......@@ -326,7 +421,7 @@ nodeType *expr_var_var(int oper, nodeType *p1, nodeType *p2)
long ngp1 = gridInqSize(p1->gridID);
long ngp2 = gridInqSize(p2->gridID);
if ( ngp1 != ngp2 ) cdoAbort("Number of grid points differ. ngp1 = %ld, ngp2 = %ld", ngp1, ngp2);
if ( ngp1 != ngp2 && ngp2 != 1 ) cdoAbort("Number of grid points differ. ngp1 = %ld, ngp2 = %ld", ngp1, ngp2);
long ngp = ngp1;
......@@ -370,85 +465,20 @@ nodeType *expr_var_var(int oper, nodeType *p1, nodeType *p2)
loff = k*ngp;
if ( nlev1 == 1 ) loff1 = 0;
else loff1 = k*ngp;
else loff1 = k*ngp1;
if ( nlev2 == 1 ) loff2 = 0;
else loff2 = k*ngp;
else loff2 = k*ngp2;
const double *restrict idat1 = p1->data+loff1;
const double *restrict idat2 = p2->data+loff2;
double *restrict odat = p->data+loff;
int nmiss = nmiss1 > 0 || nmiss2 > 0;
switch ( oper )
{
case '+':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = ADD(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] + idat2[i];
break;
case '-':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = SUB(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] - idat2[i];
break;
case '*':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MUL(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = idat1[i] * idat2[i];
break;
case '/':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = DIV(idat1[i], idat2[i]);
else
{
for ( i = 0; i < ngp; ++i )
{
if ( IS_EQUAL(idat2[i], 0.) ) odat[i] = missval1;
else odat[i] = idat1[i] / idat2[i];
}
}
break;
case '^':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = POW(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = pow(idat1[i], idat2[i]);
break;
case '<':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLT(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLT(idat1[i], idat2[i]);
break;
case '>':
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPGT(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPGT(idat1[i], idat2[i]);
break;
case LE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLE(idat1[i], idat2[i]);
break;
case GE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPGE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPGE(idat1[i], idat2[i]);
break;
case NE:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPNE(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPNE(idat1[i], idat2[i]);
break;
case EQ:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPEQ(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPEQ(idat1[i], idat2[i]);
break;
case LEG:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPLEG(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPLEG(idat1[i], idat2[i]);
break;
case AND:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPAND(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPAND(idat1[i], idat2[i]);
break;
case OR:
if ( nmiss ) for ( i=0; i<ngp; ++i ) odat[i] = MVCOMPOR(idat1[i], idat2[i]);
else for ( i=0; i<ngp; ++i ) odat[i] = COMPOR(idat1[i], idat2[i]);
break;
default:
cdoAbort("%s: operator %d (%c) unsupported!", __func__, (int)oper, oper);
break;
}
if ( ngp2 == 1 )
oper_expr_var_con(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2[0]);
else
oper_expr_var_var(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2);
}
nmiss = 0;
......
......@@ -36,10 +36,9 @@ FILE *old, *new;
static
long int_index(long n, double *x1, double x2)
{
long klo, khi, k;
klo = 0;
khi = n-1;
long k;
long klo = 0;
long khi = n-1;
while ( khi-klo > 1 )
{
......@@ -56,24 +55,21 @@ long int_index(long n, double *x1, double x2)
*/
/* printf("%d %d %g %g %g\n", klo, khi, x1[klo], x1[klo+1], x2); */
return (klo);
return klo;
}
static
double esat(double temperature)
{
double zes;
double es;
double tc;
tc = temperature-273.16;
double tc = temperature-273.16;
if ( tc < 0.0 )
zes = 21.8745584*tc / (temperature-7.66);
else
zes = 17.2693882*tc / (temperature-35.86);
es = 610.78*exp(zes);
double es = 610.78*exp(zes);
return es;
}
......@@ -103,19 +99,17 @@ void hetaeta_sc(int ltq, int lpsmod, long ij, long ngp, long nlev1, long nlev2,
long k, iv, ijk, ijk1, ijk2;
long jlev = 0, jlevr, jnop;
long klo;
long nlev1p1, nlev2p1;
long jjblt;
double zq1, zt1;
double zdff, zdffl, ztv, zb, zbb, zc, zps;
double zsump, zsumpp, zsumt, zsumtp;
double dfi, fiadj = 0, dteta = 0;
double pbl_lim, pbl_lim_need;
double rair_d_cpair;
rair_d_cpair = rair/cpair;
double rair_d_cpair = rair/cpair;
nlev1p1 = nlev1+1;
nlev2p1 = nlev2+1;
long nlev1p1 = nlev1+1;
long nlev2p1 = nlev2+1;
/******* initialise atmospheric fields in old system */
......@@ -464,8 +458,7 @@ void hetaeta(int ltq, int ngp, const int *imiss,
const double *restrict fis2, double *restrict ps2,
double *restrict t2, double *restrict q2,
int nvars, double *restrict *restrict vars1, double *restrict *restrict vars2,
double *restrict tscor, double *restrict pscor,
double *restrict secor)
double *restrict tscor, double *restrict pscor, double *restrict secor)
{
double epsm1i;
long jblt;
......@@ -473,29 +466,18 @@ void hetaeta(int ltq, int ngp, const int *imiss,
long jlev = 0;
double *zt2 = NULL, *zq2 = NULL;
double /* *etah1,*/ *ph1, *lnph1, *fi1;
double *af1, *bf1, *etaf1, *pf1, *lnpf1;
double *pf1, *lnpf1;
double *tv1, *theta1, *rh1;
double *etah2, *ph2, *lnph2, *fi2;
double *af2, *bf2, *etaf2, *pf2/*, *lnpf2*/;
double *ph2, *lnph2, *fi2;
double *pf2/*, *lnpf2*/;
double *zvar;
double *rh_pbl = NULL;
double *theta_pbl = NULL;
double *rh2;
double *w1, *w2;
double *wgt;
long *idx;
long *jl1, *jl2;
long nlev1p1;
long nlev2p1;
int lpsmod = 1;
#if defined(_OPENMP)
double **zt2_2, **zq2_2;
double **ph1_2, **lnph1_2, **fi1_2, **pf1_2, **lnpf1_2, **tv1_2, **theta1_2, **rh1_2, **zvar_2;
double **ph2_2, **lnph2_2, **fi2_2, **pf2_2;
double **rh_pbl_2, **theta_pbl_2, **rh2_2;
double **wgt_2;
double ***vars_pbl_2;
long **idx_2;
long i;
extern int ompNumThreads;
double *vars_pbl[MAX_VARS];
......@@ -509,42 +491,38 @@ void hetaeta(int ltq, int ngp, const int *imiss,
new = fopen("new.dat","w");
#endif
nlev1p1 = nlev1+1;
nlev2p1 = nlev2+1;
long nlev1p1 = nlev1+1;
long nlev2p1 = nlev2+1;
#if defined(_OPENMP)
ph1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
lnph1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
fi1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
pf1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
lnpf1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
tv1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
theta1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
rh1_2 = (double **) malloc(ompNumThreads*sizeof(double *));
zvar_2 = (double **) malloc(ompNumThreads*sizeof(double *));
ph2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
lnph2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
fi2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
pf2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
rh2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
wgt_2 = (double **) malloc(ompNumThreads*sizeof(double *));
idx_2 = (long **) malloc(ompNumThreads*sizeof(long *));
if ( ltq )
{
zt2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
zq2_2 = (double **) malloc(ompNumThreads*sizeof(double *));
rh_pbl_2 = (double **) malloc(ompNumThreads*sizeof(double *));
theta_pbl_2 = (double **) malloc(ompNumThreads*sizeof(double *));
}
if ( nvars > 0 )
{
vars_pbl_2 = (double ***) malloc(ompNumThreads*sizeof(double **));
}
double *ph1_2[ompNumThreads];
double *lnph1_2[ompNumThreads];
double *fi1_2[ompNumThreads];
double *pf1_2[ompNumThreads];
double *lnpf1_2[ompNumThreads];
double *tv1_2[ompNumThreads];
double *theta1_2[ompNumThreads];
double *rh1_2[ompNumThreads];
double *zvar_2[ompNumThreads];
double *ph2_2[ompNumThreads];
double *lnph2_2[ompNumThreads];
double *fi2_2[ompNumThreads];
double *pf2_2[ompNumThreads];
double *rh2_2[ompNumThreads];
double *wgt_2[ompNumThreads];
long *idx_2[ompNumThreads];
// if ( ltq )
double *zt2_2[ompNumThreads];
double *zq2_2[ompNumThreads];
double *rh_pbl_2[ompNumThreads];
double *theta_pbl_2[ompNumThreads];
// if ( nvars > 0 )
double **vars_pbl_2[ompNumThreads];
for ( i = 0; i < ompNumThreads; i++ )
{
......@@ -627,21 +605,20 @@ void hetaeta(int ltq, int ngp, const int *imiss,
}
#endif
af1 = (double*) malloc(nlev1*sizeof(double));
bf1 = (double*) malloc(nlev1*sizeof(double));
etaf1 = (double*) malloc(nlev1*sizeof(double));
etah2 = (double*) malloc(nlev2p1*sizeof(double));
double af1[nlev1];
double bf1[nlev1];
double etaf1[nlev1];
af2 = (double*) malloc(nlev2*sizeof(double));
bf2 = (double*) malloc(nlev2*sizeof(double));
etaf2 = (double*) malloc(nlev2*sizeof(double));
double etah2[nlev2p1];
w1 = (double*) malloc(nlev2*sizeof(double));
w2 = (double*) malloc(nlev2*sizeof(double));
jl1 = (long*) malloc(nlev2*sizeof(long));
jl2 = (long*) malloc(nlev2*sizeof(long));
double af2[nlev2];
double bf2[nlev2];
double etaf2[nlev2];
double w1[nlev2];
double w2[nlev2];
long jl1[nlev2];
long jl2[nlev2];
/******* set coordinate system ETA's, A's, B's
calculate half and full level ETA
......@@ -831,39 +808,6 @@ void hetaeta(int ltq, int ngp, const int *imiss,
free(vars_pbl_2[i]);
}
}
free(ph1_2);
free(lnph1_2);
free(fi1_2);
free(pf1_2);
free(lnpf1_2);
free(tv1_2);
free(theta1_2);
free(rh1_2);
free(zvar_2);
free(ph2_2);
free(lnph2_2);
free(fi2_2);
free(pf2_2);
free(rh2_2);
free(wgt_2);
free(idx_2);
if ( ltq )
{
free(zt2_2);
free(zq2_2);
free(rh_pbl_2);
free(theta_pbl_2);
}
if ( nvars > 0 )
{
free(vars_pbl_2);
}
#else
/* free(etah1); */
free(ph1);
......@@ -904,27 +848,17 @@ void hetaeta(int ltq, int ngp, const int *imiss,
}
#endif