From b4e05f34f7a3bbbaa62e909a2e6b88209763c4d3 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Thu, 13 Mar 2025 12:09:25 +0100
Subject: [PATCH 01/10] config/default: changed netcdf library on levante
---
config/default | 5 ++---
1 file changed, 2 insertions(+), 3 deletions(-)
diff --git a/config/default b/config/default
index b98fb1231..759af6655 100755
--- a/config/default
+++ b/config/default
@@ -369,7 +369,8 @@ case "${HOSTNAME}" in
# NETCDFPATH=/sw/spack-levante/netcdf-c-4.8.1-qk24yp
# HDF5PATH=/sw/spack-levante/hdf5-1.12.1-akf2kp
levante*)
- NETCDFPATH=/sw/spack-levante/netcdf-c-4.9.3-dp4zi4
+ #NETCDFPATH=/sw/spack-levante/netcdf-c-4.9.3-dp4zi4
+ NETCDFPATH=/sw/spack-levante/netcdf-c-4.9.3-wdybfq
UDUNITS2PATH=/sw/spack-levante/udunits-2.2.28-da6pla
FFTW3PATH=/sw/spack-levante/fftw-3.3.10-fnfhvr
ECCODESPATH=/sw/spack-levante/eccodes-2.32.5-ly6tko
@@ -378,7 +379,6 @@ case "${HOSTNAME}" in
#PROJPATH=/sw/spack-levante/proj-8.1.0-i6a6ah
CDOLIBS="--with-eccodes=$ECCODESPATH \
--with-netcdf=$NETCDFPATH \
- --with-hdf5=$HDF5PATH \
--with-udunits2=$UDUNITS2PATH \
--with-fftw3 \
--with-szlib=$SZPATH \
@@ -392,7 +392,6 @@ case "${HOSTNAME}" in
LDFLAGS="$LDFLAGS -L$SZPATH/lib64 -Wl,-rpath,$SZPATH/lib64"
LDFLAGS="$LDFLAGS -L$NETCDFPATH/lib64 -Wl,-rpath,$NETCDFPATH/lib64"
# LDFLAGS="$LDFLAGS -Wl,-rpath,/sw/spack-levante/hdf-4.2.16-2-h6elwi/lib"
- LDFLAGS="$LDFLAGS -Wl,-rpath,$HDF5PATH/lib"
LDFLAGS="$LDFLAGS -Wl,-rpath,$UDUNITS2PATH/lib"
# LDFLAGS="$LDFLAGS -Wl,-rpath,$PROJPATH/lib"
LDFLAGS="$LDFLAGS -L$FFTW3PATH/lib -Wl,-rpath,$FFTW3PATH/lib"
--
GitLab
From de8756c440d7a990fa07d212424d849e839f993b Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Fri, 14 Mar 2025 08:38:49 +0100
Subject: [PATCH 02/10] cat, mergetime: unpack input data
---
ChangeLog | 4 ++++
src/Cat.cc | 1 +
src/Merge.cc | 2 ++
src/Mergetime.cc | 1 +
4 files changed, 8 insertions(+)
diff --git a/ChangeLog b/ChangeLog
index 2d0b7415e..247c0c56f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2025-03-14 Uwe Schulzweida
+
+ * cat, mergetime: unpack input data
+
2025-03-05 Uwe Schulzweida
* Using CDI library version 2.5.1
diff --git a/src/Cat.cc b/src/Cat.cc
index 6fa6ad56d..0c1cfb9ca 100644
--- a/src/Cat.cc
+++ b/src/Cat.cc
@@ -103,6 +103,7 @@ public:
streamMode = StreamMode::CREATE;
streamID2 = cdo_open_write(numFiles);
vlistID2 = vlistDuplicate(vlistID1);
+ vlist_unpack(vlistID2);
taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
diff --git a/src/Merge.cc b/src/Merge.cc
index 0728fed89..1345bf4b8 100644
--- a/src/Merge.cc
+++ b/src/Merge.cc
@@ -189,6 +189,8 @@ public:
vlistMerge(vlistID2, vlistIDs[im]);
}
+ vlist_unpack(vlistID2);
+
int numConstVars = 0;
for (int im = 0; im < numMerge; ++im)
{
diff --git a/src/Mergetime.cc b/src/Mergetime.cc
index a3574f37b..5304e9cc9 100644
--- a/src/Mergetime.cc
+++ b/src/Mergetime.cc
@@ -288,6 +288,7 @@ public:
: ((mapFlag == MapFlag::Intersect) ? streamInfoList[vlistFileIDmin].vlistID
: streamInfoList[vlistFileIDmax].vlistID);
vlistID2 = vlistDuplicate(vlistID1);
+ vlist_unpack(vlistID2);
auto taxisID1 = vlistInqTaxis(si.vlistID);
taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
--
GitLab
From 1f32133cb138db09c139b8488f8e57819a5cc7a5 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Fri, 14 Mar 2025 18:16:23 +0100
Subject: [PATCH 03/10] enable-hirlam-extensions failed since release 2.5.1
(bug fix)
---
ChangeLog | 1 +
configure.ac | 2 +-
libcdi | 2 +-
src/cdo.cc | 4 ----
src/cdo_def_options.cc | 4 ++++
src/vertical_interp.cc | 2 +-
6 files changed, 8 insertions(+), 7 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 247c0c56f..59a1af4eb 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,7 @@
2025-03-14 Uwe Schulzweida
* cat, mergetime: unpack input data
+ * enable-hirlam-extensions failed since release 2.5.1 (bug fix)
2025-03-05 Uwe Schulzweida
diff --git a/configure.ac b/configure.ac
index 075839b1f..20a557a12 100644
--- a/configure.ac
+++ b/configure.ac
@@ -5,7 +5,7 @@
# libtool 2.4.2
AC_PREREQ([2.69])
-AC_INIT([cdo],[2.5.1],[https://mpimet.mpg.de/cdo])
+AC_INIT([cdo],[2.5.2],[https://mpimet.mpg.de/cdo])
AC_DEFINE_UNQUOTED(CDO, ["$PACKAGE_VERSION"], [CDO version])
diff --git a/libcdi b/libcdi
index b9a7b1f01..6729bee6d 160000
--- a/libcdi
+++ b/libcdi
@@ -1 +1 @@
-Subproject commit b9a7b1f01600456c195ea4b1611098e0d39cdfe7
+Subproject commit 6729bee6d0a040c99d0f6aac7d768b1adb5a81ee
diff --git a/src/cdo.cc b/src/cdo.cc
index 042d8bc7b..8e258f153 100644
--- a/src/cdo.cc
+++ b/src/cdo.cc
@@ -53,10 +53,6 @@ cdo_exit(std::string msg = "")
static bool applyDryRun = false;
-#ifdef HIRLAM_EXTENSIONS
-extern "C" void streamGrbDefDataScanningMode(int scanmode);
-#endif
-
static void
cdo_display_syntax_help(std::string const &help, FILE *p_target)
{
diff --git a/src/cdo_def_options.cc b/src/cdo_def_options.cc
index 5884f029c..28d9eb824 100644
--- a/src/cdo_def_options.cc
+++ b/src/cdo_def_options.cc
@@ -19,6 +19,10 @@
#include "cdo_zaxis.h"
#include "chunkspec.h"
+#ifdef HIRLAM_EXTENSIONS
+extern "C" void streamGrbDefDataScanningMode(int scanmode);
+#endif
+
static void
set_chunkspec_parameter(std::string const &argument)
{
diff --git a/src/vertical_interp.cc b/src/vertical_interp.cc
index 13a266fa7..111a21da8 100644
--- a/src/vertical_interp.cc
+++ b/src/vertical_interp.cc
@@ -171,7 +171,7 @@ extrapolate_Z(double pres, double halfPress, double fullPress, double geop, doub
if (tmsl > ztlim && tstar <= ztlim) tmsl = ztlim;
- if (tmsl - tstar < 0.000001 && tstar - tmsl < 0.000001)
+ if ((tmsl - tstar) < 0.000001 && (tstar - tmsl) < 0.000001)
alpha = 0.0;
else if (geop > 0.0001 || geop < -0.0001)
alpha = PlanetRD * (tmsl - tstar) / geop;
--
GitLab
From bf71e9ac0670a23b557f0cf8b4c417b319801945 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Sat, 15 Mar 2025 17:47:20 +0100
Subject: [PATCH 04/10] timpctl: Missing error handling for
CDO_PCTL_NBINS>32768 (bug fix)
---
ChangeLog | 4 +
src/Timpctl.cc | 1 -
src/percentiles_hist.cc | 168 ++++++++++++++++++++--------------------
src/percentiles_hist.h | 42 +++++-----
4 files changed, 109 insertions(+), 106 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 59a1af4eb..449802dd9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2025-03-15 Uwe Schulzweida
+
+ * timpctl: Missing error handling for CDO_PCTL_NBINS>32768 (bug fix)
+
2025-03-14 Uwe Schulzweida
* cat, mergetime: unpack input data
diff --git a/src/Timpctl.cc b/src/Timpctl.cc
index 4cf13f1d4..8dd8fcc61 100644
--- a/src/Timpctl.cc
+++ b/src/Timpctl.cc
@@ -185,7 +185,6 @@ public:
{
field1.init(var1);
cdo_read_field(streamID1, field1);
-
hset.addVarLevelValues(varID, levelID, field1);
}
}
diff --git a/src/percentiles_hist.cc b/src/percentiles_hist.cc
index 8705ee1de..4a5f01151 100644
--- a/src/percentiles_hist.cc
+++ b/src/percentiles_hist.cc
@@ -22,7 +22,7 @@
#define SHR_PTR(p) ((unsigned short *) (p))
static int
-histGetEnvNBins()
+histGetEnvNumBins()
{
constexpr int NBINS_DEFAULT = 101;
constexpr int NBINS_MINIMUM = 11;
@@ -34,47 +34,47 @@ histGetEnvNBins()
template <typename T>
static void
-hist_init_bins(int nbins, T *bins)
+hist_init_bins(int numBins, T *bins)
{
- for (int i = 0; i < nbins; ++i) bins[i] = 0;
+ for (int i = 0; i < numBins; ++i) bins[i] = 0;
}
static void
histDefBounds(HistogramEntry &hist, float a, float b)
{
- assert(hist.nbins > 0);
+ assert(hist.numBins > 0);
hist.nsamp = 0;
hist.min = std::min(a, b);
hist.max = std::max(a, b);
- hist.step = (hist.max - hist.min) / hist.nbins;
+ hist.step = (hist.max - hist.min) / hist.numBins;
if (hist.isUint32)
- hist_init_bins(hist.nbins, INT_PTR(hist.ptr));
+ hist_init_bins(hist.numBins, INT_PTR(hist.ptr));
else
- hist_init_bins(hist.nbins, SHR_PTR(hist.ptr));
+ hist_init_bins(hist.numBins, SHR_PTR(hist.ptr));
}
static inline int
-calc_bin(int nbins, float histMin, float histStep, float value)
+calc_bin(int numBins, float histMin, float histStep, float value)
{
- return (histStep > 0.0f) ? std::min((int) ((value - histMin) / histStep), nbins - 1) : 0;
+ return (histStep > 0.0f) ? std::min((int) ((value - histMin) / histStep), numBins - 1) : 0;
}
template <typename T>
static void
histBinAddValue(HistogramEntry &hist, T *bins, float value)
{
- auto bin = calc_bin(hist.nbins, hist.min, hist.step, value);
- if (bin >= 0 && bin < hist.nbins) bins[bin]++;
+ auto bin = calc_bin(hist.numBins, hist.min, hist.step, value);
+ if (bin >= 0 && bin < hist.numBins) bins[bin]++;
}
template <typename T>
static void
histBinSubValue(HistogramEntry &hist, T *bins, float value)
{
- auto bin = calc_bin(hist.nbins, hist.min, hist.step, value);
- if (bin >= 0 && bin < hist.nbins && bins[bin] > 0) bins[bin]--;
+ auto bin = calc_bin(hist.numBins, hist.min, hist.step, value);
+ if (bin >= 0 && bin < hist.numBins && bins[bin] > 0) bins[bin]--;
}
static void
@@ -88,9 +88,9 @@ histBin(HistogramEntry &hist)
for (int i = 0; i < hist.nsamp; ++i) values[i] = fltptr[i];
if (hist.isUint32)
- hist_init_bins(hist.nbins, INT_PTR(hist.ptr));
+ hist_init_bins(hist.numBins, INT_PTR(hist.ptr));
else
- hist_init_bins(hist.nbins, SHR_PTR(hist.ptr));
+ hist_init_bins(hist.numBins, SHR_PTR(hist.ptr));
if (hist.isUint32)
for (int i = 0; i < hist.nsamp; ++i) histBinAddValue(hist, INT_PTR(hist.ptr), values[i]);
@@ -101,7 +101,7 @@ histBin(HistogramEntry &hist)
static int
histReset(HistogramEntry &hist)
{
- assert(hist.nbins > 0);
+ assert(hist.numBins > 0);
if (hist.nsamp < hist.capacity)
{
@@ -110,9 +110,9 @@ histReset(HistogramEntry &hist)
else
{
if (hist.isUint32)
- for (int i = 0; i < hist.nbins; ++i) INT_PTR(hist.ptr)[i] = 0;
+ for (int i = 0; i < hist.numBins; ++i) INT_PTR(hist.ptr)[i] = 0;
else
- for (int i = 0; i < hist.nbins; ++i) SHR_PTR(hist.ptr)[i] = 0;
+ for (int i = 0; i < hist.numBins; ++i) SHR_PTR(hist.ptr)[i] = 0;
}
hist.nsamp = 0;
@@ -132,7 +132,7 @@ histCheckValue(const HistogramEntry &hist, float &value)
static int
histAddValue(HistogramEntry &hist, float value)
{
- assert(hist.nbins > 0);
+ assert(hist.numBins > 0);
// if (is_equal(hist.min, hist.max)) return 0;
@@ -178,7 +178,7 @@ histRemoveValue(HistogramEntry &hist, float value)
static int
histSubValue(HistogramEntry &hist, float value)
{
- assert(hist.nbins > 0);
+ assert(hist.numBins > 0);
if (is_equal(hist.min, hist.max)) return 0;
@@ -211,7 +211,7 @@ histGetBin(int nbins, double s, const T *ptr)
while (count < s);
assert(i > 0);
- assert(i - 1 < nbins);
+ assert((i - 1) < nbins);
assert(ptr[i - 1] > 0);
double t = (count - s) / ptr[i - 1];
@@ -223,7 +223,7 @@ static double
histGetPercentile(const HistogramEntry &hist, double p)
{
assert(hist.nsamp > 0);
- assert(hist.nbins > 0);
+ assert(hist.numBins > 0);
assert(p >= 0.0);
assert(p <= 100.0);
@@ -233,12 +233,12 @@ histGetPercentile(const HistogramEntry &hist, double p)
if (lprint && Options::cdoVerbose)
{
lprint = false;
- cdo_print("Using percentile method: histogram with %d bins", hist.nbins);
+ cdo_print("Using percentile method: histogram with %d bins", hist.numBins);
}
double s = hist.nsamp * (p / 100.0);
- auto bin = hist.isUint32 ? histGetBin(hist.nbins, s, INT_PTR(hist.ptr)) : histGetBin(hist.nbins, s, SHR_PTR(hist.ptr));
+ auto bin = hist.isUint32 ? histGetBin(hist.numBins, s, INT_PTR(hist.ptr)) : histGetBin(hist.numBins, s, SHR_PTR(hist.ptr));
// assert(hist.step > 0.0f);
@@ -248,37 +248,37 @@ histGetPercentile(const HistogramEntry &hist, double p)
}
void
-HistogramSet::createVarLevels(int varID, int nlevels, size_t nhists)
+HistogramSet::createVarLevels(int varID, int numLevels, size_t numHists)
{
- auto nbins = histGetEnvNBins();
+ auto numBins = histGetEnvNumBins();
- assert(nbins > 0);
- assert(nlevels > 0);
- assert(nhists > 0);
+ assert(numBins > 0);
+ assert(numLevels > 0);
+ assert(numHists > 0);
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- this->var_nlevels[varID] = nlevels;
- this->var_nhists[varID] = nhists;
- this->histograms[varID].resize(nlevels);
+ this->var_numLevels[varID] = numLevels;
+ this->var_numHists[varID] = numHists;
+ this->histograms[varID].resize(numLevels);
- for (int levelID = 0; levelID < nlevels; ++levelID)
+ for (int levelID = 0; levelID < numLevels; ++levelID)
{
- this->histograms[varID][levelID].resize(nhists);
+ this->histograms[varID][levelID].resize(numHists);
auto &hists = this->histograms[varID][levelID];
- for (size_t histID = 0; histID < nhists; histID++)
+ for (size_t histID = 0; histID < numHists; histID++)
{
hists[histID].min = 0.0f;
hists[histID].max = 0.0f;
hists[histID].step = 0.0f;
- hists[histID].nbins = nbins;
+ hists[histID].numBins = numBins;
hists[histID].nsamp = 0;
- auto isUint32 = (this->nsteps >= USHRT_MAX);
+ auto isUint32 = (this->numSteps >= USHRT_MAX);
auto vsize = isUint32 ? sizeof(unsigned int) : sizeof(unsigned short);
hists[histID].isUint32 = isUint32;
- hists[histID].capacity = FLT_CAPACITY(nbins, vsize);
- hists[histID].ptr = std::malloc(nbins * vsize);
+ hists[histID].capacity = FLT_CAPACITY(numBins, vsize);
+ hists[histID].ptr = std::malloc(numBins * vsize);
if (hists[histID].ptr == nullptr) cdo_abort("Not enough memory (%s)", __func__);
}
}
@@ -286,12 +286,12 @@ HistogramSet::createVarLevels(int varID, int nlevels, size_t nhists)
template <typename T1, typename T2>
static void
-def_bounds(size_t nhists, std::vector<HistogramEntry> &hists, const Varray<T1> &v1, const Varray<T2> &v2, float mv1, float mv2)
+def_bounds(size_t numHists, std::vector<HistogramEntry> &hists, Varray<T1> const &v1, Varray<T2> const &v2, float mv1, float mv2)
{
assert(!v1.empty());
assert(!v2.empty());
- for (size_t i = 0; i < nhists; ++i)
+ for (size_t i = 0; i < numHists; ++i)
{
float a = v1[i];
float b = v2[i];
@@ -306,26 +306,26 @@ def_bounds(size_t nhists, std::vector<HistogramEntry> &hists, const Varray<T1> &
void
HistogramSet::defVarLevelBounds(int varID, int levelID, Field const &field1, Field const &field2)
{
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- auto nlevels = this->var_nlevels[varID];
- if (levelID < 0 || levelID >= nlevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
+ auto numLevels = this->var_numLevels[varID];
+ if (levelID < 0 || levelID >= numLevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
- auto nhists = this->var_nhists[varID];
- if (nhists != gridInqSize(field1.grid) || nhists != gridInqSize(field2.grid)) cdo_abort("Grids are different (%s)", __func__);
+ auto numHists = this->var_numHists[varID];
+ if (numHists != gridInqSize(field1.grid) || numHists != gridInqSize(field2.grid)) cdo_abort("Grids are different (%s)", __func__);
auto &hists = this->histograms[varID][levelID];
float mv1 = (float) field1.missval;
float mv2 = (float) field2.missval;
- auto func = [&](auto const &v1, auto const &v2) { def_bounds(nhists, hists, v1, v2, mv1, mv2); };
+ auto func = [&](auto const &v1, auto const &v2) { def_bounds(numHists, hists, v1, v2, mv1, mv2); };
field_operation2(func, field1, field2);
}
template <typename T>
static int
-histAddVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<T> const &v, size_t numMissVals, double missval)
+histAddVarLevelValues(size_t numHists, std::vector<HistogramEntry> &hists, Varray<T> const &v, size_t numMissVals, double missval)
{
T mv = missval;
assert(!v.empty());
@@ -333,12 +333,12 @@ histAddVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<
int nign = 0;
if (numMissVals)
{
- for (size_t i = 0; i < nhists; ++i)
+ for (size_t i = 0; i < numHists; ++i)
if (fp_is_not_equal(v[i], mv)) nign += histAddValue(hists[i], v[i]);
}
else
{
- for (size_t i = 0; i < nhists; ++i) nign += histAddValue(hists[i], v[i]);
+ for (size_t i = 0; i < numHists; ++i) nign += histAddValue(hists[i], v[i]);
}
return nign;
@@ -346,7 +346,7 @@ histAddVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<
template <typename T>
static int
-histSubVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<T> const &v, size_t numMissVals, double missval)
+histSubVarLevelValues(size_t numHists, std::vector<HistogramEntry> &hists, Varray<T> const &v, size_t numMissVals, double missval)
{
T mv = missval;
assert(!v.empty());
@@ -354,12 +354,12 @@ histSubVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<
int nign = 0;
if (numMissVals)
{
- for (size_t i = 0; i < nhists; ++i)
+ for (size_t i = 0; i < numHists; ++i)
if (fp_is_not_equal(v[i], mv)) nign += histSubValue(hists[i], v[i]);
}
else
{
- for (size_t i = 0; i < nhists; ++i) nign += histSubValue(hists[i], v[i]);
+ for (size_t i = 0; i < numHists; ++i) nign += histSubValue(hists[i], v[i]);
}
return nign;
@@ -368,21 +368,21 @@ histSubVarLevelValues(size_t nhists, std::vector<HistogramEntry> &hists, Varray<
int
HistogramSet::addVarLevelValues(int varID, int levelID, Field const &field)
{
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- auto nlevels = this->var_nlevels[varID];
- if (levelID < 0 || levelID >= nlevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
+ auto numLevels = this->var_numLevels[varID];
+ if (levelID < 0 || levelID >= numLevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
- auto nhists = this->var_nhists[varID];
- if (nhists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
+ auto numHists = this->var_numHists[varID];
+ if (numHists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
auto &hists = this->histograms[varID][levelID];
- auto func = [&](auto &v, auto numMissVals, double mv) { return histAddVarLevelValues(nhists, hists, v, numMissVals, mv); };
+ auto func = [&](auto &v, auto numMissVals, double mv) { return histAddVarLevelValues(numHists, hists, v, numMissVals, mv); };
auto nign = field_operation(func, field, field.numMissVals, field.missval);
if (nign)
{
- cdo_warning("%d out of %d grid values are out of bounds and have been ignored (%s)", nign, nhists, __func__);
+ cdo_warning("%d out of %d grid values are out of bounds and have been ignored (%s)", nign, numHists, __func__);
return 1;
}
@@ -392,21 +392,21 @@ HistogramSet::addVarLevelValues(int varID, int levelID, Field const &field)
int
HistogramSet::subVarLevelValues(int varID, int levelID, Field const &field)
{
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- auto nlevels = this->var_nlevels[varID];
- if (levelID < 0 || levelID >= nlevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
+ auto numLevels = this->var_numLevels[varID];
+ if (levelID < 0 || levelID >= numLevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
- auto nhists = this->var_nhists[varID];
- if (nhists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
+ auto numHists = this->var_numHists[varID];
+ if (numHists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
auto &hists = this->histograms[varID][levelID];
- auto func = [&](auto &v, auto numMissVals, double mv) { return histSubVarLevelValues(nhists, hists, v, numMissVals, mv); };
+ auto func = [&](auto &v, auto numMissVals, double mv) { return histSubVarLevelValues(numHists, hists, v, numMissVals, mv); };
auto nign = field_operation(func, field, field.numMissVals, field.missval);
if (nign)
{
- cdo_warning("%d out of %d grid values are out of bounds and have been ignored (%s)", nign, nhists, __func__);
+ cdo_warning("%d out of %d grid values are out of bounds and have been ignored (%s)", nign, numHists, __func__);
return 1;
}
@@ -416,33 +416,33 @@ HistogramSet::subVarLevelValues(int varID, int levelID, Field const &field)
void
HistogramSet::reset(int varID, int levelID)
{
- assert(nvars > 0);
+ assert(numVars > 0);
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- auto nlevels = this->var_nlevels[varID];
- assert(nlevels > 0);
+ auto numLevels = this->var_numLevels[varID];
+ assert(numLevels > 0);
- if (levelID < 0 || levelID >= nlevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
+ if (levelID < 0 || levelID >= numLevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
- auto nhists = this->var_nhists[varID];
- assert(nhists > 0);
+ auto numHists = this->var_numHists[varID];
+ assert(numHists > 0);
auto &hists = this->histograms[varID][levelID];
- for (size_t i = 0; i < nhists; ++i) histReset(hists[i]);
+ for (size_t i = 0; i < numHists; ++i) histReset(hists[i]);
}
*/
template <typename T>
static size_t
-calcPercentile(size_t nhists, const std::vector<HistogramEntry> &hists, double p, Varray<T> &v, double missval)
+calcPercentile(size_t numHists, const std::vector<HistogramEntry> &hists, double p, Varray<T> &v, double missval)
{
T mv = missval;
assert(!v.empty());
size_t numMissVals = 0;
- for (size_t i = 0; i < nhists; ++i)
+ for (size_t i = 0; i < numHists; ++i)
{
if (hists[i].nsamp) { v[i] = histGetPercentile(hists[i], p); }
else
@@ -458,16 +458,16 @@ calcPercentile(size_t nhists, const std::vector<HistogramEntry> &hists, double p
void
HistogramSet::getVarLevelPercentiles(Field &field, int varID, int levelID, double p)
{
- if (varID < 0 || varID >= nvars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
+ if (varID < 0 || varID >= numVars) cdo_abort("Illegal argument: varID %d is undefined (%s)", varID, __func__);
- auto nlevels = this->var_nlevels[varID];
- if (levelID < 0 || levelID >= nlevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
+ auto numLevels = this->var_numLevels[varID];
+ if (levelID < 0 || levelID >= numLevels) cdo_abort("Illegal argument: levelID %d is undefined (%s)", levelID, __func__);
- auto nhists = this->var_nhists[varID];
- if (nhists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
+ auto numHists = this->var_numHists[varID];
+ if (numHists != gridInqSize(field.grid)) cdo_abort("Grids are different (%s)", __func__);
auto const &hists = this->histograms[varID][levelID];
- auto func = [&](auto &v, double mv) { return calcPercentile(nhists, hists, p, v, mv); };
+ auto func = [&](auto &v, double mv) { return calcPercentile(numHists, hists, p, v, mv); };
field.numMissVals = field_operation(func, field, field.missval);
}
diff --git a/src/percentiles_hist.h b/src/percentiles_hist.h
index 4f9a444c5..14dc48025 100644
--- a/src/percentiles_hist.h
+++ b/src/percentiles_hist.h
@@ -14,7 +14,7 @@ struct HistogramEntry
float step = 0.0f;
int nsamp = 0;
int capacity = 0;
- short nbins = 0;
+ int numBins = 0;
bool isUint32 = false;
};
@@ -27,57 +27,57 @@ HistogramSet
// clang-format on
{
private:
- int nvars = 0;
- int nsteps = 0;
- std::vector<int> var_nlevels;
- std::vector<size_t> var_nhists;
+ int numVars = 0;
+ int numSteps = 0;
+ std::vector<int> var_numLevels;
+ std::vector<size_t> var_numHists;
std::vector<std::vector<std::vector<HistogramEntry>>> histograms;
void
init()
{
- var_nlevels.resize(nvars, 0);
- var_nhists.resize(nvars, 0);
- histograms.resize(nvars);
+ var_numLevels.resize(numVars, 0);
+ var_numHists.resize(numVars, 0);
+ histograms.resize(numVars);
}
public:
HistogramSet() {}
- explicit HistogramSet(int _nvars) : nvars(_nvars)
+ explicit HistogramSet(int _numVars) : numVars(_numVars)
{
- assert(nvars > 0);
+ assert(numVars > 0);
init();
}
- HistogramSet(int _nvars, int _nsteps) : nvars(_nvars), nsteps(_nsteps)
+ HistogramSet(int _numVars, int _numSteps) : numVars(_numVars), numSteps(_numSteps)
{
- assert(nvars > 0);
+ assert(numVars > 0);
init();
}
void
- create(int _nvars, int _nsteps = 0)
+ create(int _numVars, int _numSteps = 0)
{
- nvars = _nvars;
- nsteps = _nsteps;
- assert(nvars > 0);
+ numVars = _numVars;
+ numSteps = _numSteps;
+ assert(numVars > 0);
init();
}
~HistogramSet()
{
- for (auto varID = nvars; varID-- > 0;)
+ for (auto varID = numVars; varID-- > 0;)
{
- auto nhists = this->var_nhists[varID];
- for (auto levelID = this->var_nlevels[varID]; levelID-- > 0;)
+ auto numHists = this->var_numHists[varID];
+ for (auto levelID = this->var_numLevels[varID]; levelID-- > 0;)
{
- for (auto histID = nhists; histID-- > 0;) std::free(this->histograms[varID][levelID][histID].ptr);
+ for (auto histID = numHists; histID-- > 0;) std::free(this->histograms[varID][levelID][histID].ptr);
}
}
}
- void createVarLevels(int varID, int nlevels, size_t nhists);
+ void createVarLevels(int varID, int numLevels, size_t numHists);
void defVarLevelBounds(int varID, int levelID, Field const &field1, Field const &field2);
int addVarLevelValues(int varID, int levelID, Field const &field);
int subVarLevelValues(int varID, int levelID, Field const &field);
--
GitLab
From 09e1eea45ea5717bd4ea0faa764740047c65cda9 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Mon, 17 Mar 2025 13:51:51 +0100
Subject: [PATCH 05/10] Detrend: failed if missing_value is between 0 and
numSteps (bug fix)
---
ChangeLog | 2 +-
NEWS | 5 +++++
src/Detrend.cc | 8 ++++----
src/Trend.cc | 4 ++--
src/field_trend.cc | 6 +++---
5 files changed, 15 insertions(+), 10 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 449802dd9..211f29c28 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,6 @@
2025-03-15 Uwe Schulzweida
- * timpctl: Missing error handling for CDO_PCTL_NBINS>32768 (bug fix)
+ * timpctl: Short integer overflow for CDO_PCTL_NBINS>32768 (bug fix)
2025-03-14 Uwe Schulzweida
diff --git a/NEWS b/NEWS
index ef15a1a5a..c27aea57a 100644
--- a/NEWS
+++ b/NEWS
@@ -3,11 +3,16 @@ CDO NEWS
Improvement
+Version 2.5.2 (15 Mar 2025):
+
New features:
* option --chunkspec to define chunkSize of t and z dimension
New operators:
* air_density
* setchunkspec - Specify chunking
+ Fixed bugs:
+ * timpctl: Short integer overflow for CDO_PCTL_NBINS>32768
+ * enable-hirlam-extensions failed since release 2.5.1
Version 2.5.1 (5 Mar 2025):
diff --git a/src/Detrend.cc b/src/Detrend.cc
index c6043e01a..89a1d300f 100644
--- a/src/Detrend.cc
+++ b/src/Detrend.cc
@@ -129,11 +129,11 @@ public:
auto const &zn = work[4][varID][levelID].vec_d;
auto trend_kernel = [&](auto i, auto is_EQ) {
- auto temp1 = SUBM(sumjx[i], DIVM(MULM(sumj[i], sumx[i]), zn[i]));
- auto temp2 = SUBM(sumjj[i], DIVM(MULM(sumj[i], sumj[i]), zn[i]));
+ auto temp1 = SUBM(sumjx[i], DIVX(MULM(sumj[i], sumx[i]), zn[i]));
+ auto temp2 = SUBM(sumjj[i], DIVX(MULM(sumj[i], sumj[i]), zn[i]));
auto temp3 = DIVM(temp1, temp2);
- paramA[i] = SUBM(DIVM(sumx[i], zn[i]), MULM(DIVM(sumj[i], zn[i]), temp3));
+ paramA[i] = SUBM(DIVX(sumx[i], zn[i]), MULM(DIVX(sumj[i], zn[i]), temp3));
paramB[i] = temp3;
};
@@ -146,7 +146,7 @@ public:
}
static void
- vars_sub_trend(FieldVector3D &work, FieldVector2D &varsData, const VarList &varList, double zj)
+ vars_sub_trend(FieldVector3D &work, FieldVector2D &varsData, VarList const &varList, double zj)
{
auto numVars = varList.numVars();
for (int varID = 0; varID < numVars; ++varID)
diff --git a/src/Trend.cc b/src/Trend.cc
index 9c9099163..ea2badeb0 100644
--- a/src/Trend.cc
+++ b/src/Trend.cc
@@ -123,7 +123,7 @@ public:
}
void
- write_output(const FieldVector3D &work)
+ write_output(FieldVector3D const &work)
{
Field field2, field3;
@@ -198,7 +198,7 @@ public:
}
static void
- fields_calc_trend_sum(FieldVector3D &work, const FieldVector2D &fields2D, const std::vector<FieldInfo> &fieldInfoList,
+ fields_calc_trend_sum(FieldVector3D &work, FieldVector2D const &fields2D, std::vector<FieldInfo> const &fieldInfoList,
double zj) noexcept
{
for (auto const &fieldInfo : fieldInfoList)
diff --git a/src/field_trend.cc b/src/field_trend.cc
index 9155952d0..a1da6d78d 100644
--- a/src/field_trend.cc
+++ b/src/field_trend.cc
@@ -99,11 +99,11 @@ calc_trend_param(const FieldVector3D &work, Field ¶mA, Field ¶mB, int va
auto const &zn = work[4][varID][levelID].vec_d;
auto trend_kernel = [&](auto i, auto is_EQ) {
- auto temp1 = SUBM(sumjx[i], DIVM(MULM(sumj[i], sumx[i]), zn[i]));
- auto temp2 = SUBM(sumjj[i], DIVM(MULM(sumj[i], sumj[i]), zn[i]));
+ auto temp1 = SUBM(sumjx[i], DIVX(MULM(sumj[i], sumx[i]), zn[i]));
+ auto temp2 = SUBM(sumjj[i], DIVX(MULM(sumj[i], sumj[i]), zn[i]));
auto temp3 = DIVM(temp1, temp2);
- paramA.vec_d[i] = SUBM(DIVM(sumx[i], zn[i]), MULM(DIVM(sumj[i], zn[i]), temp3));
+ paramA.vec_d[i] = SUBM(DIVX(sumx[i], zn[i]), MULM(DIVX(sumj[i], zn[i]), temp3));
paramB.vec_d[i] = temp3;
};
--
GitLab
From 4bc11fcc3b2a20846286bc0ed873026bafd84158 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Mon, 17 Mar 2025 14:24:36 +0100
Subject: [PATCH 06/10] Replace DIVM() by DIVX()
---
src/Timstat2.cc | 10 +++++-----
src/Timstat3.cc | 22 +++++++++++-----------
src/fieldc.cc | 4 ++--
3 files changed, 18 insertions(+), 18 deletions(-)
diff --git a/src/Timstat2.cc b/src/Timstat2.cc
index fc688dfe2..873933689 100644
--- a/src/Timstat2.cc
+++ b/src/Timstat2.cc
@@ -99,11 +99,11 @@ correlation(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varr
{
double dnvals = nvals;
auto temp0 = MULM(work[0][i], work[1][i]);
- auto temp1 = SUBM(work[4][i], DIVM(temp0, dnvals));
+ auto temp1 = SUBM(work[4][i], DIVX(temp0, dnvals));
auto temp2 = MULM(work[0][i], work[0][i]);
auto temp3 = MULM(work[1][i], work[1][i]);
- auto temp4 = SUBM(work[2][i], DIVM(temp2, dnvals));
- auto temp5 = SUBM(work[3][i], DIVM(temp3, dnvals));
+ auto temp4 = SUBM(work[2][i], DIVX(temp2, dnvals));
+ auto temp5 = SUBM(work[3][i], DIVX(temp3, dnvals));
auto temp6 = MULM(temp4, temp5);
cor = DIVM(temp1, SQRTM(temp6));
@@ -194,8 +194,8 @@ covariance(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varra
if (nvals > 0)
{
double dnvals = nvals;
- auto temp = DIVM(MULM(work[0][i], work[1][i]), dnvals * dnvals);
- covar = SUBM(DIVM(work[2][i], dnvals), temp);
+ auto temp = DIVX(MULM(work[0][i], work[1][i]), dnvals * dnvals);
+ covar = SUBM(DIVX(work[2][i], dnvals), temp);
if (fp_is_equal(covar, missval)) numMissVals++;
}
else
diff --git a/src/Timstat3.cc b/src/Timstat3.cc
index 86e954305..819978fae 100644
--- a/src/Timstat3.cc
+++ b/src/Timstat3.cc
@@ -29,8 +29,8 @@ varquot2test(double rconst, double risk, size_t gridsize, double missval, Varray
auto missval2 = missval;
auto varquot2test_kernel = [&](auto i, auto is_EQ) {
- auto temp0 = DIVM(MULM(work[0][i], work[0][i]), work[2][i]);
- auto temp1 = DIVM(MULM(work[3][i], work[3][i]), work[5][i]);
+ auto temp0 = DIVX(MULM(work[0][i], work[0][i]), work[2][i]);
+ auto temp1 = DIVX(MULM(work[3][i], work[3][i]), work[5][i]);
auto temp2 = SUBM(work[1][i], temp0);
auto temp3 = SUBM(work[4][i], temp1);
auto statistic = DIVM(temp2, ADDM(temp2, MULM(rconst, temp3)));
@@ -64,23 +64,23 @@ meandiff2test(double rconst, double risk, size_t gridsize, double missval, Varra
auto meandiff2test_kernel = [&](auto i, auto is_EQ) {
double temp0 = 0.0;
double degOfFreedom = -NIN;
- auto tmp = DIVM(MULM(work[0][i], work[0][i]), work[2][i]);
- temp0 = ADDM(temp0, DIVM(SUBM(work[1][i], tmp), varFactor[0]));
+ auto tmp = DIVX(MULM(work[0][i], work[0][i]), work[2][i]);
+ temp0 = ADDM(temp0, DIVX(SUBM(work[1][i], tmp), varFactor[0]));
degOfFreedom = ADDM(degOfFreedom, work[2][i]);
- tmp = DIVM(MULM(work[3][i], work[3][i]), work[5][i]);
- temp0 = ADDM(temp0, DIVM(SUBM(work[4][i], tmp), varFactor[1]));
+ tmp = DIVX(MULM(work[3][i], work[3][i]), work[5][i]);
+ temp0 = ADDM(temp0, DIVX(SUBM(work[4][i], tmp), varFactor[1]));
degOfFreedom = ADDM(degOfFreedom, work[5][i]);
if (fp_is_not_equal(temp0, missval1) && temp0 < 0) temp0 = 0; // This is possible because of rounding errors
- auto stddevEstimator = SQRTM(DIVM(temp0, degOfFreedom));
+ auto stddevEstimator = SQRTM(DIVX(temp0, degOfFreedom));
auto meanEstimator = -rconst;
- meanEstimator = ADDM(meanEstimator, MULM(meanFactor[0], DIVM(work[0][i], work[2][i])));
- meanEstimator = ADDM(meanEstimator, MULM(meanFactor[1], DIVM(work[3][i], work[5][i])));
+ meanEstimator = ADDM(meanEstimator, MULM(meanFactor[0], DIVX(work[0][i], work[2][i])));
+ meanEstimator = ADDM(meanEstimator, MULM(meanFactor[1], DIVX(work[3][i], work[5][i])));
double temp1 = 0.0;
- temp1 = ADDM(temp1, DIVM(factor1, work[2][i]));
- temp1 = ADDM(temp1, DIVM(factor2, work[5][i]));
+ temp1 = ADDM(temp1, DIVX(factor1, work[2][i]));
+ temp1 = ADDM(temp1, DIVX(factor2, work[5][i]));
auto norm = SQRTM(temp1);
auto temp2 = DIVM(DIVM(meanEstimator, norm), stddevEstimator);
diff --git a/src/fieldc.cc b/src/fieldc.cc
index 6c9e8d093..995219107 100644
--- a/src/fieldc.cc
+++ b/src/fieldc.cc
@@ -40,11 +40,11 @@ fieldc_div(Varray<T> &v, size_t len, size_t numMissVals, double missval, double
{
auto is_EQ = fp_is_equal;
T missval1 = missval;
- T missval2 = missval;
+ // T missval2 = missval;
if (numMissVals || is_equal(rconst, 0))
{
- for (size_t i = 0; i < len; ++i) v[i] = DIVM(v[i], rconst);
+ for (size_t i = 0; i < len; ++i) v[i] = DIVX(v[i], rconst);
if (is_equal(rconst, 0)) numMissVals = len;
}
--
GitLab
From 01671651e584b67c2cfa185cadeb9735a5bee2fc Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Mon, 17 Mar 2025 20:26:56 +0100
Subject: [PATCH 07/10] Replaced DIVX() by DIVMX()
---
src/Detrend.cc | 6 +++---
src/Timstat2.cc | 10 +++++-----
src/Timstat3.cc | 22 +++++++++++-----------
src/arithmetic.h | 2 ++
src/field_trend.cc | 6 +++---
src/fieldc.cc | 2 +-
6 files changed, 25 insertions(+), 23 deletions(-)
diff --git a/src/Detrend.cc b/src/Detrend.cc
index 89a1d300f..ce349c32c 100644
--- a/src/Detrend.cc
+++ b/src/Detrend.cc
@@ -129,11 +129,11 @@ public:
auto const &zn = work[4][varID][levelID].vec_d;
auto trend_kernel = [&](auto i, auto is_EQ) {
- auto temp1 = SUBM(sumjx[i], DIVX(MULM(sumj[i], sumx[i]), zn[i]));
- auto temp2 = SUBM(sumjj[i], DIVX(MULM(sumj[i], sumj[i]), zn[i]));
+ auto temp1 = SUBM(sumjx[i], DIVMX(MULM(sumj[i], sumx[i]), zn[i]));
+ auto temp2 = SUBM(sumjj[i], DIVMX(MULM(sumj[i], sumj[i]), zn[i]));
auto temp3 = DIVM(temp1, temp2);
- paramA[i] = SUBM(DIVX(sumx[i], zn[i]), MULM(DIVX(sumj[i], zn[i]), temp3));
+ paramA[i] = SUBM(DIVMX(sumx[i], zn[i]), MULM(DIVMX(sumj[i], zn[i]), temp3));
paramB[i] = temp3;
};
diff --git a/src/Timstat2.cc b/src/Timstat2.cc
index 873933689..fb6a10024 100644
--- a/src/Timstat2.cc
+++ b/src/Timstat2.cc
@@ -99,11 +99,11 @@ correlation(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varr
{
double dnvals = nvals;
auto temp0 = MULM(work[0][i], work[1][i]);
- auto temp1 = SUBM(work[4][i], DIVX(temp0, dnvals));
+ auto temp1 = SUBM(work[4][i], DIVMX(temp0, dnvals));
auto temp2 = MULM(work[0][i], work[0][i]);
auto temp3 = MULM(work[1][i], work[1][i]);
- auto temp4 = SUBM(work[2][i], DIVX(temp2, dnvals));
- auto temp5 = SUBM(work[3][i], DIVX(temp3, dnvals));
+ auto temp4 = SUBM(work[2][i], DIVMX(temp2, dnvals));
+ auto temp5 = SUBM(work[3][i], DIVMX(temp3, dnvals));
auto temp6 = MULM(temp4, temp5);
cor = DIVM(temp1, SQRTM(temp6));
@@ -194,8 +194,8 @@ covariance(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varra
if (nvals > 0)
{
double dnvals = nvals;
- auto temp = DIVX(MULM(work[0][i], work[1][i]), dnvals * dnvals);
- covar = SUBM(DIVX(work[2][i], dnvals), temp);
+ auto temp = DIVMX(MULM(work[0][i], work[1][i]), dnvals * dnvals);
+ covar = SUBM(DIVMX(work[2][i], dnvals), temp);
if (fp_is_equal(covar, missval)) numMissVals++;
}
else
diff --git a/src/Timstat3.cc b/src/Timstat3.cc
index 819978fae..ca8869b70 100644
--- a/src/Timstat3.cc
+++ b/src/Timstat3.cc
@@ -29,8 +29,8 @@ varquot2test(double rconst, double risk, size_t gridsize, double missval, Varray
auto missval2 = missval;
auto varquot2test_kernel = [&](auto i, auto is_EQ) {
- auto temp0 = DIVX(MULM(work[0][i], work[0][i]), work[2][i]);
- auto temp1 = DIVX(MULM(work[3][i], work[3][i]), work[5][i]);
+ auto temp0 = DIVMX(MULM(work[0][i], work[0][i]), work[2][i]);
+ auto temp1 = DIVMX(MULM(work[3][i], work[3][i]), work[5][i]);
auto temp2 = SUBM(work[1][i], temp0);
auto temp3 = SUBM(work[4][i], temp1);
auto statistic = DIVM(temp2, ADDM(temp2, MULM(rconst, temp3)));
@@ -64,23 +64,23 @@ meandiff2test(double rconst, double risk, size_t gridsize, double missval, Varra
auto meandiff2test_kernel = [&](auto i, auto is_EQ) {
double temp0 = 0.0;
double degOfFreedom = -NIN;
- auto tmp = DIVX(MULM(work[0][i], work[0][i]), work[2][i]);
- temp0 = ADDM(temp0, DIVX(SUBM(work[1][i], tmp), varFactor[0]));
+ auto tmp = DIVMX(MULM(work[0][i], work[0][i]), work[2][i]);
+ temp0 = ADDM(temp0, DIVMX(SUBM(work[1][i], tmp), varFactor[0]));
degOfFreedom = ADDM(degOfFreedom, work[2][i]);
- tmp = DIVX(MULM(work[3][i], work[3][i]), work[5][i]);
- temp0 = ADDM(temp0, DIVX(SUBM(work[4][i], tmp), varFactor[1]));
+ tmp = DIVMX(MULM(work[3][i], work[3][i]), work[5][i]);
+ temp0 = ADDM(temp0, DIVMX(SUBM(work[4][i], tmp), varFactor[1]));
degOfFreedom = ADDM(degOfFreedom, work[5][i]);
if (fp_is_not_equal(temp0, missval1) && temp0 < 0) temp0 = 0; // This is possible because of rounding errors
- auto stddevEstimator = SQRTM(DIVX(temp0, degOfFreedom));
+ auto stddevEstimator = SQRTM(DIVMX(temp0, degOfFreedom));
auto meanEstimator = -rconst;
- meanEstimator = ADDM(meanEstimator, MULM(meanFactor[0], DIVX(work[0][i], work[2][i])));
- meanEstimator = ADDM(meanEstimator, MULM(meanFactor[1], DIVX(work[3][i], work[5][i])));
+ meanEstimator = ADDM(meanEstimator, MULM(meanFactor[0], DIVMX(work[0][i], work[2][i])));
+ meanEstimator = ADDM(meanEstimator, MULM(meanFactor[1], DIVMX(work[3][i], work[5][i])));
double temp1 = 0.0;
- temp1 = ADDM(temp1, DIVX(factor1, work[2][i]));
- temp1 = ADDM(temp1, DIVX(factor2, work[5][i]));
+ temp1 = ADDM(temp1, DIVMX(factor1, work[2][i]));
+ temp1 = ADDM(temp1, DIVMX(factor2, work[5][i]));
auto norm = SQRTM(temp1);
auto temp2 = DIVM(DIVM(meanEstimator, norm), stddevEstimator);
diff --git a/src/arithmetic.h b/src/arithmetic.h
index 600117b93..2b4e9f3d4 100644
--- a/src/arithmetic.h
+++ b/src/arithmetic.h
@@ -20,6 +20,7 @@ const auto mulm = [](auto x, auto y, auto mvx, auto mvy, auto is_EQ) { return (i
const auto divm = [](auto x, auto y, auto mvx, auto mvy, auto is_EQ) { return (is_EQ(x, mvx) || is_EQ(y, mvy) || is_EQ(y, 0)) ? mvx : x / y; };
const auto powm = [](auto x, auto y, auto mvx, auto mvy, auto is_EQ) { return (is_EQ(x, mvx) || is_EQ(y, mvy)) ? mvx : std::pow(x, y); };
const auto sqrtm = [](auto x, auto mvx, auto is_EQ) { return (is_EQ(x, mvx) || x < 0) ? mvx : std::sqrt(x); };
+const auto divmx = [](auto x, auto y, auto mvx, auto is_EQ) { return (is_EQ(x, mvx) || is_EQ(y, 0)) ? mvx : x / y; };
const auto divx = [](auto x, auto y, auto mvx, auto is_EQ) { return is_EQ(y, 0) ? mvx : x / y; };
}
@@ -36,6 +37,7 @@ const auto SQRT = arithmetic::sqrt;
#define DIVM(x, y) arithmetic::divm(x, y, missval1, missval2, is_EQ)
#define POWM(x, y) arithmetic::powm(x, y, missval1, missval2, is_EQ)
#define SQRTM(x) arithmetic::sqrtm(x, missval1, is_EQ)
+#define DIVMX(x, y) arithmetic::divmx(x, y, missval1, is_EQ)
#define DIVX(x, y) arithmetic::divx(x, y, missval1, is_EQ)
// clang-format on
diff --git a/src/field_trend.cc b/src/field_trend.cc
index a1da6d78d..4800e9c76 100644
--- a/src/field_trend.cc
+++ b/src/field_trend.cc
@@ -99,11 +99,11 @@ calc_trend_param(const FieldVector3D &work, Field ¶mA, Field ¶mB, int va
auto const &zn = work[4][varID][levelID].vec_d;
auto trend_kernel = [&](auto i, auto is_EQ) {
- auto temp1 = SUBM(sumjx[i], DIVX(MULM(sumj[i], sumx[i]), zn[i]));
- auto temp2 = SUBM(sumjj[i], DIVX(MULM(sumj[i], sumj[i]), zn[i]));
+ auto temp1 = SUBM(sumjx[i], DIVMX(MULM(sumj[i], sumx[i]), zn[i]));
+ auto temp2 = SUBM(sumjj[i], DIVMX(MULM(sumj[i], sumj[i]), zn[i]));
auto temp3 = DIVM(temp1, temp2);
- paramA.vec_d[i] = SUBM(DIVX(sumx[i], zn[i]), MULM(DIVX(sumj[i], zn[i]), temp3));
+ paramA.vec_d[i] = SUBM(DIVMX(sumx[i], zn[i]), MULM(DIVMX(sumj[i], zn[i]), temp3));
paramB.vec_d[i] = temp3;
};
diff --git a/src/fieldc.cc b/src/fieldc.cc
index 995219107..ef7008490 100644
--- a/src/fieldc.cc
+++ b/src/fieldc.cc
@@ -44,7 +44,7 @@ fieldc_div(Varray<T> &v, size_t len, size_t numMissVals, double missval, double
if (numMissVals || is_equal(rconst, 0))
{
- for (size_t i = 0; i < len; ++i) v[i] = DIVX(v[i], rconst);
+ for (size_t i = 0; i < len; ++i) v[i] = DIVMX(v[i], rconst);
if (is_equal(rconst, 0)) numMissVals = len;
}
--
GitLab
From 3319914180020899cba3084ef5c3647650eb6ba8 Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Thu, 20 Mar 2025 14:08:16 +0100
Subject: [PATCH 08/10] grid_gme: cleanup
---
src/Timstat2.cc | 2 +-
src/Vertintml.cc | 6 +-
src/grid_gme.cc | 305 ++++++++++++++++--------------------
src/knndata.cc | 2 +-
src/percentiles.cc | 4 +-
src/printinfo.cc | 14 +-
src/remap_method_conserv.cc | 2 +-
src/statistic.cc | 2 +-
8 files changed, 148 insertions(+), 189 deletions(-)
diff --git a/src/Timstat2.cc b/src/Timstat2.cc
index fb6a10024..c10304d97 100644
--- a/src/Timstat2.cc
+++ b/src/Timstat2.cc
@@ -111,7 +111,7 @@ correlation(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varr
if (fp_is_equal(cor, missval)) numMissVals++;
- pvalue = (nvals <= 2) ? missval : ((fabs(cor) < 1) ? calc_pvalue(cor, nvals) : 1);
+ pvalue = (nvals <= 2) ? missval : ((std::fabs(cor) < 1) ? calc_pvalue(cor, nvals) : 1);
}
else
{
diff --git a/src/Vertintml.cc b/src/Vertintml.cc
index 377325442..8fd50db07 100644
--- a/src/Vertintml.cc
+++ b/src/Vertintml.cc
@@ -83,7 +83,7 @@ static void
check_range(double rmin, double rmax, ForwardIt first, ForwardIt last, std::string const &name)
{
const auto [min, max] = std::minmax_element(first, last);
- if (fabs(*max - *min) <= 1.0e-9 || (*min < rmin) || (*max > rmax))
+ if (std::fabs(*max - *min) <= 1.0e-9 || (*min < rmin) || (*max > rmax))
cdo_warning("%s out of range (%g - %g) min=%g max=%g", name, rmin, rmax, *min, *max);
}
@@ -459,7 +459,7 @@ pressure_level_interpolation(Varray<double> &pressureLevels, bool useHeightLevel
{
cdo_def_field(streamID2, varID, levelID);
cdo_write_field(streamID2, interpVars[varID] ? vardata2[varID] : vardata1[varID], levelID,
- varnumMissVals[varID][levelID]);
+ varnumMissVals[varID][levelID]);
}
}
}
@@ -667,7 +667,7 @@ height_level_interpolation(Varray<double> &heightLevels, bool extrapolate)
{
cdo_def_field(streamID2, varID, levelID);
cdo_write_field(streamID2, interpVars[varID] ? vardata2[varID] : vardata1[varID], levelID,
- varnumMissVals[varID][levelID]);
+ varnumMissVals[varID][levelID]);
}
}
}
diff --git a/src/grid_gme.cc b/src/grid_gme.cc
index 8d50127c6..a7ddc92a3 100644
--- a/src/grid_gme.cc
+++ b/src/grid_gme.cc
@@ -26,7 +26,7 @@ struct array
struct cart
{
- double x[3] = { 0 };
+ double x[3] = { 0, 0, 0 };
};
struct geo
@@ -42,15 +42,15 @@ struct polygon
struct geo boundary[6];
};
-static const double pid5 = 0.2 * M_PI;
-// const double pid180 = 180.0/M_PI;
+constexpr double pid5 = 0.2 * M_PI;
+// constexpr double pid180 = 180.0/M_PI;
-static const int ispokes[12] = {
+constexpr int ispokes[12] = {
1, 0, -1, -1, 0, 1, 0, 1, 1, 0, -1, -1,
};
-static const int pentagon = 5;
-static const int hexagon = 6;
+constexpr int pentagon = 5;
+constexpr int hexagon = 6;
/*****************************************************************************/
@@ -92,16 +92,11 @@ gme_factorni(int kni, int *kni2, int *kni3)
*kni3 = *kni3 + 1;
mx = mx / 3;
}
- else
- { /* error return */
- }
+ else { /* error return */ }
}
- /* kni3 must not be greater than */
-
- if (*kni3 > 1)
- { /* error return */
- }
+ // kni3 must not be greater than
+ if (*kni3 > 1) { /* error return */ }
}
/*****************************************************************************/
@@ -109,14 +104,14 @@ gme_factorni(int kni, int *kni2, int *kni3)
static int
pow_ii(int x, int n)
{
- int pow;
-
if (n <= 0)
{
if (n == 0 || x == 1) return 1;
if (x != -1) return (x == 0) ? 1 / x : 0;
n = -n;
}
+
+ int pow;
for (pow = 1;;)
{
if (n & 01) pow *= x;
@@ -139,9 +134,9 @@ circum_center(struct cart *v0, struct cart *v1, struct cart *v2)
struct cart e2;
struct cart cu;
- double *ptmp1 = ((double *) e1.x);
- double *ptmp2 = ((double *) v1->x);
- double *ptmp3 = ((double *) v0->x);
+ double *ptmp1 = e1.x;
+ double *ptmp2 = v1->x;
+ double *ptmp3 = v0->x;
for (int j = 0; j < 3; ++j)
{
@@ -153,9 +148,9 @@ circum_center(struct cart *v0, struct cart *v1, struct cart *v2)
ptmp2++;
}
- ptmp1 = ((double *) e2.x);
- ptmp2 = ((double *) v2->x);
- ptmp3 = ((double *) v0->x);
+ ptmp1 = e2.x;
+ ptmp2 = v2->x;
+ ptmp3 = v0->x;
for (int j = 0; j < 3; ++j)
{
{
@@ -170,10 +165,10 @@ circum_center(struct cart *v0, struct cart *v1, struct cart *v2)
cu.x[1] = e1.x[2] * e2.x[0] - e1.x[0] * e2.x[2];
cu.x[2] = e1.x[0] * e2.x[1] - e1.x[1] * e2.x[0];
- const auto cnorm = std::sqrt(cu.x[0] * cu.x[0] + cu.x[1] * cu.x[1] + cu.x[2] * cu.x[2]);
+ auto cnorm = std::sqrt(cu.x[0] * cu.x[0] + cu.x[1] * cu.x[1] + cu.x[2] * cu.x[2]);
- ptmp1 = ((double *) center.x);
- ptmp2 = ((double *) cu.x);
+ ptmp1 = center.x;
+ ptmp2 = cu.x;
for (int j = 0; j < 3; ++j)
{
{
@@ -191,10 +186,10 @@ circum_center(struct cart *v0, struct cart *v1, struct cart *v2)
static struct cart
gc2cc(struct geo *position)
{
- const auto sln = std::sin(position->lon);
- const auto cln = std::cos(position->lon);
- const auto slt = std::sin(position->lat);
- const auto clt = std::cos(position->lat);
+ auto sln = std::sin(position->lon);
+ auto cln = std::cos(position->lon);
+ auto slt = std::sin(position->lat);
+ auto clt = std::cos(position->lat);
struct cart x;
x.x[0] = cln * clt;
@@ -211,21 +206,21 @@ cc2gc(struct cart *x)
{
struct geo position;
- if (fabs(x->x[0]) <= 0.0) { position.lon = (x->x[1] >= 0.0) ? 0.5 * M_PI : 1.5 * M_PI; }
+ if (std::fabs(x->x[0]) <= 0.0) { position.lon = (x->x[1] >= 0.0) ? 0.5 * M_PI : 1.5 * M_PI; }
else
{
- const auto tln = x->x[1] / x->x[0];
+ auto tln = x->x[1] / x->x[0];
position.lon = std::atan(tln);
if (x->x[0] < 0.0) position.lon += M_PI;
if (position.lon < 0.0) position.lon += 2 * M_PI;
}
- const auto r = std::sqrt(x->x[0] * x->x[0] + x->x[1] * x->x[1]);
+ auto r = std::sqrt(x->x[0] * x->x[0] + x->x[1] * x->x[1]);
- if (fabs(r) <= 0.0) { position.lat = (x->x[2] > 0.0) ? 0.5 * M_PI : -0.5 * M_PI; }
+ if (std::fabs(r) <= 0.0) { position.lat = (x->x[2] > 0.0) ? 0.5 * M_PI : -0.5 * M_PI; }
else
{
- const auto tlt = x->x[2] / r;
+ auto tlt = x->x[2] / r;
position.lat = std::atan(tlt);
}
@@ -244,7 +239,6 @@ boundary(struct polygon *poly, int kip1s, int kip1e, int kip2s, int kip2e, int k
struct cart x2;
struct cart v[6];
- int j1, j2, jd;
int jm, jm1, jm2;
struct array polyinfo;
@@ -276,11 +270,11 @@ boundary(struct polygon *poly, int kip1s, int kip1e, int kip2s, int kip2e, int k
tmp5 = polyinfo.dim[2].mult;
tmp3 = polyinfo.offset;
- for (jd = 1; jd <= knd; jd++)
+ for (int jd = 1; jd <= knd; jd++)
{
- for (j2 = kip2s; j2 <= kip2e; ++j2)
+ for (int j2 = kip2s; j2 <= kip2e; ++j2)
{
- for (j1 = kip1s; j1 <= kip1e; ++j1)
+ for (int j1 = kip1s; j1 <= kip1e; ++j1)
{
ptmp1 = &poly[j1 + tmp4 * j2 + tmp5 * jd + tmp3];
c = gc2cc(&ptmp1->center);
@@ -311,9 +305,6 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
int kip1e, int kip2s, int kip2e, int knd)
{
struct polygon *ptmp1;
-
- int j1, j2, jd, jm, js1, js2;
-
struct array px1info, px2info, polyinfo;
int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9;
@@ -394,26 +385,22 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
tmp15 = polyinfo.dim[2].mult;
tmp13 = polyinfo.offset;
- for (jd = 1; jd <= kndx; jd++)
+ for (int jd = 1; jd <= kndx; jd++)
{
- for (j2 = kipx2s + 1; j2 <= kipx2e - 1; ++j2)
+ for (int j2 = kipx2s + 1; j2 <= kipx2e - 1; ++j2)
{
- for (j1 = kipx1s + 1; j1 <= kipx1e - 1; ++j1)
+ for (int j1 = kipx1s + 1; j1 <= kipx1e - 1; ++j1)
{
-
ptmp1 = &poly[j1 + tmp14 * j2 + tmp15 * jd + tmp13];
-
ptmp1->center.lon = px1[j1 + tmp4 * j2 + tmp5 * jd + tmp3];
ptmp1->center.lat = px2[j1 + tmp9 * j2 + tmp10 * jd + tmp8];
if (j1 == kipx1s + 1 && j2 == kipx2s + 1)
{
-
ptmp1->type = pentagon;
- for (jm = 1; jm <= 5; jm++)
+ for (int jm = 1; jm <= 5; jm++)
{
-
if (jd < 6)
{
ptmp1->boundary[jm - 1].lon = px1[kipx1s + 1 + tmp4 * (kipx2s + 2) + tmp5 * (jm) + tmp3];
@@ -428,7 +415,6 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
}
else if (j1 == kipx1e - 1 && j2 == kipx2s + 1)
{
-
ptmp1->type = pentagon;
ptmp1->boundary[0].lon = px1[kipx1e - 1 + tmp4 * (kipx2s + 2) + tmp5 * jd + tmp3];
@@ -444,7 +430,6 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
}
else if (j1 == kipx1s + 1 && j2 == kipx2e - 1)
{
-
ptmp1->type = pentagon;
ptmp1->boundary[0].lon = px1[kipx1s + 2 + tmp4 * (kipx2e - 2) + tmp5 * jd + tmp3];
@@ -460,7 +445,6 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
}
else if (j1 == kipx1e - 1 && j2 == kipx2e - 1)
{
-
ptmp1->type = pentagon;
ptmp1->boundary[0].lon = px1[kipx1e + tmp4 * (kipx2e) + tmp5 * jd + tmp3];
@@ -476,15 +460,12 @@ neighbours(double *px1, double *px2, int kipx1s, int kipx1e, int kipx2s, int kip
}
else
{
-
- for (jm = 1; jm <= 6; jm++)
+ for (int jm = 1; jm <= 6; jm++)
{
-
ptmp1->type = hexagon;
- js1 = j1 + ispokes[jm - 1];
- js2 = j2 + ispokes[jm + 5];
-
+ int js1 = j1 + ispokes[jm - 1];
+ int js2 = j2 + ispokes[jm + 5];
ptmp1->boundary[jm - 1].lon = px1[js1 + tmp4 * js2 + tmp5 * jd + tmp3];
ptmp1->boundary[jm - 1].lat = px2[js1 + tmp9 * js2 + tmp10 * jd + tmp8];
}
@@ -502,10 +483,6 @@ static void
xd(double *p, int kip1s, int kip1e, int kip2s, int kip2e, int knd, double *px, int kipx1s, int kipx1e, int kipx2s, int kipx2e,
int kndx)
{
- int mi1sm1, mi1ep1, mi2sm1, mi2ep1;
- int mns, mpe, mpw, maw, mae, mpp;
- int j, j1, j2, jd;
-
struct array pinfo, pxinfo;
int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10;
@@ -563,28 +540,28 @@ xd(double *p, int kip1s, int kip1e, int kip2s, int kip2e, int knd, double *px, i
// tmp2 = pxinfo.dim[1].extent;
// tmp6 = pxinfo.dim[2].extent;
- for (jd = 1; jd <= knd; jd++)
+ for (int jd = 1; jd <= knd; jd++)
{
- for (j2 = kip2s; j2 <= kip2e; ++j2)
+ for (int j2 = kip2s; j2 <= kip2e; ++j2)
{
- for (j1 = kip1s; j1 <= kip1e; ++j1) { px[j1 + tmp9 * j2 + tmp10 * jd + tmp8] = p[j1 + tmp4 * j2 + tmp5 * jd + tmp3]; }
+ for (int j1 = kip1s; j1 <= kip1e; ++j1) { px[j1 + tmp9 * j2 + tmp10 * jd + tmp8] = p[j1 + tmp4 * j2 + tmp5 * jd + tmp3]; }
}
}
- mi1sm1 = kip1s - 1;
- mi1ep1 = kip1e + 1;
- mi2sm1 = kip2s - 1;
- mi2ep1 = kip2e + 1;
+ int mi1sm1 = kip1s - 1;
+ int mi1ep1 = kip1e + 1;
+ int mi2sm1 = kip2s - 1;
+ int mi2ep1 = kip2e + 1;
- for (jd = 1; jd <= knd; jd++)
+ for (int jd = 1; jd <= knd; jd++)
{
- mns = (jd - 1) / 5;
- mpe = jd + 1 - (jd / (5 * (1 + mns))) * 5;
- mpw = jd - 1 + ((mns * 10 + 6 - jd) / (5 * (1 + mns))) * 5;
- mae = jd + 5 - 9 * mns - 5 * (jd / 10);
- maw = jd + 4 + ((6 - jd) / 5) * 5 - 9 * mns;
- mpp = jd + 3 - ((jd + 2) / 5) * 5 + 5 * mns;
- for (j = kip2s; j <= kip1e; ++j)
+ int mns = (jd - 1) / 5;
+ int mpe = jd + 1 - (jd / (5 * (1 + mns))) * 5;
+ int mpw = jd - 1 + ((mns * 10 + 6 - jd) / (5 * (1 + mns))) * 5;
+ int mae = jd + 5 - 9 * mns - 5 * (jd / 10);
+ int maw = jd + 4 + ((6 - jd) / 5) * 5 - 9 * mns;
+ int mpp = jd + 3 - ((jd + 2) / 5) * 5 + 5 * mns;
+ for (int j = kip2s; j <= kip1e; ++j)
{
px[j + tmp9 * mi2sm1 + tmp10 * jd + tmp8] = p[kip1s + 1 + tmp4 * j + tmp5 * mpw + tmp3];
px[mi1sm1 + tmp9 * (j + 1) + tmp10 * jd + tmp8] = p[j - 1 + tmp4 * (kip2s + 1) + tmp5 * mpe + tmp3];
@@ -613,10 +590,6 @@ static void
tricntr(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kjd, int kni)
{
(void) knd;
- double r1, r2, r3;
-
- int mi1, mi2;
- double zxnorm;
int id1 = kig1e - kig1s + 1;
int id2 = id1 * (kig2e - kig2s + 1);
@@ -625,8 +598,8 @@ tricntr(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kj
for (int j = 1; j <= 2; ++j)
{
- mi1 = j * kni / 3;
- mi2 = (j - 1) * kni + 1;
+ int mi1 = j * kni / 3;
+ int mi2 = (j - 1) * kni + 1;
pxn[mi1 + id1 * (mi1 + 1) + id2 * 1 + id3 * kjd + ioffset] = pxn[mi2 - 1 + id1 * (mi2) + id2 * 1 + id3 * kjd + ioffset]
+ pxn[kni + id1 * 1 + id2 * 1 + id3 * kjd + ioffset]
+ pxn[0 + id1 * (kni + 1) + id2 * 1 + id3 * kjd + ioffset];
@@ -638,10 +611,10 @@ tricntr(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kj
+ pxn[0 + id1 * (kni + 1) + id2 * 3 + id3 * kjd + ioffset];
/* Normalize to unit-sphere */
- r1 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 1 + id3 * kjd + ioffset];
- r2 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 2 + id3 * kjd + ioffset];
- r3 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 3 + id3 * kjd + ioffset];
- zxnorm = 1.0 / std::sqrt(r1 * r1 + r2 * r2 + r3 * r3);
+ double r1 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 1 + id3 * kjd + ioffset];
+ double r2 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 2 + id3 * kjd + ioffset];
+ double r3 = pxn[mi1 + id1 * (mi1 + 1) + id2 * 3 + id3 * kjd + ioffset];
+ double zxnorm = 1.0 / std::sqrt(r1 * r1 + r2 * r2 + r3 * r3);
pxn[mi1 + id1 * (mi1 + 1) + id2 * 1 + id3 * kjd + ioffset]
= zxnorm * pxn[mi1 + id1 * (mi1 + 1) + id2 * 1 + id3 * kjd + ioffset];
@@ -675,15 +648,14 @@ gcpt(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kjd,
double zchord = std::sqrt((r1 * r1) + (r2 * r2) + (r3 * r3));
- /* Calculate "ztheta", the great circle angle between x1 and x2 */
+ // Calculate "ztheta", the great circle angle between x1 and x2
double ztheta = 2.0 * std::asin(zchord * 0.5);
- /* Calculate the weighting factors which follow from the condition */
- /* that x is a point on the unit-sphere, too. */
+ // Calculate the weighting factors which follow from the condition that x is a point on the unit-sphere, too.
double zbeta = std::sin(pgamma * ztheta) / std::sin(ztheta);
double zalpha = std::sin((1.0 - pgamma) * ztheta) / std::sin(ztheta);
- /* Store the (x,y,z) coordinates of the point x into the array pxn */
+ // Store the (x,y,z) coordinates of the point x into the array pxn
pxn[ki + id1 * kj + id2 * 1 + id3 * kjd + ioffset] = zalpha * pxn[ki1 + id1 * kj1 + id2 * 1 + id3 * kjd + ioffset]
+ zbeta * pxn[ki2 + id1 * kj2 + id2 * 1 + id3 * kjd + ioffset];
pxn[ki + id1 * kj + id2 * 2 + id3 * kjd + ioffset] = zalpha * pxn[ki1 + id1 * kj1 + id2 * 2 + id3 * kjd + ioffset]
@@ -700,17 +672,14 @@ static void
glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kni2, int kni3)
{
double zsgn;
- int j1, j2;
double zrlon;
- int jb, jd, ml, mm;
+ int ml, mm;
double zgamma;
int mi1, mi2, ml2, ml3;
/*
- * Calculate the Cartesian coordinates of the gridpoints of the
- * icosahedral grid on the unit sphere. The grid resolution
- * corresponds to a subdivision of the edges of the original
- * icosahedral triangles into mni equal parts.
+ * Calculate the Cartesian coordinates of the gridpoints of the icosahedral grid on the unit sphere.
+ * The grid resolution corresponds to a subdivision of the edges of the original icosahedral triangles into mni equal parts.
*/
std::vector<int> mcosv(knd);
@@ -729,7 +698,7 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/* Compute the local array mcosv, i.e. the meridian angle locations */
- for (jd = 1; jd <= knd; ++jd)
+ for (int jd = 1; jd <= knd; ++jd)
{
if (jd % 2 == 1) { mcosv[(jd + 1) / 2 - 1] = jd - 2 - knd * ((jd - 1) / 7); }
else { mcosv[jd / 2 + 4] = jd - 2 - knd * ((jd - 1) / 7); }
@@ -740,9 +709,8 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/* coordinates and then iteratively bisecting them kni2 times. */
/* First a trisection is performed, if required (kni3=1). */
- for (jd = 1; jd <= knd; ++jd)
+ for (int jd = 1; jd <= knd; ++jd)
{
-
/* Toggle the hemisphere */
if (jd >= 6) { zsgn = -1.0; }
else { zsgn = 1.0; }
@@ -786,9 +754,9 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
ml3 = mni / 3;
/* Trisect the rows of the diamond. */
- for (j1 = 1; j1 <= 2; ++j1)
+ for (int j1 = 1; j1 <= 2; ++j1)
{
- for (j2 = 1; j2 <= 2; ++j2)
+ for (int j2 = 1; j2 <= 2; ++j2)
{
mi1 = (j1 - 1) * mni;
mi2 = j2 * ml3 + 1;
@@ -798,9 +766,9 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
}
/* Trisect the columns of the diamond. */
- for (j1 = 1; j1 <= 2; ++j1)
+ for (int j1 = 1; j1 <= 2; ++j1)
{
- for (j2 = 1; j2 <= 2; ++j2)
+ for (int j2 = 1; j2 <= 2; ++j2)
{
mi1 = j2 * ml3;
mi2 = (j1 - 1) * mni + 1;
@@ -810,7 +778,7 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
}
/* Trisect the diagonal of the diamond. */
- for (j2 = 1; j2 <= 2; ++j2)
+ for (int j2 = 1; j2 <= 2; ++j2)
{
mi1 = mni - j2 * ml3;
mi2 = j2 * ml3 + 1;
@@ -826,17 +794,16 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/* Find the coordinates of the triangle nodes by iteratively */
/* bisecting the diamond intervals. */
- for (jb = 0; jb <= kni2 - 1; ++jb)
+ for (int jb = 0; jb <= kni2 - 1; ++jb)
{
mm = pow_ii(3, kni3) * pow_ii(2, jb);
ml = mni / mm;
ml2 = ml / 2;
- /* Compute the rows of the diamond. */
-
- for (j1 = 1; j1 <= mm + 1; ++j1)
+ // Compute the rows of the diamond.
+ for (int j1 = 1; j1 <= mm + 1; ++j1)
{
- for (j2 = 1; j2 <= mm; ++j2)
+ for (int j2 = 1; j2 <= mm; ++j2)
{
mi1 = (j1 - 1) * ml;
mi2 = (j2 - 1) * ml + ml2 + 1;
@@ -844,11 +811,10 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
}
}
- /* Compute the columns of diamond. */
-
- for (j1 = 1; j1 <= mm + 1; ++j1)
+ // Compute the columns of diamond.
+ for (int j1 = 1; j1 <= mm + 1; ++j1)
{
- for (j2 = 1; j2 <= mm; ++j2)
+ for (int j2 = 1; j2 <= mm; ++j2)
{
mi1 = (j2 - 1) * ml + ml2;
mi2 = (j1 - 1) * ml + 1;
@@ -856,11 +822,10 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
}
}
- /* Compute the diagonals of the diamond. */
-
- for (j1 = 1; j1 <= mm; ++j1)
+ // Compute the diagonals of the diamond.
+ for (int j1 = 1; j1 <= mm; ++j1)
{
- for (j2 = 1; j2 <= mm; ++j2)
+ for (int j2 = 1; j2 <= mm; ++j2)
{
mi1 = (j1 - 1) * ml + ml2;
mi2 = (j2 - 1) * ml + ml2 + 1;
@@ -872,19 +837,19 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/***********************************************************************/
/* Set pxn to 0 if it is less than 2.5 e-14 to avoid round-off errors */
- for (j2 = kig2s; j2 <= kig2e; ++j2)
+ for (int j2 = kig2s; j2 <= kig2e; ++j2)
{
- for (j1 = kig1s; j1 <= kig1e; ++j1)
+ for (int j1 = kig1s; j1 <= kig1e; ++j1)
{
- if (fabs(pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset]) < 2.5e-14)
+ if (std::fabs(pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset]) < 2.5e-14)
{
pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset] = 0.0;
}
- if (fabs(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset]) < 2.5e-14)
+ if (std::fabs(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset]) < 2.5e-14)
{
pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset] = 0.0;
}
- if (fabs(pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset]) < 2.5e-14)
+ if (std::fabs(pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset]) < 2.5e-14)
{
pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset] = 0.0;
}
@@ -896,30 +861,27 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/* Calculate the longitude "prlon" and the latitude "prlat"; */
/* only for the core of the diamonds, not the extended ones. */
- for (jd = 1; jd <= knd; ++jd)
+ for (int jd = 1; jd <= knd; ++jd)
{
- for (j2 = kig2s; j2 <= kig2e; ++j2)
+ for (int j2 = kig2s; j2 <= kig2e; ++j2)
{
- for (j1 = kig1s; j1 <= kig1e; ++j1)
+ for (int j1 = kig1s; j1 <= kig1e; ++j1)
{
- prlon[j1 + id1 * j2 + id2 * jd + joffset] = atan2(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset],
- pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset] + 1.0e-20);
+ prlon[j1 + id1 * j2 + id2 * jd + joffset] = std::atan2(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset],
+ pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset] + 1.0e-20);
prlat[j1 + id1 * j2 + id2 * jd + joffset] = std::asin(pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset]);
}
}
}
return;
-} /* glo_coor */
+} // glo_coor
/*****************************************************************************/
static void
initmask(int *mask, int ni, int nd)
{
- int tmp1, tmp2, /*tmp3, tmp4, tmp5,*/ tmp6, tmp7, tmp8, tmp9;
-
- int *ptmp1;
char *ptmp2;
struct array section;
@@ -927,12 +889,12 @@ initmask(int *mask, int ni, int nd)
maskinfo.offset = 0;
maskinfo.dim[0].lower = 0;
- tmp1 = ni + 1;
+ int tmp1 = ni + 1;
if (tmp1 < 0) tmp1 = 0;
maskinfo.dim[0].extent = tmp1;
maskinfo.dim[0].mult = 1;
maskinfo.offset -= 0;
- tmp2 = tmp1;
+ int tmp2 = tmp1;
maskinfo.dim[1].lower = 1;
tmp1 = ni + 1;
if (tmp1 < 0) tmp1 = 0;
@@ -951,14 +913,14 @@ initmask(int *mask, int ni, int nd)
// tmp3 = maskinfo.offset;
tmp1 = maskinfo.dim[0].extent;
tmp2 = maskinfo.dim[1].extent;
- tmp9 = maskinfo.dim[2].extent;
+ int tmp9 = maskinfo.dim[2].extent;
- ptmp1 = mask;
- for (tmp8 = 0; tmp8 < tmp9; tmp8++)
+ int *ptmp1 = mask;
+ for (int tmp8 = 0; tmp8 < tmp9; tmp8++)
{
- for (tmp7 = 0; tmp7 < tmp2; tmp7++)
+ for (int tmp7 = 0; tmp7 < tmp2; tmp7++)
{
- for (tmp6 = 0; tmp6 < tmp1; tmp6++) *ptmp1++ = 1;
+ for (int tmp6 = 0; tmp6 < tmp1; tmp6++) *ptmp1++ = 1;
}
}
@@ -977,7 +939,7 @@ initmask(int *mask, int ni, int nd)
tmp1 *= maskinfo.dim[0].extent;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
break;
case 5:
tmp1 = 1;
@@ -987,7 +949,7 @@ initmask(int *mask, int ni, int nd)
tmp1 *= maskinfo.dim[0].extent;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 4;
section.addr = (char *) mask;
@@ -1002,7 +964,7 @@ initmask(int *mask, int ni, int nd)
section.addr += tmp1 * (jd - 1);
tmp2 = section.dim[0].extent;
ptmp2 = section.addr;
- for (tmp6 = 0; tmp6 < tmp2; tmp6++)
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++)
{
*((int *) ptmp2) = 0;
ptmp2 += section.dim[0].mult;
@@ -1018,7 +980,7 @@ initmask(int *mask, int ni, int nd)
ptmp1 += tmp1 * ni;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 4;
section.addr = (char *) mask;
@@ -1034,7 +996,7 @@ initmask(int *mask, int ni, int nd)
section.addr += tmp1 * (jd - 1);
tmp2 = section.dim[0].extent;
ptmp2 = section.addr;
- for (tmp6 = 0; tmp6 < tmp2; tmp6++)
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++)
{
*((int *) ptmp2) = 0;
ptmp2 += section.dim[0].mult;
@@ -1050,7 +1012,7 @@ initmask(int *mask, int ni, int nd)
tmp1 *= maskinfo.dim[0].extent;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 1;
ptmp1 = mask;
@@ -1060,7 +1022,7 @@ initmask(int *mask, int ni, int nd)
ptmp1 += tmp1 * ni;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 4;
section.addr = (char *) mask;
@@ -1076,7 +1038,7 @@ initmask(int *mask, int ni, int nd)
section.addr += tmp1 * (jd - 1);
tmp2 = section.dim[0].extent;
ptmp2 = section.addr;
- for (tmp6 = 0; tmp6 < tmp2; tmp6++)
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++)
{
*((int *) ptmp2) = 0;
ptmp2 += section.dim[0].mult;
@@ -1091,7 +1053,7 @@ initmask(int *mask, int ni, int nd)
tmp1 *= maskinfo.dim[0].extent;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 1;
ptmp1 = mask;
@@ -1101,7 +1063,7 @@ initmask(int *mask, int ni, int nd)
ptmp1 += tmp1 * ni;
tmp1 *= maskinfo.dim[1].extent;
ptmp1 += tmp1 * (jd - 1);
- for (tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++) *ptmp1++ = 0;
tmp1 = 4;
section.addr = (char *) mask;
@@ -1116,7 +1078,7 @@ initmask(int *mask, int ni, int nd)
section.addr += tmp1 * (jd - 1);
tmp2 = section.dim[0].extent;
ptmp2 = section.addr;
- for (tmp6 = 0; tmp6 < tmp2; tmp6++)
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++)
{
*((int *) ptmp2) = 0;
ptmp2 += section.dim[0].mult;
@@ -1136,7 +1098,7 @@ initmask(int *mask, int ni, int nd)
section.addr += tmp1 * (jd - 1);
tmp2 = section.dim[0].extent;
ptmp2 = section.addr;
- for (tmp6 = 0; tmp6 < tmp2; tmp6++)
+ for (int tmp6 = 0; tmp6 < tmp2; tmp6++)
{
*((int *) ptmp2) = 0;
ptmp2 += section.dim[0].mult;
@@ -1154,8 +1116,6 @@ template <typename T>
void
gme_grid_restore(T *p, int ni, int nd)
{
- int j;
-
struct array pinfo;
pinfo.offset = 0;
@@ -1191,30 +1151,30 @@ gme_grid_restore(T *p, int ni, int nd)
case 2:
case 3:
case 4:
- for (j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
break;
case 5:
- for (j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
- for (j = 0; j <= ni; ++j) p[tmp4 * (j + 1) + tmp5 * 5 + tmp3] = p[j + tmp4 + tmp5 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
+ for (int j = 0; j <= ni; ++j) p[tmp4 * (j + 1) + tmp5 * 5 + tmp3] = p[j + tmp4 + tmp5 + tmp3];
break;
case 6:
- for (j = 0; j <= ni; ++j) p[j + tmp4 * (ni + 1) + tmp5 * 6 + tmp3] = p[ni + tmp4 * (ni + 1 - j) + tmp5 * 2 + tmp3];
- for (j = 0; j <= ni; ++j) p[ni + tmp4 * (j + 1) + tmp5 * 6 + tmp3] = p[ni - j + tmp4 * (ni + 1) + tmp5 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 * (ni + 1) + tmp5 * 6 + tmp3] = p[ni + tmp4 * (ni + 1 - j) + tmp5 * 2 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[ni + tmp4 * (j + 1) + tmp5 * 6 + tmp3] = p[ni - j + tmp4 * (ni + 1) + tmp5 + tmp3];
break;
case 7:
case 8:
case 9:
- for (j = 0; j <= ni; ++j)
+ for (int j = 0; j <= ni; ++j)
p[j + tmp4 * (ni + 1) + tmp5 * jd + tmp3] = p[ni + tmp4 * (ni + 1 - j) + tmp5 * (jd - 4) + tmp3];
- for (j = 0; j <= ni; ++j)
+ for (int j = 0; j <= ni; ++j)
p[ni + tmp4 * (j + 1) + tmp5 * jd + tmp3] = p[ni - j + tmp4 * (ni + 1) + tmp5 * (jd - 5) + tmp3];
- for (j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * jd + tmp3] = p[tmp4 * (j + 1) + tmp5 * (jd - 1) + tmp3];
break;
case 10:
- for (j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * 10 + tmp3] = p[tmp4 * (j + 1) + tmp5 * 9 + tmp3];
- for (j = 0; j <= ni; ++j) p[tmp4 * (j + 1) + tmp5 * 10 + tmp3] = p[j + tmp4 + tmp5 * 6 + tmp3];
- for (j = 0; j <= ni; ++j) p[j + tmp4 * (ni + 1) + tmp5 * 10 + tmp3] = p[ni + tmp4 * (ni + 1 - j) + tmp5 + tmp3];
- for (j = 0; j <= ni; ++j) p[ni + tmp4 * (j + 1) + tmp5 * 10 + tmp3] = p[ni - j + tmp4 * (ni + 1) + tmp5 * 5 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 + tmp5 * 10 + tmp3] = p[tmp4 * (j + 1) + tmp5 * 9 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[tmp4 * (j + 1) + tmp5 * 10 + tmp3] = p[j + tmp4 + tmp5 * 6 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[j + tmp4 * (ni + 1) + tmp5 * 10 + tmp3] = p[ni + tmp4 * (ni + 1 - j) + tmp5 + tmp3];
+ for (int j = 0; j <= ni; ++j) p[ni + tmp4 * (j + 1) + tmp5 * 10 + tmp3] = p[ni - j + tmp4 * (ni + 1) + tmp5 * 5 + tmp3];
break;
}
}
@@ -1232,7 +1192,7 @@ void
gme_grid(int withBounds, size_t gridsize, double *rlon, double *rlat, double *blon, double *blat, int *imask, int ni, int nd,
int ni2, int ni3)
{
- /* check gridsize */
+ // check gridsize
if ((size_t) (ni + 1) * (ni + 1) * nd != gridsize)
{
fprintf(stderr, "gme_grid: Calculation of the global GME grid failed (ni=%d)!\n", ni);
@@ -1290,7 +1250,6 @@ main(int argc, char *argv[])
{
int ni2, ni3;
int im1s, im1e, im2s, im2e;
- int j1, j2, jd;
int nd = 10;
int ni = 2;
@@ -1361,11 +1320,11 @@ main(int argc, char *argv[])
exit(-1);
}
- for (jd = 1; jd <= nd; jd++)
+ for (int jd = 1; jd <= nd; jd++)
{
- for (j2 = 1; j2 <= ni + 1; ++j2)
+ for (int j2 = 1; j2 <= ni + 1; ++j2)
{
- for (j1 = 0; j1 <= ni; ++j1)
+ for (int j1 = 0; j1 <= ni; ++j1)
{
if (mask[j1 + id1 * j2 + id2 * jd + ioffset])
{
@@ -1393,11 +1352,11 @@ main(int argc, char *argv[])
double total_area = 0.0;
- for (jd = 1; jd <= nd; jd++)
+ for (int jd = 1; jd <= nd; jd++)
{
- for (j2 = 1; j2 <= ni + 1; ++j2)
+ for (int j2 = 1; j2 <= ni + 1; ++j2)
{
- for (j1 = 0; j1 <= ni; ++j1)
+ for (int j1 = 0; j1 <= ni; ++j1)
{
area[j1 + id1 * j2 + id2 * jd + ioffset] = 0.0;
if (mask[j1 + id1 * j2 + id2 * jd + ioffset])
diff --git a/src/knndata.cc b/src/knndata.cc
index 1759b509d..435eabb18 100644
--- a/src/knndata.cc
+++ b/src/knndata.cc
@@ -163,7 +163,7 @@ KnnData::compute_weights_gauss()
// If the sum of the weights is very low, which can happen in case
// the target point is very far away from the group source points.
- if (fabs(weights_sum) < 1e-9)
+ if (std::fabs(weights_sum) < 1e-9)
{
// Due to limited accuracy the exact contribution of each source
// point cannot be computed. Therefore, the normalisation would
diff --git a/src/percentiles.cc b/src/percentiles.cc
index 2fd9db8e2..d85433e0c 100644
--- a/src/percentiles.cc
+++ b/src/percentiles.cc
@@ -152,7 +152,7 @@ percentile_numpy(T *array, size_t n, double quantile)
constexpr double fuzz = 4.0 * std::numeric_limits<double>::epsilon();
size_t j = std::floor(nppn + fuzz); // m = a + probs*(1 - a - b)
double h = nppn - j;
- if (fabs(h) < fuzz) h = 0.0;
+ if (std::fabs(h) < fuzz) h = 0.0;
if (h > 0.0 && h < 1.0)
percentil = (1.0 - h) * get_nth_element(array, n, j - 1) + h * get_nth_element(array, n, j);
else
@@ -165,7 +165,7 @@ percentile_numpy(T *array, size_t n, double quantile)
size_t j = std::floor(nppm);
double h = (numpyMethod == NumpyMethod::inverted_cdf) ? (nppm > j)
: (numpyMethod == NumpyMethod::averaged_inverted_cdf) ? ((nppm > j) + 1.0) / 2.0
- : ((fabs(nppm - j) > 0.0) | ((j % 2) == 1));
+ : ((std::fabs(nppm - j) > 0.0) | ((j % 2) == 1));
if (h > 0.0 && h < 1.0)
percentil = (1.0 - h) * get_nth_element(array, n, j - 1) + h * get_nth_element(array, n, j);
else
diff --git a/src/printinfo.cc b/src/printinfo.cc
index 13e17dbc8..1965812e9 100644
--- a/src/printinfo.cc
+++ b/src/printinfo.cc
@@ -183,9 +183,9 @@ calc_curvi_xinc(int gridID)
{
Varray<double> xvals(xsize);
for (size_t i = 0; i < xsize; ++i) xvals[i] = xvals2D[i];
- xinc = fabs(xvals[xsize - 1] - xvals[0]) / (xsize - 1);
+ xinc = std::fabs(xvals[xsize - 1] - xvals[0]) / (xsize - 1);
for (size_t i = 1; i < xsize; ++i)
- if (fabs(fabs(xvals[i - 1] - xvals[i]) - xinc) > 0.005 * xinc)
+ if (std::fabs(std::fabs(xvals[i - 1] - xvals[i]) - xinc) > 0.005 * xinc)
{
xinc = 0.0;
break;
@@ -216,9 +216,9 @@ calc_curvi_yinc(int gridID)
{
Varray<double> yvals(ysize);
for (size_t i = 0; i < ysize; ++i) yvals[i] = yvals2D[i * xsize];
- yinc = fabs(yvals[ysize - 1] - yvals[0]) / (ysize - 1);
+ yinc = std::fabs(yvals[ysize - 1] - yvals[0]) / (ysize - 1);
for (size_t i = 1; i < ysize; ++i)
- if (fabs(fabs(yvals[i - 1] - yvals[i]) - yinc) > 0.005 * yinc)
+ if (std::fabs(std::fabs(yvals[i - 1] - yvals[i]) - yinc) > 0.005 * yinc)
{
yinc = 0.0;
break;
@@ -480,7 +480,7 @@ printZaxisBoundsInfo(int zaxisID, int dig, int levelsize, double zinc, std::stri
{
auto level1 = zaxisInqLbound(zaxisID, 0);
auto level2 = zaxisInqUbound(zaxisID, 0);
- if (!(levelsize == 1 && is_equal(level1, level2) && fabs(level1) <= 0))
+ if (!(levelsize == 1 && is_equal(level1, level2) && std::fabs(level1) <= 0))
{
fprintf(stdout, "%33s : ", "bounds");
fprintf(stdout, "%.*g-%.*g", dig, level1, dig, level2);
@@ -526,7 +526,7 @@ printZaxisLevelInfo(int levelsize, int zaxisID, int zaxistype, double &zinc, int
Varray<double> levels(levelsize);
zaxisInqLevels(zaxisID, levels.data());
- if (!(zaxisTypeIsSingleLayer(zaxistype) && levelsize == 1 && fabs(levels[0]) <= 0))
+ if (!(zaxisTypeIsSingleLayer(zaxistype) && levelsize == 1 && std::fabs(levels[0]) <= 0))
{
auto zfirst = levels[0];
auto zlast = levels[levelsize - 1];
@@ -534,7 +534,7 @@ printZaxisLevelInfo(int levelsize, int zaxisID, int zaxistype, double &zinc, int
{
zinc = (levels[levelsize - 1] - levels[0]) / (levelsize - 1);
for (int levelID = 2; levelID < levelsize; ++levelID)
- if (fabs(fabs(levels[levelID] - levels[levelID - 1]) - zinc) > 0.001 * zinc)
+ if (std::fabs(std::fabs(levels[levelID] - levels[levelID - 1]) - zinc) > 0.001 * zinc)
{
zinc = 0;
break;
diff --git a/src/remap_method_conserv.cc b/src/remap_method_conserv.cc
index a7c3fc508..711d4ae64 100644
--- a/src/remap_method_conserv.cc
+++ b/src/remap_method_conserv.cc
@@ -59,7 +59,7 @@ get_edge_direction(const double *ref_corner, double *corner_a, double *corner_b)
// if the reference corner is directly on the edge
// (for small angles sin(x)==x)
- if (fabs(angle) < yac_angle_tol) return 0.0;
+ if (std::fabs(angle) < yac_angle_tol) return 0.0;
return copysign(1.0, angle);
}
diff --git a/src/statistic.cc b/src/statistic.cc
index 3e13208e9..c713459c6 100644
--- a/src/statistic.cc
+++ b/src/statistic.cc
@@ -449,7 +449,7 @@ beta_distr_inv(double a, double b, double p)
x = xx;
}
#if 0
- for (x_l = 0, x_r = 1; fabs(x_l - x_r) > eps;
+ for (x_l = 0, x_r = 1; std::fabs(x_l - x_r) > eps;
x = (x_l+x_r) / 2.0, (beta_distr(a, b, x) < p) ? (x_l=x):(x_r=x));
#endif
--
GitLab
From 2256e6ce9a2955018cb1c8f04df0a86d7ba761ac Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Thu, 20 Mar 2025 20:26:10 +0100
Subject: [PATCH 09/10] gme_grid: check if calculation of coordinates failed
---
ChangeLog | 4 ++
src/grid_gme.cc | 159 +++++++++++++++++++++++-------------------------
2 files changed, 80 insertions(+), 83 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 211f29c28..4ba6c738a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2025-03-20 Uwe Schulzweida
+
+ * gme_grid: check if calculation of coordinates failed
+
2025-03-15 Uwe Schulzweida
* timpctl: Short integer overflow for CDO_PCTL_NBINS>32768 (bug fix)
diff --git a/src/grid_gme.cc b/src/grid_gme.cc
index a7ddc92a3..864dc87d6 100644
--- a/src/grid_gme.cc
+++ b/src/grid_gme.cc
@@ -635,16 +635,19 @@ gcpt(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kjd,
{
(void) knd;
- /* Calculate "zchord", the Cartesian distance between x1 and x2 */
+ // Calculate "zchord", the Cartesian distance between x1 and x2
int id1 = kig1e - kig1s + 1;
int id2 = id1 * (kig2e - kig2s + 1);
int id3 = id2 * 3;
int ioffset = -(id1 + id2 + id3);
- double r1 = (pxn[ki2 + id1 * kj2 + id2 * 1 + id3 * kjd + ioffset] - pxn[ki1 + id1 * kj1 + id2 * 1 + id3 * kjd + ioffset]);
- double r2 = (pxn[ki2 + id1 * kj2 + id2 * 2 + id3 * kjd + ioffset] - pxn[ki1 + id1 * kj1 + id2 * 2 + id3 * kjd + ioffset]);
- double r3 = (pxn[ki2 + id1 * kj2 + id2 * 3 + id3 * kjd + ioffset] - pxn[ki1 + id1 * kj1 + id2 * 3 + id3 * kjd + ioffset]);
+ int kij = ki + id1 * kj + id3 * kjd + ioffset;
+ int kij1 = ki1 + id1 * kj1 + id3 * kjd + ioffset;
+ int kij2 = ki2 + id1 * kj2 + id3 * kjd + ioffset;
+ double r1 = (pxn[kij2 + id2 * 1] - pxn[kij1 + id2 * 1]);
+ double r2 = (pxn[kij2 + id2 * 2] - pxn[kij1 + id2 * 2]);
+ double r3 = (pxn[kij2 + id2 * 3] - pxn[kij1 + id2 * 3]);
double zchord = std::sqrt((r1 * r1) + (r2 * r2) + (r3 * r3));
@@ -656,27 +659,18 @@ gcpt(double *pxn, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kjd,
double zalpha = std::sin((1.0 - pgamma) * ztheta) / std::sin(ztheta);
// Store the (x,y,z) coordinates of the point x into the array pxn
- pxn[ki + id1 * kj + id2 * 1 + id3 * kjd + ioffset] = zalpha * pxn[ki1 + id1 * kj1 + id2 * 1 + id3 * kjd + ioffset]
- + zbeta * pxn[ki2 + id1 * kj2 + id2 * 1 + id3 * kjd + ioffset];
- pxn[ki + id1 * kj + id2 * 2 + id3 * kjd + ioffset] = zalpha * pxn[ki1 + id1 * kj1 + id2 * 2 + id3 * kjd + ioffset]
- + zbeta * pxn[ki2 + id1 * kj2 + id2 * 2 + id3 * kjd + ioffset];
- pxn[ki + id1 * kj + id2 * 3 + id3 * kjd + ioffset] = zalpha * pxn[ki1 + id1 * kj1 + id2 * 3 + id3 * kjd + ioffset]
- + zbeta * pxn[ki2 + id1 * kj2 + id2 * 3 + id3 * kjd + ioffset];
+ pxn[kij + id2 * 1] = zalpha * pxn[kij1 + id2 * 1] + zbeta * pxn[kij2 + id2 * 1];
+ pxn[kij + id2 * 2] = zalpha * pxn[kij1 + id2 * 2] + zbeta * pxn[kij2 + id2 * 2];
+ pxn[kij + id2 * 3] = zalpha * pxn[kij1 + id2 * 3] + zbeta * pxn[kij2 + id2 * 3];
return;
-} /* gcpt */
+} // gcpt
/****************************************************************************/
static void
glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int kig2s, int kig2e, int knd, int kni2, int kni3)
{
- double zsgn;
- double zrlon;
- int ml, mm;
- double zgamma;
- int mi1, mi2, ml2, ml3;
-
/*
* Calculate the Cartesian coordinates of the gridpoints of the icosahedral grid on the unit sphere.
* The grid resolution corresponds to a subdivision of the edges of the original icosahedral triangles into mni equal parts.
@@ -689,15 +683,13 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
int ioffset = -(id1 + id2 + id3);
int joffset = -(id1 + id2);
- /* Compute angles associated with the icosahedron. */
-
+ // Compute angles associated with the icosahedron.
double zw = std::acos(1.0 / (std::sin(pid5) * 2.0)) * 2.0;
double zcosw = std::cos(zw);
double zsinw = std::sin(zw);
int mni = pow_ii(2, kni2) * pow_ii(3, kni3);
- /* Compute the local array mcosv, i.e. the meridian angle locations */
-
+ // Compute the local array mcosv, i.e. the meridian angle locations
for (int jd = 1; jd <= knd; ++jd)
{
if (jd % 2 == 1) { mcosv[(jd + 1) / 2 - 1] = jd - 2 - knd * ((jd - 1) / 7); }
@@ -708,105 +700,102 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/* Loop over the ten diamonds computing diamond vertices (x,y,z) */
/* coordinates and then iteratively bisecting them kni2 times. */
/* First a trisection is performed, if required (kni3=1). */
-
for (int jd = 1; jd <= knd; ++jd)
{
- /* Toggle the hemisphere */
- if (jd >= 6) { zsgn = -1.0; }
- else { zsgn = 1.0; }
+ // Toggle the hemisphere
+ double zsgn = (jd >= 6) ? -1.0 : 1.0;
- /* Compute the meridian angle for each diamond "home" vertex. */
- zrlon = mcosv[jd - 1] * pid5;
+ // Compute the meridian angle for each diamond "home" vertex.
+ double zrlon = mcosv[jd - 1] * pid5;
/* Every diamond has one vertex at a pole (N or S). */
/* Label this point (0,1,,) in each diamond, and */
/* initialize it to have the (x,y,z) coordinates of */
/* the pole point on the unit sphere. */
- pxn[0 + id1 * 1 + id2 * 1 + id3 * jd + ioffset] = 0.0;
- pxn[0 + id1 * 1 + id2 * 2 + id3 * jd + ioffset] = 0.0;
- pxn[0 + id1 * 1 + id2 * 3 + id3 * jd + ioffset] = zsgn;
+ int offset = id3 * jd + ioffset;
+ pxn[0 + id1 + id2 * 1 + offset] = 0.0;
+ pxn[0 + id1 + id2 * 2 + offset] = 0.0;
+ pxn[0 + id1 + id2 * 3 + offset] = zsgn;
/* Now initialize the (x,y,z) coordinates of the "home" vertex, */
/* which defines which diamond we are talking about, and label */
/* this point (mni,1,,). */
- pxn[mni + id1 * 1 + id2 * 1 + id3 * jd + ioffset] = zsinw * std::cos(zrlon);
- pxn[mni + id1 * 1 + id2 * 2 + id3 * jd + ioffset] = zsinw * std::sin(zrlon);
- pxn[mni + id1 * 1 + id2 * 3 + id3 * jd + ioffset] = zcosw * zsgn;
+ pxn[mni + id1 + id2 * 1 + offset] = zsinw * std::cos(zrlon);
+ pxn[mni + id1 + id2 * 2 + offset] = zsinw * std::sin(zrlon);
+ pxn[mni + id1 + id2 * 3 + offset] = zcosw * zsgn;
/* Next initialize the (x,y,z) coordinates for the corner of the */
/* diamond on the same latitude as the (mni,1,,) vertex, which */
/* is (0,mni+1,,) */
- pxn[0 + id1 * (mni + 1) + id2 * 1 + id3 * jd + ioffset] = zsinw * std::cos(zrlon + 2 * pid5);
- pxn[0 + id1 * (mni + 1) + id2 * 2 + id3 * jd + ioffset] = zsinw * std::sin(zrlon + 2 * pid5);
- pxn[0 + id1 * (mni + 1) + id2 * 3 + id3 * jd + ioffset] = zcosw * zsgn;
+ pxn[0 + id1 * (mni + 1) + id2 * 1 + offset] = zsinw * std::cos(zrlon + 2 * pid5);
+ pxn[0 + id1 * (mni + 1) + id2 * 2 + offset] = zsinw * std::sin(zrlon + 2 * pid5);
+ pxn[0 + id1 * (mni + 1) + id2 * 3 + offset] = zcosw * zsgn;
/* Initialize the last diamond vertex, which is located */
/* in the opposite hemisphere as (mni,mni+1,,) */
- pxn[mni + id1 * (mni + 1) + id2 * 1 + id3 * jd + ioffset] = zsinw * std::cos(zrlon + pid5);
- pxn[mni + id1 * (mni + 1) + id2 * 2 + id3 * jd + ioffset] = zsinw * std::sin(zrlon + pid5);
- pxn[mni + id1 * (mni + 1) + id2 * 3 + id3 * jd + ioffset] = -(zcosw * zsgn);
+ pxn[mni + id1 * (mni + 1) + id2 * 1 + offset] = zsinw * std::cos(zrlon + pid5);
+ pxn[mni + id1 * (mni + 1) + id2 * 2 + offset] = zsinw * std::sin(zrlon + pid5);
+ pxn[mni + id1 * (mni + 1) + id2 * 3 + offset] = -(zcosw * zsgn);
/***********************************************************************/
/* First a trisection is performed, if required (kni3=1). */
-
if (kni3 == 1)
{
- ml3 = mni / 3;
+ int ml3 = mni / 3;
- /* Trisect the rows of the diamond. */
+ // Trisect the rows of the diamond.
for (int j1 = 1; j1 <= 2; ++j1)
{
for (int j2 = 1; j2 <= 2; ++j2)
{
- mi1 = (j1 - 1) * mni;
- mi2 = j2 * ml3 + 1;
- zgamma = (double) j2 / 3.0;
+ int mi1 = (j1 - 1) * mni;
+ int mi2 = j2 * ml3 + 1;
+ double zgamma = (double) j2 / 3.0;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, zgamma, mi1, 1, mi1, mni + 1, mi1, mi2);
}
}
- /* Trisect the columns of the diamond. */
+ // Trisect the columns of the diamond.
for (int j1 = 1; j1 <= 2; ++j1)
{
for (int j2 = 1; j2 <= 2; ++j2)
{
- mi1 = j2 * ml3;
- mi2 = (j1 - 1) * mni + 1;
- zgamma = (double) j2 / 3.0;
+ int mi1 = j2 * ml3;
+ int mi2 = (j1 - 1) * mni + 1;
+ double zgamma = (double) j2 / 3.0;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, zgamma, 0, mi2, mni, mi2, mi1, mi2);
}
}
- /* Trisect the diagonal of the diamond. */
+ // Trisect the diagonal of the diamond.
for (int j2 = 1; j2 <= 2; ++j2)
{
- mi1 = mni - j2 * ml3;
- mi2 = j2 * ml3 + 1;
- zgamma = (double) j2 / (float) 3.;
+ int mi1 = mni - j2 * ml3;
+ int mi2 = j2 * ml3 + 1;
+ double zgamma = (double) j2 / 3.0;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, zgamma, mni, 1, 0, mni + 1, mi1, mi2);
}
- /* Compute coordinates of icosahedral triangle centers. */
+ // Compute coordinates of icosahedral triangle centers.
tricntr(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, mni);
}
/***********************************************************************/
/* Find the coordinates of the triangle nodes by iteratively */
/* bisecting the diamond intervals. */
-
for (int jb = 0; jb <= kni2 - 1; ++jb)
{
- mm = pow_ii(3, kni3) * pow_ii(2, jb);
- ml = mni / mm;
- ml2 = ml / 2;
+ int mm = pow_ii(3, kni3) * pow_ii(2, jb);
+ int ml = mni / mm;
+ int ml2 = ml / 2;
// Compute the rows of the diamond.
for (int j1 = 1; j1 <= mm + 1; ++j1)
{
for (int j2 = 1; j2 <= mm; ++j2)
{
- mi1 = (j1 - 1) * ml;
- mi2 = (j2 - 1) * ml + ml2 + 1;
+ int mi1 = (j1 - 1) * ml;
+ int mi2 = (j2 - 1) * ml + ml2 + 1;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, 0.5, mi1, mi2 - ml2, mi1, mi2 + ml2, mi1, mi2);
}
}
@@ -816,8 +805,8 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
{
for (int j2 = 1; j2 <= mm; ++j2)
{
- mi1 = (j2 - 1) * ml + ml2;
- mi2 = (j1 - 1) * ml + 1;
+ int mi1 = (j2 - 1) * ml + ml2;
+ int mi2 = (j1 - 1) * ml + 1;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, 0.5, mi1 - ml2, mi2, mi1 + ml2, mi2, mi1, mi2);
}
}
@@ -827,8 +816,8 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
{
for (int j2 = 1; j2 <= mm; ++j2)
{
- mi1 = (j1 - 1) * ml + ml2;
- mi2 = (j2 - 1) * ml + ml2 + 1;
+ int mi1 = (j1 - 1) * ml + ml2;
+ int mi2 = (j2 - 1) * ml + ml2 + 1;
gcpt(pxn, kig1s, kig1e, kig2s, kig2e, knd, jd, 0.5, mi1 - ml2, mi2 + ml2, mi1 + ml2, mi2 - ml2, mi1, mi2);
}
}
@@ -836,23 +825,14 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
/***********************************************************************/
/* Set pxn to 0 if it is less than 2.5 e-14 to avoid round-off errors */
-
for (int j2 = kig2s; j2 <= kig2e; ++j2)
{
for (int j1 = kig1s; j1 <= kig1e; ++j1)
{
- if (std::fabs(pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset]) < 2.5e-14)
- {
- pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset] = 0.0;
- }
- if (std::fabs(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset]) < 2.5e-14)
- {
- pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset] = 0.0;
- }
- if (std::fabs(pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset]) < 2.5e-14)
- {
- pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset] = 0.0;
- }
+ int j12 = j1 + id1 * j2;
+ if (std::fabs(pxn[j12 + id2 * 1 + offset]) < 2.5e-14) { pxn[j12 + id2 * 1 + offset] = 0.0; }
+ if (std::fabs(pxn[j12 + id2 * 2 + offset]) < 2.5e-14) { pxn[j12 + id2 * 2 + offset] = 0.0; }
+ if (std::fabs(pxn[j12 + id2 * 3 + offset]) < 2.5e-14) { pxn[j12 + id2 * 3 + offset] = 0.0; }
}
}
}
@@ -863,18 +843,17 @@ glo_coor(double *pxn, double *prlon, double *prlat, int kig1s, int kig1e, int ki
for (int jd = 1; jd <= knd; ++jd)
{
+ int offset = id3 * jd + ioffset;
for (int j2 = kig2s; j2 <= kig2e; ++j2)
{
for (int j1 = kig1s; j1 <= kig1e; ++j1)
{
- prlon[j1 + id1 * j2 + id2 * jd + joffset] = std::atan2(pxn[j1 + id1 * j2 + id2 * 2 + id3 * jd + ioffset],
- pxn[j1 + id1 * j2 + id2 * 1 + id3 * jd + ioffset] + 1.0e-20);
- prlat[j1 + id1 * j2 + id2 * jd + joffset] = std::asin(pxn[j1 + id1 * j2 + id2 * 3 + id3 * jd + ioffset]);
+ int j12offset = j1 + id1 * j2 + id2 * jd + joffset;
+ prlon[j12offset] = std::atan2(pxn[j1 + id1 * j2 + id2 * 2 + offset], pxn[j1 + id1 * j2 + id2 * 1 + offset] + 1.0e-20);
+ prlat[j12offset] = std::asin(pxn[j1 + id1 * j2 + id2 * 3 + offset]);
}
}
}
-
- return;
} // glo_coor
/*****************************************************************************/
@@ -1204,7 +1183,7 @@ gme_grid(int withBounds, size_t gridsize, double *rlon, double *rlat, double *bl
exit(-1);
}
- std::vector<double> xn(gridsize * 3);
+ std::vector<double> xn(gridsize * 3, 0.0);
std::vector<double> rlonx((ni + 3) * (ni + 3) * nd);
std::vector<double> rlatx((ni + 3) * (ni + 3) * nd);
@@ -1215,6 +1194,20 @@ gme_grid(int withBounds, size_t gridsize, double *rlon, double *rlat, double *bl
glo_coor(xn.data(), rlon, rlat, im1s, im1e, im2s, im2e, nd, ni2, ni3);
+ // check coordinates
+ size_t inull = 0;
+ for (size_t i = 0; i < gridsize; ++i)
+ {
+ if ((xn[i * 3 + 0] == 0.0 || std::isnan(xn[i * 3 + 0])) && (xn[i * 3 + 1] == 0.0 || std::isnan(xn[i * 3 + 1]))
+ && (xn[i * 3 + 2] == 0.0 || std::isnan(xn[i * 3 + 2])))
+ inull++;
+ }
+ if (inull > gridsize / 4)
+ {
+ fprintf(stderr, "gme_grid: Coordinates calculation of the global GME grid failed (ni=%d)!\n", ni);
+ exit(-1);
+ }
+
xd(rlon, im1s, im1e, im2s, im2e, nd, rlonx.data(), im1s - 1, im1e + 1, im2s - 1, im2e + 1, nd);
xd(rlat, im1s, im1e, im2s, im2e, nd, rlatx.data(), im1s - 1, im1e + 1, im2s - 1, im2e + 1, nd);
--
GitLab
From b7ffc33ee86b1b47978885112b5416cd79811d2b Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Thu, 20 Mar 2025 20:29:27 +0100
Subject: [PATCH 10/10] gme_grid: check if calculation of coordinates failed
---
src/grid_gme.cc | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/src/grid_gme.cc b/src/grid_gme.cc
index 864dc87d6..4696e13b9 100644
--- a/src/grid_gme.cc
+++ b/src/grid_gme.cc
@@ -1198,8 +1198,9 @@ gme_grid(int withBounds, size_t gridsize, double *rlon, double *rlat, double *bl
size_t inull = 0;
for (size_t i = 0; i < gridsize; ++i)
{
- if ((xn[i * 3 + 0] == 0.0 || std::isnan(xn[i * 3 + 0])) && (xn[i * 3 + 1] == 0.0 || std::isnan(xn[i * 3 + 1]))
- && (xn[i * 3 + 2] == 0.0 || std::isnan(xn[i * 3 + 2])))
+ if ((std::fabs(xn[i * 3 + 0]) <= 0 || std::isnan(xn[i * 3 + 0]))
+ && (std::fabs(xn[i * 3 + 1]) <= 0 || std::isnan(xn[i * 3 + 1]))
+ && (std::fabs(xn[i * 3 + 2]) <= 0 || std::isnan(xn[i * 3 + 2])))
inull++;
}
if (inull > gridsize / 4)
--
GitLab