Comp.c 8.8 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-2007 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
29
/*
   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
*/


Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <string.h>

#include "cdi.h"
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


void *Comp(void *argument)
{
  static char func[] = "Comp";
  int EQ, NE, LE, LT, GE, GT;
  int operatorID;
43
44
45
46
  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
47
48
49
50
  int nrecs, nrecs2, nvars = 0, nlev, recID;
  int tsID;
  int varID, levelID;
  int offset;
51
52
53
  int vlistIDx1, vlistIDx2, vlistID1, vlistID2, vlistID3;
  int taxisIDx1, taxisID1, taxisID2, taxisID3;
  int nmiss1, nmiss2, nmiss3;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
  int i;
  double missval1, missval2 = 0;
56
57
58
  double *missvalx1, *missvalx2;
  double *arrayx1, *arrayx2, *array1, *array2, *array3;
  double **vardata = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

  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));
  if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
  streamID2 = streamOpenRead(cdoStreamName(1));
  if ( streamID2 < 0 ) cdiError(streamID2, "Open failed on %s", cdoStreamName(1));

76
77
78
79
80
81
  streamIDx1 = streamID1;
  streamIDx2 = streamID2;

  missvalx1 = &missval1;
  missvalx2 = &missval2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
  vlistID1 = streamInqVlist(streamID1);
  vlistID2 = streamInqVlist(streamID2);
84
85
  vlistIDx1 = vlistID1;
  vlistIDx2 = vlistID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87

  taxisID1 = vlistInqTaxis(vlistID1);
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
  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
106

107
  if ( filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
    vlistCompare(vlistID1, vlistID2, func_sft);

  nospec(vlistID1);
111
  nospec(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
  gridsize = vlistGridsizeMax(vlistIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
116
117
118

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

119
120
121
  arrayx1 = array1;
  arrayx2 = array2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
  if ( cdoVerbose )
123
124
    cdoPrint("Number of timesteps: file1 %d, file2 %d",
	     vlistNtsteps(vlistID1), vlistNtsteps(vlistID2));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125

126
  if ( filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
    {
128
129
130
131
132
133
134
135
      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
136
	{
137
138
139
140
141
142
143
	  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
144
145
	}

146
147
148
149
150
151
152
153
154
155
156
      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
157
158
    }

159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  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());
  if ( streamID3 < 0 ) cdiError(streamID3, "Open failed on %s", cdoStreamName(2));

  streamDefVlist(streamID3, vlistID3);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
  tsID = 0;
178
  while ( (nrecs = streamInqTimestep(streamIDx1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
    {
180
      if ( tsID == 0 || filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
	{
182
	  nrecs2 = streamInqTimestep(streamIDx2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
184
185
186
	  if ( nrecs2 == 0 )
	    cdoAbort("Input streams have different number of timesteps!");
	}
	  
187
      taxisCopyTimestep(taxisID3, taxisIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
189
190
191
192

      streamDefTimestep(streamID3, tsID);

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

196
	  if ( tsID == 0 || filltype == FILL_NONE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
	    {
198
	      if ( recID == 0 || filltype != FILL_REC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
		{
200
201
		  streamInqRecord(streamIDx2, &varID, &levelID);
		  streamReadRecord(streamIDx2, arrayx2, &nmiss2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
203
		}

204
	      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(vardata[varID]+offset, arrayx2, gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
210
		}
	    }
211
	  else if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
	    {
213
	      gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
	      offset   = gridsize*levelID;
215
	      memcpy(arrayx2, vardata[varID]+offset, gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
217
	    }

218
219
	  gridsize1 = gridInqSize(vlistInqVarGrid(vlistIDx1, varID));
	  *missvalx1 = vlistInqVarMissval(vlistIDx1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220

221
	  if ( filltype == FILL_REC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
	    {
223
224
	      gridsize2 = gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
	      *missvalx2 = vlistInqVarMissval(vlistIDx2, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
225
	    }
226
227
228
229
230
231
232
233
234
235
	  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
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
286
287
288
289
290
291
292

	  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);

293
  if ( vardata )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
    {
295
296
      for ( varID = 0; varID < nvars; varID++ ) free(vardata[varID]);
      free(vardata);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
298
299
300
301
302
303
304
305
306
    }

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

  cdoFinish();

  return (0);
}