diff --git a/cpputil/usid/.clang-format b/cpputil/usid/.clang-format
new file mode 100644
index 0000000..09cc831
--- /dev/null
+++ b/cpputil/usid/.clang-format
@@ -0,0 +1,60 @@
+---
+Language:        Cpp
+# BasedOnStyle:  LLVM
+AccessModifierOffset: -2
+ConstructorInitializerIndentWidth: 4
+#AlignEscapedNewlinesLeft: false
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AlwaysBreakTemplateDeclarations: true
+AlwaysBreakBeforeMultilineStrings: false
+BreakBeforeBinaryOperators: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BinPackParameters: false
+BinPackArguments: false
+ColumnLimit:     120
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+DerivePointerAlignment: false
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: false
+IndentWrappedFunctionNames: false
+IndentFunctionDeclarationAfterType: false
+MaxEmptyLinesToKeep: 1
+KeepEmptyLinesAtTheStartOfBlocks: true
+NamespaceIndentation: All
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Right
+SpacesBeforeTrailingComments: 1
+Cpp11BracedListStyle: true
+Standard:        Cpp11
+IndentWidth:     4
+TabWidth:        8
+UseTab:          Never
+BreakBeforeBraces: Attach
+SpacesInParentheses: false
+SpacesInAngles:  false
+SpaceInEmptyParentheses: false
+SpacesInCStyleCastParentheses: false
+SpacesInContainerLiterals: true
+SpaceBeforeAssignmentOperators: true
+ContinuationIndentWidth: 4
+CommentPragmas:  '*'
+ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
+SpaceBeforeParens: ControlStatements
+DisableFormat:   false
+AlignAfterOpenBracket: false
+AlignEscapedNewlinesLeft: true
+...
+
diff --git a/cpputil/usid/CMakeLists.txt b/cpputil/usid/CMakeLists.txt
new file mode 100644
index 0000000..0e001dc
--- /dev/null
+++ b/cpputil/usid/CMakeLists.txt
@@ -0,0 +1,53 @@
+cmake_minimum_required(VERSION 3.14.5)
+project(usid LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 17)
+
+add_compile_options(-Wall)
+
+include(FetchContent)
+FetchContent_Declare(GridTools
+    GIT_REPOSITORY https://github.com/GridTools/gridtools.git
+    GIT_TAG        master
+)
+FetchContent_MakeAvailable(GridTools)
+
+add_library(usid_naive_helper INTERFACE)
+target_include_directories(usid_naive_helper INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include)
+target_link_libraries(usid_naive_helper INTERFACE GridTools::gridtools)
+
+include(CheckLanguage)
+check_language(CUDA)
+if(CMAKE_CUDA_COMPILER)
+    enable_language(CUDA)
+
+    add_library(usid_cuda_helper INTERFACE)
+    target_include_directories(usid_cuda_helper INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include)
+    target_link_libraries(usid_cuda_helper INTERFACE GridTools::gridtools GridTools::stencil_gpu)
+endif()
+
+set(BUILD_SHARED_LIBS OFF)
+
+FetchContent_Declare(GoogleTest
+    GIT_REPOSITORY https://github.com/google/googletest.git
+    GIT_TAG        release-1.10.0
+)
+FetchContent_MakeAvailable(GoogleTest)
+
+find_package(eckit REQUIRED)
+find_package(Atlas REQUIRED)
+
+add_library(fvm_nabla_driver INTERFACE)
+target_include_directories(fvm_nabla_driver INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include)
+target_link_libraries(fvm_nabla_driver INTERFACE atlas eckit GridTools::gridtools gtest)
+if(CMAKE_CUDA_COMPILER)
+    target_link_libraries(fvm_nabla_driver INTERFACE GridTools::stencil_gpu)
+endif()
+
+
+if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
+    include(CTest)
+    if(BUILD_TESTING)
+        add_subdirectory(tests)
+    endif()
+endif()
diff --git a/cpputil/usid/include/gridtools/usid/atlas.hpp b/cpputil/usid/include/gridtools/usid/atlas.hpp
new file mode 100644
index 0000000..ca5e7c8
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/atlas.hpp
@@ -0,0 +1,19 @@
+#pragma once
+
+#include <atlas/mesh/Connectivity.h>
+
+#include <gridtools/storage/builder.hpp>
+
+namespace atlas::mesh {
+template <class Connectivity, class MaxNeighbors>
+auto make_storage_producer(MaxNeighbors max_neighbors,
+                           Connectivity const &src) {
+  return [&src, max_neighbors](auto traits) {
+    return gridtools::storage::builder<decltype(traits)>
+            .template type<int>()
+            .dimensions(src.rows(), max_neighbors)
+            .initializer([&src](auto row, auto col) { return col < src.cols(row) ? src(row, col) : -1; })
+            .build();
+  };
+}
+} // namespace atlas::mesh
diff --git a/cpputil/usid/include/gridtools/usid/cuda_helpers.hpp b/cpputil/usid/include/gridtools/usid/cuda_helpers.hpp
new file mode 100644
index 0000000..7861453
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/cuda_helpers.hpp
@@ -0,0 +1,52 @@
+#pragma once
+
+#ifndef __CUDACC__
+#error Tried to compile CUDA code with a regular C++ compiler.
+#endif
+
+#include <gridtools/common/cuda_util.hpp>
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/common/pair.hpp>
+#include <gridtools/common/tuple_util.hpp>
+#include <gridtools/sid/allocator.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/storage/gpu.hpp>
+
+#include "dim.hpp"
+#include "helpers.hpp"
+
+namespace gridtools::usid::cuda {
+using traits_t = storage::gpu;
+
+inline auto make_allocator() {
+  return sid::device::make_cached_allocator(&cuda_util::cuda_malloc<char[]>);
+}
+
+template <class Kernel, class Ptr, class Strides, class... Neighbors>
+__global__ void kernel(int h_size, Ptr ptr_holder, Strides strides,
+                       Neighbors... neighbors) {
+  auto h = blockIdx.x * blockDim.x + threadIdx.x;
+  if (h >= h_size)
+    return;
+  Kernel()()(
+      sid::shifted(ptr_holder(), device::at_key<dim::h>(strides), h), strides,
+      tuple_util::device::make<pair>(neighbors.first(), neighbors.second)...);
+}
+
+template <class Kernel, class Size, class Sid, class... Sids>
+void call_kernel(Size size, Sid &&fields, Sids &&...neighbor_fields) {
+  int threads_per_block = 32;
+  int blocks = (size + threads_per_block - 1) / threads_per_block;
+  kernel<Kernel><<<blocks, threads_per_block>>>(
+      size, sid::get_origin(fields), sid::get_strides(fields),
+      tuple_util::make<pair>(
+          sid::get_origin(neighbor_fields),
+          at_key<dim::h>(sid::get_strides(neighbor_fields)))...);
+  GT_CUDA_CHECK(cudaGetLastError());
+}
+
+template <class Tag, class Ptr>
+__device__ decltype(auto) field(Ptr const &ptr) {
+  return *device::at_key<Tag>(ptr);
+}
+} // namespace gridtools::usid::cuda
diff --git a/cpputil/usid/include/gridtools/usid/dim.hpp b/cpputil/usid/include/gridtools/usid/dim.hpp
new file mode 100644
index 0000000..48214f8
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/dim.hpp
@@ -0,0 +1,16 @@
+#pragma once
+
+#include <gridtools/common/defs.hpp>
+#include <gridtools/common/integral_constant.hpp>
+
+namespace gridtools::usid::dim {
+using horizontal = integral_constant<int_t, 0>;
+using vertical = integral_constant<int_t, 1>;
+using neighbor = integral_constant<int_t, 1>;
+using sparse = integral_constant<int_t, 2>;
+
+using h = horizontal;
+using k = vertical;
+using n = neighbor;
+using s = sparse;
+} // namespace gridtools::usid::dim
diff --git a/cpputil/usid/include/gridtools/usid/domain.hpp b/cpputil/usid/include/gridtools/usid/domain.hpp
new file mode 100644
index 0000000..b709e93
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/domain.hpp
@@ -0,0 +1,10 @@
+#pragma once
+
+namespace gridtools::usid {
+struct domain {
+  int vertex;
+  int edge;
+  int cell;
+  int k;
+};
+} // namespace gridtools::usid
diff --git a/cpputil/usid/include/gridtools/usid/helpers.hpp b/cpputil/usid/include/gridtools/usid/helpers.hpp
new file mode 100644
index 0000000..e2f3b01
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/helpers.hpp
@@ -0,0 +1,88 @@
+#pragma once
+
+#include <limits>
+#include <type_traits>
+#include <utility>
+
+#include <gridtools/common/defs.hpp>
+#include <gridtools/common/host_device.hpp>
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/common/integral_constant.hpp>
+#include <gridtools/common/utility.hpp>
+#include <gridtools/meta/id.hpp>
+#include <gridtools/sid/composite.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/sid/contiguous.hpp>
+#include <gridtools/sid/loop.hpp>
+#include <gridtools/sid/rename_dimensions.hpp>
+
+#include "dim.hpp"
+#include "domain.hpp"
+
+namespace gridtools::usid {
+template <class T, class Alloc, class HSize, class KSize>
+auto make_simple_tmp_storage(HSize h_size, KSize k_size, Alloc &alloc) {
+  return sid::make_contiguous<T>(
+      alloc, hymap::keys<dim::h, dim::k>::values<HSize, KSize>(h_size, k_size));
+}
+
+template <class T, class Tag, int_t N, bool HasSkipValues, class F, class Init,
+          class G, class Ptr, class Strides, class Neighbors>
+GT_FUNCTION T fold_neighbors(F f, Init init, G g, Ptr &&ptr, Strides &&strides,
+                             Neighbors &&neighbors) {
+  T acc = init(meta::lazy::id<T>());
+  sid::make_loop<Tag>(
+      integral_constant<int_t, N>())([&](auto const &ptr, auto &&) {
+    auto i = *host_device::at_key<Tag>(ptr);
+    if constexpr (HasSkipValues)
+      if (i < 0)
+        return;
+    acc = f(acc, g(ptr, sid::shifted(neighbors.first, neighbors.second, i)));
+  })(wstd::forward<decltype(ptr)>(ptr), strides);
+  return acc;
+}
+
+template <class T, class Tag, int_t N, bool HasSkip = true, class F,
+          class... Args>
+GT_FUNCTION T sum_neighbors(F f, Args &&...args) {
+  return fold_neighbors<T, Tag, N, HasSkip>(
+      [](auto x, auto y) { return x + y; },
+      [](auto z) -> typename decltype(z)::type { return 0; }, f,
+      wstd::forward<Args>(args)...);
+}
+
+template <class T, class Tag, int_t N, bool HasSkip = true, class F,
+          class... Args>
+GT_FUNCTION T mul_neighbors(F f, Args &&...args) {
+  return fold_neighbors<T, Tag, N, HasSkip>(
+      [](auto x, auto y) { return x * y; },
+      [](auto z) -> typename decltype(z)::type { return 1; }, f,
+      wstd::forward<Args>(args)...);
+}
+
+template <class T, class Tag, int_t N, bool HasSkip = true, class F,
+          class... Args>
+GT_FUNCTION T min_neighbors(F f, Args &&...args) {
+  return fold_neighbors<T, Tag, N, HasSkip>(
+      [](auto x, auto y) { return x < y ? x : y; },
+      [](auto z) -> typename decltype(z)::type {
+        constexpr auto res =
+            std::numeric_limits<typename decltype(z)::type>::max();
+        return res;
+      },
+      f, wstd::forward<Args>(args)...);
+}
+
+template <class T, class Tag, int_t N, bool HasSkip = true, class F,
+          class... Args>
+GT_FUNCTION T max_neighbors(F f, Args &&...args) {
+  return fold_neighbors<T, Tag, N, HasSkip>(
+      [](auto x, auto y) { return x > y ? x : y; },
+      [](auto z) -> typename decltype(z)::type {
+        constexpr auto res =
+            std::numeric_limits<typename decltype(z)::type>::min();
+        return res;
+      },
+      f, wstd::forward<Args>(args)...);
+}
+} // namespace gridtools::usid
diff --git a/cpputil/usid/include/gridtools/usid/naive_helpers.hpp b/cpputil/usid/include/gridtools/usid/naive_helpers.hpp
new file mode 100644
index 0000000..c094bb4
--- /dev/null
+++ b/cpputil/usid/include/gridtools/usid/naive_helpers.hpp
@@ -0,0 +1,38 @@
+#pragma once
+
+#include <memory>
+#include <tuple>
+#include <utility>
+
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/sid/allocator.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/sid/loop.hpp>
+#include <gridtools/storage/cpu_ifirst.hpp>
+
+#include "dim.hpp"
+#include "helpers.hpp"
+
+namespace gridtools::usid::naive {
+using traits_t = storage::cpu_ifirst;
+
+inline auto make_allocator() {
+  return sid::make_cached_allocator(&std::make_unique<char[]>);
+}
+
+template <class Kernel, class Size, class Sid, class... Sids>
+void call_kernel(Size size, Sid &&fields, Sids &&...neighbor_fields) {
+  sid::make_loop<dim::h>(size)(
+      [params = std::make_tuple(std::make_pair(
+           sid::get_origin(neighbor_fields)(),
+           sid::get_stride<dim::h>(sid::get_strides(neighbor_fields)))...)](
+          auto &ptr, auto const &strides) {
+        std::apply(Kernel()(),
+                   std::tuple_cat(std::forward_as_tuple(ptr, strides), params));
+      })(sid::get_origin(fields)(), sid::get_strides(fields));
+}
+
+template <class Tag, class Ptr> decltype(auto) field(Ptr const &ptr) {
+  return *at_key<Tag>(ptr);
+}
+} // namespace gridtools::usid::naive
diff --git a/cpputil/usid/include/tests/fvm_nabla_driver.hpp b/cpputil/usid/include/tests/fvm_nabla_driver.hpp
new file mode 100644
index 0000000..752b42b
--- /dev/null
+++ b/cpputil/usid/include/tests/fvm_nabla_driver.hpp
@@ -0,0 +1,176 @@
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <tuple>
+
+#include <gtest/gtest.h>
+
+#include <atlas/array.h>
+#include <atlas/functionspace.h>
+#include <atlas/grid.h>
+#include <atlas/mesh.h>
+#include <atlas/mesh/actions/BuildDualMesh.h>
+#include <atlas/mesh/actions/BuildEdges.h>
+#include <atlas/meshgenerator.h>
+
+#include <gridtools/common/integral_constant.hpp>
+#include <gridtools/storage/builder.hpp>
+#include <gridtools/storage/sid.hpp>
+
+#include <gridtools/usid/atlas.hpp>
+
+#ifdef __CUDACC__
+#include <gridtools/storage/gpu.hpp>
+namespace fvm_nabla_driver_impl_ {
+    using storage_traits_t = gridtools::storage::gpu;
+}
+#else
+#include <gridtools/storage/cpu_ifirst.hpp>
+namespace fvm_nabla_driver_impl_ {
+    using storage_traits_t = gridtools::storage::cpu_ifirst;
+}
+#endif
+
+namespace fvm_nabla_driver_impl_ {
+    using namespace gridtools::literals;
+
+    inline constexpr auto storage_builder = gridtools::storage::builder<storage_traits_t>;
+
+    inline auto make_mesh() {
+        using namespace atlas;
+        using param_t = StructuredMeshGenerator::Parameters;
+        auto res = StructuredMeshGenerator(param_t("triangulate", true) | param_t("angle", -1.)).generate(Grid("O32"));
+        functionspace::EdgeColumns(res, option::levels(1) | option::halo(1));
+        functionspace::NodeColumns(res, option::levels(1) | option::halo(1));
+        mesh::actions::build_edges(res);
+        mesh::actions::build_node_to_edge_connectivity(res);
+        mesh::actions::build_median_dual_mesh(res);
+        return res;
+    }
+
+    inline const double rpi = 2 * std::asin(1);
+    inline const double radius = 6371.22e+03;
+    inline const double deg2rad = 2 * rpi / 360;
+    inline constexpr int MXX = 0;
+    inline constexpr int MYY = 1;
+    inline constexpr auto edges_per_node = 7_c;
+
+    inline auto make_vol(atlas::mesh::Nodes const &nodes) {
+        auto dual_volumes = atlas::array::make_view<double, 1>(nodes.field("dual_volumes"));
+        auto init = [&](int n) { return dual_volumes(n) * (std::pow(deg2rad, 2) * std::pow(radius, 2)); };
+        return storage_builder.type<double>().dimensions(nodes.size()).initializer(init).build();
+    }
+
+    inline auto make_sign(atlas::Mesh const &mesh) {
+        auto &&edges = mesh.edges();
+        auto &&nodes = mesh.nodes();
+        auto &&n2e = nodes.edge_connectivity();
+        auto &&e2n = edges.node_connectivity();
+        auto flags = atlas::array::make_view<int, 1>(edges.flags());
+        auto is_pole_edge = [&](auto e) {
+            using topology_t = atlas::mesh::Nodes::Topology;
+            return topology_t::check(flags(e), topology_t::POLE);
+        };
+        auto init = [&](int n, int, int e) -> double {
+            if (e >= n2e.cols(n))
+                return 0;
+            auto ee = n2e(n, e);
+            return n == e2n(ee, 0) || is_pole_edge(ee) ? 1 : -1;
+        };
+        return storage_builder.type<double>()
+            .selector<1, 0, 1>()
+            .dimensions(nodes.size(), 1, edges_per_node)
+            .initializer(init)
+            .build();
+    }
+
+    inline auto make_S_MXX(atlas::mesh::Edges const &edges) {
+        auto dual_normals = atlas::array::make_view<double, 2>(edges.field("dual_normals"));
+        return storage_builder.type<double>()
+            .dimensions(edges.size())
+            .initializer([&](int i) { return dual_normals(i, MXX) * radius * deg2rad; })
+            .build();
+    }
+
+    inline auto make_S_MYY(atlas::mesh::Edges const &edges) {
+        auto dual_normals = atlas::array::make_view<double, 2>(edges.field("dual_normals"));
+        return storage_builder.type<double>()
+            .dimensions(edges.size())
+            .initializer([&](int i) { return dual_normals(i, MYY) * radius * deg2rad; })
+            .build();
+    }
+
+    // TODO ask Christian for a proper name for this input data
+    inline auto make_pp(atlas::mesh::Nodes const &nodes) {
+        static const double zh0 = 2000;
+        static const double zrad = 3 * rpi / 4 * radius;
+        static const double zeta = rpi / 16 * radius;
+        static const double zlatc = 0;
+        static const double zlonc = 3 * rpi / 2;
+
+        auto lonlat = atlas::array::make_view<double, 2>(nodes.field("lonlat"));
+        // lonlatcr is in physical space and may differ from coords later
+        auto rlonlatcr = [&](int n, int i) { return lonlat(n, i) * deg2rad; };
+        auto init = [&](int n) {
+            double zlon = rlonlatcr(n, MXX);
+            double rcosa = std::cos(rlonlatcr(n, MYY));
+            double rsina = std::sin(rlonlatcr(n, MYY));
+            double zdist = std::sin(zlatc) * rsina + std::cos(zlatc) * rcosa * std::cos(zlon - zlonc);
+            zdist = radius * std::acos(zdist);
+            return zdist < zrad
+                       ? .5 * zh0 * (1 + std::cos(rpi * zdist / zrad)) * std::pow(std::cos(rpi * zdist / zeta), 2)
+                       : 0;
+        };
+        return storage_builder.type<double>().dimensions(nodes.size()).initializer(init).build();
+    }
+
+    template <class T>
+    auto min_max(T const &field) {
+        double min = std::numeric_limits<double>::max();
+        double max = std::numeric_limits<double>::min();
+        auto view = field->const_host_view();
+        auto lengths = field->lengths();
+        for (int i = 0; i < (int)lengths[0]; ++i)
+            for (int k = 0; k < (int)lengths[1]; ++k) {
+                double x = view(i, k);
+                min = std::min(min, x);
+                max = std::max(max, x);
+            }
+        return std::make_tuple(min, max);
+    }
+
+    template <class Nabla>
+    void fvm_nabla_driver(Nabla nabla) {
+        constexpr auto k = 1_c;
+
+        auto mesh = make_mesh();
+        auto &&edges = mesh.edges();
+        auto &&nodes = mesh.nodes();
+
+        // output
+        auto make_output = storage_builder.type<double>().dimensions(nodes.size(), k);
+        auto pnabla_MXX = make_output();
+        auto pnabla_MYY = make_output();
+
+        nabla({nodes.size(), edges.size(), mesh.cells().size(), k},
+            make_storage_producer(edges_per_node, nodes.edge_connectivity()),
+            make_storage_producer(2_c, edges.node_connectivity()))(make_S_MXX(edges),
+            make_S_MYY(edges),
+            make_pp(nodes),
+            pnabla_MXX,
+            pnabla_MYY,
+            make_vol(nodes),
+            make_sign(mesh));
+
+        auto [x_min, x_max] = min_max(pnabla_MXX);
+        auto [y_min, y_max] = min_max(pnabla_MYY);
+
+        EXPECT_DOUBLE_EQ(-3.5455427772566003E-003, x_min);
+        EXPECT_DOUBLE_EQ(3.5455427772565435E-003, x_max);
+        EXPECT_DOUBLE_EQ(-3.3540113705465301E-003, y_min);
+        EXPECT_DOUBLE_EQ(3.3540113705465301E-003, y_max);
+    } //
+} // namespace fvm_nabla_driver_impl_
+using fvm_nabla_driver_impl_::fvm_nabla_driver;
diff --git a/cpputil/usid/tests/CMakeLists.txt b/cpputil/usid/tests/CMakeLists.txt
new file mode 100644
index 0000000..5f0c7df
--- /dev/null
+++ b/cpputil/usid/tests/CMakeLists.txt
@@ -0,0 +1,11 @@
+add_executable(nabla_naive nabla_naive.cpp)
+target_link_libraries(nabla_naive PRIVATE usid_naive_helper fvm_nabla_driver gtest_main)
+add_test(NAME nabla_naive COMMAND $<TARGET_FILE:nabla_naive>)
+
+if(CMAKE_CUDA_COMPILER)
+    add_executable(nabla_cuda nabla_cuda.cu)
+    target_link_libraries(nabla_cuda PRIVATE usid_cuda_helper fvm_nabla_driver gtest_main)
+    gridtools_setup_target(nabla_cuda CUDA_ARCH sm_50)
+    target_compile_options(nabla_cuda PRIVATE "-std=c++17")
+    add_test(NAME nabla_cuda COMMAND $<TARGET_FILE:nabla_cuda>)
+endif()
diff --git a/cpputil/usid/tests/nabla_cuda.cu b/cpputil/usid/tests/nabla_cuda.cu
new file mode 100644
index 0000000..a5784f9
--- /dev/null
+++ b/cpputil/usid/tests/nabla_cuda.cu
@@ -0,0 +1,5 @@
+#include "nabla_cuda.hpp"
+#include <gtest/gtest.h>
+#include <tests/fvm_nabla_driver.hpp>
+
+TEST(fvm, nabla_cuda) { fvm_nabla_driver(nabla); }
diff --git a/cpputil/usid/tests/nabla_cuda.hpp b/cpputil/usid/tests/nabla_cuda.hpp
new file mode 100644
index 0000000..3a57185
--- /dev/null
+++ b/cpputil/usid/tests/nabla_cuda.hpp
@@ -0,0 +1,72 @@
+#pragma once
+#include <gridtools/usid/usid_cuda_helpers.hpp>
+namespace gridtools::usid::cuda::nabla_impl_ {
+    struct v2e_tag;
+    struct e2v_tag;
+    struct S_MXX_tag;
+    struct S_MYY_tag;
+    struct zavgS_MXX_tag;
+    struct zavgS_MYY_tag;
+    struct pnabla_MXX_tag;
+    struct pnabla_MYY_tag;
+    struct vol_tag;
+    struct sign_tag;
+    struct pp_tag;
+    struct kernel_0 {
+        GT_FUNCTION auto operator()() const {
+            return [](auto &&ptr, auto &&strides, auto &&neighbors) {
+                auto zavg = 0.5 * sum_neighbors<double, e2v_tag, 2, false>(
+                                      [](auto &&, auto &&n) { return field<pp_tag>(n); }, ptr, strides, neighbors);
+                field<zavgS_MXX_tag>(ptr) = field<S_MXX_tag>(ptr) * zavg;
+                field<zavgS_MYY_tag>(ptr) = field<S_MYY_tag>(ptr) * zavg;
+            };
+        }
+    };
+    struct kernel_1 {
+        GT_FUNCTION auto operator()() const {
+            return [](auto &&ptr, auto &&strides, auto &&neighbors) {
+                field<pnabla_MXX_tag>(ptr) = sum_neighbors<double, v2e_tag, 7>(
+                    [](auto &&p, auto &&n) { return field<zavgS_MXX_tag>(n) * field<sign_tag>(p); },
+                    ptr,
+                    strides,
+                    neighbors);
+                field<pnabla_MYY_tag>(ptr) = sum_neighbors<double, v2e_tag, 7>(
+                    [](auto &&p, auto &&n) { return field<zavgS_MYY_tag>(n) * field<sign_tag>(p); },
+                    ptr,
+                    strides,
+                    neighbors);
+                field<pnabla_MXX_tag>(ptr) = field<pnabla_MXX_tag>(ptr) / field<vol_tag>(ptr);
+                field<pnabla_MYY_tag>(ptr) = field<pnabla_MYY_tag>(ptr) / field<vol_tag>(ptr);
+            };
+        }
+    };
+    inline constexpr auto nabla = [](domain d, auto &&v2e, auto &&e2v) {
+        static_assert(is_sid<decltype(v2e(traits_t()))>());
+        static_assert(is_sid<decltype(e2v(traits_t()))>());
+        return
+            [d = std::move(d),
+                v2e = sid::rename_dimensions<dim::n, v2e_tag>(std::forward<decltype(v2e)>(v2e)(traits_t())),
+                e2v = sid::rename_dimensions<dim::n, e2v_tag>(std::forward<decltype(e2v)>(e2v)(traits_t()))](
+                auto &&S_MXX, auto &&S_MYY, auto &&pp, auto &&pnabla_MXX, auto &&pnabla_MYY, auto &&vol, auto &&sign) {
+                static_assert(is_sid<decltype(S_MXX)>());
+                static_assert(is_sid<decltype(S_MYY)>());
+                static_assert(is_sid<decltype(pp)>());
+                static_assert(is_sid<decltype(pnabla_MXX)>());
+                static_assert(is_sid<decltype(pnabla_MYY)>());
+                static_assert(is_sid<decltype(vol)>());
+                static_assert(is_sid<decltype(sign)>());
+                auto alloc = make_allocator();
+                auto zavgS_MXX = make_simple_tmp_storage<double>(d.edge, d.k, alloc);
+                auto zavgS_MYY = make_simple_tmp_storage<double>(d.edge, d.k, alloc);
+                call_kernel<kernel_0>(d.edge,
+                    sid::composite::make<e2v_tag, S_MXX_tag, S_MYY_tag, zavgS_MXX_tag, zavgS_MYY_tag>(
+                        e2v, S_MXX, S_MYY, zavgS_MXX, zavgS_MYY),
+                    sid::composite::make<pp_tag>(pp));
+                call_kernel<kernel_1>(d.vertex,
+                    sid::composite::make<v2e_tag, pnabla_MXX_tag, pnabla_MYY_tag, sign_tag, vol_tag>(
+                        v2e, pnabla_MXX, pnabla_MYY, sid::rename_dimensions<dim::s, v2e_tag>(sign), vol),
+                    sid::composite::make<zavgS_MXX_tag, zavgS_MYY_tag>(zavgS_MXX, zavgS_MYY));
+            };
+    };
+} // namespace gridtools::usid::cuda::nabla_impl_
+using gridtools::usid::cuda::nabla_impl_::nabla;
diff --git a/cpputil/usid/tests/nabla_naive.cpp b/cpputil/usid/tests/nabla_naive.cpp
new file mode 100644
index 0000000..a3cca65
--- /dev/null
+++ b/cpputil/usid/tests/nabla_naive.cpp
@@ -0,0 +1,5 @@
+#include "nabla_naive.hpp"
+#include <gtest/gtest.h>
+#include <tests/fvm_nabla_driver.hpp>
+
+TEST(fvm, nabla_naive) { fvm_nabla_driver(nabla); }
diff --git a/cpputil/usid/tests/nabla_naive.hpp b/cpputil/usid/tests/nabla_naive.hpp
new file mode 100644
index 0000000..8e4e115
--- /dev/null
+++ b/cpputil/usid/tests/nabla_naive.hpp
@@ -0,0 +1,66 @@
+#pragma once
+#include <gridtools/usid/naive_helpers.hpp>
+namespace gridtools::usid::naive::nabla_impl_ {
+    struct v2e_tag;
+    struct e2v_tag;
+    struct S_MXX_tag;
+    struct S_MYY_tag;
+    struct zavgS_MXX_tag;
+    struct zavgS_MYY_tag;
+    struct pnabla_MXX_tag;
+    struct pnabla_MYY_tag;
+    struct vol_tag;
+    struct sign_tag;
+    struct pp_tag;
+    struct kernel_0 {
+        GT_FUNCTION auto operator()() const {
+            return [](auto &&e, auto &&strides, auto &&v) {
+                auto zavg = 0.5 * sum_neighbors<double, e2v_tag, 2, false>(
+                                      [](auto &&e, auto &&v) { return field<pp_tag>(v); }, e, strides, v);
+                field<zavgS_MXX_tag>(e) = field<S_MXX_tag>(e) * zavg;
+                field<zavgS_MYY_tag>(e) = field<S_MYY_tag>(e) * zavg;
+            };
+        }
+    };
+    struct kernel_1 {
+        GT_FUNCTION auto operator()() const {
+            return [](auto &&v, auto &&strides, auto &&e) {
+                field<pnabla_MXX_tag>(v) = sum_neighbors<double, v2e_tag, 7>(
+                    [](auto &&v, auto &&e) { return field<zavgS_MXX_tag>(e) * field<sign_tag>(v); }, v, strides, e);
+                field<pnabla_MYY_tag>(v) = sum_neighbors<double, v2e_tag, 7>(
+                    [](auto &&v, auto &&e) { return field<zavgS_MYY_tag>(e) * field<sign_tag>(v); }, v, strides, e);
+                field<pnabla_MXX_tag>(v) = field<pnabla_MXX_tag>(v) / field<vol_tag>(v);
+                field<pnabla_MYY_tag>(v) = field<pnabla_MYY_tag>(v) / field<vol_tag>(v);
+            };
+        }
+    };
+    inline constexpr auto nabla = [](domain d, auto &&v2e, auto &&e2v) {
+        static_assert(is_sid<decltype(v2e(traits_t()))>());
+        static_assert(is_sid<decltype(e2v(traits_t()))>());
+        return
+            [d = std::move(d),
+                v2e = sid::rename_dimensions<dim::n, v2e_tag>(std::forward<decltype(v2e)>(v2e)(traits_t())),
+                e2v = sid::rename_dimensions<dim::n, e2v_tag>(std::forward<decltype(e2v)>(e2v)(traits_t()))](
+                auto &&S_MXX, auto &&S_MYY, auto &&pp, auto &&pnabla_MXX, auto &&pnabla_MYY, auto &&vol, auto &&sign) {
+                static_assert(is_sid<decltype(S_MXX)>());
+                static_assert(is_sid<decltype(S_MYY)>());
+                static_assert(is_sid<decltype(pp)>());
+                static_assert(is_sid<decltype(pnabla_MXX)>());
+                static_assert(is_sid<decltype(pnabla_MYY)>());
+                static_assert(is_sid<decltype(vol)>());
+                static_assert(is_sid<decltype(sign)>());
+                auto alloc = make_allocator();
+                auto zavgS_MXX = make_simple_tmp_storage<double>(d.edge, d.k, alloc);
+                auto zavgS_MYY = make_simple_tmp_storage<double>(d.edge, d.k, alloc);
+                call_kernel<kernel_0>(d.edge,
+                    sid::composite::make<e2v_tag, S_MXX_tag, S_MYY_tag, zavgS_MXX_tag, zavgS_MYY_tag>(
+                        e2v, S_MXX, S_MYY, zavgS_MXX, zavgS_MYY),
+                    sid::composite::make<pp_tag>(pp));
+                call_kernel<kernel_1>(d.vertex,
+                    sid::composite::make<v2e_tag, pnabla_MXX_tag, pnabla_MYY_tag, sign_tag, vol_tag>(
+                        v2e, pnabla_MXX, pnabla_MYY, sid::rename_dimensions<dim::s, v2e_tag>(sign), vol),
+                    sid::composite::make<zavgS_MXX_tag, zavgS_MYY_tag>(zavgS_MXX, zavgS_MYY));
+            };
+    };
+} // namespace gridtools::usid::naive::nabla_impl_
+using gridtools::usid::naive::nabla_impl_::nabla;