Merge develop

ChangeLog

+2022-09-30  Uwe Schulzweida
+
+	* runpctl: OpenMP parallelized
+
+2022-09-29  Uwe Schulzweida
+
+	* mastrfu: check pressure level orientation
+
 2022-09-28  Uwe Schulzweida
 
-	* Added predefined healpix grid hl<nsize>[b][_xy|_ring|_nested]
+	* Added predefined healpix grid hl<nsize>[b][_nest|_ring|_xy]
 
 2022-09-17  Uwe Schulzweida
 

libcdi @ da013596

-Subproject commit af72d9acd1be7379f35863aa4ddd7c9beb9c773a
+Subproject commit da0135962a79197bfc3ae4b1c71c9a8a165d34d4

src/Mastrfu.cc

@@ -29,6 +29,10 @@ mastrfu(int gridID, int zaxisID, const Varray2D<double> &field1, Varray2D<double
 
   cdo_zaxis_inq_levels(zaxisID, plevel.data());
 
+  const auto plevel1 = plevel[0];
+  const auto plevelN = plevel[nlev-1];
+  if (plevel1 < plevelN) cdo_abort("The 3d pressure level data is upside down! Use the operator invertlev to invert the levels.");
+
   gridInqYvals(gridID, phi.data());
 
   char units[CDI_MAX_NAME];
@@ -60,7 +64,7 @@ mastrfu(int gridID, int zaxisID, const Varray2D<double> &field1, Varray2D<double
         for (int ilev = nlev - 1; ilev >= 0; ilev--)
           for (int n = ilev; n < nlev - 1; ++n)
             {
-              if (DBL_IS_EQUAL(field1[n][ilat], missval))
+              if (dbl_is_equal(field1[n][ilat], missval) || dbl_is_equal(field1[n + 1][ilat], missval))
                 {
                   field2[ilev][ilat] = missval;
                   break;
@@ -149,6 +153,7 @@ Mastrfu(void *process)
           int varID = 0;
           int levelID = recID;
           cdo_def_record(streamID2, varID, levelID);
+          nmiss = array_num_mv(nlat, array2[levelID].data(), missval);
           cdo_write_record(streamID2, array2[levelID].data(), nmiss);
         }
 

src/Runpctl.cc

@@ -19,6 +19,63 @@
 #include "percentiles.h"
 #include "datetime.h"
 #include "field_functions.h"
+#include "cimdOmp.h"
+
+template <typename T>
+static size_t
+runpctl(double pn, int ndates, size_t gridsize, Varray<T> &v2, T missval, const FieldVector3D &vars1, int varID, int levelID, MemType memType)
+{
+  size_t nmiss = 0;
+  Varray2D<T> array_2D(Threading::ompNumThreads, Varray<T>(ndates));
+
+#ifdef _OPENMP
+#pragma omp parallel for default(shared) schedule(dynamic)
+#endif
+  for (size_t i = 0; i < gridsize; ++i)
+    {
+      const auto ompthID = cdo_omp_get_thread_num();
+      auto &array = array_2D[ompthID];
+
+      int j = 0;
+      if (memType == MemType::Float)
+        {
+          for (int inp = 0; inp < ndates; ++inp)
+            {
+              const auto val = vars1[inp][varID][levelID].vec_f[i];
+              if (!dbl_is_equal(val, missval)) array[j++] = val;
+            }
+        }
+      else
+        {
+          for (int inp = 0; inp < ndates; ++inp)
+            {
+              const auto val = vars1[inp][varID][levelID].vec_d[i];
+              if (!dbl_is_equal(val, missval)) array[j++] = val;
+            }
+        }
+
+      if (j > 0)
+        {
+          v2[i] = percentile(array.data(), j, pn);
+        }
+      else
+        {
+          v2[i] = missval;
+          nmiss++;
+        }
+    }
+
+  return nmiss;
+}
+
+static void
+runpctl(double pn, int ndates, Field &field1, const FieldVector3D &vars1, int varID, int levelID)
+{
+  if (field1.memType == MemType::Float)
+    field1.nmiss = runpctl(pn, ndates, field1.gridsize, field1.vec_f, (float)field1.missval, vars1, varID, levelID, field1.memType);
+  else
+    field1.nmiss = runpctl(pn, ndates, field1.gridsize, field1.vec_d, field1.missval, vars1, varID, levelID, field1.memType);
+}
 
 void *
 Runpctl(void *process)
@@ -96,59 +153,11 @@ Runpctl(void *process)
         {
           if (varList1[varID].timetype == TIME_CONSTANT) continue;
 
-          const auto gridsize = varList1[varID].gridsize;
           const auto nlevels = varList1[varID].nlevels;
-          const auto missval = varList1[varID].missval;
-
           for (int levelID = 0; levelID < nlevels; ++levelID)
             {
               auto &field1 = vars1[0][varID][levelID];
-              size_t nmiss = 0;
-              if (field1.memType == MemType::Float)
-                {
-                  for (size_t i = 0; i < gridsize; ++i)
-                    {
-                      int j = 0;
-                      for (int inp = 0; inp < ndates; ++inp)
-                        {
-                          const auto val = vars1[inp][varID][levelID].vec_f[i];
-                          if (!DBL_IS_EQUAL(val, (float) missval)) array_f[j++] = val;
-                        }
-
-                      if (j > 0)
-                        {
-                          field1.vec_f[i] = percentile(array_f.data(), j, pn);
-                        }
-                      else
-                        {
-                          field1.vec_f[i] = missval;
-                          nmiss++;
-                        }
-                    }
-                }
-              else
-                {
-                  for (size_t i = 0; i < gridsize; ++i)
-                    {
-                      int j = 0;
-                      for (int inp = 0; inp < ndates; ++inp)
-                        {
-                          const auto val = vars1[inp][varID][levelID].vec_d[i];
-                          if (!DBL_IS_EQUAL(val, missval)) array_d[j++] = val;
-                        }
-
-                      if (j > 0)
-                        {
-                          field1.vec_d[i] = percentile(array_d.data(), j, pn);
-                        }
-                      else
-                        {
-                          field1.vec_d[i] = missval;
-                          nmiss++;
-                        }
-                    }
-                }
-              field1.nmiss = nmiss;
+              runpctl(pn, ndates, field1, vars1, varID, levelID);
             }
         }