Consecstat.c 7.98 KB
Newer Older
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-2014 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
  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:

      Consectstep  consecsum  For each timestep, the current number of 
                              onsecutive timsteps is counted
      Consectstep  consects   For each period of consecutive timesteps, only
                              count its lenght + last contributing timesteps

   =============================================================================
   Created:  04/08/2010 11:58:01 AM
    Author:  Ralf Mueller (ram), ralf.mueller@zmaw.de
   Company:  Max-Planck-Institute for Meteorology
   =============================================================================
 */

Ralf Mueller's avatar
Ralf Mueller committed
33
#include <cdi.h>
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"

enum {CONSECSUM, CONSECTS};
#define SWITCHWARN "Hit default case!This should never happen (%s).\n"

static void selEndOfPeriod(field_t *periods, field_t history, field_t current, int isLastTimestep)
{
  long   i, len;
  double pmissval = periods->missval;
  double  *parray = periods->ptr;
  double hmissval = history.missval;
  double  *harray = history.ptr;
  double cmissval = current.missval;
  double  *carray = current.ptr;

  len = gridInqSize(periods->grid);
  if ( len != gridInqSize(current.grid) || (gridInqSize(current.grid) != gridInqSize(history.grid)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
    cdoAbort("Fields have different gridsize (%s)", __func__);
54
55
56
57
58

  if (!isLastTimestep)
  {
    if ( current.nmiss > 0 || history.nmiss > 0 )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
#if defined(_OPENMP)
60
#pragma omp parallel for default(shared)
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#endif
      for ( i = 0; i < len; i++ )
      {
        if ( !DBL_IS_EQUAL(harray[i], hmissval) )
        {
          if ( !DBL_IS_EQUAL(carray[i], cmissval) )
          {
            parray[i] = ( DBL_IS_EQUAL(carray[i], 0.0)  && IS_NOT_EQUAL(harray[i], 0.0) ) ? harray[i] : pmissval;
          }
          else /* DBL_IS_EQUAL(carray[i], cmissval) */
          {
            parray[i] = ( IS_NOT_EQUAL(harray[i], 0.0) ) ? harray[i] : pmissval;
          }
        }
        else /* DBL_IS_EQUAL(harray[i], hmissval) */
        {
          parray[i] = pmissval;
        }
      }
    }
    else
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
#if defined(_OPENMP)
84
#pragma omp parallel for default(shared)
85
86
87
88
89
90
91
92
93
#endif
      for ( i = 0; i < len; i++ )
        parray[i] = ( DBL_IS_EQUAL(carray[i], 0.0)  && IS_NOT_EQUAL(harray[i], 0.0) ) ? harray[i] : pmissval;
    }
  }
  else
  {
    if ( current.nmiss > 0 )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
#if defined(_OPENMP)
95
#pragma omp parallel for default(shared)
96
97
98
99
100
101
102
103
104
105
106
107
108
#endif
      for ( i = 0; i < len; i++ )
        if ( !DBL_IS_EQUAL(carray[i], cmissval) )
        {
          parray[i] = ( DBL_IS_EQUAL(carray[i], 0.0) ) ? pmissval : carray[i];
        }
        else /* DBL_IS_EQUAL(carray[i], cmissval) */
        {
          parray[i] = pmissval;
        }
    }
    else
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
#if defined(_OPENMP)
110
#pragma omp parallel for default(shared)
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#endif
      for ( i = 0; i < len; i++ )
        parray[i] = ( DBL_IS_EQUAL(carray[i], 0.0) ) ? pmissval : carray[i];
    }
  }

  periods->nmiss = 0;
  for ( i = 0; i < len; i++ )
    if ( DBL_IS_EQUAL(parray[i], pmissval) ) periods->nmiss++;
}

void *Consecstat (void *argument)
{
  int operatorID;
  int i;
  int istreamID, itaxisID, ivlistID, itsID;
  int ostreamID, otaxisID, ovlistID, otsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
  int vdate = 0, vtime = 0;
129
130
131
132
  int histvdate = 0, histvtime = 0;
  int recID, nrecs;
  int varID, nvars;
  int levelID, nlevels; 
Ralf Mueller's avatar
Ralf Mueller committed
133
  double refval = 0.0;
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

  field_t **vars = NULL, **hist = NULL, **periods = NULL;
  field_t field;

  cdoInitialize(argument);
  cdoOperatorAdd("consecsum",CONSECSUM, 0, "refval");
  cdoOperatorAdd("consects" ,CONSECTS , 0, NULL);
  operatorID = cdoOperatorID();

  if ( operatorID == CONSECSUM )
    if ( operatorArgc() > 0 ) refval = atof(operatorArgv()[0]);

  istreamID = streamOpenRead(cdoStreamName(0));

  ivlistID = streamInqVlist(istreamID);
  itaxisID = vlistInqTaxis(ivlistID);
  ovlistID = vlistDuplicate(ivlistID);
  otaxisID = taxisDuplicate(itaxisID);
  vlistDefTaxis(ovlistID, otaxisID);

154
  field_init(&field);
155
  field.ptr = (double *) malloc(vlistGridsizeMax(ovlistID)*sizeof(double));
156

157
  nvars     = vlistNvars(ivlistID);
158
  vars      = field_calloc(ivlistID, FIELD_PTR);
159
  if ( operatorID == CONSECTS )
160
    {
161
162
      hist      = field_malloc(ivlistID, FIELD_PTR);
      periods   = field_malloc(ivlistID, FIELD_PTR);
163
164
    }

165
  for ( varID = 0; varID < nvars; varID++ )
166
    {
167
      vlistDefVarUnits(ovlistID, varID, "steps"); /* TODO */
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

  ostreamID = streamOpenWrite(cdoStreamName(1), cdoFiletype());

  streamDefVlist(ostreamID, ovlistID);

  itsID = 0;
  otsID = 0;
  while ( (nrecs = streamInqTimestep(istreamID, itsID)) )
  {
    vdate = taxisInqVdate(itaxisID);
    vtime = taxisInqVtime(itaxisID);
    switch (operatorID)
    {
      case CONSECSUM:
        taxisDefVdate(otaxisID, vdate);
        taxisDefVtime(otaxisID, vtime);
        streamDefTimestep(ostreamID, otsID);
        break;
      case CONSECTS:
        if (itsID != 0 )
        {
          taxisDefVdate(otaxisID, histvdate);
          taxisDefVtime(otaxisID, histvtime);
          streamDefTimestep(ostreamID, otsID-1);
        }
        break;
      default:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
        printf (SWITCHWARN,__func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
        break;
198
199
200
201
202
203
204
205
    }

    for ( recID = 0; recID < nrecs; recID++ )
    {
      streamInqRecord(istreamID, &varID, &levelID);
      streamReadRecord(istreamID, field.ptr, &field.nmiss);
      field.grid    = vlistInqVarGrid(ovlistID, varID);
      field.missval = vlistInqVarMissval(ovlistID, varID);
206

207
      farsumtr(&vars[varID][levelID], field, refval);
208

209
210
211
212
213
214
215
216
217
218
219
220
221
      switch (operatorID)
      {
        case CONSECSUM:
          streamDefRecord(ostreamID, varID, levelID);
          streamWriteRecord(ostreamID, vars[varID][levelID].ptr, vars[varID][levelID].nmiss);
          break;
        case CONSECTS:
          if ( itsID != 0 )
          {
            selEndOfPeriod(&periods[varID][levelID], hist[varID][levelID], vars[varID][levelID], FALSE);
            streamDefRecord(ostreamID, varID, levelID);
            streamWriteRecord(ostreamID, periods[varID][levelID].ptr, periods[varID][levelID].nmiss);
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
#if defined(_OPENMP)
223
224
225
226
#pragma omp parallel for default(shared) schedule(static)
          for ( i = 0; i < gridInqSize(vars[varID][levelID].grid); i++ )
            hist[varID][levelID].ptr[i] = vars[varID][levelID].ptr[i];
#else
227
228
229
          memcpy(hist[varID][levelID].ptr,
                 vars[varID][levelID].ptr,
                 gridInqSize(vars[varID][levelID].grid)*sizeof(double));
230
#endif
231
232
          break;
        default:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
          printf (SWITCHWARN,__func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234
          break;
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
      }
    }
    histvdate = vdate;
    histvtime = vtime;
    itsID++;
    otsID++;
  }

  if ( operatorID == CONSECTS ) /* Save the last timestep */
  {
    taxisDefVdate(otaxisID, vdate);
    taxisDefVtime(otaxisID, vtime);
    streamDefTimestep(ostreamID, otsID-1);
    for ( varID = 0; varID < nvars; varID++ )
    {
      nlevels = zaxisInqSize(vlistInqVarZaxis(ovlistID, varID));
      for ( levelID = 0; levelID < nlevels; levelID++ )
      {
        selEndOfPeriod(&periods[varID][levelID], hist[varID][levelID], vars[varID][levelID], TRUE);
        streamDefRecord(ostreamID, varID, levelID);
        streamWriteRecord(ostreamID, periods[varID][levelID].ptr, periods[varID][levelID].nmiss);
      }
    }
  }
259
  
260
261
262
  if ( vars )    field_free(vars, ivlistID);
  if ( hist )    field_free(hist, ivlistID);
  if ( periods ) field_free(periods, ivlistID);
263
264
265
266
267
268
269
270

  streamClose(istreamID);
  streamClose(ostreamID);

  cdoFinish();

  return (0);
}