ecacore.cc 48.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2006 Brockmann Consult
  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.
*/

#include <assert.h>
#include <string.h>

Ralf Mueller's avatar
Ralf Mueller committed
21
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
22

23
#include "functs.h"
24
#include "cdo_vlist.h"
25
#include "cdo_cdi_wrapper.h"
26
#include <mpim_grid.h>
27
#include "process_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
#include "ecacore.h"
29
#include "ecautil.h"
30
#include "util_files.h"
31
#include "util_date.h"
32
#include "datetime.h"
33

34
35
#define FIRST_VAR_ID 0

Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
#define IS_NOT_SET(x) (x == nullptr)
#define IS_SET(x) (x != nullptr)
38
#define VIS_SET(x) (!x.empty())
39

40
41
void
eca1(const ECA_REQUEST_1 *request)
42
{
43
  const auto operatorID = cdoOperatorID();
44
45

  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
46
  int64_t ivdate = 0, ovdate = 0, indate21 = 0;
47
  int ivtime = 0, ovtime = 0;
48
  int nrecs;
49
  int levelID;
50
51
52
  int itsID;
  int otsID;
  long nsets;
53

54
  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);
55

56
  const auto istreamID = cdoOpenRead(0);
57

58
59
  const auto ivlistID = cdoStreamInqVlist(istreamID);
  const auto ovlistID = vlistCreate();
60

61
62
63
  const auto gridID = vlistInqVarGrid(ivlistID, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID, FIRST_VAR_ID);
  const auto missval = vlistInqVarMissval(ivlistID, FIRST_VAR_ID);
64

65
  auto varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);
66
67

  vlistDefVarMissval(ovlistID, varID, missval);
68

69
70
71
  if (IS_SET(request->var1.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var1.name);
  if (IS_SET(request->var1.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var1.longname);
  if (IS_SET(request->var1.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var1.units);
72

73
  if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
74
    {
75
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);
76

77
78
      vlistDefVarMissval(ovlistID, varID, missval);

79
80
81
      if (IS_SET(request->var2.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var2.name);
      if (IS_SET(request->var2.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var2.longname);
      if (IS_SET(request->var2.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var2.units);
82
    }
83
84

  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);
85

86
87
  const auto itaxisID = vlistInqTaxis(ivlistID);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
88
89
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
90
91
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID));
92
  taxisDefRdate(otaxisID, request->var1.refdate);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
  taxisDefRtime(otaxisID, 0);
94
95
  vlistDefTaxis(ovlistID, otaxisID);

96
  const auto ostreamID = cdoOpenWrite(1);
97
  cdoDefVlist(ostreamID, ovlistID);
98

99
  const auto gridsize = gridInqSize(gridID);
100

101
  Field field1, field2, field3;
102
  field1.resize(gridsize);
103
  field3.resize(gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
  if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3)) field2.resize(gridsize);
105

106
  const auto nlevels = zaxisInqSize(zaxisID);
107

108
109
110
111
112
  FieldVector var12(nlevels), samp1(nlevels), samp2(nlevels);
  FieldVector var13, var21, var23;
  if (IS_SET(request->var1.f3)) var13.resize(nlevels);
  if (IS_SET(request->var2.h2)) var21.resize(nlevels);
  if (IS_SET(request->var2.h3)) var23.resize(nlevels);
113
114

  for (levelID = 0; levelID < nlevels; levelID++)
115
    {
116
      var12[levelID].grid = gridID;
117
      var12[levelID].missval = missval;
118
      var12[levelID].resize(gridsize);
119

120
      samp1[levelID].grid = gridID;
121
      samp1[levelID].missval = missval;
122
      samp1[levelID].resize(gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123

124
      samp2[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
      samp2[levelID].missval = missval;
126

127
      if (IS_SET(request->var1.f3))
128
        {
129
          var13[levelID].grid = gridID;
130
          var13[levelID].missval = missval;
131
          var13[levelID].resize(gridsize);
132
        }
133
      if (IS_SET(request->var2.h2))
134
        {
135
          var21[levelID].grid = gridID;
136
          var21[levelID].missval = missval;
137
          var21[levelID].resize(gridsize);
138
        }
139
      if (IS_SET(request->var2.h3))
140
        {
141
          var23[levelID].grid = gridID;
142
          var23[levelID].missval = missval;
143
          var23[levelID].resize(gridsize);
144
145
146
        }
    }

147
148
  itsID = 0;
  otsID = 0;
149
  while (true)
150
151
    {
      nsets = 0;
152
      while ((nrecs = cdoStreamInqTimestep(istreamID, itsID)) > 0)
153
154
155
156
        {
          ivdate = taxisInqVdate(itaxisID);
          ivtime = taxisInqVtime(itaxisID);

157
158
159
160
161
          if (nsets == 0)
            {
              SET_DATE(indate2, ivdate, ivtime);
              indate21 = ivdate;
            }
162
          SET_DATE(indate1, ivdate, ivtime);
163

164
          if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;
165

166
          for (int recID = 0; recID < nrecs; recID++)
167
            {
168
              cdoInqRecord(istreamID, &varID, &levelID);
169

170
              if (varID != FIRST_VAR_ID) continue;
171

172
              if (nsets == 0)
173
                {
174
175
176
177
178
179
180
181
182
183
                  if ( cmplen != 0 && itsID == 0 && request->var1.f2 == &vfarnum2 )
                    {
                      fieldFill(var12[levelID], missval);
                      var12[levelID].nmiss = gridsize;
                    }
                  else if ( request->var1.f2 != &vfarnum2 )
                    {
                      fieldFill(var12[levelID], missval);
                      var12[levelID].nmiss = gridsize;
                    }
184
185
186
187
188
                  fieldFill(samp1[levelID], missval);
                  if (!samp2[levelID].empty()) fieldFill(samp2[levelID], 0.0);
                  if (IS_SET(request->var1.f3)) fieldFill(var13[levelID], missval);
                  if (IS_SET(request->var2.h2)) fieldFill(var21[levelID], missval);
                  if (IS_SET(request->var2.h3)) fieldFill(var23[levelID], missval);
189

Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
                  samp1[levelID].nmiss = gridsize;
191
192
193
                  if (IS_SET(request->var1.f3)) var13[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h2)) var21[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h3)) var23[levelID].nmiss = gridsize;
194
195
                }

196
              cdoReadRecord(istreamID, field1.vec_d.data(), &field1.nmiss);
197
              field1.grid = var12[levelID].grid;
198
              field1.missval = var12[levelID].missval;
199

200
              vfarnum(samp1[levelID], field1);
201

202
              if (IS_SET(request->var2.h2))
203
                {
204
                  field2.vec_d = field1.vec_d;
205
206
                  field2.nmiss = field1.nmiss;
                  field2.grid = field1.grid;
207
208
                  field2.missval = field1.missval;
                }
209

210
              if (IS_SET(request->var1.f1)) request->var1.f1(field1, request->var1.f1arg);
211

212
              if (field1.nmiss || !samp2[levelID].empty())
213
                {
214
                  if (samp2[levelID].empty())
215
                    {
216
                      samp2[levelID].resize(gridsize);
217
                      fieldFill(samp2[levelID], nsets);
218
                    }
219
                  for (size_t i = 0; i < gridsize; i++)
220
                    {
221
222
                      if (DBL_IS_EQUAL(field1.vec_d[i], field1.missval)) continue;
                      samp2[levelID].vec_d[i]++;
223
224
                    }
                }
225

226
227
              if (IS_NOT_EQUAL(request->var1.mulc, 0.0)) vfarcmul(field1, request->var1.mulc);
              if (IS_NOT_EQUAL(request->var1.addc, 0.0)) vfarcadd(field1, request->var1.addc);
228

229
230
              if (IS_SET(request->var1.f3) && request->var1.f2 == &vfarnum2)
                {
231
                  varrayCopy(gridsize, var12[levelID].vec_d, field3.vec_d);
232
233
234
235
                  field3.nmiss = var12[levelID].nmiss;
                  field3.grid = var12[levelID].grid;
                  field3.missval = var12[levelID].missval;
                }
236
              request->var1.f2(var12[levelID], field1);
237
238
239

              if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
                {
240
                  // if h2 is null, use the output of f2 as input for h1
241
                  if (IS_NOT_SET(request->var2.h2))
242
                    {
243
                      varrayCopy(gridsize, var12[levelID].vec_d, field2.vec_d);
244
245
                      field2.nmiss = var12[levelID].nmiss;
                      field2.grid = var12[levelID].grid;
246
247
                      field2.missval = var12[levelID].missval;
                    }
248

249
                  if (IS_SET(request->var2.h1)) request->var2.h1(field2, request->var2.h1arg);
250
251

                  if (IS_NOT_SET(request->var2.h2))
252
                    request->var2.h3(var23[levelID], field2);
253
254
                  else
                    {
255
256
                      request->var2.h2(var21[levelID], field2);
                      if (IS_SET(request->var2.h3)) request->var2.h3(var23[levelID], var21[levelID]);
257
258
259
                    }
                }

260
261
262
263
              if (IS_SET(request->var1.f3))
                {
                  if (cmplen != 0 && request->var1.f2 == &vfarnum2)
                    {
264
265
                      auto &array1 = field3.vec_d;
                      const auto &array2 = var12[levelID].vec_d;
266
                      const auto missval2 = field2.missval;
267
                
268
                      const auto len = field1.size;
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
                      if (len != field2.size) cdoAbort("Fields have different size (%s)", __func__);
                
                      if (field2.nmiss)
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], missval2) && !DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      else
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      request->var1.f3(var13[levelID], field3);
                    }
                  else
                    request->var1.f3(var13[levelID], var12[levelID]);
                }
292
293
294
295
296
297
298
299
            }

          ovdate = ivdate;
          ovtime = ivtime;
          nsets++;
          itsID++;
        }

300
301
      if (nrecs == 0 && nsets == 0) break;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
      if (request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME)
303
        for (levelID = 0; levelID < nlevels; levelID++)
304
          {
305
            auto &rvar = IS_SET(request->var1.f3) ? var13[levelID] : var12[levelID];
306

307
308
            if (samp2[levelID].empty())
              vfarcdiv(rvar, nsets);
309
            else
310
              vfardiv(rvar, samp2[levelID]);
311

312
            if (request->var1.epilog == PERCENT_OF_TIME) vfarcmul(rvar, 100.0);
313
          }
314
315
316
317
318
319
320
321
322
323
324
325
326
      if ( request->var1.refdate == 19550101 )
        {
          taxisDefVdate(otaxisID, ovdate);
          taxisDefVtime(otaxisID, ovtime);
          cdoDefTimestep(ostreamID, otsID);
        }
      else
        {
          int year, month, day;
          cdiDecodeDate(indate21, &year, &month, &day);
          defineMidOfTime(cdoOperatorF2(operatorID), otaxisID, year, month, 12);
          cdoDefTimestep(ostreamID, otsID);
        }
327

Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
      if (otsID && vlistInqVarTimetype(ivlistID, FIRST_VAR_ID) == TIME_CONSTANT) continue;
329
330

      varID = 0;
331
332
      for (levelID = 0; levelID < nlevels; levelID++)
        {
333
          auto &rvar = IS_SET(request->var1.f3) ? var13[levelID] : var12[levelID];
334

335
          vfarsel(rvar, samp1[levelID]);
336

337
          cdoDefRecord(ostreamID, varID, levelID);
338
          cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
339
340
341
342
343
344
345
        }

      if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
        {
          varID = 1;
          for (levelID = 0; levelID < nlevels; levelID++)
            {
346
              auto &rvar = IS_SET(request->var2.h3) ? var23[levelID] : var21[levelID];
347

348
              vfarsel(rvar, samp1[levelID]);
349

350
              cdoDefRecord(ostreamID, varID, levelID);
351
              cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
352
            }
353
354
        }

355
      if (nrecs == 0) break;
356
      otsID++;
357
358
    }

359
360
  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID);
361
362
}

363
364
void
eca2(const ECA_REQUEST_2 *request)
365
{
366
  const auto operatorID = cdoOperatorID();
367
368

  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
369
  int64_t ivdate = 0, ovdate = 0, indate21 = 0;
370
  int ivtime = 0, ovtime = 0;
371
  int nrecs;
372
  int varID, levelID;
373
374
375
  int itsID;
  int otsID;
  long nsets;
376

377
  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);
378

379
380
  const auto istreamID1 = cdoOpenRead(0);
  const auto istreamID2 = cdoOpenRead(1);
381

382
383
384
  const auto ivlistID1 = cdoStreamInqVlist(istreamID1);
  const auto ivlistID2 = cdoStreamInqVlist(istreamID2);
  const auto ovlistID = vlistCreate();
385

Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);
387

388
389
390
391
  const auto gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  const auto missval1 = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);
  const auto missval2 = vlistInqVarMissval(ivlistID2, FIRST_VAR_ID);
392

393
394
  varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

395
396
  vlistDefVarMissval(ovlistID, varID, missval1);

397
398
399
  if (IS_SET(request->var1.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var1.name);
  if (IS_SET(request->var1.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var1.longname);
  if (IS_SET(request->var1.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var1.units);
400

401
  if (IS_SET(request->var2.h2))
402
    {
403
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);
404

405
406
      vlistDefVarMissval(ovlistID, varID, missval1);

407
408
409
      if (IS_SET(request->var2.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var2.name);
      if (IS_SET(request->var2.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var2.longname);
      if (IS_SET(request->var2.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var2.units);
410
411
    }

412
  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);
413

414
415
  const auto itaxisID1 = vlistInqTaxis(ivlistID1);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
416
417
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
418
419
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
420
  taxisDefRdate(otaxisID, request->var1.refdate);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
  taxisDefRtime(otaxisID, 0);
422
423
  vlistDefTaxis(ovlistID, otaxisID);

424
  const auto ostreamID = cdoOpenWrite(2);
425
  cdoDefVlist(ostreamID, ovlistID);
426

427
  const auto gridsize = gridInqSize(gridID);
428

429
  Field field1, field2, field3;
430
431
  field1.resize(gridsize);
  field2.resize(gridsize);
432
433
434
  field3.resize(gridsize);
  constexpr int MaxDays = 373;
  FieldVector2D vars2[MaxDays];
435

436
  const auto nlevels = zaxisInqSize(zaxisID);
437

438
439
440
441
442
  FieldVector var14(nlevels), samp1(nlevels), samp2(nlevels), samp3(nlevels);
  FieldVector total, var15, var22;
  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) total.resize(nlevels);
  if (IS_SET(request->var1.f5)) var15.resize(nlevels);
  if (IS_SET(request->var2.h2)) var22.resize(nlevels);
443
444

  for (levelID = 0; levelID < nlevels; levelID++)
445
    {
446
      var14[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
      var14[levelID].missval = missval1;
448
      var14[levelID].resize(gridsize);
449
450

      samp1[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
451
      samp1[levelID].missval = missval1;
452
      samp1[levelID].resize(gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453

454
      samp2[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
      samp2[levelID].missval = missval1;
456
      samp2[levelID].resize(gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457

458
      samp3[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
459
      samp3[levelID].missval = missval1;
460
461

      if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT)
462
        {
463
          total[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
          total[levelID].missval = missval1;
465
          total[levelID].resize(gridsize);
466
        }
467
      if (IS_SET(request->var1.f5))
468
        {
469
          var15[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470
          var15[levelID].missval = missval1;
471
          var15[levelID].resize(gridsize);
472
        }
473
      if (IS_SET(request->var2.h2))
474
        {
475
          var22[levelID].grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476
          var22[levelID].missval = missval1;
477
          var22[levelID].resize(gridsize);
478
479
        }
    }
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
  itsID = 0;
  while ((nrecs = cdoStreamInqTimestep(istreamID2, itsID)))
    {
      ivdate = taxisInqVdate(vlistInqTaxis(ivlistID2));

      const auto dayoy = decodeDayOfYear(ivdate);
      if (dayoy < 0 || dayoy >= MaxDays) cdoAbort("Day %d out of range!", dayoy);

      if (!vars2[dayoy].size())
        {
          fieldsFromVlist(ivlistID2, vars2[dayoy], FIELD_VEC);
        }

      for (int recID = 0; recID < nrecs; recID++)
        {
          size_t nmiss = 0;
          cdoInqRecord(istreamID2, &varID, &levelID);
          if (varID != FIRST_VAR_ID) continue;
498
          cdoReadRecord(istreamID2, vars2[dayoy][0][levelID].vec_d.data(), &nmiss);
499
500
501
502
          vars2[dayoy][0][levelID].nmiss = nmiss;
        }
      itsID++;
    }
503

504
505
  itsID = 0;
  otsID = 0;
506
  while (true)
507
508
    {
      nsets = 0;
509
      while ((nrecs = cdoStreamInqTimestep(istreamID1, itsID)) > 0)
510
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
          ivdate = taxisInqVdate(itaxisID1);
          ivtime = taxisInqVtime(itaxisID1);
513

514
515
516
517
          const auto dayoy = decodeDayOfYear(ivdate);
          if (!vars2[dayoy].size())
            cdoAbort("Input streams have different time values!");

518
519
520
521
522
          if (nsets == 0) 
            {
              SET_DATE(indate2, ivdate, ivtime);
              indate21 = ivdate;
            }
523
524
525
          SET_DATE(indate1, ivdate, ivtime);

          if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;
526

527
          for (int recID = 0; recID < nrecs; recID++)
528
            {
529
              cdoInqRecord(istreamID1, &varID, &levelID);
530

531
              if (varID != FIRST_VAR_ID) continue;
532

533
              if (nsets == 0)
534
                {
535
536
537
538
539
540
541
542
543
544
                  if ( cmplen != 0 && itsID == 0 && request->var1.f4 == &vfarnum2 )
                    {
                      fieldFill(var14[levelID], missval1);
                      var14[levelID].nmiss = gridsize;
                    }
                  else if ( request->var1.f4 != &vfarnum2 )
                    {
                      fieldFill(var14[levelID], missval1);
                      var14[levelID].nmiss = gridsize;
                    }
545
546
547
548
549
550
                  fieldFill(samp1[levelID], missval1);
                  fieldFill(samp2[levelID], missval1);
                  if (!samp3[levelID].empty()) fieldFill(samp3[levelID], 0.0);
                  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) fieldFill(total[levelID], 0.0);
                  if (IS_SET(request->var1.f5)) fieldFill(var15[levelID], missval1);
                  if (IS_SET(request->var2.h2)) fieldFill(var22[levelID], missval1);
551

Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
553
                  samp1[levelID].nmiss = gridsize;
                  samp2[levelID].nmiss = gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554
                  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) total[levelID].nmiss = gridsize;
555
556
                  if (IS_SET(request->var1.f5)) var15[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h2)) var22[levelID].nmiss = gridsize;
557
                }
558
              cdoReadRecord(istreamID1, field1.vec_d.data(), &field1.nmiss);
559
              field1.grid = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
              field1.missval = missval1;
561
              field2.grid = vars2[dayoy][0][levelID].grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
              field2.nmiss = vars2[dayoy][0][levelID].nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
              field2.missval = missval2;
564
565
              vfarcpy(field2, vars2[dayoy][0][levelID]); 

Uwe Schulzweida's avatar
Uwe Schulzweida committed
566

567
568
              vfarnum(samp1[levelID], field1);
              vfarnum(samp2[levelID], field2);
569

570
              if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) vfarsum(total[levelID], field1);
571

572
573
              if (IS_SET(request->var1.f1)) request->var1.f1(field1, request->var1.f1arg);
              if (IS_SET(request->var1.f2)) request->var1.f2(field2, request->var1.f2arg);
574

Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
              if (field1.nmiss || !samp3[levelID].empty())
576
                {
577
                  if (samp3[levelID].empty())
578
                    {
579
                      samp3[levelID].resize(gridsize);
580
                      fieldFill(samp3[levelID], nsets);
581
                    }
582
                  for (size_t i = 0; i < gridsize; i++)
583
                    {
584
585
                      if (DBL_IS_EQUAL(field1.vec_d[i], field1.missval)) continue;
                      samp3[levelID].vec_d[i]++;
586
587
                    }
                }
588

589
590
              if (IS_SET(request->var1.f5) && request->var1.f4 == &vfarnum2)
                {
591
                  varrayCopy(gridsize, var14[levelID].vec_d, field3.vec_d);
592
593
594
595
596
                  field3.nmiss = var14[levelID].nmiss;
                  field3.grid = var14[levelID].grid;
                  field3.missval = var14[levelID].missval;
                }

597
598
              request->var1.f3(field1, field2);
              request->var1.f4(var14[levelID], field1);
599

600
              if (IS_SET(request->var2.h2))
601
                {
602
                  varrayCopy(gridsize, var14[levelID].vec_d, field2.vec_d);
603
604
                  field2.nmiss = var14[levelID].nmiss;
                  field2.grid = var14[levelID].grid;
605
606
                  field2.missval = var14[levelID].missval;

607
                  if (IS_SET(request->var2.h1)) request->var2.h1(field2, request->var2.h1arg);
608

609
                  request->var2.h2(var22[levelID], field2);
610
611
                }

612
613
614
615
              if (IS_SET(request->var1.f5))
                {
                  if (cmplen != 0 && request->var1.f4 == &vfarnum2)
                    {
616
617
                      auto &array1 = field3.vec_d;
                      const auto &array2 = var14[levelID].vec_d;
618
                      const auto missvaltemp = field1.missval;
619
                
620
                      const auto len = field1.size;
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
                      if (len != field3.size) cdoAbort("Fields have different size (%s)", __func__);
                
                      if (field1.nmiss)
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], missvaltemp) && !DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      else
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      request->var1.f5(var15[levelID], field3, request->var1.f5arg);
                    }
                  else
                    request->var1.f5(var15[levelID], var14[levelID], request->var1.f5arg);
                }
644
645
            }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
646
647
          ovdate = ivdate;
          ovtime = ivtime;
648
649
650
651
          nsets++;
          itsID++;
        }

652
653
      if (nrecs == 0 && nsets == 0) break;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
      if (request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME)
655
        for (levelID = 0; levelID < nlevels; levelID++)
656
          {
657
            auto &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];
658

659
660
            if (samp3[levelID].empty())
              vfarcdiv(rvar, nsets);
661
            else
662
              vfardiv(rvar, samp3[levelID]);
663

664
            if (request->var1.epilog == PERCENT_OF_TIME) vfarcmul(rvar, 100.0);
665
          }
666
667
      else if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT)
        for (levelID = 0; levelID < nlevels; levelID++)
668
          {
669
670
            Field &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];

671
672
            vfardiv(rvar, total[levelID]);
            vfarcmul(rvar, 100.0);
673
674
          }

675
676
677
678
679
680
681
682
683
684
685
686
687
      if ( request->var1.refdate == 19550101 )
        {
          taxisDefVdate(otaxisID, ovdate);
          taxisDefVtime(otaxisID, ovtime);
          cdoDefTimestep(ostreamID, otsID);
        }
      else
        {
          int year, month, day;
          cdiDecodeDate(indate21, &year, &month, &day);
          defineMidOfTime(cdoOperatorF2(operatorID), otaxisID, year, month, 12);
          cdoDefTimestep(ostreamID, otsID);
        }
688

Uwe Schulzweida's avatar
Uwe Schulzweida committed
689
      if (otsID && vlistInqVarTimetype(ivlistID1, FIRST_VAR_ID) == TIME_CONSTANT) continue;
690
691

      varID = 0;
692
693
      for (levelID = 0; levelID < nlevels; levelID++)
        {
694
          Field &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];
695

696
697
          vfarsel(rvar, samp1[levelID]);
          vfarsel(rvar, samp2[levelID]);
698

699
          cdoDefRecord(ostreamID, varID, levelID);
700
          cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
701
702
703
704
705
706
        }
      if (IS_SET(request->var2.h2))
        {
          varID = 1;
          for (levelID = 0; levelID < nlevels; levelID++)
            {
707
              auto &rvar = var22[levelID];
708

709
710
              vfarsel(rvar, samp1[levelID]);
              vfarsel(rvar, samp2[levelID]);
711

712
              cdoDefRecord(ostreamID, varID, levelID);
713
              cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
714
            }
715
716
        }

717
      if (nrecs == 0) break;
718
      otsID++;
719
720
    }

721
722
723
  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID2);
  cdoStreamClose(istreamID1);
724
725
}

726
727
void
eca3(const ECA_REQUEST_3 *request)
728
{
729
  const auto operatorID = cdoOperatorID();
730

731
  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
732
  int64_t ivdate1 = 0, ivdate2 = 0, indate21 = 0;
733
  int ivtime1 = 0, ivtime2 = 0;
734
735
  int64_t ovdate = 0;
  int ovtime = 0;
736
  int nrecs;
737
  int varID, levelID;
738
739
740
  int itsID;
  int otsID;
  long nsets;
741

742
  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);
743

744
745
  const auto istreamID1 = cdoOpenRead(0);
  const auto istreamID2 = cdoOpenRead(1);
746

747
748
749
  const auto ivlistID1 = cdoStreamInqVlist(istreamID1);
  const auto ivlistID2 = cdoStreamInqVlist(istreamID2);
  const auto ovlistID = vlistCreate();
750

Uwe Schulzweida's avatar
Uwe Schulzweida committed
751
  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);
752

753
754
755
  const auto gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  const auto missval = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);
756

757
  varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);
758
759

  vlistDefVarMissval(ovlistID, varID, missval);
760

761
762
763
  if (IS_SET(request->name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->name);
  if (IS_SET(request->longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->longname);
  if (IS_SET(request->units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->units);
764

765
  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);
766

767
768
769
  const auto itaxisID1 = vlistInqTaxis(ivlistID1);
  const auto itaxisID2 = vlistInqTaxis(ivlistID2);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
770
771
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
772
773
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
774
  taxisDefRdate(otaxisID, request->refdate);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
775
  taxisDefRtime(otaxisID, 0);
776
777
  vlistDefTaxis(ovlistID, otaxisID);

778
  const auto ostreamID = cdoOpenWrite(2);
779
  cdoDefVlist(ostreamID, ovlistID);
780

781
  const auto gridsize = gridInqSize(gridID);
782

783
  Field field1, field2;
784
785
  field1.resize(gridsize);
  field2.resize(gridsize);
786

787
  const int nlevels = zaxisInqSize(zaxisID);
788

789
  FieldVector var1(nlevels), var2(nlevels);
790
791

  for (levelID = 0; levelID < nlevels; levelID++)
792
    {
793
      var1[levelID].grid = gridID;
794
      var1[levelID].missval = missval;
795
      var1[levelID].resize(gridsize);
796
797

      var2[levelID].grid = gridID;
798
      var2[levelID].missval = missval;
799
      var2[levelID].resize(gridsize);
800
801
    }

802
803
  itsID = 0;
  otsID = 0;
804
  while (true)
805
806
    {
      nsets = 0;
807
      while ((nrecs = cdoStreamInqTimestep(istreamID1, itsID)) > 0)
808
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
          if (!cdoStreamInqTimestep(istreamID2, itsID)) cdoAbort("Input streams have different number of time steps!");
810