cdo-2.0.0/src/expr.cc

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

  Author: Uwe Schulzweida

*/

#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
auto expr_func_con_var = [](auto hasMissvals, auto n, auto mv, const auto cVal, const auto &vIn, auto &vOut, auto binray_operator)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = dbl_is_equal(vIn[i], mv) ? mv : binray_operator(cVal, vIn[i]);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(vIn[i], mv) ? mv : binray_operator(cVal, vIn[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binray_operator(cVal, vIn[i]);
    }
};

auto expr_mul_con_var = [](auto hasMissvals, auto n, auto mv, const auto cVal, const auto &vIn, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(cVal, 0.0) ? 0.0 : dbl_is_equal(vIn[i], mv) ? mv : binary_op_MUL(cVal, vIn[i]);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(cVal, 0.0) ? 0.0 : is_equal(vIn[i], mv) ? mv : binary_op_MUL(cVal, vIn[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binary_op_MUL(cVal, vIn[i]);
    }
};

auto expr_div_con_var = [](auto hasMissvals, auto n, auto mv, const auto cVal, const auto &vIn, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = (dbl_is_equal(vIn[i], mv) || dbl_is_equal(vIn[i], 0.0)) ? mv : binary_op_DIV(cVal, vIn[i]);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = (is_equal(vIn[i], mv) || is_equal(vIn[i], 0.0)) ? mv : binary_op_DIV(cVal, vIn[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(vIn[i], 0.0) ? mv : binary_op_DIV(cVal, vIn[i]);
    }
};

auto expr_func_var_con = [](auto hasMissvals, auto n, auto mv, const auto &vIn, const auto cVal, auto &vOut, auto binray_operator)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = dbl_is_equal(vIn[i], mv) ? mv : binray_operator(vIn[i], cVal);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(vIn[i], mv) ? mv : binray_operator(vIn[i], cVal);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binray_operator(vIn[i], cVal);
    }
};

auto expr_mul_var_con = [](auto hasMissvals, auto n, auto mv, const auto &vIn, const auto cVal, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(cVal, 0.0) ? 0.0 : dbl_is_equal(vIn[i], mv) ? mv : binary_op_MUL(vIn[i], cVal);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(cVal, 0.0) ? 0.0 : is_equal(vIn[i], mv) ? mv : binary_op_MUL(vIn[i], cVal);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binary_op_MUL(vIn[i], cVal);
    }
};

auto expr_div_var_con = [](auto hasMissvals, auto n, auto mv, const auto &vIn, const auto cVal, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv))
        for (size_t i = 0; i < n; ++i) vOut[i] = (dbl_is_equal(vIn[i], mv) || dbl_is_equal(cVal, 0.0)) ? mv : binary_op_DIV(vIn[i], cVal);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = (is_equal(vIn[i], mv) || is_equal(cVal, 0.0)) ? mv : binary_op_DIV(vIn[i], cVal);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(cVal, 0.0) ? mv : binary_op_DIV(vIn[i], cVal);
    }
};

auto expr_func_var_var = [](auto hasMissvals, auto n, auto mv1, auto mv2, const auto &vIn1, const auto &vIn2, auto &vOut, auto binray_operator)
{
  if (hasMissvals)
    {
      if (std::isnan(mv1) || std::isnan(mv2))
        for (size_t i = 0; i < n; ++i) vOut[i] = (dbl_is_equal(vIn1[i], mv1) || dbl_is_equal(vIn2[i], mv2)) ? mv1 : binray_operator(vIn1[i], vIn2[i]);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = (is_equal(vIn1[i], mv1) || is_equal(vIn2[i], mv2)) ? mv1 : binray_operator(vIn1[i], vIn2[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binray_operator(vIn1[i], vIn2[i]);
    }
};

auto expr_mul_var_var = [](auto hasMissvals, auto n, auto mv1, auto mv2, const auto &vIn1, const auto &vIn2, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv1) || std::isnan(mv2))
        for (size_t i = 0; i < n; ++i)
          vOut[i] = (dbl_is_equal(vIn1[i], 0.0) || dbl_is_equal(vIn2[i], 0.0)) ? 0.0 : (dbl_is_equal(vIn1[i], mv1) || dbl_is_equal(vIn2[i], mv2)) ? mv1 : binary_op_MUL(vIn1[i], vIn2[i]);
      else
        for (size_t i = 0; i < n; ++i)
          vOut[i] = (is_equal(vIn1[i], 0.0) || is_equal(vIn2[i], 0.0)) ? 0.0 : (is_equal(vIn1[i], mv1) || is_equal(vIn2[i], mv2)) ? mv1 : binary_op_MUL(vIn1[i], vIn2[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = binary_op_MUL(vIn1[i], vIn2[i]);
    }
};

auto expr_div_var_var = [](auto hasMissvals, auto n, auto mv1, auto mv2, const auto &vIn1, const auto &vIn2, auto &vOut)
{
  if (hasMissvals)
    {
      if (std::isnan(mv1) || std::isnan(mv2))
        for (size_t i = 0; i < n; ++i) vOut[i] = (dbl_is_equal(vIn1[i], mv1) || dbl_is_equal(vIn2[i], mv2) || dbl_is_equal(vIn2[i], 0.0)) ? mv1 : binary_op_DIV(vIn1[i], vIn2[i]);
      else
        for (size_t i = 0; i < n; ++i) vOut[i] = (is_equal(vIn1[i], mv1) || is_equal(vIn2[i], mv2) || is_equal(vIn2[i], 0.0)) ? mv1 : binary_op_DIV(vIn1[i], vIn2[i]);
    }
  else
    {
      for (size_t i = 0; i < n; ++i) vOut[i] = is_equal(vIn2[i], 0.0) ? mv1 : binary_op_DIV(vIn1[i], vIn2[i]);
    }
};
// clang-format on

// clang-format off
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.0; }
static double f_deg(const double x) { return x * 180.0 / M_PI; }
static double f_isMissval(const double x) { (void)x; return 0.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
};

// clang-format off
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)>(&field_min) },
  { FT_FLD, 0, "fldmax",  reinterpret_cast<void (*)(void)>(&field_max) },
  { FT_FLD, 0, "fldsum",  reinterpret_cast<void (*)(void)>(&field_sum) },
  { FT_FLD, 1, "fldmean", reinterpret_cast<void (*)(void)>(&field_meanw) },
  { FT_FLD, 1, "fldavg",  reinterpret_cast<void (*)(void)>(&field_avgw) },
  { FT_FLD, 1, "fldstd",  reinterpret_cast<void (*)(void)>(&field_stdw) },
  { FT_FLD, 1, "fldstd1", reinterpret_cast<void (*)(void)>(&field_std1w) },
  { FT_FLD, 1, "fldvar",  reinterpret_cast<void (*)(void)>(&field_varw) },
  { FT_FLD, 1, "fldvar1", reinterpret_cast<void (*)(void)>(&field_var1w) },

  // CDO zonal functions (Reduce grid to point)
  { FT_ZON, 0, "zonmin",  reinterpret_cast<void (*)(void)>(&zonal_min) },
  { FT_ZON, 0, "zonmax",  reinterpret_cast<void (*)(void)>(&zonal_max) },
  { FT_ZON, 0, "zonsum",  reinterpret_cast<void (*)(void)>(&zonal_sum) },
  { FT_ZON, 0, "zonmean", reinterpret_cast<void (*)(void)>(&zonal_mean) },
  { FT_ZON, 0, "zonavg",  reinterpret_cast<void (*)(void)>(&zonal_avg) },
  { FT_ZON, 0, "zonstd",  reinterpret_cast<void (*)(void)>(&zonal_std) },
  { FT_ZON, 0, "zonstd1", reinterpret_cast<void (*)(void)>(&zonal_std1) },
  { FT_ZON, 0, "zonvar",  reinterpret_cast<void (*)(void)>(&zonal_var) },
  { FT_ZON, 0, "zonvar1", reinterpret_cast<void (*)(void)>(&zonal_var1) },

  // CDO field functions (Reduce level to point)
  { FT_VERT, 0, "vertmin",  reinterpret_cast<void (*)(void)>(&field_min) },
  { FT_VERT, 0, "vertmax",  reinterpret_cast<void (*)(void)>(&field_max) },
  { FT_VERT, 0, "vertsum",  reinterpret_cast<void (*)(void)>(&field_sum) },
  { FT_VERT, 1, "vertmean", reinterpret_cast<void (*)(void)>(&field_meanw) },
  { FT_VERT, 1, "vertavg",  reinterpret_cast<void (*)(void)>(&field_avgw) },
  { FT_VERT, 1, "vertstd",  reinterpret_cast<void (*)(void)>(&field_stdw) },
  { FT_VERT, 1, "vertstd1", reinterpret_cast<void (*)(void)>(&field_std1w) },
  { FT_VERT, 1, "vertvar",  reinterpret_cast<void (*)(void)>(&field_varw) },
  { FT_VERT, 1, "vertvar1", reinterpret_cast<void (*)(void)>(&field_var1w) },

  { FT_COORD, 0, "clon",       nullptr },
  { FT_COORD, 0, "clat",       nullptr },
  { FT_COORD, 0, "clev",       nullptr },
  { FT_COORD, 0, "cdeltaz",    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},
};
// clang-format on

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) cdo_abort("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.hasMissvals = 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)
{
  auto p = (nodeType *) Calloc(1, sizeof(nodeType));

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

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

  // clang-format off
  switch (oper)
    {
    case LT:   cval1 = binary_op_LT(cval1, cval2); break;
    case GT:   cval1 = binary_op_GT(cval1, cval2); break;
    case LE:   cval1 = binary_op_LE(cval1, cval2); break;
    case GE:   cval1 = binary_op_GE(cval1, cval2); break;
    case NE:   cval1 = binary_op_NE(cval1, cval2); break;
    case EQ:   cval1 = binary_op_EQ(cval1, cval2); break;
    case LEG:  cval1 = binary_op_LEG(cval1, cval2); break;
    case AND:  cval1 = binary_op_AND(cval1, cval2); break;
    case OR:   cval1 = binary_op_OR(cval1, cval2); break;
    case '^':  cval1 = binary_op_POW(cval1, cval2); break;
    case '+':  cval1 = binary_op_ADD(cval1, cval2); break;
    case '-':  cval1 = binary_op_SUB(cval1, cval2); break;
    case '*':  cval1 = binary_op_MUL(cval1, cval2); break;
    case '/':  cval1 = binary_op_DIV(cval1, cval2); break;
    default:   cdo_abort("%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 hasMissvals, const size_t n, const double mv, double *odat,
                  const double cval, const double *idat)
{
  // clang-format off
  switch (oper)
    {
    case LT:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_LT); break;
    case GT:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_GT); break;
    case LE:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_LE); break;
    case GE:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_GE); break;
    case NE:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_NE); break;
    case EQ:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_EQ); break;
    case LEG: expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_LEG); break;
    case AND: expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_AND); break;
    case OR:  expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_OR); break;
    case '^': expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_POW); break;
    case '+': expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_ADD); break;
    case '-': expr_func_con_var(hasMissvals, n, mv, cval, idat, odat, binary_op_SUB); break;
    case '*': expr_mul_con_var(hasMissvals, n, mv, cval, idat, odat); break;
    case '/': expr_div_con_var(hasMissvals, n, mv, cval, idat, odat); break;
    default: cdo_abort("%s: operator %d unsupported!", __func__, oper); break;
    }
  // clang-format on
}

static void
oper_expr_var_con(const int oper, const bool hasMissvals, const size_t n, const double mv,
                  double *odat, const double *idat, const double cval)
{
  // clang-format off
  switch (oper)
    {
    case LT:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_LT); break;
    case GT:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_GT); break;
    case LE:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_LE); break;
    case GE:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_GE); break;
    case NE:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_NE); break;
    case EQ:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_EQ); break;
    case LEG: expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_LEG); break;
    case AND: expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_AND); break;
    case OR:  expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_OR); break;
    case '^': expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_POW); break;
    case '+': expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_ADD); break;
    case '-': expr_func_var_con(hasMissvals, n, mv, idat, cval, odat, binary_op_SUB); break;
    case '*': expr_mul_var_con(hasMissvals, n, mv, idat, cval, odat); break;
    case '/': expr_div_var_con(hasMissvals, n, mv, idat, cval, odat); break;
    default: cdo_abort("%s: operator %d unsupported!", __func__, oper); break;
    }
  // clang-format on
}

static void
oper_expr_var_var(const int oper, const bool hasMissvals, const size_t n, const double mv1, const double mv2,
                  double *odat, const double *idat1, const double *idat2)
{
  // clang-format off
  switch (oper)
    {
    case LT:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_LT); break;
    case GT:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_GT); break;
    case LE:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_LE); break;
    case GE:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_GE); break;
    case NE:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_NE); break;
    case EQ:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_EQ); break;
    case LEG: expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_LEG); break;
    case AND: expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_AND); break;
    case OR:  expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_OR); break;
    case '^': expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_POW); break;
    case '+': expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_ADD); break;
    case '-': expr_func_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat, binary_op_SUB); break;
    case '*': expr_mul_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat); break;
    case '/': expr_div_var_var(hasMissvals, n, mv1, mv2, idat1, idat2, odat); break;
    default: cdo_abort("%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 auto ngp = (p2->param.ngp > 0) ? p2->param.ngp : 1;
  const auto nlev = (p2->param.nlev > 0) ? p2->param.nlev : 1;
  const auto nmiss = p2->param.nmiss;
  const auto datatype = p2->param.datatype;
  const auto missval2 = p2->param.missval;

  const auto n = ngp * nlev;

  auto 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));
      auto odat = p->param.data;
      const auto idat = p2->param.data;
      auto cval = p1->u.con.value;
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;

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

      p->param.nmiss = array_num_mv(n, odat, missval2);
    }

  return p;
}

static nodeType *
expr_var_con(const int init, const int oper, const nodeType *p1, const nodeType *p2)
{
  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  const auto nmiss = p1->param.nmiss;
  const auto datatype = p1->param.datatype;
  const auto missval1 = p1->param.missval;

  const auto n = ngp * nlev;

  auto 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));
      auto odat = p->param.data;
      const auto idat = p1->param.data;
      auto cval = p2->u.con.value;
      if (datatype == CDI_DATATYPE_FLT32 && isCompare(oper)) cval = (float) cval;

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

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

  return p;
}

static nodeType *
expr_var_var(int init, int oper, nodeType *p1, nodeType *p2)
{
  auto px = p1;
  const auto nmiss1 = p1->param.nmiss;
  const auto nmiss2 = p2->param.nmiss;
  const auto missval1 = p1->param.missval;
  const auto missval2 = p2->param.missval;

  const auto ngp1 = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto ngp2 = (p2->param.ngp > 0) ? p2->param.ngp : 1;
  auto ngp = ngp1;

  const auto nlat1 = p1->param.nlat;
  const auto nlat2 = p2->param.nlat;
  auto 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
        {
          cdo_abort("%s: Number of grid points differ (%s[%zu] <-> %s[%zu])", __func__, p1->param.name, ngp1, p2->param.name, ngp2);
        }
    }

  const auto nlev1 = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  const auto nlev2 = (p2->param.nlev > 0) ? p2->param.nlev : 1;

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

  auto 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 auto steptype1 = p1->param.steptype;
      const auto 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++)
        {
          const auto loff = k * ngp;
          const size_t loff1 = (nlev1 > 1) ? k * ngp1 : 0;
          const size_t loff2 = (nlev2 > 1) ? k * ngp2 : 0;

          const auto idat1 = p1->param.data + loff1;
          const auto idat2 = p2->param.data + loff2;
          auto odat = p->param.data + loff;
          const size_t nmiss = (nmiss1 || nmiss2);

          if (ngp1 != ngp2)
            {
              if (lzonal)
                {
                  const auto nlon = ngp1 / nlat1;
                  for (size_t j = 0; j < nlat1; ++j)
                    oper_expr_var_con(oper, nmiss, nlon, missval1, odat + j * nlon, idat1 + j * nlon, idat2[j]);
                }
              else
                {
                  if (ngp2 == 1)
                    oper_expr_var_con(oper, nmiss, ngp, missval1, odat, idat1, idat2[0]);
                  else
                    oper_expr_con_var(oper, nmiss, ngp, missval2, odat, idat1[0], idat2);
                }
            }
          else
            {
              oper_expr_var_var(oper, nmiss, ngp, missval1, missval2, odat, idat1, idat2);
            }
        }

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

  return p;
}

static void
ex_copy_var(int init, nodeType *p2, const nodeType *p1)
{
  const auto copyConstValue1 = (p1->param.ngp == 0 && p1->param.nlev == 0);
  const auto ngp1 = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev1 = (p1->param.nlev > 0) ? p1->param.nlev : 1;

  if (Options::cdoVerbose)
    {
      if (copyConstValue1)
        cdo_print("\t%s\tcopy\t%s[N%zu][L%zu] = %s", ExIn[init], p2->param.name, p2->param.ngp, p2->param.nlev, p1->param.name);
      else
        cdo_print("\t%s\tcopy\t%s[N%zu][L%zu] = %s[N%zu][L%zu]", ExIn[init], p2->param.name, p2->param.ngp, p2->param.nlev,
                  p1->param.name, p1->param.ngp, p1->param.nlev);
    }

  const auto ngp2 = (p2->param.ngp > 0) ? p2->param.ngp : 1;
  const auto nlev2 = (p2->param.nlev > 0) ? p2->param.nlev : 1;

  if (!copyConstValue1 && ngp2 != ngp1)
    cdo_abort("%s: Number of grid points differ (%s[N%zu] = %s[N%zu])", __func__, p2->param.name, ngp2, p1->param.name, ngp1);

  if (!copyConstValue1 && nlev2 != nlev1)
    cdo_abort("%s: Number of levels differ (%s[L%zu] = %s[L%zu])", __func__, p2->param.name, nlev2, p1->param.name, nlev1);

  if (!init)
    {
      if (copyConstValue1)
        {
          varray_fill(ngp2 * nlev2, p2->param.data, p1->param.data[0]);
        }
      else
        {
          array_copy(ngp2 * nlev2, 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 auto cval = p1->u.con.value;

  auto ngp = p2->param.ngp;
  auto nlev = p2->param.nlev;

  if (Options::cdoVerbose)
    {
      if (ngp == 0 && nlev == 0)
        cdo_print("\t%s\tcopy\t%s = %g", ExIn[init], p2->param.name, cval);
      else
        cdo_print("\t%s\tcopy\t%s[N%zu][L%zu] = %g", ExIn[init], p2->param.name, ngp, nlev, cval);
    }

  if (ngp == 0 && nlev == 0)
    {
      ngp = 1;
      nlev = 1;
    }

  assert(ngp > 0);
  assert(nlev > 0);

  if (!init)
    {
      assert(p2->param.data != nullptr);
      varray_fill(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: cdo_abort("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)
        cdo_print("\t%s\tarith\t%s[N%zu][L%zu] = %s %s %s", ExIn[init], p->u.var.nm, p->param.ngp, p->param.nlev, 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)
        cdo_print("\t%s\tarith\t%s[N%zu][L%zu] = %s %s %g", ExIn[init], p->u.var.nm, p->param.ngp, p->param.nlev, 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)
        cdo_print("\t%s\tarith\t%s[N%zu][L%zu] = %g %s %s", ExIn[init], p->u.var.nm, p->param.ngp, p->param.nlev, p1->u.con.value,
                  coper, p2->u.var.nm);
    }
  else if (p1->type == typeCon && p2->type == typeCon)
    {
      p = expr_con_con(oper, p1, p2);
      if (Options::cdoVerbose)
        cdo_print("\t%s\tarith\t%g = %g %s %g", ExIn[init], p->u.con.value, p1->u.con.value, coper, p2->u.con.value);
    }
  else
    cdo_abort("Internal problem!");

  if (p1->ltmpobj) node_delete(p1);
  if (p2->ltmpobj) node_delete(p2);

  return p;
}

static nodeType *
ex_fun0_con(const int init, const int funcID, double *data)
{
  auto functype = fun_sym_tbl[funcID].type;
  if (functype != FT_0) cdo_abort("Function %s not available for constant values!", fun_sym_tbl[funcID].name);

  if (Options::cdoVerbose) cdo_print("\t%s\tfunc\t%s", ExIn[init], fun_sym_tbl[funcID].name);

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

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

  if (!init)
    {
      double (*exprfunc)(const double *) = (double (*)(const double *)) fun_sym_tbl[funcID].func;
      p->u.con.value = exprfunc(data);
    }

  return p;
}

static nodeType *
ex_fun_con(const int funcID, nodeType *p1)
{
  const auto functype = fun_sym_tbl[funcID].type;
  if (functype != FT_STD) cdo_abort("Function %s not available for constant values!", fun_sym_tbl[funcID].name);

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

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

  double (*exprfunc)(double) = (double (*)(double)) fun_sym_tbl[funcID].func;
  p->u.con.value = exprfunc(p1->u.con.value);

  if (p1->ltmpobj)
    node_delete(p1);
  else
    Free(p1);

  return p;
}

static void
ex_fun_std(const int funcID, size_t ngp, size_t nlev, size_t nmiss, double missval, double *p1data, double *pdata)
{
  double (*exprfunc)(double) = (double (*)(double)) fun_sym_tbl[funcID].func;
  if (nmiss)
    {
      if (cdo_cmpstr(fun_sym_tbl[funcID].name, "isMissval"))
        {
          for (size_t i = 0; i < ngp * nlev; i++)
            {
              pdata[i] = DBL_IS_EQUAL(p1data[i], missval);
            }
        }
      else
        {
          for (size_t i = 0; i < ngp * nlev; i++)
            {
              errno = -1;
              pdata[i] = DBL_IS_EQUAL(p1data[i], missval) ? missval : exprfunc(p1data[i]);
              if (errno == EDOM || errno == ERANGE || std::isnan(pdata[i])) pdata[i] = missval;
            }
        }
    }
  else
    {
      for (size_t i = 0; i < ngp * nlev; i++)
        {
          errno = -1;
          pdata[i] = exprfunc(p1data[i]);
          if (errno == EDOM || errno == ERANGE || std::isnan(pdata[i])) pdata[i] = missval;
        }
    }
}

static nodeType *
ex_fun_var(const int init, const int funcID, nodeType *p1)
{
  const auto funcname = fun_sym_tbl[funcID].name;
  const auto functype = fun_sym_tbl[funcID].type;
  const auto funcflag = fun_sym_tbl[funcID].flag;

  const auto gridID = p1->param.gridID;
  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  const auto nlat = p1->param.nlat;
  const auto nmiss = p1->param.nmiss;
  const auto missval = p1->param.missval;

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

  if (functype == FT_CONST)
    {
      p->type = typeCon;
      double (*exprfunc)(const paramType *) = (double (*)(const paramType *)) fun_sym_tbl[funcID].func;
      p->u.con.value = exprfunc(&p1->param);
    }
  else if (functype == FT_FLD)
    {
      p->param.gridID = pointID;
      p->param.ngp = 1;
    }
  else if (functype == FT_ZON)
    {
      if (zonalID == -1) cdo_abort("Function %s() is only available for regular 2D grids!", fun_sym_tbl[funcID].name);
      p->param.gridID = zonalID;
      p->param.ngp = nlat;
    }
  else if (functype == FT_VERT)
    {
      p->param.zaxisID = surfaceID;
      p->param.nlev = 1;
    }

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

      if (functype == FT_STD)
        {
          ex_fun_std(funcID, ngp, nlev, nmiss, missval, p1data, pdata);
        }
      else if (functype == FT_FLD)
        {
          Field field;
          field.resize(ngp);
          if (funcflag == 1)
            {
              assert(p1->param.weight != nullptr);
              field.weightv.resize(ngp);
            }

          double (*exprfunc)(const Field &) = (double (*)(const Field &)) fun_sym_tbl[funcID].func;
          for (size_t k = 0; k < nlev; k++)
            {
              fld_field_init(field, nmiss, missval, ngp, p1data + k * ngp, p1->param.weight);
              pdata[k] = exprfunc(field);
            }
        }
      else if (functype == FT_ZON)
        {
          Field field1, field2;
          field1.resize(ngp);
          field2.resize(nlat);
          void (*exprfunc)(const Field &, Field &) = (void (*)(const Field &, Field &)) fun_sym_tbl[funcID].func;
          for (size_t k = 0; k < nlev; k++)
            {
              fld_field_init(field1, nmiss, missval, ngp, &p1data[k * ngp], nullptr);
              field1.grid = gridID;
              fld_field_init(field2, nmiss, missval, nlat, &pdata[k * nlat], nullptr);
              exprfunc(field1, field2);
              array_copy(nlat, field2.vec_d.data(), &pdata[k * nlat]);
            }
        }
      else if (functype == FT_VERT)
        {
          Field field;
          field.resize(nlev);
          if (funcflag == 1) vert_weights(p1->param.zaxisID, nlev, field.weightv);
          double (*exprfunc)(const Field &) = (double (*)(const Field &)) fun_sym_tbl[funcID].func;
          for (size_t i = 0; i < ngp; i++)
            {
              for (size_t k = 0; k < nlev; k++) field.vec_d[k] = p1data[k * ngp + i];
              fld_field_init(field, nmiss, missval, nlev, nullptr, nullptr);
              pdata[i] = exprfunc(field);
            }
        }
      else if (functype == FT_CONST)
        {
        }
      else
        cdo_abort("Intermal error, wrong function type (%d) for %s()!", functype, funcname);

      p->param.nmiss = array_num_mv(p->param.ngp * p->param.nlev, pdata, missval);
    }

  if (p1->ltmpobj)
    node_delete(p1);
  else
    Free(p1);

  return p;
}

static nodeType *
ex_fun(const int init, const int funcID, nodeType *p1)
{
  nodeType *p = nullptr;

  if (p1->type == typeVar)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tfunc\t%s (%s)", ExIn[init], fun_sym_tbl[funcID].name, p1->u.var.nm);
      p = ex_fun_var(init, funcID, p1);
    }
  else if (p1->type == typeCon)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tfunc\t%s (%g)", ExIn[init], fun_sym_tbl[funcID].name, p1->u.con.value);
      p = ex_fun_con(funcID, p1);
    }
  else
    cdo_abort("Internal problem!");

  return p;
}

static size_t
get_levidx(const size_t nlev, const double *data, const double value, const char *funcname)
{
  size_t levidx;

  for (levidx = 0; levidx < nlev; ++levidx)
    if (IS_EQUAL(data[levidx], value)) break;
  if (levidx == nlev) cdo_abort("%s(): level %g not found!", funcname, value);

  return levidx;
}

static void
get_levidxrange(size_t nlev, const double *data, double value1, double value2, const char *funcname, size_t &levidx1,
                size_t &levidx2)
{
  long i, n = nlev;
  if (data[0] <= data[nlev - 1])
    {
      for (i = 0; i < n; ++i)
        if (data[i] >= value1) break;
      if (i == n) cdo_abort("%s(): lower level %g not found!", funcname, value1);
      levidx1 = i;

      for (i = n - 1; i >= 0; --i)
        if (data[i] <= value2) break;
      if (i < 0) cdo_abort("%s(): upper level %g not found!", funcname, value2);
      if (i < (long) levidx1) cdo_abort("%s(): level range %g to %g not found!", funcname, value1, value2);
      levidx2 = i;
    }
  else
    {
      for (i = 0; i < n; ++i)
        if (data[i] <= value2) break;
      if (i == n) cdo_abort("%s(): upper level %g not found!", funcname, value1);
      levidx1 = i;

      for (i = n - 1; i >= 0; --i)
        if (data[i] >= value1) break;
      if (i < 0) cdo_abort("%s(): lower level %g not found!", funcname, value2);
      if (i < (long) levidx1) cdo_abort("%s(): level range %g to %g not found!", funcname, value1, value2);
      levidx2 = i;
    }
}

static nodeType *
fun1c(int init, int funcID, nodeType *p1, double value, parseParamType *parse_arg)
{
  const auto funcname = fun_sym_tbl[funcID].name;
  if (p1->type != typeVar) cdo_abort("Parameter of function %s() needs to be a variable!", funcname);
  if (p1->ltmpobj) cdo_abort("Temporary objects not allowed in function %s()!", funcname);

  if (parse_arg->debug)
    cdo_print("\t%s\tfunc\t%s=%s(%s[N%zu][L%zu], %g)", ExIn[init], tmpvnm, funcname, p1->param.name, p1->param.ngp, p1->param.nlev,
              value);

  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  auto nmiss = p1->param.nmiss;
  const auto missval = p1->param.missval;

  auto 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)
    {
      if (p1->param.longname) p->param.longname = strdup(p1->param.longname);
      if (p1->param.units) p->param.units = strdup(p1->param.units);
    }

  p->param.nlev = 1;

  const auto zaxisID = p1->param.zaxisID;
  const auto coordID = params_get_coord_ID(parse_arg, 'z', zaxisID);

  std::vector<double> data;
  double *pdata = nullptr;

  if (init)
    {
      parse_arg->coords[coordID].needed = true;

      data.resize(nlev);
      pdata = data.data();
      cdo_zaxis_inq_levels(zaxisID, pdata);
    }
  else
    {
      pdata = parse_arg->coords[coordID].data.data();
    }

  size_t levidx = 0;
  if (cdo_cmpstr(funcname, "sellevidx"))
    {
      const auto ilevidx = std::lround(value);
      if (ilevidx < 1 || ilevidx > (long) nlev)
        cdo_abort("%s(): level index %ld out of range (range: 1-%zu)!", funcname, ilevidx, nlev);
      levidx = (size_t) ilevidx - 1;
    }
  else if (cdo_cmpstr(funcname, "sellevel"))
    {
      levidx = get_levidx(nlev, pdata, value, funcname);
    }
  else
    cdo_abort("Function %s() not implemented!", funcname);

  if (init)
    {
      const auto level = data[levidx];
      const auto zaxisID2 = zaxisCreate(zaxisInqType(zaxisID), 1);
      zaxisDefLevels(zaxisID2, &level);
      p->param.zaxisID = zaxisID2;
    }

  if (!init)
    {
      p->param.data = (double *) Malloc(ngp * sizeof(double));
      pdata = p->param.data;
      const auto p1data = p1->param.data + ngp * levidx;
      array_copy(ngp, p1data, pdata);
      if (nmiss) nmiss = array_num_mv(ngp, pdata, missval);
      p->param.nmiss = nmiss;
    }

  if (p1->ltmpobj) node_delete(p1);

  return p;
}

static nodeType *
fun2c(int init, int funcID, nodeType *p1, double value1, double value2, parseParamType *parse_arg)
{
  const auto funcname = fun_sym_tbl[funcID].name;
  if (p1->type != typeVar) cdo_abort("Parameter of function %s() needs to be a variable!", funcname);
  if (p1->ltmpobj) cdo_abort("Temporary objects not allowed in function %s()!", funcname);

  if (parse_arg->debug)
    cdo_print("\t%s\tfunc\t%s=%s(%s[N%zu][L%zu], %g, %g)", ExIn[init], tmpvnm, funcname, p1->param.name, p1->param.ngp,
              p1->param.nlev, value1, value2);

  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  auto nmiss = p1->param.nmiss;
  const auto missval = p1->param.missval;

  auto 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)
    {
      if (p1->param.longname) p->param.longname = strdup(p1->param.longname);
      if (p1->param.units) p->param.units = strdup(p1->param.units);
    }

  const auto zaxisID = p1->param.zaxisID;
  const auto coordID = params_get_coord_ID(parse_arg, 'z', zaxisID);
  size_t levidx1 = 0;
  size_t levidx2 = 0;

  double *data = nullptr;
  if (init)
    {
      parse_arg->coords[coordID].needed = true;

      data = (double *) Malloc(nlev * sizeof(double));
      cdo_zaxis_inq_levels(zaxisID, data);
    }
  else
    {
      data = parse_arg->coords[coordID].data.data();
    }

  if (cdo_cmpstr(funcname, "sellevidxrange"))
    {
      if (value2 < value1) cdo_abort("%s(): first level index is greater than last level index!", funcname);
      const auto ilevidx1 = std::lround(value1);
      if (ilevidx1 < 1 || ilevidx1 > (long) nlev)
        cdo_abort("%s(): level index1 %ld out of range (range: 1-%zu)!", funcname, ilevidx1, nlev);
      levidx1 = (size_t) ilevidx1 - 1;
      const auto ilevidx2 = std::lround(value2);
      if (ilevidx2 < 1 || ilevidx2 > (long) nlev)
        cdo_abort("%s(): level index2 %ld out of range (range: 1-%zu)!", funcname, ilevidx2, nlev);
      levidx2 = (size_t) ilevidx2 - 1;
    }
  else if (cdo_cmpstr(funcname, "sellevelrange"))
    {
      if (value2 < value1) cdo_abort("%s(): first level is greater than last level!", funcname);
      get_levidxrange(nlev, data, value1, value2, funcname, levidx1, levidx2);
    }
  else
    cdo_abort("Function %s() not implemented!", funcname);

  const int nlevout = levidx2 - levidx1 + 1;
  p->param.nlev = nlevout;

  if (init)
    {
      const auto zaxisID2 = zaxisCreate(zaxisInqType(zaxisID), nlevout);
      zaxisDefLevels(zaxisID2, &data[levidx1]);
      if (zaxisInqLbounds(zaxisID, nullptr) && zaxisInqUbounds(zaxisID, nullptr))
        {
          std::vector<double> bounds(nlev);
          zaxisInqLbounds(zaxisID, bounds.data());
          zaxisDefLbounds(zaxisID2, &bounds[levidx1]);
          zaxisInqUbounds(zaxisID, bounds.data());
          zaxisDefUbounds(zaxisID2, &bounds[levidx1]);
        }
      p->param.zaxisID = zaxisID2;
    }

  if (init) Free(data);

  if (!init)
    {
      p->param.data = (double *) Malloc(ngp * nlevout * sizeof(double));
      auto pdata = p->param.data;
      const auto p1data = p1->param.data + ngp * levidx1;
      array_copy(ngp * nlevout, p1data, pdata);
      if (nmiss) nmiss = array_num_mv(ngp * nlevout, pdata, missval);
      p->param.nmiss = nmiss;
    }

  if (p1->ltmpobj) node_delete(p1);

  return p;
}

static nodeType *
coord_fun(const int init, const int funcID, nodeType *p1, parseParamType *parse_arg)
{
  const auto funcname = fun_sym_tbl[funcID].name;
  if (p1->type != typeVar) cdo_abort("Parameter of function %s() needs to be a variable!", funcname);
  if (p1->ltmpobj) cdo_abort("Temporary objects not allowed in function %s()!", funcname);

  const auto len = 3 + strlen(p1->u.var.nm);
  auto cname = (char *) Calloc(len, 1);
  strcpy(cname, p1->u.var.nm);

  // clang-format off
  if      (cdo_cmpstr(funcname, "clon"))       strcat(cname, ".x");
  else if (cdo_cmpstr(funcname, "clat"))       strcat(cname, ".y");
  else if (cdo_cmpstr(funcname, "clev"))       strcat(cname, ".z");
  else if (cdo_cmpstr(funcname, "cdeltaz"))    strcat(cname, ".d");
  else if (cdo_cmpstr(funcname, "gridarea"))   strcat(cname, ".a");
  else if (cdo_cmpstr(funcname, "gridweight")) strcat(cname, ".w");
  else cdo_abort("Implementation missing for function %s!", funcname);
  // clang-format on

  Free(p1->u.var.nm);
  p1->u.var.nm = cname;

  auto p = expr_run(p1, parse_arg);
  p->param.hasMissvals = false;

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

      for (size_t i = 0; i < ngp*nlev; i++) pdata[i] = p1data[i];
      */
    }
  /*
  Free(cname);
  Free(p1);
  */
  return p;
}

static nodeType *
ex_uminus_var(const int init, nodeType *p1)
{
  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  const auto nmiss = p1->param.nmiss;
  const auto missval = p1->param.missval;

  auto 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(ngp * nlev * sizeof(double));
      double *pdata = p->param.data;
      const double *p1data = p1->param.data;

      if (nmiss)
        {
          for (size_t i = 0; i < ngp * nlev; ++i) pdata[i] = DBL_IS_EQUAL(p1data[i], missval) ? missval : -(p1data[i]);
        }
      else
        {
          for (size_t i = 0; i < ngp * nlev; ++i) pdata[i] = -(p1data[i]);
        }

      p->param.nmiss = nmiss;
    }

  if (p1->ltmpobj) node_delete(p1);

  return p;
}

static nodeType *
ex_uminus_con(nodeType *p1)
{
  auto p = (nodeType *) Calloc(1, sizeof(nodeType));

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

  p->u.con.value = -(p1->u.con.value);

  return p;
}

static nodeType *
ex_uminus(const int init, nodeType *p1)
{
  nodeType *p = nullptr;

  if (p1->type == typeVar)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tneg\t- (%s)", ExIn[init], p1->u.var.nm);
      p = ex_uminus_var(init, p1);
    }
  else if (p1->type == typeCon)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tneg\t- (%g)", ExIn[init], p1->u.con.value);
      p = ex_uminus_con(p1);
    }
  else
    cdo_abort("Internal problem!");

  return p;
}

static nodeType *
ex_not_var(const int init, nodeType *p1)
{
  const auto ngp = (p1->param.ngp > 0) ? p1->param.ngp : 1;
  const auto nlev = (p1->param.nlev > 0) ? p1->param.nlev : 1;
  const auto nmiss = p1->param.nmiss;
  const auto missval = p1->param.missval;

  auto 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(ngp * nlev * sizeof(double));
      double *pdata = p->param.data;
      const double *p1data = p1->param.data;

      if (nmiss)
        {
          for (size_t i = 0; i < ngp * nlev; ++i) pdata[i] = DBL_IS_EQUAL(p1data[i], missval) ? missval : unary_op_not(p1data[i]);
        }
      else
        {
          for (size_t i = 0; i < ngp * nlev; ++i) pdata[i] = unary_op_not(p1data[i]);
        }

      p->param.nmiss = nmiss;
    }

  if (p1->ltmpobj) node_delete(p1);

  return p;
}

static nodeType *
ex_not_con(const nodeType *p1)
{
  auto p = (nodeType *) Calloc(1, sizeof(nodeType));

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

  p->u.con.value = unary_op_not(p1->u.con.value);

  return p;
}

static nodeType *
ex_not(const int init, nodeType *p1)
{
  nodeType *p = nullptr;

  if (p1->type == typeVar)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tnot\t! (%s)", ExIn[init], p1->u.var.nm);
      p = ex_not_var(init, p1);
    }
  else if (p1->type == typeCon)
    {
      if (Options::cdoVerbose) cdo_print("\t%s\tnot\t! (%g)", ExIn[init], p1->u.con.value);
      p = ex_not_con(p1);
    }
  else
    cdo_abort("Internal problem!");

  return p;
}

static void
strAddNodeInfo(char *string, size_t stringlen, nodeType *p, const char *ext)
{
  const auto len = strlen(string);
  if (p->type == typeCon)
    snprintf(string + len, stringlen - len, "%g%s", p->u.con.value, ext);
  else
    snprintf(string + len, stringlen - len, "%s[N%zu][L%zu]%s", p->u.var.nm, p->param.ngp, p->param.nlev, ext);
}

static nodeType *
ex_ifelse(const int init, nodeType *p1, nodeType *p2, nodeType *p3)
{
  if (Options::cdoVerbose)
    {
      char strbuffer[1024];
      snprintf(strbuffer, sizeof(strbuffer), "\t%s\tifelse\t", ExIn[init]);
      strAddNodeInfo(strbuffer, sizeof(strbuffer), p1, " ? ");
      strAddNodeInfo(strbuffer, sizeof(strbuffer), p2, " : ");
      strAddNodeInfo(strbuffer, sizeof(strbuffer), p3, "");
      cdo_print(strbuffer);
    }

  nodeType *p0 = nullptr;
  nodeType *pnodes[3] = { p1, p2, p3 };
  for (unsigned i = 0; i < 3; ++i)
    if (pnodes[i]->type != typeCon)
      {
        p0 = pnodes[i];
        break;
      }

  if (p0 == nullptr) cdo_abort("expr?expr:expr: no data variable found!");

  auto missval = p0->param.missval;
  auto ngp = (p0->param.ngp > 0) ? p0->param.ngp : 1;
  auto nlev = (p0->param.nlev > 0) ? p0->param.nlev : 1;
  auto px = p0;

  size_t nmiss1 = 0;
  auto missval1 = missval;
  double *pdata1 = nullptr;
  size_t ngp1 = 1;
  size_t nlev1 = 1;

  if (p1->type == typeCon)
    {
      pdata1 = &p1->u.con.value;
    }
  else
    {
      nmiss1 = p1->param.nmiss;
      ngp1 = (p1->param.ngp > 0) ? p1->param.ngp : 1;
      nlev1 = (p1->param.nlev > 0) ? p1->param.nlev : 1;
      missval1 = p1->param.missval;
      pdata1 = p1->param.data;

      if (ngp1 > 1 && ngp1 != ngp)
        {
          if (ngp != 1) cdo_abort("expr?expr:expr: Number of grid points differ (ngp = %zu, ngp1 = %zu)", ngp, ngp1);

          ngp = ngp1;
          px = p1;
        }

      if (nlev1 > 1 && nlev1 != nlev)
        {
          if (nlev != 1) cdo_abort("expr?expr:expr: Number of levels differ (nlev = %zu, nlev1 = %zu)", nlev, nlev1);

          nlev = nlev1;
          px = p1;
        }
    }

  auto missval2 = missval1;
  double *pdata2 = nullptr;
  size_t ngp2 = 1;
  size_t nlev2 = 1;

  if (p2->type == typeCon)
    {
      pdata2 = &p2->u.con.value;
    }
  else
    {
      ngp2 = (p2->param.ngp > 0) ? p2->param.ngp : 1;
      nlev2 = (p2->param.nlev > 0) ? p2->param.nlev : 1;
      missval2 = p2->param.missval;
      pdata2 = p2->param.data;

      if (ngp2 > 1 && ngp2 != ngp)
        {
          if (ngp != 1) cdo_abort("expr?expr:expr: Number of grid points differ (ngp = %zu, ngp2 = %zu)", ngp, ngp2);

          ngp = ngp2;
          px = p2;
        }

      if (nlev2 > 1 && nlev2 != nlev)
        {
          if (nlev != 1) cdo_abort("expr?expr:expr: Number of levels differ (nlev = %zu, nlev2 = %zu)", nlev, nlev2);

          nlev = nlev2;
          px = p2;
        }
    }

  auto missval3 = missval1;
  double *pdata3 = nullptr;
  size_t ngp3 = 1;
  size_t nlev3 = 1;

  if (p3->type == typeCon)
    {
      pdata3 = &p3->u.con.value;
    }
  else
    {
      ngp3 = (p3->param.ngp > 0) ? p3->param.ngp : 1;
      nlev3 = (p3->param.nlev > 0) ? p3->param.nlev : 1;
      missval3 = p3->param.missval;
      pdata3 = p3->param.data;

      if (ngp3 > 1 && ngp3 != ngp)
        {
          if (ngp != 1) cdo_abort("expr?expr:expr: Number of grid points differ (ngp = %zu, ngp3 = %zu)", ngp, ngp3);

          ngp = ngp3;
          px = p3;
        }

      if (nlev3 > 1 && nlev3 != nlev)
        {
          if (nlev != 1) cdo_abort("expr?expr:expr: Number of levels differ (nlev = %zu, nlev3 = %zu)", nlev, nlev3);

          nlev = nlev3;
          px = p3;
        }
    }

  auto 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);
  p->param.name = p->u.var.nm;

  if (!init)
    {
      size_t nmiss = 0;

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

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

          const auto idat1 = pdata1 + loff1;
          const auto idat2 = pdata2 + loff2;
          const auto idat3 = pdata3 + loff3;
          auto odat = p->param.data + loff;

          for (size_t i = 0; i < ngp; ++i)
            {
              const auto ival1 = idat1[(ngp1 > 1) ? i : 0];
              const auto ival2 = idat2[(ngp2 > 1) ? i : 0];
              const auto ival3 = idat3[(ngp3 > 1) ? i : 0];

              if (nmiss1 && DBL_IS_EQUAL(ival1, missval1))
                odat[i] = missval1;
              else if (IS_NOT_EQUAL(ival1, 0.0))
                odat[i] = DBL_IS_EQUAL(ival2, missval2) ? missval1 : ival2;
              else
                odat[i] = DBL_IS_EQUAL(ival3, missval3) ? missval1 : ival3;
            }

          nmiss += array_num_mv(ngp, odat, missval1);
        }

      p->param.nmiss = nmiss;
    }

  if (p1->ltmpobj) node_delete(p1);
  if (p2->ltmpobj) node_delete(p2);
  if (p3->ltmpobj) node_delete(p3);

  return p;
}
/*
static int
exNode(nodeType *p, parseParamType *parse_arg)
{
  if (!p) return 0;

  // node is leaf
  if (p->type == typeCon || p->type == typeVar || p->u.opr.nops == 0) return 0;

  // node has children
  for (int k = 0; k < p->u.opr.nops; k++) exNode(p->u.opr.op[k], parse_arg);

  return 0;
}
*/

static int
param_search_name(const int nparam, const std::vector<paramType> &params, const char *name)
{
  for (int varID = 0; varID < nparam; ++varID)
    {
      if (params[varID].isValid && cdo_cmpstr(params[varID].name, name)) return varID;
    }

  return -1;
}

static int
param_search_name_size(const int nparam, std::vector<paramType> &params, const char *name, size_t ngp, size_t nlev)
{
  for (int varID = 0; varID < nparam; ++varID)
    {
      if (cdo_cmpstr(params[varID].name, name))
        {
          if (ngp == params[varID].ngp && nlev == params[varID].nlev)
            return varID;
          else if (ngp == 0 && nlev == 0)
            return varID;
          else
            params[varID].isValid = false;
        }
    }

  return -1;
}

static void
param_print(const char *vname, const paramType &param, const long tsID)
{
  constexpr size_t maxout = 100;
  const auto data = param.data;
  for (size_t k = 0; k < param.nlev; ++k)
    for (size_t i = 0; i < param.ngp; ++i)
      {
        if (i < maxout || (param.ngp > maxout && i >= (param.ngp - maxout)))
          {
            if (param.steptype == TIME_CONSTANT)
              fprintf(stdout, "   %s[lev=%zu:gp=%zu] = %g\n", vname, k + 1, i + 1, data[k * param.ngp + i]);
            else
              fprintf(stdout, "   %s[ts=%ld:lev=%zu:gp=%zu] = %g\n", vname, tsID, k + 1, i + 1, data[k * param.ngp + i]);
          }
        else if (i == maxout)
          {
            fprintf(stdout, "   .......\n");
          }
      }
}

static void
add_new_constant(const char *varname, parseParamType *parse_arg, std::vector<paramType> &params, const paramType &param)
{
  const auto varID = parse_arg->nparams;
  if (varID >= parse_arg->maxparams) cdo_abort("Too many parameter (limit=%d)", parse_arg->maxparams);

  param_meta_copy(params[varID], param);
  params[varID].type = PARAM_CONST;
  params[varID].isValid = true;
  params[varID].ngp = 0;
  params[varID].nlat = 0;
  params[varID].nlev = 0;
  params[varID].missval = -9.e33;
  params[varID].nmiss = 0;
  params[varID].name = strdup(varname);
  parse_arg->nparams++;
  parse_arg->cnparams++;
}

static void
add_new_param(const char *varname, parseParamType *parse_arg, std::vector<paramType> &params, const paramType &param)
{
  const auto varID = parse_arg->nparams;
  if (varID >= parse_arg->maxparams) cdo_abort("Too many parameter (limit=%d)", parse_arg->maxparams);

  param_meta_copy(params[varID], param);
  params[varID].isValid = true;
  params[varID].hasMissvals = param.hasMissvals;
  params[varID].name = strdup(varname);
  params[varID].nmiss = param.nmiss;
  if (param.units) params[varID].units = strdup(param.units);
  if (param.longname) params[varID].longname = strdup(param.longname);
  parse_arg->nparams++;
  parse_arg->cnparams++;
}

static nodeType *
expr_run_type_com(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;
  auto &params = parse_arg->params;

  const auto cname = p->u.com.cname;
  const auto vname = p->u.com.vname;
  if (parse_arg->debug) cdo_print("\tstatement\t\t%s(%s)", cname, vname);

  auto varID = param_search_name(parse_arg->nparams, params, vname);
  if (varID == -1) cdo_abort("Variable %s not found, needed for statement %s(%s)!", vname, cname, vname);

  if (init)
    {
      if (cdo_cmpstr(cname, "remove")) params[varID].remove = true;
    }
  else
    {
      if (cdo_cmpstr(cname, "print")) param_print(vname, params[varID], std::lround(params[parse_arg->tsID].data[0]));
    }

  return nullptr;
}

static nodeType *
expr_run_type_con(nodeType *p, parseParamType *parse_arg)
{
  if (parse_arg->debug) cdo_print("\tpush\tconst\t%g", p->u.con.value);
  return p;
}

static int
expr_run_type_var_grid(const char *vnm, int coord, parseParamType *parse_arg)
{
  auto &params = parse_arg->params;

  const auto len = strlen(vnm);
  auto varname = strdup(vnm);
  varname[len - 2] = 0;
  const auto varID = param_search_name(parse_arg->nparams, params, varname);
  if (varID == -1) cdo_abort("Coordinate %c: variable >%s< not found!", coord, varname);

  const auto nvarID = parse_arg->nparams;
  if (nvarID >= parse_arg->maxparams) cdo_abort("Too many parameter (limit=%d)", parse_arg->maxparams);

  const auto coordID = params_get_coord_ID(parse_arg, coord, params[varID].gridID);
  parse_arg->coords[coordID].needed = true;
  const auto units = parse_arg->coords[coordID].units.data();
  const auto longname = parse_arg->coords[coordID].longname.data();

  params[nvarID].isValid = true;
  params[nvarID].coord = coord;
  params[nvarID].hasMissvals = false;
  params[nvarID].name = strdup(vnm);
  params[nvarID].missval = params[varID].missval;
  params[nvarID].gridID = params[varID].gridID;
  params[nvarID].zaxisID = parse_arg->surfaceID;
  params[nvarID].steptype = TIME_CONSTANT;
  params[nvarID].ngp = params[varID].ngp;
  params[nvarID].nlat = params[varID].nlat;
  params[nvarID].nlev = 1;
  if (units) params[nvarID].units = strdup(units);
  if (longname) params[nvarID].longname = strdup(longname);
  parse_arg->nparams++;
  parse_arg->cnparams++;

  free(varname);

  return nvarID;
}

static int
expr_run_type_var_zaxis(const char *vnm, int coord, parseParamType *parse_arg)
{
  auto &params = parse_arg->params;

  const auto len = strlen(vnm);
  auto varname = strdup(vnm);
  varname[len - 2] = 0;
  const auto varID = param_search_name(parse_arg->nparams, params, varname);
  if (varID == -1) cdo_abort("Coordinate %c: variable >%s< not found!", coord, varname);

  const auto nvarID = parse_arg->nparams;
  if (nvarID >= parse_arg->maxparams) cdo_abort("Too many parameter (limit=%d)", parse_arg->maxparams);

  const auto coordID = params_get_coord_ID(parse_arg, coord, params[varID].zaxisID);
  parse_arg->coords[coordID].needed = true;
  const auto units = parse_arg->coords[coordID].units.data();
  const auto longname = parse_arg->coords[coordID].longname.data();

  params[nvarID].isValid = true;
  params[nvarID].coord = coord;
  params[nvarID].hasMissvals = false;
  params[nvarID].name = strdup(vnm);
  params[nvarID].missval = params[varID].missval;
  params[nvarID].gridID = parse_arg->pointID;
  params[nvarID].zaxisID = params[varID].zaxisID;
  params[nvarID].steptype = TIME_CONSTANT;
  params[nvarID].ngp = 1;
  params[nvarID].nlev = params[varID].nlev;
  if (units) params[nvarID].units = strdup(units);
  if (longname) params[nvarID].longname = strdup(longname);
  parse_arg->nparams++;
  parse_arg->cnparams++;

  free(varname);

  return nvarID;
}

static nodeType *
expr_run_type_var(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;
  const auto &params = parse_arg->params;

  const auto vnm = p->u.var.nm;
  auto varID = param_search_name(parse_arg->nparams, params, vnm);
  if (varID == -1 && init)
    {
      const auto len = strlen(vnm);
      if (len > 2 && vnm[len - 2] == '.')
        {
          const auto coord = vnm[len - 1];
          if (coord == 'x' || coord == 'y' || coord == 'a' || coord == 'w')
            {
              varID = expr_run_type_var_grid(vnm, coord, parse_arg);
            }
          else if (coord == 'z' || coord == 'd')
            {
              varID = expr_run_type_var_zaxis(vnm, coord, parse_arg);
            }
        }
    }

  if (varID == -1)
    {
      cdo_abort("Variable >%s< not found!", vnm);
    }
  else if (init)
    {
      if (varID < parse_arg->nvars1 && !parse_arg->needed[varID]) parse_arg->needed[varID] = true;
    }

  param_meta_copy(p->param, params[varID]);
  p->param.coord = params[varID].coord;
  p->param.hasMissvals = params[varID].hasMissvals;
  p->param.name = params[varID].name;
  p->param.longname = params[varID].longname;
  p->param.units = params[varID].units;
  p->ltmpobj = false;

  if (!init)
    {
      p->param.data = params[varID].data;
      p->param.nmiss = params[varID].nmiss;
    }

  if (parse_arg->debug)
    {
      if (p->param.ngp == 0 && p->param.nlev == 0)
        cdo_print("\tpush\tvar\t%s", vnm);
      else
        cdo_print("\tpush\tvar\t%s[N%zu][L%zu]", vnm, p->param.ngp, p->param.nlev);
    }

  return p;
}

static nodeType *
expr_run_type_fun1c(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;

  const auto funcID = get_funcID(p->u.fun1c.name);
  const auto functype = fun_sym_tbl[funcID].type;

  if (functype == FT_1C)
    {
      auto fnode = expr_run(p->u.fun1c.op, parse_arg);
      const auto value = p->u.fun1c.value;
      return fun1c(init, funcID, fnode, value, parse_arg);
    }
  else
    {
      cdo_abort("Syntax error in call to %s(), check number of parameter!\n", p->u.fun1c.name);
    }

  return nullptr;
}

static nodeType *
expr_run_type_fun2c(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;

  const auto funcID = get_funcID(p->u.fun2c.name);
  const auto functype = fun_sym_tbl[funcID].type;

  if (functype == FT_2C)
    {
      auto fnode = expr_run(p->u.fun2c.op, parse_arg);
      const auto value1 = p->u.fun2c.value1;
      const auto value2 = p->u.fun2c.value2;
      return fun2c(init, funcID, fnode, value1, value2, parse_arg);
    }
  else
    {
      cdo_abort("Syntax error in call to %s(), check number of parameter!\n", p->u.fun2c.name);
    }

  return nullptr;
}

static nodeType *
expr_run_type_fun(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;
  const auto &params = parse_arg->params;

  const auto funcID = get_funcID(p->u.fun.name);
  auto functype = fun_sym_tbl[funcID].type;

  auto fnode = expr_run(p->u.fun.op, parse_arg);

  if (functype == FT_COORD)
    {
      return coord_fun(init, funcID, fnode, parse_arg);
    }
  else if (functype == FT_0)
    {
      const auto vartsID = parse_arg->tsID;
      return ex_fun0_con(init, funcID, params[vartsID].data);
    }
  else
    {
      functype = fun_sym_tbl[funcID].type;
      const auto funcflag = fun_sym_tbl[funcID].flag;
      if (functype == FT_FLD && funcflag == 1)
        {
          const auto coordID = params_get_coord_ID(parse_arg, 'w', fnode->param.gridID);
          if (init)
            parse_arg->coords[coordID].needed = true;
          else
            fnode->param.weight = parse_arg->coords[coordID].data.data();
        }
      auto rnode = ex_fun(init, funcID, fnode);
      // if (fnode->ltmpobj) node_delete(fnode);
      // Free(fnode);
      return rnode;
    }
}

static nodeType *
expr_run_type_opr(nodeType *p, parseParamType *parse_arg)
{
  const auto init = parse_arg->init;
  auto &params = parse_arg->params;

  // clang-format off
  switch (p->u.opr.oper)
    {
    case '=':
      {
        auto rnode = expr_run(p->u.opr.op[1], parse_arg);

        const auto varname2 = p->u.opr.op[0]->u.var.nm;

        if (parse_arg->debug)
          {
            if (rnode && rnode->type == typeVar && rnode->param.ngp > 0 && rnode->param.nlev > 0)
              cdo_print("\tpop\tvar\t%s[N%zu][L%zu]", varname2, rnode->param.ngp, rnode->param.nlev);
            else
              cdo_print("\tpop\tconst\t%s", varname2);
          }

        // auto varID = param_search_name(parse_arg->nparams, params, varname2);
        auto varID = param_search_name_size(parse_arg->nparams, params, varname2, rnode->param.ngp, rnode->param.nlev);
        if (init)
          {
            if (varID >= 0)
              {
                // printf("  found %s\n", varname2);
                if (varID < parse_arg->nvars1)
                  {
                    if (rnode->param.nlev > 0 && params[varID].nlev != rnode->param.nlev)
                      cdo_abort("The number of layers must not change (name=%s layers: in=%zu out=%zu)!", params[varID].name,
                               params[varID].nlev, rnode->param.nlev);

                    params[varID].select = true;
                    parse_arg->needed[varID] = true;
                  }
                else if (params[varID].coord)
                  cdo_abort("Coordinate variable %s is read only!", varname2);
                /*
                  else
                  cdo_warning("Variable %s already defined!", varname2);
                */
              }
            else if (rnode && rnode->type == typeCon)
              {
                add_new_constant(varname2, parse_arg, params, rnode->param);
              }
            else if (p->u.opr.op[1]->type != typeCon)
              {
                add_new_param(varname2, parse_arg, params, rnode->param);
              }
          }
        else
          {
            if (varID < 0) cdo_abort("Variable >%s< not found!", varname2);

            if (params[varID].coord) cdo_abort("Coordinate variable %s is read only!", varname2);
            param_meta_copy(p->param, params[varID]);
            p->param.name = params[varID].name;
            p->param.data = params[varID].data;
            p->ltmpobj = false;

            ex_copy(init, p, rnode);
            params[varID].nmiss = p->param.nmiss;
          }

        if (rnode && rnode->ltmpobj)
          {
            node_delete(rnode);
            rnode = nullptr;
          }
        // else Free(rnode);
        break;
      }
    case UMINUS: return ex_uminus(init, expr_run(p->u.opr.op[0], parse_arg));
    case NOT:    return ex_not(init, expr_run(p->u.opr.op[0], parse_arg));
    case '?':    return ex_ifelse(init, expr_run(p->u.opr.op[0], parse_arg), expr_run(p->u.opr.op[1], parse_arg),
                                  expr_run(p->u.opr.op[2], parse_arg));
    default:     return expr(init, p->u.opr.oper, expr_run(p->u.opr.op[0], parse_arg), expr_run(p->u.opr.op[1], parse_arg));
    }
  // clang-format on

  return nullptr;
}

nodeType *
expr_run(nodeType *p, parseParamType *parse_arg)
{
  pointID = parse_arg->pointID;
  zonalID = parse_arg->zonalID;
  surfaceID = parse_arg->surfaceID;

  if (!p) return nullptr;

  // clang-format off
  switch (p->type)
    {
    case typeCom:     return expr_run_type_com(p, parse_arg);
    case typeCon:     return expr_run_type_con(p, parse_arg);
    case typeVar:     return expr_run_type_var(p, parse_arg);
    case typeFun1c:   return expr_run_type_fun1c(p, parse_arg);
    case typeFun2c:   return expr_run_type_fun2c(p, parse_arg);
    case typeFun:     return expr_run_type_fun(p, parse_arg);
    case typeOpr:     return expr_run_type_opr(p, parse_arg);
    }
  // clang-format on

  return nullptr;
}