Detrend: failed if missing_value is between 0 and numSteps (bug fix)

ChangeLog

@@ -4,6 +4,10 @@
 	* Using YAC version 3.4.0
 	* Version 2.5.1 release
 
+2025-03-04  Uwe Schulzweida
+
+	* Detrend: failed if missing_value is between 0 and numSteps (bug fix)
+
 2025-03-03  Uwe Schulzweida
 
 	* New operator: air_density

libcdi @ 8dbf0a5b

-Subproject commit 34bec0e7aae27fe52d28d879baedd052e7e437fe
+Subproject commit 8dbf0a5b3427ca421f488041bbfa032e2e4d21dd

src/Detrend.cc

@@ -124,9 +124,9 @@ public:
             auto &paramB = work[1][varID][levelID].vec_d;
             auto &sumj = work[0][varID][levelID].vec_d;
             auto &sumjj = work[1][varID][levelID].vec_d;
-            const auto &sumjx = work[2][varID][levelID].vec_d;
-            const auto &sumx = work[3][varID][levelID].vec_d;
-            const auto &zn = work[4][varID][levelID].vec_d;
+            auto const &sumjx = work[2][varID][levelID].vec_d;
+            auto const &sumx = work[3][varID][levelID].vec_d;
+            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]));
@@ -151,13 +151,13 @@ public:
     auto numVars = varList.numVars();
     for (int varID = 0; varID < numVars; ++varID)
       {
-        const auto &var = varList.vars[varID];
+        auto const &var = varList.vars[varID];
         if (var.isConstant) continue;
         for (int levelID = 0; levelID < var.nlevels; ++levelID)
           {
             auto &field = varsData[varID][levelID];
-            const auto &paramA = work[0][varID][levelID];
-            const auto &paramB = work[1][varID][levelID];
+            auto const &paramA = work[0][varID][levelID];
+            auto const &paramB = work[1][varID][levelID];
             sub_trend(zj, field, paramA, paramB);
           }
       }

src/Trendarith.cc

@@ -16,6 +16,7 @@
 
 #include "process_int.h"
 #include "cdo_vlist.h"
+#include "cdo_omp.h"
 #include "datetime.h"
 #include "pmlist.h"
 #include "param_conversion.h"
@@ -30,12 +31,18 @@ add_trend(double zj, Varray<T> &v1, const Varray<double> &v2, const Varray<doubl
   auto missval1 = mv;
   auto missval2 = mv;
 
-  auto add_kernel = [&](auto i, auto is_EQ) { return ADDM(v1[i], ADDM(v2[i], MULM(v3[i], zj))); };
+  auto add_kernel = [&](auto is_EQ) {
+#ifdef _OPENMP
+#pragma omp parallel for if (n > cdoMinLoopSize) default(shared) schedule(static)
+#endif
+    for (size_t i = 0; i < n; ++i)
+      {
+        auto tmp = (is_EQ(v2[i], mv) || is_EQ(v3[i], mv)) ? mv : (v2[i] + v3[i] * zj);
+        v1[i] = ADDM(v1[i], tmp);
+      }
+  };
 
-  if (std::isnan(missval1))
-    for (size_t i = 0; i < n; ++i) v1[i] = add_kernel(i, fp_is_equal);
-  else
-    for (size_t i = 0; i < n; ++i) v1[i] = add_kernel(i, is_equal);
+  std::isnan(mv) ? add_kernel(fp_is_equal) : add_kernel(is_equal);
 }
 
 static void
@@ -52,12 +59,18 @@ sub_trend(double zj, Varray<T> &v1, const Varray<double> &v2, const Varray<doubl
   auto missval1 = mv;
   auto missval2 = mv;
 
-  auto sub_kernel = [&](auto i, auto is_EQ) { return SUBM(v1[i], ADDM(v2[i], MULM(v3[i], zj))); };
+  auto sub_kernel = [&](auto is_EQ) {
+#ifdef _OPENMP
+#pragma omp parallel for if (n > cdoMinLoopSize) default(shared) schedule(static)
+#endif
+    for (size_t i = 0; i < n; ++i)
+      {
+        auto tmp = (is_EQ(v2[i], mv) || is_EQ(v3[i], mv)) ? mv : (v2[i] + v3[i] * zj);
+        v1[i] = SUBM(v1[i], tmp);
+      }
+  };
 
-  if (std::isnan(missval1))
-    for (size_t i = 0; i < n; ++i) v1[i] = sub_kernel(i, fp_is_equal);
-  else
-    for (size_t i = 0; i < n; ++i) v1[i] = sub_kernel(i, is_equal);
+  std::isnan(mv) ? sub_kernel(fp_is_equal) : sub_kernel(is_equal);
 }
 
 static void
@@ -125,7 +138,6 @@ public:
   bool tstepIsEqual;
 
   VarList varList1;
-  FieldVector2D vars2, vars3;
 
 public:
   void
@@ -161,7 +173,12 @@ public:
 
     streamID4 = cdo_open_write(3);
     cdo_def_vlist(streamID4, vlistID4);
+  }
 
+  void
+  run() override
+  {
+    FieldVector2D vars2, vars3;
     field2D_init(vars2, varList1, FIELD_VEC);
     field2D_init(vars3, varList1, FIELD_VEC);
 
@@ -182,11 +199,7 @@ public:
           cdo_read_field(streamID3, vars3[varID][levelID]);
         }
     }
-  }
 
-  void
-  run() override
-  {
     auto calendar = taxisInqCalendar(taxisID1);
     CheckTimeIncr checkTimeIncr;
     JulianDate julianDate0;

src/field_trend.cc

@@ -64,22 +64,18 @@ sub_trend(double zj, Varray<T> &v1, const Varray<double> &v2, const Varray<doubl
   auto missval1 = mv;
   auto missval2 = mv;
 
-  auto sub_kernel = [&](auto i, auto is_EQ) { return SUBM(v1[i], ADDM(v2[i], MULM(v3[i], zj))); };
-
-  if (std::isnan(missval1))
-    {
-#ifdef _OPENMP
-#pragma omp parallel for if (len > cdoMinLoopSize) default(shared) schedule(static)
-#endif
-      for (size_t i = 0; i < len; ++i) { v1[i] = sub_kernel(i, fp_is_equal); }
-    }
-  else
-    {
+  auto sub_kernel = [&](auto is_EQ) {
 #ifdef _OPENMP
 #pragma omp parallel for if (len > cdoMinLoopSize) default(shared) schedule(static)
 #endif
-      for (size_t i = 0; i < len; ++i) { v1[i] = sub_kernel(i, is_equal); }
-    }
+    for (size_t i = 0; i < len; ++i)
+      {
+        auto tmp = (is_EQ(v2[i], mv) || is_EQ(v3[i], mv)) ? mv : (v2[i] + v3[i] * zj);
+        v1[i] = SUBM(v1[i], tmp);
+      }
+  };
+
+  std::isnan(mv) ? sub_kernel(fp_is_equal) : sub_kernel(is_equal);
 }
 
 void