diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6e19fee71..09676f01e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,10 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
+## [0.32.1] - 2023-02-09
+### Added
+- Added (lon, lat) to (alpha, beta) transforms to cubed sphere projection
+
 ## [0.32.0] - 2023-01-23
 ### Added
 - Added BlockStructuredColumns FunctionSpace
@@ -428,6 +432,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[0.32.1]: https://github.com/ecmwf/atlas/compare/0.32.0...0.32.1
 [0.32.0]: https://github.com/ecmwf/atlas/compare/0.31.1...0.32.0
 [0.31.1]: https://github.com/ecmwf/atlas/compare/0.31.0...0.31.1
 [0.31.0]: https://github.com/ecmwf/atlas/compare/0.30.0...0.31.0
diff --git a/VERSION b/VERSION
index 9eb2aa3f1..fd9620c08 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-0.32.0
+0.32.1
diff --git a/src/atlas/interpolation/method/cubedsphere/CellFinder.cc b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc
index 2c9c4477b..c4c1e4d09 100644
--- a/src/atlas/interpolation/method/cubedsphere/CellFinder.cc
+++ b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc
@@ -26,105 +26,100 @@ CellFinder::CellFinder(const Mesh& mesh, const util::Config& config): mesh_{mesh
     projection_ = &(csGrid.cubedSphereProjection());
 
     // Get views to cell data.
-    const auto xyView   = array::make_view<double, 2>(mesh_.cells().field("xy"));
-    const auto tijView  = array::make_view<idx_t, 2>(mesh_.cells().field("tij"));
+    const auto lonlatView   = array::make_view<double, 2>(mesh_.cells().field("lonlat"));
     const auto haloView = array::make_view<int, 1>(mesh_.cells().halo());
 
     // make points and payloads vectors.
-    auto points   = std::array<std::vector<Point2>, 6>{};
-    auto payloads = std::array<std::vector<idx_t>, 6>{};
+    auto points   = std::vector<PointLonLat>{};
+    auto payloads = std::vector<idx_t>{};
 
     // Iterate over cells.
     auto halo = config.getInt("halo", 0);
     for (idx_t i = 0; i < mesh_.cells().size(); ++i) {
         if (haloView(i) <= halo) {
-            const auto t  = tijView(i, 0);
-            const auto xy = PointXY(xyView(i, XX), xyView(i, YY));
-
-            points[size_t(t)].push_back(projection_->xy2alphabeta(xy, t));
-            payloads[size_t(t)].push_back(i);
+            points.push_back(PointLonLat(lonlatView(i, LON), lonlatView(i, LAT)));
+            payloads.push_back(i);
         }
     }
 
-    // build cell kd-trees.
-    for (size_t t = 0; t < 6; ++t) {
-        trees_[t].build(points[t], payloads[t]);
-    }
+    // build cell kd-tree.
+    tree_.build(points, payloads);
 }
 
-CellFinder::Cell CellFinder::getCell(const PointXY& xy, size_t listSize, double edgeEpsilon, double epsilon) const {
+CellFinder::Cell CellFinder::getCell(const PointLonLat& lonlat, size_t listSize, double edgeEpsilon, double epsilon) const {
     // Convert xy to alphabeta and t;
     const auto& tiles    = projection_->getCubedSphereTiles();
-    const auto t         = tiles.indexFromXY(xy.data());
-    const auto alphabeta = projection_->xy2alphabeta(xy, t);
-
 
     // Get mesh nodes and connectivity table.
-    const auto nodeXyView        = array::make_view<double, 2>(mesh_.nodes().xy());
-    const auto& nodeConnectivity = mesh_.cells().node_connectivity();
+    const auto nodeXyView            = array::make_view<double, 2>(mesh_.nodes().xy());
+    const auto& nodeConnectivity     = mesh_.cells().node_connectivity();
 
     // Get four nearest cell-centres to xy.
-    const auto& tree = trees_[size_t(t)];
-    if (tree.size() == 0) {
-        return Cell{{}, Intersect{}.fail()};
+    if (tree_.size() == 0) {
+        return Cell{-1, {}, Intersect{}.fail()};
     }
 
-    const auto valueList = tree.closestPoints(alphabeta, std::min(listSize, tree.size()));
+    const auto valueList = tree_.closestPoints(lonlat, std::min(listSize, tree_.size()));
 
     // Get view to cell flags.
     const auto cellFlagsView = array::make_view<int, 1>(mesh_.cells().flags());
+    const auto cellTijView   = array::make_view<idx_t, 2>(mesh_.cells().field("tij"));
 
     for (const auto& value : valueList) {
         const auto i = value.payload();
 
         // Ignore invalid cells.
-        if (!(cellFlagsView(i) & Topology::INVALID)) {
-            const auto row = nodeConnectivity.row(i);
+        if ((cellFlagsView(i) & Topology::INVALID)) {
+            break;
+        }
 
-            if (row.size() == 4) {
-                auto quadAlphabeta = std::array<Vector2D, 4>{};
-                for (size_t k = 0; k < 4; ++k) {
-                    const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
-                    quadAlphabeta[k]  = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
-                }
+        const auto t   = cellTijView(i);
+        const auto row = nodeConnectivity.row(i);
 
-                const auto quad =
-                    element::Quad2D(quadAlphabeta[0], quadAlphabeta[1], quadAlphabeta[2], quadAlphabeta[3]);
+        if (row.size() == 4) {
+            auto quadAlphabeta = std::array<Vector2D, 4>{};
+            for (size_t k = 0; k < 4; ++k) {
+                const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
+                quadAlphabeta[k]  = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
+            }
+
+            const auto quad =
+                element::Quad2D(quadAlphabeta[0], quadAlphabeta[1], quadAlphabeta[2], quadAlphabeta[3]);
 
-                const auto isect = quad.localRemap(alphabeta, edgeEpsilon, epsilon);
+            const auto alphabeta = projection_->lonlat2alphabeta(lonlat, t);
+            const auto isect = quad.localRemap(alphabeta, edgeEpsilon, epsilon);
 
-                if (isect) {
-                    return Cell{{row(0), row(1), row(2), row(3)}, isect};
-                }
+            if (isect) {
+                return Cell{i, {row(0), row(1), row(2), row(3)}, isect};
+            }
+        }
+        else {
+            // Cell is triangle.
+            auto triagAlphabeta = std::array<Vector2D, 3>{};
+            for (size_t k = 0; k < 3; ++k) {
+                const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
+                triagAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
             }
-            else {
-                // Cell is triangle.
-                auto triagAlphabeta = std::array<Vector2D, 3>{};
-                for (size_t k = 0; k < 3; ++k) {
-                    const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
-                    triagAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
-                }
 
-                const auto triag = element::Triag2D(triagAlphabeta[0], triagAlphabeta[1], triagAlphabeta[2]);
+            const auto triag = element::Triag2D(triagAlphabeta[0], triagAlphabeta[1], triagAlphabeta[2]);
 
-                const auto isect = triag.intersects(alphabeta, edgeEpsilon, epsilon);
+            const auto alphabeta = projection_->lonlat2alphabeta(lonlat, t);
+            const auto isect = triag.intersects(alphabeta, edgeEpsilon, epsilon);
 
-                if (isect) {
-                    return Cell{{row(0), row(1), row(2)}, isect};
-                }
+            if (isect) {
+                return Cell{i, {row(0), row(1), row(2)}, isect};
             }
         }
     }
 
     // Couldn't find a cell.
-    return Cell{{}, Intersect{}.fail()};
+    return Cell{-1, {}, Intersect{}.fail()};
 }
 
-CellFinder::Cell CellFinder::getCell(const PointLonLat& lonlat, size_t listSize, double edgeEpsilon,
+CellFinder::Cell CellFinder::getCell(const PointXY& xy, size_t listSize, double edgeEpsilon,
                                      double epsilon) const {
-    auto xy = PointXY(lonlat.data());
-    projection_->lonlat2xy(xy.data());
-    return getCell(xy, listSize, edgeEpsilon, epsilon);
+    const auto lonlat = projection_->lonlat(xy);
+    return getCell(lonlat, listSize, edgeEpsilon, epsilon);
 }
 
 
diff --git a/src/atlas/interpolation/method/cubedsphere/CellFinder.h b/src/atlas/interpolation/method/cubedsphere/CellFinder.h
index 76932ba99..74424765c 100644
--- a/src/atlas/interpolation/method/cubedsphere/CellFinder.h
+++ b/src/atlas/interpolation/method/cubedsphere/CellFinder.h
@@ -12,6 +12,7 @@
 #include "atlas/interpolation/method/Intersect.h"
 #include "atlas/mesh/Mesh.h"
 #include "atlas/util/Config.h"
+#include "atlas/util/Geometry.h"
 #include "atlas/util/KDTree.h"
 #include "atlas/util/Point.h"
 
@@ -36,6 +37,7 @@ using namespace util;
 class CellFinder {
 public:
     struct Cell {
+        idx_t idx;
         std::vector<idx_t> nodes;
         Intersect isect;
     };
@@ -57,7 +59,7 @@ class CellFinder {
 private:
     Mesh mesh_{};
     const projection::detail::CubedSphereProjectionBase* projection_{};
-    std::array<util::IndexKDTree2D, 6> trees_{};
+    util::IndexKDTree tree_{Geometry{}};
 };
 
 }  // namespace cubedsphere
diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc
index ac0d76353..0dd6527d6 100644
--- a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc
+++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc
@@ -53,6 +53,10 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
     auto weights          = std::vector<Triplet>{};
     const auto ghostView  = array::make_view<int, 1>(target_.ghost());
     const auto lonlatView = array::make_view<double, 2>(target_.lonlat());
+    const auto tijView    = array::make_view<idx_t, 2>(ncSource.mesh().cells().field("tij"));
+
+    // Make vector of tile indices for each target point (needed for vector field interpolation).
+    std::vector<idx_t> tileIndex{};
 
     for (idx_t i = 0; i < target_.size(); ++i) {
         if (!ghostView(i)) {
@@ -66,6 +70,7 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
                     std::to_string(i) + ".");
             }
 
+            tileIndex.push_back(tijView(cell.idx, 0));
             const auto& isect = cell.isect;
             const auto& j     = cell.nodes;
 
@@ -95,6 +100,9 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
     // fill sparse matrix and return.
     Matrix A(target_.size(), source_.size(), weights);
     setMatrix(A);
+
+    // Add tile index metadata to target.
+    target_->metadata().set("tile index", tileIndex);
 }
 
 void CubedSphereBilinear::print(std::ostream&) const {
diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h
index 3cf75c87f..166ec8c61 100644
--- a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h
+++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h
@@ -20,7 +20,7 @@ namespace method {
 ///          accross cubed sphere tiles in (alpha, beta) coordinates.
 ///          Adding int "halo" (default 0) to the config controls the amount of
 ///          halo to consider when seraching for interpolation polygons. Adding
-///          int "list_size" (defualt 4) will change the number of cells
+///          int "list_size" (defualt 8) will change the number of cells
 ///          returned by the internal kd-tree search. Increasing both values
 ///          may resolve errors if setup method fails to find cells.
 ///          The automatic halo exchange in the execute method can be disabled
@@ -46,7 +46,7 @@ class CubedSphereBilinear : public Method {
     FunctionSpace target_;
 
     int halo_{0};
-    int listSize_{4};
+    int listSize_{8};
     bool halo_exchange_{true};
 };
 
diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc
index 74a893e8b..0540ced8f 100644
--- a/src/atlas/output/detail/GmshIO.cc
+++ b/src/atlas/output/detail/GmshIO.cc
@@ -221,7 +221,8 @@ std::vector<int> get_levels(int nlev, const Metadata& gmsh_options) {
 std::string field_lev(const Field& field, int jlev) {
     if (field.levels()) {
         char str[6] = {0, 0, 0, 0, 0, 0};
-        std::sprintf(str, "[%03d]", jlev);
+        auto str_len = std::snprintf(str, sizeof(str), "[%03d]", jlev);
+        ATLAS_ASSERT(str_len == 5);
         return std::string(str);
     }
     else {
@@ -388,8 +389,9 @@ void write_field_nodes(const Metadata& gmsh_options, const functionspace::NoFunc
 
 // ----------------------------------------------------------------------------
 
-void print_field_lev(char field_lev[], int jlev) {
-    std::sprintf(field_lev, "[%03d]", jlev);
+void print_field_lev(char field_lev[], size_t size, int jlev) {
+    ATLAS_ASSERT(size > 5);
+    auto str_len = std::snprintf(field_lev, size, "[%03d]", jlev);
 }
 
 /* unused
@@ -433,7 +435,7 @@ void write_field_nodes(const Metadata& gmsh_options, const functionspace::Struct
         char field_lev[6] = {0, 0, 0, 0, 0, 0};
 
         if (field.levels()) {
-            print_field_lev(field_lev, jlev);
+            print_field_lev(field_lev, sizeof(field_lev), jlev);
         }
 
         out << "$NodeData\n";
diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc
index eddcd9f88..0eb17ddc3 100644
--- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc
+++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc
@@ -22,6 +22,12 @@
 #include "atlas/util/Constants.h"
 #include "atlas/util/CoordinateEnums.h"
 
+
+
+namespace atlas {
+namespace projection {
+namespace detail {
+
 namespace {
 
 static constexpr bool debug     = false;  // constexpr so compiler can optimize `if ( debug ) { ... }` out
@@ -54,14 +60,24 @@ bool greaterEqual(double a, double b) {
     return a > b || equal(a, b);
 }
 
+// Sine/cosine helper struct.
+struct SineCosine {
+    double sinLambda{};
+    double cosLambda{};
+    double sinPhi{};
+    double cosPhi{};
+    explicit SineCosine(const PointLonLat& lonlat) {
+        const double lambda = lonlat.lon() * deg2rad;
+        const double phi    = lonlat.lat() * deg2rad;
+        sinLambda           = std::sin(lambda);
+        cosLambda           = std::cos(lambda);
+        sinPhi              = std::sin(phi);
+        cosPhi              = std::sqrt(1. - sinPhi * sinPhi); // Safe, as phi is between +/- 90 degrees
+    }
+};
 
 }  // namespace
 
-namespace atlas {
-namespace projection {
-namespace detail {
-
-
 // -------------------------------------------------------------------------------------------------
 
 CubedSphereEquiAnglProjection::CubedSphereEquiAnglProjection(const eckit::Parametrisation& params):
@@ -85,10 +101,10 @@ void CubedSphereEquiAnglProjection::xy2alphabeta(double crd[], idx_t t) const {
     }
 
     // Get alphaBeta Jacobian.
-    const auto alphabetaJacobian = getCubedSphereTiles().tileJacobian(static_cast<size_t>(t)).inverse();
+    const auto abJacobian = getCubedSphereTiles().tileJacobian(static_cast<size_t>(t)).inverse();
 
     // Set (alpha, beta) coord.
-    const Point2 alphabeta = alphabetaJacobian * (Point2(crd) - xyCentre);
+    const Point2 alphabeta = abJacobian * (Point2(crd) - xyCentre);
     crd[0]                 = alphabeta[0];
     crd[1]                 = alphabeta[1];
 
@@ -159,17 +175,73 @@ void CubedSphereEquiAnglProjection::alphabeta2xy(double crd[], idx_t t) const {
 
 // -------------------------------------------------------------------------------------------------
 
+void CubedSphereEquiAnglProjection::lonlat2alphabeta(double crd[], idx_t t) const {
+
+    const auto angles = SineCosine(PointLonLat(crd));
+
+    // Convert lonlat to xyz on unit sphere.
+    auto xyz = PointXYZ{angles.cosLambda * angles.cosPhi, angles.sinLambda * angles.cosPhi, -angles.sinPhi};
+
+    const auto& tiles = getCubedSphereTiles();
+    tiles.unrotate(t, xyz.data());
+
+    // project into xy and rotate into alphabeta coordinates.
+    auto alphaBeta = tiles.tileJacobian(t).inverse() *
+                     Point2{std::atan2(xyz.y(), xyz.x()), -std::atan2(xyz.z(), xyz.x())};
+
+    alphaBeta = alphaBeta * rad2deg;
+
+    crd[0] = alphaBeta[0];
+    crd[1] = alphaBeta[1];
+}
+
+
+// -------------------------------------------------------------------------------------------------
+
+void CubedSphereEquiAnglProjection::alphabeta2lonlat(double crd[], idx_t t) const {
+
+    // Rotate alphabeta into xy direction.
+    const auto& tiles = getCubedSphereTiles();
+    auto alphaBeta = tiles.tileJacobian(t) * Point2{crd};
+    alphaBeta = alphaBeta * deg2rad;
+
+    // project alphabeta into xyz space.
+    auto xyz = PointXYZ{-1., -std::tan(alphaBeta[0]), std::tan(alphaBeta[1])};
+
+    tiles.rotate(t, xyz.data());
+
+    // Honestly, I figured out some of the -1 factors by trial and error.
+    const auto r = PointXYZ::norm(xyz);
+
+    auto lonlat = PointLonLat{};
+    if (equal(xyz.x(), 0.) && equal(xyz.y(), 0.)) {
+        lonlat.lon() = 0.;
+    }
+    else {
+        lonlat.lon() = std::atan2(-xyz.y(), -xyz.x());
+    }
+    lonlat.lat() = std::asin(xyz.z() / r);
+
+    lonlat = lonlat * rad2deg;
+    lonlat.normalise();
+
+    crd[0] = lonlat.lon();
+    crd[1] = lonlat.lat();
+}
+
+// -------------------------------------------------------------------------------------------------
+
 Jacobian CubedSphereEquiAnglProjection::jacobian(const PointLonLat& lonlat, idx_t t) const {
     // Note: angular units cancel, so we leave all values in radians.
 
+    const auto angles = SineCosine(lonlat);
+
     // Convert lonlat to xyz on unit sphere.
-    const double lambda = lonlat.lon() * deg2rad;
-    const double phi    = lonlat.lat() * deg2rad;
-    auto xyz            = PointXYZ{std::cos(lambda) * std::cos(phi), std::sin(lambda) * std::cos(phi), -std::sin(phi)};
+    auto xyz = PointXYZ{angles.cosLambda * angles.cosPhi, angles.sinLambda * angles.cosPhi, -angles.sinPhi};
 
     // Get derivatives of xyz with respect to lambda and phi.
-    auto dxyz_by_dlambda = PointXYZ{-std::sin(lambda) * std::cos(phi), std::cos(lambda) * std::cos(phi), 0.};
-    auto dxyz_by_dphi = PointXYZ{-std::cos(lambda) * std::sin(phi), -std::sin(lambda) * std::sin(phi), -std::cos(phi)};
+    auto dxyz_by_dlambda = PointXYZ{-angles.sinLambda * angles.cosPhi, angles.cosLambda * angles.cosPhi, 0.};
+    auto dxyz_by_dphi = PointXYZ{-angles.cosLambda * angles.sinPhi, -angles.sinLambda * angles.sinPhi, -angles.cosPhi};
 
     // Rotate vectors.
     const auto& tiles = getCubedSphereTiles();
@@ -218,14 +290,14 @@ void CubedSphereEquiAnglProjection::lonlat2xy(double crd[]) const {
     // left coordinate system
 
     if (debug) {
-        Log::info() << "equiangular lonlat2xy xyz ab : " << xyz[XX] << " " << xyz[YY] << " " << xyz[ZZ] << " "
+        Log::debug() << "equiangular lonlat2xy xyz ab : " << xyz[XX] << " " << xyz[YY] << " " << xyz[ZZ] << " "
                     << ab[LON] << " " << ab[LAT] << std::endl;
     }
 
     CubedSphereProjectionBase::alphabetat2xy(t, ab, crd);
 
     if (debug) {
-        Log::info() << "equiangular lonlat2xy end : xy = " << crd[LON] << " " << crd[LAT] << std::endl;
+        Log::debug() << "equiangular lonlat2xy end : xy = " << crd[LON] << " " << crd[LAT] << std::endl;
     }
 }
 
@@ -234,7 +306,7 @@ void CubedSphereEquiAnglProjection::lonlat2xy(double crd[]) const {
 //
 void CubedSphereEquiAnglProjection::xy2lonlat(double crd[]) const {
     if (debug) {
-        Log::info() << "xy2lonlat start xy = " << crd[LON] << " " << crd[LAT] << std::endl;
+        Log::debug() << "xy2lonlat start xy = " << crd[LON] << " " << crd[LAT] << std::endl;
     }
 
     static const double rsq3 = 1.0 / std::sqrt(3.0);
@@ -257,7 +329,7 @@ void CubedSphereEquiAnglProjection::xy2lonlat(double crd[]) const {
     CubedSphereProjectionBase::xy2lonlat_post(xyz, t, crd);
 
     if (debug) {
-        Log::info() << "end of equiangular xy2lonlat lonlat = " << crd[LON] << " " << crd[LAT] << std::endl;
+        Log::debug() << "end of equiangular xy2lonlat lonlat = " << crd[LON] << " " << crd[LAT] << std::endl;
     }
 }
 
diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h
index c297919df..7c2631f74 100644
--- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h
+++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h
@@ -34,6 +34,12 @@ class CubedSphereEquiAnglProjection final : public CubedSphereProjectionBase {
     /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t.
     void alphabeta2xy(double crd[], idx_t t) const override;
 
+    /// @brief Convert (lon, lat) coordinate to (alpha, beta) on tile t.
+    void lonlat2alphabeta(double crd[], idx_t t) const override;
+
+    /// @brief Convert (alpha, beta) coordinate to (lon, lat) on tile t.
+    void alphabeta2lonlat(double crd[], idx_t t) const override;
+
     /// @brief Jacobian of (x, y) with respect to (lon, lat) on tile t
     Jacobian jacobian(const PointLonLat& lonlat, idx_t t) const override;
 
diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc
index 0ff91b55c..95c3fdc4d 100644
--- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc
+++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc
@@ -47,6 +47,18 @@ void CubedSphereEquiDistProjection::alphabeta2xy(double crd[], idx_t t) const {
 
 // -------------------------------------------------------------------------------------------------
 
+void CubedSphereEquiDistProjection::lonlat2alphabeta(double crd[], idx_t t) const {
+    throw_NotImplemented("lonlat2alphabeta not implemented for CubedSphereEquiDistProjection", Here());
+}
+
+// -------------------------------------------------------------------------------------------------
+
+void CubedSphereEquiDistProjection::alphabeta2lonlat(double crd[], idx_t t) const {
+    throw_NotImplemented("alphabeta2lonlat not implemented for CubedSphereEquiDistProjection", Here());
+}
+
+// -------------------------------------------------------------------------------------------------
+
 Jacobian CubedSphereEquiDistProjection::jacobian(const PointLonLat& lonlat, idx_t t) const {
     throw_NotImplemented("jacobian not implemented for CubedSphereEquiDistProjection", Here());
 }
diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h
index 6d7ea0b7f..60480b909 100644
--- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h
+++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h
@@ -33,6 +33,12 @@ class CubedSphereEquiDistProjection final : public CubedSphereProjectionBase {
     /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t.
     void alphabeta2xy(double crd[], idx_t t) const override;
 
+    /// @brief Convert (lon, lat) coordinate to (alpha, beta) on tile t.
+    void lonlat2alphabeta(double crd[], idx_t t) const override;
+
+    /// @brief Convert (alpha, beta) coordinate to (lon, lat) on tile t.
+    void alphabeta2lonlat(double crd[], idx_t t) const override;
+
     /// @brief Jacobian of (alpha, beta) with respect to (lon, lat) on tile t
     Jacobian alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const override;
 
diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.cc b/src/atlas/projection/detail/CubedSphereProjectionBase.cc
index a8295be4b..be98ff3f4 100644
--- a/src/atlas/projection/detail/CubedSphereProjectionBase.cc
+++ b/src/atlas/projection/detail/CubedSphereProjectionBase.cc
@@ -248,6 +248,22 @@ Point2 CubedSphereProjectionBase::alphabeta2xy(const Point2& alphabeta, idx_t t)
 
 // -------------------------------------------------------------------------------------------------
 
+Point2 CubedSphereProjectionBase::lonlat2alphabeta(const Point2& lonlat, idx_t t) const {
+    auto alphabeta = Point2(lonlat);
+    lonlat2alphabeta(alphabeta.data(), t);
+    return alphabeta;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+Point2 CubedSphereProjectionBase::alphabeta2lonlat(const Point2& alphabeta, idx_t t) const {
+    auto lonlat = Point2(alphabeta);
+    alphabeta2lonlat(lonlat.data(), t);
+    return lonlat;
+}
+
+// -------------------------------------------------------------------------------------------------
+
 }  // namespace detail
 }  // namespace projection
 }  // namespace atlas
diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.h b/src/atlas/projection/detail/CubedSphereProjectionBase.h
index b71154853..0989dc6a2 100644
--- a/src/atlas/projection/detail/CubedSphereProjectionBase.h
+++ b/src/atlas/projection/detail/CubedSphereProjectionBase.h
@@ -58,6 +58,23 @@ class CubedSphereProjectionBase : public ProjectionImpl {
     virtual void alphabeta2xy(double crd[], idx_t) const = 0;
     ///@}
 
+    /// @brief   Convert (lon, lat) coordinate to (alpha, beta) on tile t
+    ///
+    /// @details Converts the lon lat coordinates to the angular coordinates
+    ///          described of tile t, described by by Ronchi et al. (1996,
+    ///          Journal of Computational Physics, 124, 93).
+    /// @{
+    Point2 lonlat2alphabeta(const Point2& lonlat, idx_t t) const;
+    virtual void lonlat2alphabeta(double crd[], idx_t) const = 0;
+    /// @}
+
+    /// @brief   Convert (lon, lat) coordinate to (alpha, beta) on tile t
+    ///
+    /// @details Performs the inverse of lonlat2alphabeta.
+    ///
+    Point2 alphabeta2lonlat(const Point2& alphabeta, idx_t t) const;
+    virtual void alphabeta2lonlat(double crd[], idx_t) const = 0;
+
     /// @brief   Jacobian of (x, y) with respect to (lon, lat) on tile t
     ///
     /// @details Returns the Jacobian
diff --git a/src/tests/interpolation/test_interpolation_cubedsphere.cc b/src/tests/interpolation/test_interpolation_cubedsphere.cc
index 04c3b1cb5..8d77fd461 100644
--- a/src/tests/interpolation/test_interpolation_cubedsphere.cc
+++ b/src/tests/interpolation/test_interpolation_cubedsphere.cc
@@ -262,10 +262,12 @@ CASE("cubedsphere_wind_interpolation") {
         const auto vBeta  = array::make_view<double, 1>(targetFieldSet["v_beta"]);
         auto error0       = array::make_view<double, 1>(targetFieldSet["error_field_0"]);
         auto error1       = array::make_view<double, 1>(targetFieldSet["error_field_1"]);
+        const auto& tVec  = interp.target()->metadata().getIntVector("tile index");
+
         for (idx_t idx = 0; idx < targetFunctionspace.size(); ++idx) {
             if (!ghost(idx)) {
                 const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT));
-                const idx_t t = proj.getCubedSphereTiles().indexFromLonLat(ll);
+                const idx_t t = tVec[idx];
 
                 // Get inverse wind transform jacobian.
                 auto jac = windTransform(ll, t).inverse();
@@ -322,12 +324,13 @@ CASE("cubedsphere_wind_interpolation") {
         const auto vAdj   = array::make_view<double, 1>(targetFieldSet["v_adjoint"]);
         auto vAlphaAdj    = array::make_view<double, 1>(targetFieldSet["v_alpha_adjoint"]);
         auto vBetaAdj     = array::make_view<double, 1>(targetFieldSet["v_beta_adjoint"]);
+        const auto& tVec  = interp.target()->metadata().getIntVector("tile index");
         vAlphaAdj.assign(0.);
         vBetaAdj.assign(0.);
         for (idx_t idx = 0; idx < targetFunctionspace.size(); ++idx) {
             if (!ghost(idx)) {
                 const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT));
-                const idx_t t = proj.getCubedSphereTiles().indexFromLonLat(ll);
+                const idx_t t = tVec[idx];
 
                 // Get adjoint of inverse wind transform jacobian.
                 auto jac = windTransform(ll, t).inverse().transpose();
diff --git a/src/tests/projection/test_cubedsphere_projection.cc b/src/tests/projection/test_cubedsphere_projection.cc
index 24b9b405b..59d836f11 100644
--- a/src/tests/projection/test_cubedsphere_projection.cc
+++ b/src/tests/projection/test_cubedsphere_projection.cc
@@ -20,7 +20,7 @@ namespace atlas {
 namespace test {
 
 // Small number relative to 360.
-constexpr double epsilon = std::numeric_limits<double>::epsilon() * 360.;
+constexpr double epsilon = 2. * std::numeric_limits<double>::epsilon() * 360.;
 
 void testProjection(const std::string& gridType, const std::string& meshType, const std::string& outputID) {
     // Create grid.
@@ -52,6 +52,7 @@ void testProjection(const std::string& gridType, const std::string& meshType, co
     // "Hack" mesh xy coordinates.
     auto xyView        = array::make_view<double, 2>(mesh.nodes().xy());
     const auto tijView = array::make_view<idx_t, 2>(mesh.nodes().field("tij"));
+    const auto lonlatView = array::make_view<double, 2>(mesh.nodes().lonlat());
     for (idx_t i = 0; i < mesh.nodes().size(); ++i) {
         const auto t = tijView(i, 0);
 
@@ -73,6 +74,15 @@ void testProjection(const std::string& gridType, const std::string& meshType, co
         // overwrite mesh xy field.
         xyView(i, XX) = alphabeta[0];
         xyView(i, YY) = alphabeta[1];
+
+        // Check lonlat2alphabeta methods.
+        const auto lonlat       = csProjection.alphabeta2lonlat(alphabeta, t);
+        const auto newAlphabeta = csProjection.lonlat2alphabeta(lonlat, t);
+
+        EXPECT_APPROX_EQ(alphabeta[0], newAlphabeta[0], epsilon);
+        EXPECT_APPROX_EQ(alphabeta[1], newAlphabeta[1], epsilon);
+        EXPECT_APPROX_EQ(lonlat[0], lonlatView(i, 0), epsilon);
+        EXPECT_APPROX_EQ(lonlat[1], lonlatView(i, 1), epsilon);
     }
 
     // Output mesh updated mesh.