Timstat.cc 17.5 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

18
19
20
/*
   This module contains the following operators:

21
      Timstat    timrange        Time range
22
23
24
25
26
      Timstat    timmin          Time minimum
      Timstat    timmax          Time maximum
      Timstat    timsum          Time sum
      Timstat    timmean         Time mean
      Timstat    timavg          Time average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
      Timstat    timvar          Time variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
      Timstat    timvar1         Time variance [Normalize by (n-1)]
29
      Timstat    timstd          Time standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
      Timstat    timstd1         Time standard deviation [Normalize by (n-1)]
31
      Hourstat   hourrange       Hourly range
32
33
34
35
36
      Hourstat   hourmin         Hourly minimum
      Hourstat   hourmax         Hourly maximum
      Hourstat   hoursum         Hourly sum
      Hourstat   hourmean        Hourly mean
      Hourstat   houravg         Hourly average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
      Hourstat   hourvar         Hourly variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
      Hourstat   hourvar1        Hourly variance [Normalize by (n-1)]
39
      Hourstat   hourstd         Hourly standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
      Hourstat   hourstd1        Hourly standard deviation [Normalize by (n-1)]
41
      Daystat    dayrange        Daily range
42
43
44
45
46
      Daystat    daymin          Daily minimum
      Daystat    daymax          Daily maximum
      Daystat    daysum          Daily sum
      Daystat    daymean         Daily mean
      Daystat    dayavg          Daily average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
      Daystat    dayvar          Daily variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
      Daystat    dayvar1         Daily variance [Normalize by (n-1)]
49
      Daystat    daystd          Daily standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
      Daystat    daystd1         Daily standard deviation [Normalize by (n-1)]
51
      Monstat    monrange        Monthly range
52
53
54
55
56
      Monstat    monmin          Monthly minimum
      Monstat    monmax          Monthly maximum
      Monstat    monsum          Monthly sum
      Monstat    monmean         Monthly mean
      Monstat    monavg          Monthly average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
      Monstat    monvar          Monthly variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
      Monstat    monvar1         Monthly variance [Normalize by (n-1)]
59
      Monstat    monstd          Monthly standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
      Monstat    monstd1         Monthly standard deviation [Normalize by (n-1)]
61
      Yearstat   yearrange       Yearly range
62
63
64
65
66
      Yearstat   yearmin         Yearly minimum
      Yearstat   yearmax         Yearly maximum
      Yearstat   yearsum         Yearly sum
      Yearstat   yearmean        Yearly mean
      Yearstat   yearavg         Yearly average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
      Yearstat   yearvar         Yearly variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
      Yearstat   yearvar1        Yearly variance [Normalize by (n-1)]
69
      Yearstat   yearstd         Yearly standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
      Yearstat   yearstd1        Yearly standard deviation [Normalize by (n-1)]
71
72
73
*/


Ralf Mueller's avatar
Ralf Mueller committed
74
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
78
79
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
80
81
82
83
84
85
typedef struct {
  short varID;
  short levelID;
} recinfo_t;


Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
void *Timstat(void *argument)
{
88
  enum {HOUR_LEN=4, DAY_LEN=6, MON_LEN=8, YEAR_LEN=10};
89
  int timestat_date = TIMESTAT_MEAN;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
  int vdate0 = 0, vtime0 = 0;
91
  int nrecs;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
92
  int varID, levelID;
93
94
  int streamID3 = -1;
  int vlistID3, taxisID3 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
  int nmiss;
96
  bool lvfrac = false;
97
  int nwpv; // number of words per value; real:1  complex:2
98
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
  double vfrac = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
102

  cdoInitialize(argument);

103
  cdoOperatorAdd("timrange",  func_range, DATE_LEN, NULL);
104
105
106
  cdoOperatorAdd("timmin",    func_min,   DATE_LEN, NULL);
  cdoOperatorAdd("timmax",    func_max,   DATE_LEN, NULL);
  cdoOperatorAdd("timsum",    func_sum,   DATE_LEN, NULL);
107
  cdoOperatorAdd("timmean",   func_mean,  DATE_LEN, NULL);
108
109
110
111
112
  cdoOperatorAdd("timavg",    func_avg,   DATE_LEN, NULL);
  cdoOperatorAdd("timvar",    func_var,   DATE_LEN, NULL);
  cdoOperatorAdd("timvar1",   func_var1,  DATE_LEN, NULL);
  cdoOperatorAdd("timstd",    func_std,   DATE_LEN, NULL);
  cdoOperatorAdd("timstd1",   func_std1,  DATE_LEN, NULL);
113
  cdoOperatorAdd("yearrange", func_range, YEAR_LEN, NULL);
114
115
116
  cdoOperatorAdd("yearmin",   func_min,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearmax",   func_max,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearsum",   func_sum,   YEAR_LEN, NULL);
117
  cdoOperatorAdd("yearmean",  func_mean,  YEAR_LEN, NULL);
118
119
120
121
122
  cdoOperatorAdd("yearavg",   func_avg,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearvar",   func_var,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearvar1",  func_var1,  YEAR_LEN, NULL);
  cdoOperatorAdd("yearstd",   func_std,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearstd1",  func_std1,  YEAR_LEN, NULL);
123
  cdoOperatorAdd("monrange",  func_range, MON_LEN, NULL);
124
125
126
  cdoOperatorAdd("monmin",    func_min,   MON_LEN, NULL);
  cdoOperatorAdd("monmax",    func_max,   MON_LEN, NULL);
  cdoOperatorAdd("monsum",    func_sum,   MON_LEN, NULL);
127
  cdoOperatorAdd("monmean",   func_mean,  MON_LEN, NULL);
128
129
130
131
132
  cdoOperatorAdd("monavg",    func_avg,   MON_LEN, NULL);
  cdoOperatorAdd("monvar",    func_var,   MON_LEN, NULL);
  cdoOperatorAdd("monvar1",   func_var1,  MON_LEN, NULL);
  cdoOperatorAdd("monstd",    func_std,   MON_LEN, NULL);
  cdoOperatorAdd("monstd1",   func_std1,  MON_LEN, NULL);
133
  cdoOperatorAdd("dayrange",  func_range, DAY_LEN, NULL);
134
135
136
  cdoOperatorAdd("daymin",    func_min,   DAY_LEN, NULL);
  cdoOperatorAdd("daymax",    func_max,   DAY_LEN, NULL);
  cdoOperatorAdd("daysum",    func_sum,   DAY_LEN, NULL);
137
  cdoOperatorAdd("daymean",   func_mean,  DAY_LEN, NULL);
138
139
140
141
142
  cdoOperatorAdd("dayavg",    func_avg,   DAY_LEN, NULL);
  cdoOperatorAdd("dayvar",    func_var,   DAY_LEN, NULL);
  cdoOperatorAdd("dayvar1",   func_var1,  DAY_LEN, NULL);
  cdoOperatorAdd("daystd",    func_std,   DAY_LEN, NULL);
  cdoOperatorAdd("daystd1",   func_std1,  DAY_LEN, NULL);
143
  cdoOperatorAdd("hourrange", func_range, HOUR_LEN, NULL);
144
145
146
  cdoOperatorAdd("hourmin",   func_min,   HOUR_LEN, NULL);
  cdoOperatorAdd("hourmax",   func_max,   HOUR_LEN, NULL);
  cdoOperatorAdd("hoursum",   func_sum,   HOUR_LEN, NULL);
147
  cdoOperatorAdd("hourmean",  func_mean,  HOUR_LEN, NULL);
148
149
150
151
152
  cdoOperatorAdd("houravg",   func_avg,   HOUR_LEN, NULL);
  cdoOperatorAdd("hourvar",   func_var,   HOUR_LEN, NULL);
  cdoOperatorAdd("hourvar1",  func_var1,  HOUR_LEN, NULL);
  cdoOperatorAdd("hourstd",   func_std,   HOUR_LEN, NULL);
  cdoOperatorAdd("hourstd1",  func_std1,  HOUR_LEN, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153

154
  int operatorID = cdoOperatorID();
155
156
  int operfunc   = cdoOperatorF1(operatorID);
  int comparelen = cdoOperatorF2(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157

158
  bool lrange  = operfunc == func_range;
159
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
160
161
  bool lstd    = operfunc == func_std || operfunc == func_std1;
  bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
162
  int  divisor = operfunc == func_std1 || operfunc == func_var1;
163

164
  if ( operfunc == func_mean )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
167
168
169
170
    {
      int oargc = operatorArgc();
      char **oargv = operatorArgv();

      if ( oargc == 1 )
	{
171
	  lvfrac = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174
175
176
177
178
179
	  vfrac = atof(oargv[0]);
	  if ( cdoVerbose ) cdoPrint("Set vfrac to %g", vfrac);
	  if ( vfrac < 0 || vfrac > 1 ) cdoAbort("vfrac out of range!");
	}
      else if ( oargc > 1 )
	cdoAbort("Too many arguments!");
    }

180
  int cmplen = DATE_LEN - comparelen;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
  int streamID1 = streamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183

184
185
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186

187
  if ( cmplen == 0 ) vlistDefNtsteps(vlistID2, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188

189
190
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
191
  if ( taxisInqType(taxisID2) == TAXIS_FORECAST ) taxisDefType(taxisID2, TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
193
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
194
195
  int nvars   = vlistNvars(vlistID1);
  int maxrecs = vlistNrecs(vlistID1);
196
197
198
199
200

  if ( cmplen == 0 && CDO_Reduce_Dim )
    for ( varID = 0; varID < nvars; ++varID )
      vlistDefVarTsteptype(vlistID2, varID, TSTEP_CONSTANT);

201
202
203
204
  const char *freq = NULL;
  if      ( comparelen == DAY_LEN )  freq = "day";
  else if ( comparelen == MON_LEN )  freq = "mon";
  else if ( comparelen == YEAR_LEN ) freq = "year";
205
  if ( freq ) cdiDefAttTxt(vlistID2, CDI_GLOBAL, "frequency", (int)strlen(freq), freq);
206

207
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
209
210

  streamDefVlist(streamID2, vlistID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
212
  if ( cdoDiag )
    {
213
      char filename[8192];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214

Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
      strcpy(filename, cdoOperatorName(operatorID));
      strcat(filename, "_");
217
      strcat(filename, cdoStreamName(1)->args);
218
219
220
      argument_t *fileargument = file_argument_new(filename);
      streamID3 = streamOpenWrite(fileargument, cdoFiletype());
      file_argument_free(fileargument);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
222
223

      vlistID3 = vlistDuplicate(vlistID1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
225
      for ( varID = 0; varID < nvars; ++varID )
	{
226
	  vlistDefVarDatatype(vlistID3, varID, CDI_DATATYPE_INT32);
227
	  vlistDefVarMissval(vlistID3, varID, -1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
230
231
232
	  vlistDefVarUnits(vlistID3, varID, "");
	  vlistDefVarAddoffset(vlistID3, varID, 0);
	  vlistDefVarScalefactor(vlistID3, varID, 1);
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
234
235
      taxisID3 = taxisDuplicate(taxisID1);
      vlistDefTaxis(vlistID3, taxisID3);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
      streamDefVlist(streamID3, vlistID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
238
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
239
  recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
240
241
242
243

  dtlist_type *dtlist = dtlist_new();
  dtlist_set_stat(dtlist, timestat_date);
  dtlist_set_calendar(dtlist, taxisInqCalendar(taxisID1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244

Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
246
  int gridsizemax = vlistGridsizeMax(vlistID1);
  if ( vlistNumber(vlistID1) != CDI_REAL ) gridsizemax *= 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
248
249
250
  int FIELD_MEMTYPE = 0;
  if ( CDO_Memtype == MEMTYPE_FLOAT ) FIELD_MEMTYPE = MEMTYPE_FLOAT;

251
  field_type field;
252
  field_init(&field);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
253
254
  field.memtype = FIELD_MEMTYPE;
  if ( FIELD_MEMTYPE == MEMTYPE_FLOAT )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
    field.ptrf = (float*) Malloc(gridsizemax*sizeof(float));
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
256
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
257
    field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258

259
  field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
  field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
261
  field_type **vars2 = NULL;
262
  if ( lvarstd || lrange ) vars2 = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263

264
265
  int tsID  = 0;
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
  while ( TRUE )
    {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
268
      int nsets = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
270
      while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
	{
271
	  dtlist_taxisInqTimestep(dtlist, taxisID1, nsets);
272
273
	  int vdate = dtlist_get_vdate(dtlist, nsets);
	  int vtime = dtlist_get_vtime(dtlist, nsets);
274

275
276
	  if ( nsets == 0 ) SET_DATE(indate2, vdate, vtime);
	  SET_DATE(indate1, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277

278
	  if ( DATE_IS_NEQ(indate1, indate2, cmplen) ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
280
	  for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
282
283
284
285
	    {
	      streamInqRecord(streamID1, &varID, &levelID);

	      if ( tsID == 0 )
		{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
286
287
                  recinfo[recID].varID   = varID;
                  recinfo[recID].levelID = levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
289
		}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
              field_type *psamp1 = &samp1[varID][levelID];
291
292
293
294
              field_type *pvars1 = &vars1[varID][levelID];
              field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;

	      nwpv     = pvars1->nwpv;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
	      int gridsize = pvars1->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
297
298

	      if ( nsets == 0 )
		{
299
300
		  streamReadRecord(streamID1, pvars1->ptr, &nmiss);
		  pvars1->nmiss = (size_t)nmiss;
301
302
                  if ( lrange )
                    {
303
                      pvars2->nmiss = (size_t)nmiss;
304
		      for ( int i = 0; i < nwpv*gridsize; i++ )
305
                        pvars2->ptr[i] = pvars1->ptr[i];
306
307
                    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
		  if ( nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310
311
		      if ( psamp1->ptr == NULL )
			psamp1->ptr = (double*) Malloc(nwpv*gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312

313
		      for ( int i = 0; i < nwpv*gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314
                        psamp1->ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
318
		    }
		}
	      else
		{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
319
                  if ( CDO_Memtype == MEMTYPE_FLOAT )
320
                    streamReadRecordF(streamID1, field.ptrf, &nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
321
                  else
322
323
                    streamReadRecord(streamID1, field.ptr, &nmiss);
                  field.nmiss   = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
		  field.size    = gridsize;
325
326
		  field.grid    = pvars1->grid;
		  field.missval = pvars1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
		  if ( field.nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
		      if ( psamp1->ptr == NULL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
			{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
			  psamp1->ptr = (double*) Malloc(nwpv*gridsize*sizeof(double));
332
			  for ( int i = 0; i < nwpv*gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
			    psamp1->ptr[i] = nsets;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
			}

336
		      for ( int i = 0; i < nwpv*gridsize; i++ )
337
			if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
			  psamp1->ptr[i]++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
340
		    }

341
		  if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
		    {
343
344
		      farsumq(pvars2, field);
		      farsum(pvars1, field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
		    }
346
347
                  else if ( lrange )
                    {
348
349
                      farmin(pvars2, field);
                      farmax(pvars1, field);
350
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
352
		  else
		    {
353
		      farfun(pvars1, field, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
354
355
356
357
		    }
		}
	    }

358
	  if ( nsets == 0 && lvarstd )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
359
360
361
362
            for ( int recID = 0; recID < maxrecs; recID++ )
              {
                int varID   = recinfo[recID].varID;
                int levelID = recinfo[recID].levelID;
363
364
                field_type *pvars1 = &vars1[varID][levelID];
                field_type *pvars2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
365

Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
		if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
367

368
                farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
370
371
372
373
374
375
376
377
378
	      }

	  vdate0 = vdate;
	  vtime0 = vtime;
	  nsets++;
	  tsID++;
	}

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
      for ( int recID = 0; recID < maxrecs; recID++ )
        {
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
          field_type *psamp1 = &samp1[varID][levelID];
          field_type *pvars1 = &vars1[varID][levelID];
          field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;

          if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;

          if ( lmean )
            {
              if ( psamp1->ptr ) fardiv(pvars1, *psamp1);
              else               farcdiv(pvars1, (double)nsets);
            }
          else if ( lvarstd )
            {
              if ( psamp1->ptr )
                {
                  if ( lstd ) farstd(pvars1, *pvars2, *psamp1, divisor);
                  else        farvar(pvars1, *pvars2, *psamp1, divisor);
                }
              else
                {
                  if ( lstd ) farcstd(pvars1, *pvars2, nsets, divisor);
                  else        farcvar(pvars1, *pvars2, nsets, divisor);
                }
            }
          else if ( lrange )
            {
              farsub(pvars1, *pvars2);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412

413
414
      if ( cdoVerbose )
	{
415
	  char vdatestr[32], vtimestr[32];
416
417
418
419
	  date2str(vdate0, vdatestr, sizeof(vdatestr));
	  time2str(vtime0, vtimestr, sizeof(vtimestr));
	  cdoPrint("%s %s  vfrac = %g, nsets = %d", vdatestr, vtimestr, vfrac, nsets);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420

421
      if ( lvfrac && operfunc == func_mean )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
422
423
424
425
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
            field_type *psamp1 = &samp1[varID][levelID];
427
            field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
428

Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
430

Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
432
433
            int nwpv     = pvars1->nwpv;
            int gridsize = gridInqSize(pvars1->grid);
            double missval = pvars1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
            if ( psamp1->ptr )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
435
436
437
438
              {
                int irun = 0;
                for ( int i = 0; i < nwpv*gridsize; ++i )
                  {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
                    if ( (psamp1->ptr[i] / nsets) < vfrac )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
440
                      {
441
                        pvars1->ptr[i] = missval;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
442
443
444
445
446
447
448
449
                        irun++;
                      }
                  }

                if ( irun )
                  {
                    nmiss = 0;
                    for ( int i = 0; i < nwpv*gridsize; ++i )
450
451
                      if ( DBL_IS_EQUAL(pvars1->ptr[i], missval) ) nmiss++;
                    pvars1->nmiss = (size_t)nmiss;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
452
                  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
454
455
	      }
	  }

456
      dtlist_stat_taxisDefTimestep(dtlist, taxisID2, nsets);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
458
459
460
      streamDefTimestep(streamID2, otsID);

      if ( cdoDiag )
	{
461
	  dtlist_stat_taxisDefTimestep(dtlist, taxisID3, nsets);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
463
464
	  streamDefTimestep(streamID3, otsID);
	}

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
465
      for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
466
	{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
467
468
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
          field_type *psamp1 = &samp1[varID][levelID];
470
          field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471

Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
	  if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
473
474

	  streamDefRecord(streamID2, varID, levelID);
475
	  streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
476
          
477
	  if ( cdoDiag )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
	    {
479
              double *sampptr = field.ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
	      if ( psamp1->ptr ) sampptr = psamp1->ptr;
481
482
              else
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
                  int gridsize = pvars1->size;
484
485
486
487
488
                  for ( int i = 0; i < gridsize; ++i ) sampptr[i] = nsets;
                }

              streamDefRecord(streamID3, varID, levelID);
              streamWriteRecord(streamID3, sampptr, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
490
491
492
	    }
	}

      if ( nrecs == 0 ) break;
493
      otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
495
496
    }


497
498
  field_free(vars1, vlistID1);
  field_free(samp1, vlistID1);
499
  if ( lvarstd || lrange ) field_free(vars2, vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500

501
502
  dtlist_delete(dtlist);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
503
  if ( cdoDiag ) streamClose(streamID3);
504
505
506
  streamClose(streamID2);
  streamClose(streamID1);

507
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
508
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
509
510
511

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
512
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513
}