Add function read_source_coordinate()

src/Intlevel3d.cc

@@ -37,7 +37,7 @@ bool levelDirDown(const int nlev, const double *const lev);
  * 3d version of vert_gen_weights() (src/Intlevel.cc)
  */
 static void
-vert_gen_weights3d(bool expol, size_t gridsize, int nlev1, const Varray<float> &xlev1, int nlev2, const Varray<double> &xlev2,
+vert_gen_weights3d(bool expol, size_t gridsize, int nlev1, const Varray<float> &xlev1, int nlev2, const Varray<float> &xlev2,
                    Varray<int> &xlev_idx, Varray<float> &xlev_wgt)
 {
   auto nthreads = Threading::ompNumThreads;
@@ -82,6 +82,36 @@ vlist_copy_var_attributes(const VarList &varList0, int varID0, int vlistID3, int
   if (var0.units.size()) cdiDefKeyString(vlistID3, oz3dvarID, CDI_KEY_UNITS, var0.units.c_str());
 }
 
+static void
+read_source_coordinate(int streamNumber, VarList &varList2, Varray<float> &zlevelsIn)
+{
+  auto streamID2 = cdo_open_read(streamNumber);  // 3d vertical source coordinate
+  auto vlistID2 = cdo_stream_inq_vlist(streamID2);
+
+  varListInit(varList2, vlistID2);
+
+  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
+
+  zlevelsIn.resize(gridsize * nlev);
+  auto nrecs = cdo_stream_inq_timestep(streamID2, 0);
+  if (Options::cdoVerbose) cdo_print("%d records input 3d vertical height", nrecs);
+
+  for (int recID = 0; recID < nrecs; ++recID)
+    {
+      int varID, levelID;
+      cdo_inq_record(streamID2, &varID, &levelID);
+      auto offset = gridsize * levelID;
+      size_t nmiss;
+      cdo_read_record_f(streamID2, &zlevelsIn[offset], &nmiss);
+      if (0 != nmiss) cdo_abort("Input vertical coordinate variables are not allowed to contain missing values.");
+    }
+
+  cdo_stream_close(streamID2);
+}
+
 void *
 Intlevel3d(void *process)
 {
@@ -118,36 +148,13 @@ Intlevel3d(void *process)
   auto memType = varList1[0].memType;
 
   // Read 3d source coordinate (streamID2)
-  size_t gridsizei;
-  int nlevi;
-  Varray<float> zlevels_in;
-  size_t zlevels_in_miss = 0;
+  VarList varList2;
+  Varray<float> zlevelsIn;
 
-  {
-    auto streamID2 = cdo_open_read(1);  // 3d vertical source coordinate
-    auto vlistID2 = cdo_stream_inq_vlist(streamID2);
-
-    if (vlistNvars(vlistID2) != 1) cdo_abort("infile2: Only one single variable is allowed!");
-
-    auto gridID = vlistInqVarGrid(vlistID2, 0);
-    auto zaxisID = vlistInqVarZaxis(vlistID2, 0);
-    gridsizei = gridInqSize(gridID);  // horizontal gridsize of input z coordinate
-    nlevi = zaxisInqSize(zaxisID);    // number of input levels for later use
+  read_source_coordinate(1, varList2, zlevelsIn);
 
-    zlevels_in.resize(gridsizei * nlevi);
-    auto nrecs = cdo_stream_inq_timestep(streamID2, 0);
-    if (Options::cdoVerbose) cdo_print("%d records input 3d vertical height", nrecs);
-
-    for (int recID = 0; recID < nrecs; ++recID)
-      {
-        int varID, levelID;
-        cdo_inq_record(streamID2, &varID, &levelID);
-        auto offset = gridsizei * levelID;
-        cdo_read_record_f(streamID2, &zlevels_in[offset], &zlevels_in_miss);
-      }
-
-    cdo_stream_close(streamID2);
-  }
+  auto gridsizei = varList2[0].gridsize;  // horizontal gridsize of input z coordinate
+  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
@@ -157,7 +164,7 @@ Intlevel3d(void *process)
   VarList varList0;
   varListInit(varList0, vlistID0);
 
-  Varray<double> zlevelsOut;
+  Varray<float> zlevelsOut;
 
   size_t gridsizeo;
   int nlevo;
@@ -184,16 +191,13 @@ Intlevel3d(void *process)
         streamInqRecord(streamID0, &varID, &levelID);
         auto offset = gridsize * levelID;
         size_t nmiss;
-        streamReadRecord(streamID0, &zlevelsOut[offset], &nmiss);
+        streamReadRecordF(streamID0, &zlevelsOut[offset], &nmiss);
         if (0 != nmiss) cdo_abort("Output vertical coordinate variables are not allowd to contain missing values.");
       }
 
     streamClose(streamID0);
   }
 
-  // Missing values are not allowed for coordinate variables
-  if (0 != zlevels_in_miss) cdo_abort("Input vertical coordinate variables are not allowed to contain missing values.");
-
   // 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!");
 
@@ -231,8 +235,8 @@ Intlevel3d(void *process)
   Varray<int> lev_idx(nlevo * gridSize);
   Varray<float> lev_wgt(nlevo * gridSize);
 
-  vert_gen_weights3d(expol, gridSize, nlevi, zlevels_in, nlevo, zlevelsOut, lev_idx, lev_wgt);
-  varray_free(zlevels_in);
+  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)
@@ -358,7 +362,7 @@ Intlevel3d(void *process)
         {
           auto offset = varList3[oz3dvarID].gridsize * levelID;
           cdo_def_record(streamID3, oz3dvarID, levelID);
-          cdo_write_record(streamID3, &zlevelsOut[offset], 0);
+          cdo_write_record_f(streamID3, &zlevelsOut[offset], 0);
         }
 
       tsID++;