Merstat.cc 5.19 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:

      Merstat    mermin          Meridional minimum
      Merstat    mermax          Meridional maximum
      Merstat    mersum          Meridional sum
      Merstat    mermean         Meridional mean
      Merstat    meravg          Meridional average
      Merstat    merstd          Meridional standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
      Merstat    merstd          Meridional standard deviation [Normalize by (n-1)]
28
      Merstat    mervar          Meridional variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
      Merstat    mervar          Meridional variance [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
      Merstat    merpctl         Meridional percentiles
31
32
33
*/


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


void *Merstat(void *argument)
{
43
  int gridID1, gridID2 = -1, lastgrid = -1;
44
  int wstatus = FALSE;
45
  int index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
  int nmiss;
47
  int nrecs;
48
  int varID, levelID;
49
  char varname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52

  cdoInitialize(argument);

53
54
55
  cdoOperatorAdd("mermin",  func_min,   0, NULL);
  cdoOperatorAdd("mermax",  func_max,   0, NULL);
  cdoOperatorAdd("mersum",  func_sum,   0, NULL);
56
57
58
59
60
61
  cdoOperatorAdd("mermean", func_meanw, 1, NULL);
  cdoOperatorAdd("meravg",  func_avgw,  1, NULL);
  cdoOperatorAdd("mervar",  func_varw,  1, NULL);
  cdoOperatorAdd("mervar1", func_var1w, 1, NULL);
  cdoOperatorAdd("merstd",  func_stdw,  1, NULL);
  cdoOperatorAdd("merstd1", func_std1w, 1, NULL);
62
  cdoOperatorAdd("merpctl", func_pctl,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
 
64
65
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
66
  bool needWeights = cdoOperatorF2(operatorID) != 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67

68
  double pn = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
71
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
72
73
      pn = parameter2double(operatorArgv()[0]);
      percentile_check_number(pn);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
    }

76
  int streamID1 = streamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77

78
79
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80

81
82
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
  vlistDefTaxis(vlistID2, taxisID2);

85
86
  int ngrids = vlistNgrids(vlistID1);
  int ndiffgrids = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
89
90
91
92
  for ( index = 1; index < ngrids; index++ )
    if ( vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index))
      ndiffgrids++;

  if ( ndiffgrids > 0 ) cdoAbort("Too many different grids!");

Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
95
  index = 0;
  gridID1 = vlistGrid(vlistID1, index);

96
97
98
99
100
101
102
103
104
105
106
  if ( gridInqType(gridID1) == GRID_LONLAT   ||
       gridInqType(gridID1) == GRID_GAUSSIAN ||
       gridInqType(gridID1) == GRID_GENERIC )
    {
      gridID2 = gridToMeridional(gridID1);
    }
  else
    {
      cdoAbort("Unsupported gridtype: %s", gridNamePtr(gridInqType(gridID1)));
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
  vlistChangeGridIndex(vlistID2, index, gridID2);

109
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
111
112
113

  streamDefVlist(streamID2, vlistID2);

  gridID1 = vlistInqVarGrid(vlistID1, 0);
114
115
  int nlonmax = gridInqXsize(gridID1); /* max nlon ? */
  int lim = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116

117
  field_type field1, field2;
118
119
120
  field_init(&field1);
  field_init(&field2);

121
  field1.ptr    = (double*) Malloc(lim*sizeof(double));
122
123
  field1.weight = NULL;
  if ( needWeights )
124
    field1.weight = (double*) Malloc(lim*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125

126
  field2.ptr  = (double*) Malloc(nlonmax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
  field2.grid = gridID2;

129
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132
133
134
135
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

136
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
	{
	  streamInqRecord(streamID1, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
	  streamReadRecord(streamID1, field1.ptr, &nmiss);
          field1.nmiss = (size_t) nmiss;
141
	  field1.grid = vlistInqVarGrid(vlistID1, varID);
142
	  if ( needWeights && field1.grid != lastgrid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
144
	    {
	      lastgrid = field1.grid;
145
	      wstatus = gridWeights(field1.grid, field1.weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
	    }
147
148
149
150
151
	  if ( wstatus != 0 && tsID == 0 && levelID == 0 )
	    {
	      vlistInqVarName(vlistID1, varID, varname);
	      cdoWarning("Using constant grid cell area weights for variable %s!", varname);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
153
154
	  field1.missval = vlistInqVarMissval(vlistID1, varID);
	  field2.missval = vlistInqVarMissval(vlistID1, varID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
155
156
157
158
	  if ( operfunc == func_pctl )
	    merpctl(field1, & field2, pn);
	  else  
	    merfun(field1, &field2, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160

	  streamDefRecord(streamID2, varID,  levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
	  streamWriteRecord(streamID2, field2.ptr, (int)field2.nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
163
164
165
166
167
168
	}
      tsID++;
    }

  streamClose(streamID2);
  streamClose(streamID1);

169
170
171
  if ( field1.ptr )    Free(field1.ptr);
  if ( field1.weight ) Free(field1.weight);
  if ( field2.ptr )    Free(field2.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
}