Add function read_target_coordinate()

src/Intlevel3d.cc

@@ -92,10 +92,11 @@ read_source_coordinate(int streamNumber, VarList &varList2, Varray<float> &zleve
 
   if (vlistNvars(vlistID2) != 1) cdo_abort("infile2: Only one single variable is allowed!");
 
-  auto gridsize = varList2[0].gridsize;  // horizontal gridsize of input z coordinate
-  auto nlev = varList2[0].nlevels;       // number of input levels for later use
+  auto gridsize = varList2[0].gridsize;
+  auto nlevels = varList2[0].nlevels;
+
+  zlevelsIn.resize(gridsize * nlevels);
 
-  zlevelsIn.resize(gridsize * nlev);
   auto nrecs = cdo_stream_inq_timestep(streamID2, 0);
   if (Options::cdoVerbose) cdo_print("%d records input 3d vertical height", nrecs);
 
@@ -112,6 +113,37 @@ read_source_coordinate(int streamNumber, VarList &varList2, Varray<float> &zleve
   cdo_stream_close(streamID2);
 }
 
+static void
+read_target_coordinate(const std::string &fileName, VarList &varList0, Varray<float> &zlevelsOut)
+{
+  auto streamID0 = stream_open_read_locked(fileName.c_str());  // 3d vertical target coordinate
+  auto vlistID0 = streamInqVlist(streamID0);
+
+  varListInit(varList0, vlistID0);
+
+  if (vlistNvars(vlistID0) != 1) cdo_abort("tgtcoordinate: Only one single variable is allowed!");
+
+  auto gridsize = varList0[0].gridsize;
+  auto nlevels = varList0[0].nlevels;
+
+  zlevelsOut.resize(gridsize * nlevels);
+
+  auto nrecs = streamInqTimestep(streamID0, 0);
+  if (Options::cdoVerbose) cdo_print("%d records target 3d vertical height and gridsize %zu", nrecs, gridsize);
+
+  for (int recID = 0; recID < nrecs; ++recID)
+    {
+      int varID, levelID;
+      streamInqRecord(streamID0, &varID, &levelID);
+      auto offset = gridsize * levelID;
+      size_t nmiss;
+      streamReadRecordF(streamID0, &zlevelsOut[offset], &nmiss);
+      if (0 != nmiss) cdo_abort("Output vertical coordinate variables are not allowd to contain missing values.");
+    }
+
+  streamClose(streamID0);
+}
+
 void *
 Intlevel3d(void *process)
 {
@@ -157,46 +189,15 @@ Intlevel3d(void *process)
   auto nlevi = varList2[0].nlevels;       // number of input levels for later use
 
   // Read 3d target coordinate (streamID0)
-  auto streamID0 = stream_open_read_locked(cdo_operator_argv(0).c_str());  // 3d vertical target coordinate
-  auto vlistID0 = streamInqVlist(streamID0);
-  if (vlistNvars(vlistID0) != 1) cdo_abort("tgtcoordinate: Only one single variable is allowed!");
-
   VarList varList0;
-  varListInit(varList0, vlistID0);
-
   Varray<float> zlevelsOut;
 
-  size_t gridsizeo;
-  int nlevo;
-  int gridID3 = -1;
-  int zaxisID0 = -1;
+  read_target_coordinate(cdo_operator_argv(0), varList0, zlevelsOut);
 
-  {
-    auto gridID = varList0[0].gridID;
-    gridID3 = gridID;
-    auto zaxisID = varList0[0].zaxisID;
-    auto gridsize = varList0[0].gridsize;
-    auto nlevel = varList0[0].nlevels;
-
-    zlevelsOut.resize(gridsize * nlevel);
-    zaxisID0 = zaxisID;
-    nlevo = nlevel;        // number of output levels for later use
-    gridsizeo = gridsize;  // horizontal gridsize of output z coordinate
-    auto nrecs = streamInqTimestep(streamID0, 0);
-    if (Options::cdoVerbose) cdo_print("%d records target 3d vertical height and gridsize %zu", nrecs, gridsize);
-
-    for (int recID = 0; recID < nrecs; ++recID)
-      {
-        int varID, levelID;
-        streamInqRecord(streamID0, &varID, &levelID);
-        auto offset = gridsize * levelID;
-        size_t nmiss;
-        streamReadRecordF(streamID0, &zlevelsOut[offset], &nmiss);
-        if (0 != nmiss) cdo_abort("Output vertical coordinate variables are not allowd to contain missing values.");
-      }
-
-    streamClose(streamID0);
-  }
+  auto gridsizeo = varList0[0].gridsize;  // horizontal gridsize of output z coordinate
+  auto nlevo = varList0[0].nlevels;       // number of output levels for later use
+  auto gridID3 = varList0[0].gridID;
+  auto zaxisID3 = varList0[0].zaxisID;
 
   // gridsize of input and output vertical coordinate must be equal (later use of gridsizeo ONLY)
   if (gridsizei != gridsizeo) cdo_abort("Input and output vertical coordinate must have the same gridsize!");
@@ -238,7 +239,6 @@ Intlevel3d(void *process)
   vert_gen_weights3d(expol, gridSize, nlevi, zlevelsIn, nlevo, zlevelsOut, lev_idx, lev_wgt);
   varray_free(zlevelsIn);
 
-  auto zaxisID3 = zaxisID0;
   for (int index = 0; index < nzaxis; ++index)
     if (zaxisID1 == vlistZaxis(vlistID1, index)) vlistChangeZaxisIndex(vlistID3, index, zaxisID3);
 
@@ -256,14 +256,11 @@ Intlevel3d(void *process)
   auto maxlev = std::max(nlevi, nlevo);
   auto nvars = vlistNvars(vlistID1);
 
-  std::vector<bool> vars(nvars);
-  std::vector<bool> varinterp(nvars);                // marker for variables to be interpolated
+  std::vector<bool> processVars(nvars);
+  std::vector<bool> varinterp(nvars, false);         // marker for variables to be interpolated
   std::vector<std::vector<size_t>> varnmiss(nvars);  // can for missing values of arbitrary variables
   Field3DVector vardata1(nvars), vardata2(nvars);
 
-  // by default no variable should be interpolated
-  for (int varID = 0; varID < nvars; ++varID) varinterp[varID] = false;
-
   for (int varID = 0; varID < nvars; ++varID)
     {
       auto zaxisID = varList1[varID].zaxisID;
@@ -304,7 +301,7 @@ Intlevel3d(void *process)
       auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
       if (nrecs == 0) break;
 
-      for (int varID = 0; varID < nvars; ++varID) vars[varID] = false;
+      for (int varID = 0; varID < nvars; ++varID) processVars[varID] = false;
 
       cdo_taxis_copy_timestep(taxisID3, taxisID1);
       cdo_def_timestep(streamID3, tsID);
@@ -315,13 +312,13 @@ Intlevel3d(void *process)
           int varID, levelID;
           cdo_inq_record(streamID1, &varID, &levelID);
           cdo_read_record(streamID1, vardata1[varID], levelID, &varnmiss[varID][levelID]);
-          vars[varID] = true;
+          processVars[varID] = true;
         }
 
       // Perform the interpolation on all valid data variables
       for (int varID = 0; varID < nvars; ++varID)
         {
-          if (vars[varID] && varinterp[varID])
+          if (processVars[varID] && varinterp[varID])
             {
               auto gridsize = varList1[varID].gridsize;
               auto missval = varList1[varID].missval;
@@ -346,7 +343,7 @@ Intlevel3d(void *process)
       // write the output
       for (int varID = 0; varID < nvars; ++varID)
         {
-          if (vars[varID])
+          if (processVars[varID])
             {
               for (int levelID = 0; levelID < varList3[varID].nlevels; ++levelID)
                 {