expr.cc

/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  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.
*/

#include <cstdlib>
#include <cassert>
#include <cdi.h>

#include "cdo_options.h"
#include "dmemory.h"
#include "process_int.h"
#include "field.h"
#include "expr.h"
#include "expr_fun.h"
#include "expr_yacc.h"
#include "cdo_zaxis.h"

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;

enum
{
  FT_STD,
  FT_CONST,
  FT_FLD,
  FT_ZON,
  FT_VERT,
  FT_COORD,
  FT_1C,
  FT_2C,
  FT_0
};

// clang-format off
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); }
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); }

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); }
static double f_rand(const double x) { (void)x; return ((double) std::rand()) / ((double) RAND_MAX); }
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; }
static double f_isMissval(const double x) { (void)x; return 0; }
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]; }
// clang-format on

struct func_t
{
  int type;
  int flag;
  const char *name;    // function name
  void (*func)(void);  // pointer to function
};

static const func_t fun_sym_tbl[] = {
  // scalar functions
  { 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) },
  { FT_STD, 0, "rand", reinterpret_cast<void (*)(void)>(&f_rand) },
  { 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) },
  { FT_STD, 0, "isMissval", reinterpret_cast<void (*)(void)>(&f_isMissval) },

  // constant functions
  { 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

  // CDO field functions (Reduce grid to point)
  { FT_FLD, 0, "fldmin", reinterpret_cast<void (*)(void)>(&fieldMin) },
  { FT_FLD, 0, "fldmax", reinterpret_cast<void (*)(void)>(&fieldMax) },
  { FT_FLD, 0, "fldsum", reinterpret_cast<void (*)(void)>(&fieldSum) },
  { FT_FLD, 1, "fldmean", reinterpret_cast<void (*)(void)>(&fieldMeanw) },
  { FT_FLD, 1, "fldavg", reinterpret_cast<void (*)(void)>(&fieldAvgw) },
  { 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) },

  // CDO zonal functions (Reduce grid to point)
  { 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) },

  // CDO field functions (Reduce level to point)
  { FT_VERT, 0, "vertmin", reinterpret_cast<void (*)(void)>(&fieldMin) },
  { FT_VERT, 0, "vertmax", reinterpret_cast<void (*)(void)>(&fieldMax) },
  { FT_VERT, 0, "vertsum", reinterpret_cast<void (*)(void)>(&fieldSum) },
  { FT_VERT, 1, "vertmean", reinterpret_cast<void (*)(void)>(&fieldMeanw) },
  { FT_VERT, 1, "vertavg", reinterpret_cast<void (*)(void)>(&fieldAvgw) },
  { 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) },

  { 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 },

  { 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) },

  { 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},
};

static const int NumFunc = sizeof(fun_sym_tbl) / sizeof(fun_sym_tbl[0]);

static void
node_data_delete(nodeType *p)
{
  if (p)
    {
      if (p->param.data)
        {
          Free(p->param.data);
          p->param.data = nullptr;
        }
    }
}

static void
node_delete(nodeType *p)
{
  if (p)
    {
      if (p->type == typeVar) node_data_delete(p);
      Free(p);
    }
}

static int
get_funcID(const char *fun)
{
  int funcID = -1;
  for (int i = 0; i < NumFunc; i++)
    if (cdo_cmpstr(fun, fun_sym_tbl[i].name))
      {
        funcID = i;
        break;
      }

  if (funcID == -1) cdoAbort("Function >%s< not available!", fun);

  return funcID;
}

static constexpr bool
isCompare(const int oper) noexcept
{
  return oper == LEG || oper == GE || oper == LE || oper == EQ || oper == NE || oper == GT || oper == LT;
}

static void
param_meta_copy(paramType &out, const paramType &in)
{
  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;
  out.name = nullptr;
  out.longname = nullptr;
  out.units = nullptr;
  out.data = nullptr;
}

static nodeType *
expr_con_con(const int oper, const nodeType *p1, const nodeType *p2)
{
  nodeType *p = (nodeType *) Calloc(1, sizeof(nodeType));

  p->type = typeCon;
  p->ltmpobj = true;

  double cval1 = p1->u.con.value;
  const double cval2 = p2->u.con.value;

  // clang-format off
  switch (oper)
    {
    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;
    case '+':  cval1 = cval1 + cval2; break;
    case '-':  cval1 = cval1 - cval2; break;
    case '*':  cval1 = cval1 * cval2; break;
    case '/':  cval1 = cval1 / cval2; break;
    case '^':  cval1 = std::pow(cval1, cval2); break;
    default:   cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
    }
  // clang-format on

  p->u.con.value = cval1;

  return p;
}

static void
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)
{
  size_t i;

  // clang-format off
  if (nmiss)
    {
      switch (oper)
        {
        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;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  else
    {
      switch (oper)
        {
        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;
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(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] = cval * idat[i]; break;
        case '/': for (i=0; i<n; ++i) odat[i] = DIVMN(cval, idat[i]); break;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  // clang-format on
}

static void
oper_expr_var_con(const int oper, const bool lmiss, const size_t n, const double missval1, const double missval2,
                  double *restrict odat, const double *restrict idat, const double cval)
{
  size_t i;

  // clang-format off
  if (lmiss)
    {
      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;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  else
    {
      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;
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(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 '*': 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;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  // clang-format on
}

static void
oper_expr_var_var(const int oper, const bool lmiss, const size_t n, const double missval1, const double missval2,
                  double *restrict odat, const double *restrict idat1, const double *restrict idat2)
{
  size_t i;

  // clang-format off
  if (lmiss)
    {
      switch (oper)
        {
        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;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  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;
        case '^': for (i=0; i<n; ++i) odat[i] = std::pow(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;
        case '*': for (i=0; i<n; ++i) odat[i] = idat1[i] * idat2[i]; break;
        case '/': for (i=0; i<n; ++i) odat[i] = IS_EQUAL(idat2[i], 0.) ? missval1 : idat1[i] / idat2[i]; break;
        default: cdoAbort("%s: operator %d unsupported!", __func__, oper); break;
        }
    }
  // clang-format on
}

static nodeType *
expr_con_var(const int init, const int oper, const nodeType *p1, const nodeType *p2)
{
  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;

  const size_t n = ngp * nlev;

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

  p->type = typeVar;
  p->ltmpobj = true;
  p->u.var.nm = strdup(tmpvnm);
  param_meta_copy(p->param, p2->param);
  p->param.name = p->u.var.nm;

  if (!init)
    {
      p->param.data = (double *) Malloc(n * sizeof(double));
      double *restrict odat = p->param.data;
      const double *restrict idat = p2->param.data;
      double cval = p1->u.con.value;
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;

      oper_expr_con_var(oper, nmiss, n, missval1, missval2, odat, cval, idat);

      p->param.nmiss = arrayNumMV(n, odat, missval1);
    }

  return p;
}

static nodeType *
expr_var_con(const int init, const int oper, const nodeType *p1, const nodeType *p2)
{
  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;

  const size_t n = ngp * nlev;

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

  p->type = typeVar;
  p->ltmpobj = true;
  p->u.var.nm = strdup(tmpvnm);
  param_meta_copy(p->param, p1->param);
  p->param.name = p->u.var.nm;

  if (!init)
    {
      p->param.data = (double *) Malloc(n * sizeof(double));
      double *restrict odat = p->param.data;
      const double *restrict idat = p1->param.data;
      double cval = p2->u.con.value;
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;

      oper_expr_var_con(oper, nmiss, n, missval1, missval2, odat, idat, cval);

      p->param.nmiss = arrayNumMV(n, odat, missval1);
    }

  return p;
}

static nodeType *
expr_var_var(int init, int oper, nodeType *p1, nodeType *p2)
{
  nodeType *px = p1;
  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;

  const size_t ngp1 = p1->param.ngp;
  const size_t ngp2 = p2->param.ngp;
  size_t ngp = ngp1;

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

  if (ngp1 != ngp2)
    {
      if (ngp1 == 1 || ngp2 == 1)
        {
          if (ngp1 == 1)
            {
              ngp = ngp2;
              px = p2;
            }
        }
      else if (nlat1 == nlat2 && ngp1 > ngp2)
        {
          lzonal = true;
        }
      else
        {
          cdoAbort("%s: Number of grid points differ (%s[%zu] <-> %s[%zu])", __func__, p1->param.name, ngp1, p2->param.name, ngp2);
        }
    }

  const size_t nlev1 = p1->param.nlev;
  const size_t nlev2 = p2->param.nlev;

  size_t nlev = nlev1;
  if (nlev1 != nlev2)
    {
      if (nlev1 == 1 || nlev2 == 1)
        {
          if (nlev1 == 1)
            {
              nlev = nlev2;
              px = p2;
            }
        }
      else
        {
          cdoAbort("%s: Number of levels differ (%s[%zu] <-> %s[%zu])", __func__, p1->param.name, nlev1, p2->param.name, nlev2);
        }
    }

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

  p->type = typeVar;
  p->ltmpobj = true;
  p->u.var.nm = strdup(tmpvnm);

  param_meta_copy(p->param, px->param);

  if (p->param.steptype == TIME_CONSTANT)
    {
      const int steptype1 = p1->param.steptype;
      const int steptype2 = p2->param.steptype;
      if (steptype1 != TIME_CONSTANT)
        p->param.steptype = steptype1;
      else if (steptype2 != TIME_CONSTANT)
        p->param.steptype = steptype2;
    }

  p->param.name = p->u.var.nm;
  // printf("%s %s nmiss %zu %zu\n", p->u.var.nm, px->param.name, nmiss1, nmiss2);

  if (!init)
    {
      p->param.data = (double *) Malloc(ngp * nlev * sizeof(double));

      for (size_t k = 0; k < nlev; k++)
        {
          size_t loff = k * ngp;
          const size_t loff1 = (nlev1 > 1) ? k * ngp1 : 0;
          const size_t loff2 = (nlev2 > 1) ? k * ngp2 : 0;

          const double *restrict idat1 = p1->param.data + loff1;
          const double *restrict idat2 = p2->param.data + loff2;
          double *restrict odat = p->param.data + loff;
          const size_t nmiss = nmiss1 || nmiss2;

          if (ngp1 != ngp2)
            {
              if (lzonal)
                {
                  const size_t nlon = ngp1 / nlat1;
                  for (size_t j = 0; j < nlat1; ++j)
                    oper_expr_var_con(oper, nmiss, nlon, missval1, missval2, odat + j * nlon, idat1 + j * nlon, idat2[j]);
                }
              else
                {
                  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);
                }
            }
          else
            {
              oper_expr_var_var(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2);
            }
        }

      p->param.nmiss = arrayNumMV(ngp * nlev, p->param.data, missval1);
    }

  return p;
}

static void
ex_copy_var(int init, nodeType *p2, const nodeType *p1)
{
  if (Options::cdoVerbose)
    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);

  const size_t ngp = p1->param.ngp;
  assert(ngp > 0);

  if (ngp != p2->param.ngp)
    cdoAbort("%s: Number of grid points differ (%s[%zu] = %s[%zu])", __func__, p2->param.name, p2->param.ngp, p1->param.name, ngp);

  const size_t nlev = p1->param.nlev;
  assert(nlev > 0);

  if (nlev != p2->param.nlev)
    cdoAbort("%s: Number of levels differ (%s[%zu] = %s[%zu])", __func__, p2->param.name, p2->param.nlev, p1->param.name, nlev);

  if (!init)
    {
      arrayCopy(ngp * nlev, p1->param.data, p2->param.data);
      p2->param.missval = p1->param.missval;
      p2->param.nmiss = p1->param.nmiss;
    }
}

static void
ex_copy_con(int init, nodeType *p2, const nodeType *p1)
{
  const double cval = p1->u.con.value;

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

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

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

  if (!init)
    {
      assert(p2->param.data != nullptr);
      varrayFill(ngp * nlev, p2->param.data, cval);
    }
}

static void
ex_copy(int init, nodeType *p2, const nodeType *p1)
{
  if (p1->type == typeCon)
    ex_copy_con(init, p2, p1);
  else
    ex_copy_var(init, p2, p1);
}

static nodeType *
expr(const int init, const int oper, nodeType *p1, nodeType *p2)
{
  if (p1 == nullptr || p2 == nullptr) return nullptr;

  const char *coper = "???";

  if (Options::cdoVerbose)
    {
      switch (oper)
        {
        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;
        case '^': coper = "^"; break;
        case '+': coper = "+"; break;
        case '-': coper = "-"; break;
        case '*': coper = "*"; break;
        case '/': coper = "/"; break;
        default: cdoAbort("Internal error, expr operator %d not implemented!\n", oper);
        }
    }

  nodeType *p = nullptr;

  if (p1->type == typeVar && p2->type == typeVar)
    {
      p = expr_var_var(init, oper, p1, p2);
      if (Options::cdoVerbose)
        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);
    }
  else if (p1->type == typeVar && p2->type == typeCon)
    {
      p = expr_var_con(init, oper, p1, p2);
      if (Options::cdoVerbose)
        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);
    }
  else if (p1->type == typeCon && p2->type == typeVar)
    {
      p = expr_con_var(init, oper, p1, p2);
      if (Options::cdoVerbose)
        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);
    }
  else if (p1->type == typeCon && p2->