Comp.c 8.5 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.

5
  Copyright (C) 2003-2011 Uwe Schulzweida, Uwe.Schulzweida@zmaw.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
27
28
/*
   This module contains the following operators:

      Comp       eq              Equal
      Comp       ne              Not equal
      Comp       le              Less equal
      Comp       lt              Less than
      Comp       ge              Greater equal
      Comp       gt              Greater than
*/

Ralf Mueller's avatar
Ralf Mueller committed
29
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
31
32
33
34
35
36
37
38
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


void *Comp(void *argument)
{
  int EQ, NE, LE, LT, GE, GT;
  int operatorID;
39
40
41
42
  enum {FILL_NONE, FILL_TS, FILL_REC};
  int filltype = FILL_NONE;
  int streamIDx1, streamIDx2, streamID1, streamID2, streamID3;
  int gridsize, gridsize1, gridsize2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
46
  int nrecs, nrecs2, nvars = 0, nlev, recID;
  int tsID;
  int varID, levelID;
  int offset;
47
48
49
  int vlistIDx1, vlistIDx2, vlistID1, vlistID2, vlistID3;
  int taxisIDx1, taxisID1, taxisID2, taxisID3;
  int nmiss1, nmiss2, nmiss3;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
  int i;
  double missval1, missval2 = 0;
52
53
54
  double *missvalx1, *missvalx2;
  double *arrayx1, *arrayx2, *array1, *array2, *array3;
  double **vardata = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

  cdoInitialize(argument);

  EQ = cdoOperatorAdd("eq", 0, 0, NULL);
  NE = cdoOperatorAdd("ne", 0, 0, NULL);
  LE = cdoOperatorAdd("le", 0, 0, NULL);
  LT = cdoOperatorAdd("lt", 0, 0, NULL);
  GE = cdoOperatorAdd("ge", 0, 0, NULL);
  GT = cdoOperatorAdd("gt", 0, 0, NULL);

  operatorID = cdoOperatorID();

  streamID1 = streamOpenRead(cdoStreamName(0));
  streamID2 = streamOpenRead(cdoStreamName(1));

70
71
72
73
74
75
  streamIDx1 = streamID1;
  streamIDx2 = streamID2;

  missvalx1 = &missval1;
  missvalx2 = &missval2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
77
  vlistID1 = streamInqVlist(streamID1);
  vlistID2 = streamInqVlist(streamID2);
78
79
  vlistIDx1 = vlistID1;
  vlistIDx2 = vlistID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81

  taxisID1 = vlistInqTaxis(vlistID1);
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
  taxisID2 = vlistInqTaxis(vlistID2);
  taxisIDx1 = taxisID1;

  if ( vlistNrecs(vlistID1) != 1 && vlistNrecs(vlistID2) == 1 )
    {
      filltype = FILL_REC;
      cdoPrint("Filling up stream2 >%s< by copying the first record.", cdoStreamName(1));
    }
  else if ( vlistNrecs(vlistID1) == 1 && vlistNrecs(vlistID2) != 1 )
    {
      filltype = FILL_REC;
      cdoPrint("Filling up stream1 >%s< by copying the first record.", cdoStreamName(0));
      streamIDx1 = streamID2;
      streamIDx2 = streamID1;
      vlistIDx1 = vlistID2;
      vlistIDx2 = vlistID1;
      taxisIDx1 = taxisID2;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100

101
  if ( filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
    vlistCompare(vlistID1, vlistID2, CMP_SFT);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
104

  nospec(vlistID1);
105
  nospec(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106

107
  gridsize = vlistGridsizeMax(vlistIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
111
112

  array1 = (double *) malloc(gridsize*sizeof(double));
  array2 = (double *) malloc(gridsize*sizeof(double));
  array3 = (double *) malloc(gridsize*sizeof(double));

113
114
115
  arrayx1 = array1;
  arrayx2 = array2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
  if ( cdoVerbose )
117
118
    cdoPrint("Number of timesteps: file1 %d, file2 %d",
	     vlistNtsteps(vlistID1), vlistNtsteps(vlistID2));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119

120
  if ( filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
    {
122
123
124
125
126
127
128
129
      if ( vlistNtsteps(vlistID1) != 1 && vlistNtsteps(vlistID1) != 0 &&
	  (vlistNtsteps(vlistID2) == 1 || vlistNtsteps(vlistID2) == 0) )
	{
	  filltype = FILL_TS;
	  cdoPrint("Filling up stream2 >%s< by copying the first timestep.", cdoStreamName(1));
	}
      else if ( (vlistNtsteps(vlistID1) == 1 || vlistNtsteps(vlistID1) == 0) &&
		 vlistNtsteps(vlistID2) != 1 && vlistNtsteps(vlistID2) != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
	{
131
132
133
134
135
136
137
	  filltype = FILL_TS;
	  cdoPrint("Filling up stream1 >%s< by copying the first timestep.", cdoStreamName(0));
	  streamIDx1 = streamID2;
          streamIDx2 = streamID1;
	  vlistIDx1 = vlistID2;
	  vlistIDx2 = vlistID1;
	  taxisIDx1 = taxisID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
	}

140
141
142
143
144
145
146
147
148
149
150
      if ( filltype == FILL_TS )
	{
	  nvars  = vlistNvars(vlistIDx2);
	  vardata  = (double **) malloc(nvars*sizeof(double *));
	  for ( varID = 0; varID < nvars; varID++ )
	    {
	      gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
	      nlev     = zaxisInqSize(vlistInqVarZaxis(vlistIDx2, varID));
	      vardata[varID]  = (double *) malloc(nlev*gridsize*sizeof(double));
	    }
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
152
    }

153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
  if ( filltype != FILL_NONE && vlistNtsteps(vlistID1) == 1 )
    {
      arrayx1 = array2;
      arrayx2 = array1;
      missvalx1 = &missval2;
      missvalx2 = &missval1;
    }

  vlistID3 = vlistDuplicate(vlistIDx1);

  taxisID3 = taxisDuplicate(taxisIDx1);
  vlistDefTaxis(vlistID3, taxisID3);

  streamID3 = streamOpenWrite(cdoStreamName(2), cdoFiletype());

  streamDefVlist(streamID3, vlistID3);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
  tsID = 0;
171
  while ( (nrecs = streamInqTimestep(streamIDx1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
    {
173
      if ( tsID == 0 || filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
	{
175
	  nrecs2 = streamInqTimestep(streamIDx2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
177
178
179
	  if ( nrecs2 == 0 )
	    cdoAbort("Input streams have different number of timesteps!");
	}
	  
180
      taxisCopyTimestep(taxisID3, taxisIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
184
185

      streamDefTimestep(streamID3, tsID);

      for ( recID = 0; recID < nrecs; recID++ )
	{
186
187
	  streamInqRecord(streamIDx1, &varID, &levelID);
	  streamReadRecord(streamIDx1, arrayx1, &nmiss1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188

189
	  if ( tsID == 0 || filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
	    {
191
	      if ( recID == 0 || filltype != FILL_REC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
		{
193
194
		  streamInqRecord(streamIDx2, &varID, &levelID);
		  streamReadRecord(streamIDx2, arrayx2, &nmiss2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
196
		}

197
	      if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
		{
199
		  gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
		  offset   = gridsize*levelID;
201
		  memcpy(vardata[varID]+offset, arrayx2, gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
203
		}
	    }
204
	  else if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
	    {
206
	      gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
	      offset   = gridsize*levelID;
208
	      memcpy(arrayx2, vardata[varID]+offset, gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
210
	    }

211
212
	  gridsize1 = gridInqSize(vlistInqVarGrid(vlistIDx1, varID));
	  *missvalx1 = vlistInqVarMissval(vlistIDx1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213

214
	  if ( filltype == FILL_REC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
	    {
216
217
	      gridsize2 = gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
	      *missvalx2 = vlistInqVarMissval(vlistIDx2, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
	    }
219
220
221
222
223
224
225
226
227
228
	  else
	    {
	      gridsize2 = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
	      *missvalx2 = vlistInqVarMissval(vlistIDx2, varID);
	    }

	  if ( gridsize1 != gridsize2 ) cdoAbort("Streams have different gridsize (gridsize1 = %d; gridsize2 = %d!",
						 gridsize1, gridsize2);

	  gridsize = gridsize1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
230
231
232
233
234
235
236
237
238
239
240
241
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
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285

	  if ( operatorID == EQ )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : DBL_IS_EQUAL(array1[i], array2[i]));
	    }
	  else if ( operatorID == NE )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : !DBL_IS_EQUAL(array1[i], array2[i]));
	    }
	  else if ( operatorID == LE )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : array1[i] <= array2[i]);
	    }
	  else if ( operatorID == LT )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : array1[i] < array2[i]);
	    }
	  else if ( operatorID == GE )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : array1[i] >= array2[i]);
	    }
	  else if ( operatorID == GT )
	    {
	      for ( i = 0; i < gridsize; i++ )
		array3[i] = (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) ?
			     missval1 : array1[i] > array2[i]);
	    }
	  else
	    {
	      cdoAbort("Operator not implemented!");
	    }

	  nmiss3 = 0;
	  for ( i = 0; i < gridsize; i++ )
	    if ( DBL_IS_EQUAL(array3[i], missval1) ) nmiss3++;

	  streamDefRecord(streamID3, varID, levelID);
	  streamWriteRecord(streamID3, array3, nmiss3);
	}

      tsID++;
    }

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

286
  if ( vardata )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
287
    {
288
289
      for ( varID = 0; varID < nvars; varID++ ) free(vardata[varID]);
      free(vardata);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
291
292
293
294
295
296
297
298
299
    }

  if ( array3 ) free(array3);
  if ( array2 ) free(array2);
  if ( array1 ) free(array1);

  cdoFinish();

  return (0);
}