change_hybrid_zaxis: check number of levels.

src/Vertintml.cc

@@ -39,7 +39,7 @@ zaxis_is_hybrid(int zaxistype)
 }
 
 static void
-change_hybrid_zaxis(int vlistID1, int vlistID2, int nvct, double *vct, int zaxisID2)
+change_hybrid_zaxis(int vlistID1, int vlistID2, int nvct, double *vct, int zaxisID2, int nhlevf, int nhlevh)
 {
   int nzaxis = vlistNzaxis(vlistID1);
   for (int iz = 0; iz < nzaxis; ++iz)
@@ -48,7 +48,7 @@ change_hybrid_zaxis(int vlistID1, int vlistID2, int nvct, double *vct, int zaxis
       int nlevel = zaxisInqSize(zaxisID);
       int zaxistype = zaxisInqType(zaxisID);
 
-      if (zaxis_is_hybrid(zaxistype) && nlevel > 1)
+      if (zaxis_is_hybrid(zaxistype) && (nlevel == nhlevh || nlevel == nhlevf))
         {
           int nvct2 = zaxisInqVctSize(zaxisID);
           if (nvct2 == nvct && memcmp(vct, zaxisInqVctPtr(zaxisID), nvct * sizeof(double)) == 0)
@@ -163,7 +163,7 @@ Vertintml(void *process)
   int nhlev = 0, nhlevf = 0, nhlevh = 0;
   double *vct = vlist_read_vct(vlistID1, &zaxisIDh, &nvct, &nhlev, &nhlevf, &nhlevh);
 
-  change_hybrid_zaxis(vlistID1, vlistID2, nvct, vct, zaxisIDp);
+  change_hybrid_zaxis(vlistID1, vlistID2, nvct, vct, zaxisIDp, nhlevf, nhlevh);
 
   int psvarID = -1;
   bool linvertvct = false;
@@ -293,8 +293,7 @@ Vertintml(void *process)
           //  KNMI: HIRLAM model version 7.2 uses tableNum=1    (LAMH_D11*)
           //  KNMI: HARMONIE model version 36 uses tableNum=1   (grib*)
           //  (opreational NWP version) KNMI: HARMONIE model version 38 uses
-          //  tableNum=253 (grib,grib_md) and tableNum=1 (grib_sfx) (research
-          //  version)
+          //  tableNum=253 (grib,grib_md) and tableNum=1 (grib_sfx) (research version)
           else if (tableNum == 1 || tableNum == 253)
             {
               mode = ModelMode::HIRLAM;
@@ -326,45 +325,36 @@ Vertintml(void *process)
           if (code == -1)
             {
               //                                       ECHAM                         ECMWF
-              if (sgeopotID == -1 && (strcmp(varname, "geosp") == 0 || strcmp(varname, "z") == 0))
-                code = gribcodes.geopot;
-              else if (tempID == -1 && (strcmp(varname, "st") == 0 || strcmp(varname, "t") == 0))
-                code = gribcodes.temp;
-              else if (psID == -1 && (strcmp(varname, "aps") == 0 || strcmp(varname, "sp") == 0))
-                code = gribcodes.ps;
-              else if (lnpsID == -1 && (strcmp(varname, "lsp") == 0 || strcmp(varname, "lnsp") == 0))
-                code = gribcodes.lsp;
-              else if (geopotID == -1 && strcmp(stdname, "geopotential_full") == 0)
-                code = gribcodes.geopot;
+              // clang-format off
+              if      (sgeopotID == -1 && (strcmp(varname, "geosp") == 0 || strcmp(varname, "z") == 0)) code = gribcodes.geopot;
+              else if (tempID == -1 && (strcmp(varname, "st") == 0 || strcmp(varname, "t") == 0)) code = gribcodes.temp;
+              else if (psID == -1 && (strcmp(varname, "aps") == 0 || strcmp(varname, "sp") == 0)) code = gribcodes.ps;
+              else if (lnpsID == -1 && (strcmp(varname, "lsp") == 0 || strcmp(varname, "lnsp") == 0)) code = gribcodes.lsp;
+              else if (geopotID == -1 && strcmp(stdname, "geopotential_full") == 0) code = gribcodes.geopot;
               // else if ( strcmp(varname, "geopoth") == 0 ) code = 156;
+              // clang-format on
             }
         }
 
       if (mode == ModelMode::ECHAM)
         {
-          if (code == gribcodes.geopot && nlevel == 1)
-            sgeopotID = varID;
-          else if (code == gribcodes.geopot && nlevel == nhlevf)
-            geopotID = varID;
-          else if (code == gribcodes.temp && nlevel == nhlevf)
-            tempID = varID;
-          else if (code == gribcodes.ps && nlevel == 1)
-            psID = varID;
-          else if (code == gribcodes.lsp && nlevel == 1)
-            lnpsID = varID;
-          else if (code == gribcodes.gheight && nlevel == nhlevf)
-            gheightID = varID;
+          // clang-format off
+          if      (code == gribcodes.geopot && nlevel == 1) sgeopotID = varID;
+          else if (code == gribcodes.geopot && nlevel == nhlevf) geopotID = varID;
+          else if (code == gribcodes.temp && nlevel == nhlevf) tempID = varID;
+          else if (code == gribcodes.ps && nlevel == 1) psID = varID;
+          else if (code == gribcodes.lsp && nlevel == 1) lnpsID = varID;
+          else if (code == gribcodes.gheight && nlevel == nhlevf) gheightID = varID;
+          // clang-format on
         }
       else if (mode == ModelMode::WMO || mode == ModelMode::HIRLAM)
         {
-          if (code == gribcodes.geopot && nlevel == 1)
-            sgeopotID = varID;
-          else if (code == gribcodes.geopot && nlevel == nhlevf)
-            geopotID = varID;
-          else if (code == gribcodes.temp && nlevel == nhlevf)
-            tempID = varID;
-          else if (code == gribcodes.ps && nlevel == 1)
-            psID = varID;
+          // clang-format off
+          if      (code == gribcodes.geopot && nlevel == 1) sgeopotID = varID;
+          else if (code == gribcodes.geopot && nlevel == nhlevf) geopotID = varID;
+          else if (code == gribcodes.temp && nlevel == nhlevf) tempID = varID;
+          else if (code == gribcodes.ps && nlevel == 1) psID = varID;
+          // clang-format on
         }
 
       if (gridInqType(gridID) == GRID_SPECTRAL && zaxis_is_hybrid(zaxistype)) cdoAbort("Spectral data on model level unsupported!");
@@ -376,8 +366,7 @@ Vertintml(void *process)
       else
         vardata1[varID] = (double *) Malloc(gridsize * nlevel * sizeof(double));
 
-      /* if ( zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevel == nhlev )
-       */
+      // if ( zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevel == nhlev )
       if (zaxisID == zaxisIDh || (zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && (nlevel == nhlevh || nlevel == nhlevf)))
         {
           varinterp[varID] = true;