Add extended cubed-sphere dual halo with a correct alpha-beta projection.

src/atlas/CMakeLists.txt

@@ -710,6 +710,7 @@ util/KDTree.cc
 util/KDTree.h
 util/PolygonXY.cc
 util/PolygonXY.h
+util/Matrix.h
 util/Metadata.cc
 util/Metadata.h
 util/Point.cc
@@ -908,6 +909,7 @@ target_link_libraries( atlas PUBLIC
   eckit
   eckit_geometry
   eckit_linalg
+  eckit_maths
   eckit_mpi
   eckit_option
 )

src/atlas/functionspace/CubedSphereColumns.cc

@@ -150,6 +150,11 @@ idx_t CubedSphereColumns<BaseFunctionSpace>::index(idx_t t, idx_t i, idx_t j) co
     return cubedSphereColumnsHandle_.get()->index(t, i, j);
 }
 
+template <typename BaseFunctionSpace>
+bool CubedSphereColumns<BaseFunctionSpace>::is_valid_index(idx_t t, idx_t i, idx_t j) const {
+    return cubedSphereColumnsHandle_.get()->is_valid_index(t, i, j);
+}
+
 template <typename BaseFunctionSpace>
 Field CubedSphereColumns<BaseFunctionSpace>::tij() const {
     return cubedSphereColumnsHandle_.get()->tij();

src/atlas/functionspace/CubedSphereColumns.h

@@ -54,6 +54,9 @@ class CubedSphereColumns : public BaseFunctionSpace {
     /// Return array_view index for (t, i, j).
     idx_t index(idx_t t, idx_t i, idx_t j) const;
 
+    /// Return true if (t, i, j) is a valid index.
+    bool is_valid_index(idx_t t, idx_t i, idx_t j) const;
+
     /// Return tij field.
     Field tij() const;
 

src/atlas/functionspace/detail/CubedSphereStructure.cc

@@ -120,6 +120,29 @@ idx_t CubedSphereStructure::index(idx_t t, idx_t i, idx_t j) const {
     return tijToIdx_[static_cast<size_t>(t)][vecIndex(t, i, j)];
 }
 
+bool CubedSphereStructure::is_valid_index(idx_t t, idx_t i, idx_t j) const {
+
+
+    // Check if t is in range.
+    if (t < 0 || t > 5) {
+        return false;
+    }
+
+    // Check if i and j are in range in index method.
+    if ( i < i_begin(t) || i >= i_end(t) ||
+         j < j_begin(t) || j >= j_end(t)) {
+        return false;
+    }
+
+    // Check if (t, i, j) is a valid index.
+    if (index(t, i, j) == invalid_index()) {
+        return false;
+    }
+
+    return true;
+
+}
+
 Field CubedSphereStructure::tij() const {
     return tij_;
 }

src/atlas/functionspace/detail/CubedSphereStructure.h

@@ -53,6 +53,10 @@ class CubedSphereStructure : public util::Object {
     /// Return array_view index for (t, i, j).
     idx_t index(idx_t t, idx_t i, idx_t j) const;
 
+    /// Return true if (t, i, j) is a valid index.
+    bool is_valid_index(idx_t t, idx_t i, idx_t j) const;
+
+
     /// Return ijt field.
     Field tij() const;
 

src/atlas/grid/CubedSphereGrid.h

@@ -295,9 +295,19 @@ class CubedSphereGrid : public Grid {
     // Return the size of the cubed sphere grid, where N is the number of grid boxes along the edge of a tile
     inline int N() const { return grid_->N(); }
 
-    // Return the number of tiles
+    /// @brief return tiles object.
     inline atlas::grid::CubedSphereTiles tiles() const { return grid_->tiles(); }
 
+    /// @brief return cubed sphere projection object.
+    inline const projection::detail::CubedSphereProjectionBase& cubedSphereProjection() const {
+
+        const auto projPtr =
+            dynamic_cast<const projection::detail::CubedSphereProjectionBase*>
+            (projection().get());
+
+        return *projPtr;
+    };
+
     temporary::IterateTIJ tij() const { return temporary::IterateTIJ(*grid_); }
 
     const std::string& stagger() const { return grid_->stagger(); }

src/atlas/grid/Tiles.cc

@@ -72,5 +72,13 @@ std::ostream& operator<<(std::ostream& os, const CubedSphereTiles& t) {
     return os;
 }
 
+const PointXY& CubedSphereTiles::tileCentre(size_t t) const {
+    return get()->tileCentre(t);
+}
+
+const util::Matrix22& CubedSphereTiles::tileJacobian(size_t t) const {
+    return get()->tileJacobian(t);
+}
+
 }  // namespace grid
 }  // namespace atlas

src/atlas/grid/Tiles.h

@@ -16,6 +16,7 @@
 
 #include "atlas/library/config.h"
 #include "atlas/util/ObjectHandle.h"
+#include "atlas/util/Matrix.h"
 
 //---------------------------------------------------------------------------------------------------------------------
 
@@ -85,6 +86,13 @@ class CubedSphereTiles : DOXYGEN_HIDE(public util::ObjectHandle<atlas::grid::det
     // a "diagonal" extension over corners of the cube.
     atlas::PointXY tileCubePeriodicity(const atlas::PointXY& xyExtended, const atlas::idx_t tile) const;
 
+    /// @brief Return the position of the tile centre in xy space.
+    const PointXY& tileCentre(size_t t) const;
+
+    /// @brief Return the Jacobian of xy with respect to the curvilinear
+    ///        coordinates of the tile.
+    const util::Matrix22& tileJacobian(size_t t) const;
+
 private:
     /// Output to stream
     void print(std::ostream&) const;

src/atlas/grid/detail/tiles/FV3Tiles.cc

@@ -453,6 +453,15 @@ void FV3CubedSphereTiles::print(std::ostream& os) const {
 }
 
 
+const PointXY& FV3CubedSphereTiles::tileCentre(size_t t) const {
+    throw_NotImplemented("tileCentre not implemented for FV3Tiles", Here());
+}
+
+const util::Matrix22& FV3CubedSphereTiles::tileJacobian(size_t t) const {
+    throw_NotImplemented("tileJacobian not implemented for FV3Tiles", Here());;
+}
+
+
 namespace {
 static CubedSphereTilesBuilder<FV3CubedSphereTiles> register_builder(FV3CubedSphereTiles::static_type());
 }

src/atlas/grid/detail/tiles/FV3Tiles.h

@@ -46,6 +46,11 @@ class FV3CubedSphereTiles : public CubedSphereTiles {
 
     virtual void print(std::ostream&) const override;
 
+    virtual const PointXY& tileCentre(size_t t) const override;
+
+    virtual const util::Matrix22& tileJacobian(size_t t) const override;
+
+
 private:
 };