Ymonstat.cc 9.9 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.
*/

Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20
/*
   This module contains the following operators:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
      Ymonstat   ymonrange       Multi-year monthly range
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
23
      Ymonstat   ymonmin         Multi-year monthly minimum
      Ymonstat   ymonmax         Multi-year monthly maximum
24
      Ymonstat   ymonsum         Multi-year monthly sum
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
      Ymonstat   ymonmean        Multi-year monthly mean
      Ymonstat   ymonavg         Multi-year monthly average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
      Ymonstat   ymonvar         Multi-year monthly variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
      Ymonstat   ymonvar1        Multi-year monthly variance [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
      Ymonstat   ymonstd         Multi-year monthly standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
      Ymonstat   ymonstd1        Multi-year monthly standard deviation [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
33
*/


Ralf Mueller's avatar
Ralf Mueller committed
34
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
#include "cdo.h"
36
#include "cdo_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
38
39
40
41
#include "pstream.h"


#define  NMONTH     17

42
43
44
45
46
47
typedef struct {
  short varID;
  short levelID;
} recinfo_t;


48
/*
49
50
51
52
static
int cmpint(const void *s1, const void *s2)
{
  int cmp = 0;
53
54
  const int *x = s1;
  const int *y = s2;
55
56
57
58

  if      ( *x < *y ) cmp = -1;
  else if ( *x > *y ) cmp =  1;

59
  return cmp;
60
}
61
*/
62

Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
65
void *Ymonstat(void *argument)
{
  int varID;
66
  int year, month, day;
67
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
  int levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
  int month_nsets[NMONTH];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
  int nmiss;
  int vdates[NMONTH], vtimes[NMONTH];
72
73
  int mon[NMONTH];
  int nmon = 0;
74
75
  field_type **vars1[NMONTH], **vars2[NMONTH], **samp1[NMONTH];
  field_type field;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
77
78

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
81
82
83
84
85
86
87
88
  cdoOperatorAdd("ymonrange", func_range, 0, NULL);
  cdoOperatorAdd("ymonmin",   func_min,   0, NULL);
  cdoOperatorAdd("ymonmax",   func_max,   0, NULL);
  cdoOperatorAdd("ymonsum",   func_sum,   0, NULL);
  cdoOperatorAdd("ymonmean",  func_mean,  0, NULL);
  cdoOperatorAdd("ymonavg",   func_avg,   0, NULL);
  cdoOperatorAdd("ymonvar",   func_var,   0, NULL);
  cdoOperatorAdd("ymonvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("ymonstd",   func_std,   0, NULL);
  cdoOperatorAdd("ymonstd1",  func_std1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89

90
91
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
  bool lrange  = operfunc == func_range;
94
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
95
96
  bool lstd    = operfunc == func_std || operfunc == func_std1;
  bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
  int  divisor = operfunc == func_std1 || operfunc == func_var1;
98

Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
102
  for ( month = 0; month < NMONTH; month++ )
    {
      vars1[month] = NULL;
      vars2[month] = NULL;
103
      samp1[month] = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
      month_nsets[month] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
    }

107
  int streamID1 = streamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108

109
110
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111

112
113
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
  if ( taxisHasBounds(taxisID2) ) taxisDeleteBounds(taxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
116
  vlistDefTaxis(vlistID2, taxisID2);

117
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
119
120

  streamDefVlist(streamID2, vlistID2);

121
  int maxrecs = vlistNrecs(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122

123
  recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124

125
  int gridsizemax = vlistGridsizeMax(vlistID1);
126
  field_init(&field);
127
  field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128

129
130
  int tsID = 0;
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
133
134
      int vdate = taxisInqVdate(taxisID1);
      int vtime = taxisInqVtime(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
137

      if ( cdoVerbose ) cdoPrint("process timestep: %d %d %d", tsID+1, vdate, vtime);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
      cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
141
142
143
      if ( month < 0 || month >= NMONTH )
	cdoAbort("month %d out of range!", month);

      vdates[month] = vdate;
      vtimes[month] = vtime;
144
      // mon[month] = vdate;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146
147

      if ( vars1[month] == NULL )
	{
148
	  mon[nmon++] = month;
149
150
	  vars1[month] = field_malloc(vlistID1, FIELD_PTR);
	  samp1[month] = field_malloc(vlistID1, FIELD_NONE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
	  if ( lvarstd || lrange )
152
	    vars2[month] = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
	}

155
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
158
	{
	  streamInqRecord(streamID1, &varID, &levelID);

159
160
	  if ( tsID == 0 )
	    {
161
162
              recinfo[recID].varID   = varID;
              recinfo[recID].levelID = levelID;
163
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164

Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
          field_type *psamp1 = &samp1[month][varID][levelID];
166
167
          field_type *pvars1 = &vars1[month][varID][levelID];
          field_type *pvars2 = vars2[month] ? &vars2[month][varID][levelID] : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
          int nsets = month_nsets[month];
169
170

          int gridsize = pvars1->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171

Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
	  if ( nsets == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
	    {
174
175
	      streamReadRecord(streamID1, pvars1->ptr, &nmiss);
	      pvars1->nmiss = (size_t) nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
177
178
179
180
181
182
183
              if ( lrange )
                {
                  pvars2->nmiss = pvars1->nmiss;
                  for ( int i = 0; i < gridsize; i++ )
                    pvars2->ptr[i] = pvars1->ptr[i];
                }

	      if ( nmiss > 0 || psamp1->ptr )
184
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
		  if ( psamp1->ptr == NULL )
		    psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
187

188
		  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
                    psamp1->ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
190
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
191
192
193
	    }
	  else
	    {
194
195
	      streamReadRecord(streamID1, field.ptr, &nmiss);
              field.nmiss   = (size_t) nmiss;
196
197
	      field.grid    = pvars1->grid;
	      field.missval = pvars1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
	      if ( field.nmiss > 0 || psamp1->ptr )
200
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
		  if ( psamp1->ptr == NULL )
202
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
		      psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
204
		      for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
			psamp1->ptr[i] = nsets;
206
207
		    }
		  
208
209
		  for ( int i = 0; i < gridsize; i++ )
		    if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
		      psamp1->ptr[i]++;
211
212
		}

213
	      if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
		{
215
216
		  farsumq(pvars2, field);
		  farsum(pvars1, field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
220
221
222
              else if ( lrange )
                {
                  farmin(pvars2, field);
                  farmax(pvars1, field);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224
	      else
		{
225
		  farfun(pvars1, field, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
228
229
		}
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
      if ( month_nsets[month] == 0 && lvarstd )
231
232
233
234
235
236
237
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[month][varID][levelID];
            field_type *pvars2 = &vars2[month][varID][levelID];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
239
240

            farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
      month_nsets[month]++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
246
      tsID++;
    }

247
248
249
  if ( nmon == 12 )
    {
      int smon = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
      for ( month = 1; month <= 12; month++ ) if ( month_nsets[month] ) smon++;
251
252
253
      if ( smon == 12 ) for ( month = 1; month <= 12; month++ ) mon[month-1] = month;
    }

254
255
256
257
258
  /* sort output time steps */
  /*
  nmon = 0;
  for ( month = 0; month < NMONTH; month++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
      if ( month_nsets[month] == 0 )
260
261
262
263
264
265
266
267
268
269
270
271
272
	for ( i = month+1; i < NMONTH; i++ ) mon[i-1] = mon[i];
      else
	nmon++;
    }

  qsort(mon, nmon, sizeof(int), cmpint);
	      
  for ( i = 0; i < nmon; i++ )
    {
      cdiDecodeDate(mon[i], &year, &month, &day);
      mon[i] = month;
    }
  */
273
  for ( int i = 0; i < nmon; i++ )
274
275
    {
      month = mon[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276
      if ( month_nsets[month] == 0 ) cdoAbort("Internal problem, nsets[%d] not defined!", month);
277

Uwe Schulzweida's avatar
Uwe Schulzweida committed
278
279
280
281
282
283
284
285
      int nsets = month_nsets[month];
      for ( int recID = 0; recID < maxrecs; recID++ )
        {
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
          field_type *psamp1 = &samp1[month][varID][levelID];
          field_type *pvars1 = &vars1[month][varID][levelID];
          field_type *pvars2 = vars2[month] ? &vars2[month][varID][levelID] : NULL;
286

Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
          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);
            }
        }
312
313
314

      taxisDefVdate(taxisID2, vdates[month]);
      taxisDefVtime(taxisID2, vtimes[month]);
315
      streamDefTimestep(streamID2, otsID);
316

317
      for ( int recID = 0; recID < maxrecs; recID++ )
318
	{
319
320
321
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
          field_type *pvars1 = &vars1[month][varID][levelID];
322
	  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
	  if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
324
325

	  streamDefRecord(streamID2, varID, levelID);
326
	  streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
327
	}
328
329

      otsID++;
330
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
332
333
334
335

  for ( month = 0; month < NMONTH; month++ )
    {
      if ( vars1[month] != NULL )
	{
336
337
338
	  field_free(vars1[month], vlistID1);
	  field_free(samp1[month], vlistID1);
	  if ( lvarstd ) field_free(vars2[month], vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
340
341
	}
    }

342
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343

344
345
346
347
  Free(recinfo);

  streamClose(streamID2);
  streamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
349
350

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
}