diff --git a/src/Timstat2.cc b/src/Timstat2.cc
index fb6a10024fbcdc4ba60bff8233af4ad55c278b58..c10304d97248decad10a9ca010248c926d49fafc 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 377325442823f3135d78ccd70f84936c511394cf..8fd50db07f21d90f409a85dc06458e8c5cf4f250 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 8d50127c6c59c334b13c80a100fc52cbc61607f8..a7ddc92a34a1d8c866d008b8f38e1429d0ac67a3 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 1759b509d6d356c47c37f16803934edbbb3f3728..435eabb18a51907b0076b959f98fe59192f2e879 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 2fd9db8e2d68fd7fdaa5ce4d69d1cae1977a042f..d85433e0cad7c6be817b17555ca78ee295f0b334 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 13e17dbc86911464902aad4cae191abdb87fa62f..1965812e9ae57c47fab638d2a5684f3f7d874525 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 a7c3fc5081f77284d3a07e6423e535c74791822a..711d4ae64d8c5f9894af8d83a2c48fa645632897 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 3e13208e98ac04fe5c0244eb94b7264d75823c42..c713459c684c98382aad3ef0b09f6e3659dba3c9 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