src/horizontal/CMakeLists.txt

@@ -11,7 +11,7 @@
 
 add_library(
   iconmath-horizontal
-  lib_divrot.cpp
+  mo_lib_divrot.cpp
   mo_lib_divrot.F90
   mo_lib_laplace.F90
   mo_lib_gradients.F90)

src/horizontal/lib_divrot.cpp → src/horizontal/mo_lib_divrot.cpp

@@ -12,7 +12,7 @@
 #include <iostream>
 #include <vector>
 
-#include <horizontal/lib_divrot.hpp>
+#include <horizontal/mo_lib_divrot.hpp>
 #include <support/mo_lib_loopindices.hpp>
 
 template <typename T>

test/c/CMakeLists.txt

@@ -27,7 +27,9 @@ endif()
 
 set(SOURCES
   main.cpp
-  test_horizontal_divrot.cpp
+  test_horizontal_div.cpp
+  test_horizontal_recon.cpp
+  test_horizontal_rot.cpp
   test_tdma_solver.cpp
   test_interpolation_vector.cpp
   test_intp_rbf.cpp
@@ -36,6 +38,8 @@ set(SOURCES
 # Create the test executable from your test files, including main.cpp.
 add_executable(iconmath_test_c ${SOURCES})
 
+target_include_directories(iconmath_test_c PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+
 # Link the test executable with GoogleTest and Kokkos.
 target_link_libraries(iconmath_test_c
   PUBLIC

test/c/dim_helper.hpp

+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2025, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#pragma once
+
+// Template function for computing array size.
+// For example, we get the array size of a 4-dimensional array A(2, 3, 4, 5) by
+//    dim_combine(2, 3, 4, 5).
+// Which will automatically instantiate
+//    dim_combine<int, int, int, int>(2, 3, 4, 5).
+// The function then call dim_combine recursively
+//    dim_combine<int, int, int, int>(2, 3, 4, 5) {
+//      return static_cast<size_t>(2) * dim_combine<int, int, int>(3, 4, 5);
+//    }
+//    dim_combine<int, int, int>(3, 4, 5) {
+//      return static_cast<size_t>(3) * dim_combine<int, int>(4, 5);
+//    }
+//    dim_combine<int, int>(4, 5) {
+//      return static_cast<size_t>(4) * dim_combine<int>(5);
+//    }
+// Where the last dim_combine is specialized as
+//    dim_combine<int>(5) {
+//      return static_cast<size_t>(5);
+//    }
+// Which gives
+//    dim_combine<int, int, int, int>(2, 3, 4, 5) =
+//      static_cast<size_t>(2) * static_cast<size_t>(3) *
+//      static_cast<size_t>(4) * static_cast<size_t>(5)
+/// Template helpers for combining multiple dimension array sizes.
+/// The base function of dimension combine. Should not be used.
+template <typename... Ts> size_t dim_combine(Ts... dims) { return 0; }
+/// Template specialization of only one dimension, returns the dimension itself.
+template <typename T> size_t dim_combine(T dim) {
+  return static_cast<size_t>(dim);
+}
+/// Template specialization of picking out the first dimension. The combined
+/// dimension is the first dimension times the combined dimension of the rest.
+template <typename T, typename... Ts> size_t dim_combine(T dim, Ts... dims) {
+  return static_cast<size_t>(dim) * dim_combine(dims...);
+}
+
+// Template function for LayoutLeft ID access in compile time.
+// For example, a multi-dimensional array A of dimensions <2, 3, 4, 5> gets its
+// corresponding vector id (LayoutLeft) by
+//    at<2, 3, 4, 5>(id1, id2, id3, id4).
+// The at_impl then adds the id from beginning to the end and pass the id prefix
+// to the next recursive at_impl function. In this example,
+//    at<2, 3, 4, 5>(id1, id2, id3, id4) {
+//      return id1 + at_impl<3, 4, 5>(2, id2, id3, id4);
+//    }
+//    at_impl<3, 4, 5>(2, id2, id3, id4) {
+//      return id2 * 2 + at_impl<4, 5>(2 * 3, id3, id4);
+//    }
+//    at_impl<4, 5>(2 * 3, id3, id4) {
+//      return id3 * 2 * 3 + at_impl<5>(2 * 3 * 4, id4);
+//    }
+//    at_impl<5>(2 * 3 * 4, id4) {
+//      return id4 * 2 * 3 * 4;
+//    }
+// Which gives
+//    at<2, 3, 4, 5>(id1, id2, id3, id4) = id1         + id2 * 2         +
+//                                         id3 * 2 * 3 + id4 * 2 * 3 * 4
+/// Helper type converting integer numbers to int
+template <class T, auto> using always_t = T;
+/// Base function of at_impl. Should not be used.
+template <int... Dims> int at_impl(always_t<int, Dims>... ids) { return 0; }
+/// Template specialization of the last ID
+template <int LastDim> int at_impl(int prefix, int id) { return id * prefix; }
+/// Template specialization of at_impl, accumulate the return value using the
+/// first id and pass the prefix to the next recursive at_impl function.
+template <int FirstDim, int... Dims>
+int at_impl(int prefix, int id, always_t<int, Dims>... ids) {
+  return id * prefix + at_impl<Dims...>(prefix * FirstDim, ids...);
+}
+/// at<dim1, dim2, ...>(id1, id2, ...) gets its memory index in vector assuming
+/// LayoutLeft. Use this function instead of at_impl.
+template <int FirstDim, int... Dims>
+int at(int id, always_t<int, Dims>... ids) {
+  return id + at_impl<Dims...>(FirstDim, ids...);
+}

test/c/test_horizontal_rot.cpp

+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2025, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include <iostream>
+#include <random>
+#include <vector>
+
+#include <Kokkos_Core.hpp>
+#include <gtest/gtest.h>
+#include <dim_helper.hpp>
+#include <horizontal/mo_lib_divrot.hpp>
+#include <support/mo_lib_loopindices.hpp>
+
+/// Test class for the horizontal rotation tests. Templated for the ValueType.
+template <typename ValueType>
+class HorizontalRotVertexTest : public ::testing::Test {
+protected:
+  static constexpr int nproma = 3;  // inner loop length
+  static constexpr int nlev = 2;    // number of vertical levels
+  static constexpr int nblks_e = 1; // number of edge blocks
+  static constexpr int nblks_v = 1; // number of vertex blocks
+  static constexpr int dim4d = 2;   // 4th dimension size
+
+  int i_startblk = 0;
+  int i_endblk = nblks_v; // Test blocks [0 .. nblks_v-1]
+  int i_startidx_in = 0;
+  int i_endidx_in = nproma; // Full range: 0 .. nproma-1
+  std::vector<int> slev;
+  std::vector<int> elev;
+  bool lacc = false;      // Not using ACC-specific behavior.
+  bool acc_async = false; // Not using ACC-specific behavior.
+
+  std::vector<ValueType> vec_e;
+  std::vector<int> vert_edge_idx;
+  std::vector<int> vert_edge_blk;
+  std::vector<ValueType> geofac_rot;
+  std::vector<ValueType> rot_vec;
+  std::vector<ValueType> f4din;
+  std::vector<ValueType> f4dout;
+
+  HorizontalRotVertexTest() {
+    slev.resize(dim4d, 0);
+    elev.resize(dim4d, nlev); // Full vertical range (0 .. nlev-1)
+
+    vec_e.resize(dim_combine(nproma, nlev, nblks_e));
+    vert_edge_idx.resize(dim_combine(nproma, nblks_v, 6));
+    vert_edge_blk.resize(dim_combine(nproma, nblks_v, 6));
+    geofac_rot.resize(dim_combine(nproma, 6, nblks_v));
+    rot_vec.resize(dim_combine(nproma, nlev, nblks_v));
+    f4din.resize(dim_combine(nproma, nlev, nblks_e, dim4d));
+    f4dout.resize(dim_combine(nproma, nlev, nblks_v, dim4d));
+  }
+};
+
+/// ValueTypes which the divrot tests should run with
+typedef ::testing::Types<float, double> ValueTypes;
+
+TYPED_TEST_SUITE(HorizontalRotVertexTest, ValueTypes);
+
+TYPED_TEST(HorizontalRotVertexTest, TestRotVertexAtmosSpecific) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_e = this->nblks_e;
+  constexpr int nblks_v = this->nblks_v;
+
+  const auto &vec_e_at = at<nproma, nlev, nblks_e>;
+  const auto &vert_edge_at = at<nproma, nblks_v, 6>;
+  const auto &geofac_rot_at = at<nproma, 6, nblks_v>;
+  const auto &rot_vec_at = at<nproma, nlev, nblks_v>;
+
+  // Initialization with specific values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->vec_e[vec_e_at(i, k, 0)] = (i + 1) * (k + 1); // Simple pattern
+    }
+
+    // Set edge indices to point to specific edges
+    for (int j = 0; j < 6; ++j) {
+      this->vert_edge_idx[vert_edge_at(i, 0, j)] = (i + j) % nproma;
+      // All edges are in the same block for this test
+      this->vert_edge_blk[vert_edge_at(i, 0, j)] = 0;
+    }
+
+    // Geometric factors for rotation
+    this->geofac_rot[geofac_rot_at(i, 0, 0)] = 0.3;
+    this->geofac_rot[geofac_rot_at(i, 1, 0)] = 0.2;
+    this->geofac_rot[geofac_rot_at(i, 2, 0)] = 0.1;
+    this->geofac_rot[geofac_rot_at(i, 3, 0)] = 0.2;
+    this->geofac_rot[geofac_rot_at(i, 4, 0)] = 0.1;
+    this->geofac_rot[geofac_rot_at(i, 5, 0)] = 0.1;
+
+    // Initialize rot_vec to zero
+    for (int k = 0; k < nlev; ++k) {
+      this->rot_vec[rot_vec_at(i, k, 0)] = 0.0;
+    }
+  }
+
+  // Call the rot_vertex_atmos function
+  rot_vertex_atmos<TypeParam>(
+      this->vec_e.data(), this->vert_edge_idx.data(),
+      this->vert_edge_blk.data(), this->geofac_rot.data(), this->rot_vec.data(),
+      this->i_startblk, this->i_endblk, this->i_startidx_in, this->i_endidx_in,
+      this->slev[0], this->elev[0], this->nproma, this->lacc, this->nlev,
+      this->nblks_e, this->nblks_v);
+
+  // Expected values based on the initialization pattern
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(0, 0, 0)], 1.7, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(0, 1, 0)], 3.4, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(1, 0, 0)], 2.1, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(1, 1, 0)], 4.2, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(2, 0, 0)], 2.2, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(2, 1, 0)], 4.4, 1e-6);
+}
+
+TYPED_TEST(HorizontalRotVertexTest, TestRotVertexAtmosRandom) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_e = this->nblks_e;
+  constexpr int nblks_v = this->nblks_v;
+
+  const auto &vec_e_at = at<nproma, nlev, nblks_e>;
+  const auto &vert_edge_at = at<nproma, nblks_v, 6>;
+  const auto &geofac_rot_at = at<nproma, 6, nblks_v>;
+  const auto &rot_vec_at = at<nproma, nlev, nblks_v>;
+
+  // Set up random number generators
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<int> int_distrib(0, nproma - 1);
+  std::uniform_real_distribution<TypeParam> real_distrib(-10.0, 10.0);
+
+  // Initialization with random values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->vec_e[vec_e_at(i, k, 0)] = real_distrib(gen);
+    }
+
+    // Set random edge indices
+    for (int j = 0; j < 6; ++j) {
+      this->vert_edge_idx[vert_edge_at(i, 0, j)] = int_distrib(gen);
+      this->vert_edge_blk[vert_edge_at(i, 0, j)] =
+          0; // Keep in same block for simplicity
+    }
+
+    // Random geometric factors
+    for (int j = 0; j < 6; ++j) {
+      this->geofac_rot[geofac_rot_at(i, j, 0)] = real_distrib(gen);
+    }
+
+    // Initialize rot_vec to random values
+    for (int k = 0; k < nlev; ++k) {
+      this->rot_vec[rot_vec_at(i, k, 0)] = real_distrib(gen);
+    }
+  }
+
+  // Call the rot_vertex_atmos function
+  rot_vertex_atmos<TypeParam>(
+      this->vec_e.data(), this->vert_edge_idx.data(),
+      this->vert_edge_blk.data(), this->geofac_rot.data(), this->rot_vec.data(),
+      this->i_startblk, this->i_endblk, this->i_startidx_in, this->i_endidx_in,
+      this->slev[0], this->elev[0], this->nproma, this->lacc, this->nlev,
+      this->nblks_e, this->nblks_v);
+
+  // Calculate reference values separately and verify results
+  std::vector<TypeParam> ref_rot_vec(nproma * nlev * nblks_v, 0.0);
+
+  for (int jb = this->i_startblk; jb < this->i_endblk; ++jb) {
+    int i_startidx, i_endidx;
+    get_indices_v_lib(this->i_startidx_in, this->i_endidx_in, nproma, jb,
+                      this->i_startblk, this->i_endblk, i_startidx, i_endidx);
+
+    for (int jk = this->slev[0]; jk < this->elev[0]; ++jk) {
+      for (int jv = i_startidx; jv < i_endidx; ++jv) {
+        ref_rot_vec[rot_vec_at(jv, jk, jb)] =
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 0)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 0)])] *
+                this->geofac_rot[geofac_rot_at(jv, 0, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 1)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 1)])] *
+                this->geofac_rot[geofac_rot_at(jv, 1, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 2)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 2)])] *
+                this->geofac_rot[geofac_rot_at(jv, 2, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 3)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 3)])] *
+                this->geofac_rot[geofac_rot_at(jv, 3, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 4)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 4)])] *
+                this->geofac_rot[geofac_rot_at(jv, 4, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 5)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 5)])] *
+                this->geofac_rot[geofac_rot_at(jv, 5, jb)];
+      }
+    }
+  }
+
+  // Verify results
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      EXPECT_NEAR(this->rot_vec[rot_vec_at(i, k, 0)],
+                  ref_rot_vec[rot_vec_at(i, k, 0)], 1e-5)
+          << "Results differ at i=" << i << ", k=" << k;
+    }
+  }
+}
+
+TYPED_TEST_SUITE(HorizontalRotVertexTest, ValueTypes);
+
+TYPED_TEST(HorizontalRotVertexTest, TestRotVertexRISpecific) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_e = this->nblks_e;
+  constexpr int nblks_v = this->nblks_v;
+
+  const auto &vec_e_at = at<nproma, nlev, nblks_e>;
+  const auto &vert_edge_at = at<nproma, nblks_v, 6>;
+  const auto &geofac_rot_at = at<nproma, 6, nblks_v>;
+  const auto &rot_vec_at = at<nproma, nlev, nblks_v>;
+
+  // Initialization with specific values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->vec_e[vec_e_at(i, k, 0)] = (i + 1) * (k + 1); // Simple pattern
+    }
+
+    // Set edge indices to point to specific edges
+    for (int j = 0; j < 6; ++j) {
+      this->vert_edge_idx[vert_edge_at(i, 0, j)] = (i + j) % nproma;
+      // All edges are in the same block for this test
+      this->vert_edge_blk[vert_edge_at(i, 0, j)] = 0;
+    }
+
+    // Geometric factors for rotation
+    this->geofac_rot[geofac_rot_at(i, 0, 0)] = 0.3;
+    this->geofac_rot[geofac_rot_at(i, 1, 0)] = 0.2;
+    this->geofac_rot[geofac_rot_at(i, 2, 0)] = 0.1;
+    this->geofac_rot[geofac_rot_at(i, 3, 0)] = 0.2;
+    this->geofac_rot[geofac_rot_at(i, 4, 0)] = 0.1;
+    this->geofac_rot[geofac_rot_at(i, 5, 0)] = 0.1;
+
+    // Initialize rot_vec to zero
+    for (int k = 0; k < nlev; ++k) {
+      this->rot_vec[rot_vec_at(i, k, 0)] = 0.0;
+    }
+  }
+
+  // Call the rot_vertex_ri function
+  rot_vertex_ri<TypeParam>(
+      this->vec_e.data(), this->vert_edge_idx.data(),
+      this->vert_edge_blk.data(), this->geofac_rot.data(), this->rot_vec.data(),
+      this->i_startblk, this->i_endblk, this->i_startidx_in, this->i_endidx_in,
+      this->slev[0], this->elev[0], this->nproma, this->lacc, this->acc_async,
+      this->nlev, this->nblks_e, this->nblks_v);
+
+  // Expected values based on the initialization pattern
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(0, 0, 0)], 1.7, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(0, 1, 0)], 3.4, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(1, 0, 0)], 2.1, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(1, 1, 0)], 4.2, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(2, 0, 0)], 2.2, 1e-6);
+  EXPECT_NEAR(this->rot_vec[rot_vec_at(2, 1, 0)], 4.4, 1e-6);
+}
+
+TYPED_TEST(HorizontalRotVertexTest, TestRotVertexRIRandom) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_e = this->nblks_e;
+  constexpr int nblks_v = this->nblks_v;
+
+  const auto &vec_e_at = at<nproma, nlev, nblks_e>;
+  const auto &vert_edge_at = at<nproma, nblks_v, 6>;
+  const auto &geofac_rot_at = at<nproma, 6, nblks_v>;
+  const auto &rot_vec_at = at<nproma, nlev, nblks_v>;
+
+  // Set up random number generators
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<int> int_distrib(0, nproma - 1);
+  std::uniform_real_distribution<TypeParam> real_distrib(-10.0, 10.0);
+
+  // Initialization with random values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->vec_e[vec_e_at(i, k, 0)] = real_distrib(gen);
+    }
+
+    // Set random edge indices
+    for (int j = 0; j < 6; ++j) {
+      this->vert_edge_idx[vert_edge_at(i, 0, j)] = int_distrib(gen);
+      this->vert_edge_blk[vert_edge_at(i, 0, j)] =
+          0; // Keep in same block for simplicity
+    }
+
+    // Random geometric factors
+    for (int j = 0; j < 6; ++j) {
+      this->geofac_rot[geofac_rot_at(i, j, 0)] = real_distrib(gen);
+    }
+
+    // Initialize rot_vec to random values
+    for (int k = 0; k < nlev; ++k) {
+      this->rot_vec[rot_vec_at(i, k, 0)] = real_distrib(gen);
+    }
+  }
+
+  // Call the rot_vertex_ri function
+  rot_vertex_ri<TypeParam>(
+      this->vec_e.data(), this->vert_edge_idx.data(),
+      this->vert_edge_blk.data(), this->geofac_rot.data(), this->rot_vec.data(),
+      this->i_startblk, this->i_endblk, this->i_startidx_in, this->i_endidx_in,
+      this->slev[0], this->elev[0], this->nproma, this->lacc, this->acc_async,
+      this->nlev, this->nblks_e, this->nblks_v);
+
+  // Ensure computation is complete for both modes
+  Kokkos::fence();
+
+  // Calculate reference values separately and verify results
+  std::vector<TypeParam> ref_rot_vec(nproma * nlev * nblks_v, 0.0);
+
+  for (int jb = this->i_startblk; jb < this->i_endblk; ++jb) {
+    int i_startidx, i_endidx;
+    get_indices_v_lib(this->i_startidx_in, this->i_endidx_in, nproma, jb,
+                      this->i_startblk, this->i_endblk, i_startidx, i_endidx);
+
+    for (int jk = this->slev[0]; jk < this->elev[0]; ++jk) {
+      for (int jv = i_startidx; jv < i_endidx; ++jv) {
+        ref_rot_vec[rot_vec_at(jv, jk, jb)] =
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 0)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 0)])] *
+                this->geofac_rot[geofac_rot_at(jv, 0, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 1)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 1)])] *
+                this->geofac_rot[geofac_rot_at(jv, 1, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 2)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 2)])] *
+                this->geofac_rot[geofac_rot_at(jv, 2, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 3)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 3)])] *
+                this->geofac_rot[geofac_rot_at(jv, 3, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 4)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 4)])] *
+                this->geofac_rot[geofac_rot_at(jv, 4, jb)] +
+            this->vec_e[vec_e_at(
+                this->vert_edge_idx[vert_edge_at(jv, jb, 5)], jk,
+                this->vert_edge_blk[vert_edge_at(jv, jb, 5)])] *
+                this->geofac_rot[geofac_rot_at(jv, 5, jb)];
+      }
+    }
+  }
+
+  // Verify results
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      EXPECT_NEAR(this->rot_vec[rot_vec_at(i, k, 0)],
+                  ref_rot_vec[rot_vec_at(i, k, 0)], 1e-5)
+          << "Results differ at i=" << i << ", k=" << k << ")";
+    }
+  }
+}
+