Commit 3a4b3a96 authored by Uwe Schulzweida's avatar Uwe Schulzweida

Merge develop.

parents bf801b6d 3969ca9f
2021-10-25 Uwe Schulzweida
* Using CDI library version 2.0.0
* Using spherepart and clipping from YAC2
* Version 2.0.0 release
2021-02-25 Uwe Schulzweida
* Changed to YAC2
* Added operator fldint: field integral
2021-02-22 Uwe Schulzweida
* Merstat: Added memory support for 32-bit float data.
* Zonstat: Added memory support for 32-bit float data.
* Fldstat2: Added memory support for 32-bit float data.
2021-02-18 Uwe Schulzweida
* Seasstat: Added memory support for 32-bit float data.
2021-02-17 Uwe Schulzweida
* Yearmonstat: Added memory support for 32-bit float data.
* Ydaystat: Added memory support for 32-bit float data.
* Pack: added missing value support
2021-02-16 Uwe Schulzweida
* Runstat: Added memory support for 32-bit float data.
......
......@@ -7,6 +7,9 @@ Version 2.0.0 (25 October 2021):
New features:
* Changed to C++14
New operators:
* fldint: field integral
Fixed bugs:
Version 1.9.10 (25 January 2021):
......
......@@ -284,6 +284,7 @@ Operator catalog:
Fldstat fldmax Field maximum
Fldstat fldrange Field range
Fldstat fldsum Field sum
Fldstat fldint Field integral
Fldstat fldmean Field mean
Fldstat fldavg Field average
Fldstat fldstd Field standard deviation
......
......@@ -318,7 +318,7 @@ case "${HOSTNAME}" in
elif test "$COMP" = pgi ; then
${CONFPATH}configure --without-fortran \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,/sw/rhel6-x64/eccodes/eccodes-2.6.0-gcc64/lib" \
LDFLAGS="-Wl,-rpath,/sw/rhel6-x64/eccodes/eccodes-2.6.0-gcc64/lib -lpgf90 -lpgf90_rpm1 -lpgf902 -lpgf90rtl -lpgftnrtl -lnspgc -lpgc -lrt -lm -lgcc -lc -lgcc" \
F77=pgf90 \
CXX=pgc++ CXXFLAGS="-g -fast -D_REENTRANT" \
CC=pgcc CFLAGS="-g -fast -D_REENTRANT"
......
......@@ -212,7 +212,7 @@ colorlinks=true}
\end{picture}
\begin{flushright}
\large \textbf{Climate Data Operator \\ Version 1.9.9 \\ October 2020}
\large \textbf{Climate Data Operator \\ Version 2.0.0 \\ October 2021}
\end{flushright}
\vfill
......
......@@ -204,7 +204,7 @@ keepaspectratio]{cdo_libdep.pdf}}%
\end{picture}
\begin{flushright}
\large \textbf{Climate Data Operator \\ Version 1.9.9 \\ October 2020}
\large \textbf{Climate Data Operator \\ Version 2.0.0 \\ October 2021}
\end{flushright}
\vfill
......
......@@ -14,7 +14,7 @@
\put(0,0.0){\line(1,0){3.95}}
\end{picture}
\begin{flushright}
{\small{Climate Data Operator \\ Version 1.9.9 \\ October 2020}}
{\small{Climate Data Operator \\ Version 2.0.0 \\ October 2021}}
\end{flushright}
\vspace*{0mm}
......
......@@ -5,11 +5,11 @@
@Section = Statistical values
@Class = Statistic
@Arguments = infile outfile
@Operators = fldmin fldmax fldrange fldsum fldmean fldavg fldstd fldstd1 fldvar fldvar1 fldpctl
@Operators = fldmin fldmax fldrange fldsum fldint fldmean fldavg fldstd fldstd1 fldvar fldvar1 fldpctl
@BeginDescription
This module computes statistical values of all input fields. A field is a horizontal layer of a data variable.
According to the chosen operator the field minimum, maximum, range, sum, average, variance, standard deviation
According to the chosen operator the field minimum, maximum, range, sum, integral, average, variance, standard deviation
or a certain percentile is written to @file{outfile}.
@EndDescription
@EndModule
......@@ -95,6 +95,26 @@ o(t,1) = \mbox{\textbf{sum}}\{i(t,x'), x_1 < x' \leq x_n\}
@EndOperator
@BeginOperator_fldint
@Title = Field integral
@Parameter = weights
@BeginDescription
@IfMan
For every gridpoint x_1, ..., x_n of the same field it is:
o(t,1) = sum{i(t,x')*cellarea(x'), x_1<x'<=x_n}
@EndifMan
@IfDoc
For every gridpoint \begin{math}x_1, ..., x_n\end{math} of the same field it is: \\
@BeginMath
o(t,1) = \mbox{\textbf{sum}}\{i(t,x')*cellarea(x'), x_1 < x' \leq x_n\}
@EndMath
@EndifDoc
@EndDescription
@EndOperator
@BeginOperator_fldmean
@Title = Field mean
@Parameter = [weights]
......
......@@ -81,8 +81,6 @@ The following options are available for all operators:
\> NetCDF4 chunk type: auto, grid or lines. \\
\makebox[1.5in][l]{\hspace*{1cm}\textsl{-L}}
\> Lock I/O (sequential access). \\
\makebox[1.5in][l]{\hspace*{1cm}\textsl{-M}}
\> Switch to indicate that the I/O streams have missing values. \\
\makebox[1.5in][l]{\hspace*{1cm}\textsl{-m $<$missval$>$}}
\> Set the missing value of non NetCDF files (default: \texttt{-9e+33}). \\
\makebox[1.5in][l]{\hspace*{1cm}\textsl{-O}}
......
......@@ -77,7 +77,6 @@ static int DataType = -1;
static char *filename;
static const char **ifiles;
static char *ifile = nullptr;
static const char *ofile = nullptr;
static int ofileidx = 0;
......@@ -1888,7 +1887,6 @@ Afterburner::run(void *process)
ofileidx = nfiles;
ifile = strdup(cdoGetStreamName(0));
ofile = cdoGetStreamName(nfiles);
globs.Nfiles = nfiles - 1;
if (globs.Nfiles > 0)
......
......@@ -42,6 +42,8 @@
#include "pmlist.h"
#include "cdo_zaxis.h"
void grid_cell_area(int gridID, double *array);
template <typename T>
static void
print_location_LL(int operfunc, int vlistID, int varID, int levelID, int gridID, T sglval, const Varray<T> &fieldvec,
......@@ -85,6 +87,39 @@ print_location_LL(int operfunc, int vlistID, int varID, int levelID, int gridID,
}
}
template <typename T>
static void
field_mul_weights(size_t len, Varray<T> &v1, const Varray<double> &v2, size_t nmiss, T missval)
{
if (nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(v1[i], missval)) v1[i] *= v2[i];
}
else
{
for (size_t i = 0; i < len; ++i) v1[i] *= v2[i];
}
}
static void
field_mul_weights(Field &field)
{
if (field.memType == MemType::Float)
field_mul_weights(field.size, field.vec_f, field.weightv, field.nmiss, (float)field.missval);
else
field_mul_weights(field.size, field.vec_d, field.weightv, field.nmiss, field.missval);
}
static void
printWeightsWarning(int ngrids, const char *varname)
{
if (ngrids == 1)
cdoWarning("Grid cell bounds not available, using constant grid cell area weights!");
else
cdoWarning("Grid cell bounds not available, using constant grid cell area weights for variable %s!", varname);
}
static void
fldstatGetParameter(bool &useweights)
{
......@@ -113,41 +148,33 @@ fldstatGetParameter(bool &useweights)
}
}
static void
addOperators(void)
{
// clang-format off
cdoOperatorAdd("fldrange", func_range, 0, nullptr);
cdoOperatorAdd("fldmin", func_min, 0, nullptr);
cdoOperatorAdd("fldmax", func_max, 0, nullptr);
cdoOperatorAdd("fldsum", func_sum, 0, nullptr);
cdoOperatorAdd("fldmean", func_meanw, 1, nullptr);
cdoOperatorAdd("fldavg", func_avgw, 1, nullptr);
cdoOperatorAdd("fldstd", func_stdw, 1, nullptr);
cdoOperatorAdd("fldstd1", func_std1w, 1, nullptr);
cdoOperatorAdd("fldvar", func_varw, 1, nullptr);
cdoOperatorAdd("fldvar1", func_var1w, 1, nullptr);
cdoOperatorAdd("fldskew", func_skew, 0, nullptr);
cdoOperatorAdd("fldkurt", func_kurt, 0, nullptr);
cdoOperatorAdd("fldpctl", func_pctl, 0, nullptr);
// clang-format on
}
void *
Fldstat::run(void *process)
{
int lastgrid = -1;
int nrecs;
int varID, levelID;
cdoInitialize(process);
addOperators();
// clang-format off
cdoOperatorAdd("fldrange", func_range, 0, nullptr);
cdoOperatorAdd("fldmin", func_min, 0, nullptr);
cdoOperatorAdd("fldmax", func_max, 0, nullptr);
cdoOperatorAdd("fldsum", func_sum, 0, nullptr);
auto FLDINT = cdoOperatorAdd("fldint", func_sum, 0, nullptr);
cdoOperatorAdd("fldmean", func_meanw, 1, nullptr);
cdoOperatorAdd("fldavg", func_avgw, 1, nullptr);
cdoOperatorAdd("fldstd", func_stdw, 1, nullptr);
cdoOperatorAdd("fldstd1", func_std1w, 1, nullptr);
cdoOperatorAdd("fldvar", func_varw, 1, nullptr);
cdoOperatorAdd("fldvar1", func_var1w, 1, nullptr);
cdoOperatorAdd("fldskew", func_skew, 0, nullptr);
cdoOperatorAdd("fldkurt", func_kurt, 0, nullptr);
cdoOperatorAdd("fldpctl", func_pctl, 0, nullptr);
// clang-format on
const auto operatorID = cdoOperatorID();
const auto operfunc = cdoOperatorF1(operatorID);
const bool needWeights = cdoOperatorF2(operatorID) != 0;
const auto needWeights = (cdoOperatorF2(operatorID) != 0);
bool useweights = true;
const auto needCellarea = (operatorID == FLDINT);
double pn = 0;
if (operfunc == func_pctl)
......@@ -164,7 +191,6 @@ Fldstat::run(void *process)
operatorCheckArgc(0);
}
const auto streamID1 = cdoOpenRead(0);
const auto vlistID1 = cdoStreamInqVlist(streamID1);
......@@ -192,10 +218,10 @@ Fldstat::run(void *process)
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
if (needWeights)
if (needWeights || needCellarea)
{
field.weightv.resize(gridsizemax);
if (!useweights)
if (needWeights && !useweights)
{
cdoPrint("Using constant grid cell area weights!");
for (size_t i = 0; i < gridsizemax; ++i) field.weightv[i] = 1;
......@@ -205,9 +231,13 @@ Fldstat::run(void *process)
VarList varList1;
varListInit(varList1, vlistID1);
int lastgrid = -1;
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID);
......@@ -216,6 +246,7 @@ Fldstat::run(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
int varID, levelID;
cdoInqRecord(streamID1, &varID, &levelID);
field.init(varList1[varID]);
cdoReadRecord(streamID1, field);
......@@ -226,18 +257,19 @@ Fldstat::run(void *process)
field.weightv[0] = 1;
if (useweights && field.size > 1)
{
const bool wstatus = gridWeights(field.grid, field.weightv.data()) != 0;
if (wstatus && tsID == 0 && levelID == 0)
{
if (ngrids == 1)
cdoWarning("Grid cell bounds not available, using constant grid cell area weights!");
else
cdoWarning("Grid cell bounds not available, using constant grid cell area weights for variable %s!",
varList1[varID].name);
}
const auto wstatus = (gridWeights(field.grid, field.weightv.data()) != 0);
if (wstatus && tsID == 0 && levelID == 0) printWeightsWarning(ngrids, varList1[varID].name);
}
}
if (needCellarea && field.grid != lastgrid)
{
lastgrid = field.grid;
grid_cell_area(field.grid, field.weightv.data());
}
if (needCellarea) field_mul_weights(field);
auto singleValue = (operfunc == func_pctl) ? fieldPctl(field, pn) : fieldFunction(field, operfunc);
if (Options::cdoVerbose && (operfunc == func_min || operfunc == func_max))
......
......@@ -31,21 +31,24 @@
// routine corr copied from PINGO
// correclation in space
template <typename T>
static double
correlation_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<double> &weight, double missval1,
correlation_s(const Varray<T> &v1, const Varray<T> &v2, const Varray<double> &weight, double missval1,
double missval2, size_t gridsize)
{
double sum0 = 0.0, sum1 = 0.0, sum00 = 0.0, sum01 = 0.0, sum11 = 0.0, wsum0 = 0.0;
for (size_t i = 0; i < gridsize; ++i)
{
if (IS_NOT_EQUAL(weight[i], missval1) && IS_NOT_EQUAL(in0[i], missval1) && IS_NOT_EQUAL(in1[i], missval2))
double vv1 = v1[i];
double vv2 = v2[i];
if (IS_NOT_EQUAL(weight[i], missval1) && IS_NOT_EQUAL(vv1, missval1) && IS_NOT_EQUAL(vv2, missval2))
{
sum0 += weight[i] * in0[i];
sum1 += weight[i] * in1[i];
sum00 += weight[i] * in0[i] * in0[i];
sum01 += weight[i] * in0[i] * in1[i];
sum11 += weight[i] * in1[i] * in1[i];
sum0 += weight[i] * vv1;
sum1 += weight[i] * vv2;
sum00 += weight[i] * vv1 * vv1;
sum01 += weight[i] * vv1 * vv2;
sum11 += weight[i] * vv2 * vv2;
wsum0 += weight[i];
}
}
......@@ -57,20 +60,32 @@ correlation_s(const Varray<double> &in0, const Varray<double> &in1, const Varray
return out;
}
static double
correlation_s(const Field &field1, const Field &field2, const Varray<double> &weight)
{
if (field1.memType == MemType::Float)
return correlation_s(field1.vec_f, field2.vec_f, weight, field1.missval, field2.missval, field1.size);
else
return correlation_s(field1.vec_d, field2.vec_d, weight, field1.missval, field2.missval, field1.size);
}
// covariance in space
template <typename T>
static double
covariance_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<double> &weight, double missval1,
covariance_s(const Varray<T> &v1, const Varray<T> &v2, const Varray<double> &weight, double missval1,
double missval2, size_t gridsize)
{
double sum0 = 0.0, sum1 = 0.0, sum01 = 0.0, wsum0 = 0.0, wsum00 = 0.0;
for (size_t i = 0; i < gridsize; ++i)
{
if (IS_NOT_EQUAL(weight[i], missval1) && IS_NOT_EQUAL(in0[i], missval1) && IS_NOT_EQUAL(in1[i], missval2))
double vv1 = v1[i];
double vv2 = v2[i];
if (IS_NOT_EQUAL(weight[i], missval1) && IS_NOT_EQUAL(vv1, missval1) && IS_NOT_EQUAL(vv2, missval2))
{
sum0 += weight[i] * in0[i];
sum1 += weight[i] * in1[i];
sum01 += weight[i] * in0[i] * in1[i];
sum0 += weight[i] * vv1;
sum1 += weight[i] * vv2;
sum01 += weight[i] * vv1 * vv2;
wsum0 += weight[i];
wsum00 += weight[i] * weight[i];
}
......@@ -81,13 +96,18 @@ covariance_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<
return out;
}
static double
covariance_s(const Field &field1, const Field &field2, const Varray<double> &weight)
{
if (field1.memType == MemType::Float)
return covariance_s(field1.vec_f, field2.vec_f, weight, field1.missval, field2.missval, field1.size);
else
return covariance_s(field1.vec_d, field2.vec_d, weight, field1.missval, field2.missval, field1.size);
}
void *
Fldstat2::run(void *process)
{
int lastgridID = -1;
int nrecs;
int varID, levelID;
size_t nmiss1, nmiss2;
bool wstatus = false;
bool needWeights = true;
......@@ -132,17 +152,20 @@ Fldstat2::run(void *process)
const auto streamID3 = cdoOpenWrite(2);
cdoDefVlist(streamID3, vlistID3);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Field field1, field2;
Varray<double> array1(gridsizemax), array2(gridsizemax);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> weight;
if (needWeights) weight.resize(gridsizemax);
int lastgridID = -1;
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs2 = cdoStreamInqTimestep(streamID2, tsID);
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
const auto nrecs2 = cdoStreamInqTimestep(streamID2, tsID);
if (nrecs2 == 0)
{
cdoWarning("Input streams have different number of time steps!");
......@@ -154,13 +177,15 @@ Fldstat2::run(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
int varID, levelID;
cdoInqRecord(streamID1, &varID, &levelID);
field1.init(varList1[varID]);
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID1, array1.data(), &nmiss1);
cdoReadRecord(streamID2, array2.data(), &nmiss2);
field2.init(varList2[varID]);
cdoReadRecord(streamID1, field1);
cdoReadRecord(streamID2, field2);
const auto gridID = varList1[varID].gridID;
const auto gridsize = varList1[varID].gridsize;
if (needWeights && gridID != lastgridID)
{
lastgridID = gridID;
......@@ -169,16 +194,13 @@ Fldstat2::run(void *process)
if (wstatus && tsID == 0 && levelID == 0)
cdoWarning("Using constant grid cell area weights for variable %s!", varList1[varID].name);
const auto missval1 = varList1[varID].missval;
const auto missval2 = varList2[varID].missval;
double sglval = 0.0;
if (operfunc == func_cor)
sglval = correlation_s(array1, array2, weight, missval1, missval2, gridsize);
sglval = correlation_s(field1, field2, weight);
else if (operfunc == func_covar)
sglval = covariance_s(array1, array2, weight, missval1, missval2, gridsize);
sglval = covariance_s(field1, field2, weight);
const auto nmiss3 = DBL_IS_EQUAL(sglval, missval1) ? 1 : 0;
const auto nmiss3 = DBL_IS_EQUAL(sglval, varList1[varID].missval) ? 1 : 0;
cdoDefRecord(streamID3, varID, levelID);
cdoWriteRecord(streamID3, &sglval, nmiss3);
......
......@@ -387,7 +387,6 @@ intlevelGetParameter(std::vector<double> &lev2, std::string &zdescription, std::
void *
Intlevel::run(void *process)
{
int nrecs;
int varID, levelID;
int zaxisID1 = -1;
int nlevel = 0;
......@@ -578,8 +577,11 @@ Intlevel::run(void *process)
}
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
for (varID = 0; varID < nvars; ++varID) vars[varID] = false;
taxisCopyTimestep(taxisID2, taxisID1);
......
......@@ -87,7 +87,6 @@ void *
Intlevel3d::run(void *process)
{
size_t gridsizeo;
int nrecs;
int i;
int varID = -1, levelID = -1;
int gridID3 = -1;
......@@ -147,7 +146,7 @@ Intlevel3d::run(void *process)
nlevi = zaxisInqSize(zaxisID); // number of input levels for later use
zlevels_in.resize(gridsizei * nlevi);
nrecs = cdoStreamInqTimestep(streamID2, 0);
const auto nrecs = cdoStreamInqTimestep(streamID2, 0);
if (Options::cdoVerbose) cdoPrint("%d records input 3d vertical height", nrecs);
for (int recID = 0; recID < nrecs; recID++)
......@@ -178,7 +177,7 @@ Intlevel3d::run(void *process)
zlevels_out.resize(gridsize * nlevel);
nlevo = nlevel; // number of output levels for later use
gridsizeo = gridsize; // horizontal gridsize of output z coordinate
nrecs = streamInqTimestep(streamID0, 0);
const auto nrecs = streamInqTimestep(streamID0, 0);
if (Options::cdoVerbose) cdoPrint("%d records target 3d vertical height and gridsize %zu", nrecs, gridsize);
for (int recID = 0; recID < nrecs; recID++)
......@@ -314,8 +313,11 @@ Intlevel3d::run(void *process)
if (varID == nvars) cdoAbort("No processable variable found!");
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
for (varID = 0; varID < nvars; ++varID) vars[varID] = false;
taxisCopyTimestep(taxisID3, taxisID1);
......
......@@ -64,8 +64,6 @@ Merstat::run(void *process)
{
int gridID1, gridID2 = -1, lastgrid = -1;
int index;
int nrecs;
int varID, levelID;
cdoInitialize(process);
......@@ -124,29 +122,32 @@ Merstat::run(void *process)
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Field field1, field2;
field1.resize(gridsizemax);
if (needWeights) field1.weightv.resize(gridsizemax);
field2.resize(nlonmax);
field2.grid = gridID2;
field2.memType = MemType::Double;
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID);