Setgrid.cc 17 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1 2 3 4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6 7 8 9 10 11 12 13 14 15 16 17
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

18 19 20 21 22
/*
   This module contains the following operators:

      Setgrid    setgrid         Set grid
      Setgrid    setgridtype     Set grid type
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23 24
      Setgrid    setgridarea     Set grid area
      Setgrid    setgridmask     Set grid mask
25 26
*/

Ralf Müller's avatar
Ralf Müller committed
27
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
28

29
#include "cdo_options.h"
30
#include "process_int.h"
31
#include "cdi_lockedIO.h"
32
#include "param_conversion.h"
33
#include <mpim_grid.h>
34
#include "griddes.h"
35
#include "gridreference.h"
36 37
#include "util_files.h"

38 39 40
static void
changeGrid(const int vlistID1, const int vlistID2, const int gridID2)
{
41
  const auto gridsize2 = gridInqSize(gridID2);
42
  int found = 0;
43
  const auto ngrids = vlistNgrids(vlistID1);
44 45
  for (int index = 0; index < ngrids; index++)
    {
46
      if (gridsize2 == gridInqSize(vlistGrid(vlistID1, index)))
47 48 49 50 51
        {
          vlistChangeGridIndex(vlistID2, index, gridID2);
          found++;
        }
    }
52
  if (!found) cdoWarning("No horizontal grid with %zu cells found!", gridsize2);
53 54 55
}

static void
56
setGridtype(const int vlistID1, const int vlistID2, const int gridtype, const std::string &gridname, const bool lregular,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
            const bool lregularnn, const bool ldereference, const bool lbounds, std::vector<int> &grid2_vgpm)
58 59
{
  int gridID2;
60
  auto lrgrid = false;
61
  const auto ngrids = vlistNgrids(vlistID1);
62 63
  for (int index = 0; index < ngrids; index++)
    {
64 65
      const auto gridID1 = vlistGrid(vlistID1, index);
      const auto gridtype1 = gridInqType(gridID1);
66 67 68 69 70 71 72 73 74 75
      gridID2 = -1;

      if (gridtype1 == GRID_GENERIC && gridInqSize(gridID1) == 1) continue;

      if (lregular || lregularnn)
        {
          if (gridtype1 == GRID_GAUSSIAN_REDUCED) gridID2 = gridToRegular(gridID1);
        }
      else if (ldereference)
        {
76 77 78
          const auto reference = dereferenceGrid(gridID1);
          if (reference.isValid) gridID2 = reference.gridID;
          if (reference.notFound) cdoAbort("Reference to horizontal grid not found!");
79 80 81 82 83 84 85 86 87 88
        }
      else
        {
          if (gridtype == GRID_CURVILINEAR)
            {
              gridID2 = (gridtype1 == GRID_CURVILINEAR) ? gridID1 : gridToCurvilinear(gridID1, lbounds);
            }
          else if (gridtype == GRID_UNSTRUCTURED)
            {
              gridID2 = gridToUnstructured(gridID1, 1);
89
              if (gridtype1 == GRID_GME)
90
                {
91
                  const auto grid2_nvgp = gridInqSize(gridID2);
92 93 94 95 96 97 98 99
                  grid2_vgpm.resize(grid2_nvgp);
                  gridInqMaskGME(gridID2, grid2_vgpm.data());
                  gridCompress(gridID2);
                }
            }
          else if (gridtype == GRID_LONLAT && gridtype1 == GRID_CURVILINEAR)
            {
              gridID2 = gridCurvilinearToRegular(gridID1);
100
              if (gridID2 == -1) cdoWarning("Conversion of curvilinear grid to regular grid failed!");
101 102 103 104
            }
          else if (gridtype == GRID_LONLAT && gridtype1 == GRID_UNSTRUCTURED)
            {
              gridID2 = -1;
105
              if (gridID2 == -1) cdoWarning("Conversion of unstructured grid to regular grid failed!");
106 107 108 109
            }
          else if (gridtype == GRID_LONLAT && gridtype1 == GRID_GENERIC)
            {
              gridID2 = -1;
110
              if (gridID2 == -1) cdoWarning("Conversion of generic grid to regular grid failed!");
111 112 113 114 115
            }
          else if (gridtype == GRID_LONLAT && gridtype1 == GRID_LONLAT)
            {
              gridID2 = gridID1;
            }
116 117 118 119 120 121 122
          else if (gridtype == GRID_PROJECTION)
            {
              if (gridInqType(gridID1) == GRID_PROJECTION)
                gridID2 = gridID1;
              else
                {
                  const int projID = gridInqProj(gridID1);
123
                  if (projID != CDI_UNDEFID && gridInqType(projID) == GRID_PROJECTION) gridID2 = projID;
124 125 126
                }
              if (gridID2 == -1) cdoWarning("Projection not found!");
            }
127
          else
128
            cdoAbort("Unsupported grid name: %s", gridname.c_str());
129 130 131 132 133 134 135 136 137 138 139 140 141 142
        }

      if (gridID2 == -1)
        {
          if (!(lregular || lregularnn)) cdoAbort("Unsupported grid type!");
        }

      if (gridID2 != -1)
        {
          if (lregular || lregularnn) lrgrid = true;
          vlistChangeGridIndex(vlistID2, index, gridID2);
        }
    }

143
  if ((lregular || lregularnn) && !lrgrid) cdoWarning("No reduced Gaussian grid found!");
144 145 146
}

static void
147
setGridcellArea(const int vlistID1, const int vlistID2, Varray<double> &gridcellArea)
148
{
149
  const auto areasize = gridcellArea.size();
150
  int numGridsChanges = 0;
151
  const auto ngrids = vlistNgrids(vlistID1);
152 153
  for (int index = 0; index < ngrids; index++)
    {
154 155
      const auto gridID1 = vlistGrid(vlistID1, index);
      const auto gridsize = gridInqSize(gridID1);
156 157
      if (gridsize == areasize)
        {
158
          const auto gridID2 = gridDuplicate(gridID1);
159 160 161 162 163
          gridDefArea(gridID2, gridcellArea.data());
          vlistChangeGridIndex(vlistID2, index, gridID2);
          numGridsChanges++;
        }
    }
164
  if (!numGridsChanges) cdoWarning("No grid with %zu cells found!", areasize);
165 166 167
}

static void
168
setGridMask(const int vlistID1, const int vlistID2, Varray<double> &gridmask)
169
{
170 171
  const auto masksize = gridmask.size();
  const auto ngrids = vlistNgrids(vlistID1);
172 173
  for (int index = 0; index < ngrids; index++)
    {
174 175
      const auto gridID1 = vlistGrid(vlistID1, index);
      const auto gridsize = gridInqSize(gridID1);
176 177 178 179 180
      if (gridsize == masksize)
        {
          std::vector<int> mask(masksize);
          for (size_t i = 0; i < masksize; i++)
            {
181
              mask[i] = (gridmask[i] < 0 || gridmask[i] > 255) ? 0 : (int) std::lround(gridmask[i]);
182
            }
183
          const auto gridID2 = gridDuplicate(gridID1);
184 185 186 187 188 189 190 191 192
          gridDefMask(gridID2, mask.data());
          vlistChangeGridIndex(vlistID2, index, gridID2);
        }
    }
}

static void
unsetGridMask(const int vlistID1, const int vlistID2)
{
193
  const auto ngrids = vlistNgrids(vlistID1);
194 195
  for (int index = 0; index < ngrids; index++)
    {
196 197
      const auto gridID1 = vlistGrid(vlistID1, index);
      const auto gridID2 = gridDuplicate(gridID1);
198 199 200 201 202
      gridDefMask(gridID2, nullptr);
      vlistChangeGridIndex(vlistID2, index, gridID2);
    }
}

203
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
setProjParams(const int vlistID1, const int vlistID2, const char *projparams)
205
{
206
  const auto ngrids = vlistNgrids(vlistID1);
207 208
  for (int index = 0; index < ngrids; index++)
    {
209
      const auto gridID1 = vlistGrid(vlistID1, index);
210 211
      if (gridInqType(gridID1) == GRID_PROJECTION)
        {
212
          const auto gridID2 = gridDuplicate(gridID1);
213
          cdiDefAttTxt(gridID2, CDI_GLOBAL, "proj4_params", (int) (strlen(projparams)), projparams);
214 215 216 217 218
          vlistChangeGridIndex(vlistID2, index, gridID2);
        }
    }
}

219
static void
220
readAreaFromFile(const std::string &areafile, Varray<double> &gridcellArea)
221
{
222
  auto searchName = false;
223 224 225
  std::string filename = areafile;
  std::string varname;

226
  if (!fileExists(areafile.c_str()))
227
    {
228
      const auto pos = filename.find_last_of(':');
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
      if (pos > 1 && pos < (filename.size() - 1))
230
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
          varname = filename.substr(pos + 1);
232 233 234 235 236
          filename = filename.substr(0, pos);
          searchName = true;
        }
    }

237 238
  const auto streamID = streamOpenReadLocked(filename.c_str());
  const auto vlistID = streamInqVlist(streamID);
239 240 241 242 243

  int svarID = 0;
  if (searchName)
    {
      char name[CDI_MAX_NAME];
244
      const auto nvars = vlistNvars(vlistID);
245 246 247 248
      int varID;
      for (varID = 0; varID < nvars; ++varID)
        {
          vlistInqVarName(vlistID, varID, name);
249
          if (cdo_cmpstr(name, varname.c_str()))
250 251 252 253 254 255 256 257
            {
              svarID = varID;
              break;
            }
        }
      if (varID == nvars) cdoAbort("Variable %s not found in %s\n", varname.c_str(), filename.c_str());
    }

258
  const auto nrecs = streamInqTimestep(streamID, 0);
259 260
  for (int recID = 0; recID < nrecs; ++recID)
    {
261
      int varID, levelID;
262 263 264
      streamInqRecord(streamID, &varID, &levelID);
      if (varID == svarID)
        {
265 266
          const auto gridID = vlistInqVarGrid(vlistID, varID);
          const auto areasize = gridInqSize(gridID);
267
          gridcellArea.resize(areasize);
268 269

          size_t nmiss;
270
          streamReadRecord(streamID, gridcellArea.data(), &nmiss);
271 272 273 274 275 276
          break;
        }
    }

  streamClose(streamID);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277

278 279
void *
Setgrid(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280 281
{
  int nrecs;
282 283
  int varID, levelID;
  int gridID2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
  int gridtype = -1;
285
  size_t nmiss;
286 287 288
  bool lregular = false;
  bool lregularnn = false;
  bool ldereference = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
  int number = 0, position = 0;
290
  bool lbounds = true;
291
  std::string gridname;
292 293
  const char *griduri = nullptr;
  const char *projparams = nullptr;
294
  std::vector<int> grid2_vgpm;
295
  Varray<double> gridmask, gridcellArea;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296

297
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298

Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
  // clang-format off
300 301 302 303 304 305 306 307 308
  const auto SETGRID       = cdoOperatorAdd("setgrid",       0, 0, "grid description file or name");
  const auto SETGRIDTYPE   = cdoOperatorAdd("setgridtype",   0, 0, "grid type");
  const auto SETGRIDAREA   = cdoOperatorAdd("setgridarea",   0, 0, "filename with area weights");
  const auto SETGRIDMASK   = cdoOperatorAdd("setgridmask",   0, 0, "filename with grid mask");
  const auto UNSETGRIDMASK = cdoOperatorAdd("unsetgridmask", 0, 0, nullptr);
  const auto SETGRIDNUMBER = cdoOperatorAdd("setgridnumber", 0, 0, "grid number and optionally grid position");
  const auto SETGRIDURI    = cdoOperatorAdd("setgriduri",    0, 0, "reference URI of the horizontal grid");
  const auto USEGRIDNUMBER = cdoOperatorAdd("usegridnumber", 0, 0, "use existing grid identified by grid number");
  const auto SETPROJPARAMS = cdoOperatorAdd("setprojparams", 0, 0, "proj library parameter (e.g.: +init=EPSG:3413)");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310

311
  const auto operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312

Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
  if (operatorID != UNSETGRIDMASK) operatorInputArg(cdoOperatorEnter(operatorID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314

315
  if (operatorID == SETGRID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317
      operatorCheckArgc(1);
318
      gridID2 = cdoDefineGrid(cdoOperatorArgv(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
    }
320
  else if (operatorID == SETGRIDTYPE)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322
      operatorCheckArgc(1);
323
      gridname = cdoOperatorArgv(0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324

325
      // clang-format off
326 327 328 329 330 331 332 333 334 335 336 337 338
      if      (gridname == "curvilinear0")  {gridtype = GRID_CURVILINEAR; lbounds = false;}
      else if (gridname == "curvilinear")   {gridtype = GRID_CURVILINEAR; lbounds = true;}
      else if (gridname == "cell")           gridtype = GRID_UNSTRUCTURED;
      else if (gridname == "unstructured0") {gridtype = GRID_UNSTRUCTURED; lbounds = false;}
      else if (gridname == "unstructured")  {gridtype = GRID_UNSTRUCTURED; lbounds = true;}
      else if (gridname == "generic")        gridtype = GRID_GENERIC;
      else if (gridname == "dereference")    ldereference = true;
      else if (gridname == "lonlat")         gridtype = GRID_LONLAT;
      else if (gridname == "gaussian")       gridtype = GRID_GAUSSIAN;
      else if (gridname == "regularnn")     {gridtype = GRID_GAUSSIAN; lregularnn = true;}
      else if (gridname == "regular")       {gridtype = GRID_GAUSSIAN; lregular = true;}
      else if (gridname == "projection")    {gridtype = GRID_PROJECTION;}
      else cdoAbort("Unsupported grid name: %s", gridname.c_str());
339
      // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
    }
341
  else if (operatorID == SETGRIDAREA)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
      operatorCheckArgc(1);
344

345
      readAreaFromFile(cdoOperatorArgv(0), gridcellArea);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
346

347
      if (Options::cdoVerbose)
348
        {
349
          const auto areasize = gridcellArea.size();
350 351
          auto mmm = varrayMinMaxMean(areasize, gridcellArea);
          cdoPrint("gridcellAreas: %zu %#12.5g%#12.5g%#12.5g", areasize, mmm.min, mmm.mean, mmm.max);
352
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
    }
354
  else if (operatorID == SETGRIDMASK)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
      operatorCheckArgc(1);
357
      const auto maskfile = cdoOperatorArgv(0).c_str();
358 359
      const auto streamID = streamOpenReadLocked(maskfile);
      const auto vlistID = streamInqVlist(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360

361
      (void) streamInqTimestep(streamID, 0);
362
      streamInqRecord(streamID, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363

364 365 366
      const auto missval = vlistInqVarMissval(vlistID, varID);
      const auto gridID = vlistInqVarGrid(vlistID, varID);
      const auto masksize = gridInqSize(gridID);
367
      gridmask.resize(masksize);
368

369
      streamReadRecord(streamID, gridmask.data(), &nmiss);
370
      streamClose(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371

372 373
      for (size_t i = 0; i < masksize; i++)
        if (DBL_IS_EQUAL(gridmask[i], missval)) gridmask[i] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
    }
375
  else if (operatorID == USEGRIDNUMBER)
376 377
    {
      operatorCheckArgc(1);
378
      number = parameter2int(cdoOperatorArgv(0));
379
    }
380
  else if (operatorID == SETGRIDNUMBER)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
381
    {
382 383
      if (operatorArgc() >= 1 && operatorArgc() <= 2)
        {
384 385
          number = parameter2int(cdoOperatorArgv(0));
          if (operatorArgc() == 2) position = parameter2int(cdoOperatorArgv(1));
386
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
      else
388 389 390
        {
          operatorCheckArgc(1);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
    }
392
  else if (operatorID == SETGRIDURI)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393 394
    {
      operatorCheckArgc(1);
395
      griduri = cdoOperatorArgv(0).c_str();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
    }
397 398 399
  else if (operatorID == SETPROJPARAMS)
    {
      operatorCheckArgc(1);
400
      projparams = cdoOperatorArgv(0).c_str();
401
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402

403
  const auto streamID1 = cdoOpenRead(0);
404

405 406
  const auto vlistID1 = cdoStreamInqVlist(streamID1);
  const auto vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407

408 409
  const auto taxisID1 = vlistInqTaxis(vlistID1);
  const auto taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410 411
  vlistDefTaxis(vlistID2, taxisID2);

412
  if (operatorID == SETGRID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
    {
414
      changeGrid(vlistID1, vlistID2, gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
  else if (operatorID == SETGRIDNUMBER || operatorID == SETGRIDURI || operatorID == USEGRIDNUMBER)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
      int gridID2x = -1;
419 420
      if (operatorID == SETGRIDNUMBER)
        {
421
          const auto gridID1 = vlistGrid(vlistID1, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
422 423 424
          gridID2x = gridCreate(GRID_UNSTRUCTURED, gridInqSize(gridID1));
          cdiDefKeyInt(gridID2x, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, number);
          cdiDefKeyInt(gridID2x, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, position);
425 426
        }
      else if (operatorID == USEGRIDNUMBER)
427
        {
428
          if (number < 1 || number > vlistNgrids(vlistID1))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
            cdoAbort("Invalid grid number: %d (max = %d)!", number, vlistNgrids(vlistID1));
430

Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
          gridID2x = vlistGrid(vlistID1, number - 1);
432
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
      else
434
        {
435
          const auto gridID1 = vlistGrid(vlistID1, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436 437
          gridID2x = gridDuplicate(gridID1);
          cdiDefKeyString(gridID2x, CDI_GLOBAL, CDI_KEY_REFERENCEURI, griduri);
438
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439

Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
      changeGrid(vlistID1, vlistID2, gridID2x);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
    }
442
  else if (operatorID == SETGRIDTYPE)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443
    {
444
      setGridtype(vlistID1, vlistID2, gridtype, gridname, lregular, lregularnn, ldereference, lbounds, grid2_vgpm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
445
    }
446
  else if (operatorID == SETGRIDAREA)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
    {
448
      setGridcellArea(vlistID1, vlistID2, gridcellArea);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449
    }
450
  else if (operatorID == SETGRIDMASK)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
451
    {
452
      setGridMask(vlistID1, vlistID2, gridmask);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
    }
454
  else if (operatorID == UNSETGRIDMASK)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
    {
456
      unsetGridMask(vlistID1, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
    }
458 459 460 461
  else if (operatorID == SETPROJPARAMS)
    {
      setProjParams(vlistID1, vlistID2, projparams);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462

463
  const auto streamID2 = cdoOpenWrite(1);
464

465
  cdoDefVlist(streamID2, vlistID2);
466
  // vlistPrint(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467

468
  auto gridsizemax = (lregular || lregularnn) ? vlistGridsizeMax(vlistID2) : vlistGridsizeMax(vlistID1);
469

470
  if (vlistNumber(vlistID1) != CDI_REAL) gridsizemax *= 2;
471
  Varray<double> array(gridsizemax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472

473
  int tsID = 0;
474
  while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475 476
    {
      taxisCopyTimestep(taxisID2, taxisID1);
477
      cdoDefTimestep(streamID2, tsID);
478 479 480

      for (int recID = 0; recID < nrecs; recID++)
        {
481 482
          cdoInqRecord(streamID1, &varID, &levelID);
          cdoDefRecord(streamID2, varID, levelID);
483

484
          cdoReadRecord(streamID1, array.data(), &nmiss);
485

486
          const auto gridID1 = vlistInqVarGrid(vlistID1, varID);
487 488 489 490
          if (lregular || lregularnn)
            {
              if (gridInqType(gridID1) == GRID_GAUSSIAN_REDUCED)
                {
491
                  const auto missval = vlistInqVarMissval(vlistID1, varID);
492
                  const int lnearst = lregularnn ? 1 : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
                  field2regular(gridID1, vlistInqVarGrid(vlistID2, varID), missval, array.data(), nmiss, lnearst);
494 495 496 497
                }
            }
          else if (gridInqType(gridID1) == GRID_GME)
            {
498
              const auto gridsize = gridInqSize(gridID1);
499 500 501 502 503
              size_t j = 0;
              for (size_t i = 0; i < gridsize; i++)
                if (grid2_vgpm[i]) array[j++] = array[i];
            }

504
          cdoWriteRecord(streamID2, array.data(), nmiss);
505
        }
506

Uwe Schulzweida's avatar
Uwe Schulzweida committed
507 508 509
      tsID++;
    }

510 511
  cdoStreamClose(streamID2);
  cdoStreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
512 513 514

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
515
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
516
}