Ensstat.c 8.13 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-2015 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
*/

Uwe Schulzweida's avatar
Uwe Schulzweida committed
33

Ralf Mueller's avatar
Ralf Mueller 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"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
#include "util.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
41
42
43


void *Ensstat(void *argument)
{
  int i;
44
  int varID = 0, recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
46
  int gridID;
  int nrecs, nrecs0;
47
  int levelID = 0;
48
  int streamID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
54
  int nmiss;
  double missval;
  typedef struct
  {
    int streamID;
    int vlistID;
55
    double missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
58
59
60
61
62
63
64
65
66
    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);
  cdoOperatorAdd("ensmean", func_mean, 0, NULL);
  cdoOperatorAdd("ensavg",  func_avg,  0, NULL);
  cdoOperatorAdd("ensstd",  func_std,  0, NULL);
67
  cdoOperatorAdd("ensstd1", func_std1, 0, NULL);
68
69
  cdoOperatorAdd("ensvar",  func_var,  0, NULL);
  cdoOperatorAdd("ensvar1", func_var1, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
  cdoOperatorAdd("enspctl", func_pctl, 0, NULL);
71

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

75
76
  int argc = operatorArgc();
  int nargc = argc;
77
78

  double pn = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
81
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
82
      pn = parameter2int(operatorArgv()[0]);
83
      percentile_check_number(pn);
84
85
86
87
88
89
90
91
      argc--;
    }

  int count_data = FALSE;
  if ( argc == 1 )
    {
      if ( strcmp("count", operatorArgv()[nargc-1]) == 0 ) count_data = TRUE;
      else cdoAbort("Unknown parameter: >%s<", operatorArgv()[nargc-1]); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
93
    }
    
94
  int nfiles = cdoStreamCnt() - 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

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

98
  const char *ofilename = cdoStreamName(nfiles)->args;
99

100
  if ( !cdoSilentMode && !cdoOverwriteMode )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
    if ( fileExists(ofilename) )
102
      if ( !userFileOverwrite(ofilename) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
	cdoAbort("Outputfile %s already exists!", ofilename);
104

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

107
  field_t *field = (field_t *) Malloc(ompNumThreads*sizeof(field_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
  for ( i = 0; i < ompNumThreads; i++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
      field_init(&field[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
      field[i].size   = nfiles;
112
113
      field[i].ptr    = (double*) Malloc(nfiles*sizeof(double));
      field[i].weight = (double*) Malloc(nfiles*sizeof(double));
114
      for ( int fileID = 0; fileID < nfiles; fileID++ ) field[i].weight[fileID] = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
116
    }

117
  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
    {
119
120
      ef[fileID].streamID = streamOpenRead(cdoStreamName(fileID));
      ef[fileID].vlistID  = streamInqVlist(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
    }

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

127
128
129
130
  int vlistID1 = ef[0].vlistID;
  int vlistID2 = vlistDuplicate(vlistID1);
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
  vlistDefTaxis(vlistID2, taxisID2);

133
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134

135
  for ( int fileID = 0; fileID < nfiles; fileID++ )
136
    ef[fileID].array = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137

138
  double *array2 = (double *) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139

140
141
142
143
  int nvars = vlistNvars(vlistID2);
  double *count2 = NULL;
  if ( count_data )
    {
144
      count2 = (double *) Malloc(gridsize*sizeof(double));
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
      for ( varID = 0; varID < nvars; ++varID )
	{
	  char name[CDI_MAX_NAME];
	  vlistInqVarName(vlistID2, varID, name);
	  strcat(name, "_count");
	  gridID = vlistInqVarGrid(vlistID2, varID);
	  int zaxisID = vlistInqVarZaxis(vlistID2, varID);
	  int tsteptype = vlistInqVarTsteptype(vlistID2, varID);
	  int cvarID = vlistDefVar(vlistID2, gridID, zaxisID, tsteptype);
	  vlistDefVarName(vlistID2, cvarID, name);
	  vlistDefVarDatatype(vlistID2, cvarID, DATATYPE_INT16);
	  if ( cvarID != (varID+nvars) ) cdoAbort("Internal error, varIDs do not match!");
	}
    }

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

  streamDefVlist(streamID2, vlistID2);

164
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
167
  do
    {
      nrecs0 = streamInqTimestep(ef[0].streamID, tsID);
168
      for ( int fileID = 1; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
171
172
	{
	  streamID = ef[fileID].streamID;
	  nrecs = streamInqTimestep(streamID, tsID);
	  if ( nrecs != nrecs0 )
173
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
	      if ( nrecs == 0 )
175
		cdoAbort("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(fileID)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
177
	      else
		cdoAbort("Inconsistent ensemble file, number of records at time step %d of %s and %s differ!",
178
			 tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
179
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181
	}

182
183
184
185
186
      if ( nrecs0 > 0 )
	{
	  taxisCopyTimestep(taxisID2, taxisID1);
	  streamDefTimestep(streamID2, tsID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
189

      for ( recID = 0; recID < nrecs0; recID++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
#if defined(_OPENMP)
191
192
#pragma omp parallel for default(none) shared(ef, nfiles)         \
                                         private(streamID, nmiss) \
193
194
                                     lastprivate(varID, levelID)
#endif
195
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
197
198
199
	    {
	      streamID = ef[fileID].streamID;
	      streamInqRecord(streamID, &varID, &levelID);
	      streamReadRecord(streamID, ef[fileID].array, &nmiss);
200
              ef[fileID].missval = vlistInqVarMissval(ef[fileID].vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
202
	    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
	  gridID   = vlistInqVarGrid(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
	  gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
	  missval  = vlistInqVarMissval(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207

	  nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
#if defined(_OPENMP)
209
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
#endif
211
	  for ( i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
	    {
213
214
	      int ompthID = cdo_omp_get_thread_num();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
	      field[ompthID].missval = missval;
	      field[ompthID].nmiss = 0;
217
	      for ( int fileID = 0; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
		{
		  field[ompthID].ptr[fileID] = ef[fileID].array[i];
220
221
222
223
224
225
		  if ( DBL_IS_EQUAL(field[ompthID].ptr[fileID], ef[fileID].missval) )
                    {
                      field[ompthID].ptr[fileID] = missval;
                      field[ompthID].nmiss++;
                    }
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226

Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
228
229
230
	      if ( operfunc == func_pctl )
	        array2[i] = fldpctl(field[ompthID], pn);
	      else  
	        array2[i] = fldfun(field[ompthID], operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231

232
233
	      if ( DBL_IS_EQUAL(array2[i], field[ompthID].missval) )
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234
#if defined(_OPENMP)
235
#include "pragma_omp_atomic_update.h"
236
237
238
#endif
		  nmiss++;
		}
239
240

	      if ( count_data ) count2[i] = nfiles - field[ompthID].nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
243
244
	    }

	  streamDefRecord(streamID2, varID, levelID);
	  streamWriteRecord(streamID2, array2, nmiss);
245
246
247
248
249
250

	  if ( count_data )
	    {
	      streamDefRecord(streamID2, varID+nvars, levelID);
	      streamWriteRecord(streamID2, count2, 0);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
252
253
254
255
256
	}

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

257
  for ( int fileID = 0; fileID < nfiles; fileID++ )
258
    streamClose(ef[fileID].streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
260
261

  streamClose(streamID2);

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

265
266
267
  if ( ef ) Free(ef);
  if ( array2 ) Free(array2);
  if ( count2 ) Free(count2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
268
269

  for ( i = 0; i < ompNumThreads; i++ )
270
    {
271
272
      if ( field[i].ptr    ) Free(field[i].ptr);
      if ( field[i].weight ) Free(field[i].weight);
273
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

275
  if ( field ) Free(field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276
277
278

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280
}