XTimstat.cc 19.6 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
18
19
20
21
22
23
24
25
26
  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.
*/

/*
   This module contains the following operators:

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


#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
72
#include "cdo_task.h"
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
73
74
#include "pstream.h"
//#include "pstream_write.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75

Uwe Schulzweida's avatar
Uwe Schulzweida committed
76

Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79
80
81
typedef struct {
  short varID;
  short levelID;
} recinfo_t;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
  int tsIDnext;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
  int streamID, nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
  recinfo_t *recinfo;
86
  field_type **vars;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
}
readarg_t;
89

Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
91
92
93
94
95
96
static int num_recs = 0;

static
void *cdoReadTimestep(void *rarg)
{
  int varID, levelID, nmiss;
  readarg_t *readarg = (readarg_t *) rarg;
97
  field_type **input_vars = readarg->vars;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
  recinfo_t *recinfo = readarg->recinfo;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
  int streamID = readarg->streamID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
  int tsIDnext = readarg->tsIDnext;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
102
  int nrecs = readarg->nrecs;

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
103
  // timer_start(timer_read);
104

Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
107
  for ( int recID = 0; recID < nrecs; ++recID )
    {
      streamInqRecord(streamID, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
111
112
113
114

      if ( tsIDnext == 1 && recinfo )
        {
          recinfo[recID].varID   = varID;
          recinfo[recID].levelID = levelID;
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
      if ( CDO_Memtype == MEMTYPE_FLOAT )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
        streamReadRecordF(streamID, (float*)input_vars[varID][levelID].ptr2, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
        streamReadRecord(streamID, (double*)input_vars[varID][levelID].ptr2, &nmiss);
119
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
      input_vars[varID][levelID].nmiss2 = nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
    }

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
123
  // timer_stop(timer_read);
124

Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
  num_recs = streamInqTimestep(streamID, tsIDnext);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
129
  return ((void *) &num_recs);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
static
131
void cdoUpdateVars(int nvars, int vlistID, field_type **vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
{
133
  void *tmp = NULL;
134

Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
137
138
139
  for ( int varID = 0; varID < nvars; varID++ )
    {
      int nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID, varID));
      for ( int levelID = 0; levelID < nlevels; levelID++ )
        {
140
141
142
          if ( CDO_Memtype == MEMTYPE_FLOAT )
            {
              tmp = vars[varID][levelID].ptrf;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
              vars[varID][levelID].ptrf = (float*) vars[varID][levelID].ptr2;
144
145
146
147
            }
          else
            {
              tmp = vars[varID][levelID].ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
              vars[varID][levelID].ptr = (double*) vars[varID][levelID].ptr2;
149
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150
151
152
153
154
155
          vars[varID][levelID].ptr2  = tmp;
          vars[varID][levelID].nmiss = vars[varID][levelID].nmiss2;
        }
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
156

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
158
void *XTimstat(void *argument)
{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
159
  enum {HOUR_LEN=4, DAY_LEN=6, MON_LEN=8, YEAR_LEN=10};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
160
161
162
  int timestat_date = TIMESTAT_MEAN;
  int vdate = 0, vtime = 0;
  int vdate0 = 0, vtime0 = 0;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
163
  int varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
165
166
  int streamID3 = -1;
  int vlistID3, taxisID3 = -1;
  int nmiss;
167
  bool lvfrac = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
169
170
171
172
173
  int nwpv; // number of words per value; real:1  complex:2
  char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
  double vfrac = 1;

  cdoInitialize(argument);

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
174
175
176
  cdoOperatorAdd("xtimmin",    func_min,   DATE_LEN, NULL);
  cdoOperatorAdd("xtimmax",    func_max,   DATE_LEN, NULL);
  cdoOperatorAdd("xtimsum",    func_sum,   DATE_LEN, NULL);
177
  cdoOperatorAdd("xtimmean",   func_mean,  DATE_LEN, NULL);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
178
179
180
181
182
183
184
185
  cdoOperatorAdd("xtimavg",    func_avg,   DATE_LEN, NULL);
  cdoOperatorAdd("xtimvar",    func_var,   DATE_LEN, NULL);
  cdoOperatorAdd("xtimvar1",   func_var1,  DATE_LEN, NULL);
  cdoOperatorAdd("xtimstd",    func_std,   DATE_LEN, NULL);
  cdoOperatorAdd("xtimstd1",   func_std1,  DATE_LEN, NULL);
  cdoOperatorAdd("xyearmin",   func_min,   YEAR_LEN, NULL);
  cdoOperatorAdd("xyearmax",   func_max,   YEAR_LEN, NULL);
  cdoOperatorAdd("xyearsum",   func_sum,   YEAR_LEN, NULL);
186
  cdoOperatorAdd("xyearmean",  func_mean,  YEAR_LEN, NULL);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
187
188
189
190
191
192
193
194
  cdoOperatorAdd("xyearavg",   func_avg,   YEAR_LEN, NULL);
  cdoOperatorAdd("xyearvar",   func_var,   YEAR_LEN, NULL);
  cdoOperatorAdd("xyearvar1",  func_var1,  YEAR_LEN, NULL);
  cdoOperatorAdd("xyearstd",   func_std,   YEAR_LEN, NULL);
  cdoOperatorAdd("xyearstd1",  func_std1,  YEAR_LEN, NULL);
  cdoOperatorAdd("xmonmin",    func_min,   MON_LEN, NULL);
  cdoOperatorAdd("xmonmax",    func_max,   MON_LEN, NULL);
  cdoOperatorAdd("xmonsum",    func_sum,   MON_LEN, NULL);
195
  cdoOperatorAdd("xmonmean",   func_mean,  MON_LEN, NULL);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
196
197
198
199
200
  cdoOperatorAdd("xmonavg",    func_avg,   MON_LEN, NULL);
  cdoOperatorAdd("xmonvar",    func_var,   MON_LEN, NULL);
  cdoOperatorAdd("xmonvar1",   func_var1,  MON_LEN, NULL);
  cdoOperatorAdd("xmonstd",    func_std,   MON_LEN, NULL);
  cdoOperatorAdd("xmonstd1",   func_std1,  MON_LEN, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
202
203
204
205

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

206
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
207
208
209
  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;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210

211
  if ( operfunc == func_mean )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
213
214
215
216
217
    {
      int oargc = operatorArgc();
      char **oargv = operatorArgv();

      if ( oargc == 1 )
	{
218
	  lvfrac = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
221
222
223
224
225
226
227
228
	  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!");
    }

  int cmplen = DATE_LEN - comparelen;

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
229
230
  int streamID1 = streamOpenRead(cdoStreamName(0));
  //int streamID1 = streamOpenRead(cdoStreamName(0)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
233
234
235
236
237
238
239
240
241

  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);

  if ( cmplen == 0 ) vlistDefNtsteps(vlistID2, 1);

  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
  if ( taxisInqType(taxisID2) == TAXIS_FORECAST ) taxisDefType(taxisID2, TAXIS_RELATIVE);
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
  int nvars = vlistNvars(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
245
246
247

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

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
248
249
250
251
  const char *freq = NULL;
  if      ( comparelen == DAY_LEN )  freq = "day";
  else if ( comparelen == MON_LEN )  freq = "mon";
  else if ( comparelen == YEAR_LEN ) freq = "year";
252
  if ( freq ) cdiDefAttTxt(vlistID2, CDI_GLOBAL, "frequency", (int)strlen(freq), freq);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
253

Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());

  streamDefVlist(streamID2, vlistID2);

  if ( cdoDiag )
    {
      char filename[8192];
      strcpy(filename, cdoOperatorName(operatorID));
      strcat(filename, "_");
      strcat(filename, cdoStreamName(1)->args);
      argument_t *fileargument = file_argument_new(filename);
      streamID3 = streamOpenWrite(fileargument, cdoFiletype());
      file_argument_free(fileargument);

      vlistID3 = vlistDuplicate(vlistID1);

      for ( varID = 0; varID < nvars; ++varID )
	{
272
	  vlistDefVarDatatype(vlistID3, varID, CDI_DATATYPE_INT32);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
	  vlistDefVarMissval(vlistID3, varID, -1);
	  vlistDefVarUnits(vlistID3, varID, "");
	  vlistDefVarAddoffset(vlistID3, varID, 0);
	  vlistDefVarScalefactor(vlistID3, varID, 1);
	}

      taxisID3 = taxisDuplicate(taxisID1);
      vlistDefTaxis(vlistID3, taxisID3);

      streamDefVlist(streamID3, vlistID3);
    }

  dtlist_type *dtlist = dtlist_new();
  dtlist_set_stat(dtlist, timestat_date);
  dtlist_set_calendar(dtlist, taxisInqCalendar(taxisID1));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
  int gridsizemax = vlistGridsizeMax(vlistID1);
  if ( vlistNumber(vlistID1) != CDI_REAL ) gridsizemax *= 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
291

292
293
  int FIELD_MEMTYPE = 0;
  if ( CDO_Memtype == MEMTYPE_FLOAT ) FIELD_MEMTYPE = FIELD_FLT;
294
295
296
297
  field_type **input_vars = field_malloc(vlistID1, FIELD_PTR | FIELD_PTR2 | FIELD_MEMTYPE);
  field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
  field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
  field_type **vars2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
299
  if ( lvarstd ) vars2 = field_malloc(vlistID1, FIELD_PTR);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
300
301
302
303
  readarg_t readarg;
  readarg.streamID = streamID1;
  readarg.vars = input_vars;

304
  int lparallelread = CDO_Parallel_Read;
305
  bool ltsfirst = true;
306
307
  void *read_task = NULL;
  void *readresult = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
309
310

  if ( lparallelread )
    {
311
312
313
314
315
316
      read_task = cdo_task_new();
      if ( read_task == NULL )
        {
          lparallelread = FALSE;
          cdoWarning("CDO tasks not available!");
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318

Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
320
  int tsID  = 0;
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
322
  int nrecs = streamInqTimestep(streamID1, tsID);
  int maxrecs = nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
  recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
  tsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
326
327
  while ( TRUE )
    {
328
      int nsets = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
      while ( nrecs > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
331
332
333
334
335
336
337
338
339
	{
	  dtlist_taxisInqTimestep(dtlist, taxisID1, nsets);
	  vdate = dtlist_get_vdate(dtlist, nsets);
	  vtime = dtlist_get_vtime(dtlist, nsets);

	  if ( nsets == 0 ) SET_DATE(indate2, vdate, vtime);
	  SET_DATE(indate1, vdate, vtime);

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
          readarg.tsIDnext = tsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
341
342
          readarg.nrecs    = nrecs;
          readarg.recinfo  = recinfo;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343

Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
          if ( ltsfirst || lparallelread == FALSE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
            {
346
              ltsfirst = false;
347
348
              readresult = cdoReadTimestep(&readarg);
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349
350
          else
            {
351
              readresult = cdo_task_wait(read_task);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
            }
353
354
          
          nrecs = *(int *)readresult;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355
356
357
358
359
360

          cdoUpdateVars(nvars, vlistID1, input_vars);

          if ( nrecs && lparallelread )
            {
              readarg.tsIDnext = tsID+1;
361
              cdo_task_start(read_task, cdoReadTimestep, &readarg);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
362
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363

364
          if ( nsets == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
            {
366
#if defined(_OPENMP)
367
#pragma omp parallel for default(none) shared(maxrecs, recinfo, input_vars, vars1, samp1) if(maxrecs>1)
368
369
#endif
              for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
                {
371
372
373
                  int varID    = recinfo[recID].varID;
                  int levelID  = recinfo[recID].levelID;

374
                  field_type *pvars1 = &vars1[varID][levelID];
375
                  field_type *pinput_var = &input_vars[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
376

377
378
                  int nwpv     = pvars1->nwpv;
                  int gridsize = pvars1->size;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
379
380
                  int nmiss    = pinput_var->nmiss;

381
382
                  farcpy(pvars1, *pinput_var);
                  pvars1->nmiss = nmiss;
383
                  if ( nmiss > 0 || samp1[varID][levelID].ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384
                    {
385
                      if ( samp1[varID][levelID].ptr == NULL )
386
                        samp1[varID][levelID].ptr = (double*) malloc(nwpv*gridsize*sizeof(double));
387
388
                      
                      for ( int i = 0; i < nwpv*gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389
                        samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
                    }
391
392
393
394
395
                }
            }
          else
            {
#if defined(_OPENMP)
396
#pragma omp parallel for default(none) shared(lvarstd, nsets, maxrecs, recinfo, input_vars, vars1, samp1, vars2, operfunc) if(maxrecs>1)
397
398
399
400
401
402
#endif
              for ( int recID = 0; recID < maxrecs; recID++ )
                {
                  int varID    = recinfo[recID].varID;
                  int levelID  = recinfo[recID].levelID;
                  
403
                  field_type *pvars1 = &vars1[varID][levelID];
404
                  field_type *pinput_var = &input_vars[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
405

406
407
                  int nwpv     = pvars1->nwpv;
                  int gridsize = pvars1->size;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
408
409
                  int nmiss    = pinput_var->nmiss;

410
                  if ( nmiss > 0 || samp1[varID][levelID].ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
                    {
412
                      if ( samp1[varID][levelID].ptr == NULL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
                        {
414
                          samp1[varID][levelID].ptr = (double*) malloc(nwpv*gridsize*sizeof(double));
415
416
                          for ( int i = 0; i < nwpv*gridsize; i++ )
                            samp1[varID][levelID].ptr[i] = nsets;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
                        }
418
419
                          
                      for ( int i = 0; i < nwpv*gridsize; i++ )
420
                        if ( !DBL_IS_EQUAL(pinput_var->ptr[i], pvars1->missval) )
421
422
                          samp1[varID][levelID].ptr[i]++;
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423
                  
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
424
425
		  if ( lvarstd )
		    {
426
427
428
                      field_type *pvars2 = &vars2[varID][levelID];
		      farsumq(pvars2, *pinput_var);
		      farsum(pvars1, *pinput_var);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
429
430
431
		    }
		  else
		    {
432
		      farfun(pvars1, *pinput_var, operfunc);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
433
		    }
434
435
                }
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
437

	  if ( nsets == 0 && lvarstd )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
438
439
440
441
            for ( int recID = 0; recID < maxrecs; recID++ )
              {
                int varID   = recinfo[recID].varID;
                int levelID = recinfo[recID].levelID;
442
443
                field_type *pvars1 = &vars1[varID][levelID];
                field_type *pvars2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
444

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

447
                farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
449
450
451
452
453
454
455
456
457
458
	      }

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

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

      if ( lmean )
459
        {
460
#if defined(_OPENMP)
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
461
#pragma omp parallel for default(none) shared(vlistID1, nsets, maxrecs, recinfo, vars1, samp1) if(maxrecs>1)
462
#endif
463
464
465
466
          for ( int recID = 0; recID < maxrecs; recID++ )
            {
              int varID   = recinfo[recID].varID;
              int levelID = recinfo[recID].levelID;
467
              field_type *pvars1 = &vars1[varID][levelID];
468

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

471
              if ( samp1[varID][levelID].ptr == NULL )
472
                farcdiv(pvars1, (double)nsets);
473
              else
474
                fardiv(pvars1, samp1[varID][levelID]);
475
476
            }
        }
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
477
      else if ( lvarstd )
478
479
480
481
482
        {
          for ( int recID = 0; recID < maxrecs; recID++ )
            {
              int varID   = recinfo[recID].varID;
              int levelID = recinfo[recID].levelID;
483
484
              field_type *pvars1 = &vars1[varID][levelID];
              field_type *pvars2 = &vars2[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
485

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

488
489
              if ( samp1[varID][levelID].ptr == NULL )
                {
490
491
                  if ( lstd ) farcstd(pvars1, *pvars2, nsets, divisor);
                  else        farcvar(pvars1, *pvars2, nsets, divisor);
492
493
494
                }
              else
                {
495
496
                  if ( lstd ) farstd(pvars1, *pvars2, samp1[varID][levelID], divisor);
                  else        farvar(pvars1, *pvars2, samp1[varID][levelID], divisor);
497
498
499
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
501
502
503
504
505
506
507
508

      if ( cdoVerbose )
	{
	  char vdatestr[32], vtimestr[32];
	  date2str(vdate0, vdatestr, sizeof(vdatestr));
	  time2str(vtime0, vtimestr, sizeof(vtimestr));
	  cdoPrint("%s %s  vfrac = %g, nsets = %d", vdatestr, vtimestr, vfrac, nsets);
	}

509
      if ( lvfrac && operfunc == func_mean )
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
510
511
512
513
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
514
            field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
515

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

518
            nwpv     = pvars1->nwpv;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
520
            int gridsize = pvars1->size;
            double missval = pvars1->missval;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
521
522
523
524
525
526
527
            if ( samp1[varID][levelID].ptr )
              {
                int irun = 0;
                for ( int i = 0; i < nwpv*gridsize; ++i )
                  {
                    if ( (samp1[varID][levelID].ptr[i] / nsets) < vfrac )
                      {
528
                        pvars1->ptr[i] = missval;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
529
530
531
532
533
534
535
536
                        irun++;
                      }
                  }

                if ( irun )
                  {
                    nmiss = 0;
                    for ( int i = 0; i < nwpv*gridsize; ++i )
537
538
                      if ( DBL_IS_EQUAL(pvars1->ptr[i], missval) ) nmiss++;
                    pvars1->nmiss = nmiss;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
539
                  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
541
542
543
544
545
546
547
548
549
550
551
	      }
	  }

      dtlist_stat_taxisDefTimestep(dtlist, taxisID2, nsets);
      streamDefTimestep(streamID2, otsID);

      if ( cdoDiag )
	{
	  dtlist_stat_taxisDefTimestep(dtlist, taxisID3, nsets);
	  streamDefTimestep(streamID3, otsID);
	}

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
552
      for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
553
	{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
554
555
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
556
          field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
557

Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
559
	  if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
560
          streamDefRecord(streamID2, varID, levelID);
561
	  streamWriteRecord(streamID2, pvars1->ptr,  pvars1->nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
562
563
              
          if ( cdoDiag )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
            {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
565
              if ( samp1[varID][levelID].ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
                {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
567
568
                  streamDefRecord(streamID3, varID, levelID);
                  streamWriteRecord(streamID3, samp1[varID][levelID].ptr, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
569
570
                }
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
573
574
575
576
577
	}

      if ( nrecs == 0 ) break;
      otsID++;
    }


Uwe Schulzweida's avatar
Uwe Schulzweida committed
578
  field_free(input_vars, vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
579
580
581
582
583
  field_free(vars1, vlistID1);
  field_free(samp1, vlistID1);
  if ( lvarstd ) field_free(vars2, vlistID1);

  dtlist_delete(dtlist);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
584
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585

Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
587
  if ( cdoDiag ) pstreamClose(streamID3);
  pstreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
589
590
591
592
593
  streamClose(streamID1);

  cdoFinish();

  return 0;
}