ecacore.cc 51.6 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>
22
23
#include "cdo.h"
#include "cdo_int.h"
24
#include "grid.h"
25
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
#include "ecacore.h"
27
28
29
30
31
32
33
34
35
36
37
38
#include "ecautil.h"

#define FIRST_VAR_ID 0

#define IS_NOT_SET(x) (x == NULL)
#define IS_SET(x)     (x != NULL)


void eca1(const ECA_REQUEST_1 *request)
{
  const int operatorID = cdoOperatorID();
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
  int nmiss;
40
41
  int cmplen;
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
42
43
44
45
  int gridsize;
  int ivdate = 0, ivtime = 0;
  int ovdate = 0, ovtime = 0;
  int nrecs, nrecords;
46
  int gridID, zaxisID, varID, levelID;
47
48
49
50
51
52
53
54
55
  int itsID;
  int otsID;
  long nsets;
  int i;
  int istreamID, ostreamID;
  int ivlistID, ovlistID, itaxisID, otaxisID;
  int nlevels;
  int *recVarID, *recLevelID;
  double missval;
56
57
  field_type *var12 = NULL, *samp1 = NULL, *samp2 = NULL, *var13 = NULL, *var21 = NULL, *var23 = NULL, *var;
  field_type field1, field2;
58
  
59
  cmplen = DATE_LEN - cdoOperatorF2(operatorID);
60
61
62
63
64
65
66
67

  istreamID = streamOpenRead(cdoStreamName(0));

  ivlistID = streamInqVlist(istreamID);
  ovlistID = vlistCreate();
  
  gridID  = vlistInqVarGrid(ivlistID, FIRST_VAR_ID);
  zaxisID = vlistInqVarZaxis(ivlistID, FIRST_VAR_ID);
68
69
  missval = vlistInqVarMissval(ivlistID, FIRST_VAR_ID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
  varID   = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
71
72

  vlistDefVarMissval(ovlistID, varID, missval);
73
74
75
76
77
78
79
80
81
82
  
  if ( IS_SET(request->var1.name) )
    vlistDefVarName(ovlistID, varID, request->var1.name);
  if ( IS_SET(request->var1.longname) ) 
    vlistDefVarLongname(ovlistID, varID, request->var1.longname);
  if ( IS_SET(request->var1.units) ) 
    vlistDefVarUnits(ovlistID, varID, request->var1.units);

  if ( IS_SET(request->var2.h2) || IS_SET(request->var2.h3) )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
84
  
85
86
      vlistDefVarMissval(ovlistID, varID, missval);

87
88
89
90
      if ( IS_SET(request->var2.name) ) 
        vlistDefVarName(ovlistID, varID, request->var2.name);
      if ( IS_SET(request->var2.longname) ) 
        vlistDefVarLongname(ovlistID, varID, request->var2.longname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
      if ( IS_SET(request->var2.units) ) 
        vlistDefVarUnits(ovlistID, varID, request->var2.units);
93
94
    }
    
95
  if ( cdoOperatorF2(operatorID) == 16 ) vlistDefNtsteps(ovlistID, 1);
96
97

  itaxisID = vlistInqTaxis(ivlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
  otaxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
  taxisDefTunit(otaxisID, TUNIT_MINUTE);
100
101
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
  taxisDefRdate(otaxisID, 19550101);
  taxisDefRtime(otaxisID, 0);
104
105
106
107
108
109
110
  vlistDefTaxis(ovlistID, otaxisID);

  ostreamID = streamOpenWrite(cdoStreamName(1), cdoFiletype());

  streamDefVlist(ostreamID, ovlistID);

  nrecords   = vlistNrecs(ivlistID);
111
112
  recVarID   = (int*) Malloc(nrecords*sizeof(int));
  recLevelID = (int*) Malloc(nrecords*sizeof(int));
113
114
115

  gridsize = gridInqSize(gridID);

116
117
118
  field_init(&field1);
  field_init(&field2);

119
  field1.ptr = (double*) Malloc(gridsize*sizeof(double));
120
  if ( IS_SET(request->var2.h2) || IS_SET(request->var2.h3) ) 
121
    field2.ptr = (double*) Malloc(gridsize*sizeof(double));
122
123
124
125
126
  else
    field2.ptr = NULL;

  nlevels = zaxisInqSize(zaxisID);

127
128
129
  var12 = (field_type*) Malloc(nlevels*sizeof(field_type));
  samp1 = (field_type*) Malloc(nlevels*sizeof(field_type));
  samp2 = (field_type*) Malloc(nlevels*sizeof(field_type));
130
  if ( IS_SET(request->var1.f3) ) 
131
    var13 = (field_type*) Malloc(nlevels*sizeof(field_type));
132
133
    
  if ( IS_SET(request->var2.h2) ) 
134
    var21 = (field_type*) Malloc(nlevels*sizeof(field_type));
135
  if ( IS_SET(request->var2.h3) ) 
136
    var23 = (field_type*) Malloc(nlevels*sizeof(field_type));
137
138
139
      
  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
140
      field_init(&var12[levelID]);
141
142
143
      var12[levelID].grid    = gridID;
      var12[levelID].nmiss   = 0;
      var12[levelID].missval = missval;
144
      var12[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
145

146
      field_init(&samp1[levelID]);
147
148
149
      samp1[levelID].grid    = gridID;
      samp1[levelID].nmiss   = 0;
      samp1[levelID].missval = missval;
150
      samp1[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151

152
      field_init(&samp2[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
155
156
      samp2[levelID].grid    = gridID;
      samp2[levelID].nmiss   = 0;
      samp2[levelID].missval = missval;
      samp2[levelID].ptr     = NULL;
157
158
159

      if ( IS_SET(request->var1.f3) )
        {
160
	  field_init(&var13[levelID]);
161
162
163
          var13[levelID].grid    = gridID;
          var13[levelID].nmiss   = 0;
          var13[levelID].missval = missval;
164
          var13[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
165
166
167
        }
      if ( IS_SET(request->var2.h2) )
        {
168
	  field_init(&var21[levelID]);
169
170
171
          var21[levelID].grid    = gridID;
          var21[levelID].nmiss   = 0;
          var21[levelID].missval = missval;
172
          var21[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
173
174
175
        }
      if ( IS_SET(request->var2.h3) )
        {
176
	  field_init(&var23[levelID]);
177
178
179
          var23[levelID].grid    = gridID;
          var23[levelID].nmiss   = 0;
          var23[levelID].missval = missval;
180
          var23[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
181
182
183
184
185
186
187
188
189
190
191
192
193
        }
    }

  itsID   = 0;
  otsID   = 0;
  while ( TRUE )
    {
      nsets = 0;
      while ( (nrecs = streamInqTimestep(istreamID, itsID)) > 0 )
        {
          ivdate = taxisInqVdate(itaxisID);
          ivtime = taxisInqVtime(itaxisID);

194
195
	  if ( nsets == 0 ) SET_DATE(indate2, ivdate, ivtime);
	  SET_DATE(indate1, ivdate, ivtime);
196

197
	  if ( DATE_IS_NEQ(indate1, indate2, cmplen) ) break;
198

199
          for ( int recID = 0; recID < nrecs; recID++ )
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
            {
              streamInqRecord(istreamID, &varID, &levelID);

              if ( itsID == 0 )
                {
                  recVarID[recID]   = varID;
                  recLevelID[recID] = levelID;
                }
              if ( varID != FIRST_VAR_ID ) continue;

              if ( nsets == 0 )
                {
                  for ( i = 0; i < gridsize; i++ )
                    {
                      var12[levelID].ptr[i] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
                      samp1[levelID].ptr[i] = missval;
216
217
218
219
                      if ( IS_SET(samp2[levelID].ptr) ) samp2[levelID].ptr[i] = 0.0;
                      if ( IS_SET(request->var1.f3) ) var13[levelID].ptr[i] = missval;
                      if ( IS_SET(request->var2.h2) ) var21[levelID].ptr[i] = missval;
                      if ( IS_SET(request->var2.h3) ) var23[levelID].ptr[i] = missval;
220
221
                    }
                  var12[levelID].nmiss = gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
                  samp1[levelID].nmiss = gridsize;
223
224
225
                  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; 
226
227
                }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
              streamReadRecord(istreamID, field1.ptr, &nmiss);
              field1.nmiss   = (size_t)nmiss;
230
231
              field1.grid    = var12[levelID].grid;
              field1.missval = var12[levelID].missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
233
              
              farnum(&samp1[levelID], field1);
234
235
236
237
238
239
240
241
242
243
244

              if ( IS_SET(request->var2.h2) )
                {
                  memcpy(field2.ptr, field1.ptr, gridsize*sizeof(double));
                  field2.nmiss   = field1.nmiss;
                  field2.grid    = field1.grid;
                  field2.missval = field1.missval;
                }
		
              if ( IS_SET(request->var1.f1) ) 
                request->var1.f1(&field1, request->var1.f1arg);
245

Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
              if ( field1.nmiss > 0 || IS_SET(samp2[levelID].ptr) )
247
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
                  if ( IS_NOT_SET(samp2[levelID].ptr) )
249
                    {
250
                      samp2[levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
251
                      for ( i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
                        samp2[levelID].ptr[i] = nsets;
253
254
255
256
257
                    }
                  for ( i = 0; i < gridsize; i++ )
                    {
                      if ( DBL_IS_EQUAL(field1.ptr[i], field1.missval) )
                        continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
                      samp2[levelID].ptr[i]++;
259
260
261
                    }
                }
              
262
              if ( IS_NOT_EQUAL(request->var1.mulc, 0.0) )
263
                farcmul(&field1, request->var1.mulc);
264
              if ( IS_NOT_EQUAL(request->var1.addc, 0.0) )
265
266
267
268
269
270
271
272
273
274
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
300
301
302
303
304
305
306
307
308
309
310
311
312
                farcadd(&field1, request->var1.addc);
                
              request->var1.f2(&var12[levelID], field1);
              
              if ( IS_SET(request->var2.h2) || IS_SET(request->var2.h3) )
                { 
                  /* if h2 is null, use the output of f2 as input for h1 */
                  if ( IS_NOT_SET(request->var2.h2) )
                    {
                      memcpy(field2.ptr, var12[levelID].ptr, gridsize*sizeof(double));
                      field2.nmiss   = var12[levelID].nmiss;
                      field2.grid    = var12[levelID].grid;
                      field2.missval = var12[levelID].missval;
                    }
                  
                  if ( IS_SET(request->var2.h1) ) 
                    request->var2.h1(&field2, request->var2.h1arg);
                        
                  if ( IS_NOT_SET(request->var2.h2) )    
                    request->var2.h3(&var23[levelID], field2);
                  else
                    {
                      request->var2.h2(&var21[levelID], field2);
                      if ( IS_SET(request->var2.h3) )
                        request->var2.h3(&var23[levelID], var21[levelID]);
                    }
                }

              if ( IS_SET(request->var1.f3) )
                request->var1.f3(&var13[levelID], var12[levelID]);
            }

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

      if ( nrecs == 0 && nsets == 0 ) break;
      
      if ( request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME )
        for ( levelID = 0; levelID < nlevels; levelID++ )
          {
            if ( IS_SET(request->var1.f3) )
              var = &var13[levelID];
            else 
              var = &var12[levelID];
              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
            if ( IS_NOT_SET(samp2[levelID].ptr) )
314
315
              farcdiv(var, nsets);
            else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
              fardiv(var, samp2[levelID]);
317
318
319
320
321
322
323
              
            if ( request->var1.epilog == PERCENT_OF_TIME )
              farcmul(var, 100.0);
          }

      taxisDefVdate(otaxisID, ovdate);
      taxisDefVtime(otaxisID, ovtime);
324
      streamDefTimestep(ostreamID, otsID);
325

Uwe Schulzweida's avatar
Uwe Schulzweida committed
326
      if ( otsID && vlistInqVarTsteptype(ivlistID, FIRST_VAR_ID) == TSTEP_CONSTANT ) continue;
327
328
329
330
331
332
333
334

      varID = 0;
      for ( levelID = 0; levelID < nlevels; levelID++ )
	{
	  if ( IS_SET(request->var1.f3) )
	    var = &var13[levelID];
	  else 
	    var = &var12[levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335

336
	  farsel(var, samp1[levelID]);
337

338
339
340
	  streamDefRecord(ostreamID, varID, levelID);
	  streamWriteRecord(ostreamID, var->ptr, var->nmiss);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
341

342
343
344
345
346
347
348
349
350
      if ( IS_SET(request->var2.h2) || IS_SET(request->var2.h3) )
	{
	  varID = 1;
	  for ( levelID = 0; levelID < nlevels; levelID++ )
	    {
	      if ( IS_SET(request->var2.h3) )
		var = &var23[levelID];
	      else 
		var = &var21[levelID];
351
                  
352
	      farsel(var, samp1[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
                  
354
355
356
	      streamDefRecord(ostreamID, varID, levelID);
	      streamWriteRecord(ostreamID, var->ptr, var->nmiss);
	    }
357
358
359
        }

      if ( nrecs == 0 ) break;
360
      otsID++;
361
362
363
364
    }

  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
365
366
367
      Free(var12[levelID].ptr);
      Free(samp1[levelID].ptr);
      if ( IS_SET(samp2[levelID].ptr) ) Free(samp2[levelID].ptr);
368
    }  
369
370
371
  Free(var12);
  Free(samp1);
  Free(samp2);
372
373
374
375
  
  if ( IS_SET(var13) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
376
377
        Free(var13[levelID].ptr);
      Free(var13);
378
379
380
381
    }
  if ( IS_SET(var21) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
382
383
        Free(var21[levelID].ptr);
      Free(var21);
384
385
386
387
    }
  if ( IS_SET(var23) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
388
389
        Free(var23[levelID].ptr);
      Free(var23);
390
391
    }
    
392
393
  if ( IS_SET(field1.ptr) ) Free(field1.ptr);
  if ( IS_SET(field2.ptr) ) Free(field2.ptr);
394

395
396
  if ( IS_SET(recVarID) )   Free(recVarID);
  if ( IS_SET(recLevelID) ) Free(recLevelID);
397
398
399
400
401
402
403
404
405
406

  streamClose(ostreamID);
  streamClose(istreamID);
}


void eca2(const ECA_REQUEST_2 *request)
{
  const int operatorID = cdoOperatorID();
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407
  int nmiss;
408
409
  int cmplen;
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
410
  int gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
  int ivdate = 0, ivtime = 0;
412
413
  int ovdate = 0, ovtime = 0;
  int nrecs, nrecords;
414
  int gridID, zaxisID, varID, levelID;
415
416
417
418
419
  int itsID;
  int otsID;
  long nsets;
  int i;
  int istreamID1, istreamID2, ostreamID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
  int ivlistID1, ivlistID2, ovlistID, itaxisID1, otaxisID;
421
422
  int nlevels;
  int *recVarID, *recLevelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423
  double missval1, missval2;
424
425
  field_type *var14 = NULL, *samp1 = NULL, *samp2 = NULL, *samp3 = NULL, *total = NULL, *var15 = NULL, *var22 = NULL, *var;
  field_type field1, field2;
426
  
427
  cmplen = DATE_LEN - cdoOperatorF2(operatorID);
428
429
430
431
432
433
434
435

  istreamID1 = streamOpenRead(cdoStreamName(0));
  istreamID2 = streamOpenRead(cdoStreamName(1));

  ivlistID1 = streamInqVlist(istreamID1);
  ivlistID2 = streamInqVlist(istreamID2);
  ovlistID  = vlistCreate();
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);
437
438
439
  
  gridID  = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
440
441
442
  missval1 = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);
  missval2 = vlistInqVarMissval(ivlistID2, FIRST_VAR_ID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
443
  varID   = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
444
  
445
446
  vlistDefVarMissval(ovlistID, varID, missval1);

447
448
449
450
451
452
453
454
455
  if ( IS_SET(request->var1.name) ) 
    vlistDefVarName(ovlistID, varID, request->var1.name);
  if ( IS_SET(request->var1.longname) ) 
    vlistDefVarLongname(ovlistID, varID, request->var1.longname);
  if ( IS_SET(request->var1.units) ) 
    vlistDefVarUnits(ovlistID, varID, request->var1.units);

  if ( IS_SET(request->var2.h2) )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
457
  
458
459
      vlistDefVarMissval(ovlistID, varID, missval1);

460
461
462
463
      if ( IS_SET(request->var2.name) ) 
        vlistDefVarName(ovlistID, varID, request->var2.name);
      if ( IS_SET(request->var2.longname) ) 
        vlistDefVarLongname(ovlistID, varID, request->var2.longname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
465
      if ( IS_SET(request->var2.units) ) 
        vlistDefVarUnits(ovlistID, varID, request->var2.units);
466
467
    }

468
  if ( cdoOperatorF2(operatorID) == 16 ) vlistDefNtsteps(ovlistID, 1);
469
470

  itaxisID1 = vlistInqTaxis(ivlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
  otaxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
  taxisDefTunit(otaxisID, TUNIT_MINUTE);
473
474
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
476
  taxisDefRdate(otaxisID, 19550101);
  taxisDefRtime(otaxisID, 0);
477
478
479
480
481
482
483
  vlistDefTaxis(ovlistID, otaxisID);

  ostreamID = streamOpenWrite(cdoStreamName(2), cdoFiletype());

  streamDefVlist(ostreamID, ovlistID);

  nrecords   = vlistNrecs(ivlistID1);
484
485
  recVarID   = (int*) Malloc(nrecords*sizeof(int));
  recLevelID = (int*) Malloc(nrecords*sizeof(int));
486
487
488

  gridsize = gridInqSize(gridID);

489
490
  field_init(&field1);
  field_init(&field2);
491
492
  field1.ptr = (double*) Malloc(gridsize*sizeof(double));
  field2.ptr = (double*) Malloc(gridsize*sizeof(double));
493
494

  nlevels = zaxisInqSize(zaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
  
496
497
498
499
  var14 = (field_type*) Malloc(nlevels*sizeof(field_type));
  samp1 = (field_type*) Malloc(nlevels*sizeof(field_type));
  samp2 = (field_type*) Malloc(nlevels*sizeof(field_type));
  samp3 = (field_type*) Malloc(nlevels*sizeof(field_type));
500
501
  
  if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
502
    total = (field_type*) Malloc(nlevels*sizeof(field_type));
503
  if ( IS_SET(request->var1.f5) )
504
    var15 = (field_type*) Malloc(nlevels*sizeof(field_type));
505
  if ( IS_SET(request->var2.h2) )
506
    var22 = (field_type*) Malloc(nlevels*sizeof(field_type));
507
508
509
      
  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
510
      field_init(&var14[levelID]);
511
512
      var14[levelID].grid    = gridID;
      var14[levelID].nmiss   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513
      var14[levelID].missval = missval1;
514
      var14[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
515
      
516
      field_init(&samp1[levelID]);
517
518
      samp1[levelID].grid    = gridID;
      samp1[levelID].nmiss   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
      samp1[levelID].missval = missval1;
520
      samp1[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
521

522
      field_init(&samp2[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
524
525
      samp2[levelID].grid    = gridID;
      samp2[levelID].nmiss   = 0;
      samp2[levelID].missval = missval1;
526
      samp2[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527

528
      field_init(&samp3[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
530
531
532
      samp3[levelID].grid    = gridID;
      samp3[levelID].nmiss   = 0;
      samp3[levelID].missval = missval1;
      samp3[levelID].ptr     = NULL;
533
534
535
      
      if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
        {
536
	  field_init(&total[levelID]);
537
538
          total[levelID].grid    = gridID;
          total[levelID].nmiss   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
539
          total[levelID].missval = missval1;
540
          total[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
541
542
543
        }
      if ( IS_SET(request->var1.f5) )
        {
544
	  field_init(&var15[levelID]);
545
546
          var15[levelID].grid    = gridID;
          var15[levelID].nmiss   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547
          var15[levelID].missval = missval1;
548
          var15[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
549
550
551
        }
      if ( IS_SET(request->var2.h2) )
        {
552
	  field_init(&var22[levelID]);
553
554
          var22[levelID].grid    = gridID;
          var22[levelID].nmiss   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555
          var22[levelID].missval = missval1;
556
          var22[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
557
558
559
560
561
562
563
564
565
566
567
568
        }
    }

  itsID   = 0;
  otsID   = 0;
  while ( TRUE )
    {
      nsets = 0;
      while ( (nrecs = streamInqTimestep(istreamID1, itsID)) > 0 )
        {
	  if ( !streamInqTimestep(istreamID2, itsID) )
	    cdoAbort("Input streams have different number of time steps!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
569
570
571
572

          ivdate = taxisInqVdate(itaxisID1);
          ivtime = taxisInqVtime(itaxisID1);
          
573
574
	  if ( nsets == 0 ) SET_DATE(indate2, ivdate, ivtime);
	  SET_DATE(indate1, ivdate, ivtime);
575

576
	  if ( DATE_IS_NEQ(indate1, indate2, cmplen) ) break;
577

578
          for ( int recID = 0; recID < nrecs; recID++ )
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
            {
              streamInqRecord(istreamID1, &varID, &levelID);
              streamInqRecord(istreamID2, &varID, &levelID);

              if ( itsID == 0 )
                {
                  recVarID[recID]   = varID;
                  recLevelID[recID] = levelID;
                }
              if ( varID != FIRST_VAR_ID ) continue;

              if ( nsets == 0 )
                {
                  for ( i = 0; i < gridsize; i++ )
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
594
595
596
597
598
                      var14[levelID].ptr[i] = missval1;
                      samp1[levelID].ptr[i] = missval1;
                      samp2[levelID].ptr[i] = missval1;
                      if ( IS_SET(samp3[levelID].ptr) )
                        samp3[levelID].ptr[i] = 0.0;
599
600
601
                      if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
                        total[levelID].ptr[i] = 0.0;
                      if ( IS_SET(request->var1.f5) ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
                        var15[levelID].ptr[i] = missval1;
603
                      if ( IS_SET(request->var2.h2) ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604
                        var22[levelID].ptr[i] = missval1;
605
606
                    }
                  var14[levelID].nmiss = gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
                  samp1[levelID].nmiss = gridsize;
                  samp2[levelID].nmiss = gridsize;
609
610
611
612
613
614
615
616
                  if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
                    total[levelID].nmiss = gridsize;
                  if ( IS_SET(request->var1.f5) )
                    var15[levelID].nmiss = gridsize;
                  if ( IS_SET(request->var2.h2) )
                    var22[levelID].nmiss = gridsize;
                }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
618
              streamReadRecord(istreamID1, field1.ptr, &nmiss);
              field1.nmiss   = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
620
621
              field1.grid    = gridID;
              field1.missval = missval1;
              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
623
              streamReadRecord(istreamID2, field2.ptr, &nmiss);
              field2.nmiss   = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
624
625
626
627
628
              field2.grid    = gridID;
              field2.missval = missval2;

              farnum(&samp1[levelID], field1);
              farnum(&samp2[levelID], field2);
629
630
631
632
633
634
635
636
637
638

              if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
                farsum(&total[levelID], field1);
                             
              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);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
639
              if ( field1.nmiss > 0 || IS_SET(samp3[levelID].ptr) )
640
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
641
                  if ( IS_NOT_SET(samp3[levelID].ptr) )
642
                    {
643
                      samp3[levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
644
                      for ( i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
645
                        samp3[levelID].ptr[i] = nsets;
646
647
648
649
650
                    }
                  for ( i = 0; i < gridsize; i++ )
                    {
                      if ( DBL_IS_EQUAL(field1.ptr[i], field1.missval) )
                        continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
651
                      samp3[levelID].ptr[i]++;
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
                    }
                }
              
              request->var1.f3(&field1, field2);
              request->var1.f4(&var14[levelID], field1);

              if ( IS_SET(request->var2.h2) )
                {
                  memcpy(field2.ptr, var14[levelID].ptr, gridsize*sizeof(double));
                  field2.nmiss   = var14[levelID].nmiss;
                  field2.grid    = var14[levelID].grid;
                  field2.missval = var14[levelID].missval;

                  if ( IS_SET(request->var2.h1) )
                    request->var2.h1(&field2, request->var2.h1arg);
                    
                  request->var2.h2(&var22[levelID], field2);
                }

              if ( IS_SET(request->var1.f5) )
                request->var1.f5(&var15[levelID], var14[levelID], request->var1.f5arg);              
            }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
675
676
          ovdate = ivdate;
          ovtime = ivtime;
677
678
679
680
681
682
683
684
685
686
687
688
689
690
          nsets++;
          itsID++;
        }

      if ( nrecs == 0 && nsets == 0 ) break;
      
      if ( request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME )
        for ( levelID = 0; levelID < nlevels; levelID++ )
          {
	    if ( IS_SET(request->var1.f5) )
	      var = &var15[levelID];
	    else
	      var = &var14[levelID];
	  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
691
            if ( IS_NOT_SET(samp3[levelID].ptr) )
692
693
              farcdiv(var, nsets);
            else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
              fardiv(var, samp3[levelID]);
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712

	    if ( request->var1.epilog == PERCENT_OF_TIME )
	      farcmul(var, 100.0);
          }
      else if ( request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT )
        for ( levelID = 0; levelID < nlevels; levelID++ )
          {
	    if ( IS_SET(request->var1.f5) )
	      var = &var15[levelID];
	    else
	      var = &var14[levelID];
	  
            fardiv(var, total[levelID]);
	    farcmul(var, 100.0);
          }

      taxisDefVdate(otaxisID, ovdate);
      taxisDefVtime(otaxisID, ovtime);
713
      streamDefTimestep(ostreamID, otsID);
714

Uwe Schulzweida's avatar
Uwe Schulzweida committed
715
      if ( otsID && vlistInqVarTsteptype(ivlistID1, FIRST_VAR_ID) == TSTEP_CONSTANT ) continue;
716
717
718
719
720
721
722
723

      varID = 0;
      for ( levelID = 0; levelID < nlevels; levelID++ )
	{
	  if ( IS_SET(request->var1.f5) )
	    var = &var15[levelID];
	  else 
	    var = &var14[levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
724
              
725
726
	  farsel(var, samp1[levelID]);
	  farsel(var, samp2[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
              
728
729
730
731
732
733
734
735
736
737
738
739
	  streamDefRecord(ostreamID, varID, levelID);
	  streamWriteRecord(ostreamID, var->ptr, var->nmiss);
	}
      if ( IS_SET(request->var2.h2) )
	{
	  varID = 1;
	  for ( levelID = 0; levelID < nlevels; levelID++ )
	    {
	      var = &var22[levelID];
              
	      farsel(var, samp1[levelID]);
	      farsel(var, samp2[levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
740
                  
741
742
743
	      streamDefRecord(ostreamID, varID, levelID);
	      streamWriteRecord(ostreamID, var->ptr, var->nmiss);
	    }
744
745
746
        }

      if ( nrecs == 0 ) break;
747
      otsID++;
748
749
750
751
    }

  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
752
753
754
755
      Free(var14[levelID].ptr);
      Free(samp1[levelID].ptr);
      Free(samp2[levelID].ptr);
      if ( IS_SET(samp3[levelID].ptr) ) Free(samp3[levelID].ptr);
756
    }
757
758
759
760
  Free(var14);
  Free(samp1);
  Free(samp2);
  Free(samp3);
761
762
763
764
  
  if ( IS_SET(total) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
765
766
        Free(total[levelID].ptr);
      Free(total);
767
768
769
770
    }
  if ( IS_SET(var15) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
771
772
        Free(var15[levelID].ptr);
      Free(var15);
773
774
775
776
    }
  if ( IS_SET(var22) ) 
    {
      for ( levelID = 0; levelID < nlevels; levelID++ )
777
778
        Free(var22[levelID].ptr);
      Free(var22);
779
780
    }
  
781
782
  if ( IS_SET(field1.ptr) ) Free(field1.ptr);
  if ( IS_SET(field2.ptr) ) Free(field2.ptr);
783

784
785
  if ( IS_SET(recVarID) )   Free(recVarID);
  if ( IS_SET(recLevelID) ) Free(recLevelID);
786
787
788
789
790
791
792
793
794
795

  streamClose(ostreamID);
  streamClose(istreamID2);
  streamClose(istreamID1);
}


void eca3(const ECA_REQUEST_3 *request)
{
  const int operatorID = cdoOperatorID();
796

Uwe Schulzweida's avatar
Uwe Schulzweida committed
797
  int nmiss;
798
799
  int cmplen;
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
800
801
802
803
804
  int gridsize;
  int ivdate1 = 0, ivtime1 = 0;
  int ivdate2 = 0, ivtime2 = 0;
  int ovdate = 0, ovtime = 0;
  int nrecs, nrecords;
805
  int gridID, zaxisID, varID, levelID;
806
807
808
809
810
811
812
813
814
  int itsID;
  int otsID;
  long nsets;
  int i;
  int istreamID1, istreamID2, ostreamID;
  int ivlistID1, ivlistID2, ovlistID, itaxisID1, itaxisID2, otaxisID;
  int nlevels;
  int *recVarID, *recLevelID;
  double missval;
815
816
  field_type *var1 = NULL, *var2 = NULL;
  field_type field1, field2;
817
  
818
  cmplen = DATE_LEN - cdoOperatorF2(operatorID);
819
820
821
822
823
824
825
826

  istreamID1 = streamOpenRead(cdoStreamName(0));
  istreamID2 = streamOpenRead(cdoStreamName(1));

  ivlistID1 = streamInqVlist(istreamID1);
  ivlistID2 = streamInqVlist(istreamID2);
  ovlistID  = vlistCreate();
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827
  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);
828
829
830
  
  gridID  = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
831
832
  missval = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
  varID   = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
834
835

  vlistDefVarMissval(ovlistID, varID, missval);
836
837
838
839
840
841
842
843
  
  if ( IS_SET(request->name) ) 
    vlistDefVarName(ovlistID, varID, request->name);
  if ( IS_SET(request->longname) ) 
    vlistDefVarLongname(ovlistID, varID, request->longname);
  if ( IS_SET(request->units) ) 
    vlistDefVarUnits(ovlistID, varID, request->units);

844
  if ( cdoOperatorF2(operatorID) == 16 ) vlistDefNtsteps(ovlistID, 1);
845
846
847

  itaxisID1 = vlistInqTaxis(ivlistID1);
  itaxisID2 = vlistInqTaxis(ivlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
848
  otaxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
849
  taxisDefTunit(otaxisID, TUNIT_MINUTE);
850
851
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
852
853
  taxisDefRdate(otaxisID, 19550101);
  taxisDefRtime(otaxisID, 0);
854
855
856
857
858
859
860
  vlistDefTaxis(ovlistID, otaxisID);

  ostreamID = streamOpenWrite(cdoStreamName(2), cdoFiletype());

  streamDefVlist(ostreamID, ovlistID);

  nrecords   = vlistNrecs(ivlistID1);
861
862
  recVarID   = (int*) Malloc(nrecords*sizeof(int));
  recLevelID = (int*) Malloc(nrecords*sizeof(int));
863
864
865

  gridsize = gridInqSize(gridID);

866
867
  field_init(&field1);
  field_init(&field2);
868
869
  field1.ptr = (double*) Malloc(gridsize*sizeof(double));
  field2.ptr = (double*) Malloc(gridsize*sizeof(double));
870
871
872

  nlevels = zaxisInqSize(zaxisID);

873
874
  var1 = (field_type*) Malloc(nlevels*sizeof(field_type));
  var2 = (field_type*) Malloc(nlevels*sizeof(field_type));
875
876
877
        
  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
878
      field_init(&var1[levelID]);
879
880
881
      var1[levelID].grid    = gridID;
      var1[levelID].nmiss   = 0;
      var1[levelID].missval = missval;
882
      var1[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
883
            
884
      field_init(&var2[levelID]);
885
886
887
      var2[levelID].grid    = gridID;
      var2[levelID].nmiss   = 0;
      var2[levelID].missval = missval;
888
      var2[levelID].ptr     = (double*) Malloc(gridsize*sizeof(double));
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
    }

  itsID   = 0;
  otsID   = 0;
  while ( TRUE )
    {
      nsets = 0;
      while ( (nrecs = streamInqTimestep(istreamID1, itsID)) > 0 )
        {
	  if ( !streamInqTimestep(istreamID2, itsID) )
	    cdoAbort("Input streams have different number of time steps!");
            
          ivdate1 = taxisInqVdate(itaxisID1);
          ivdate2 = taxisInqVdate(itaxisID2);
          if ( ivdate1 != ivdate2 )
904
            cdoAbort("Input streams have different verification dates at time step %d!", itsID+1);
905
906
907
908
            
          ivtime1 = taxisInqVtime(itaxisID1);
          ivtime2 = taxisInqVtime(itaxisID2);
          if ( ivtime1 != ivtime2 )
909
            cdoAbort("Input streams have different verification times at time step %d!", itsID+1);
910
911
912
   
	  if ( nsets == 0 ) SET_DATE(indate2, ivdate1, ivtime1);
	  SET_DATE(indate1, ivdate1, ivtime1);
913

914
	  if ( DATE_IS_NEQ(indate1, indate2, cmplen) ) break;
915

916
          for ( int recID = 0; recID < nrecs; recID++ )
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
            {
              streamInqRecord(istreamID1, &varID, &levelID);
              streamInqRecord(istreamID2, &varID, &levelID);

              if ( itsID == 0 )
                {
                  recVarID[recID]   = varID;
                  recLevelID[recID] = levelID;
                }
              if ( varID != FIRST_VAR_ID ) continue;

              if ( nsets == 0 )
                {
                  for ( i = 0; i < gridsize; i++ )
                    {
                      var1[levelID].ptr[i] = missval;
                      var2[levelID].ptr[i] = missval;
                    }
                  var1[levelID].nmiss = gridsize;
                  var2[levelID].nmiss = gridsize;
                }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
939
940
              streamReadRecord(istreamID1, field1.ptr, &nmiss);
              field1.nmiss   = (size_t)nmiss;
941
942
943
              field1.grid    = var1[levelID].grid;
              field1.missval = var1[levelID].missval;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
944
945
              streamReadRecord(istreamID2, field2.ptr, &nmiss);
              field2.nmiss   = (size_t)nmiss;
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
              field2.grid    = var1[levelID].grid;
              field2.missval = var1[levelID].missval;

              request->f1(&var1[levelID], field1);
              request->f2(&var2[levelID], field2);
            }

          ovdate = ivdate1;
          ovtime = ivtime1;
          nsets++;
          itsID++;
        }

      if ( nrecs == 0 && nsets == 0 ) break;
      
      for ( levelID = 0; levelID < nlevels; levelID++ )
        request->f3(&var1[levelID], var2[levelID]);

      taxisDefVdate(otaxisID, ovdate);
      taxisDefVtime(otaxisID, ovtime);
966
      streamDefTimestep(ostreamID, otsID);
967

Uwe Schulzweida's avatar
Uwe Schulzweida committed
968
      if ( otsID && vlistInqVarTsteptype(ivlistID1, FIRST_VAR_ID) == TSTEP_CONSTANT ) continue;
969
970
971
972
973
974
975

      varID = 0;
      for ( levelID = 0; levelID < nlevels; levelID++ )
	{
	  streamDefRecord(ostreamID, varID, levelID);
	  streamWriteRecord(ostreamID, var1[levelID].ptr, var1[levelID].nmiss);
	}
976
977

      if ( nrecs == 0 ) break;
978
      otsID++;
979
980
981
982
    }

  for ( levelID = 0; levelID < nlevels; levelID++ )
    {
983
984
      Free(var1[levelID].ptr);
      Free(var2[levelID].ptr);
985
    }
986
987
  Free(var1);
  Free(var2);
988
    
989
990
  if ( IS_SET(field1.ptr) ) Free(field1.ptr);
  if ( IS_SET(field2.ptr) ) Free(field2.ptr);
991

992
993
  if ( IS_SET(recVarID) )   Free(recVarID);
  if ( IS_SET(recLevelID) ) Free(recLevelID);
994
995
996
997
998
999
1000
1001
1002
1003

  streamClose(ostreamID);
  streamClose(istreamID2);
  streamClose(istreamID1);
}


void eca4(const ECA_REQUEST_4 *request)
{
  const int operatorID = cdoOperatorID();
1004

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1005
  int nmiss;
1006
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
1007
1008
  int ivdate = 0, ivtime = 0;
  int ovdate = 0, ovtime = 0;
1009
1010
1011
  int yearcnt = 0;
  int nrecs;
  int varID, levelID;
1012
  field_type *startDateWithHist[2], *endDateWithHist[2];
1013
  int resetAtJan = FALSE, resetAtJul = FALSE;
1014
1015
  int isFirstYear = TRUE;

1016
  int cmplen = DATE_LEN - cdoOperatorF2(operatorID);
1017

1018
1019
  int istreamID1 = streamOpenRead(cdoStreamName(0));
  int istreamID2 = streamOpenRead(cdoStreamName(1));
1020

1021
1022
1023
  int ivlistID1 = streamInqVlist(istreamID1);
  int ivlistID2 = streamInqVlist(istreamID2);
  int ovlistID = vlistCreate();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024

1025
  int gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1026
  if ( gridID != vlistInqVarGrid(ivlistID2, FIRST_VAR_ID) ) cdoAbort("Grid sizes of the input fields do not match!");
1027

1028
1029
  int zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  double missval = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);
1030

1031
  int ovarID1 = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
1032

1033
1034
  vlistDefVarMissval(ovlistID, ovarID1, missval);

1035
1036
1037
  if ( IS_SET(request->name) )      vlistDefVarName(ovlistID, ovarID1, request->name);
  if ( IS_SET(request->longname) )  vlistDefVarLongname(ovlistID, ovarID1, request->longname);
  if ( IS_SET(request->units) )     vlistDefVarUnits(ovlistID, ovarID1, request->units);
1038

1039
  int ovarID2 = vlistDefVar(ovlistID, gridID, zaxisID, TSTEP_INSTANT);
1040
  
1041
1042
  vlistDefVarMissval(ovlistID, ovarID2, missval);

1043
1044
1045
  if ( IS_SET(request->name2) )      vlistDefVarName(ovlistID, ovarID2, request->name2);
  if ( IS_SET(request->longname2) )  vlistDefVarLongname(ovlistID, ovarID2, request->longname2);
  if ( IS_SET(request->units2) )     vlistDefVarUnits(ovlistID, ovarID2, request->units2);
1046

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

1049
1050
  int itaxisID = vlistInqTaxis(ivlistID1);
  int otaxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1051
  taxisDefTunit(otaxisID, TUNIT_MINUTE);
1052
1053
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1054
1055
  taxisDefRdate(otaxisID, 19550101);
  taxisDefRtime(otaxisID, 0);
1056
1057
  vlistDefTaxis(ovlistID, otaxisID);

1058
  int ostreamID = streamOpenWrite(cdoStreamName(2), cdoFiletype());
1059
1060
1061

  streamDefVlist(ostreamID, ovlistID);

1062
1063
1064
  int nrecords    = vlistNrecs(ivlistID1);
  int *recVarID    = (int*) Malloc(nrecords*sizeof(int));
  int *recLevelID  = (int*) Malloc(nrecords*sizeof(int));
1065

1066
  int