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

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

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

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

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

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
103
  double missval;

  cdoInitialize(argument);

104
  cdoOperatorAdd("timrange",  func_range, DATE_LEN, NULL);
105
106
107
  cdoOperatorAdd("timmin",    func_min,   DATE_LEN, NULL);
  cdoOperatorAdd("timmax",    func_max,   DATE_LEN, NULL);
  cdoOperatorAdd("timsum",    func_sum,   DATE_LEN, NULL);
108
  cdoOperatorAdd("timmean",   func_mean,  DATE_LEN, NULL);
109
110
111
112
113
  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);
114
  cdoOperatorAdd("yearrange", func_range, YEAR_LEN, NULL);
115
116
117
  cdoOperatorAdd("yearmin",   func_min,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearmax",   func_max,   YEAR_LEN, NULL);
  cdoOperatorAdd("yearsum",   func_sum,   YEAR_LEN, NULL);
118
  cdoOperatorAdd("yearmean",  func_mean,  YEAR_LEN, NULL);
119
120
121
122
123
  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);
124
  cdoOperatorAdd("monrange",  func_range, MON_LEN, NULL);
125
126
127
  cdoOperatorAdd("monmin",    func_min,   MON_LEN, NULL);
  cdoOperatorAdd("monmax",    func_max,   MON_LEN, NULL);
  cdoOperatorAdd("monsum",    func_sum,   MON_LEN, NULL);
128
  cdoOperatorAdd("monmean",   func_mean,  MON_LEN, NULL);
129
130
131
132
133
  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);
134
  cdoOperatorAdd("dayrange",  func_range, DAY_LEN, NULL);
135
136
137
  cdoOperatorAdd("daymin",    func_min,   DAY_LEN, NULL);
  cdoOperatorAdd("daymax",    func_max,   DAY_LEN, NULL);
  cdoOperatorAdd("daysum",    func_sum,   DAY_LEN, NULL);
138
  cdoOperatorAdd("daymean",   func_mean,  DAY_LEN, NULL);
139
140
141
142
143
  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);
144
  cdoOperatorAdd("hourrange", func_range, HOUR_LEN, NULL);
145
146
147
  cdoOperatorAdd("hourmin",   func_min,   HOUR_LEN, NULL);
  cdoOperatorAdd("hourmax",   func_max,   HOUR_LEN, NULL);
  cdoOperatorAdd("hoursum",   func_sum,   HOUR_LEN, NULL);
148
  cdoOperatorAdd("hourmean",  func_mean,  HOUR_LEN, NULL);
149
150
151
152
153
  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
154

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

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

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

      if ( oargc == 1 )
	{
172
	  lvfrac = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
174
175
176
177
178
179
180
	  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!");
    }

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

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

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

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

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

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

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

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

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

  streamDefVlist(streamID2, vlistID2);

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

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

      vlistID3 = vlistDuplicate(vlistID1);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

291
              field_type *pvar1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
292
293
294
              
	      nwpv     = pvar1->nwpv;
	      gridsize = gridInqSize(pvar1->grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
296
297

	      if ( nsets == 0 )
		{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
298
		  streamReadRecord(streamID1, pvar1->ptr, &nmiss);
299
		  pvar1->nmiss = (size_t)nmiss;
300
301
                  if ( lrange )
                    {
302
                      field_type *pvar2 = &vars2[varID][levelID];
303
304
305
306
307
		      for ( int i = 0; i < nwpv*gridsize; i++ )
                        pvar2->ptr[i] = pvar1->ptr[i];
                      pvar2->nmiss = (size_t)nmiss;
                    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
309
310
		  if ( nmiss > 0 || samp1[varID][levelID].ptr )
		    {
		      if ( samp1[varID][levelID].ptr == NULL )
311
			samp1[varID][levelID].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
cleanup    
Uwe Schulzweida committed
314
			if ( DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
318
319
320
321
			  samp1[varID][levelID].ptr[i] = 0;
			else
			  samp1[varID][levelID].ptr[i] = 1;
		    }
		}
	      else
		{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
322
                  if ( CDO_Memtype == MEMTYPE_FLOAT )
323
                    streamReadRecordF(streamID1, field.ptrf, &nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
324
                  else
325
326
                    streamReadRecord(streamID1, field.ptr, &nmiss);
                  field.nmiss   = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
		  field.size    = gridsize;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
328
329
		  field.grid    = pvar1->grid;
		  field.missval = pvar1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
331
332
333
		  if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
		    {
		      if ( samp1[varID][levelID].ptr == NULL )
			{
334
			  samp1[varID][levelID].ptr = (double*) Malloc(nwpv*gridsize*sizeof(double));
335
			  for ( int i = 0; i < nwpv*gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
337
338
			    samp1[varID][levelID].ptr[i] = nsets;
			}

339
		      for ( int i = 0; i < nwpv*gridsize; i++ )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
340
			if ( !DBL_IS_EQUAL(field.ptr[i], pvar1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
341
342
343
			  samp1[varID][levelID].ptr[i]++;
		    }

344
		  if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
		    {
346
                      field_type *pvar2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
347
348
		      farsumq(pvar2, field);
		      farsum(pvar1, field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349
		    }
350
351
                  else if ( lrange )
                    {
352
                      field_type *pvar2 = &vars2[varID][levelID];
353
354
355
                      farmin(pvar2, field);
                      farmax(pvar1, field);
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
357
		  else
		    {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
358
		      farfun(pvar1, field, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
360
361
362
		    }
		}
	    }

363
	  if ( nsets == 0 && lvarstd )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
364
365
366
367
            for ( int recID = 0; recID < maxrecs; recID++ )
              {
                int varID   = recinfo[recID].varID;
                int levelID = recinfo[recID].levelID;
368
369
                field_type *pvar1 = &vars1[varID][levelID];
                field_type *pvar2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
370

Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
		if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
372
373

                farmoq(pvar2, *pvar1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
375
376
377
378
379
380
381
382
383
	      }

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

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

384
      if ( lmean )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
385
386
387
388
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
389
            field_type *pvar1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
390

Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
392
393
394
395
396

            if ( samp1[varID][levelID].ptr == NULL )
              farcdiv(pvar1, (double)nsets);
            else
              fardiv(pvar1, samp1[varID][levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
	  }
398
      else if ( lvarstd )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
399
400
401
402
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
403
404
            field_type *pvar1 = &vars1[varID][levelID];
            field_type *pvar2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
405
406
407
408
409
410
411
412
413
414
415
416

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

            if ( samp1[varID][levelID].ptr == NULL )
              {
                if ( lstd ) farcstd(pvar1, *pvar2, nsets, divisor);
                else        farcvar(pvar1, *pvar2, nsets, divisor);
              }
            else
              {
                if ( lstd ) farstd(pvar1, *pvar2, samp1[varID][levelID], divisor);
                else        farvar(pvar1, *pvar2, samp1[varID][levelID], divisor);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
418
	      }
	  }
419
420
421
422
423
      else if ( lrange )
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
424
425
            field_type *pvar1 = &vars1[varID][levelID];
            field_type *pvar2 = &vars2[varID][levelID];
426
427
428
429
430

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

            farsub(pvar1, *pvar2);
	  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431

432
433
      if ( cdoVerbose )
	{
434
	  char vdatestr[32], vtimestr[32];
435
436
437
438
	  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
439

440
      if ( lvfrac && operfunc == func_mean )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
441
442
443
444
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
445
            field_type *pvar1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
446

Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468

            nwpv     = pvar1->nwpv;
            gridsize = gridInqSize(pvar1->grid);
            missval  = pvar1->missval;
            if ( samp1[varID][levelID].ptr )
              {
                int irun = 0;
                for ( int i = 0; i < nwpv*gridsize; ++i )
                  {
                    if ( (samp1[varID][levelID].ptr[i] / nsets) < vfrac )
                      {
                        pvar1->ptr[i] = missval;
                        irun++;
                      }
                  }

                if ( irun )
                  {
                    nmiss = 0;
                    for ( int i = 0; i < nwpv*gridsize; ++i )
                      if ( DBL_IS_EQUAL(pvar1->ptr[i], missval) ) nmiss++;
469
                    pvar1->nmiss = (size_t)nmiss;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
470
                  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
472
473
	      }
	  }

474
      dtlist_stat_taxisDefTimestep(dtlist, taxisID2, nsets);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
476
477
478
      streamDefTimestep(streamID2, otsID);

      if ( cdoDiag )
	{
479
	  dtlist_stat_taxisDefTimestep(dtlist, taxisID3, nsets);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
481
482
	  streamDefTimestep(streamID3, otsID);
	}

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
483
      for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484
	{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
485
486
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
487
          field_type *pvar1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488

Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
	  if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
490
491

	  streamDefRecord(streamID2, varID, levelID);
492
	  streamWriteRecord(streamID2, pvar1->ptr, (int)pvar1->nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
493
          
494
	  if ( cdoDiag )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
	    {
496
497
498
499
500
501
502
503
504
505
              double *sampptr = field.ptr;
	      if ( samp1[varID][levelID].ptr ) sampptr = samp1[varID][levelID].ptr;
              else
                {
                  gridsize = gridInqSize(pvar1->grid);
                  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
506
507
508
509
	    }
	}

      if ( nrecs == 0 ) break;
510
      otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
513
    }


514
515
  field_free(vars1, vlistID1);
  field_free(samp1, vlistID1);
516
  if ( lvarstd || lrange ) field_free(vars2, vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517

518
519
  dtlist_delete(dtlist);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
  if ( cdoDiag ) streamClose(streamID3);
521
522
523
  streamClose(streamID2);
  streamClose(streamID1);

524
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
525
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526
527
528

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
}