Arithc.cc 4.07 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
18
  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.
*/

/*
19
   This module contains the following operators:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20

21
22
23
24
      Arithc     addc            Add by constant
      Arithc     subc            Subtract by constant
      Arithc     mulc            Multiply by constant
      Arithc     divc            Divide by constant
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
      Arithc     mod             Modulo operator
26
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
28


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


Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
37
38
int *fill_vars(int vlistID)
{
  int varID;
  int nvars = vlistNvars(vlistID);
39
  int *vars = (int*) Malloc(nvars*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

  if ( cdoNumVarnames )
    {
      int nfound = 0;
      char varname[CDI_MAX_NAME];

      for ( varID = 0; varID < nvars; ++varID )
	{
	  vars[varID] = 0;

	  vlistInqVarName(vlistID, varID, varname);

	  for ( int i = 0; i < cdoNumVarnames; ++i )
	    if ( strcmp(varname, cdoVarnames[i]) == 0 )
	      {
		vars[varID] = 1;
		nfound++;
		break;
	      }
	}

      if ( nfound == 0 )
	cdoAbort("Variable -n %s%s not found!", cdoVarnames[0], cdoNumVarnames > 1 ? ",..." : "");
    }
  else
    {
      for ( varID = 0; varID < nvars; ++varID ) vars[varID] = 1;
    }

69
  return vars;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
72
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
74
void *Arithc(void *argument)
{
75
  int nmiss;
76
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79
80
  int varID, levelID;

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
84
85
86
  cdoOperatorAdd("addc", func_add, 0, "constant value");
  cdoOperatorAdd("subc", func_sub, 0, "constant value");
  cdoOperatorAdd("mulc", func_mul, 0, "constant value");
  cdoOperatorAdd("divc", func_div, 0, "constant value");
  cdoOperatorAdd("mod",  func_mod, 0, "divisor");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88

89
90
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91

Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
  operatorInputArg(cdoOperatorEnter(operatorID));
93
  operatorCheckArgc(1);
94
  double rconst = parameter2double(operatorArgv()[0]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

96
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97

98
  int vlistID1 = pstreamInqVlist(streamID1);
99
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100

101
  int *vars = fill_vars(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102

103
104
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
  vlistDefTaxis(vlistID2, taxisID2);

107
108
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109

110
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111

112
  field_type field;
113
  field_init(&field);
114
  field.ptr    = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
116
  field.weight = NULL;

117
  int tsID = 0;
118
  while ( (nrecs = pstreamInqTimestep(streamID1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
121
    {
      taxisCopyTimestep(taxisID2, taxisID1);

122
      pstreamDefTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123

124
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
	{
126
127
	  pstreamInqRecord(streamID1, &varID, &levelID);
	  pstreamReadRecord(streamID1, field.ptr, &nmiss);
128
          field.nmiss = (size_t) nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129

Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132
133
	  if ( vars[varID] )
	    {
	      field.grid    = vlistInqVarGrid(vlistID1, varID);
	      field.missval = vlistInqVarMissval(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134

Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
	      farcfun(&field, rconst, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136

Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
139
	      /* recalculate number of missing values */
	      gridsize = gridInqSize(field.grid);
	      field.nmiss = 0;
140
	      for ( int i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
		if ( DBL_IS_EQUAL(field.ptr[i], field.missval) ) field.nmiss++;
	    }
143

144
          nmiss = (int) field.nmiss;
145
146
	  pstreamDefRecord(streamID2, varID, levelID);
	  pstreamWriteRecord(streamID2, field.ptr, nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
148
149
150
	}
      tsID++;
    }

151
152
  pstreamClose(streamID2);
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153

154
155
  if ( field.ptr ) Free(field.ptr);
  if ( vars ) Free(vars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156

157
158
  vlistDestroy(vlistID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160
  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
}