expr.cc 71 KB
Newer Older
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-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
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
#include <cstdlib>
Oliver Heidmann's avatar
Oliver Heidmann committed
19
#include <cassert>
20
#include <cdi.h>
21

22
#include "cdo_options.h"
23
#include "dmemory.h"
24
#include "process_int.h"
25
#include "field.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
#include "expr.h"
27
#include "expr_fun.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
#include "expr_yacc.h"
29
#include "cdo_zaxis.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
33
34
35
static const char *const ExIn[] = { "expr", "init" };
static const char *const tmpvnm = "_tmp_";
static int pointID = -1;
static int zonalID = -1;
static int surfaceID = -1;
36

37
38
39
40
41
enum
{
  FT_STD,
  FT_CONST,
  FT_FLD,
42
  FT_ZON,
43
44
45
  FT_VERT,
  FT_COORD,
  FT_1C,
46
  FT_2C,
47
48
49
  FT_0
};

50
// clang-format off
51
52
53
54
55
56
57
58
59
60
static constexpr double COMPLT(const double x, const double y) noexcept { return x < y; }
static constexpr double COMPGT(const double x, const double y) noexcept { return x > y; }
static constexpr double COMPLE(const double x, const double y) noexcept { return x <= y; }
static constexpr double COMPGE(const double x, const double y) noexcept { return x >= y; }
static constexpr double COMPNE(const double x, const double y) noexcept { return IS_NOT_EQUAL(x, y); }
static constexpr double COMPEQ(const double x, const double y) noexcept { return IS_EQUAL(x, y); }
static constexpr double COMPLEG(const double x, const double y) noexcept { return (x < y) ? -1. : (x > y); }
static constexpr double COMPAND(const double x, const double y) noexcept { return IS_NOT_EQUAL(x, 0) && IS_NOT_EQUAL(y, 0); }
static constexpr double COMPOR(const double x, const double y) noexcept { return IS_NOT_EQUAL(x, 0) || IS_NOT_EQUAL(y, 0); }
static constexpr double COMPNOT(const double x) noexcept { return IS_EQUAL(x, 0); }
61
62
63
64
65
66
67
68
69
static inline    double MVCOMPLT(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPLT(x, y); }
static inline    double MVCOMPGT(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPGT(x, y); }
static inline    double MVCOMPLE(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPLE(x, y); }
static inline    double MVCOMPGE(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPGE(x, y); }
static inline    double MVCOMPNE(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPNE(x, y); }
static inline    double MVCOMPEQ(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPEQ(x, y); }
static inline    double MVCOMPLEG(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPLEG(x, y); }
static inline    double MVCOMPAND(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPAND(x, y); }
static inline    double MVCOMPOR(const double x, const double y, const double m) noexcept { return DBL_IS_EQUAL(x, m) ? m : COMPOR(x, y); }
70

71
72
73
static double f_float(const double x) { return (float) (x); }
static double f_int(const double x) { return (int) (x); }
static double f_nint(const double x) { return std::round(x); }
74
static double f_rand(const double x) { (void)x; return ((double) std::rand()) / ((double) RAND_MAX); }
75
76
77
static double f_sqr(const double x) { return x * x; }
static double f_rad(const double x) { return x * M_PI / 180.; }
static double f_deg(const double x) { return x * 180. / M_PI; }
78
static double f_isMissval(const double x) { (void)x; return 0; }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
static double pt_ngp(const paramType *const p) { return p->ngp; }
static double pt_nlev(const paramType *const p) { return p->nlev; }
static double pt_size(const paramType *const p) { return p->ngp * p->nlev; }
static double pt_missval(const paramType *const p) { return p->missval; }
static double ts_ctimestep(const double *const data) { return std::lround(data[CTIMESTEP]); }
static double ts_cdate(const double *const data) { return std::lround(data[CDATE]); }
static double ts_ctime(const double *const data) { return std::lround(data[CTIME]); }
static double ts_cdeltat(const double *const data) { return data[CDELTAT]; }
static double ts_cday(const double *const data) { return data[CDAY]; }
static double ts_cmonth(const double *const data) { return data[CMONTH]; }
static double ts_cyear(const double *const data) { return data[CYEAR]; }
static double ts_csecond(const double *const data) { return data[CSECOND]; }
static double ts_cminute(const double *const data) { return data[CMINUTE]; }
static double ts_chour(const double *const data) { return data[CHOUR]; }
93
// clang-format on
94

Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
struct func_t
96
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
  int type;
98
  int flag;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
  const char *name;    // function name
100
  void (*func)(void);  // pointer to function
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102

Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
static const func_t fun_sym_tbl[] = {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
  // scalar functions
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
  { FT_STD, 0, "abs", (void (*)(void))((double (*)(double)) fabs) },  // math functions could be inlined
  { FT_STD, 0, "floor", (void (*)(void))((double (*)(double)) floor) },
  { FT_STD, 0, "ceil", (void (*)(void))((double (*)(double)) ceil) },
  { FT_STD, 0, "sqrt", (void (*)(void))((double (*)(double)) sqrt) },
  { FT_STD, 0, "exp", (void (*)(void))((double (*)(double)) exp) },
  { FT_STD, 0, "erf", (void (*)(void))((double (*)(double)) erf) },
  { FT_STD, 0, "log", (void (*)(void))((double (*)(double)) log) },
  { FT_STD, 0, "ln", (void (*)(void))((double (*)(double)) log) },
  { FT_STD, 0, "log10", (void (*)(void))((double (*)(double)) log10) },
  { FT_STD, 0, "sin", (void (*)(void))((double (*)(double)) sin) },
  { FT_STD, 0, "cos", (void (*)(void))((double (*)(double)) cos) },
  { FT_STD, 0, "tan", (void (*)(void))((double (*)(double)) tan) },
  { FT_STD, 0, "sinh", (void (*)(void))((double (*)(double)) sinh) },
  { FT_STD, 0, "cosh", (void (*)(void))((double (*)(double)) cosh) },
  { FT_STD, 0, "tanh", (void (*)(void))((double (*)(double)) tanh) },
  { FT_STD, 0, "asin", (void (*)(void))((double (*)(double)) asin) },
  { FT_STD, 0, "acos", (void (*)(void))((double (*)(double)) acos) },
  { FT_STD, 0, "atan", (void (*)(void))((double (*)(double)) atan) },
  { FT_STD, 0, "asinh", (void (*)(void))((double (*)(double)) asinh) },
  { FT_STD, 0, "acosh", (void (*)(void))((double (*)(double)) acosh) },
  { FT_STD, 0, "atanh", (void (*)(void))((double (*)(double)) atanh) },
  { FT_STD, 0, "gamma", (void (*)(void))((double (*)(double)) tgamma) },
  { FT_STD, 0, "float", reinterpret_cast<void (*)(void)>(&f_float) },
  { FT_STD, 0, "int", reinterpret_cast<void (*)(void)>(&f_int) },
  { FT_STD, 0, "nint", reinterpret_cast<void (*)(void)>(&f_nint) },
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
  { FT_STD, 0, "rand", reinterpret_cast<void (*)(void)>(&f_rand) },
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
133
  { FT_STD, 0, "sqr", reinterpret_cast<void (*)(void)>(&f_sqr) },
  { FT_STD, 0, "rad", reinterpret_cast<void (*)(void)>(&f_rad) },
  { FT_STD, 0, "deg", reinterpret_cast<void (*)(void)>(&f_deg) },
134
  { FT_STD, 0, "isMissval", reinterpret_cast<void (*)(void)>(&f_isMissval) },
135
136

  // constant functions
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
139
140
  { FT_CONST, 0, "ngp", reinterpret_cast<void (*)(void)>(&pt_ngp) },          // number of horizontal grid points
  { FT_CONST, 0, "nlev", reinterpret_cast<void (*)(void)>(&pt_nlev) },        // number of vertical levels
  { FT_CONST, 0, "size", reinterpret_cast<void (*)(void)>(&pt_size) },        // ngp*nlev
  { FT_CONST, 0, "missval", reinterpret_cast<void (*)(void)>(&pt_missval) },  // Returns the missing value of a variable
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141

142
  // CDO field functions (Reduce grid to point)
143
144
  { FT_FLD, 0, "fldmin", reinterpret_cast<void (*)(void)>(&fieldMin) },
  { FT_FLD, 0, "fldmax", reinterpret_cast<void (*)(void)>(&fieldMax) },
145
  { FT_FLD, 0, "fldsum", reinterpret_cast<void (*)(void)>(&fieldSum) },
146
147
  { FT_FLD, 1, "fldmean", reinterpret_cast<void (*)(void)>(&fieldMeanw) },
  { FT_FLD, 1, "fldavg", reinterpret_cast<void (*)(void)>(&fieldAvgw) },
148
149
150
151
  { FT_FLD, 1, "fldstd", reinterpret_cast<void (*)(void)>(&fieldStdw) },
  { FT_FLD, 1, "fldstd1", reinterpret_cast<void (*)(void)>(&fieldStd1w) },
  { FT_FLD, 1, "fldvar", reinterpret_cast<void (*)(void)>(&fieldVarw) },
  { FT_FLD, 1, "fldvar1", reinterpret_cast<void (*)(void)>(&fieldVar1w) },
152

153
  // CDO zonal functions (Reduce grid to point)
154
155
156
157
158
159
160
161
162
  { FT_ZON, 0, "zonmin", reinterpret_cast<void (*)(void)>(&fieldZonMin) },
  { FT_ZON, 0, "zonmax", reinterpret_cast<void (*)(void)>(&fieldZonMax) },
  { FT_ZON, 0, "zonsum", reinterpret_cast<void (*)(void)>(&fieldZonSum) },
  { FT_ZON, 0, "zonmean", reinterpret_cast<void (*)(void)>(&fieldZonMean) },
  { FT_ZON, 0, "zonavg", reinterpret_cast<void (*)(void)>(&fieldZonAvg) },
  { FT_ZON, 0, "zonstd", reinterpret_cast<void (*)(void)>(&fieldZonStd) },
  { FT_ZON, 0, "zonstd1", reinterpret_cast<void (*)(void)>(&fieldZonStd1) },
  { FT_ZON, 0, "zonvar", reinterpret_cast<void (*)(void)>(&fieldZonVar) },
  { FT_ZON, 0, "zonvar1", reinterpret_cast<void (*)(void)>(&fieldZonVar1) },
163
164

  // CDO field functions (Reduce level to point)
165
166
  { FT_VERT, 0, "vertmin", reinterpret_cast<void (*)(void)>(&fieldMin) },
  { FT_VERT, 0, "vertmax", reinterpret_cast<void (*)(void)>(&fieldMax) },
167
  { FT_VERT, 0, "vertsum", reinterpret_cast<void (*)(void)>(&fieldSum) },
168
169
  { FT_VERT, 1, "vertmean", reinterpret_cast<void (*)(void)>(&fieldMeanw) },
  { FT_VERT, 1, "vertavg", reinterpret_cast<void (*)(void)>(&fieldAvgw) },
170
171
172
173
  { FT_VERT, 1, "vertstd", reinterpret_cast<void (*)(void)>(&fieldStdw) },
  { FT_VERT, 1, "vertstd1", reinterpret_cast<void (*)(void)>(&fieldStd1w) },
  { FT_VERT, 1, "vertvar", reinterpret_cast<void (*)(void)>(&fieldVarw) },
  { FT_VERT, 1, "vertvar1", reinterpret_cast<void (*)(void)>(&fieldVar1w) },
174

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
176
177
178
179
  { FT_COORD, 0, "clon", nullptr },
  { FT_COORD, 0, "clat", nullptr },
  { FT_COORD, 0, "clev", nullptr },
  { FT_COORD, 0, "gridarea", nullptr },
  { FT_COORD, 0, "gridweight", nullptr },
180

Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
184
185
186
187
188
189
190
  { FT_0, 0, "ctimestep", reinterpret_cast<void (*)(void)>(&ts_ctimestep) },
  { FT_0, 0, "cdate", reinterpret_cast<void (*)(void)>(&ts_cdate) },
  { FT_0, 0, "ctime", reinterpret_cast<void (*)(void)>(&ts_ctime) },
  { FT_0, 0, "cdeltat", reinterpret_cast<void (*)(void)>(&ts_cdeltat) },
  { FT_0, 0, "cday", reinterpret_cast<void (*)(void)>(&ts_cday) },
  { FT_0, 0, "cmonth", reinterpret_cast<void (*)(void)>(&ts_cmonth) },
  { FT_0, 0, "cyear", reinterpret_cast<void (*)(void)>(&ts_cyear) },
  { FT_0, 0, "csecond", reinterpret_cast<void (*)(void)>(&ts_csecond) },
  { FT_0, 0, "cminute", reinterpret_cast<void (*)(void)>(&ts_cminute) },
  { FT_0, 0, "chour", reinterpret_cast<void (*)(void)>(&ts_chour) },
191

Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
193
194
195
196
  { FT_1C, 0, "sellevel", nullptr },
  { FT_1C, 0, "sellevidx", nullptr },
  { FT_2C, 0, "sellevelrange", nullptr },
  { FT_2C, 0, "sellevidxrange", nullptr },
  // {FT_1C, 0, "gridindex", nullptr},
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
198
};

Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
static const int NumFunc = sizeof(fun_sym_tbl) / sizeof(fun_sym_tbl[0]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200

201
202
static void
node_data_delete(nodeType *p)
203
{
204
  if (p)
205
    {
206
207
208
      if (p->param.data)
        {
          Free(p->param.data);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
          p->param.data = nullptr;
210
        }
211
212
213
    }
}

214
215
static void
node_delete(nodeType *p)
216
{
217
  if (p)
218
    {
219
      if (p->type == typeVar) node_data_delete(p);
220
221
222
223
      Free(p);
    }
}

224
225
static int
get_funcID(const char *fun)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
{
  int funcID = -1;
228
  for (int i = 0; i < NumFunc; i++)
229
    if (cdo_cmpstr(fun, fun_sym_tbl[i].name))
230
231
232
      {
        funcID = i;
        break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
234
      }

235
  if (funcID == -1) cdoAbort("Function >%s< not available!", fun);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
237
238
239

  return funcID;
}

240
static constexpr bool
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
isCompare(const int oper) noexcept
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
  return oper == LEG || oper == GE || oper == LE || oper == EQ || oper == NE || oper == GT || oper == LT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
}

246
static void
247
param_meta_copy(paramType &out, const paramType &in)
248
{
249
250
251
252
253
254
255
256
257
258
259
260
  out.type = in.type;
  out.gridID = in.gridID;
  out.zaxisID = in.zaxisID;
  out.datatype = in.datatype;
  out.steptype = in.steptype;
  out.ngp = in.ngp;
  out.nlat = in.nlat;
  out.nlev = in.nlev;
  out.missval = in.missval;
  out.nmiss = 0;
  out.coord = 0;
  out.lmiss = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261
262
263
264
  out.name = nullptr;
  out.longname = nullptr;
  out.units = nullptr;
  out.data = nullptr;
265
266
}

267
static nodeType *
268
expr_con_con(const int oper, const nodeType *p1, const nodeType *p2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
{
270
  nodeType *p = (nodeType *) Calloc(1, sizeof(nodeType));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271

272
  p->type = typeCon;
273
  p->ltmpobj = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
  double cval1 = p1->u.con.value;
276
  const double cval2 = p2->u.con.value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277

278
  // clang-format off
279
  switch (oper)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280
    {
281
282
283
284
285
286
287
288
289
    case LT:   cval1 = COMPLT(cval1, cval2); break;
    case GT:   cval1 = COMPGT(cval1, cval2); break;
    case LE:   cval1 = COMPLE(cval1, cval2); break;
    case GE:   cval1 = COMPGE(cval1, cval2); break;
    case NE:   cval1 = COMPNE(cval1, cval2); break;
    case EQ:   cval1 = COMPEQ(cval1, cval2); break;
    case LEG:  cval1 = COMPLEG(cval1, cval2); break;
    case AND:  cval1 = COMPAND(cval1, cval2); break;
    case OR:   cval1 = COMPOR(cval1, cval2); break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
291
292
293
    case '+':  cval1 = cval1 + cval2; break;
    case '-':  cval1 = cval1 - cval2; break;
    case '*':  cval1 = cval1 * cval2; break;
    case '/':  cval1 = cval1 / cval2; break;
294
    case '^':  cval1 = std::pow(cval1, cval2); break;
295
    default:   cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
    }
297
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298

299
300
  p->u.con.value = cval1;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
  return p;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
303
}

304
static void
305
306
oper_expr_con_var(const int oper, const bool nmiss, const size_t n, const double missval1, const double missval2,
                  double *restrict odat, const double cval, const double *restrict idat)
307
{
308
  size_t i;
309

310
  // clang-format off
311
312
  if (nmiss)
    {
313
      switch (oper)
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
        {
        case LT:  for (i=0; i<n; ++i) odat[i] = MVCOMPLT(cval, idat[i], missval1); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = MVCOMPGT(cval, idat[i], missval1); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = MVCOMPLE(cval, idat[i], missval1); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = MVCOMPGE(cval, idat[i], missval1); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = MVCOMPNE(cval, idat[i], missval1); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = MVCOMPEQ(cval, idat[i], missval1); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = MVCOMPLEG(cval, idat[i], missval1); break;
        case AND: for (i=0; i<n; ++i) odat[i] = MVCOMPAND(cval, idat[i], missval1); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = MVCOMPOR(cval, idat[i], missval1); break;
        case '^': for (i=0; i<n; ++i) odat[i] = POWMN(cval, idat[i]); break;
        case '+': for (i=0; i<n; ++i) odat[i] = ADDMN(cval, idat[i]); break;
        case '-': for (i=0; i<n; ++i) odat[i] = SUBMN(cval, idat[i]); break;
        case '*': for (i=0; i<n; ++i) odat[i] = MULMN(cval, idat[i]); break;
        case '/': for (i=0; i<n; ++i) odat[i] = DIVMN(cval, idat[i]); break;
329
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
330
331
332
        }
    }
  else
333
    {
334
      switch (oper)
335
336
337
338
339
340
341
342
343
344
        {
        case LT:  for (i=0; i<n; ++i) odat[i] = COMPLT(cval, idat[i]); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = COMPGT(cval, idat[i]); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = COMPLE(cval, idat[i]); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = COMPGE(cval, idat[i]); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = COMPNE(cval, idat[i]); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = COMPEQ(cval, idat[i]); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = COMPLEG(cval, idat[i]); break;
        case AND: for (i=0; i<n; ++i) odat[i] = COMPAND(cval, idat[i]); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = COMPOR(cval, idat[i]); break;
345
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(cval, idat[i]); break;
346
347
348
349
        case '+': for (i=0; i<n; ++i) odat[i] = cval + idat[i]; break;
        case '-': for (i=0; i<n; ++i) odat[i] = cval - idat[i]; break;
        case '*': for (i=0; i<n; ++i) odat[i] = cval * idat[i]; break;
        case '/': for (i=0; i<n; ++i) odat[i] = DIVMN(cval, idat[i]); break;
350
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
351
        }
352
    }
353
  // clang-format on
354
355
}

356
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
357
oper_expr_var_con(const int oper, const bool lmiss, const size_t n, const double missval1, const double missval2,
358
                  double *restrict odat, const double *restrict idat, const double cval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
{
360
  size_t i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361

362
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363
  if (lmiss)
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
    {
      switch (oper)
        {
        case LT:  for (i=0; i<n; ++i) odat[i] = MVCOMPLT(idat[i], cval, missval1); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = MVCOMPGT(idat[i], cval, missval1); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = MVCOMPLE(idat[i], cval, missval1); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = MVCOMPGE(idat[i], cval, missval1); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = MVCOMPNE(idat[i], cval, missval1); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = MVCOMPEQ(idat[i], cval, missval1); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = MVCOMPLEG(idat[i], cval, missval1); break;
        case AND: for (i=0; i<n; ++i) odat[i] = MVCOMPAND(idat[i], cval, missval1); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = MVCOMPOR(idat[i], cval, missval1); break;
        case '^': for (i=0; i<n; ++i) odat[i] = POWMN(idat[i], cval); break;
        case '+': for (i=0; i<n; ++i) odat[i] = ADDMN(idat[i], cval); break;
        case '-': for (i=0; i<n; ++i) odat[i] = SUBMN(idat[i], cval); break;
        case '*': for (i=0; i<n; ++i) odat[i] = MULMN(idat[i], cval); break;
        case '/': for (i=0; i<n; ++i) odat[i] = DIVMN(idat[i], cval); break;
381
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
382
383
384
        }
    }
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
385
    {
386
387
388
389
390
391
392
393
394
395
396
      switch (oper)
        {
        case LT:  for (i=0; i<n; ++i) odat[i] = COMPLT(idat[i], cval); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = COMPGT(idat[i], cval); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = COMPLE(idat[i], cval); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = COMPGE(idat[i], cval); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = COMPNE(idat[i], cval); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = COMPEQ(idat[i], cval); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = COMPLEG(idat[i], cval); break;
        case AND: for (i=0; i<n; ++i) odat[i] = COMPAND(idat[i], cval); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = COMPOR(idat[i], cval); break;
397
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(idat[i], cval); break;
398
399
400
401
402
403
404
        case '+': for (i=0; i<n; ++i) odat[i] = idat[i] + cval; break;
        case '-': for (i=0; i<n; ++i) odat[i] = idat[i] - cval; break;
        case '*': for (i=0; i<n; ++i) odat[i] = idat[i] * cval; break;
        case '/':
          if (IS_EQUAL(cval, 0)) for (i=0; i<n; ++i) odat[i] = missval1;
          else                   for (i=0; i<n; ++i) odat[i] = idat[i] / cval;
          break;
405
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
406
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407
    }
408
  // clang-format on
409
410
}

411
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
oper_expr_var_var(const int oper, const bool lmiss, const size_t n, const double missval1, const double missval2,
413
                  double *restrict odat, const double *restrict idat1, const double *restrict idat2)
414
{
415
  size_t i;
416

417
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
  if (lmiss)
419
    {
420
      switch (oper)
421
        {
422
423
424
425
426
427
428
429
430
431
432
433
434
435
        case LT:  for (i=0; i<n; ++i) odat[i] = MVCOMPLT(idat1[i], idat2[i], missval1); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = MVCOMPGT(idat1[i], idat2[i], missval1); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = MVCOMPLE(idat1[i], idat2[i], missval1); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = MVCOMPGE(idat1[i], idat2[i], missval1); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = MVCOMPNE(idat1[i], idat2[i], missval1); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = MVCOMPEQ(idat1[i], idat2[i], missval1); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = MVCOMPLEG(idat1[i], idat2[i], missval1); break;
        case AND: for (i=0; i<n; ++i) odat[i] = MVCOMPAND(idat1[i], idat2[i], missval1); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = MVCOMPOR(idat1[i], idat2[i], missval1); break;
        case '^': for (i=0; i<n; ++i) odat[i] = POWMN(idat1[i], idat2[i]); break;
        case '+': for (i=0; i<n; ++i) odat[i] = ADDMN(idat1[i], idat2[i]); break;
        case '-': for (i=0; i<n; ++i) odat[i] = SUBMN(idat1[i], idat2[i]); break;
        case '*': for (i=0; i<n; ++i) odat[i] = MULMN(idat1[i], idat2[i]); break;
        case '/': for (i=0; i<n; ++i) odat[i] = DIVMN(idat1[i], idat2[i]); break;
436
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
        }
    }
  else
    {
      switch (oper)
        {
        case LT:  for (i=0; i<n; ++i) odat[i] = COMPLT(idat1[i], idat2[i]); break;
        case GT:  for (i=0; i<n; ++i) odat[i] = COMPGT(idat1[i], idat2[i]); break;
        case LE:  for (i=0; i<n; ++i) odat[i] = COMPLE(idat1[i], idat2[i]); break;
        case GE:  for (i=0; i<n; ++i) odat[i] = COMPGE(idat1[i], idat2[i]); break;
        case NE:  for (i=0; i<n; ++i) odat[i] = COMPNE(idat1[i], idat2[i]); break;
        case EQ:  for (i=0; i<n; ++i) odat[i] = COMPEQ(idat1[i], idat2[i]); break;
        case LEG: for (i=0; i<n; ++i) odat[i] = COMPLEG(idat1[i], idat2[i]); break;
        case AND: for (i=0; i<n; ++i) odat[i] = COMPAND(idat1[i], idat2[i]); break;
        case OR:  for (i=0; i<n; ++i) odat[i] = COMPOR(idat1[i], idat2[i]); break;
452
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(idat1[i], idat2[i]); break;
453
454
455
        case '+': for (i=0; i<n; ++i) odat[i] = idat1[i] + idat2[i]; break;
        case '-': for (i=0; i<n; ++i) odat[i] = idat1[i] - idat2[i]; break;
        case '*': for (i=0; i<n; ++i) odat[i] = idat1[i] * idat2[i]; break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
        case '/': for (i=0; i<n; ++i) odat[i] = IS_EQUAL(idat2[i], 0.) ? missval1 : idat1[i] / idat2[i]; break;
457
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
458
459
        }
    }
460
  // clang-format on
461
462
}

463
static nodeType *
464
expr_con_var(const int init, const int oper, const nodeType *p1, const nodeType *p2)
465
{
466
467
468
469
470
471
  const size_t ngp = p2->param.ngp;
  const size_t nlev = p2->param.nlev;
  const size_t nmiss = p2->param.nmiss;
  const int datatype = p2->param.datatype;
  const double missval1 = p2->param.missval;
  const double missval2 = p2->param.missval;
472

473
  const size_t n = ngp * nlev;
474

475
  nodeType *p = (nodeType *) Calloc(1, sizeof(nodeType));
476

477
478
  p->type = typeVar;
  p->ltmpobj = true;
479
  p->u.var.nm = strdup(tmpvnm);
480
  param_meta_copy(p->param, p2->param);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
  p->param.name = p->u.var.nm;
482

483
  if (!init)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484
    {
485
      p->param.data = (double *) Malloc(n * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
487
488
      double *restrict odat = p->param.data;
      const double *restrict idat = p2->param.data;
      double cval = p1->u.con.value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;
490

Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
      oper_expr_con_var(oper, nmiss, n, missval1, missval2, odat, cval, idat);
492

Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
      p->param.nmiss = arrayNumMV(n, odat, missval1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
    }
495

Uwe Schulzweida's avatar
Uwe Schulzweida committed
496
  return p;
497
498
}

499
static nodeType *
500
expr_var_con(const int init, const int oper, const nodeType *p1, const nodeType *p2)
501
{
502
503
504
505
506
507
  const size_t ngp = p1->param.ngp;
  const size_t nlev = p1->param.nlev;
  const size_t nmiss = p1->param.nmiss;
  const int datatype = p1->param.datatype;
  const double missval1 = p1->param.missval;
  const double missval2 = p1->param.missval;
508

509
  const size_t n = ngp * nlev;
510

511
  nodeType *p = (nodeType *) Calloc(1, sizeof(nodeType));
512

513
514
  p->type = typeVar;
  p->ltmpobj = true;
515
  p->u.var.nm = strdup(tmpvnm);
516
  param_meta_copy(p->param, p1->param);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
  p->param.name = p->u.var.nm;
518

519
  if (!init)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
    {
521
      p->param.data = (double *) Malloc(n * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
522
523
524
      double *restrict odat = p->param.data;
      const double *restrict idat = p1->param.data;
      double cval = p2->u.con.value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526

Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
      oper_expr_var_con(oper, nmiss, n, missval1, missval2, odat, idat, cval);
528

Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
      p->param.nmiss = arrayNumMV(n, odat, missval1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
    }
531

Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
  return p;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
534
}

535
536
static nodeType *
expr_var_var(int init, int oper, nodeType *p1, nodeType *p2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
{
538
  nodeType *px = p1;
539
540
541
542
  const size_t nmiss1 = p1->param.nmiss;
  const size_t nmiss2 = p2->param.nmiss;
  const double missval1 = p1->param.missval;
  const double missval2 = p2->param.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543

544
545
  const size_t ngp1 = p1->param.ngp;
  const size_t ngp2 = p2->param.ngp;
546
  size_t ngp = ngp1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547

548
549
  const size_t nlat1 = p1->param.nlat;
  const size_t nlat2 = p2->param.nlat;
550
551
  bool lzonal = false;

552
  if (ngp1 != ngp2)
553
    {
554
      if (ngp1 == 1 || ngp2 == 1)
555
        {
556
557
558
559
560
          if (ngp1 == 1)
            {
              ngp = ngp2;
              px = p2;
            }
561
        }
562
563
564
565
      else if (nlat1 == nlat2 && ngp1 > ngp2)
        {
          lzonal = true;
        }
566
      else
567
        {
568
          cdoAbort("%s: Number of grid points differ (%s[%zu] <-> %s[%zu])", __func__, p1->param.name, ngp1, p2->param.name, ngp2);
569
570
571
        }
    }

572
573
  const size_t nlev1 = p1->param.nlev;
  const size_t nlev2 = p2->param.nlev;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574

575
  size_t nlev = nlev1;
576
  if (nlev1 != nlev2)
577
    {
578
      if (nlev1 == 1 || nlev2 == 1)
579
        {
580
581
582
583
584
          if (nlev1 == 1)
            {
              nlev = nlev2;
              px = p2;
            }
585
        }
586
      else
587
        {
588
          cdoAbort("%s: Number of levels differ (%s[%zu] <-> %s[%zu])", __func__, p1->param.name, nlev1, p2->param.name, nlev2);
589
590
591
        }
    }

592
  nodeType *p = (nodeType *) Calloc(1, sizeof(nodeType));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
593

594
595
  p->type = typeVar;
  p->ltmpobj = true;
596
  p->u.var.nm = strdup(tmpvnm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
597

598
  param_meta_copy(p->param, px->param);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
599

600
  if (p->param.steptype == TIME_CONSTANT)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
601
    {
602
603
      const int steptype1 = p1->param.steptype;
      const int steptype2 = p2->param.steptype;
604
      if (steptype1 != TIME_CONSTANT)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
        p->param.steptype = steptype1;
606
      else if (steptype2 != TIME_CONSTANT)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
609
        p->param.steptype = steptype2;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
610
  p->param.name = p->u.var.nm;
611
  // printf("%s %s nmiss %zu %zu\n", p->u.var.nm, px->param.name, nmiss1, nmiss2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612

613
  if (!init)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614
    {
615
      p->param.data = (double *) Malloc(ngp * nlev * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
616

617
      for (size_t k = 0; k < nlev; k++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
        {
619
          size_t loff = k * ngp;
620
621
          const size_t loff1 = (nlev1 > 1) ? k * ngp1 : 0;
          const size_t loff2 = (nlev2 > 1) ? k * ngp2 : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622

623
624
625
          const double *restrict idat1 = p1->param.data + loff1;
          const double *restrict idat2 = p2->param.data + loff2;
          double *restrict odat = p->param.data + loff;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
          const size_t nmiss = nmiss1 || nmiss2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
627

628
          if (ngp1 != ngp2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
629
            {
630
631
              if (lzonal)
                {
632
                  const size_t nlon = ngp1 / nlat1;
633
                  for (size_t j = 0; j < nlat1; ++j)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
634
                    oper_expr_var_con(oper, nmiss, nlon, missval1, missval2, odat + j * nlon, idat1 + j * nlon, idat2[j]);
635
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
636
              else
637
638
639
640
641
642
                {
                  if (ngp2 == 1)
                    oper_expr_var_con(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2[0]);
                  else
                    oper_expr_con_var(oper, nmiss, ngp, missval1, missval2, odat, idat1[0], idat2);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643
            }
644
          else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
645
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
646
              oper_expr_var_var(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
647
            }
648
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
649

650
      p->param.nmiss = arrayNumMV(ngp * nlev, p->param.data, missval1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
651
    }
652

Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
  return p;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
655
}

656
static void
657
ex_copy_var(int init, nodeType *p2, const nodeType *p1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
658
{
659
  if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
660
661
    cdoPrint("\t%s\tcopy\t%s[L%zu][N%zu] = %s[L%zu][N%zu]", ExIn[init], p2->param.name, p2->param.nlev, p2->param.ngp,
             p1->param.name, p2->param.nlev, p2->param.ngp);
662

663
  const size_t ngp = p1->param.ngp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
664
665
  assert(ngp > 0);

666
  if (ngp != p2->param.ngp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
667
    cdoAbort("%s: Number of grid points differ (%s[%zu] = %s[%zu])", __func__, p2->param.name, p2->param.ngp, p1->param.name, ngp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
668

669
  const size_t nlev = p1->param.nlev;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
670
  assert(nlev > 0);
671

672
  if (nlev != p2->param.nlev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
673
    cdoAbort("%s: Number of levels differ (%s[%zu] = %s[%zu])", __func__, p2->param.name, p2->param.nlev, p1->param.name, nlev);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674

675
  if (!init)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
    {
677
      arrayCopy(ngp * nlev, p1->param.data, p2->param.data);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
678
      p2->param.missval = p1->param.missval;
679
      p2->param.nmiss = p1->param.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
680
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
681
682
}

683
static void
684
ex_copy_con(int init, nodeType *p2, const nodeType *p1)
685
{
686
  const double cval = p1->u.con.value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
687

Uwe Schulzweida's avatar
Uwe Schulzweida committed
688
689
  if (Options::cdoVerbose)
    cdoPrint("\t%s\tcopy\t%s[L%zu][N%zu] = %g", ExIn[init], p2->param.name, p2->param.nlev, p2->param.ngp, cval);
690

691
  const size_t ngp = p2->param.ngp;
692
693
  assert(ngp > 0);

694
  const size_t nlev = p2->param.nlev;
695
696
  assert(nlev > 0);

697
  if (!init)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
698
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699
      assert(p2->param.data != nullptr);
700
      varrayFill(ngp * nlev, p2->param.data, cval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
701
    }
702
703
}

704
static void
705
ex_copy(int init, nodeType *p2, const nodeType *p1)
706
{
707
  if (p1->type == typeCon)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
708
    ex_copy_con(init, p2, p1);
709
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
710
    ex_copy_var(init, p2, p1);
711
712
}

713
static nodeType *
714
expr(const int init, const int oper, nodeType *p1, nodeType *p2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
715
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
716
  if (p1 == nullptr || p2 == nullptr) return nullptr;
717

Uwe Schulzweida's avatar
Uwe Schulzweida committed
718
  const char *coper = "???";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
719

720
  if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
    {
722
      switch (oper)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
723
        {
724
725
726
727
728
729
730
731
732
        case LT: coper = "<"; break;
        case GT: coper = ">"; break;
        case LE: coper = "<="; break;
        case GE: coper = ">="; break;
        case NE: coper = "!="; break;
        case EQ: coper = "=="; break;
        case LEG: coper = "<=>"; break;
        case AND: coper = "&&"; break;
        case OR: coper = "||"; break;
733
734
735
736
737
        case '^': coper = "^"; break;
        case '+': coper = "+"; break;
        case '-': coper = "-"; break;
        case '*': coper = "*"; break;
        case '/': coper = "/"; break;
738
        default: cdoAbort("Internal error, expr operator %d not implemented!\n", oper);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
740
741
        }
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
742
  nodeType *p = nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
743

744
  if (p1->type == typeVar && p2->type == typeVar)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
746
      p = expr_var_var(init, oper, p1, p2);
747
      if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
748
749
        cdoPrint("\t%s\tarith\t%s[L%zu][N%zu] = %s %s %s", ExIn[init], p->u.var.nm, p->param.nlev, p->param.ngp, p1->u.var.nm,
                 coper, p2->u.var.nm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
750
    }
751
  else if (p1->type == typeVar && p2->type == typeCon)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
753
      p = expr_var_con(init, oper, p1, p2);
754
      if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
756
        cdoPrint("\t%s\tarith\t%s[L%zu][N%zu] = %s %s %g", ExIn[init], p->u.var.nm, p->param.nlev, p->param.ngp, p1->u.var.nm,
                 coper, p2->u.con.value);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757
    }
758
  else if (p1->type == typeCon && p2->type == typeVar)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
760
      p = expr_con_var(init, oper, p1, p2);
761
      if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
763
        cdoPrint("\t%s\tarith\t%s[L%zu][N%zu] = %g %s %s", ExIn[init], p->u.var.nm, p->param.nlev, p->param.ngp, p1->u.con.value,
                 coper, p2->u.var.nm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
    }
765
  else if (p1->type == typeCon && p2->