Arith.cc 10.3 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:

      Arith      add             Add two fields
      Arith      sub             Subtract two fields
      Arith      mul             Multiply two fields
      Arith      div             Divide two fields
      Arith      min             Minimum of two fields
      Arith      max             Maximum of two fields
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
      Arith      atan2           Arc tangent of two fields
28
29
*/

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


void *Arith(void *argument)
{
38
  enum {FILL_NONE, FILL_TS, FILL_VAR, FILL_VARTS, FILL_FILE};
39
  int filltype = FILL_NONE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
  int nmiss;
41
  int nrecs, nrecs2, nvars = 0;
42
  int nlevels2 = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
  int varID, levelID;
44
  int levelID2;
45
  int *varnmiss2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
  int **varnmiss = NULL;
47
  double *vardata2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
51
  double **vardata = NULL;

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
54
55
56
57
58
59
  cdoOperatorAdd("add",   func_add,   0, NULL);
  cdoOperatorAdd("sub",   func_sub,   0, NULL);
  cdoOperatorAdd("mul",   func_mul,   0, NULL);
  cdoOperatorAdd("div",   func_div,   0, NULL);
  cdoOperatorAdd("min",   func_min,   0, NULL);
  cdoOperatorAdd("max",   func_max,   0, NULL);
  cdoOperatorAdd("atan2", func_atan2, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61

62
63
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
64
  operatorCheckArgc(0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65

66
67
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
  int streamID2 = pstreamOpenRead(cdoStreamName(1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68

69
70
  int streamIDx1 = streamID1;
  int streamIDx2 = streamID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71

72
73
74
  field_type field1, field2;
  field_type *fieldx1 = &field1;
  field_type *fieldx2 = &field2;
75

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

81
  if ( cdoVerbose ) vlistPrint(vlistID1);
82
  if ( cdoVerbose ) vlistPrint(vlistID2);
83

84
85
86
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = vlistInqTaxis(vlistID2);
  int taxisIDx1 = taxisID1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87

88
89
  int ntsteps1 = vlistNtsteps(vlistID1);
  int ntsteps2 = vlistNtsteps(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
91
  if ( ntsteps1 == 0 ) ntsteps1 = 1;
  if ( ntsteps2 == 0 ) ntsteps2 = 1;
92

93
  bool lfill1, lfill2;
94
95
96
97
98
99
100
101
102
103
104
105
  if ( vlistNvars(vlistID1) == 1 && vlistNvars(vlistID2) == 1 )
    {
      lfill2 = vlistNrecs(vlistID1) != 1 && vlistNrecs(vlistID2) == 1;
      lfill1 = vlistNrecs(vlistID1) == 1 && vlistNrecs(vlistID2) != 1;
    }
  else
    {
      lfill2 = vlistNvars(vlistID1) != 1 && vlistNvars(vlistID2) == 1;
      lfill1 = vlistNvars(vlistID1) == 1 && vlistNvars(vlistID2) != 1;
    }

  if ( lfill2 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
    {
107
108
      nlevels2 = vlistCompareX(vlistID1, vlistID2, CMP_DIM);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
      if ( ntsteps1 != 1 && ntsteps2 == 1 )
110
	{
111
	  filltype = FILL_VAR;
112
	  cdoPrint("Filling up stream2 >%s< by copying the first variable.", cdoStreamName(1)->args);
113
114
115
	}
      else
	{
116
	  filltype = FILL_VARTS;
117
	  cdoPrint("Filling up stream2 >%s< by copying the first variable of each timestep.", cdoStreamName(1)->args);
118
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
    }
120
  else if ( lfill1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
    {
122
123
      nlevels2 = vlistCompareX(vlistID2, vlistID1, CMP_DIM);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
      if ( ntsteps1 == 1 && ntsteps2 != 1 )
125
	{
126
	  filltype = FILL_VAR;
127
	  cdoPrint("Filling up stream1 >%s< by copying the first variable.", cdoStreamName(0)->args);
128
129
130
	}
      else
	{
131
	  filltype = FILL_VARTS;
132
	  cdoPrint("Filling up stream1 >%s< by copying the first variable of each timestep.", cdoStreamName(0)->args);
133
	}
134
135
136
137
138
139
140
      streamIDx1 = streamID2;
      streamIDx2 = streamID1;
      vlistIDx1 = vlistID2;
      vlistIDx2 = vlistID1;
      taxisIDx1 = taxisID2;
      fieldx1 = &field2;
      fieldx2 = &field1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
    }

143
  if ( filltype == FILL_NONE ) vlistCompare(vlistID1, vlistID2, CMP_ALL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144

145
  int gridsize = vlistGridsizeMax(vlistIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146

147
148
  field_init(&field1);
  field_init(&field2);
149
150
  field1.ptr = (double*) Malloc(gridsize*sizeof(double));
  field2.ptr = (double*) Malloc(gridsize*sizeof(double));
151
152
  if ( filltype == FILL_VAR || filltype == FILL_VARTS )
    {
153
154
      vardata2 = (double*) Malloc(gridsize*nlevels2*sizeof(double));
      varnmiss2 = (int*) Malloc(nlevels2*sizeof(int));
155
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156

157
  if ( cdoVerbose ) cdoPrint("Number of timesteps: file1 %d, file2 %d", ntsteps1, ntsteps2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
159
160

  if ( filltype == FILL_NONE )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
      if ( ntsteps1 != 1 && ntsteps2 == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
	{
163
	  filltype = FILL_TS;
164
	  cdoPrint("Filling up stream2 >%s< by copying the first timestep.", cdoStreamName(1)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
      else if ( ntsteps1 == 1 && ntsteps2 != 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
	{
168
	  filltype = FILL_TS;
169
	  cdoPrint("Filling up stream1 >%s< by copying the first timestep.", cdoStreamName(0)->args);
170
171
172
173
174
175
176
	  streamIDx1 = streamID2;
          streamIDx2 = streamID1;
	  vlistIDx1 = vlistID2;
	  vlistIDx2 = vlistID1;
	  taxisIDx1 = taxisID2;
	  fieldx1 = &field2;
	  fieldx2 = &field1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
178
	}

179
      if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
	{
181
	  nvars  = vlistNvars(vlistIDx2);
182
183
	  vardata  = (double **) Malloc(nvars*sizeof(double *));
	  varnmiss = (int **) Malloc(nvars*sizeof(int *));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
185
	  for ( varID = 0; varID < nvars; varID++ )
	    {
186
187
	      int gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
	      int nlev     = zaxisInqSize(vlistInqVarZaxis(vlistIDx2, varID));
188
189
	      vardata[varID]  = (double*) Malloc(nlev*gridsize*sizeof(double));
	      varnmiss[varID] = (int*) Malloc(nlev*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
192
193
	    }
	}
    }

194
  int vlistID3 = vlistDuplicate(vlistIDx1);
195
196
197
198
199
200
  if ( filltype == FILL_TS && vlistIDx1 != vlistID1 )
    {
      nvars  = vlistNvars(vlistID1);
      for ( varID = 0; varID < nvars; varID++ )
	vlistDefVarMissval(vlistID3, varID, vlistInqVarMissval(vlistID1, varID));
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201

202
  int taxisID3 = taxisDuplicate(taxisIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
  vlistDefTaxis(vlistID3, taxisID3);

205
206
  int streamID3 = pstreamOpenWrite(cdoStreamName(2), cdoFiletype());
  pstreamDefVlist(streamID3, vlistID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207

208
209
  int tsID = 0;
  int tsID2 = 0;
210
  while ( (nrecs = pstreamInqTimestep(streamIDx1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
    {
212
      if ( tsID == 0 || filltype == FILL_NONE || filltype == FILL_FILE || filltype == FILL_VARTS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
	{
214
	  nrecs2 = pstreamInqTimestep(streamIDx2, tsID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
	  if ( nrecs2 == 0 )
216
217
218
219
	    {
	      if ( filltype == FILL_NONE && streamIDx2 == streamID2 )
		{
		  filltype = FILL_FILE;
220
		  cdoPrint("Filling up stream2 >%s< by copying all timesteps.", cdoStreamName(1)->args);
221
222
223
224
225
		}

	      if ( filltype == FILL_FILE )
		{
		  tsID2 = 0;
226
227
		  pstreamClose(streamID2);
		  streamID2 = pstreamOpenRead(cdoStreamName(1));
228
229
		  streamIDx2 = streamID2;

230
		  vlistID2 = pstreamInqVlist(streamID2);
231
232
		  vlistIDx2 = vlistID2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
		  vlistCompare(vlistID1, vlistID2, CMP_DIM);
234

235
		  nrecs2 = pstreamInqTimestep(streamIDx2, tsID2);
236
		  if ( nrecs2 == 0 )
237
		    cdoAbort("Empty input stream %s!", cdoStreamName(1)->args);
238
239
240
241
		}
	      else
		cdoAbort("Input streams have different number of timesteps!");
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
	}

244
      taxisCopyTimestep(taxisID3, taxisIDx1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245

246
      pstreamDefTimestep(streamID3, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247

248
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
	{
250
251
	  pstreamInqRecord(streamIDx1, &varID, &levelID);
	  pstreamReadRecord(streamIDx1, fieldx1->ptr, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
          fieldx1->nmiss = (size_t) nmiss;
253
          int varID2 = varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
          
255
	  if ( tsID == 0 || filltype == FILL_NONE || filltype == FILL_FILE || filltype == FILL_VARTS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
256
	    {
257
258
259
	      int lstatus = nlevels2 > 1 ? varID == 0 : recID == 0;

	      if ( lstatus || (filltype != FILL_VAR && filltype != FILL_VARTS) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
		{
261
262
		  pstreamInqRecord(streamIDx2, &varID2, &levelID2);
		  pstreamReadRecord(streamIDx2, fieldx2->ptr, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
                  fieldx2->nmiss = (size_t) nmiss;
264
265
                  if ( varID   != varID2 ) cdoAbort("Internal error, varIDs of input streams differ!");
                  if ( levelID != levelID2 ) cdoAbort("Internal error, levelIDs of input streams differ!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
		}

268
	      if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
		{
270
271
		  int gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
		  int offset   = gridsize*levelID;
272
273
		  memcpy(vardata[varID]+offset, fieldx2->ptr, gridsize*sizeof(double));
		  varnmiss[varID][levelID] = fieldx2->nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274
		}
275
	      else if ( lstatus && (filltype == FILL_VAR || filltype == FILL_VARTS) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276
		{
277
278
		  int gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
		  int offset   = gridsize*levelID2;
279
280
		  memcpy(vardata2+offset, fieldx2->ptr, gridsize*sizeof(double));
		  varnmiss2[levelID2] = fieldx2->nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
282
		}
	    }
283
	  else if ( filltype == FILL_TS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
	    {
285
286
	      int gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, varID2));
	      int offset   = gridsize*levelID;
287
288
	      memcpy(fieldx2->ptr, vardata[varID]+offset, gridsize*sizeof(double));
	      fieldx2->nmiss = varnmiss[varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
	    }

291
292
	  fieldx1->grid    = vlistInqVarGrid(vlistIDx1, varID);
	  fieldx1->missval = vlistInqVarMissval(vlistIDx1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293

294
	  if ( filltype == FILL_VAR || filltype == FILL_VARTS )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
	    {
296
	      int gridsize = gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
297
298
	      levelID2 = 0;
	      if ( nlevels2 > 1 ) levelID2 = levelID;
299
	      int offset   = gridsize*levelID2;
300
301
	      memcpy(fieldx2->ptr, vardata2+offset, gridsize*sizeof(double));
	      fieldx2->nmiss   = varnmiss2[levelID2];
302
303
	      fieldx2->grid    = vlistInqVarGrid(vlistIDx2, 0);
	      fieldx2->missval = vlistInqVarMissval(vlistIDx2, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
305
306
	    }
	  else
	    {
307
308
	      fieldx2->grid    = vlistInqVarGrid(vlistIDx2, varID2);
	      fieldx2->missval = vlistInqVarMissval(vlistIDx2, varID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312
	    }

	  farfun(&field1, field2, operfunc);

313
314
	  pstreamDefRecord(streamID3, varID, levelID);
	  pstreamWriteRecord(streamID3, field1.ptr, (int)field1.nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
	}

      tsID++;
318
      tsID2++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
320
    }

321
322
323
  pstreamClose(streamID3);
  pstreamClose(streamID2);
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324

Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
326
  vlistDestroy(vlistID3);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
328
329
330
  if ( vardata )
    {
      for ( varID = 0; varID < nvars; varID++ )
	{
331
332
	  Free(vardata[varID]);
	  Free(varnmiss[varID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
334
	}

335
336
      Free(vardata);
      Free(varnmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
338
    }

339
340
341
342
  if ( field1.ptr ) Free(field1.ptr);
  if ( field2.ptr ) Free(field2.ptr);
  if ( vardata2   ) Free(vardata2);
  if ( varnmiss2  ) Free(varnmiss2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
344
345

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
346
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
}