diff --git a/src/expr.cc b/src/expr.cc
index 941154b..6e44547 100644
--- a/src/expr.cc
+++ b/src/expr.cc
@@ -49,6 +49,40 @@ static double f_rad(double x) { return x * M_PI / 180.0; }
 static double f_deg(double x) { return x * 180.0 / M_PI; }
 static double f_isMissval(double x) { (void)x; return 0.0; }
 
+static double trimrel(double x, double keepbits) { // https://gmd.copernicus.org/articles/14/377/2021/
+  int maskbits;
+  int64_t  *b;
+  int64_t  mask, half_quantum1;
+  //
+  // We might wish to add something here that will check that 0 <= keepbits <= 52
+  b = (int64_t *) &x;
+  maskbits = 52 - (int) keepbits; // explicit bits in double precision
+  mask = (0xFFFFFFFFFFFFFFFF>>maskbits)<<maskbits;
+  half_quantum1 = ( 1<<(maskbits-1) ) - 1;
+  *b += ((*b>>maskbits) & 1) +half_quantum1;
+  *b &= mask;
+  return x;
+}
+
+static double trimabs(double x, double max_abs_err) {//Sec. A6 https://gmd.copernicus.org/articles/14/377/2021/
+  int maskbits;
+  double quantum;
+  int64_t  *b;
+  const int64_t zerofrac = ~0ULL << 52; // zeros in place of mantissa in IEEE754 double
+  const double maxval = double(1ULL << 52); // max quantizable value in quanta
+  //
+  // We might wish to add something here that will check that max_abs_err > 0
+  //
+  if (abs(x) < maxval*max_abs_err) { //quantisation is meaningful, no integer overflow
+    b = (int64_t *) (& max_abs_err); // access it as int64
+    *b &= zerofrac; // zero fraction bits of xmax_abs_err: reduce it to benefit from compression
+    quantum = max_abs_err * 2; //quantum -- twice of that
+    x = nearbyint (x / quantum) * quantum;
+  }
+
+  return x;
+}
+
 static double pt_definegrid(const ParamEntry *const p) { printf("pt_definegrid: >%s<\n", p->name.c_str()); return -1;/*return cdo_define_grid(xxx);*/ }
 static double pt_grid(const ParamEntry *const p) { return p->gridID; }
 static double pt_ngp(const ParamEntry *const p) { return p->ngp; }
@@ -112,6 +146,8 @@ static const ExprFuncEntry funcSymbolTable[] = {
   { FT_STD, 0, "pow",   (void (*)(void))((double (*)(double, double)) std::pow) },
   { FT_STD, 0, "hypot", (void (*)(void))((double (*)(double, double)) std::hypot) },
   { FT_STD, 0, "atan2", (void (*)(void))((double (*)(double, double)) std::atan2) },
+  { FT_STD, 0, "trimrel", (void (*)(void))((double (*)(double, double)) (&trimrel)) },
+  { FT_STD, 0, "trimabs", (void (*)(void))((double (*)(double, double)) (&trimabs)) },
   // scalar functions   y=func(x)
   { FT_STD, 0, "float",     reinterpret_cast<void (*)(void)>(&f_float) },
   { FT_STD, 0, "int",       reinterpret_cast<void (*)(void)>(&f_int) },
diff --git a/src/operator_help.cc b/src/operator_help.cc
index 1533aad..149a485 100644
--- a/src/operator_help.cc
+++ b/src/operator_help.cc
@@ -1989,6 +1989,13 @@ const CdoHelp ExprHelp = {
     "    pow(x,y)    "    "    Power function",
     "    hypot(x,y)  "    "    Euclidean distance function, sqrt(x*x + y*y)",
     "    atan2(x,y)  "    "    Arc tangent function of y/x, using signs to determine quadrants",
+    "   trimrel(x,kb)"    "    Trim relative precision to kb keep-bits. Max relative error to a value",
+    "                "    "    | result - x | / x   < 2^(-1-kb) for any finite x. ",
+    "                "    "    trimrel is a non-decreasing funciorn of x ",
+    "                "    "    Loosely follows A5 of https://gmd.copernicus.org/articles/14/377/2021/",
+    "  trimabs(x,err)"    "    Trim absolute precision introducing given max. absolute error",
+    "                "    "    | result - x |  < err, trimabs is a non-decreasing funciorn of x",
+    "                "    "    Loosely follows A6 of https://gmd.copernicus.org/articles/14/377/2021/",
     "    ",
     "    Coordinates:",
     "    ",