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

vertvar, vertstd: changed to weighted var/std if layer bounds are available

parent 82417228
......@@ -3,6 +3,10 @@
* using CDI library version 1.7.0
* Version 1.7.0 released
2015-06-01 Uwe Schulzweida
* vertvar, vertstd: changed to weighted var/std if layer bounds are available
2015-05-29 Uwe Schulzweida
* yseaspctl: check of verification date failed [Bug #5810]
......@@ -14,7 +18,7 @@
2015-05-26 Uwe Schulzweida
* New operator: vertstd1 - Vertical standard deviation [Divisor is (n-1)]
* New operator: vertvar1 - Vartical variance [Divisor is (n-1)]
* New operator: vertvar1 - Vertical variance [Divisor is (n-1)]
* New operator: gridboxstd1 - Gridbox standard deviation [Divisor is (n-1)]
* New operator: gridboxvar1 - Gridbox variance [Divisor is (n-1)]
* New operator: merstd1 - Meridional standard deviation [Divisor is (n-1)]
......
......@@ -4,10 +4,14 @@ CDO NEWS
Version 1.7.0 (28 October 2015):
New operators:
* vertstd1: Vertical standard deviation [Divisor is (n-1)]
* vertvar1: Vertical variance [Divisor is (n-1)]
* seasvar1: Seasonal variance [Divisor is (n-1)]
* seasstd1: Seasonal standard deviation [Divisor is (n-1)]
* yseasvar1: Multi-year seasonally variance [Divisor is (n-1)]
* yseasstd1: Multi-year seasonally standard deviation [Divisor is (n-1)]
Changed operators:
* vertvar, vertstd: changed to weighted var/std if layer bounds are available
Fixed bugs:
* env. CDO_TIMESTAT_DATE does not work [Bug #5758]
* splityear*: support for constant fields is missing [Bug #5759]
......
......@@ -179,10 +179,10 @@ void *Vertstat(void *argument)
int VERTINT = cdoOperatorAdd("vertint", func_sum, 1, NULL);
cdoOperatorAdd("vertmean", func_mean, 1, NULL);
cdoOperatorAdd("vertavg", func_avg, 1, NULL);
cdoOperatorAdd("vertvar", func_var, 0, NULL);
cdoOperatorAdd("vertvar1", func_var1, 0, NULL);
cdoOperatorAdd("vertstd", func_std, 0, NULL);
cdoOperatorAdd("vertstd1", func_std1, 0, NULL);
cdoOperatorAdd("vertvar", func_var, 1, NULL);
cdoOperatorAdd("vertvar1", func_var1, 1, NULL);
cdoOperatorAdd("vertstd", func_std, 1, NULL);
cdoOperatorAdd("vertstd1", func_std1, 1, NULL);
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
......@@ -193,7 +193,7 @@ void *Vertstat(void *argument)
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
double divisor = operfunc == func_std1 || operfunc == func_var1;
int applyWeights = lmean;
//int applyWeights = lmean;
int streamID1 = streamOpenRead(cdoStreamName(0));
......@@ -340,10 +340,20 @@ void *Vertstat(void *argument)
vars1[varID].nmiss = nmiss;
if ( operatorID == VERTINT && IS_NOT_EQUAL(layer_thickness, 1.0) ) farcmul(&vars1[varID], layer_thickness);
if ( applyWeights && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&vars1[varID], layer_weight);
if ( lmean && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&vars1[varID], layer_weight);
if ( lvarstd )
farmoq(&vars2[varID], vars1[varID]);
{
if ( IS_NOT_EQUAL(layer_weight, 1.0) )
{
farmoqw(&vars2[varID], vars1[varID], layer_weight);
farcmul(&vars1[varID], layer_weight);
}
else
{
farmoq(&vars2[varID], vars1[varID]);
}
}
if ( nmiss > 0 || samp1[varID].ptr || needWeights )
{
......@@ -364,7 +374,7 @@ void *Vertstat(void *argument)
field.missval = vars1[varID].missval;
if ( operatorID == VERTINT && IS_NOT_EQUAL(layer_thickness, 1.0) ) farcmul(&field, layer_thickness);
if ( applyWeights && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&field, layer_weight);
if ( lmean && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&field, layer_weight);
if ( field.nmiss > 0 || samp1[varID].ptr )
{
......@@ -382,8 +392,16 @@ void *Vertstat(void *argument)
if ( lvarstd )
{
farsumq(&vars2[varID], field);
farsum(&vars1[varID], field);
if ( IS_NOT_EQUAL(layer_weight, 1.0) )
{
farsumqw(&vars2[varID], field, layer_weight);
farsumw(&vars1[varID], field, layer_weight);
}
else
{
farsumq(&vars2[varID], field);
farsum(&vars1[varID], field);
}
}
else
{
......
......@@ -145,7 +145,9 @@ void farfun(field_t *field1, field_t field2, const int function);
void faradd(field_t *field1, field_t field2);
void farsum(field_t *field1, field_t field2);
void farsumw(field_t *field1, field_t field2, double w);
void farsumq(field_t *field1, field_t field2);
void farsumqw(field_t *field1, field_t field2, double w);
void farsumtr(field_t *field1, field_t field2, const double refval);
void farsub(field_t *field1, field_t field2);
void farmul(field_t *field1, field_t field2);
......@@ -157,6 +159,7 @@ void farstd(field_t *field1, field_t field2, field_t field3, const double diviso
void farcvar(field_t *field1, field_t field2, const double rconst1, const double divisor);
void farcstd(field_t *field1, field_t field2, const double rconst1, const double divisor);
void farmoq(field_t *field1, field_t field2);
void farmoqw(field_t *field1, field_t field2, double w);
void faratan2(field_t *field1, field_t field2);
/* RQ */
......
......@@ -68,6 +68,12 @@ void arradd(const size_t n, double * restrict a, const double * restrict b)
#endif
}
static
void arraddw(const size_t n, double * restrict a, const double * restrict b, double w)
{
for ( size_t i = 0; i < n; i++ ) a[i] += w*b[i];
}
void faradd(field_t *field1, field_t field2)
{
......@@ -148,6 +154,49 @@ void farsum(field_t *field1, field_t field2)
}
}
void farsumw(field_t *field1, field_t field2, double w)
{
size_t i, len;
int nwpv = field1->nwpv;
const int grid1 = field1->grid;
const int nmiss1 = field1->nmiss;
const double missval1 = field1->missval;
double *array1 = field1->ptr;
/* double *weight1 = field1->weight; */
const int grid2 = field2.grid;
const int nmiss2 = field2.nmiss;
const double missval2 = field2.missval;
double *array2 = field2.ptr;
if ( nwpv != 2 ) nwpv = 1;
len = (size_t) (nwpv*gridInqSize(grid1));
if ( len != (size_t) (nwpv*gridInqSize(grid2)) )
cdoAbort("Fields have different gridsize (%s)", __func__);
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array2[i], missval2) )
{
if ( !DBL_IS_EQUAL(array1[i], missval1) )
array1[i] += w*array2[i];
else
array1[i] = w*array2[i];
}
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) ) field1->nmiss++;
}
else
{
arraddw(len, array1, array2, w);
}
}
/*
* Compute the occurrence of values in field, if they do not equal refval.
* This can be used to compute the lengths of multiple periods in a timeseries.
......@@ -247,6 +296,50 @@ void farsumq(field_t *field1, field_t field2)
}
void farsumqw(field_t *field1, field_t field2, double w)
{
size_t i, len;
int nwpv = field1->nwpv;
const int grid1 = field1->grid;
const int nmiss1 = field1->nmiss;
const double missval1 = field1->missval;
double *array1 = field1->ptr;
/* double *weight1 = field1->weight; */
const int grid2 = field2.grid;
const int nmiss2 = field2.nmiss;
const double missval2 = field2.missval;
double *array2 = field2.ptr;
if ( nwpv != 2 ) nwpv = 1;
len = (size_t) (nwpv*gridInqSize(grid1));
if ( len != (size_t) (nwpv*gridInqSize(grid2)) )
cdoAbort("Fields have different gridsize (%s)", __func__);
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array2[i], missval2) )
{
if ( !DBL_IS_EQUAL(array1[i], missval1) )
array1[i] += w*array2[i]*array2[i];
else
array1[i] = w*array2[i]*array2[i];
}
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) ) field1->nmiss++;
}
else
{
for ( i = 0; i < len; i++ )
array1[i] += w*array2[i]*array2[i];
}
}
void farsub(field_t *field1, field_t field2)
{
size_t i, len;
......@@ -807,6 +900,45 @@ void farmoq(field_t *field1, field_t field2)
}
void farmoqw(field_t *field1, field_t field2, double w)
{
size_t i, len;
int nwpv = field1->nwpv;
const int grid1 = field1->grid;
const double missval1 = field1->missval;
double *array1 = field1->ptr;
const int grid2 = field2.grid;
const int nmiss2 = field2.nmiss;
double missval2 = field2.missval;
double *array2 = field2.ptr;
if ( nwpv != 2 ) nwpv = 1;
len = (size_t) (nwpv*gridInqSize(grid1));
if ( len != (size_t) (nwpv*gridInqSize(grid2)) )
cdoAbort("Fields have different gridsize (%s)", __func__);
if ( nmiss2 > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array2[i], missval2) )
array1[i] = w*array2[i]*array2[i];
else
array1[i] = missval1;
field1->nmiss = 0;
for ( i = 0; i < len; i++ )
if ( DBL_IS_EQUAL(array1[i], missval1) ) field1->nmiss++;
}
else
{
for ( i = 0; i < len; i++ )
array1[i] = w*array2[i]*array2[i];
}
}
/* RQ */
/**
* Counts the number of nonmissing values. The result of the operation
......
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