Ensstat.c 8.81 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
21
22
23
24
25
26
/*
   This module contains the following operators:

      Ensstat    ensmin          Ensemble minimum
      Ensstat    ensmax          Ensemble maximum
      Ensstat    enssum          Ensemble sum
      Ensstat    ensmean         Ensemble mean
      Ensstat    ensavg          Ensemble average
      Ensstat    ensstd          Ensemble standard deviation
27
      Ensstat    ensstd1         Ensemble standard deviation
28
      Ensstat    ensvar          Ensemble variance
29
      Ensstat    ensvar1         Ensemble variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
      Ensstat    enspctl         Ensemble percentiles
31
32
*/

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


void *Ensstat(void *argument)
{
42
  int nrecs0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
46
  typedef struct
  {
    int streamID;
    int vlistID;
47
    int nmiss;
48
    double missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
54
55
56
    double *array;
  } ens_file_t;

  cdoInitialize(argument);

  cdoOperatorAdd("ensmin",  func_min,  0, NULL);
  cdoOperatorAdd("ensmax",  func_max,  0, NULL);
  cdoOperatorAdd("enssum",  func_sum,  0, NULL);
57
  cdoOperatorAdd("ensmean", func_mean, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
  cdoOperatorAdd("ensavg",  func_avg,  0, NULL);
59
60
61
62
  cdoOperatorAdd("ensstd",  func_stdw,  0, NULL);
  cdoOperatorAdd("ensstd1", func_std1w, 0, NULL);
  cdoOperatorAdd("ensvar",  func_varw,  0, NULL);
  cdoOperatorAdd("ensvar1", func_var1w, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
  cdoOperatorAdd("enspctl", func_pctl, 0, NULL);
64

65
66
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67

68
69
  bool lpctl = operfunc == func_pctl;

70
71
  int argc = operatorArgc();
  int nargc = argc;
72
73

  double pn = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
76
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
77
      pn = parameter2int(operatorArgv()[0]);
78
      percentile_check_number(pn);
79
80
81
      argc--;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
  bool count_data = false;
83
84
  if ( argc == 1 )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
      if ( strcmp("count", operatorArgv()[nargc-1]) == 0 ) count_data = true;
86
      else cdoAbort("Unknown parameter: >%s<", operatorArgv()[nargc-1]); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
    }
    
89
  int nfiles = cdoStreamCnt() - 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

91
  if ( cdoVerbose ) cdoPrint("Ensemble over %d files.", nfiles);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

93
  const char *ofilename = cdoStreamName(nfiles)->args;
94

95
96
  if ( !cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename) )
    cdoAbort("Outputfile %s already exists!", ofilename);
97

98
  ens_file_t *ef = (ens_file_t *) Malloc(nfiles*sizeof(ens_file_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99

100
  field_type *field = (field_type *) Malloc(ompNumThreads*sizeof(field_type));
101
  for ( int i = 0; i < ompNumThreads; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
      field_init(&field[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
      field[i].size   = nfiles;
105
106
      field[i].ptr    = (double*) Malloc(nfiles*sizeof(double));
      field[i].weight = (double*) Malloc(nfiles*sizeof(double));
107
      for ( int fileID = 0; fileID < nfiles; fileID++ ) field[i].weight[fileID] = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
    }

110
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
    {
112
113
      ef[fileID].streamID = streamOpenRead(cdoStreamName(fileID));
      ef[fileID].vlistID  = streamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
  /* check that the contents is always the same */
117
  for ( int fileID = 1; fileID < nfiles; fileID++ )
118
    vlistCompare(ef[0].vlistID, ef[fileID].vlistID, CMP_ALL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119

120
121
122
123
  int vlistID1 = ef[0].vlistID;
  int vlistID2 = vlistDuplicate(vlistID1);
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
  vlistDefTaxis(vlistID2, taxisID2);

126
  int gridsizemax = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127

128
  for ( int fileID = 0; fileID < nfiles; fileID++ )
129
    ef[fileID].array = (double*) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130

131
  double *array2 = (double *) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132

133
134
135
136
  int nvars = vlistNvars(vlistID2);
  double *count2 = NULL;
  if ( count_data )
    {
137
      count2 = (double *) Malloc(gridsizemax*sizeof(double));
138
      for ( int varID = 0; varID < nvars; ++varID )
139
140
141
142
	{
	  char name[CDI_MAX_NAME];
	  vlistInqVarName(vlistID2, varID, name);
	  strcat(name, "_count");
143
	  int gridID = vlistInqVarGrid(vlistID2, varID);
144
145
146
147
	  int zaxisID = vlistInqVarZaxis(vlistID2, varID);
	  int tsteptype = vlistInqVarTsteptype(vlistID2, varID);
	  int cvarID = vlistDefVar(vlistID2, gridID, zaxisID, tsteptype);
	  vlistDefVarName(vlistID2, cvarID, name);
148
	  vlistDefVarDatatype(vlistID2, cvarID, CDI_DATATYPE_INT16);
149
150
151
152
153
154
155
156
	  if ( cvarID != (varID+nvars) ) cdoAbort("Internal error, varIDs do not match!");
	}
    }

  int streamID2 = streamOpenWrite(cdoStreamName(nfiles), cdoFiletype());

  streamDefVlist(streamID2, vlistID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
158
  bool lwarning = false;
  bool lerror = false;
159
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
160
161
162
  do
    {
      nrecs0 = streamInqTimestep(ef[0].streamID, tsID);
163
      for ( int fileID = 1; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
	{
165
166
	  int streamID = ef[fileID].streamID;
	  int nrecs = streamInqTimestep(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
	  if ( nrecs != nrecs0 )
168
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
	      if ( nrecs == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
171
172
173
174
175
176
177
178
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(fileID)->args);
                }
	      else if ( nrecs0 == 0 )
                {
                  lwarning = true;
                  cdoWarning("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(0)->args);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
	      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181
182
183
184
185
                {
                  lerror = true;
                  cdoWarning("Inconsistent ensemble file, number of records at time step %d of %s and %s differ!",
                             tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
                }
              goto CLEANUP;
186
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
	}

189
190
191
192
193
      if ( nrecs0 > 0 )
	{
	  taxisCopyTimestep(taxisID2, taxisID1);
	  streamDefTimestep(streamID2, tsID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194

195
      for ( int recID = 0; recID < nrecs0; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
	{
197
          int varID = 0, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

199
200
201
202
203
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
	      streamInqRecord(ef[fileID].streamID, &varID, &levelID);
              ef[fileID].missval = vlistInqVarMissval(ef[fileID].vlistID, varID);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
          //#pragma omp parallel for default(none) shared(ef, nfiles)
205
206
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
	    {
207
	      streamReadRecord(ef[fileID].streamID, ef[fileID].array, &ef[fileID].nmiss);
208
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209

210
211
212
          bool lmiss = false;
          for ( int fileID = 0; fileID < nfiles; fileID++ ) if ( ef[fileID].nmiss > 0 ) lmiss = true;

213
	  int gridID = vlistInqVarGrid(vlistID1, varID);
214
	  int gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
	  double missval = vlistInqVarMissval(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216

Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
	  int nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
#if defined(_OPENMP)
219
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220
#endif
221
	  for ( int i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
	    {
223
224
	      int ompthID = cdo_omp_get_thread_num();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
225
226
	      field[ompthID].missval = missval;
	      field[ompthID].nmiss = 0;
227
	      for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
		{
		  field[ompthID].ptr[fileID] = ef[fileID].array[i];
230
		  if ( lmiss && DBL_IS_EQUAL(field[ompthID].ptr[fileID], ef[fileID].missval) )
231
232
233
234
235
                    {
                      field[ompthID].ptr[fileID] = missval;
                      field[ompthID].nmiss++;
                    }
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236

237
              array2[i] = lpctl ? fldpctl(field[ompthID], pn) : fldfun(field[ompthID], operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238

239
240
	      if ( DBL_IS_EQUAL(array2[i], field[ompthID].missval) )
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
#if defined(_OPENMP)
242
#include "pragma_omp_atomic_update.h"
243
244
245
#endif
		  nmiss++;
		}
246
247

	      if ( count_data ) count2[i] = nfiles - field[ompthID].nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
	    }

	  streamDefRecord(streamID2, varID, levelID);
	  streamWriteRecord(streamID2, array2, nmiss);
252
253
254
255
256
257

	  if ( count_data )
	    {
	      streamDefRecord(streamID2, varID+nvars, levelID);
	      streamWriteRecord(streamID2, count2, 0);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
259
260
261
262
263
	}

      tsID++;
    }
  while ( nrecs0 > 0 );

Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
266
267
268
 CLEANUP:

  if ( lwarning ) cdoWarning("Inconsistent ensemble, processed only the first %d timesteps!", tsID);
  if ( lerror ) cdoAbort("Inconsistent ensemble, processed only the first %d timesteps!", tsID);

269
  for ( int fileID = 0; fileID < nfiles; fileID++ )
270
    streamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271
272
273

  streamClose(streamID2);

274
  for ( int fileID = 0; fileID < nfiles; fileID++ )
275
    if ( ef[fileID].array ) Free(ef[fileID].array);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276

277
278
279
  if ( ef ) Free(ef);
  if ( array2 ) Free(array2);
  if ( count2 ) Free(count2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280

281
  for ( int i = 0; i < ompNumThreads; i++ )
282
    {
283
284
      if ( field[i].ptr    ) Free(field[i].ptr);
      if ( field[i].weight ) Free(field[i].weight);
285
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286

287
  if ( field ) Free(field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
289
290

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
291
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
292
}