Collgrid.c 15.5 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-2017 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

Ralf Mueller's avatar
Ralf Mueller committed
19
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
21
22
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
23
#include "grid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
25
26
27
28
29
30
31
#include "util.h"


typedef struct
{
  int streamID;
  int vlistID;
  int gridID;
32
  int nmiss;
33
34
  int gridsize;
  int *gridindex;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
  double *array;
} ens_file_t;


typedef struct
{
  double x, y;
  int id;
}
xyinfo_t;

static
int cmpx(const void *s1, const void *s2)
{
  int cmp = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
  const xyinfo_t *xy1 = (const xyinfo_t *) s1;
  const xyinfo_t *xy2 = (const xyinfo_t *) s2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
53
54
55

  if      ( xy1->x < xy2->x ) cmp = -1;
  else if ( xy1->x > xy2->x ) cmp =  1;

56
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
58
59
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
int cmpxy_lt(const void *s1, const void *s2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
62
{
  int cmp = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
  const xyinfo_t *xy1 = (const xyinfo_t *) s1;
  const xyinfo_t *xy2 = (const xyinfo_t *) s2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65

Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
67
  if      ( xy1->y < xy2->y || (!(fabs(xy1->y - xy2->y) > 0) && xy1->x < xy2->x) ) cmp = -1;
  else if ( xy1->y > xy2->y || (!(fabs(xy1->y - xy2->y) > 0) && xy1->x > xy2->x) ) cmp =  1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
73
74
75
static
int cmpxy_gt(const void *s1, const void *s2)
{
  int cmp = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
77
  const xyinfo_t *xy1 = (const xyinfo_t *) s1;
  const xyinfo_t *xy2 = (const xyinfo_t *) s2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
80
81

  if      ( xy1->y > xy2->y || (!(fabs(xy1->y - xy2->y) > 0) && xy1->x < xy2->x) ) cmp = -1;
  else if ( xy1->y < xy2->y || (!(fabs(xy1->y - xy2->y) > 0) && xy1->x > xy2->x) ) cmp =  1;

82
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
static
86
int genGrid(int ngrids, int nfiles, ens_file_t *ef, bool ginit, int igrid, int nxblocks)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
{
88
89
90
  bool lsouthnorth = true;
  bool lregular = false;
  bool lcurvilinear = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
  int gridID2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
  double *xvals2 = NULL, *yvals2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93

94
95
96
  int nx = -1;
  if ( nxblocks != -1 ) nx = nxblocks;

97
98
  int gridID   = vlistGrid(ef[0].vlistID, igrid);
  int gridtype = gridInqType(gridID);
99
  if ( ngrids > 1 && gridtype == GRID_GENERIC && gridInqXsize(gridID) == 0 && gridInqYsize(gridID) == 0 )
100
    return gridID2;
101

102
103
104
105
106
  int *xsize = (int*) Malloc(nfiles*sizeof(int));
  int *ysize = (int*) Malloc(nfiles*sizeof(int));
  xyinfo_t *xyinfo = (xyinfo_t*) Malloc(nfiles*sizeof(xyinfo_t));
  double **xvals = (double**) Malloc(nfiles*sizeof(double*));
  double **yvals = (double**) Malloc(nfiles*sizeof(double*));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107

108
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
      gridID   = vlistGrid(ef[fileID].vlistID, igrid);
111
      int gridtype = gridInqType(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
      if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN )
113
        lregular = true;
114
      else if ( gridtype == GRID_CURVILINEAR )
115
        lcurvilinear = true;
116
      else if ( gridtype == GRID_GENERIC /*&& gridInqXsize(gridID) > 0 && gridInqYsize(gridID) > 0*/ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
        ;
118
      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
	cdoAbort("Unsupported grid type: %s!", gridNamePtr(gridtype));

121
122
      xsize[fileID] = gridInqXsize(gridID);
      ysize[fileID] = gridInqYsize(gridID);
123
124
      if ( xsize[fileID] == 0 ) xsize[fileID] = 1;
      if ( ysize[fileID] == 0 ) ysize[fileID] = 1;
125

126
127
      if ( lregular )
        {
128
129
          xvals[fileID] = (double*) Malloc(xsize[fileID]*sizeof(double));
          yvals[fileID] = (double*) Malloc(ysize[fileID]*sizeof(double));
130
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
      else if ( lcurvilinear )
132
        {
133
134
          xvals[fileID] = (double*) Malloc(xsize[fileID]*ysize[fileID]*sizeof(double));
          yvals[fileID] = (double*) Malloc(xsize[fileID]*ysize[fileID]*sizeof(double));
135
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
137
138
139
140
      else
        {
          xvals[fileID] = NULL;
          yvals[fileID] = NULL;
        }
141
        
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
144
145
146
      if ( lregular || lcurvilinear )
        {
          gridInqXvals(gridID, xvals[fileID]);
          gridInqYvals(gridID, yvals[fileID]);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
      // printf("fileID %d, gridID %d\n", fileID, gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148

Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
150
151
152
153
154
155
156
      if ( lregular )
        {
          xyinfo[fileID].x  = xvals[fileID][0];
          xyinfo[fileID].y  = yvals[fileID][0];
          xyinfo[fileID].id = fileID;

          if ( ysize[fileID] > 1 )
            {
157
              if ( yvals[fileID][0] > yvals[fileID][ysize[fileID]-1] ) lsouthnorth = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
159
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
160
161
162
163
164
165
      else
        {
          xyinfo[fileID].x  = 0;
          xyinfo[fileID].y  = 0;
          xyinfo[fileID].id = fileID;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
    }
167

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

172
173
174
  if ( lregular )
    {
      qsort(xyinfo, nfiles, sizeof(xyinfo_t), cmpx);  	      
175

176
      if ( cdoVerbose )
177
        for ( int fileID = 0; fileID < nfiles; fileID++ )
178
179
180
181
182
183
184
185
          printf("2 %d %g %g \n",  xyinfo[fileID].id, xyinfo[fileID].x, xyinfo[fileID].y);

      if ( lsouthnorth )
        qsort(xyinfo, nfiles, sizeof(xyinfo_t), cmpxy_lt);  
      else
        qsort(xyinfo, nfiles, sizeof(xyinfo_t), cmpxy_gt);  	      

      if ( cdoVerbose )
186
        for ( int fileID = 0; fileID < nfiles; fileID++ )
187
188
          printf("3 %d %g %g \n",  xyinfo[fileID].id, xyinfo[fileID].x, xyinfo[fileID].y);

189
      if ( nx <= 0 )
190
        {
191
192
193
194
195
196
          nx = 1;
          for ( int fileID = 1; fileID < nfiles; fileID++ )
            {
              if ( DBL_IS_EQUAL(xyinfo[0].y, xyinfo[fileID].y) ) nx++;
              else break;
            }
197
198
199
        }
    }
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
    {
201
      if ( nx <= 0 ) nx = nfiles;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
    }
203

204
  int ny = nfiles/nx;
205
  if ( nx*ny != nfiles ) cdoAbort("Number of input files (%d) and number of blocks (%dx%d) differ!", nfiles, nx, ny);
206
 
207
  int xsize2 = 0;
208
  for ( int i = 0; i < nx; ++i ) xsize2 += xsize[xyinfo[i].id];
209
  int ysize2 = 0;
210
  for ( int j = 0; j < ny; ++j ) ysize2 += ysize[xyinfo[j*nx].id];
211
  if ( cdoVerbose ) cdoPrint("xsize2 %d  ysize2 %d", xsize2, ysize2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212

213
214
  if ( lregular )
    {
215
216
      xvals2 = (double*) Malloc(xsize2*sizeof(double));
      yvals2 = (double*) Malloc(ysize2*sizeof(double));
217
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
  else if ( lcurvilinear )
219
    {
220
221
      xvals2 = (double*) Malloc(xsize2*ysize2*sizeof(double));
      yvals2 = (double*) Malloc(xsize2*ysize2*sizeof(double));
222
223
    }
    
224
225
  int *xoff = (int*) Malloc((nx+1)*sizeof(int));
  int *yoff = (int*) Malloc((ny+1)*sizeof(int));
226
227

  xoff[0] = 0;
228
  for ( int i = 0; i < nx; ++i )
229
    {
230
      int idx = xyinfo[i].id;
231
      if ( lregular ) memcpy(xvals2+xoff[i], xvals[idx], xsize[idx]*sizeof(double));
232
233
      xoff[i+1] = xoff[i] + xsize[idx];
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234

235
  yoff[0] = 0;
236
  for ( int j = 0; j < ny; ++j )
237
    {
238
      int idx = xyinfo[j*nx].id;
239
      if ( lregular ) memcpy(yvals2+yoff[j], yvals[idx], ysize[idx]*sizeof(double));
240
241
      yoff[j+1] = yoff[j] + ysize[idx];
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242

243
  if ( ginit == false )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
    {
245
      for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
	{
247
248
249
	  int idx = xyinfo[fileID].id;
	  int iy = fileID/nx;
	  int ix = fileID - iy*nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250

251
          int offset = yoff[iy]*xsize2 + xoff[ix];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
253
254
255
	  /*
	  printf("fileID %d %d, iy %d, ix %d, offset %d\n",
		 fileID, xyinfo[fileID].id, iy, ix, offset);
	  */
256
257
258
	  int ij = 0;
	  for ( int j = 0; j < ysize[idx]; ++j )
	    for ( int i = 0; i < xsize[idx]; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
	      {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
                if ( lcurvilinear )
261
262
263
264
                  {
                    xvals2[offset+j*xsize2+i] = xvals[idx][ij];
                    yvals2[offset+j*xsize2+i] = yvals[idx][ij];
                  }
265
                ef[idx].gridindex[ij++] = offset+j*xsize2+i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
268
269
270
271
272
	      }
	}
    }

  gridID2 = gridCreate(gridtype, xsize2*ysize2);
  gridDefXsize(gridID2, xsize2);
  gridDefYsize(gridID2, ysize2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
275
276
277
  if ( lregular || lcurvilinear )
    {
      gridDefXvals(gridID2, xvals2);
      gridDefYvals(gridID2, yvals2);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
278

279
280
281
282
283
284
  Free(xoff);
  Free(yoff);
  Free(xsize);
  Free(ysize);
  if ( xvals2 ) Free(xvals2);
  if ( yvals2 ) Free(yvals2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285

Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
287
288
289
290
291
292
293
294
  for ( int fileID = 0; fileID < nfiles; fileID++ )
    {
      if ( xvals[fileID] ) Free(xvals[fileID]);
      if ( yvals[fileID] ) Free(yvals[fileID]);
    }
  Free(xvals);
  Free(yvals);
  Free(xyinfo);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
  gridID = vlistGrid(ef[0].vlistID, igrid);
296
297

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

299
300
  if ( gridtype == GRID_PROJECTION ) grid_copy_mapping(gridID, gridID2);

301
  return gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
303
304
}


305
void *Collgrid(void *argument)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
{
307
  int nxblocks = -1;
308
  int varID, levelID;
309
  int nrecs0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310

Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
312
  cdoInitialize(argument);
    
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
314
  int nfiles = cdoStreamCnt() - 1;
  const char *ofilename = cdoStreamName(nfiles)->args;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315

316
317
  if ( !cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename) )
    cdoAbort("Outputfile %s already exists!", ofilename);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318

319
  ens_file_t *ef = (ens_file_t*) Malloc(nfiles*sizeof(ens_file_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320

321
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
      ef[fileID].streamID = streamOpenRead(cdoStreamName(fileID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
      ef[fileID].vlistID  = streamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
326
    }

327
328
  int vlistID1 = ef[0].vlistID;
  vlistClearFlag(vlistID1);
329
330

  /* check that the contents is always the same */
331
  for ( int fileID = 1; fileID < nfiles; fileID++ )
332
333
334
    vlistCompare(vlistID1, ef[fileID].vlistID, CMP_NAME | CMP_NLEVEL);

  int nvars = vlistNvars(vlistID1);
335
  bool *vars = (bool*) Malloc(nvars*sizeof(bool));
336
  for ( varID = 0; varID < nvars; varID++ ) vars[varID] = false;
337
  bool *vars1 = (bool*) Malloc(nvars*sizeof(bool));
338
  for ( varID = 0; varID < nvars; varID++ ) vars1[varID] = false;
339
340

  int nsel = operatorArgc();
341
  int noff = 0;
342
343
344
345
346
347
348
349
350

  if ( nsel > 0 )
    {
      int len = (int) strlen(operatorArgv()[0]);
      while ( --len >= 0 && isdigit(operatorArgv()[0][len]) ) ;

      if ( len == -1 )
        {
          nsel--;
351
          noff++;
352
353
354
355
          nxblocks = parameter2int(operatorArgv()[0]);
        }
    }

356
357
  if ( nsel == 0 )
    {
358
      for ( varID = 0; varID < nvars; varID++ ) vars1[varID] = true;
359
360
361
    }
  else
    {
362
      char **argnames = operatorArgv() + noff;
363
364
365
366
367

      if ( cdoVerbose )
	for ( int i = 0; i < nsel; i++ )
	  fprintf(stderr, "name %d = %s\n", i+1, argnames[i]);

368
369
      bool *selfound = (bool*) Malloc(nsel*sizeof(bool));
      for ( int i = 0; i < nsel; i++ ) selfound[i] = false;
370
371
372
373
374
375
376
377
378
379

      char varname[CDI_MAX_NAME];
      for ( varID = 0; varID < nvars; varID++ )
	{
	  vlistInqVarName(vlistID1, varID, varname);

	  for ( int isel = 0; isel < nsel; isel++ )
	    {
	      if ( strcmp(argnames[isel], varname) == 0 )
		{
380
381
		  selfound[isel] = true;
		  vars1[varID] = true;
382
383
384
385
386
		}
	    }
	}

      for ( int isel = 0; isel < nsel; isel++ )
387
	if ( selfound[isel] == false )
388
389
	  cdoAbort("Variable name %s not found!", argnames[isel]);

390
      Free(selfound);
391
392
393
394
    }

  for ( varID = 0; varID < nvars; varID++ )
    {
395
      if ( vars1[varID] )
396
397
398
399
400
401
402
403
	{
	  int zaxisID  = vlistInqVarZaxis(vlistID1, varID);
	  int nlevs    = zaxisInqSize(zaxisID);
	  for ( int levID = 0; levID < nlevs; levID++ )
	    vlistDefFlag(vlistID1, varID, levID, TRUE);
	}
    }

404
  for ( int fileID = 0; fileID < nfiles; fileID++ )
405
    {
406
407
408
409
      int gridsize = vlistGridsizeMax(ef[fileID].vlistID);
      ef[fileID].gridsize = gridsize;
      ef[fileID].gridindex = (int*) Malloc(gridsize*sizeof(int));
      ef[fileID].array = (double*) Malloc(gridsize*sizeof(double));
410
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411

412
413
  int vlistID2 = vlistCreate();
  vlistCopyFlag(vlistID2, vlistID1);
414
  /*
415
416
417
  if ( cdoVerbose )
    {
      vlistPrint(vlistID1);
418
      vlistPrint(vlistID2);
419
    }
420
  */
421
422
  //int vlistID2 = vlistDuplicate(vlistID1);
  int nvars2 = vlistNvars(vlistID2);
423
  // int *vars  = (int*) Malloc(nvars*sizeof(int));
424
  //for ( varID = 0; varID < nvars; varID++ ) vars[varID] = false;
425

426
427
428
  int ngrids1 = vlistNgrids(vlistID1);
  int ngrids2 = vlistNgrids(vlistID2);

429
  int *gridIDs = (int*) Malloc(ngrids2*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430

431
  bool ginit = false;
432
  for ( int i2 = 0; i2 < ngrids2; ++i2 )
433
    {
434
435
436
437
438
439
      int i1;
      for ( i1 = 0; i1 < ngrids1; ++i1 )
	if ( vlistGrid(vlistID1, i1) == vlistGrid(vlistID2, i2) ) break;

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

440
      if ( !ginit )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
	{
442
	  gridIDs[i2] = genGrid(ngrids2, nfiles, ef, ginit, i1, nxblocks);
443
	  if ( gridIDs[i2] != -1 ) ginit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
	}
445
      else
446
	gridIDs[i2] = genGrid(ngrids2, nfiles, ef, ginit, i1, nxblocks);
447
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448

449

450
451
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452
453
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
  int gridsize2 = 0;
455
  for ( int i = 0; i < ngrids2; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
    {
457
458
459
460
461
462
      if ( gridIDs[i] != -1 ) 
	{
	  if ( gridsize2 == 0 ) gridsize2 = gridInqSize(gridIDs[i]);
	  if ( gridsize2 != gridInqSize(gridIDs[i]) ) cdoAbort("gridsize differ!");
	  vlistChangeGridIndex(vlistID2, i, gridIDs[i]);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463
    }
464

465
  for ( varID = 0; varID < nvars2; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
466
    {
467
      int gridID = vlistInqVarGrid(vlistID2, varID);
468

469
      for ( int i = 0; i < ngrids2; ++i )
470
471
472
	{
	  if ( gridIDs[i] != -1 ) 
	    {
473
	      if ( gridID == vlistGrid(vlistID2, i) ) vars[varID] = true;
474
475
476
	      break;
	    }
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
478
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
  int streamID2 = streamOpenWrite(cdoStreamName(nfiles), cdoFiletype());
480
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
482
  streamDefVlist(streamID2, vlistID2);
	  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
  double *array2 = (gridsize2 > 0) ? (double*) Malloc(gridsize2*sizeof(double)) : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484

Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
487
488
  do
    {
      nrecs0 = streamInqTimestep(ef[0].streamID, tsID);
489
      for ( int fileID = 1; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490
	{
491
	  int nrecs = streamInqTimestep(ef[fileID].streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
492
	  if ( nrecs != nrecs0 )
493
	    cdoAbort("Number of records at time step %d of %s and %s differ!", tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
495
496
497
498
499
	}

      taxisCopyTimestep(taxisID2, taxisID1);

      if ( nrecs0 > 0 ) streamDefTimestep(streamID2, tsID);
      
500
      for ( int recID = 0; recID < nrecs0; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
	  streamInqRecord(ef[0].streamID, &varID, &levelID);
503
	  if ( cdoVerbose && tsID == 0 ) printf(" tsID, recID, varID, levelID %d %d %d %d\n", tsID, recID, varID, levelID);
504

505
	  for ( int fileID = 1; fileID < nfiles; fileID++ )
506
507
	    {
	      int varIDx, levelIDx;
508
	      streamInqRecord(ef[fileID].streamID, &varIDx, &levelIDx);
509
510
	    }

511
512
513
514
515
	  if ( vlistInqFlag(vlistID1, varID, levelID) == TRUE )
	    {
	      int varID2   = vlistFindVar(vlistID2, varID);
	      int levelID2 = vlistFindLevel(vlistID2, varID, levelID);
	      if ( cdoVerbose && tsID == 0 ) printf("varID %d %d levelID %d %d\n", varID, varID2, levelID, levelID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
516

517
	      double missval = vlistInqVarMissval(vlistID2, varID2);
518
	      for ( int i = 0; i < gridsize2; i++ ) array2[i] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519

Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
#if defined(_OPENMP)
521
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
522
#endif
523
	      for ( int fileID = 0; fileID < nfiles; fileID++ )
524
		{
525
		  streamReadRecord(ef[fileID].streamID, ef[fileID].array, &ef[fileID].nmiss);
526
527
528

		  if ( vars[varID2] )
		    {
529
                      int gridsize = ef[fileID].gridsize;
530
		      for ( int i = 0; i < gridsize; ++i )
531
			array2[ef[fileID].gridindex[i]] = ef[fileID].array[i];
532
		    }
533
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534

535
	      streamDefRecord(streamID2, varID2, levelID2);
536

537
538
	      if ( vars[varID2] )
		{
539
		  int nmiss = 0;
540
541
542
543
544
545
		  for ( int i = 0; i < gridsize2; i++ )
		    if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;

		  streamWriteRecord(streamID2, array2, nmiss);
		}
	      else
546
		streamWriteRecord(streamID2, ef[0].array, ef[0].nmiss);
547
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
549
550
551
552
553
	}

      tsID++;
    }
  while ( nrecs0 > 0 );

554
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555
    streamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
557
558

  streamClose(streamID2);

559
  for ( int fileID = 0; fileID < nfiles; fileID++ )
560
561
562
563
    {
      if ( ef[fileID].gridindex ) Free(ef[fileID].gridindex);
      if ( ef[fileID].array ) Free(ef[fileID].array);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564

565
566
  if ( ef ) Free(ef);
  if ( array2 ) Free(array2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567

568
569
570
  Free(gridIDs);
  if ( vars   ) Free(vars);
  if ( vars1  ) Free(vars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
573

  cdoFinish();

574
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
}