Seasstat.cc 9.74 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:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
24
25
26
27
28
29
30
      Seasstat   seasrange       Seasonal range
      Seasstat   seasmin         Seasonal minimum
      Seasstat   seasmax         Seasonal maximum
      Seasstat   seassum         Seasonal sum
      Seasstat   seasmean        Seasonal mean
      Seasstat   seasavg         Seasonal average
      Seasstat   seasvar         Seasonal variance
      Seasstat   seasvar1        Seasonal variance [Normalize by (n-1)]
      Seasstat   seasstd         Seasonal standard deviation
      Seasstat   seasstd1        Seasonal standard deviation [Normalize by (n-1)]
31
32
33
*/


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


typedef struct {
  short varID;
  short levelID;
} recinfo_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
45
46
47


void *Seasstat(void *argument)
{
48
  int timestat_date = TIMESTAT_MEAN;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
  int vdate0 = 0, vtime0 = 0;
50
  int vdate1 = 0, vtime1 = 0;
51
  int nrecs;
52
53
  int varID, levelID;
  int year, month, day, seas0 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
  int nmiss;
55
  int oldmon = 0;
56
  int nseason = 0;
57
  const char *seas_name[4];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
59
60

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
  // clang-format off
62
63
64
65
66
67
68
69
70
71
  cdoOperatorAdd("seasrange", func_range, 0, NULL);
  cdoOperatorAdd("seasmin",   func_min,   0, NULL);
  cdoOperatorAdd("seasmax",   func_max,   0, NULL);
  cdoOperatorAdd("seassum",   func_sum,   0, NULL);
  cdoOperatorAdd("seasmean",  func_mean,  0, NULL);
  cdoOperatorAdd("seasavg",   func_avg,   0, NULL);
  cdoOperatorAdd("seasvar",   func_var,   0, NULL);
  cdoOperatorAdd("seasvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("seasstd",   func_std,   0, NULL);
  cdoOperatorAdd("seasstd1",  func_std1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72

73
74
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75

76
77
78
79
80
  bool lrange  = operfunc == func_range;
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
  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
81
82
83
84
  // clang-format on

  int season_start = get_season_start();
  get_season_name(seas_name);
85

86
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87

88
  int vlistID1 = pstreamInqVlist(streamID1);
89
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

91
92
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
93
  if ( taxisInqType(taxisID2) == TAXIS_FORECAST ) taxisDefType(taxisID2, TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
95
  vlistDefTaxis(vlistID2, taxisID2);

96
97
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98

99
  int maxrecs = vlistNrecs(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100

101
  recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
102
103
104
105

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
  int gridsizemax = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108

109
  field_type field;
110
  field_init(&field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
  field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
114
115
  field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
  field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
  field_type **vars2 = NULL;
116
  if ( lvarstd || lrange ) vars2 = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117

118
119
  int tsID  = 0;
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
  while ( TRUE )
    {
122
123
      long nsets = 0;
      bool newseas = false;
124
      while ( (nrecs = pstreamInqTimestep(streamID1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
	{
126
	  dtlist_taxisInqTimestep(dtlist, taxisID1, nsets);
127
128
	  int vdate = dtlist_get_vdate(dtlist, nsets);
	  int vtime = dtlist_get_vtime(dtlist, nsets);
129

130
	  cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131

132
	  int newmon = month;
133

134
	  if ( season_start == START_DEC && newmon == 12 ) newmon = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135

136
          int seas = month_to_season(month);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137

138
139
	  if ( nsets == 0 )
	    {
140
141
142
143
	      nseason++;
	      vdate0 = vdate;
	      vtime0 = vtime;
	      seas0  = seas;
144
145
146
	      oldmon = newmon;
	    }

147
	  if ( newmon < oldmon ) newseas = true;
148
149
150
151

	  if ( (seas != seas0) || newseas ) break;

	  oldmon = newmon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
              field_type *psamp1 = &samp1[varID][levelID];
164
165
              field_type *pvars1 = &vars1[varID][levelID];
              field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;
166

Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
	      int gridsize = pvars1->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
169
170

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
		  if ( nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
183
		      if ( psamp1->ptr == NULL )
			psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
		  if ( field.nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
		      if ( psamp1->ptr == NULL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
			{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
			  psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
201
			  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
			    psamp1->ptr[i] = nsets;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
			}

205
		      for ( int i = 0; i < gridsize; i++ )
206
			if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
			  psamp1->ptr[i]++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
209
		    }

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

227
	  if ( nsets == 0 && lvarstd )
228
229
230
231
            for ( int recID = 0; recID < maxrecs; recID++ )
              {
                int varID   = recinfo[recID].varID;
                int levelID = recinfo[recID].levelID;
232
233
                field_type *pvars1 = &vars1[varID][levelID];
                field_type *pvars2 = &vars2[varID][levelID];
234

Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
		if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
236

237
                farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
	      }

240
241
	  vdate1 = vdate;
	  vtime1 = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
244
245
246
247
	  nsets++;
	  tsID++;
	}

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
      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
281

282
      if ( cdoVerbose )
283
	{
284
285
286
287
288
289
	  char vdatestr0[32], vtimestr0[32];
	  char vdatestr1[32], vtimestr1[32];
	  date2str(vdate0, vdatestr0, sizeof(vdatestr0));
	  time2str(vtime0, vtimestr0, sizeof(vtimestr0));
	  date2str(vdate1, vdatestr1, sizeof(vdatestr1));
	  time2str(vtime1, vtimestr1, sizeof(vtimestr1));
290
	  cdoPrint("season: %3d %3s  start: %s %s  end: %s %s ntimesteps: %d", 
291
		   nseason, seas_name[seas0], vdatestr0, vtimestr0, vdatestr1, vtimestr1, nsets);
292
	}
293

294
      dtlist_stat_taxisDefTimestep(dtlist, taxisID2, nsets);
295
      pstreamDefTimestep(streamID2, otsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296

297
298
      if ( nsets < 3 )
	{
299
300
301
	  char vdatestr[32];
	  date2str(vdate0, vdatestr, sizeof(vdatestr));
	  cdoWarning("Season %3d (%s) has only %d input time step%s!", 
302
		     otsID+1, vdatestr, nsets, nsets == 1 ? "" : "s");
303
304
	}

305
      for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
	{
307
308
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
309
          field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310

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

313
314
	  pstreamDefRecord(streamID2, varID, levelID);
	  pstreamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
	}

      if ( nrecs == 0 ) break;
318
      otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
320
321
    }


322
323
  field_free(vars1, vlistID1);
  field_free(samp1, vlistID1);
324
  if ( lvarstd ) field_free(vars2, vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325

326
  dtlist_delete(dtlist);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327

328
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329

330
  if ( field.ptr ) Free(field.ptr);
331

332
333
  pstreamClose(streamID2);
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
336

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
}