Collgrid.cc 18.3 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.
*/

Ralf Mueller's avatar
Ralf Mueller committed
18
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
19

20
#include "dmemory.h"
21
#include "functs.h"
22
#include "process_int.h"
23
#include "cdo_vlist.h"
24
#include "param_conversion.h"
25
#include <mpim_grid.h>
26
#include "util_files.h"
Oliver Heidmann's avatar
Oliver Heidmann committed
27
#include "cdo_options.h"
28
#include "cdi_lockedIO.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29

30
31
static int globalGridType = CDI_UNDEFID;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
struct EnsfileType
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
{
34
  CdoStreamID streamID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
  int vlistID;
36
  size_t nmiss;
37
  size_t gridsize;
38
  std::vector<long> gridindex;
39
  Varray<double> array;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41

Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
struct xyinfoType
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
{
  double x, y;
  int id;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47

48
static int
49
cmpx(const void *a, const void *b)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
52
  const auto x = ((const xyinfoType *) a)->x;
  const auto y = ((const xyinfoType *) b)->x;
53
  return ((x > y) - (x < y)) * 2 + (x > y) - (x < y);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
}

56
static int
57
cmpxy_lt(const void *a, const void *b)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
{
59
60
61
62
  const auto x1 = ((const xyinfoType *) a)->x;
  const auto x2 = ((const xyinfoType *) b)->x;
  const auto y1 = ((const xyinfoType *) a)->y;
  const auto y2 = ((const xyinfoType *) b)->y;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63

64
  int cmp = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
  if (y1 < y2 || (std::fabs(y1 - y2) <= 0 && x1 < x2))
66
    cmp = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
  else if (y1 > y2 || (std::fabs(y1 - y2) <= 0 && x1 > x2))
68
    cmp = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69

70
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
72
}

73
static int
74
cmpxy_gt(const void *a, const void *b)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
{
76
77
78
79
  const auto x1 = ((const xyinfoType *) a)->x;
  const auto x2 = ((const xyinfoType *) b)->x;
  const auto y1 = ((const xyinfoType *) a)->y;
  const auto y2 = ((const xyinfoType *) b)->y;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80

81
  int cmp = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
  if (y1 > y2 || (std::fabs(y1 - y2) <= 0 && x1 < x2))
83
    cmp = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
  else if (y1 < y2 || (std::fabs(y1 - y2) <= 0 && x1 > x2))
85
    cmp = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86

87
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
89
}

90
static int
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
genGrid(int ngrids, int nfiles, std::vector<EnsfileType> &ef, bool ginit, int igrid, long nxblocks)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
{
93
94
95
  bool lsouthnorth = true;
  bool lregular = false;
  bool lcurvilinear = false;
96
  bool lunstructured = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
  int gridID2 = -1;

99
  long nx = -1;
100
  if (nxblocks != -1) nx = nxblocks;
101

102
  auto gridID = vlistGrid(ef[0].vlistID, igrid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
104
105
  auto gridtype0 = (globalGridType != CDI_UNDEFID) ? globalGridType : gridInqType(gridID);
  if (ngrids > 1 && gridtype0 == GRID_GENERIC && gridInqXsize(gridID) == 0 && gridInqYsize(gridID) == 0) return gridID2;
  if (gridtype0 == GRID_UNSTRUCTURED) lunstructured = true;
106
107
  const int nv = lunstructured ? gridInqNvertex(gridID) : 0;
  const bool lcenter = globalGridType == CDI_UNDEFID && gridInqXvals(gridID, nullptr) > 0 && gridInqYvals(gridID, nullptr) > 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
  const bool lbounds = lunstructured && globalGridType == CDI_UNDEFID && gridInqXbounds(gridID, nullptr) > 0
                       && gridInqYbounds(gridID, nullptr) > 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110

Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
  std::vector<long> xsize(nfiles), ysize(nfiles);
112
  xyinfoType *xyinfo = (xyinfoType *) Malloc(nfiles * sizeof(xyinfoType));
113
114
  Varray2D<double> xvals(nfiles), yvals(nfiles);
  Varray2D<double> xbounds(nfiles), ybounds(nfiles);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115

116
  for (int fileID = 0; fileID < nfiles; fileID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
    {
118
      gridID = vlistGrid(ef[fileID].vlistID, igrid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
      const auto gridtype = (globalGridType != CDI_UNDEFID) ? globalGridType : gridInqType(gridID);
120
      if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN)
121
        lregular = true;
122
      else if (gridtype == GRID_CURVILINEAR)
123
        lcurvilinear = true;
124
125
      else if (gridtype == GRID_UNSTRUCTURED)
        lunstructured = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
      else if (gridtype == GRID_GENERIC /*&& gridInqXsize(gridID) > 0 && gridInqYsize(gridID) > 0*/)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
        ;
128
      else
129
        cdoAbort("Unsupported grid type: %s!", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130

131
      xsize[fileID] = lunstructured ? gridInqSize(gridID) : gridInqXsize(gridID);
132
      ysize[fileID] = lunstructured ? 1 : gridInqYsize(gridID);
133
134
      if (xsize[fileID] == 0) xsize[fileID] = 1;
      if (ysize[fileID] == 0) ysize[fileID] = 1;
135

136
      if (lregular)
137
        {
138
139
          xvals[fileID].resize(xsize[fileID]);
          yvals[fileID].resize(ysize[fileID]);
140
        }
141
      else if (lcurvilinear || lunstructured)
142
        {
143
144
145
146
147
          if (lcenter)
            {
              xvals[fileID].resize(xsize[fileID] * ysize[fileID]);
              yvals[fileID].resize(xsize[fileID] * ysize[fileID]);
            }
148
149
150
151
152
          if (lbounds)
            {
              xbounds[fileID].resize(nv * xsize[fileID] * ysize[fileID]);
              ybounds[fileID].resize(nv * xsize[fileID] * ysize[fileID]);
            }
153
        }
154

155
      if (lregular || lcurvilinear || lunstructured)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
        {
157
158
159
160
161
          if (lcenter)
            {
              gridInqXvals(gridID, xvals[fileID].data());
              gridInqYvals(gridID, yvals[fileID].data());
            }
162
163
164
165
166
          if (lbounds)
            {
              gridInqXbounds(gridID, xbounds[fileID].data());
              gridInqYbounds(gridID, ybounds[fileID].data());
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
      // printf("fileID %d, gridID %d\n", fileID, gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169

170
      if (lregular)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
        {
172
173
          xyinfo[fileID].x = xvals[fileID][0];
          xyinfo[fileID].y = yvals[fileID][0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
          xyinfo[fileID].id = fileID;

176
          if (ysize[fileID] > 1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
              if (yvals[fileID][0] > yvals[fileID][ysize[fileID] - 1]) lsouthnorth = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
180
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
      else
        {
183
184
          xyinfo[fileID].x = 0;
          xyinfo[fileID].y = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
          xyinfo[fileID].id = fileID;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
    }
188

189
  if (Options::cdoVerbose && lregular)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
    for (int fileID = 0; fileID < nfiles; fileID++) printf("1 %d %g %g \n", xyinfo[fileID].id, xyinfo[fileID].x, xyinfo[fileID].y);
191

192
  if (lregular)
193
    {
194
      std::qsort(xyinfo, nfiles, sizeof(xyinfoType), cmpx);
195

196
      if (Options::cdoVerbose)
197
        for (int fileID = 0; fileID < nfiles; fileID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
          printf("2 %d %g %g \n", xyinfo[fileID].id, xyinfo[fileID].x, xyinfo[fileID].y);
199

200
      if (lsouthnorth)
201
        std::qsort(xyinfo, nfiles, sizeof(xyinfoType), cmpxy_lt);
202
      else
203
        std::qsort(xyinfo, nfiles, sizeof(xyinfoType), cmpxy_gt);
204

205
      if (Options::cdoVerbose)
206
        for (int fileID = 0; fileID < nfiles; fileID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
          printf("3 %d %g %g \n", xyinfo[fileID].id, xyinfo[fileID].x, xyinfo[fileID].y);
208

209
      if (nx <= 0)
210
        {
211
          nx = 1;
212
          for (int fileID = 1; fileID < nfiles; fileID++)
213
            {
214
215
216
217
              if (DBL_IS_EQUAL(xyinfo[0].y, xyinfo[fileID].y))
                nx++;
              else
                break;
218
            }
219
220
221
        }
    }
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
    {
223
      if (nx <= 0) nx = nfiles;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
    }
225

226
  const long ny = nfiles / nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
  if (nx * ny != nfiles) cdoAbort("Number of input files (%ld) and number of blocks (%ldx%ld) differ!", nfiles, nx, ny);
228

229
  long xsize2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
  for (long i = 0; i < nx; ++i) xsize2 += xsize[xyinfo[i].id];
231
  long ysize2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
  for (long j = 0; j < ny; ++j) ysize2 += ysize[xyinfo[j * nx].id];
233
  if (Options::cdoVerbose) cdoPrint("xsize2 %ld  ysize2 %ld", xsize2, ysize2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234

Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
  {  // verify size of data
236
    const long xs = xsize[xyinfo[0].id];
237
    for (long j = 1; j < ny; ++j)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
      if (xsize[xyinfo[j * nx].id] != xs) cdoAbort("xsize=%ld differ from first file (xsize=%ld)!", xsize[xyinfo[j * nx].id], xs);
239
    const long ys = ysize[xyinfo[0].id];
240
    for (long i = 1; i < nx; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
      if (ysize[xyinfo[i].id] != ys) cdoAbort("ysize=%ld differ from first file (ysize=%ld)!", ysize[xyinfo[i].id], ys);
242
243
  }

244
245
  Varray<double> xvals2, yvals2;
  Varray<double> xbounds2, ybounds2;
246
  if (lregular)
247
    {
248
249
      xvals2.resize(xsize2);
      yvals2.resize(ysize2);
250
    }
251
  else if (lcurvilinear || lunstructured)
252
    {
253
254
255
256
257
      if (lcenter)
        {
          xvals2.resize(xsize2 * ysize2);
          yvals2.resize(xsize2 * ysize2);
        }
258
259
260
261
262
      if (lbounds)
        {
          xbounds2.resize(nv * xsize2 * ysize2);
          ybounds2.resize(nv * xsize2 * ysize2);
        }
263
    }
264

Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
  std::vector<long> xoff(nx + 1), yoff(ny + 1);
266
267

  xoff[0] = 0;
268
  for (long i = 0; i < nx; ++i)
269
    {
270
      const long idx = xyinfo[i].id;
271
      if (lregular) arrayCopy(xsize[idx], xvals[idx].data(), &xvals2[xoff[i]]);
272
      xoff[i + 1] = xoff[i] + xsize[idx];
273
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

275
  yoff[0] = 0;
276
  for (long j = 0; j < ny; ++j)
277
    {
278
      const long idx = xyinfo[j * nx].id;
279
      if (lregular) arrayCopy(ysize[idx], yvals[idx].data(), &yvals2[yoff[j]]);
280
      yoff[j + 1] = yoff[j] + ysize[idx];
281
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282

Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
  if (!ginit)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
    {
285
286
      for (int fileID = 0; fileID < nfiles; fileID++)
        {
287
288
289
          const long idx = xyinfo[fileID].id;
          const long iy = fileID / nx;
          const long ix = fileID - iy * nx;
290

291
          const long offset = yoff[iy] * xsize2 + xoff[ix];
292
293
          // printf("fileID %d %d, iy %d, ix %d, offset %d\n", fileID, xyinfo[fileID].id, iy, ix, offset);

294
295
296
297
          long ij = 0;
          for (long j = 0; j < ysize[idx]; ++j)
            for (long i = 0; i < xsize[idx]; ++i)
              {
298
                if (lcurvilinear || lunstructured)
299
                  {
300
301
302
303
304
                    if (lcenter)
                      {
                        xvals2[offset + j * xsize2 + i] = xvals[idx][ij];
                        yvals2[offset + j * xsize2 + i] = yvals[idx][ij];
                      }
305
306
307
308
                    if (lbounds)
                      {
                        for (long k = 0; k < nv; ++k)
                          {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
                            xbounds2[(offset + j * xsize2 + i) * nv + k] = xbounds[idx][ij * nv + k];
                            ybounds2[(offset + j * xsize2 + i) * nv + k] = ybounds[idx][ij * nv + k];
311
312
                          }
                      }
313
                  }
314
315
316
                ef[idx].gridindex[ij++] = offset + j * xsize2 + i;
              }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317
318
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
  gridID2 = gridCreate(gridtype0, xsize2 * ysize2);
320
321
322
323
324
  if (!lunstructured)
    {
      gridDefXsize(gridID2, xsize2);
      gridDefYsize(gridID2, ysize2);
    }
325
326
327
328
  else if (nv > 0)
    {
      gridDefNvertex(gridID2, nv);
    }
329

330
  if (lregular || lcurvilinear || lunstructured)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
    {
332
333
334
335
336
      if (lcenter)
        {
          gridDefXvals(gridID2, xvals2.data());
          gridDefYvals(gridID2, yvals2.data());
        }
337
338
339
340
341
      if (lbounds)
        {
          gridDefXbounds(gridID2, xbounds2.data());
          gridDefYbounds(gridID2, ybounds2.data());
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343

Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
345
  Free(xyinfo);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
346
  gridID = vlistGrid(ef[0].vlistID, igrid);
347
348

  grid_copy_attributes(gridID, gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349

Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
  if (gridtype0 == GRID_PROJECTION) grid_copy_mapping(gridID, gridID2);
351

352
  return gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
354
}

355
356
void *
Collgrid(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
357
{
358
  int nxblocks = -1;
359
  int varID, levelID;
360
  int nrecs0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361

362
  cdoInitialize(process);
363

364
365
  const auto nfiles = cdoStreamCnt() - 1;
  const auto ofilename = cdoGetStreamName(nfiles);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366

367
  if (!Options::cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename))
368
    cdoAbort("Outputfile %s already exists!", ofilename);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369

Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
  std::vector<EnsfileType> ef(nfiles);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371

372
  for (int fileID = 0; fileID < nfiles; fileID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
    {
374
      ef[fileID].streamID = cdoOpenRead(fileID);
375
      ef[fileID].vlistID = cdoStreamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
376
377
    }

378
  const auto vlistID1 = ef[0].vlistID;
379
  vlistClearFlag(vlistID1);
380
381

  /* check that the contents is always the same */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
382
  for (int fileID = 1; fileID < nfiles; fileID++) vlistCompare(vlistID1, ef[fileID].vlistID, CMP_NAME | CMP_NLEVEL);
383

384
  const auto nvars = vlistNvars(vlistID1);
385
  std::vector<bool> vars(nvars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
  for (varID = 0; varID < nvars; varID++) vars[varID] = false;
387
  std::vector<bool> vars1(nvars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
388
  for (varID = 0; varID < nvars; varID++) vars1[varID] = false;
389

390
  auto nsel = operatorArgc();
391
  int noff = 0;
392

393
  if (nsel > 0)
394
    {
395
396
      const char *argument = cdoOperatorArgv(0).c_str();
      if (strcmp(argument, "gridtype=unstructured") == 0)
397
398
        {
          nsel--;
399
          noff++;
400
401
402
403
          globalGridType = GRID_UNSTRUCTURED;
        }
      else
        {
404
405
          int len = (int) strlen(argument);
          while (--len >= 0 && isdigit(argument[len]))
406
407
408
409
410
411
            ;

          if (len == -1)
            {
              nsel--;
              noff++;
412
              nxblocks = parameter2int(argument);
413
            }
414
415
416
        }
    }

417
  if (nsel == 0)
418
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
      for (varID = 0; varID < nvars; varID++) vars1[varID] = true;
420
421
422
    }
  else
    {
423
      if (Options::cdoVerbose)
424
        for (int i = 0; i < nsel; i++) cdoPrint("name %d = %s\n", i + 1, cdoOperatorArgv(noff + i).c_str());
425

426
      std::vector<bool> selfound(nsel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
      for (int i = 0; i < nsel; i++) selfound[i] = false;
428
429

      char varname[CDI_MAX_NAME];
430
431
432
433
434
435
      for (varID = 0; varID < nvars; varID++)
        {
          vlistInqVarName(vlistID1, varID, varname);

          for (int isel = 0; isel < nsel; isel++)
            {
436
              if (cdoOperatorArgv(noff + isel).compare(varname) == 0)
437
438
439
440
441
442
443
                {
                  selfound[isel] = true;
                  vars1[varID] = true;
                }
            }
        }

444
      int err = 0;
445
      for (int isel = 0; isel < nsel; isel++)
446
447
448
449
        {
          if (selfound[isel] == false)
            {
              err++;
450
              cdoWarning("Variable name %s not found!", cdoOperatorArgv(noff + isel).c_str());
451
452
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
      if (err) cdoAbort("Could not find all requested variables: (%d/%d)", nsel - err, nsel);
454
455
    }

456
  for (varID = 0; varID < nvars; varID++)
457
    {
458
459
      if (vars1[varID])
        {
460
461
          const auto zaxisID = vlistInqVarZaxis(vlistID1, varID);
          const auto nlevs = zaxisInqSize(zaxisID);
462
          for (int levID = 0; levID < nlevs; levID++) vlistDefFlag(vlistID1, varID, levID, true);
463
        }
464
465
    }

466
  for (int fileID = 0; fileID < nfiles; fileID++)
467
    {
468
      const auto gridsize = vlistGridsizeMax(ef[fileID].vlistID);
469
      ef[fileID].gridsize = gridsize;
470
471
      ef[fileID].gridindex.resize(gridsize);
      ef[fileID].array.resize(gridsize);
472
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473

474
  const auto vlistID2 = vlistCreate();
475
  cdoVlistCopyFlag(vlistID2, vlistID1);
476
  /*
477
  if ( Options::cdoVerbose )
478
479
    {
      vlistPrint(vlistID1);
480
      vlistPrint(vlistID2);
481
    }
482
  */
483
484
  // auto vlistID2 = vlistDuplicate(vlistID1);
  const auto nvars2 = vlistNvars(vlistID2);
485

486
487
  const auto ngrids1 = vlistNgrids(vlistID1);
  const auto ngrids2 = vlistNgrids(vlistID2);
488

489
  std::vector<int> gridIDs(ngrids2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490

491
  bool ginit = false;
492
  for (int i2 = 0; i2 < ngrids2; ++i2)
493
    {
494
      int i1;
495
496
      for (i1 = 0; i1 < ngrids1; ++i1)
        if (vlistGrid(vlistID1, i1) == vlistGrid(vlistID2, i2)) break;
497
498
499

      //   printf("i1 %d i2 %d\n", i1, i2);

500
501
502
503
504
      if (!ginit)
        {
          gridIDs[i2] = genGrid(ngrids2, nfiles, ef, ginit, i1, nxblocks);
          if (gridIDs[i2] != -1) ginit = true;
        }
505
      else
506
        gridIDs[i2] = genGrid(ngrids2, nfiles, ef, ginit, i1, nxblocks);
507
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508

509
510
  const auto taxisID1 = vlistInqTaxis(vlistID1);
  const auto taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
  vlistDefTaxis(vlistID2, taxisID2);

513
  size_t gridsize2 = 0;
514
  for (int i = 0; i < ngrids2; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
515
    {
516
517
518
      if (gridIDs[i] != -1)
        {
          if (gridsize2 == 0) gridsize2 = gridInqSize(gridIDs[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
520
          if (gridsize2 != gridInqSize(gridIDs[i]))
            cdoAbort("Target gridsize differ (gridID=%d  gridsize=%zu)!", i + 1, gridInqSize(gridIDs[i]));
521
522
          vlistChangeGridIndex(vlistID2, i, gridIDs[i]);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
    }
524

525
  for (varID = 0; varID < nvars2; varID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526
    {
527
      const auto gridID = vlistInqVarGrid(vlistID2, varID);
528

529
530
531
532
      for (int i = 0; i < ngrids2; ++i)
        {
          if (gridIDs[i] != -1)
            {
533
              if (gridID == vlistGrid(vlistID2, i) || gridsize2 == gridInqSize(gridID)) vars[varID] = true;
534
535
536
              break;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
538
    }

539
  const auto streamID2 = cdoOpenWrite(nfiles);
540
  cdoDefVlist(streamID2, vlistID2);
541

542
  Varray<double> array2;
543
  if (gridsize2 > 0) array2.resize(gridsize2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
544

Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547
  do
    {
548
      nrecs0 = cdoStreamInqTimestep(ef[0].streamID, tsID);
549
550
      for (int fileID = 1; fileID < nfiles; fileID++)
        {
551
          const auto nrecs = cdoStreamInqTimestep(ef[fileID].streamID, tsID);
552
          if (nrecs != nrecs0)
553
554
            cdoAbort("Number of records at time step %d of %s and %s differ!", tsID + 1, cdoGetStreamName(0),
                     cdoGetStreamName(fileID));
555
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
557
558

      taxisCopyTimestep(taxisID2, taxisID1);

559
      if (nrecs0 > 0) cdoDefTimestep(streamID2, tsID);
560
561
562

      for (int recID = 0; recID < nrecs0; recID++)
        {
563
          cdoInqRecord(ef[0].streamID, &varID, &levelID);
564
          if (Options::cdoVerbose && tsID == 0) printf(" tsID, recID, varID, levelID %d %d %d %d\n", tsID, recID, varID, levelID);
565
566
567
568

          for (int fileID = 1; fileID < nfiles; fileID++)
            {
              int varIDx, levelIDx;
569
              cdoInqRecord(ef[fileID].streamID, &varIDx, &levelIDx);
570
571
            }

572
          if (vlistInqFlag(vlistID1, varID, levelID) == true)
573
            {
574
575
              const auto varID2 = vlistFindVar(vlistID2, varID);
              const auto levelID2 = vlistFindLevel(vlistID2, varID, levelID);
576
              if (Options::cdoVerbose && tsID == 0) printf("varID %d %d levelID %d %d\n", varID, varID2, levelID, levelID2);
577

578
              const auto missval = vlistInqVarMissval(vlistID2, varID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
579
              for (size_t i = 0; i < gridsize2; i++) array2[i] = missval;
580
581

#ifdef _OPENMP
582
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
583
#endif
584
585
              for (int fileID = 0; fileID < nfiles; fileID++)
                {
586
                  cdoReadRecord(ef[fileID].streamID, &ef[fileID].array[0], &ef[fileID].nmiss);
587

588
589
                  if (vars[varID2])
                    {
590
                      const auto gridsize = ef[fileID].gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
                      for (size_t i = 0; i < gridsize; ++i) array2[ef[fileID].gridindex[i]] = ef[fileID].array[i];
592
593
594
                    }
                }

595
              cdoDefRecord(streamID2, varID2, levelID2);
596
597
598

              if (vars[varID2])
                {
599
600
                  const auto nmiss = varrayNumMV(gridsize2, array2, missval);
                  cdoWriteRecord(streamID2, array2.data(), nmiss);
601
602
                }
              else
603
                cdoWriteRecord(streamID2, &ef[0].array[0], ef[0].nmiss);
604
605
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
606
607
608

      tsID++;
    }
609
  while (nrecs0 > 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
610

611
  for (int fileID = 0; fileID < nfiles; fileID++) cdoStreamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612

613
  cdoStreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614
615
616

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
}