Merge branch 'release/1.24.0'

CMakeLists.txt

@@ -21,7 +21,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 ecbuild_add_option(FEATURE BUILD_TOOLS        DEFAULT OFF DESCRIPTION "build the command line tools (install)")
 ecbuild_add_option(FEATURE MIR_DOWNLOAD_MASKS DEFAULT ON  DESCRIPTION "download LSM files")
-ecbuild_add_option(FEATURE NETCDF             DEFAULT ON  DESCRIPTION "support for netCDF" REQUIRED_PACKAGES "NetCDF COMPONENTS C")
+ecbuild_add_option(FEATURE NETCDF             DEFAULT ON  DESCRIPTION "support for netCDF" REQUIRED_PACKAGES "NetCDF COMPONENTS C CXX")
 ecbuild_add_option(FEATURE PNG                DEFAULT OFF DESCRIPTION "support for PNG encoding" REQUIRED_PACKAGES PNG)
 ecbuild_add_option(FEATURE OMP                DEFAULT OFF DESCRIPTION "support for OpenMP shared memory parallelism" REQUIRED_PACKAGES "OpenMP COMPONENTS CXX")
 

VERSION

@@ -1 +1 @@
-1.23.0
+1.24.0

etc/mir/grib-input.yaml

@@ -19,6 +19,10 @@ md5GridSection=f5dc74ec36353f4c83f7de3bf46e1aef, centre:i=98, gridName=O640, edi
 md5GridSection=6ae0076681a2541abe053ab52fa60f30, centre:i=84/85, gridType=lambert, edition=1:
   - first_point_bottom_left: true
 
+# gridType=polar_stereographic does not support the default area method
+gridType=polar_stereographic:
+  - area-mode: mask
+
 # gridType=space_view does not support the default interpolation method
 gridType=space_view:
   - interpolation: k-nearest

etc/mir/parameter-class.yaml

@@ -1952,6 +1952,7 @@ intensity:
   - 261020  # Physiological equivalent temperature
   - 261021  # Saturation water vapor pressure
   - 261022  # Wet-bulb potential temperature
+  - 261023  # Wet-bulb temperature
   - 262000  # Sea ice thickness
   - 262002  # Snow thickness over sea ice
   - 262003  # Eastward sea ice velocity

src/mir/CMakeLists.txt

@@ -41,6 +41,8 @@ list(APPEND mir_srcs
     action/context/Context.h
     action/filter/AddRandomFilter.cc
     action/filter/AddRandomFilter.h
+    action/filter/AdjustWindsAtPoles.cc
+    action/filter/AdjustWindsAtPoles.h
     action/filter/AdjustWindsDirections.cc
     action/filter/AdjustWindsDirections.h
     action/filter/AdjustWindsScaleCosLatitude.cc
@@ -246,6 +248,8 @@ list(APPEND mir_srcs
     input/GribMemoryInput.h
     input/GribStreamInput.cc
     input/GribStreamInput.h
+    input/GriddefInput.cc
+    input/GriddefInput.h
     input/MIRInput.cc
     input/MIRInput.h
     input/MultiDimensionalGribFileInput.cc
@@ -457,6 +461,8 @@ list(APPEND mir_srcs
     output/GribOutput.h
     output/GribStreamOutput.cc
     output/GribStreamOutput.h
+    output/GriddefOutput.cc
+    output/GriddefOutput.h
     output/MIROutput.cc
     output/MIROutput.h
     output/MultiDimensionalOutput.cc

src/mir/action/area/Area.h

@@ -38,7 +38,7 @@ class Area : public Action {
 
     // -- Destructor
 
-    virtual ~Area() override;
+    ~Area() override;
 
     // -- Convertors
     // None
@@ -72,7 +72,7 @@ class Area : public Action {
 
     // -- Overridden methods
 
-    virtual bool sameAs(const Action&) const override;
+    bool sameAs(const Action&) const override;
 
     // -- Class members
     // None

src/mir/action/filter/AdjustWindsAtPoles.cc

@@ -0,0 +1,114 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#include "mir/action/filter/AdjustWindsAtPoles.h"
+
+#include <limits>
+#include <memory>
+#include <ostream>
+
+#include "mir/action/context/Context.h"
+#include "mir/data/MIRField.h"
+#include "mir/repres/Iterator.h"
+#include "mir/repres/Representation.h"
+#include "mir/util/Exceptions.h"
+
+
+namespace mir::action {
+
+
+bool AdjustWindsAtPoles::sameAs(const Action& other) const {
+    const auto* o = dynamic_cast<const AdjustWindsAtPoles*>(&other);
+    return (o != nullptr);
+}
+
+
+void AdjustWindsAtPoles::print(std::ostream& out) const {
+    out << "AdjustWindsAtPoles[]";
+}
+
+
+void AdjustWindsAtPoles::execute(context::Context& ctx) const {
+    auto& field = ctx.field();
+
+    ASSERT(field.dimensions() > 0);
+
+    const auto missingValue = field.hasMissing() ? field.missingValue() : std::numeric_limits<double>::quiet_NaN();
+
+    const auto N = field.values(0).size();
+    ASSERT(N > 0);
+
+    // find points at pole(s)
+    const auto* representation(field.representation());
+    ASSERT(representation);
+
+    std::vector<size_t> north;
+    std::vector<size_t> south;
+
+    for (const std::unique_ptr<repres::Iterator> it(representation->iterator()); it->next();) {
+        const auto lat = it->pointUnrotated().lat();
+        if (lat == Latitude::NORTH_POLE) {
+            north.push_back(it->index());
+        }
+        else if (lat == Latitude::SOUTH_POLE) {
+            south.push_back(it->index());
+        }
+    }
+
+    if (field.dimensions() % 2 == 0) {
+        // two-component mode: set v to zero at pole(s)
+        for (size_t i = 0; i < field.dimensions(); i += 2) {
+            auto& u = field.direct(i);
+            auto& v = field.direct(i + 1);
+            ASSERT(u.size() == N);
+            ASSERT(v.size() == N);
+
+            for (const auto& indices : {north, south}) {
+                for (const auto& idx : indices) {
+                    if (u[idx] == missingValue || v[idx] == missingValue) {
+                        u[idx] = missingValue;
+                        v[idx] = missingValue;
+                    }
+                    else {
+                        v[idx] = 0;
+                    }
+                }
+            }
+        }
+    }
+    else {
+        // single-component mode: set u, v to zero at pole(s)
+        for (size_t i = 0; i < field.dimensions(); ++i) {
+            auto& v = field.direct(i);
+            ASSERT(v.size() == N);
+
+            for (const auto& indices : {north, south}) {
+                for (const auto& idx : indices) {
+                    if (v[idx] != missingValue) {
+                        v[idx] = 0.;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+const char* AdjustWindsAtPoles::name() const {
+    return "AdjustWindsAtPoles";
+}
+
+
+static const ActionBuilder<AdjustWindsAtPoles> __action("filter.adjust-winds-at-poles");
+
+
+}  // namespace mir::action

src/mir/action/filter/AdjustWindsAtPoles.h

@@ -0,0 +1,37 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#pragma once
+
+#include "mir/action/plan/Action.h"
+
+
+namespace mir::action {
+
+
+class AdjustWindsAtPoles : public Action {
+public:
+    // -- Constructors
+
+    using Action::Action;
+
+private:
+    // -- Overridden methods
+
+    void execute(context::Context&) const override;
+    bool sameAs(const Action&) const override;
+    const char* name() const override;
+    void print(std::ostream&) const override;
+};
+
+
+}  // namespace mir::action

src/mir/action/filter/AdjustWindsScaleCosLatitude.cc

@@ -14,7 +14,6 @@
 
 #include <algorithm>
 #include <cmath>
-#include <limits>
 #include <memory>
 #include <ostream>
 
@@ -45,35 +44,33 @@ void AdjustWindsScaleCosLatitude::print(std::ostream& out) const {
 
 
 void AdjustWindsScaleCosLatitude::execute(context::Context& ctx) const {
-    data::MIRField& field = ctx.field();
-    ASSERT(field.dimensions() > 0);
+    auto& field = ctx.field();
 
+    ASSERT(field.dimensions() > 0);
+    ASSERT(!field.hasMissing());
 
     // determine scaling vector (all fields share the same representation)
     const size_t N = field.values(0).size();
     ASSERT(N > 0);
 
-    const repres::Representation* representation(field.representation());
+    const auto* representation(field.representation());
     ASSERT(representation);
 
-    std::vector<double> scale(N, std::numeric_limits<double>::quiet_NaN());
+    std::vector<double> scale(N);
     for (const std::unique_ptr<repres::Iterator> it(representation->iterator()); it->next();) {
-        const auto& p      = it->pointUnrotated();
-        scale[it->index()] = (p.lat() == Latitude::SOUTH_POLE) ? 0.
-                             : (p.lat() == Latitude::NORTH_POLE)
+        const auto lat     = it->pointUnrotated().lat();
+        scale[it->index()] = (lat == Latitude::SOUTH_POLE || lat == Latitude::NORTH_POLE)
                                  ? 0.
-                                 : 1. / std::cos(util::degree_to_radian(p.lat().value()));
+                                 : 1. / std::cos(util::degree_to_radian(lat.value()));
     }
 
     // apply scaling to each field component directly
     for (size_t i = 0; i < field.dimensions(); ++i) {
         auto& values = field.direct(i);
         ASSERT(values.size() == N);
 
-        std::transform(
-            values.begin(), values.end(), scale.begin(), values.begin(), [](const double& a, const double& b) {
-                return std::isfinite(a) && std::isfinite(b) ? a * b : std::numeric_limits<double>::quiet_NaN();
-            });
+        std::transform(values.begin(), values.end(), scale.begin(), values.begin(),
+                       [](const double& a, const double& b) { return a * b; });
     }
 }
 

src/mir/action/io/EndAction.h

@@ -35,7 +35,7 @@ class EndAction : public Action {
 
     // -- Destructor
 
-    virtual ~EndAction() override;
+    ~EndAction() override;
 
     // -- Convertors
     // None
@@ -87,8 +87,8 @@ class EndAction : public Action {
     bool sameAs(const Action&) const override;
     bool isEndAction() const override;
 
-    virtual void print(std::ostream&) const override;
-    virtual void custom(std::ostream&) const override;
+    void print(std::ostream&) const override;
+    void custom(std::ostream&) const override;
 
     // -- Class members
     // None

src/mir/api/mir_config.h.in

@@ -1,11 +1,36 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
 #pragma once
 
 #include "mir/mir_ecbuild_config.h"
 
 #include "mir/api/mir_version.h"
 
+
 #cmakedefine01 mir_HAVE_ATLAS
 #cmakedefine01 mir_HAVE_NETCDF
 #cmakedefine01 mir_HAVE_PNG
 #cmakedefine01 mir_HAVE_OMP
 
+
+constexpr bool _to_bool(int x = 0) {
+    return x != 0;
+}
+
+constexpr bool HAVE_ATLAS  = _to_bool(@mir_HAVE_ATLAS@);
+constexpr bool HAVE_NETCDF = _to_bool(@mir_HAVE_NETCDF@);
+constexpr bool HAVE_PNG    = _to_bool(@mir_HAVE_PNG@);
+constexpr bool HAVE_OMP    = _to_bool(@mir_HAVE_OMP@);
+
+constexpr bool HAVE_TESSELATION = _to_bool(@atlas_HAVE_TESSELATION@);
+constexpr bool HAVE_PROJ        = _to_bool(@atlas_HAVE_PROJ@);

src/mir/caching/LegendreCache.cc

@@ -12,6 +12,8 @@
 
 #include "mir/caching/LegendreCache.h"
 
+#include "eckit/config/Resource.h"
+
 #include "mir/config/LibMir.h"
 
 

src/mir/caching/MeshCache.cc

@@ -12,6 +12,8 @@
 
 #include "mir/caching/MeshCache.h"
 
+#include "eckit/config/Resource.h"
+
 #include "mir/config/LibMir.h"
 
 

src/mir/caching/WeightCache.cc

@@ -14,6 +14,8 @@
 
 #include <sstream>
 
+#include "eckit/config/Resource.h"
+
 #include "mir/caching/matrix/MatrixLoader.h"
 #include "mir/config/LibMir.h"
 #include "mir/method/WeightMatrix.h"
@@ -81,10 +83,7 @@ void WeightCacheTraits::load(const eckit::CacheManagerBase& manager, value_type&
     value_type tmp(matrix::MatrixLoaderFactory::build(manager.loader(), path));
     w.swap(tmp);
 
-    static bool matrixValidate = eckit::Resource<bool>("$MIR_MATRIX_VALIDATE", false);
-    if (matrixValidate) {
-        w.validate("fromCache");
-    }
+    w.validate("fromCache");  // check matrix structure (only)
 }
 
 

src/mir/caching/legendre/MappedMemoryLoader.h

@@ -30,9 +30,9 @@ class MappedMemoryLoader : public LegendreLoader {
     void print(std::ostream&) const override;
 
 private:
-    virtual const void* address() const override;
-    virtual size_t size() const override;
-    virtual bool inSharedMemory() const override;
+    const void* address() const override;
+    size_t size() const override;
+    bool inSharedMemory() const override;
 
 private:
     int fd_;