Ensstat.cc 17.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-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
19
20
/*
   This module contains the following operators:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
      Ensstat    ensrange        Ensemble range
22
23
24
25
26
27
      Ensstat    ensmin          Ensemble minimum
      Ensstat    ensmax          Ensemble maximum
      Ensstat    enssum          Ensemble sum
      Ensstat    ensmean         Ensemble mean
      Ensstat    ensavg          Ensemble average
      Ensstat    ensstd          Ensemble standard deviation
28
      Ensstat    ensstd1         Ensemble standard deviation
29
      Ensstat    ensvar          Ensemble variance
30
      Ensstat    ensvar1         Ensemble variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
      Ensstat    enspctl         Ensemble percentiles
32
33
*/

Ralf Mueller's avatar
Ralf Mueller committed
34
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
#include "cdo.h"
#include "cdo_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
#include "cdo_task.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
#include "util.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
41


Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
43
44
45
typedef struct
{
  int streamID;
  int vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
47
48
  int nmiss[2];
  double missval[2];
  double *array[2];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
} ens_file_t;


typedef struct
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
56
  int t;
  int varID[2];
  int levelID[2];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
  int vlistID1;
  int streamID2;
  int nfiles;
  ens_file_t *ef;
  double *array2;
  double *count2;
  field_type *field;
  int operfunc;
  double pn;
  bool lpctl;
  bool count_data;
  int nvars;
} ensstat_arg_t;


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
void *ensstat_func(void *ensarg)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
  if ( CDO_task ) cdo_omp_set_num_threads(ompNumThreads);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76

Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
  ensstat_arg_t *arg = (ensstat_arg_t*) ensarg;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
  int t = arg->t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
  int nfiles = arg->nfiles;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
  ens_file_t *ef = arg->ef;
  field_type *field = arg->field;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
  bool lmiss = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
  for ( int fileID = 0; fileID < nfiles; fileID++ ) if ( ef[fileID].nmiss[t] > 0 ) lmiss = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85

Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
  int gridID = vlistInqVarGrid(arg->vlistID1, arg->varID[t]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
  int gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
  double missval = vlistInqVarMissval(arg->vlistID1, arg->varID[t]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
91
92
93
94
95
96
97

  int nmiss = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(shared)
#endif
  for ( int i = 0; i < gridsize; ++i )
    {
      int ompthID = cdo_omp_get_thread_num();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
99
      field[ompthID].missval = missval;
      field[ompthID].nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
      for ( int fileID = 0; fileID < nfiles; fileID++ )
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
          field[ompthID].ptr[fileID] = ef[fileID].array[t][i];
          if ( lmiss && DBL_IS_EQUAL(field[ompthID].ptr[fileID], ef[fileID].missval[t]) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
              field[ompthID].ptr[fileID] = missval;
              field[ompthID].nmiss++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
109
            }
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
      arg->array2[i] = arg->lpctl ? fldpctl(field[ompthID], arg->pn) : fldfun(field[ompthID], arg->operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111

Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
      if ( DBL_IS_EQUAL(arg->array2[i], field[ompthID].missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
114
115
116
117
118
119
        {
#if defined(_OPENMP)
#include "pragma_omp_atomic_update.h"
#endif
          nmiss++;
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
      if ( arg->count_data ) arg->count2[i] = nfiles - field[ompthID].nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
    }

123
124
  pstreamDefRecord(arg->streamID2, arg->varID[t], arg->levelID[t]);
  pstreamWriteRecord(arg->streamID2, arg->array2, nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
127

  if ( arg->count_data )
    {
128
129
      pstreamDefRecord(arg->streamID2, arg->varID[t]+arg->nvars, arg->levelID[t]);
      pstreamWriteRecord(arg->streamID2, arg->count2, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132

  return NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
134
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
137
138
139
140
141
142
143
144
//#define TEST_TASK
#ifdef TEST_TASK
void *Ensstat(void *argument)
{
  void *task = CDO_task ? cdo_task_new() : NULL;
  ensstat_arg_t ensstat_arg;
  int nrecs0;

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
147
148
149
150
151
152
153
154
155
156
  cdoOperatorAdd("ensrange", func_range, 0, NULL);
  cdoOperatorAdd("ensmin",   func_min,   0, NULL);
  cdoOperatorAdd("ensmax",   func_max,   0, NULL);
  cdoOperatorAdd("enssum",   func_sum,   0, NULL);
  cdoOperatorAdd("ensmean",  func_mean,  0, NULL);
  cdoOperatorAdd("ensavg",   func_avg,   0, NULL);
  cdoOperatorAdd("ensstd",   func_std,   0, NULL);
  cdoOperatorAdd("ensstd1",  func_std1,  0, NULL);
  cdoOperatorAdd("ensvar",   func_var,   0, NULL);
  cdoOperatorAdd("ensvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("enspctl",  func_pctl,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192

  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);

  bool lpctl = operfunc == func_pctl;

  int argc = operatorArgc();
  int nargc = argc;

  double pn = 0;
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
      pn = parameter2int(operatorArgv()[0]);
      percentile_check_number(pn);
      argc--;
    }

  bool count_data = false;
  if ( argc == 1 )
    {
      if ( strcmp("count", operatorArgv()[nargc-1]) == 0 ) count_data = true;
      else cdoAbort("Unknown parameter: >%s<", operatorArgv()[nargc-1]); 
    }
    
  int nfiles = cdoStreamCnt() - 1;

  if ( cdoVerbose ) cdoPrint("Ensemble over %d files.", nfiles);

  const char *ofilename = cdoStreamName(nfiles)->args;

  if ( !cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename) )
    cdoAbort("Outputfile %s already exists!", ofilename);

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
196
197
198
199
200
201
202
203
  field_type *field = (field_type *) Malloc(ompNumThreads*sizeof(field_type));
  for ( int i = 0; i < ompNumThreads; i++ )
    {
      field_init(&field[i]);
      field[i].size = nfiles;
      field[i].ptr  = (double*) Malloc(nfiles*sizeof(double));
    }

  for ( int fileID = 0; fileID < nfiles; fileID++ )
    {
204
205
      ef[fileID].streamID = pstreamOpenRead(cdoStreamName(fileID));
      ef[fileID].vlistID  = pstreamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
    }

  /* check that the contents is always the same */
  for ( int fileID = 1; fileID < nfiles; fileID++ )
    vlistCompare(ef[0].vlistID, ef[fileID].vlistID, CMP_ALL);

  int vlistID1 = ef[0].vlistID;
  int vlistID2 = vlistDuplicate(vlistID1);
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
  vlistDefTaxis(vlistID2, taxisID2);

  int gridsizemax = vlistGridsizeMax(vlistID1);

  for ( int fileID = 0; fileID < nfiles; fileID++ )
    {
      ef[fileID].array[0] = (double*) Malloc(gridsizemax*sizeof(double));
      ef[fileID].array[1] = task ? (double*) Malloc(gridsizemax*sizeof(double)) : NULL;
    }

  double *array2 = (double *) Malloc(gridsizemax*sizeof(double));

  int nvars = vlistNvars(vlistID2);
  double *count2 = NULL;
  if ( count_data )
    {
      count2 = (double *) Malloc(gridsizemax*sizeof(double));
      for ( int varID = 0; varID < nvars; ++varID )
	{
	  char name[CDI_MAX_NAME];
	  vlistInqVarName(vlistID2, varID, name);
	  strcat(name, "_count");
	  int gridID = vlistInqVarGrid(vlistID2, varID);
	  int zaxisID = vlistInqVarZaxis(vlistID2, varID);
	  int tsteptype = vlistInqVarTsteptype(vlistID2, varID);
	  int cvarID = vlistDefVar(vlistID2, gridID, zaxisID, tsteptype);
	  vlistDefVarName(vlistID2, cvarID, name);
	  vlistDefVarDatatype(vlistID2, cvarID, CDI_DATATYPE_INT16);
	  if ( cvarID != (varID+nvars) ) cdoAbort("Internal error, varIDs do not match!");
	}
    }

248
249
  int streamID2 = pstreamOpenWrite(cdoStreamName(nfiles), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

  ensstat_arg.vlistID1 = vlistID1;
  ensstat_arg.streamID2 = streamID2;
  ensstat_arg.nfiles = nfiles;
  ensstat_arg.array2 = array2;
  ensstat_arg.count2 = count2;
  ensstat_arg.field = field;
  ensstat_arg.operfunc = operfunc;
  ensstat_arg.pn = pn;
  ensstat_arg.lpctl = lpctl;
  ensstat_arg.count_data = count_data;
  ensstat_arg.nvars = nvars;
  ensstat_arg.t = 0;

  bool lwarning = false;
  bool lerror = false;
  int t = 0;
  int tsID = 0;
  do
    {
270
      nrecs0 = pstreamInqTimestep(ef[0].streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271
272
273
      for ( int fileID = 1; fileID < nfiles; fileID++ )
	{
	  int streamID = ef[fileID].streamID;
274
	  int nrecs = pstreamInqTimestep(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
	  if ( nrecs != nrecs0 )
	    {
	      if ( nrecs == 0 )
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(fileID)->args);
                }
	      else if ( nrecs0 == 0 )
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(0)->args);
                }
	      else
                {
                  lerror = true;
                  cdoWarning("Inconsistent ensemble file, number of records at time step %d of %s and %s differ!",
                             tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
                }
              goto CLEANUP;
	    }
	}

      if ( nrecs0 > 0 )
	{
	  taxisCopyTimestep(taxisID2, taxisID1);
300
	  pstreamDefTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
302
303
304
305
306
307
308
	}

      for ( int recID = 0; recID < nrecs0; recID++ )
	{
          int varID = 0, levelID;

	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
309
	      pstreamInqRecord(ef[fileID].streamID, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310
311
312
313
314
              ef[fileID].missval[t] = vlistInqVarMissval(ef[fileID].vlistID, varID);
	    }
          //#pragma omp parallel for default(none) shared(ef, nfiles)
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
315
	      pstreamReadRecord(ef[fileID].streamID, ef[fileID].array[t], &ef[fileID].nmiss[t]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
	    }

          ensstat_arg.ef = ef;
          ensstat_arg.varID[t] = varID;
          ensstat_arg.levelID[t] = levelID;
          if ( task )
            {
              cdo_task_start(task, ensstat_func, &ensstat_arg);
              cdo_task_wait(task);
              t = !t;
            }
          else
            {
              ensstat_func(&ensstat_arg);
            }
        }

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

 CLEANUP:

  if ( lwarning ) cdoWarning("Inconsistent ensemble, processed only the first %d timesteps!", tsID);
  if ( lerror ) cdoAbort("Inconsistent ensemble, processed only the first %d timesteps!", tsID);

  if ( task ) cdo_task_wait(task);

  for ( int fileID = 0; fileID < nfiles; fileID++ )
345
    pstreamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
346

347
  pstreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368

  for ( int fileID = 0; fileID < nfiles; fileID++ )
    {
      if ( ef[fileID].array[0] ) Free(ef[fileID].array[0]);
      if ( ef[fileID].array[1] ) Free(ef[fileID].array[1]);
    }

  if ( ef ) Free(ef);
  if ( array2 ) Free(array2);
  if ( count2 ) Free(count2);

  for ( int i = 0; i < ompNumThreads; i++ ) if ( field[i].ptr ) Free(field[i].ptr);
  if ( field ) Free(field);

  if ( task ) cdo_task_delete(task);

  cdoFinish();

  return 0;
}
#else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
370
void *Ensstat(void *argument)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
  void *task = CDO_task ? cdo_task_new() : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
  ensstat_arg_t ensstat_arg;
373
  int nrecs0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
375
376

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
377
378
379
380
381
382
383
384
385
386
387
  cdoOperatorAdd("ensrange", func_range, 0, NULL);
  cdoOperatorAdd("ensmin",   func_min,   0, NULL);
  cdoOperatorAdd("ensmax",   func_max,   0, NULL);
  cdoOperatorAdd("enssum",   func_sum,   0, NULL);
  cdoOperatorAdd("ensmean",  func_mean,  0, NULL);
  cdoOperatorAdd("ensavg",   func_avg,   0, NULL);
  cdoOperatorAdd("ensstd",   func_std,   0, NULL);
  cdoOperatorAdd("ensstd1",  func_std1,  0, NULL);
  cdoOperatorAdd("ensvar",   func_var,   0, NULL);
  cdoOperatorAdd("ensvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("enspctl",  func_pctl,  0, NULL);
388

389
390
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391

392
393
  bool lpctl = operfunc == func_pctl;

394
395
  int argc = operatorArgc();
  int nargc = argc;
396
397

  double pn = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
399
400
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
401
      pn = parameter2int(operatorArgv()[0]);
402
      percentile_check_number(pn);
403
404
405
      argc--;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
  bool count_data = false;
407
408
  if ( argc == 1 )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
      if ( strcmp("count", operatorArgv()[nargc-1]) == 0 ) count_data = true;
410
      else cdoAbort("Unknown parameter: >%s<", operatorArgv()[nargc-1]); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
412
    }
    
413
  int nfiles = cdoStreamCnt() - 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
414

415
  if ( cdoVerbose ) cdoPrint("Ensemble over %d files.", nfiles);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416

417
  const char *ofilename = cdoStreamName(nfiles)->args;
418

419
420
  if ( !cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename) )
    cdoAbort("Outputfile %s already exists!", ofilename);
421

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

424
  field_type *field = (field_type *) Malloc(ompNumThreads*sizeof(field_type));
425
  for ( int i = 0; i < ompNumThreads; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
      field_init(&field[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
429
      field[i].size = nfiles;
      field[i].ptr  = (double*) Malloc(nfiles*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430
431
    }

432
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
    {
434
435
      ef[fileID].streamID = pstreamOpenRead(cdoStreamName(fileID));
      ef[fileID].vlistID  = pstreamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
437
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
  /* check that the contents is always the same */
439
  for ( int fileID = 1; fileID < nfiles; fileID++ )
440
    vlistCompare(ef[0].vlistID, ef[fileID].vlistID, CMP_ALL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441

442
443
444
445
  int vlistID1 = ef[0].vlistID;
  int vlistID2 = vlistDuplicate(vlistID1);
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446
447
  vlistDefTaxis(vlistID2, taxisID2);

448
  int gridsizemax = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449

450
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
451
452
453
454
    {
      ef[fileID].array[0] = (double*) Malloc(gridsizemax*sizeof(double));
      ef[fileID].array[1] = task ? (double*) Malloc(gridsizemax*sizeof(double)) : NULL;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455

456
  double *array2 = (double *) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457

458
459
460
461
  int nvars = vlistNvars(vlistID2);
  double *count2 = NULL;
  if ( count_data )
    {
462
      count2 = (double *) Malloc(gridsizemax*sizeof(double));
463
      for ( int varID = 0; varID < nvars; ++varID )
464
465
466
467
	{
	  char name[CDI_MAX_NAME];
	  vlistInqVarName(vlistID2, varID, name);
	  strcat(name, "_count");
468
	  int gridID = vlistInqVarGrid(vlistID2, varID);
469
470
471
472
	  int zaxisID = vlistInqVarZaxis(vlistID2, varID);
	  int tsteptype = vlistInqVarTsteptype(vlistID2, varID);
	  int cvarID = vlistDefVar(vlistID2, gridID, zaxisID, tsteptype);
	  vlistDefVarName(vlistID2, cvarID, name);
473
	  vlistDefVarDatatype(vlistID2, cvarID, CDI_DATATYPE_INT16);
474
475
476
477
	  if ( cvarID != (varID+nvars) ) cdoAbort("Internal error, varIDs do not match!");
	}
    }

478
479
  int streamID2 = pstreamOpenWrite(cdoStreamName(nfiles), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
480

Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
482
483
484
485
486
487
488
489
490
491
  ensstat_arg.vlistID1 = vlistID1;
  ensstat_arg.streamID2 = streamID2;
  ensstat_arg.nfiles = nfiles;
  ensstat_arg.array2 = array2;
  ensstat_arg.count2 = count2;
  ensstat_arg.field = field;
  ensstat_arg.operfunc = operfunc;
  ensstat_arg.pn = pn;
  ensstat_arg.lpctl = lpctl;
  ensstat_arg.count_data = count_data;
  ensstat_arg.nvars = nvars;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
492
  ensstat_arg.t = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
493

Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
495
  bool lwarning = false;
  bool lerror = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
496
  int t = 0;
497
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
498
499
  do
    {
500
      nrecs0 = pstreamInqTimestep(ef[0].streamID, tsID);
501
      for ( int fileID = 1; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
	{
503
	  int streamID = ef[fileID].streamID;
504
	  int nrecs = pstreamInqTimestep(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505
	  if ( nrecs != nrecs0 )
506
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
	      if ( nrecs == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508
509
510
511
512
513
514
515
516
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(fileID)->args);
                }
	      else if ( nrecs0 == 0 )
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(0)->args);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
	      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
518
519
520
521
522
523
                {
                  lerror = true;
                  cdoWarning("Inconsistent ensemble file, number of records at time step %d of %s and %s differ!",
                             tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
                }
              goto CLEANUP;
524
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
526
	}

527
528
529
      if ( nrecs0 > 0 )
	{
	  taxisCopyTimestep(taxisID2, taxisID1);
530
	  pstreamDefTimestep(streamID2, tsID);
531
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532

533
      for ( int recID = 0; recID < nrecs0; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
	{
535
          int varID = 0, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536

537
538
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
539
	      pstreamInqRecord(ef[fileID].streamID, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
              ef[fileID].missval[t] = vlistInqVarMissval(ef[fileID].vlistID, varID);
541
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
          //#pragma omp parallel for default(none) shared(ef, nfiles)
543
544
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
545
	      pstreamReadRecord(ef[fileID].streamID, ef[fileID].array[t], &ef[fileID].nmiss[t]);
546
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547

Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
          ensstat_arg.ef = ef;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
          ensstat_arg.varID[t] = varID;
          ensstat_arg.levelID[t] = levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
552
553
554
555
556
557
558
559
          if ( task )
            {
              cdo_task_start(task, ensstat_func, &ensstat_arg);
              cdo_task_wait(task);
            }
          else
            {
              ensstat_func(&ensstat_arg);
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561
562
563
564
565

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
567
568
569
570
 CLEANUP:

  if ( lwarning ) cdoWarning("Inconsistent ensemble, processed only the first %d timesteps!", tsID);
  if ( lerror ) cdoAbort("Inconsistent ensemble, processed only the first %d timesteps!", tsID);

571
  for ( int fileID = 0; fileID < nfiles; fileID++ )
572
    pstreamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
573

574
  pstreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575

576
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
579
580
    {
      if ( ef[fileID].array[0] ) Free(ef[fileID].array[0]);
      if ( ef[fileID].array[1] ) Free(ef[fileID].array[1]);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
581

582
583
584
  if ( ef ) Free(ef);
  if ( array2 ) Free(array2);
  if ( count2 ) Free(count2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585

Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
  for ( int i = 0; i < ompNumThreads; i++ ) if ( field[i].ptr ) Free(field[i].ptr);
587
  if ( field ) Free(field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588

Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
  if ( task ) cdo_task_delete(task);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
590

Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
592
  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
594
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
595
#endif