From 33107354deea6766c44311508dbd9a85dad7b61b Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Thu, 27 Mar 2025 11:28:32 +0100
Subject: [PATCH] sortlevel: added NetCDF support

---
 ChangeLog       |   4 ++
 libcdi          |   2 +-
 src/Intlevel.cc |   4 +-
 src/Sort.cc     | 105 ++++++++++++++++++++++++++----------------------
 4 files changed, 65 insertions(+), 50 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 17ce50d70..cd4a172b8 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2025-03-27  Uwe Schulzweida
+
+	* sortlevel: added NetCDF support
+
 2025-03-25  Uwe Schulzweida
 
 	* Input: failed with more than one record (bug fix)
diff --git a/libcdi b/libcdi
index 3b65d721f..1549d6e7f 160000
--- a/libcdi
+++ b/libcdi
@@ -1 +1 @@
-Subproject commit 3b65d721f75aed9057376f2b271ebd129a49e909
+Subproject commit 1549d6e7f8f506b9837e31d6ea3572516e5eda01
diff --git a/src/Intlevel.cc b/src/Intlevel.cc
index 26c4caaf6..55440e61e 100644
--- a/src/Intlevel.cc
+++ b/src/Intlevel.cc
@@ -536,8 +536,8 @@ public:
     for (int index = 0; index < numZaxes; ++index)
       {
         auto zaxisID = vlistZaxis(vlistID1, index);
-        auto nlevels = zaxisInqSize(zaxisID);
-        if (zaxisID == zaxisID1 || nlevels == nlev1) vlistChangeZaxisIndex(vlistID2, index, zaxisID2);
+        auto numLevels = zaxisInqSize(zaxisID);
+        if (zaxisID == zaxisID1 || numLevels == nlev1) vlistChangeZaxisIndex(vlistID2, index, zaxisID2);
       }
 
     varList2 = VarList(vlistID2);
diff --git a/src/Sort.cc b/src/Sort.cc
index 9ffa4e770..cbff75530 100644
--- a/src/Sort.cc
+++ b/src/Sort.cc
@@ -28,7 +28,7 @@ struct LevInfo
 struct VarInfo
 {
   int varID;
-  int nlevels;
+  int numLevels;
   int code;
   std::string param;
   std::string name;
@@ -36,22 +36,22 @@ struct VarInfo
 };
 
 static void
-setNmiss(int varID, int levelID, int nvars, std::vector<VarInfo> &varsInfo, size_t numMissVals)
+setNmiss(int varID, int levelID, int numVars, std::vector<VarInfo> &varsInfo, size_t numMissVals)
 {
-  int vindex, lindex;
+  int varIndex, levIndex;
 
-  for (vindex = 0; vindex < nvars; vindex++)
-    if (varsInfo[vindex].varID == varID) break;
+  for (varIndex = 0; varIndex < numVars; varIndex++)
+    if (varsInfo[varIndex].varID == varID) break;
 
-  if (vindex == nvars) cdo_abort("Internal problem; varID not found!");
+  if (varIndex == numVars) cdo_abort("Internal problem; varID not found!");
 
-  auto nlevels = varsInfo[vindex].nlevels;
-  for (lindex = 0; lindex < nlevels; lindex++)
-    if (varsInfo[vindex].levInfo[lindex].levelID == levelID) break;
+  auto numLevels = varsInfo[varIndex].numLevels;
+  for (levIndex = 0; levIndex < numLevels; levIndex++)
+    if (varsInfo[varIndex].levInfo[levIndex].levelID == levelID) break;
 
-  if (lindex == nlevels) cdo_abort("Internal problem; levelID not found!");
+  if (levIndex == numLevels) cdo_abort("Internal problem; levelID not found!");
 
-  varsInfo[vindex].levInfo[lindex].numMissVals = numMissVals;
+  varsInfo[varIndex].levInfo[levIndex].numMissVals = numMissVals;
 }
 
 class Sort : public Process
@@ -75,7 +75,7 @@ public:
   CdoStreamID streamID2;
   int taxisID1;
   int taxisID2;
-  int nvars;
+  int numVars;
 
   VarList varList1;
   std::vector<VarInfo> varsInfo;
@@ -117,24 +117,40 @@ public:
      else if ( operatorID == SORTLEVEL )
         ;
     */
+    if (operatorID == SORTLEVEL)
+      {
+        auto numZaxes = vlistNumZaxis(vlistID2);
+        for (int index = 0; index < numZaxes; ++index)
+          {
+            auto zaxisID1 = vlistZaxis(vlistID2, index);
+            auto zaxisID2 = zaxisDuplicate(zaxisID1);
+            auto numLevels = zaxisInqSize(zaxisID2);
+            Varray<double> levels(numLevels);
+            cdo_zaxis_inq_levels(zaxisID2, levels.data());
+
+            compareLess ? ranges::sort(levels, ranges::less()) : ranges::sort(levels, ranges::greater());
+
+            zaxisDefLevels(zaxisID2, levels.data());
+            vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
+          }
+      }
 
     streamID2 = cdo_open_write(1);
     cdo_def_vlist(streamID2, vlistID2);
 
     varList1 = VarList(vlistID1);
+    numVars = varList1.numVars();
 
-    nvars = vlistNvars(vlistID1);
-
-    varsInfo = std::vector<VarInfo>(nvars);
-    for (int varID = 0; varID < nvars; ++varID)
+    varsInfo = std::vector<VarInfo>(numVars);
+    for (int varID = 0; varID < numVars; ++varID)
       {
         auto const &var = varList1.vars[varID];
-        varsInfo[varID].nlevels = var.nlevels;
+        varsInfo[varID].numLevels = var.nlevels;
         varsInfo[varID].levInfo.resize(var.nlevels);
       }
 
-    vardata = Varray2D<double>(nvars);
-    for (int varID = 0; varID < nvars; ++varID)
+    vardata = Varray2D<double>(numVars);
+    for (int varID = 0; varID < numVars; ++varID)
       {
         auto const &var = varList1.vars[varID];
         vardata[varID].resize(var.gridsize * var.nlevels);
@@ -175,60 +191,55 @@ public:
             size_t numMissVals;
             cdo_read_field(streamID1, single, &numMissVals);
 
-            setNmiss(varID, levelID, nvars, varsInfo, numMissVals);
+            setNmiss(varID, levelID, numVars, varsInfo, numMissVals);
             // varsInfo[varID].levInfo[levelID].numMissVals = numMissVals;
           }
 
         if (tsID == 0)
           {
             if (Options::cdoVerbose)
-              for (int vindex = 0; vindex < nvars; vindex++)
+              for (int varID = 0; varID < numVars; varID++)
                 {
-                  auto const &varInfo = varsInfo[vindex];
-                  for (int lindex = 0; lindex < varInfo.nlevels; ++lindex)
-                    printf("sort in: %d %s %d %d %g\n", vindex, varInfo.name.c_str(), varInfo.code, varInfo.nlevels,
-                           varInfo.levInfo[lindex].level);
+                  auto const &varInfo = varsInfo[varID];
+                  for (int levelID = 0; levelID < varInfo.numLevels; ++levelID)
+                    printf("sort in: %d %s %d %d %d %g\n", varID, varInfo.name.c_str(), varInfo.code, varInfo.numLevels,
+                           varInfo.levInfo[levelID].levelID, varInfo.levInfo[levelID].level);
                 }
 
-            if (operatorID == SORTCODE)
-              ranges::sort(varsInfo, {}, &VarInfo::code);
-            else if (operatorID == SORTPARAM)
-              ranges::sort(varsInfo, {}, &VarInfo::param);
-            else if (operatorID == SORTNAME)
-              ranges::sort(varsInfo, {}, &VarInfo::name);
+            if (operatorID == SORTCODE) { ranges::sort(varsInfo, {}, &VarInfo::code); }
+            else if (operatorID == SORTPARAM) { ranges::sort(varsInfo, {}, &VarInfo::param); }
+            else if (operatorID == SORTNAME) { ranges::sort(varsInfo, {}, &VarInfo::name); }
             else if (operatorID == SORTLEVEL)
               {
-                for (int vindex = 0; vindex < nvars; vindex++)
+                for (int varID = 0; varID < numVars; varID++)
                   {
-                    if (compareLess)
-                      ranges::sort(varsInfo[vindex].levInfo, ranges::less(), &LevInfo::level);
-                    else
-                      ranges::sort(varsInfo[vindex].levInfo, ranges::greater(), &LevInfo::level);
+                    if (compareLess) { ranges::sort(varsInfo[varID].levInfo, ranges::less(), &LevInfo::level); }
+                    else { ranges::sort(varsInfo[varID].levInfo, ranges::greater(), &LevInfo::level); }
                   }
               }
 
             if (Options::cdoVerbose)
-              for (int vindex = 0; vindex < nvars; vindex++)
+              for (int varID = 0; varID < numVars; varID++)
                 {
-                  auto const &varInfo = varsInfo[vindex];
-                  for (int lindex = 0; lindex < varInfo.nlevels; ++lindex)
-                    printf("sort out: %d %s %d %d %g\n", vindex, varInfo.name.c_str(), varInfo.code, varInfo.nlevels,
-                           varInfo.levInfo[lindex].level);
+                  auto const &varInfo = varsInfo[varID];
+                  for (int levelID = 0; levelID < varInfo.numLevels; ++levelID)
+                    printf("sort out: %d %s %d %d %d %g\n", varID, varInfo.name.c_str(), varInfo.code, varInfo.numLevels,
+                           varInfo.levInfo[levelID].levelID, varInfo.levInfo[levelID].level);
                 }
           }
 
-        for (int vindex = 0; vindex < nvars; vindex++)
+        for (int varID = 0; varID < numVars; varID++)
           {
-            auto const &varInfo = varsInfo[vindex];
+            auto const &varInfo = varsInfo[varID];
             auto const &var = varList1.vars[varInfo.varID];
-            for (int lindex = 0; lindex < var.nlevels; ++lindex)
+            for (int levelID = 0; levelID < var.nlevels; ++levelID)
               {
-                auto levelID = varInfo.levInfo[lindex].levelID;
-                auto numMissVals = varInfo.levInfo[lindex].numMissVals;
+                auto levelID2 = varInfo.levInfo[levelID].levelID;
+                auto numMissVals = varInfo.levInfo[levelID].numMissVals;
 
                 if (tsID == 0 || !var.isConstant)
                   {
-                    auto offset = var.gridsize * levelID;
+                    auto offset = var.gridsize * levelID2;
                     auto single = &vardata[varInfo.varID][offset];
 
                     cdo_def_field(streamID2, varInfo.varID, levelID);
-- 
GitLab