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

Merge develop.

parents 5e6ca0ff 875136a5
2021-01-29 Uwe Schulzweida
* Using CDI library version 1.9.10
* Version 1.9.10 release
2021-01-11 Uwe Schulzweida
* EOF: data race, wrong result with multiple OpenMP threads (bug fix)
2021-01-05 Uwe Schulzweida
* Added option --ignore_time_bounds to ignore time bounds for time range statistics
* Isosurface: Added memory support for 32-bit float data.
2020-12-17 Uwe Schulzweida
* Added warning message if a non-thread-safe NetCDF4/HDF5 library is used
2020-12-15 Uwe Schulzweida
* Exprf: added support for same variable name with different number of levels
2020-11-23 Uwe Schulzweida
* Ymonarith: failed with variables on different grids
2020-11-20 Uwe Schulzweida
* timselmean: failed with variables on different grids [Bug #9978]
2020-11-18 Uwe Schulzweida
* Detrend: wrong result with parameter equal=false [Bug #9961]
* subtrend: added parameter equal=false
2020-11-12 Uwe Schulzweida
* Fldstat: optional parameter weights failed (bug fix)
* Wind: check that numLPE is > 0 (bug fix)
2020-10-29 Uwe Schulzweida
* Using CDI library version 1.9.9
......
......@@ -3,6 +3,18 @@ CDO NEWS
Improvement
Version 1.9.10 (29 January 2021):
New features:
* Added option --ignore_time_bounds to ignore time bounds for time range statistics
Fixed bugs:
* EOF: fix wrong result with multiple OpenMP threads (data race)
* timselmean: failed with variables on different grids [Bug #9978]
* Ymonarith: failed with variables on different grids
* Detrend: wrong result with parameter equal=false [Bug #9961]
* Fldstat: optional parameter weights failed
* Wind: check that numLPE is > 0
Version 1.9.9 (29 October 2020):
New features:
......
......@@ -5,7 +5,7 @@
# libtool 2.4.2
AC_PREREQ([2.68])
AC_INIT([cdo], [1.9.9], [https://mpimet.mpg.de/cdo])
AC_INIT([cdo], [1.9.10], [https://mpimet.mpg.de/cdo])
AC_DEFINE_UNQUOTED(CDO, ["$PACKAGE_VERSION"], [CDO version])
......@@ -168,7 +168,7 @@ AS_IF([test x$acx_cv_cfortran_works = 'xno'],[AC_SUBST([FORTRAN_WORKS],[no])],[A
# ----------------------------------------------------------------------
CFLAGS="$CFLAGS ${OPENMP_CFLAGS}"
CXXFLAGS="$CXXFLAGS ${OPENMP_CFLAGS}"
----------------------------------------------------------------------
# ----------------------------------------------------------------------
# Add configure options
ACX_CDO_OPTIONS
ACX_CDI_OPTIONS
......@@ -215,6 +215,7 @@ AC_CONFIG_FILES([
test/Filter.test
test/Fldpctl.test
test/Fldstat.test
test/Fldstat2.test
test/Genweights.test
test/Gradsdes.test
test/Gridarea.test
......@@ -222,6 +223,7 @@ AC_CONFIG_FILES([
test/Inttime.test
test/Isosurface.test
test/MapReduce.test
test/Merge.test
test/Mergetime.test
test/Merstat.test
test/Monarith.test
......
......@@ -82,7 +82,7 @@ INTEGER Comma-separated list or first/last[/inc] range of years.
@Item = dom
STRING Comma-separated list of the day of month (e.g. 29feb).
@Item = timestep
INTEGER Comma-separated list or first/last[/inc] range of timesteps. Negative values selects timesteps from the end (NetCDF only).
INTEGER Comma-separated list or first/last[/inc] range of timesteps. Negative values select timesteps from the end (NetCDF only).
@Item = timestep_of_year
INTEGER Comma-separated list or first/last[/inc] range of timesteps of year.
@Item = timestepmask
......
......@@ -107,7 +107,7 @@ Selects a month and optional an arbitrary number of timesteps before and after t
@BeginParameter
@Item = timesteps
INTEGER Comma-separated list or first/last[/inc] range of timesteps. Negative values selects timesteps from the end (NetCDF only).
INTEGER Comma-separated list or first/last[/inc] range of timesteps. Negative values select timesteps from the end (NetCDF only).
@Item = times
STRING Comma-separated list of times (format hh:mm:ss).
@Item = hours
......
@BeginModule
@NewPage
@Name = Trend
@Title = Trend of time series
@Section = Regression
......
......@@ -15,6 +15,7 @@ This module is for adding or subtracting a trend computed by the operator @mod{t
@BeginOperator_addtrend
@Title = Add trend
@Parameter = [equal]
@BeginDescription
@IfMan
......@@ -35,6 +36,7 @@ where t is the timesteps.
@BeginOperator_subtrend
@Title = Subtract trend
@Parameter = [equal]
@BeginDescription
@IfMan
......@@ -53,6 +55,12 @@ where t is the timesteps.
@EndOperator
@BeginParameter
@Item = equal
BOOL Set to false for unequal distributed timesteps (default: true)
@EndParameter
@BeginExample
The typical call for detrending the data in @file{infile} and storing the detrended data in @file{outfile} is:
@BeginVerbatim
......
......@@ -25,8 +25,9 @@ Supported parameter tables are: WMO standard table number 2 and ECMWF local tabl
geopotential_height & m & 156
@EndTable
Use the alias @bold{ml2plx}/@bold{ml2hlx} or the environment variable @env{EXTRAPOLATE}
to extrapolate missing values. This operator requires all variables on the same horizontal grid.
Use the alias @bold{ml2plx}/@bold{ml2hlx} or the environment variable @env{EXTRAPOLATE} to extrapolate
missing values. This operator requires all variables on the same horizontal grid.
Missing values in the input data are not supported.
@EndDescription
@EndModule
......
......@@ -48,6 +48,7 @@
#include "cdo_options.h"
#include "mpmo_color.h"
#include "cdi_lockedIO.h"
#include "gaussian_latitudes.h"
int scan_par_obsolete(char *namelist, const char *name, int def);
int scan_par(int verbose, char *namelist, const char *name, int def);
......@@ -1018,7 +1019,7 @@ after_defineGrid(const AfterControl &globs, struct Variable *vars)
const size_t nlat = globs.Latitudes;
std::vector<double> lats(nlat), latw(nlat);
gaussianLatitudes(lats.data(), latw.data(), nlat);
gaussian_latitudes(nlat, lats.data(), latw.data());
for (size_t j = 0; j < nlat; ++j) lats[j] = 180. / M_PI * std::asin(lats[j]);
gridDefYvals(gaussGridID, lats.data());
}
......
......@@ -98,6 +98,9 @@ Arithdays(void *process)
const auto streamID2 = cdoOpenWrite(1);
cdoDefVlist(streamID2, vlistID2);
VarList varList1;
varListInit(varList1, vlistID1);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
......@@ -126,9 +129,9 @@ Arithdays(void *process)
{
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, field.vec_d.data(), &field.nmiss);
field.grid = vlistInqVarGrid(vlistID1, varID);
field.missval = vlistInqVarMissval(vlistID1, varID);
field.size = varList1[varID].gridsize;
field.grid = varList1[varID].gridID;
field.missval = varList1[varID].missval;
vfarcfun(field, rconst, operfunc);
......
......@@ -54,7 +54,7 @@ detrend(const long nts, const Varray<double> &deltaTS0, const double missval1, c
const auto work1 = DIVMN(SUBMN(sumjx, DIVMN(MULMN(sumj, sumx), n)), SUBMN(sumjj, DIVMN(MULMN(sumj, sumj), n)));
const auto work2 = SUBMN(DIVMN(sumx, n), MULMN(DIVMN(sumj, n), work1));
for (long j = 0; j < nts; j++) array2[j] = SUBMN(array1[j], ADDMN(work2, MULMN(j, work1)));
for (long j = 0; j < nts; j++) array2[j] = SUBMN(array1[j], ADDMN(work2, MULMN(work1, deltaTS0[j])));
}
static void
......
......@@ -44,7 +44,7 @@
#include "datetime.h"
#include "eof_mode.h"
// NO MISSING VALUE SUPPORT ADDED SO FAR
// No missing value support added so far!
static void
scale_eigvec_grid(Varray<double> &out, int tsID, size_t npack, const std::vector<size_t> &pack, const Varray<double> &weight,
......@@ -169,9 +169,9 @@ EOFs(void *process)
for (int index = 1; index < ngrids; index++)
if (vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index)) cdoAbort("Too many different grids!");
/* eigenvalues */
// eigenvalues
/* COUNT NUMBER OF TIMESTEPS if EOF_ or EOF_TIME */
// Count number of timesteps if EOF_ or EOF_TIME
if (operfunc == EOF_ || operfunc == EOF_TIME)
{
if (Options::cdoVerbose) cdoPrint("Counting timesteps in ifile");
......@@ -471,13 +471,13 @@ EOFs(void *process)
eofdata[varID][levelID].eig_val.resize(nts);
#ifdef _OPENMP
#pragma omp parallel for default(shared) schedule(dynamic)
#endif
for (int j1 = 0; j1 < nts; j1++)
{
const auto &df1p = data[j1];
for (int j2 = 0; j2 < j1; j2++) covar[j1][j2] = covar[j2][j1];
#ifdef _OPENMP
#pragma omp parallel for default(shared) schedule(dynamic)
#endif
for (int j2 = j1; j2 < nts; j2++)
{
const auto &df2p = data[j2];
......@@ -492,7 +492,7 @@ EOFs(void *process)
if (Options::Timer) timer_stop(timer_cov);
// SOLVE THE EIGEN PROBLEM
// Solve the eigen problem
if (Options::Timer) timer_start(timer_eig);
auto &eig_val = eofdata[varID][levelID].eig_val;
......
......@@ -144,7 +144,7 @@ static void calculateOuterPeriod(Field &field, int MaxMonths, int recentYear, in
if (cdoAssertFilesOnly() == false)
cdoAbort("infile1 cannot be a pipe");
const auto maxrecs = vlistNrecs(vlistID1);
CdiStreamID cdiStream = streamOpenRead(cdoGetStreamName(0));
auto cdiStream = streamOpenRead(cdoGetStreamName(0));
const auto cdiVlistID = streamInqVlist(cdiStream);
const auto cdiTaxisID = vlistInqTaxis(cdiVlistID);
......@@ -152,14 +152,17 @@ static void calculateOuterPeriod(Field &field, int MaxMonths, int recentYear, in
int tempdpy = 0;
for ( int i = 0; i<MaxMonths; i++)
tempdpm[i] = 0;
int year, month, day, tsID = 0, nrecs = 0, varID, levelID;
int year, month, day, tsID = 0, varID, levelID;
bool lHasStarted = false;
if (Options::cdoVerbose) cdoPrint("Start to process variables");
while ( ( nrecs = streamInqTimestep(cdiStream, tsID++) ) )
while (true)
{
int64_t vdate = taxisInqVdate(cdiTaxisID);
const auto nrecs = streamInqTimestep(cdiStream, tsID++);
if (nrecs == 0) break;
auto vdate = taxisInqVdate(cdiTaxisID);
cdiDecodeDate(vdate, &year, &month, &day);
if ( !lHasStarted && year != recentYear )
continue;
......@@ -180,14 +183,13 @@ static void calculateOuterPeriod(Field &field, int MaxMonths, int recentYear, in
fieldFill(cei[loopmonth][0][0], 0.);
}
}
if ( year == endOfCalc && func2 == func_avg )
break;
if (year == endOfCalc && func2 == func_avg) break;
tempdpy++;
int dayoy = (month >= 1 && month <= 12) ? (month - 1) * 31 + day : 0;
tempdpm[month-1]++;
if ( func2 == func_sum )
dayoy = 1;
if ( func2 == func_sum ) dayoy = 1;
for (int recID = 0; recID < nrecs; recID++)
{
streamInqRecord(cdiStream, &varID, &levelID);
......@@ -217,7 +219,7 @@ static void calculateOuterPeriod(Field &field, int MaxMonths, int recentYear, in
fieldFill(cei[0][0][0], 0.);
if ( frequency == 8 )
for ( int loopmonth = 1; loopmonth < MaxMonths; loopmonth++)
for (int loopmonth = 1; loopmonth < MaxMonths; loopmonth++)
{
tempdpm[loopmonth] = 0;
fieldFill(cei[loopmonth][0][0], 0.);
......@@ -232,7 +234,6 @@ etccdi_op(ETCCDI_REQUEST *request)
constexpr int MaxDays = 373;
constexpr int MaxMonths = 12;
int varID;
int nrecs;
int levelID;
size_t nmiss;
int year, month, day, dayoy;
......@@ -336,8 +337,11 @@ etccdi_op(ETCCDI_REQUEST *request)
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID2, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID2, tsID);
if (nrecs == 0) break;
if (nrecs != cdoStreamInqTimestep(streamID3, tsID))
cdoAbort("Number of records at time step %d of %s and %s differ!", tsID + 1, cdoGetStreamName(1),
cdoGetStreamName(2));
......@@ -396,6 +400,7 @@ etccdi_op(ETCCDI_REQUEST *request)
else
cdoReadRecord(streamID3, field.vec_d.data(), &nmiss);
field.nmiss = nmiss;
field.size = vars2[dayoy][varID][levelID].size;
field.grid = vars2[dayoy][varID][levelID].grid;
field.missval = vars2[dayoy][varID][levelID].missval;
......@@ -414,8 +419,11 @@ etccdi_op(ETCCDI_REQUEST *request)
tsID = 0;
bool lOnlyRefPeriod = true;
int firstYear = 0, lastYear = 0;
while ( (nrecs = cdoStreamInqTimestep(streamID1, tsID++) ) )
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID++);
if (nrecs == 0) break;
vdate = taxisInqVdate(taxisID1);
vtime = taxisInqVtime(taxisID1);
......@@ -434,7 +442,7 @@ etccdi_op(ETCCDI_REQUEST *request)
if (dayoy < 0 || dayoy >= MaxDays) cdoAbort("Day %d out of range!", dayoy);
if (wdaysSrc[dayoy] || request->func2 == func_sum )
{
/* Variable independent ? */
/* Variable independent ? */
wdaysRead[dayoy+(year-request->startboot)*(MaxDays-1)] = dayoy+(year-request->startboot)*(MaxDays-1);
dpy[year-request->startboot]++;
dpm[(year-request->startboot)*MaxMonths+(int)((dayoy-1)/31.)]++;
......@@ -792,6 +800,7 @@ etccdi_op(ETCCDI_REQUEST *request)
}
field.resize(gridsizemax);
field.missval = vars1[1][0][0].missval;
field.size = vars1[1][0][0].size;
field.grid = vars1[1][0][0].grid;
calculateOuterPeriod(field, MaxMonths, request->endboot+1, lastYear+1, cei, varsPtemp, frequency, taxisID4, streamID4, &otsID, vlistID1, recinfo, selection, request->func2);
}
......
......@@ -104,7 +104,7 @@ exprsFromFile(const std::vector<std::string> &exprArgv)
if (exprArgc != 1) operatorCheckArgc(1);
auto exprf = exprArgv[0].c_str();
/* Open expr script file for reading */
// Open expr script file for reading
auto fp = fopen(exprf, "r");
if (!fp) cdoAbort("Open failed on %s", exprf);
......@@ -136,7 +136,7 @@ str_replace(char *target, const char *needle, const char *replacement)
const auto needle_len = strlen(needle);
const auto repl_len = strlen(replacement);
while (1)
while (true)
{
char *p = strstr(tmp, needle);
......@@ -222,6 +222,7 @@ params_init(std::vector<paramType> &params, const VarList &varList, int vlistID)
vlistInqVarUnits(vlistID, varID, units);
params[varID].type = PARAM_VAR;
params[varID].isValid = true;
params[varID].select = false;
params[varID].remove = false;
params[varID].lmiss = true;
......@@ -329,11 +330,9 @@ params_add_coordinates(int vlistID, parseParamType &parse_arg)
static int
params_add_ts(parseParamType &parse_arg)
{
int varID = -1;
auto params = parse_arg.params;
if (params)
{
varID = parse_arg.nparams;
auto &params = parse_arg.params;
auto varID = parse_arg.nparams;
if (varID >= parse_arg.maxparams) cdoAbort("Too many parameter (limit=%d)", parse_arg.maxparams);
params[varID].name = strdup("_ts");
......@@ -345,13 +344,12 @@ params_add_ts(parseParamType &parse_arg)
parse_arg.nparams++;
parse_arg.cnparams++;
}
return varID;
}
static void
parseParamInit(parseParamType &parse_arg, int vlistID, int pointID, int zonalID, int surfaceID, paramType *params)
parseParamInit(parseParamType &parse_arg, int vlistID, int pointID, int zonalID, int surfaceID)
{
const auto nvars = vlistNvars(vlistID);
const auto ngrids = vlistNgrids(vlistID);
......@@ -359,12 +357,12 @@ parseParamInit(parseParamType &parse_arg, int vlistID, int pointID, int zonalID,
const auto maxcoords = ngrids * 4 + nzaxis;
parse_arg.maxparams = MAX_PARAMS;
parse_arg.params.resize(MAX_PARAMS);
parse_arg.nparams = nvars;
parse_arg.cnparams = nvars;
parse_arg.nvars1 = nvars;
parse_arg.init = true;
parse_arg.debug = Options::cdoVerbose != 0;
parse_arg.params = params;
parse_arg.pointID = pointID;
parse_arg.zonalID = zonalID;
parse_arg.surfaceID = surfaceID;
......@@ -395,6 +393,47 @@ genZonalID(int vlistID)
return zonalID;
}
static
void setDateAndTime(paramType &varts, int calendar, int tsID, int64_t vdate0, int vtime0, int64_t vdate, int vtime)
{
double jdelta = 0.0;
if (tsID)
{
const auto juldate0 = julianDateEncode(calendar, vdate0, vtime0);
const auto juldate = julianDateEncode(calendar, vdate, vtime);
jdelta = julianDateToSeconds(julianDateSub(juldate, juldate0));
}
varts.data[CTIMESTEP] = tsID + 1;
varts.data[CDATE] = vdate;
varts.data[CTIME] = vtime;
varts.data[CDELTAT] = jdelta;
int year, mon, day;
int hour, minute, second;
cdiDecodeDate(vdate, &year, &mon, &day);
cdiDecodeTime(vtime, &hour, &minute, &second);
varts.data[CDAY] = day;
varts.data[CMONTH] = mon;
varts.data[CYEAR] = year;
varts.data[CSECOND] = second;
varts.data[CMINUTE] = minute;
varts.data[CHOUR] = hour;
}
static void
addOperators(void)
{
// clang-format off
cdoOperatorAdd("expr", 1, 1, "expressions");
cdoOperatorAdd("exprf", 1, 0, "expr script filename");
cdoOperatorAdd("aexpr", 0, 1, "expressions");
cdoOperatorAdd("aexprf", 0, 0, "expr script filename");
// clang-format on
}
void *
Expr(void *process)
{
......@@ -406,12 +445,7 @@ Expr(void *process)
parseParamType parse_arg = {};
yyset_extra(&parse_arg, scanner);
// clang-format off
cdoOperatorAdd("expr", 1, 1, "expressions");
cdoOperatorAdd("exprf", 1, 0, "expr script filename");
cdoOperatorAdd("aexpr", 0, 1, "expressions");
cdoOperatorAdd("aexprf", 0, 0, "expr script filename");
// clang-format on
addOperators();
const auto operatorID = cdoOperatorID();
const bool replacesVariables = cdoOperatorF1(operatorID);
......@@ -438,10 +472,10 @@ Expr(void *process)
const auto zonalID = genZonalID(vlistID1);
const auto surfaceID = getSurfaceID(vlistID1);
std::vector<paramType> params(MAX_PARAMS);
params_init(params, varList1, vlistID1);
parseParamInit(parse_arg, vlistID1, pointID, zonalID, surfaceID);
parseParamInit(parse_arg, vlistID1, pointID, zonalID, surfaceID, params.data());
auto &params = parse_arg.params;
params_init(params, varList1, vlistID1);
// Set all input variables to 'needed' if replacing is switched off
for (int varID = 0; varID < nvars1; varID++) parse_arg.needed[varID] = !replacesVariables;
......@@ -551,7 +585,7 @@ Expr(void *process)
{
if (parse_arg.coords[i].needed)
{
const int coord = parse_arg.coords[i].coord;
const auto coord = parse_arg.coords[i].coord;
if (coord == 'x' || coord == 'y' || coord == 'a' || coord == 'w')
{
auto gridID = parse_arg.coords[i].cdiID;
......@@ -607,7 +641,7 @@ Expr(void *process)
for (int varID = parse_arg.nvars1; varID < parse_arg.nparams; varID++)
{
const int coord = params[varID].coord;
const auto coord = params[varID].coord;
if (coord)
{
if (coord == 'x' || coord == 'y' || coord == 'a' || coord == 'w')
......@@ -650,46 +684,26 @@ Expr(void *process)
int vtime0 = 0;
const auto calendar = taxisInqCalendar(taxisID1);
int nrecs;
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
const auto vdate = taxisInqVdate(taxisID1);
const auto vtime = taxisInqVtime(taxisID1);
double jdelta = 0;
if (tsID)
{
const auto juldate0 = julianDateEncode(calendar, vdate0, vtime0);
const auto juldate = julianDateEncode(calendar, vdate, vtime);
jdelta = julianDateToSeconds(julianDateSub(juldate, juldate0));
}
setDateAndTime(params[vartsID], calendar, tsID, vdate0, vtime0, vdate, vtime);
vdate0 = vdate;
vtime0 = vtime;
params[vartsID].data[CTIMESTEP] = tsID + 1;
params[vartsID].data[CDATE] = vdate;
params[vartsID].data[CTIME] = vtime;
params[vartsID].data[CDELTAT] = jdelta;
int year, mon, day;
int hour, minute, second;
cdiDecodeDate(vdate, &year, &mon, &day);
cdiDecodeTime(vtime, &hour, &minute, &second);
params[vartsID].data[CDAY] = day;
params[vartsID].data[CMONTH] = mon;
params[vartsID].data[CYEAR] = year;
params[vartsID].data[CSECOND] = second;
params[vartsID].data[CMINUTE] = minute;
params[vartsID].data[CHOUR] = hour;
taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID);
//for (int varID = 0; varID < nvars1; varID++) printf(">>> %s %d\n", params[varID].name, params[varID].isValid);
for (int varID = 0; varID < nvars1; varID++) params[varID].isValid = true;
for (int varID = 0; varID < nvars1; varID++)
if (tsID == 0 || params[varID].steptype != TIME_CONSTANT) params[varID].nmiss = 0;
......@@ -713,11 +727,9 @@ Expr(void *process)
{
const auto pidx = varIDmap[varID];
if (pidx < nvars1) continue;
const auto ngp = params[pidx].ngp;
const auto nlev = params[pidx].nlev;
params[pidx].nmiss = 0;
varrayFill(ngp * nlev, params[pidx].data, 0.0);
varrayFill(params[pidx].ngp * params[pidx].nlev, params[pidx].data, 0.0);
}
parse_arg.cnparams = vartsID + 1;
......@@ -733,7 +745,7 @@ Expr(void *process)
const auto missval = vlistInqVarMissval(vlistID2, varID);
const auto ngp = params[pidx].ngp;
const int nlev = (int) params[pidx].nlev;
const auto nlev = (int) params[pidx].nlev;
for (int levelID = 0; levelID < nlev; levelID++)
{
const auto offset = ngp * levelID;
......
......@@ -85,7 +85,7 @@ print_location_LL(int operfunc, int vlistID, int varID, int levelID, int gridID,
}