diff --git a/ChangeLog b/ChangeLog
index 211f29c2866fbb7919698a59996556e1ff3e04d4..4ba6c738ac8142af457348c19b19a9bd5bd4c104 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 a7ddc92a34a1d8c866d008b8f38e1429d0ac67a3..864dc87d6a0308d0ebb23ab0b429323d80a8374b 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);