Ydrunpctl.c 9.94 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) 2006 Brockmann Consult
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
  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:

      Ydrunpctl    ydrunpctl         Multi-year daily running percentiles
*/

Ralf Mueller's avatar
Ralf Mueller committed
24
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
28
#include "percentiles_hist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29

Uwe Schulzweida's avatar
Uwe Schulzweida committed
30

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
#define NDAY 373

33
34
35
36
37
38
39
int getmonthday(int date)
{
  int year, month, day;
  cdiDecodeDate(date, &year, &month, &day);
  return month*100+day;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
41
42
43
44
45
46

void *Ydrunpctl(void *argument)
{
  int gridsize;
  int varID;
  int recID;
  int gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
  int levelID;
  int tsID;
  int otsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
  int inp, its;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
  int nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
  int nlevels;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
  field_t ***vars1 = NULL, **vars2[NDAY];
55
  datetime_t *datetime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
  int vdate, vtime;
  int vdates1[NDAY], vtimes1[NDAY];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
  int vdates2[NDAY] /*, vtimes2[NDAY]*/;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
  int nsets[NDAY];
  int year, month, day, dayoy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
  field_t field;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
63
64
65
66
67
68
  HISTOGRAM_SET *hsets[NDAY];
    
  cdoInitialize(argument);
  cdoOperatorAdd("ydrunpctl", func_pctl, 0, NULL);

  operatorInputArg("percentile number, number of timesteps");
  operatorCheckArgc(2);
69
  double pn  = parameter2double(operatorArgv()[0]);
70
  percentile_check_number(pn);
71
  int ndates = parameter2int(operatorArgv()[1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
73
74
75
76
77
78
79
  
  for ( dayoy = 0; dayoy < NDAY; dayoy++ )
    {
      vars2[dayoy] = NULL;
      hsets[dayoy] = NULL;
      nsets[dayoy] = 0;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
  int streamID1 = streamOpenRead(cdoStreamName(0));
  int streamID2 = streamOpenRead(cdoStreamName(1));
  int streamID3 = streamOpenRead(cdoStreamName(2));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83

Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
85
86
87
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = streamInqVlist(streamID2);
  int vlistID3 = streamInqVlist(streamID3);
  int vlistID4 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88

Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
  vlistCompare(vlistID1, vlistID2, CMP_ALL);
  vlistCompare(vlistID1, vlistID3, CMP_ALL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91

Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
93
94
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = vlistInqTaxis(vlistID2);
  int taxisID3 = vlistInqTaxis(vlistID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
  /* TODO - check that time axes 2 and 3 are equal */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
  int taxisID4 = taxisDuplicate(taxisID1);
98
  if ( taxisHasBounds(taxisID4) ) taxisDeleteBounds(taxisID4);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
  vlistDefTaxis(vlistID4, taxisID4);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
102
  int calendar = taxisInqCalendar(taxisID1);
  int dpy      = calendar_dpy(calendar);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103

Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
  int streamID4 = streamOpenWrite(cdoStreamName(3), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
107

  streamDefVlist(streamID4, vlistID4);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
  int nvars    = vlistNvars(vlistID1);
  int nrecords = vlistNrecs(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110

111
112
  int *recVarID   = (int*) Malloc(nrecords*sizeof(int));
  int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
114

  gridsize = vlistGridsizeMax(vlistID1);
115
  field_init(&field);
116
  field.ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117

118
  datetime = (datetime_t*) Malloc((ndates+1)*sizeof(datetime_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
  
120
  vars1 = (field_t ***) Malloc((ndates+1)*sizeof(field_t **));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
123
  
  for ( its = 0; its < ndates; its++ )
    {
124
      vars1[its] = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
127
    }

  tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
  while ( (nrecs = streamInqTimestep(streamID2, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
    {
      if ( nrecs != streamInqTimestep(streamID3, tsID) )
131
        cdoAbort("Number of records at time step %d of %s and %s differ!", tsID+1, cdoStreamName(1)->args, cdoStreamName(2)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
133
134
135
136
      
      vdate = taxisInqVdate(taxisID2);
      vtime = taxisInqVtime(taxisID2);
      
      if ( vdate != taxisInqVdate(taxisID3) || vtime != taxisInqVtime(taxisID3) )
137
        cdoAbort("Verification dates at time step %d of %s and %s differ!", tsID+1, cdoStreamName(1)->args, cdoStreamName(2)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
140
        
      if ( cdoVerbose ) cdoPrint("process timestep: %d %d %d", tsID+1, vdate, vtime);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
      cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
144
145
146
147
148
149
150
151

      if ( month >= 1 && month <= 12 )
	dayoy = (month-1)*31 + day;
      else
	dayoy = 0;

      if ( dayoy < 0 || dayoy >= NDAY )
	cdoAbort("Day %d out of range!", dayoy);

      vdates2[dayoy] = vdate;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
      //vtimes2[dayoy] = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
155

      if ( vars2[dayoy] == NULL )
	{
156
	  vars2[dayoy] = field_malloc(vlistID2, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
          hsets[dayoy] = hsetCreate(nvars);

	  for ( varID = 0; varID < nvars; varID++ )
	    {
	      gridID   = vlistInqVarGrid(vlistID2, varID);
	      nlevels  = zaxisInqSize(vlistInqVarZaxis(vlistID2, varID));

              hsetCreateVarLevels(hsets[dayoy], varID, nlevels, gridID);
	    }
	}
      
      for ( recID = 0; recID < nrecs; recID++ )
        {
          streamInqRecord(streamID2, &varID, &levelID);
	  streamReadRecord(streamID2, vars2[dayoy][varID][levelID].ptr, &nmiss);
          vars2[dayoy][varID][levelID].nmiss = nmiss;
        }
      for ( recID = 0; recID < nrecs; recID++ )
        {
          streamInqRecord(streamID3, &varID, &levelID);
	  streamReadRecord(streamID3, field.ptr, &nmiss);
          field.nmiss   = nmiss;
          field.grid    = vars2[dayoy][varID][levelID].grid;
	  field.missval = vars2[dayoy][varID][levelID].missval;
	  
	  hsetDefVarLevelBounds(hsets[dayoy], varID, levelID, &vars2[dayoy][varID][levelID], &field);
        }
      
      tsID++;
    }
  
  for ( tsID = 0; tsID < ndates; tsID++ )
    {
      nrecs = streamInqTimestep(streamID1, tsID);
      if ( nrecs == 0 )
	cdoAbort("File has less then %d timesteps!", ndates);

      datetime[tsID].date = taxisInqVdate(taxisID1);
      datetime[tsID].time = taxisInqVtime(taxisID1);
	
      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);

	  if ( tsID == 0 )
	    {
	      recVarID[recID]   = varID;
	      recLevelID[recID] = levelID;
	    }
	  
	  streamReadRecord(streamID1, vars1[tsID][varID][levelID].ptr, &nmiss);
	  vars1[tsID][varID][levelID].nmiss = nmiss;
	}
    }
  
  while ( TRUE )
    {
      datetime_avg(dpy, ndates, datetime);
      
      vdate = datetime[ndates].date;
      vtime = datetime[ndates].time;
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
      cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220
221
222
223
224
225
226
227
228
229
230
231
232

      if ( month >= 1 && month <= 12 )
	dayoy = (month-1)*31 + day;
      else
	dayoy = 0;

      if ( dayoy < 0 || dayoy >= NDAY )
	cdoAbort("Day %d out of range!", dayoy);

      vdates1[dayoy] = vdate;
      vtimes1[dayoy] = vtime;
      
      if ( vars2[dayoy] == NULL )
233
        cdoAbort("No data for day %d in %s and %s", dayoy, cdoStreamName(1)->args, cdoStreamName(2)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234
235
236

      for ( varID = 0; varID < nvars; varID++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
	  if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
240
	  nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
	      
	  for ( levelID = 0; levelID < nlevels; levelID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
	    for ( inp = 0; inp < ndates; inp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
	      hsetAddVarLevelValues(hsets[dayoy], varID, levelID, &vars1[inp][varID][levelID]);
	}
        
      datetime[ndates] = datetime[0];
      vars1[ndates] = vars1[0];

      for ( inp = 0; inp < ndates; inp++ )
	{
	  datetime[inp] = datetime[inp+1];
	  vars1[inp] = vars1[inp+1];
	}

      nrecs = streamInqTimestep(streamID1, tsID);
      if ( nrecs == 0 ) break;

      datetime[ndates-1].date = taxisInqVdate(taxisID1);
      datetime[ndates-1].time = taxisInqVtime(taxisID1);

      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);
	  
	  streamReadRecord(streamID1, vars1[ndates-1][varID][levelID].ptr, &nmiss);
	  vars1[ndates-1][varID][levelID].nmiss = nmiss;
	}

      nsets[dayoy] += ndates;
      tsID++;
    }
271
272
273
274
275
276
277
278
279
  /*
  int outyear = 1e9;
  for ( dayoy = 0; dayoy < NDAY; dayoy++ )
    if ( nsets[dayoy] )
      {
	int year, month, day;
	cdiDecodeDate(vdates1[dayoy], &year, &month, &day);
	if ( year < outyear ) outyear = year;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280

281
282
283
284
285
286
287
288
  for ( dayoy = 0; dayoy < NDAY; dayoy++ )
    if ( nsets[dayoy] )
      {
	int year, month, day;
	cdiDecodeDate(vdates1[dayoy], &year, &month, &day);
	vdates1[dayoy] = cdiEncodeDate(outyear, month, day);
      }
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
291
292
  otsID = 0;
  for ( dayoy = 0; dayoy < NDAY; dayoy++ )
    if ( nsets[dayoy] )
      {
293
        if ( getmonthday(vdates1[dayoy]) != getmonthday(vdates2[dayoy]) )
294
295
296
          cdoAbort("Verification dates for day %d of %s, %s and %s are different!",
                   dayoy, cdoStreamName(1)->args, cdoStreamName(2)->args, cdoStreamName(3)->args);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
298
	for ( varID = 0; varID < nvars; varID++ )
	  {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
300
301
302
303
304
305
306
307
	    nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
	      
	    for ( levelID = 0; levelID < nlevels; levelID++ )
	      hsetGetVarLevelPercentiles(&vars2[dayoy][varID][levelID], hsets[dayoy], varID, levelID, pn);
	  }

	taxisDefVdate(taxisID4, vdates1[dayoy]);
	taxisDefVtime(taxisID4, vtimes1[dayoy]);
308
	streamDefTimestep(streamID4, otsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312
313
314

	for ( recID = 0; recID < nrecords; recID++ )
	  {
	    varID    = recVarID[recID];
	    levelID  = recLevelID[recID];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
	    if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
316
317
318
319

	    streamDefRecord(streamID4, varID, levelID);
	    streamWriteRecord(streamID4, vars2[dayoy][varID][levelID].ptr,
			      vars2[dayoy][varID][levelID].nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320
	  }
321
322

	otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
324
325
326
327
328
329
330
      }
  
  for ( its = 0; its < ndates; its++ )
    {
      for ( varID = 0; varID < nvars; varID++ )
	{
	  nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
	  for ( levelID = 0; levelID < nlevels; levelID++ )
331
332
	    Free(vars1[its][varID][levelID].ptr);
	  Free(vars1[its][varID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
	}
334
      Free(vars1[its]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
    }
336
  Free(vars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
338
339
340
341

  for ( dayoy = 0; dayoy < NDAY; dayoy++ )
    {
      if ( vars2[dayoy] != NULL )
	{
342
	  field_free(vars2[dayoy], vlistID2); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
344
345
346
	  hsetDestroy(hsets[dayoy]);
	}
    }

347
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
  
349
350
  if ( recVarID   ) Free(recVarID);
  if ( recLevelID ) Free(recLevelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
352
353
354
355
356
357
358

  streamClose(streamID4);
  streamClose(streamID3);
  streamClose(streamID2);
  streamClose(streamID1);

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360
}