diff --git a/model/atmosphere/advection/tests/advection/integration_tests/test_advection.py b/model/atmosphere/advection/tests/advection/integration_tests/test_advection.py
index fc17a97d3c..98f9f28929 100644
--- a/model/atmosphere/advection/tests/advection/integration_tests/test_advection.py
+++ b/model/atmosphere/advection/tests/advection/integration_tests/test_advection.py
@@ -21,6 +21,7 @@
 from icon4py.model.common.utils import data_allocation as data_alloc
 from icon4py.model.testing import (
     definitions,
+    exchange_utils,
     grid_utils,
     grid_utils as gridtest_utils,
     serialbox as sb,
@@ -47,7 +48,6 @@
     construct_least_squares_state,
     construct_metric_state,
     construct_prep_adv,
-    dummy_exchange,
     log_serialized,
     verify_advection_fields,
 )
@@ -162,7 +162,7 @@ def test_advection_run_single_step(
         ),
         min_rlcell_int=icon_grid.end_index(h_grid.domain(dims.CellDim)(h_grid.Zone.LOCAL)),
         geometry_type=icon_grid.geometry_type,
-        exchange=dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
     )
 
     least_squares_state = construct_least_squares_state(least_squares_coeffs, backend=backend)
@@ -270,7 +270,7 @@ def test_compute_lsq_coeffs(
         start_idx,
         min_rlcell_int,
         icon_grid.geometry_type,
-        exchange=dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
     )
 
     assert test_helpers.dallclose(
diff --git a/model/atmosphere/advection/tests/advection/utils.py b/model/atmosphere/advection/tests/advection/utils.py
index c06148e4f0..4cf2b66b3d 100644
--- a/model/atmosphere/advection/tests/advection/utils.py
+++ b/model/atmosphere/advection/tests/advection/utils.py
@@ -197,7 +197,3 @@ def verify_advection_fields(
         p_tracer_new_ref.asnumpy()[p_tracer_new_range, :],
         atol=1e-16,
     )
-
-
-def dummy_exchange(*field: data_alloc.NDArray) -> None:
-    return None
diff --git a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
index 011800753a..080ffac929 100644
--- a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
+++ b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
@@ -1094,7 +1094,7 @@ def test_compute_perturbed_quantities_and_interpolation(
     reference_theta_at_cells_on_half_levels = metrics_savepoint.theta_ref_ic()
     d2dexdz2_fac1_mc = metrics_savepoint.d2dexdz2_fac1_mc()
     d2dexdz2_fac2_mc = metrics_savepoint.d2dexdz2_fac2_mc()
-    wgtfacq_c = metrics_savepoint.wgtfacq_c_dsl()
+    wgtfacq_c = metrics_savepoint.wgtfacq_c()
     wgtfac_c = metrics_savepoint.wgtfac_c()
     exner_w_explicit_weight_parameter = metrics_savepoint.vwind_expl_wgt()
     ddz_of_reference_exner_at_cells_on_half_levels = metrics_savepoint.d_exner_dz_ref_ic()
@@ -1166,17 +1166,24 @@ def test_compute_perturbed_quantities_and_interpolation(
     assert test_utils.dallclose(
         perturbed_theta_v_at_cells_on_model_levels.asnumpy(), z_rth_pr_2_ref.asnumpy()
     )
+    # `z_exner_ex_pr` is only computed in a subset of the whole domain, reference may contain garbage outside this range
     assert test_utils.dallclose(
-        temporal_extrapolation_of_perturbed_exner.asnumpy(), z_exner_ex_pr_ref.asnumpy()
+        temporal_extrapolation_of_perturbed_exner.asnumpy()[
+            start_cell_lateral_boundary_level_3:end_cell_halo, :
+        ],
+        z_exner_ex_pr_ref.asnumpy()[start_cell_lateral_boundary_level_3:end_cell_halo, :],
     )
     assert test_utils.dallclose(
         perturbed_exner_at_cells_on_model_levels.asnumpy(), exner_pr_ref.asnumpy()
     )
     assert test_utils.dallclose(rho_at_cells_on_half_levels.asnumpy(), rho_ic_ref.asnumpy())
 
+    # `exner_at_cells_on_half_levels` is only computed in a subset of the whole domain, reference may contain garbage outside this range
     assert test_utils.dallclose(
-        exner_at_cells_on_half_levels.asnumpy()[:, nflatlev:],
-        z_exner_ic_ref.asnumpy()[:, nflatlev:],
+        exner_at_cells_on_half_levels.asnumpy()[
+            start_cell_lateral_boundary_level_3:end_cell_halo, nflatlev:
+        ],
+        z_exner_ic_ref.asnumpy()[start_cell_lateral_boundary_level_3:end_cell_halo, nflatlev:],
         rtol=1e-11,
     )
 
@@ -1768,7 +1775,7 @@ def test_compute_horizontal_velocity_quantities_and_fluxes(
     ddxn_z_full = metrics_savepoint.ddxn_z_full()
     ddxt_z_full = metrics_savepoint.ddxt_z_full()
     wgtfac_e = metrics_savepoint.wgtfac_e()
-    wgtfacq_e = metrics_savepoint.wgtfacq_e_dsl()
+    wgtfacq_e = metrics_savepoint.wgtfacq_e()
     rbf_vec_coeff_e = interpolation_savepoint.rbf_vec_coeff_e()
     geofac_grdiv = interpolation_savepoint.geofac_grdiv()
     nflatlev = vertical_params.nflatlev
@@ -2152,7 +2159,7 @@ def test_vertically_implicit_solver_at_predictor_step(
         reference_exner_at_cells_on_model_levels=metrics_savepoint.exner_ref_mc(),
         e_bln_c_s=interpolation_savepoint.e_bln_c_s(),
         wgtfac_c=metrics_savepoint.wgtfac_c(),
-        wgtfacq_c=metrics_savepoint.wgtfacq_c_dsl(),
+        wgtfacq_c=metrics_savepoint.wgtfacq_c(),
         iau_wgt_dyn=iau_wgt_dyn,
         dtime=savepoint_nonhydro_init.get_metadata("dtime").get("dtime"),
         is_iau_active=is_iau_active,
diff --git a/model/atmosphere/dycore/tests/dycore/integration_tests/test_velocity_advection.py b/model/atmosphere/dycore/tests/dycore/integration_tests/test_velocity_advection.py
index cacfda444c..ee02eda2b5 100644
--- a/model/atmosphere/dycore/tests/dycore/integration_tests/test_velocity_advection.py
+++ b/model/atmosphere/dycore/tests/dycore/integration_tests/test_velocity_advection.py
@@ -494,7 +494,7 @@ def test_compute_diagnostics_from_normal_wind(
     ddxn_z_full = metrics_savepoint.ddxn_z_full()
     ddxt_z_full = metrics_savepoint.ddxt_z_full()
     contravariant_correction_at_edges_on_model_levels = savepoint_velocity_init.z_w_concorr_me()
-    wgtfacq_e = metrics_savepoint.wgtfacq_e_dsl()
+    wgtfacq_e = metrics_savepoint.wgtfacq_e()
     nflatlev = grid_savepoint.nflatlev()
     c_intp = interpolation_savepoint.c_intp()
     inv_dual_edge_length = grid_savepoint.inv_dual_edge_length()
diff --git a/model/atmosphere/dycore/tests/dycore/utils.py b/model/atmosphere/dycore/tests/dycore/utils.py
index efe09d1b39..dccb5d755b 100644
--- a/model/atmosphere/dycore/tests/dycore/utils.py
+++ b/model/atmosphere/dycore/tests/dycore/utils.py
@@ -51,7 +51,7 @@ def construct_metric_state(
         time_extrapolation_parameter_for_exner=metrics_savepoint.exner_exfac(),
         reference_exner_at_cells_on_model_levels=metrics_savepoint.exner_ref_mc(),
         wgtfac_c=metrics_savepoint.wgtfac_c(),
-        wgtfacq_c=metrics_savepoint.wgtfacq_c_dsl(),
+        wgtfacq_c=metrics_savepoint.wgtfacq_c(),
         inv_ddqz_z_full=metrics_savepoint.inv_ddqz_z_full(),
         reference_rho_at_cells_on_model_levels=metrics_savepoint.rho_ref_mc(),
         reference_theta_at_cells_on_model_levels=metrics_savepoint.theta_ref_mc(),
@@ -71,7 +71,7 @@ def construct_metric_state(
         ddqz_z_full_e=metrics_savepoint.ddqz_z_full_e(),
         ddxt_z_full=metrics_savepoint.ddxt_z_full(),
         wgtfac_e=metrics_savepoint.wgtfac_e(),
-        wgtfacq_e=metrics_savepoint.wgtfacq_e_dsl(),
+        wgtfacq_e=metrics_savepoint.wgtfacq_e(),
         exner_w_implicit_weight_parameter=metrics_savepoint.vwind_impl_wgt(),
         horizontal_mask_for_3d_divdamp=metrics_savepoint.hmask_dd3d(),
         scaling_factor_for_3d_divdamp=metrics_savepoint.scalfac_dd3d(),
diff --git a/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_full_muphys.py b/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_full_muphys.py
index 3b9cb91672..da069257aa 100644
--- a/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_full_muphys.py
+++ b/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_full_muphys.py
@@ -16,6 +16,7 @@
 
 from icon4py.model.atmosphere.subgrid_scale_physics.muphys.driver import common, run_full_muphys
 from icon4py.model.common import dimension as dims, model_backends
+from icon4py.model.testing import test_utils
 from icon4py.model.testing.fixtures.datatest import backend_like
 
 from . import utils
@@ -114,10 +115,10 @@ def test_full_muphys(
     rtol = 1e-14
     atol = 1e-16
 
-    np.testing.assert_allclose(ref.qv.asnumpy(), out.qv.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qc.asnumpy(), out.qc.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qi.asnumpy(), out.qi.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qr.asnumpy(), out.qr.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qs.asnumpy(), out.qs.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qg.asnumpy(), out.qg.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.t.asnumpy(), out.t.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qv.asnumpy(), out.qv.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qc.asnumpy(), out.qc.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qi.asnumpy(), out.qi.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qr.asnumpy(), out.qr.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qs.asnumpy(), out.qs.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qg.asnumpy(), out.qg.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.t.asnumpy(), out.t.asnumpy(), atol=atol, rtol=rtol)
diff --git a/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_graupel_only.py b/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_graupel_only.py
index 19e228de60..101184eed0 100644
--- a/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_graupel_only.py
+++ b/model/atmosphere/subgrid_scale_physics/muphys/tests/muphys/integration_tests/test_graupel_only.py
@@ -16,6 +16,7 @@
 
 from icon4py.model.atmosphere.subgrid_scale_physics.muphys.driver import common, run_graupel_only
 from icon4py.model.common import dimension as dims, model_backends
+from icon4py.model.testing import test_utils
 from icon4py.model.testing.fixtures.datatest import backend_like
 
 from . import utils
@@ -106,10 +107,10 @@ def test_graupel_only(
     rtol = 1e-14
     atol = 1e-16
 
-    np.testing.assert_allclose(ref.qv.asnumpy(), out.qv.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qc.asnumpy(), out.qc.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qi.asnumpy(), out.qi.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qr.asnumpy(), out.qr.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qs.asnumpy(), out.qs.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.qg.asnumpy(), out.qg.asnumpy(), atol=atol, rtol=rtol)
-    np.testing.assert_allclose(ref.t.asnumpy(), out.t.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qv.asnumpy(), out.qv.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qc.asnumpy(), out.qc.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qi.asnumpy(), out.qi.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qr.asnumpy(), out.qr.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qs.asnumpy(), out.qs.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.qg.asnumpy(), out.qg.asnumpy(), atol=atol, rtol=rtol)
+    test_utils.assert_dallclose(ref.t.asnumpy(), out.t.asnumpy(), atol=atol, rtol=rtol)
diff --git a/model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp_dsl.py b/model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp.py
similarity index 96%
rename from model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp_dsl.py
rename to model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp.py
index dff144be43..645d4c9d47 100644
--- a/model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp_dsl.py
+++ b/model/common/src/icon4py/model/common/metrics/compute_zdiff_gradp.py
@@ -15,7 +15,7 @@
 from icon4py.model.common.utils import data_allocation as data_alloc
 
 
-def compute_zdiff_gradp_dsl(  # noqa: PLR0912 [too-many-branches]
+def compute_zdiff_gradp(  # noqa: PLR0912 [too-many-branches]
     e2c,
     z_mc: data_alloc.NDArray,
     c_lin_e: data_alloc.NDArray,
@@ -129,6 +129,9 @@ def compute_zdiff_gradp_dsl(  # noqa: PLR0912 [too-many-branches]
                         break
 
     vertoffset_gradp = vertidx_gradp - vertoffset_gradp
-    vertoffset_gradp[:horizontal_start_1, :, :] = 0.0
+
+    # TODO(havogt): NumpyDataProvider needs to be extended to support implict exchange.
+    exchange(zdiff_gradp[:, 0, :])
+    exchange(zdiff_gradp[:, 1, :])
 
     return zdiff_gradp, vertoffset_gradp
diff --git a/model/common/src/icon4py/model/common/metrics/metrics_factory.py b/model/common/src/icon4py/model/common/metrics/metrics_factory.py
index 0f9100ae76..2e3422b7f6 100644
--- a/model/common/src/icon4py/model/common/metrics/metrics_factory.py
+++ b/model/common/src/icon4py/model/common/metrics/metrics_factory.py
@@ -36,7 +36,7 @@
     compute_advection_metrics,
     compute_coeff_gradekin,
     compute_diffusion_metrics,
-    compute_zdiff_gradp_dsl,
+    compute_zdiff_gradp,
     metric_fields as mf,
     metrics_attributes as attrs,
     reference_atmosphere,
@@ -762,9 +762,9 @@ def _register_computed_fields(self) -> None:  # noqa: PLR0915 [too-many-statemen
         )
         self.register_provider(compute_horizontal_mask_for_3d_divdamp)
 
-        compute_zdiff_gradp_dsl_np = factory.NumpyDataProvider(
+        compute_zdiff_gradp_np = factory.NumpyDataProvider(
             func=functools.partial(
-                compute_zdiff_gradp_dsl.compute_zdiff_gradp_dsl,
+                compute_zdiff_gradp.compute_zdiff_gradp,
                 array_ns=self._xp,
                 exchange=functools.partial(
                     self._exchange.exchange, dims.EdgeDim, stream=decomposition.BLOCK
@@ -793,7 +793,7 @@ def _register_computed_fields(self) -> None:  # noqa: PLR0915 [too-many-statemen
                 ),
             },
         )
-        self.register_provider(compute_zdiff_gradp_dsl_np)
+        self.register_provider(compute_zdiff_gradp_np)
 
         coeff_gradekin = factory.NumpyDataProvider(
             func=functools.partial(
diff --git a/model/common/src/icon4py/model/common/utils/data_allocation.py b/model/common/src/icon4py/model/common/utils/data_allocation.py
index 6f3ee8d6a4..eee4ee10e3 100644
--- a/model/common/src/icon4py/model/common/utils/data_allocation.py
+++ b/model/common/src/icon4py/model/common/utils/data_allocation.py
@@ -223,14 +223,6 @@ def list2field(
     return gtx.as_field(domain, arr, allocator=allocator)
 
 
-def kflip_wgtfacq(
-    arr: NDArray,
-    domain: gtx.Domain,
-    allocator: gtx_typing.Allocator,
-) -> gtx.Field:
-    return gtx.as_field(domain, arr[:, ::-1], allocator=allocator)  # type: ignore [arg-type] # type "ndarray[Any, dtype[Any] | Any"; expected "NDArrayObject"
-
-
 def adjust_fortran_indices(inp: NDArray) -> NDArray:
     """For some Fortran arrays we need to subtract 1 to be compatible with Python indexing."""
     return inp - 1
diff --git a/model/common/src/icon4py/model/common/utils/field_utils.py b/model/common/src/icon4py/model/common/utils/field_utils.py
new file mode 100644
index 0000000000..6179ca7099
--- /dev/null
+++ b/model/common/src/icon4py/model/common/utils/field_utils.py
@@ -0,0 +1,74 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022-2024, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
+
+import array_api_compat
+from gt4py import next as gtx
+from gt4py.next import typing as gtx_typing
+
+
+def flip(field: gtx.Field, dim: gtx.Dimension, allocator: gtx_typing.Allocator) -> gtx.Field:
+    """Flip a field along a given dimension.
+
+    Args:
+        field: The field to flip.
+        dim: The dimension along which to flip the field.
+        allocator: Allocator to use for the output field.
+    """
+    # Note: `allocator` needs to be passed explicitly since GT4Py fields currently don't persist how they were allocated.
+    xp = array_api_compat.array_namespace(field.ndarray)
+    flipped_array = xp.flip(field.ndarray, axis=field.domain.dims.index(dim))
+    return gtx.as_field(field.domain, flipped_array, allocator=allocator)
+
+
+def index2offset(
+    index_field: gtx.Field, dim: gtx.Dimension, allocator: gtx_typing.Allocator
+) -> gtx.Field:
+    """Convert an index field to an offset field.
+
+    Note: Additionally clips negative indices to become zero offset as Fortran initializes some indices with `0` (which corresponds to `-1` in Python) to indicate that they are not used.
+    As GT4Py's unstructured domain inference is incomplete and runs over the whole domain we might use out-of-bounds offsets in intermediate computations.
+
+    Args:
+        index_field: Index field in Python indexing (0-based).
+        dim: The dimension along which to convert indices to offsets.
+        allocator: Allocator to use for the output field.
+
+    Example:
+        >>> import numpy as np
+        >>> from gt4py import next as gtx
+        >>> from icon4py.model.common import dimension as dims
+        >>> index_field = gtx.as_field(
+        ...     {dims.CellDim: range(2, 6)}, np.array([5, 3, 4, 2], dtype=gtx.int32)
+        ... )
+        >>> result = index2offset(index_field, dims.CellDim, allocator=np)
+        >>> result.domain
+        Domain(dims=(Dimension(value='Cell', kind=<DimensionKind.HORIZONTAL: 'horizontal'>),), ranges=(UnitRange(2, 6),))
+        >>> result.ndarray
+        array([ 3,  0,  0, -3], dtype=int32)
+    """
+    # Note: `allocator` needs to be passed explicitly since GT4Py fields currently don't persist how they were allocated.
+    xp = array_api_compat.array_namespace(index_field.ndarray)
+
+    current_index = gtx.as_field(
+        gtx.Domain(index_field.domain[dim]),
+        xp.arange(
+            index_field.domain[dim].unit_range.start,
+            index_field.domain[dim].unit_range.stop,
+            dtype=index_field.ndarray.dtype,
+        ),
+        allocator=allocator,
+    )
+    # use GT4Py's broadcasting and field arithmetic (includes clipping)
+    offset_field = gtx.where(  # type: ignore[attr-defined]
+        index_field >= 0,  # type: ignore[operator]
+        index_field - current_index,
+        0,
+    )
+
+    # if GT4Py embedded would propagate the allocator, we could avoid this extra conversion.
+    return gtx.as_field(offset_field.domain, offset_field.ndarray, allocator=allocator)
diff --git a/model/common/tests/common/decomposition/unit_tests/test_halo.py b/model/common/tests/common/decomposition/unit_tests/test_halo.py
index 420f1c20ba..567965e4b1 100644
--- a/model/common/tests/common/decomposition/unit_tests/test_halo.py
+++ b/model/common/tests/common/decomposition/unit_tests/test_halo.py
@@ -130,7 +130,7 @@ def test_no_halo():
     decomposition = decomp.SingleNodeDecomposer()
     decomposition_info = halo_generator(decomposition(np.arange(grid_size.num_cells), 1))
     # cells
-    np.testing.assert_allclose(
+    test_utils.assert_dallclose(
         np.arange(grid_size.num_cells), decomposition_info.global_index(dims.CellDim)
     )
     assert np.all(decomposition_info.owner_mask(dims.CellDim))
@@ -138,7 +138,7 @@ def test_no_halo():
         decomposition_info.halo_levels(dims.CellDim) == definitions.DecompositionFlag.OWNED
     )
     # edges
-    np.testing.assert_allclose(
+    test_utils.assert_dallclose(
         np.arange(grid_size.num_edges), decomposition_info.global_index(dims.EdgeDim)
     )
     assert np.all(
@@ -146,7 +146,7 @@ def test_no_halo():
     )
     assert np.all(decomposition_info.owner_mask(dims.EdgeDim))
     # vertices
-    np.testing.assert_allclose(
+    test_utils.assert_dallclose(
         np.arange(grid_size.num_vertices), decomposition_info.global_index(dims.VertexDim)
     )
     assert np.all(
diff --git a/model/common/tests/common/grid/mpi_tests/test_parallel_grid_manager.py b/model/common/tests/common/grid/mpi_tests/test_parallel_grid_manager.py
index da17c4d217..a59dec248d 100644
--- a/model/common/tests/common/grid/mpi_tests/test_parallel_grid_manager.py
+++ b/model/common/tests/common/grid/mpi_tests/test_parallel_grid_manager.py
@@ -123,7 +123,7 @@ def check_local_global_field(
     check_halos: bool,
 ) -> None:
     if dim == dims.KDim:
-        np.testing.assert_allclose(global_reference_field, local_field)
+        test_utils.assert_dallclose(global_reference_field, local_field)
         return
 
     _log.info(
@@ -138,7 +138,7 @@ def check_local_global_field(
 
     # Compare halo against global reference field
     if check_halos:
-        np.testing.assert_allclose(
+        test_utils.assert_dallclose(
             global_reference_field[
                 decomposition_info.global_index(dim, decomp_defs.DecompositionInfo.EntryType.HALO)
             ],
@@ -180,7 +180,7 @@ def check_local_global_field(
             f" rank = {processor_props.rank}: SHAPES: global reference field {global_reference_field.shape}, gathered = {gathered_field.shape}"
         )
 
-        np.testing.assert_allclose(sorted_, global_reference_field, atol=1e-9, verbose=True)
+        test_utils.assert_dallclose(sorted_, global_reference_field, atol=1e-9, verbose=True)
 
 
 # These fields can't be computed with the embedded backend for one reason or
diff --git a/model/common/tests/common/grid/unit_tests/test_topography.py b/model/common/tests/common/grid/unit_tests/test_topography.py
index 0ed305d658..3982095682 100644
--- a/model/common/tests/common/grid/unit_tests/test_topography.py
+++ b/model/common/tests/common/grid/unit_tests/test_topography.py
@@ -13,11 +13,9 @@
 
 from icon4py.model.common.grid import topography as topo
 from icon4py.model.common.utils import data_allocation as data_alloc
-from icon4py.model.testing import definitions, test_utils
+from icon4py.model.testing import exchange_utils, test_utils
 from icon4py.model.testing.fixtures import *  # noqa: F403
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
@@ -53,7 +51,7 @@ def test_topography_smoothing_with_serialized_data(
         c2e2co=icon_grid.get_connectivity("C2E2CO").ndarray,
         num_iterations=num_iterations,
         array_ns=xp,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
     )
 
     assert test_utils.dallclose(topography_smoothed_ref, topography_smoothed, atol=1.0e-14)
diff --git a/model/common/tests/common/grid/unit_tests/test_vertical.py b/model/common/tests/common/grid/unit_tests/test_vertical.py
index 78bbeef5b5..19adb890c6 100644
--- a/model/common/tests/common/grid/unit_tests/test_vertical.py
+++ b/model/common/tests/common/grid/unit_tests/test_vertical.py
@@ -17,7 +17,7 @@
 from icon4py.model.common import dimension as dims, type_alias as ta
 from icon4py.model.common.grid import vertical as v_grid
 from icon4py.model.common.utils import data_allocation as data_alloc, device_utils
-from icon4py.model.testing import definitions, test_utils
+from icon4py.model.testing import definitions, exchange_utils, test_utils
 from icon4py.model.testing.fixtures import (
     backend,
     damping_height,
@@ -39,8 +39,6 @@
     topography_savepoint,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     from collections.abc import Iterator
@@ -430,7 +428,7 @@ def test_compute_vertical_coordinate(
         SLEVE_minimum_relative_layer_thickness_2=0.5,
         lowest_layer_thickness=vertical_config.lowest_layer_thickness,
         array_ns=xp,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
     )
 
     assert test_utils.dallclose(
diff --git a/model/common/tests/common/interpolation/unit_tests/test_interpolation_fields.py b/model/common/tests/common/interpolation/unit_tests/test_interpolation_fields.py
index e33277a4bb..95b3495608 100644
--- a/model/common/tests/common/interpolation/unit_tests/test_interpolation_fields.py
+++ b/model/common/tests/common/interpolation/unit_tests/test_interpolation_fields.py
@@ -32,7 +32,7 @@
     compute_pos_on_tplane_e_x_y_torus,
 )
 from icon4py.model.common.utils import data_allocation as data_alloc
-from icon4py.model.testing import definitions, serialbox as sb
+from icon4py.model.testing import definitions, exchange_utils, serialbox as sb
 from icon4py.model.testing.fixtures.datatest import (
     backend,
     data_provider,
@@ -44,8 +44,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 cell_domain = h_grid.domain(dims.CellDim)
 edge_domain = h_grid.domain(dims.EdgeDim)
@@ -61,7 +59,9 @@ def test_compute_c_lin_e(
     backend: gtx_typing.Backend,
 ) -> None:
     xp = data_alloc.import_array_ns(backend)
-    func = functools.partial(compute_c_lin_e, array_ns=xp, exchange=utils.dummy_exchange)
+    func = functools.partial(
+        compute_c_lin_e, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
+    )
     inv_dual_edge_length = grid_savepoint.inv_dual_edge_length()
     edge_cell_length = grid_savepoint.edge_cell_length()
     edge_owner_mask = grid_savepoint.e_owner_mask()
@@ -74,7 +74,7 @@ def test_compute_c_lin_e(
         inv_dual_edge_length.asnumpy(),
         edge_owner_mask.asnumpy(),
         horizontal_start,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=xp,
     )
     assert test_helpers.dallclose(c_lin_e, c_lin_e_ref.asnumpy())
@@ -162,7 +162,9 @@ def test_compute_geofac_n2s(
     e2c = icon_grid.get_connectivity(dims.E2C).ndarray
     c2e2c = icon_grid.get_connectivity(dims.C2E2C).ndarray
     horizontal_start = icon_grid.start_index(cell_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2))
-    geofac_n2s = functools.partial(compute_geofac_n2s, array_ns=xp, exchange=utils.dummy_exchange)(
+    geofac_n2s = functools.partial(
+        compute_geofac_n2s, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
+    )(
         dual_edge_length.ndarray,
         geofac_div.ndarray,
         c2e,
@@ -194,7 +196,7 @@ def test_compute_geofac_grg(
     horizontal_start = icon_grid.start_index(cell_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2))
 
     geofac_grg_0, geofac_grg_1 = functools.partial(
-        compute_geofac_grg, array_ns=xp, exchange=utils.dummy_exchange
+        compute_geofac_grg, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
     )(
         primal_normal_cell_x,
         primal_normal_cell_y,
@@ -238,7 +240,7 @@ def test_compute_geofac_grdiv(
     e2c2e = icon_grid.get_connectivity(dims.E2C2E).ndarray
     horizontal_start = icon_grid.start_index(edge_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2))
     geofac_grdiv = functools.partial(
-        compute_geofac_grdiv, array_ns=xp, exchange=utils.dummy_exchange
+        compute_geofac_grdiv, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
     )(
         geofac_div.ndarray,
         inv_dual_edge_length.ndarray,
@@ -284,7 +286,7 @@ def test_compute_c_bln_avg(
                 divergence_averaging_central_cell_weight,
                 horizontal_start,
                 horizontal_start_p2,
-                exchange=utils.dummy_exchange,
+                exchange=exchange_utils.dummy_exchange_with_bound_dim,
                 array_ns=xp,
             )
         case base_grid.GeometryType.TORUS:
@@ -295,7 +297,7 @@ def test_compute_c_bln_avg(
                 divergence_averaging_central_cell_weight,
                 horizontal_start,
                 horizontal_start_p2,
-                exchange=utils.dummy_exchange,
+                exchange=exchange_utils.dummy_exchange_with_bound_dim,
                 array_ns=xp,
             )
 
@@ -325,7 +327,9 @@ def test_compute_e_flx_avg(
     horizontal_start_1 = icon_grid.start_index(edge_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_4))
     horizontal_start_2 = icon_grid.start_index(edge_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_5))
 
-    e_flx_avg = functools.partial(compute_e_flx_avg, array_ns=xp, exchange=utils.dummy_exchange)(
+    e_flx_avg = functools.partial(
+        compute_e_flx_avg, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
+    )(
         c_bln_avg,
         geofac_div,
         owner_mask,
@@ -364,7 +368,7 @@ def test_compute_cells_aw_verts(
     )
 
     cells_aw_verts = functools.partial(
-        compute_cells_aw_verts, array_ns=xp, exchange=utils.dummy_exchange
+        compute_cells_aw_verts, array_ns=xp, exchange=exchange_utils.dummy_exchange_with_bound_dim
     )(
         dual_area=dual_area,
         edge_vert_length=edge_vert_length,
@@ -446,14 +450,14 @@ def test_compute_pos_on_tplane_e(
                 owner_mask,
                 e2c,
                 horizontal_start,
-                exchange=utils.dummy_exchange,
+                exchange=exchange_utils.dummy_exchange_with_bound_dim,
                 array_ns=xp,
             )
         case base_grid.GeometryType.TORUS:
             pos_on_tplane_e_x, pos_on_tplane_e_y = compute_pos_on_tplane_e_x_y_torus(
                 dual_edge_length,
                 e2c,
-                exchange=utils.dummy_exchange,
+                exchange=exchange_utils.dummy_exchange_with_bound_dim,
                 array_ns=xp,
             )
     assert test_helpers.dallclose(pos_on_tplane_e_x, pos_on_tplane_e_x_ref, atol=1e-8, rtol=1e-9)
diff --git a/model/common/tests/common/interpolation/unit_tests/test_rbf_interpolation.py b/model/common/tests/common/interpolation/unit_tests/test_rbf_interpolation.py
index c6bfb0000b..190852b068 100644
--- a/model/common/tests/common/interpolation/unit_tests/test_rbf_interpolation.py
+++ b/model/common/tests/common/interpolation/unit_tests/test_rbf_interpolation.py
@@ -20,6 +20,7 @@
 from icon4py.model.common.utils import data_allocation as data_alloc
 from icon4py.model.testing import (
     definitions,
+    exchange_utils,
     grid_utils as gridtest_utils,
     test_utils as test_helpers,
 )
@@ -34,8 +35,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
@@ -205,7 +204,7 @@ def test_rbf_interpolation_coeffs_cell(
         horizontal_end,
         grid.global_properties.domain_length,
         grid.global_properties.domain_height,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=data_alloc.import_array_ns(backend),
     )
 
@@ -282,7 +281,7 @@ def test_rbf_interpolation_coeffs_vertex(
         horizontal_end,
         grid.global_properties.domain_length,
         grid.global_properties.domain_height,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=data_alloc.import_array_ns(backend),
     )
 
@@ -361,7 +360,7 @@ def test_rbf_interpolation_coeffs_edge(
         horizontal_end,
         grid.global_properties.domain_length,
         grid.global_properties.domain_height,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=data_alloc.import_array_ns(backend),
     )
 
diff --git a/model/common/tests/common/metrics/mpi_tests/test_parallel_metrics.py b/model/common/tests/common/metrics/mpi_tests/test_parallel_metrics.py
index b5f6f5ddf5..f008add026 100644
--- a/model/common/tests/common/metrics/mpi_tests/test_parallel_metrics.py
+++ b/model/common/tests/common/metrics/mpi_tests/test_parallel_metrics.py
@@ -8,11 +8,13 @@
 
 from __future__ import annotations
 
+from types import EllipsisType
 from typing import TYPE_CHECKING
 
 import pytest
 
 from icon4py.model.common.decomposition import definitions as decomposition, mpi_decomposition
+from icon4py.model.common.grid import base as base_grid, horizontal as h_grid
 from icon4py.model.common.metrics import metrics_attributes as attrs, metrics_factory
 from icon4py.model.testing import definitions as test_defs, parallel_helpers, test_utils
 
@@ -34,7 +36,8 @@
 
 
 if TYPE_CHECKING:
-    import gt4py.next.typing as gtx_typing
+    from gt4py import next as gtx
+    from gt4py.next import typing as gtx_typing
 
     from icon4py.model.testing import serialbox as sb
 
@@ -43,21 +46,48 @@
     pytest.skip("Skipping parallel tests on single node installation", allow_module_level=True)
 
 
+def _get_slice_tuple_from_horizontal_range(
+    grid: base_grid.Grid,
+    horizontal_dim: gtx.Dimension,
+    horizontal_range: tuple[h_grid.Zone | None, h_grid.Zone | None],
+) -> tuple[slice | None | EllipsisType, ...]:
+    # TODO(havogt): Ideally we refactor the factories to only construct fields on the domain where they matter,
+    # then this function disappears as we get the verification range directly from the constructed field.
+    start_zone, end_zone = horizontal_range
+    horizontal_start = (
+        grid.start_index(h_grid.domain(horizontal_dim)(start_zone))
+        if start_zone is not None
+        else None
+    )
+    horizontal_end = (
+        grid.end_index(h_grid.domain(horizontal_dim)(end_zone)) if end_zone is not None else None
+    )
+    return (slice(horizontal_start, horizontal_end), ...)
+
+
 @pytest.mark.datatest
 @pytest.mark.mpi
 @pytest.mark.parametrize("processor_props", [True], indirect=True)
 @pytest.mark.parametrize(
-    "attrs_name, metrics_name",
+    "attrs_name, metrics_name, horizontal_range",
     [
-        (attrs.CELL_HEIGHT_ON_HALF_LEVEL, "z_ifc"),
-        (attrs.DDQZ_Z_FULL_E, "ddqz_z_full_e"),
-        (attrs.ZDIFF_GRADP, "zdiff_gradp"),
-        (attrs.VERTOFFSET_GRADP, "vertoffset_gradp"),
-        (attrs.Z_MC, "z_mc"),
-        (attrs.DDQZ_Z_HALF, "ddqz_z_half"),
-        (attrs.SCALING_FACTOR_FOR_3D_DIVDAMP, "scalfac_dd3d"),
-        (attrs.RAYLEIGH_W, "rayleigh_w"),
-        (attrs.COEFF_GRADEKIN, "coeff_gradekin"),
+        (attrs.CELL_HEIGHT_ON_HALF_LEVEL, "z_ifc", None),
+        (attrs.DDQZ_Z_FULL_E, "ddqz_z_full_e", None),
+        (
+            attrs.ZDIFF_GRADP,
+            "zdiff_gradp",
+            (h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2, None),
+        ),
+        (
+            attrs.VERTOFFSET_GRADP,
+            "vertoffset_gradp",
+            (h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2, None),
+        ),
+        (attrs.Z_MC, "z_mc", None),
+        (attrs.DDQZ_Z_HALF, "ddqz_z_half", None),
+        (attrs.SCALING_FACTOR_FOR_3D_DIVDAMP, "scalfac_dd3d", None),
+        (attrs.RAYLEIGH_W, "rayleigh_w", None),
+        (attrs.COEFF_GRADEKIN, "coeff_gradekin", None),
     ],
 )
 def test_distributed_metrics_attrs(
@@ -69,6 +99,7 @@ def test_distributed_metrics_attrs(
     metrics_factory_from_savepoint: metrics_factory.MetricsFieldsFactory,
     attrs_name: str,
     metrics_name: str,
+    horizontal_range: tuple[h_grid.Zone | None, h_grid.Zone | None] | None,
     experiment: test_defs.Experiment,
 ) -> None:
     if test_utils.is_embedded(backend) and metrics_name == "ddqz_z_half":
@@ -81,7 +112,15 @@ def test_distributed_metrics_attrs(
 
     field = factory.get(attrs_name).asnumpy()
     field_ref = metrics_savepoint.__getattribute__(metrics_name)().asnumpy()
-    assert test_utils.dallclose(field, field_ref, rtol=1e-8, atol=1.0e-8)
+    if horizontal_range is not None:
+        # We assume that the horizontal dimension exists and is the first one.
+        slicer = _get_slice_tuple_from_horizontal_range(
+            factory.grid, attrs.attrs[attrs_name]["dims"][0], horizontal_range
+        )
+        field = field[slicer]
+        field_ref = field_ref[slicer]
+
+    test_utils.assert_dallclose(field, field_ref, rtol=1e-8, atol=1.0e-8)
 
 
 @pytest.mark.datatest
@@ -188,7 +227,7 @@ def test_distributed_metrics_wgtfacq_e(
     factory = metrics_factory_from_savepoint
 
     field = factory.get(attrs.WGTFACQ_E).asnumpy()
-    field_ref = metrics_savepoint.wgtfacq_e_dsl().asnumpy()
+    field_ref = metrics_savepoint.wgtfacq_e().asnumpy()
     assert test_utils.dallclose(field, field_ref)
 
 
diff --git a/model/common/tests/common/metrics/unit_tests/test_compute_coeff_gradekin.py b/model/common/tests/common/metrics/unit_tests/test_compute_coeff_gradekin.py
index 24640fbad6..3a7b33790e 100644
--- a/model/common/tests/common/metrics/unit_tests/test_compute_coeff_gradekin.py
+++ b/model/common/tests/common/metrics/unit_tests/test_compute_coeff_gradekin.py
@@ -14,7 +14,7 @@
 import icon4py.model.common.grid.horizontal as h_grid
 from icon4py.model.common import dimension as dims
 from icon4py.model.common.metrics.compute_coeff_gradekin import compute_coeff_gradekin
-from icon4py.model.testing import test_utils
+from icon4py.model.testing import exchange_utils, test_utils
 from icon4py.model.testing.fixtures.datatest import (
     backend,
     data_provider,
@@ -26,8 +26,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     from icon4py.model.common.grid import base as base_grid
@@ -52,6 +50,6 @@ def test_compute_coeff_gradekin(
         edge_cell_length,
         inv_dual_edge_length,
         horizontal_start,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
     )
     assert test_utils.dallclose(coeff_gradekin_ref.asnumpy(), coeff_gradekin_full)
diff --git a/model/common/tests/common/metrics/unit_tests/test_compute_weight_factors.py b/model/common/tests/common/metrics/unit_tests/test_compute_weight_factors.py
index 746307ae84..7fcfc10064 100644
--- a/model/common/tests/common/metrics/unit_tests/test_compute_weight_factors.py
+++ b/model/common/tests/common/metrics/unit_tests/test_compute_weight_factors.py
@@ -14,7 +14,7 @@
 from icon4py.model.common import dimension as dims, type_alias as ta
 from icon4py.model.common.metrics import compute_weight_factors as weight_factors
 from icon4py.model.common.utils import data_allocation as data_alloc
-from icon4py.model.testing import test_utils
+from icon4py.model.testing import exchange_utils, test_utils
 from icon4py.model.testing.fixtures.datatest import (
     backend,
     data_provider,
@@ -27,8 +27,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
@@ -79,22 +77,22 @@ def test_compute_wgtfacq_e_dsl(
     icon_grid: base_grid.Grid,
     backend: gtx_typing.Backend | None,
 ) -> None:
-    wgtfacq_e_dsl_ref = metrics_savepoint.wgtfacq_e_dsl()
-    wgtfacq_c_dsl_ref = metrics_savepoint.wgtfacq_c_dsl()
+    wgtfacq_e_ref = metrics_savepoint.wgtfacq_e()
+    wgtfacq_c_ref = metrics_savepoint.wgtfacq_c()
 
     xp = data_alloc.import_array_ns(backend)
     wgtfacq_e_dsl = weight_factors.compute_wgtfacq_e_dsl(
         e2c=icon_grid.get_connectivity("E2C").ndarray,
         z_ifc=metrics_savepoint.z_ifc().ndarray,
-        wgtfacq_c_dsl=wgtfacq_c_dsl_ref.ndarray,
+        wgtfacq_c_dsl=wgtfacq_c_ref.ndarray,
         c_lin_e=interpolation_savepoint.c_lin_e().ndarray,
         n_edges=icon_grid.num_edges,
         nlev=icon_grid.num_levels,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=xp,
     )
 
-    assert test_utils.dallclose(data_alloc.as_numpy(wgtfacq_e_dsl), wgtfacq_e_dsl_ref.asnumpy())
+    assert test_utils.dallclose(data_alloc.as_numpy(wgtfacq_e_dsl), wgtfacq_e_ref.asnumpy())
 
 
 @pytest.mark.datatest
@@ -103,7 +101,7 @@ def test_compute_wgtfacq_c_dsl(
     metrics_savepoint: sb.MetricSavepoint,
     backend: gtx_typing.Backend | None,
 ) -> None:
-    wgtfacq_c_dsl_ref = metrics_savepoint.wgtfacq_c_dsl()
+    wgtfacq_c_ref = metrics_savepoint.wgtfacq_c()
 
     xp = data_alloc.import_array_ns(backend)
     wgtfacq_c_dsl = weight_factors.compute_wgtfacq_c_dsl(
@@ -111,4 +109,4 @@ def test_compute_wgtfacq_c_dsl(
         nlev=icon_grid.num_levels,
         array_ns=xp,
     )
-    assert test_utils.dallclose(data_alloc.as_numpy(wgtfacq_c_dsl), wgtfacq_c_dsl_ref.asnumpy())
+    assert test_utils.dallclose(data_alloc.as_numpy(wgtfacq_c_dsl), wgtfacq_c_ref.asnumpy())
diff --git a/model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp_dsl.py b/model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp.py
similarity index 90%
rename from model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp_dsl.py
rename to model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp.py
index 1c19a7ae09..f788c5497e 100644
--- a/model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp_dsl.py
+++ b/model/common/tests/common/metrics/unit_tests/test_compute_zdiff_gradp.py
@@ -14,10 +14,10 @@
 
 import icon4py.model.common.grid.horizontal as h_grid
 from icon4py.model.common import dimension as dims
-from icon4py.model.common.metrics.compute_zdiff_gradp_dsl import compute_zdiff_gradp_dsl
+from icon4py.model.common.metrics.compute_zdiff_gradp import compute_zdiff_gradp
 from icon4py.model.common.metrics.metric_fields import compute_flat_max_idx
 from icon4py.model.common.utils import data_allocation as data_alloc
-from icon4py.model.testing import definitions, test_utils
+from icon4py.model.testing import definitions, exchange_utils, test_utils
 from icon4py.model.testing.fixtures.datatest import (
     backend,
     data_provider,
@@ -30,8 +30,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
@@ -42,7 +40,7 @@
 
 @pytest.mark.level("unit")
 @pytest.mark.datatest
-def test_compute_zdiff_gradp_dsl(
+def test_compute_zdiff_gradp(
     icon_grid: base_grid.Grid,
     metrics_savepoint: sb.MetricSavepoint,
     interpolation_savepoint: sb.InterpolationSavepoint,
@@ -69,11 +67,11 @@ def test_compute_zdiff_gradp_dsl(
         c_lin_e=c_lin_e.ndarray,
         z_ifc=z_ifc.ndarray,
         k_lev=k_lev.ndarray,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=xp,
     )
 
-    zdiff_gradp_full_field, vertoffset_gradp_full_field = compute_zdiff_gradp_dsl(
+    zdiff_gradp_full_field, vertoffset_gradp_full_field = compute_zdiff_gradp(
         e2c=icon_grid.get_connectivity("E2C").ndarray,
         z_mc=z_mc.ndarray,
         c_lin_e=c_lin_e.ndarray,
@@ -83,7 +81,7 @@ def test_compute_zdiff_gradp_dsl(
         nlev=icon_grid.num_levels,
         horizontal_start=horizontal_start_edge,
         horizontal_start_1=start_nudging,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=xp,
     )
 
diff --git a/model/common/tests/common/metrics/unit_tests/test_metric_fields.py b/model/common/tests/common/metrics/unit_tests/test_metric_fields.py
index 9caf3fa81c..4b807246af 100644
--- a/model/common/tests/common/metrics/unit_tests/test_metric_fields.py
+++ b/model/common/tests/common/metrics/unit_tests/test_metric_fields.py
@@ -17,7 +17,7 @@
 from icon4py.model.common.grid import grid_refinement as refinement, horizontal
 from icon4py.model.common.metrics import metric_fields as mf
 from icon4py.model.common.utils import data_allocation as data_alloc
-from icon4py.model.testing import definitions, test_utils as testing_helpers
+from icon4py.model.testing import definitions, exchange_utils, test_utils as testing_helpers
 from icon4py.model.testing.definitions import construct_metrics_config
 from icon4py.model.testing.fixtures.datatest import (
     backend,
@@ -31,8 +31,6 @@
     processor_props,
 )
 
-from ... import utils
-
 
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
@@ -429,7 +427,7 @@ def test_compute_pressure_gradient_downward_extrapolation_mask_distance(
         c_lin_e=c_lin_e.ndarray,
         z_ifc=z_ifc.ndarray,
         k_lev=k.ndarray,
-        exchange=utils.dummy_exchange,
+        exchange=exchange_utils.dummy_exchange_with_bound_dim,
         array_ns=xp,
     )
     # TODO (nfarabullini): fix type ignore
diff --git a/model/common/tests/common/metrics/unit_tests/test_metrics_factory.py b/model/common/tests/common/metrics/unit_tests/test_metrics_factory.py
index 516f7d2d5a..7dc2bf014f 100644
--- a/model/common/tests/common/metrics/unit_tests/test_metrics_factory.py
+++ b/model/common/tests/common/metrics/unit_tests/test_metrics_factory.py
@@ -13,9 +13,11 @@
 if TYPE_CHECKING:
     import gt4py.next.typing as gtx_typing
 import pytest
+from numpy import testing as np_testing
 
+from icon4py.model.common import dimension as dims
 from icon4py.model.common.decomposition import definitions as decomposition
-from icon4py.model.common.grid import vertical as v_grid
+from icon4py.model.common.grid import horizontal as h_grid, vertical as v_grid
 from icon4py.model.common.interpolation import interpolation_attributes, interpolation_factory
 from icon4py.model.common.metrics import metrics_attributes as attrs, metrics_factory
 from icon4py.model.common.utils import data_allocation as data_alloc
@@ -501,11 +503,21 @@ def test_factory_zdiff_gradp(
     )
     field_1 = factory.get(attrs.ZDIFF_GRADP)
     field_2 = factory.get(attrs.VERTOFFSET_GRADP)
-    assert test_helpers.dallclose(
-        zdiff_gradp_ref.asnumpy(), field_1.asnumpy(), atol=1.0e-10, rtol=1e-9
+
+    # on the Fortran side, the vertidx_gradp is not initialized below start_lat_level2
+    start_lat_level2 = factory._grid.start_index(
+        h_grid.domain(dims.EdgeDim)(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_2)
+    )
+    test_helpers.assert_dallclose(
+        zdiff_gradp_ref.asnumpy()[start_lat_level2:, :, :],
+        field_1.asnumpy()[start_lat_level2:, :, :],
+        atol=1.0e-10,
+        rtol=1e-9,
     )
-    assert test_helpers.dallclose(
-        vertoffset_gradp_ref.asnumpy(), field_2.asnumpy(), atol=1.0e-10, rtol=1e-9
+
+    np_testing.assert_array_equal(
+        vertoffset_gradp_ref.asnumpy()[start_lat_level2:, :, :],
+        field_2.asnumpy()[start_lat_level2:, :, :],
     )
 
 
@@ -545,7 +557,7 @@ def test_factory_wgtfacq_e(
         topography_savepoint=topography_savepoint,
     )
     field = factory.get(attrs.WGTFACQ_E)
-    field_ref = metrics_savepoint.wgtfacq_e_dsl()
+    field_ref = metrics_savepoint.wgtfacq_e()
     # TODO: upgrade the dallclose such that it verifies the domain ranges.
     # This field is defined on k (nlev-3, nlev) an converting to numpy
     # doesn't know if it's there or (whatever-3, whatever)
diff --git a/model/common/tests/common/utils/__init__.py b/model/common/tests/common/utils/__init__.py
new file mode 100644
index 0000000000..de9850de36
--- /dev/null
+++ b/model/common/tests/common/utils/__init__.py
@@ -0,0 +1,7 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022-2024, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
diff --git a/model/common/tests/common/utils/unit_tests/__init__.py b/model/common/tests/common/utils/unit_tests/__init__.py
new file mode 100644
index 0000000000..de9850de36
--- /dev/null
+++ b/model/common/tests/common/utils/unit_tests/__init__.py
@@ -0,0 +1,7 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022-2024, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
diff --git a/model/common/tests/common/utils/unit_tests/test_field_utils.py b/model/common/tests/common/utils/unit_tests/test_field_utils.py
new file mode 100644
index 0000000000..30070a29db
--- /dev/null
+++ b/model/common/tests/common/utils/unit_tests/test_field_utils.py
@@ -0,0 +1,145 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022-2024, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
+
+import gt4py.next as gtx
+import numpy as np
+import pytest
+
+from icon4py.model.common import dimension as dims
+from icon4py.model.common.utils import field_utils
+
+
+@pytest.fixture
+def allocator():
+    return np
+
+
+class TestFlip:
+    def test_flip_1d(self, allocator):
+        field = gtx.as_field({dims.CellDim: range(0, 4)}, np.array([10.0, 20.0, 30.0, 40.0]))
+        result = field_utils.flip(field, dims.CellDim, allocator=allocator)
+
+        assert np.array_equal(result.ndarray, [40.0, 30.0, 20.0, 10.0])
+        assert result.domain == field.domain
+
+    def test_flip_1d_nonzero_start(self, allocator):
+        field = gtx.as_field({dims.CellDim: range(3, 7)}, np.array([10.0, 20.0, 30.0, 40.0]))
+        result = field_utils.flip(field, dims.CellDim, allocator=allocator)
+
+        assert np.array_equal(result.ndarray, [40.0, 30.0, 20.0, 10.0])
+        assert result.domain == field.domain
+
+    def test_flip_2d_along_first_dim(self, allocator):
+        data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
+        field = gtx.as_field({dims.CellDim: range(0, 2), dims.KDim: range(0, 3)}, data)
+        result = field_utils.flip(field, dims.CellDim, allocator=allocator)
+
+        expected = np.array([[4.0, 5.0, 6.0], [1.0, 2.0, 3.0]])
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == field.domain
+
+    def test_flip_2d_along_second_dim(self, allocator):
+        data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
+        field = gtx.as_field({dims.CellDim: range(0, 2), dims.KDim: range(0, 3)}, data)
+        result = field_utils.flip(field, dims.KDim, allocator=allocator)
+
+        expected = np.array([[3.0, 2.0, 1.0], [6.0, 5.0, 4.0]])
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == field.domain
+
+    def test_flip_2d_nonzero_start_along_first_dim(self, allocator):
+        data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
+        field = gtx.as_field({dims.CellDim: range(5, 7), dims.KDim: range(2, 5)}, data)
+        result = field_utils.flip(field, dims.CellDim, allocator=allocator)
+
+        expected = np.array([[4.0, 5.0, 6.0], [1.0, 2.0, 3.0]])
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == field.domain
+
+    def test_flip_2d_nonzero_start_along_second_dim(self, allocator):
+        data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
+        field = gtx.as_field({dims.CellDim: range(5, 7), dims.KDim: range(2, 5)}, data)
+        result = field_utils.flip(field, dims.KDim, allocator=allocator)
+
+        expected = np.array([[3.0, 2.0, 1.0], [6.0, 5.0, 4.0]])
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == field.domain
+
+    def test_flip_preserves_dtype(self, allocator):
+        field = gtx.as_field({dims.KDim: range(3, 6)}, np.array([1, 2, 3], dtype=np.int32))
+        result = field_utils.flip(field, dims.KDim, allocator=allocator)
+
+        assert result.ndarray.dtype == np.int32
+        assert np.array_equal(result.ndarray, [3, 2, 1])
+
+
+class TestIndex2Offset:
+    def test_1d_zero_start(self, allocator):
+        # indices: [2, 1, 3, 0], positions: [0, 1, 2, 3]
+        # offsets: [2-0, 1-1, 3-2, 0-3] = [2, 0, 1, -3]  # noqa: ERA001
+        index_field = gtx.as_field(
+            {dims.CellDim: range(0, 4)}, np.array([2, 1, 3, 0], dtype=np.int32)
+        )
+        result = field_utils.index2offset(index_field, dims.CellDim, allocator=allocator)
+
+        assert np.array_equal(result.ndarray, [2, 0, 1, -3])
+        assert result.domain == index_field.domain
+
+    def test_1d_nonzero_start(self, allocator):
+        # indices: [5, 3, 4, 2], positions: [2, 3, 4, 5]
+        # offsets: [5-2, 3-3, 4-4, 2-5] = [3, 0, 0, -3]  # noqa: ERA001
+        index_field = gtx.as_field(
+            {dims.CellDim: range(2, 6)}, np.array([5, 3, 4, 2], dtype=np.int32)
+        )
+        result = field_utils.index2offset(index_field, dims.CellDim, allocator=allocator)
+
+        assert np.array_equal(result.ndarray, [3, 0, 0, -3])
+        assert result.domain == index_field.domain
+
+    def test_1d_identity_permutation(self, allocator):
+        # indices == positions => all offsets are 0
+        index_field = gtx.as_field({dims.KDim: range(3, 7)}, np.array([3, 4, 5, 6], dtype=np.int32))
+        result = field_utils.index2offset(index_field, dims.KDim, allocator=allocator)
+
+        assert np.array_equal(result.ndarray, [0, 0, 0, 0])
+
+    def test_2d_along_second_dim(self, allocator):
+        # shape (3, 4), apply along KDim (axis=1)
+        # positions for KDim with range(0, 4): [0, 1, 2, 3]
+        data = np.array([[0, 3, 1, 2], [1, 0, 2, 3], [2, 1, 3, 0]], dtype=np.int32)
+        index_field = gtx.as_field({dims.CellDim: range(0, 3), dims.KDim: range(0, 4)}, data)
+        result = field_utils.index2offset(index_field, dims.KDim, allocator=allocator)
+
+        expected = np.array([[0, 2, -1, -1], [1, -1, 0, 0], [2, 0, 1, -3]], dtype=np.int32)
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == index_field.domain
+
+    def test_2d_nonzero_start_along_second_dim(self, allocator):
+        # shape (2, 3), KDim starts at 5
+        # positions for KDim with range(5, 8): [5, 6, 7]
+        data = np.array([[7, 5, 6], [6, 7, 5]], dtype=np.int32)
+        index_field = gtx.as_field({dims.CellDim: range(2, 4), dims.KDim: range(5, 8)}, data)
+        result = field_utils.index2offset(index_field, dims.KDim, allocator=allocator)
+
+        # offsets: [[7-5, 5-6, 6-7], [6-5, 7-6, 5-7]] = [[2, -1, -1], [1, 1, -2]]  # noqa: ERA001
+        expected = np.array([[2, -1, -1], [1, 1, -2]], dtype=np.int32)
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == index_field.domain
+
+    def test_2d_along_first_dim_nonzero_start(self, allocator):
+        # shape (3, 2), apply along CellDim (axis=0)
+        # positions for CellDim with range(10, 13): [10, 11, 12]
+        # arange is 1D along CellDim axis, broadcast subtracted from (3, 2) data
+        data = np.array([[12, 11], [10, 12], [11, 10]], dtype=np.int32)
+        index_field = gtx.as_field({dims.CellDim: range(10, 13), dims.KDim: range(0, 2)}, data)
+        result = field_utils.index2offset(index_field, dims.CellDim, allocator=allocator)
+
+        # offsets: [[12-10, 11-10], [10-11, 12-11], [11-12, 10-12]] = [[2, 1], [-1, 1], [-1, -2]]  # noqa: ERA001
+        expected = np.array([[2, 1], [-1, 1], [-1, -2]], dtype=np.int32)
+        assert np.array_equal(result.ndarray, expected)
+        assert result.domain == index_field.domain
diff --git a/model/driver/src/icon4py/model/driver/initialization_utils.py b/model/driver/src/icon4py/model/driver/initialization_utils.py
index a78be8f6c3..9188fad450 100644
--- a/model/driver/src/icon4py/model/driver/initialization_utils.py
+++ b/model/driver/src/icon4py/model/driver/initialization_utils.py
@@ -491,7 +491,7 @@ def read_static_fields(
             time_extrapolation_parameter_for_exner=metrics_savepoint.exner_exfac(),
             reference_exner_at_cells_on_model_levels=metrics_savepoint.exner_ref_mc(),
             wgtfac_c=metrics_savepoint.wgtfac_c(),
-            wgtfacq_c=metrics_savepoint.wgtfacq_c_dsl(),
+            wgtfacq_c=metrics_savepoint.wgtfacq_c(),
             inv_ddqz_z_full=metrics_savepoint.inv_ddqz_z_full(),
             reference_rho_at_cells_on_model_levels=metrics_savepoint.rho_ref_mc(),
             reference_theta_at_cells_on_model_levels=metrics_savepoint.theta_ref_mc(),
@@ -511,7 +511,7 @@ def read_static_fields(
             ddqz_z_full_e=metrics_savepoint.ddqz_z_full_e(),
             ddxt_z_full=metrics_savepoint.ddxt_z_full(),
             wgtfac_e=metrics_savepoint.wgtfac_e(),
-            wgtfacq_e=metrics_savepoint.wgtfacq_e_dsl(),
+            wgtfacq_e=metrics_savepoint.wgtfacq_e(),
             exner_w_implicit_weight_parameter=metrics_savepoint.vwind_impl_wgt(),
             horizontal_mask_for_3d_divdamp=metrics_savepoint.hmask_dd3d(),
             scaling_factor_for_3d_divdamp=metrics_savepoint.scalfac_dd3d(),
diff --git a/model/driver/tests/driver/integration_tests/test_icon4py.py b/model/driver/tests/driver/integration_tests/test_icon4py.py
index d5c6572884..1d686c154b 100644
--- a/model/driver/tests/driver/integration_tests/test_icon4py.py
+++ b/model/driver/tests/driver/integration_tests/test_icon4py.py
@@ -227,7 +227,7 @@ def test_run_timeloop_single_step(
         time_extrapolation_parameter_for_exner=metrics_savepoint.exner_exfac(),
         reference_exner_at_cells_on_model_levels=metrics_savepoint.exner_ref_mc(),
         wgtfac_c=metrics_savepoint.wgtfac_c(),
-        wgtfacq_c=metrics_savepoint.wgtfacq_c_dsl(),
+        wgtfacq_c=metrics_savepoint.wgtfacq_c(),
         inv_ddqz_z_full=metrics_savepoint.inv_ddqz_z_full(),
         reference_rho_at_cells_on_model_levels=metrics_savepoint.rho_ref_mc(),
         reference_theta_at_cells_on_model_levels=metrics_savepoint.theta_ref_mc(),
@@ -247,7 +247,7 @@ def test_run_timeloop_single_step(
         ddqz_z_full_e=metrics_savepoint.ddqz_z_full_e(),
         ddxt_z_full=metrics_savepoint.ddxt_z_full(),
         wgtfac_e=metrics_savepoint.wgtfac_e(),
-        wgtfacq_e=metrics_savepoint.wgtfacq_e_dsl(),
+        wgtfacq_e=metrics_savepoint.wgtfacq_e(),
         exner_w_implicit_weight_parameter=metrics_savepoint.vwind_impl_wgt(),
         horizontal_mask_for_3d_divdamp=metrics_savepoint.hmask_dd3d(),
         scaling_factor_for_3d_divdamp=metrics_savepoint.scalfac_dd3d(),
diff --git a/model/testing/src/icon4py/model/testing/definitions.py b/model/testing/src/icon4py/model/testing/definitions.py
index 3e63e40fa9..41f2f0ff76 100644
--- a/model/testing/src/icon4py/model/testing/definitions.py
+++ b/model/testing/src/icon4py/model/testing/definitions.py
@@ -226,7 +226,7 @@ class Experiment:
     description: str
     grid: GridDescription
     num_levels: int
-    version: int = 2
+    version: int = 3
 
 
 class Experiments:
diff --git a/model/common/tests/common/utils.py b/model/testing/src/icon4py/model/testing/exchange_utils.py
similarity index 78%
rename from model/common/tests/common/utils.py
rename to model/testing/src/icon4py/model/testing/exchange_utils.py
index 602e187dc4..1a43c4bf39 100644
--- a/model/common/tests/common/utils.py
+++ b/model/testing/src/icon4py/model/testing/exchange_utils.py
@@ -11,5 +11,5 @@
 from icon4py.model.common.utils import data_allocation as data_alloc
 
 
-def dummy_exchange(*field: data_alloc.NDArray, stream: Any = None) -> None:
+def dummy_exchange_with_bound_dim(*field: data_alloc.NDArray, stream: Any = None) -> None:
     return None
diff --git a/model/testing/src/icon4py/model/testing/serialbox.py b/model/testing/src/icon4py/model/testing/serialbox.py
index d015b62e58..15e86d11e3 100644
--- a/model/testing/src/icon4py/model/testing/serialbox.py
+++ b/model/testing/src/icon4py/model/testing/serialbox.py
@@ -7,6 +7,7 @@
 # SPDX-License-Identifier: BSD-3-Clause
 import functools
 import logging
+from collections.abc import Sequence
 from typing import Final, Literal, TypeAlias
 
 import gt4py.next as gtx
@@ -20,7 +21,7 @@
 from icon4py.model.common import dimension as dims, model_backends, type_alias
 from icon4py.model.common.grid import base, horizontal as h_grid, icon, utils as grid_utils
 from icon4py.model.common.states import prognostic_state
-from icon4py.model.common.utils import data_allocation as data_alloc
+from icon4py.model.common.utils import data_allocation as data_alloc, field_utils
 
 
 log = logging.getLogger(__name__)
@@ -84,8 +85,20 @@ def wrapper(self, *args, **kwargs):
     def log_meta_info(self):
         self.log.info(self.savepoint.metainfo)
 
-    def _get_field(self, name, *dimensions, dtype=float):
+    def _get_field(
+        self,
+        name,
+        *dimensions,
+        dtype=float,
+        slice_: int | slice | tuple[int | slice, ...] | None = None,
+        transpose: None | Sequence[int] = None,
+    ):
+        # Note: slice is applied before transpose!
         buffer = np.squeeze(self.serializer.read(name, self.savepoint).astype(dtype))
+        if slice_ is not None:
+            buffer = buffer[slice_]
+        if transpose is not None:
+            buffer = np.transpose(buffer, axes=transpose)
         buffer = self._reduce_to_dim_size(buffer, dimensions)
 
         self.log.debug(f"{name} {buffer.shape}")
@@ -673,31 +686,33 @@ def pos_on_tplane_e_y(self):
         return self._get_field("pos_on_tplane_e_y", dims.EdgeDim, dims.E2CDim)[:, 0:2]
 
     def rbf_vec_coeff_e(self):
-        return self._get_field("rbf_vec_coeff_e", dims.EdgeDim, dims.E2C2EDim)
+        return self._get_field("rbf_vec_coeff_e", dims.EdgeDim, dims.E2C2EDim, transpose=(1, 0))
 
     @IconSavepoint.optionally_registered()
     def rbf_vec_coeff_c1(self):
-        dimensions = (dims.CellDim, dims.C2E2C2EDim)
-        buffer = np.squeeze(
-            self.serializer.read("rbf_vec_coeff_c1", self.savepoint).astype(float)
-        ).transpose()
-        buffer = self._reduce_to_dim_size(buffer, dimensions)
-        return gtx.as_field(dimensions, buffer, allocator=self.backend)
+        return self._get_field("rbf_vec_coeff_c1", dims.CellDim, dims.C2E2C2EDim, transpose=(1, 0))
 
     @IconSavepoint.optionally_registered()
     def rbf_vec_coeff_c2(self):
-        dimensions = (dims.CellDim, dims.C2E2C2EDim)
-        buffer = np.squeeze(
-            self.serializer.read("rbf_vec_coeff_c2", self.savepoint).astype(float)
-        ).transpose()
-        buffer = self._reduce_to_dim_size(buffer, dimensions)
-        return gtx.as_field(dimensions, buffer, allocator=self.backend)
+        return self._get_field("rbf_vec_coeff_c2", dims.CellDim, dims.C2E2C2EDim, transpose=(1, 0))
 
     def rbf_vec_coeff_v1(self):
-        return self._get_field("rbf_vec_coeff_v1", dims.VertexDim, dims.V2EDim)
+        return self._get_field(
+            "rbf_vec_coeff_v",
+            dims.VertexDim,
+            dims.V2EDim,
+            slice_=(slice(None), 0, slice(None)),
+            transpose=(1, 0),
+        )
 
     def rbf_vec_coeff_v2(self):
-        return self._get_field("rbf_vec_coeff_v2", dims.VertexDim, dims.V2EDim)
+        return self._get_field(
+            "rbf_vec_coeff_v",
+            dims.VertexDim,
+            dims.V2EDim,
+            slice_=(slice(None), 1, slice(None)),
+            transpose=(1, 0),
+        )
 
     def rbf_vec_idx_v(self):
         return self._get_field("rbf_vec_idx_v", dims.VertexDim, dims.V2EDim)
@@ -806,34 +821,27 @@ def vwind_impl_wgt(self):
         return self._get_field("vwind_impl_wgt", dims.CellDim)
 
     def wgtfacq_c(self):
-        return self._get_field("wgtfacq_c", dims.CellDim, dims.KDim)
-
-    def wgtfacq_c_dsl(self):
-        ar = self.wgtfacq_c().ndarray
-        k = ar.shape[1]
-        wgtfac_c = self.wgtfac_c()
-        cell_range = wgtfac_c.domain[dims.CellDim].unit_range
-        nlev = wgtfac_c.domain[dims.KDim].unit_range.stop - 1
-        k_range = (nlev - k, nlev)
-        cell_kflip_domain = gtx.domain(
-            {
-                dims.CellDim: cell_range,
-                dims.KDim: k_range,
-            }
-        )
-        return data_alloc.kflip_wgtfacq(
-            arr=ar,
-            domain=cell_kflip_domain,
+        # The Fortran array stores the surface levels in reversed order.
+        wgtfacq_c_fortran = self._get_field("wgtfacq_c", dims.CellDim, dims.KDim)
+        assert len(wgtfacq_c_fortran.domain[dims.KDim].unit_range) == 3
+        nlev = self.sizes[dims.KDim]
+        return field_utils.flip(
+            wgtfacq_c_fortran(dims.KDim - (nlev - 3)),  # GT4Py embedded shift
+            dims.KDim,
             allocator=model_backends.get_allocator(self.backend),
         )
 
     def zdiff_gradp(self):
-        return self._get_field("zdiff_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim)
+        return self._get_field("zdiff_gradp", dims.EdgeDim, dims.E2CDim, dims.KDim)
 
     def vertoffset_gradp(self):
-        return self._get_field(
-            "vertoffset_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32
+        # In Fortran `vertidx_gradp` contains `0`s in areas where the array is not used.
+        # When we translate to offsets we just subtract the current index, therefore these values will be negative.
+        # Since in Fortran accessing index `0` would be out-of-bounds, we should be safe.
+        vertidx_gradp = data_alloc.adjust_fortran_indices(
+            self._get_field("vertidx_gradp", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32)
         )
+        return field_utils.index2offset(vertidx_gradp, dims.KDim, self.backend)
 
     def coeff1_dwdz(self):
         return self._get_field("coeff1_dwdz", dims.CellDim, dims.KDim)
@@ -866,24 +874,13 @@ def wgtfac_e(self):
         return self._get_field("wgtfac_e", dims.EdgeDim, dims.KDim)
 
     def wgtfacq_e(self):
-        return self._get_field("wgtfacq_e", dims.EdgeDim, dims.KDim)
-
-    def wgtfacq_e_dsl(self):
-        ar = self.wgtfacq_e().ndarray
-        k = ar.shape[1]
-        wgtfac_e = self.wgtfac_e()
-        edge_range = wgtfac_e.domain[dims.EdgeDim].unit_range
-        nlev = wgtfac_e.domain[dims.KDim].unit_range.stop - 1
-        k_range = (nlev - k, nlev)
-        edge_kflip_domain = gtx.domain(
-            {
-                dims.EdgeDim: edge_range,
-                dims.KDim: k_range,
-            }
-        )
-        return data_alloc.kflip_wgtfacq(
-            arr=ar,
-            domain=edge_kflip_domain,
+        # The Fortran array stores the surface levels in reversed order.
+        wgtfacq_e_fortran = self._get_field("wgtfacq_e", dims.EdgeDim, dims.KDim)
+        assert len(wgtfacq_e_fortran.domain[dims.KDim].unit_range) == 3
+        nlev = self.sizes[dims.KDim]
+        return field_utils.flip(
+            wgtfacq_e_fortran(dims.KDim - (nlev - 3)),  # GT4Py embedded shift
+            dims.KDim,
             allocator=model_backends.get_allocator(self.backend),
         )
 
diff --git a/model/testing/src/icon4py/model/testing/stencil_tests.py b/model/testing/src/icon4py/model/testing/stencil_tests.py
index d19a11ebdf..444c0a8f4e 100644
--- a/model/testing/src/icon4py/model/testing/stencil_tests.py
+++ b/model/testing/src/icon4py/model/testing/stencil_tests.py
@@ -31,6 +31,7 @@
 from icon4py.model.common import model_backends, model_options
 from icon4py.model.common.grid import base
 from icon4py.model.common.utils import device_utils
+from icon4py.model.testing import test_utils
 
 
 def allocate_data(
@@ -266,7 +267,7 @@ def _verify_stencil_test(
             relative_tolerance = 3e-6
             if isinstance(input_data_name, tuple):
                 for i_out_field, out_field in enumerate(input_data_name):
-                    np.testing.assert_allclose(
+                    test_utils.assert_dallclose(
                         out_field.asnumpy()[gtslice],
                         reference_outputs[name][i_out_field][refslice],
                         equal_nan=True,
@@ -276,7 +277,7 @@ def _verify_stencil_test(
             else:
                 reference_outputs_name = reference_outputs[name]  # for mypy
                 assert isinstance(reference_outputs_name, np.ndarray)
-                np.testing.assert_allclose(
+                test_utils.assert_dallclose(
                     input_data_name.asnumpy()[gtslice],
                     reference_outputs_name[refslice],
                     equal_nan=True,
diff --git a/model/testing/src/icon4py/model/testing/test_utils.py b/model/testing/src/icon4py/model/testing/test_utils.py
index 40cceb2462..312356861a 100644
--- a/model/testing/src/icon4py/model/testing/test_utils.py
+++ b/model/testing/src/icon4py/model/testing/test_utils.py
@@ -11,6 +11,7 @@
 
 import gt4py.next.typing as gtx_typing
 import numpy as np
+import numpy.testing as np_testing
 import numpy.typing as npt
 import pytest
 from typing_extensions import Buffer
@@ -25,9 +26,35 @@ def dallclose(
     atol: float = 0.0,
     equal_nan: bool = False,
 ) -> bool:
+    """
+    'numpy.allclose', but with double precision default tolerances.
+    """
     return np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan)
 
 
+def assert_dallclose(
+    actual: npt.ArrayLike,
+    desired: npt.ArrayLike,
+    rtol: float = 1.0e-12,
+    atol: float = 0.0,
+    equal_nan: bool = False,
+    err_msg: str = "",
+    verbose: bool = True,
+) -> None:
+    """
+    'numpy.testing.assert_allclose', but with double precision default tolerances.
+    """
+    np_testing.assert_allclose(
+        actual,  # type: ignore[arg-type]
+        desired,  # type: ignore[arg-type]
+        rtol=rtol,
+        atol=atol,
+        equal_nan=equal_nan,
+        err_msg=err_msg,
+        verbose=verbose,
+    )
+
+
 def is_sorted(array: np.ndarray) -> bool:
     return bool((array[:-1] <= array[1:]).all())
 
diff --git a/tools/src/icon4py/tools/py2fgen/wrappers/common.py b/tools/src/icon4py/tools/py2fgen/wrappers/common.py
index 79785b991e..f6d908926c 100644
--- a/tools/src/icon4py/tools/py2fgen/wrappers/common.py
+++ b/tools/src/icon4py/tools/py2fgen/wrappers/common.py
@@ -97,6 +97,16 @@
     ),
 ]
 
+Float64Array3D: TypeAlias = Annotated[
+    data_alloc.NDArray,
+    py2fgen.ArrayParamDescriptor(
+        rank=3,
+        dtype=ts.ScalarKind.FLOAT64,
+        memory_space=py2fgen.MemorySpace.MAYBE_DEVICE,
+        is_optional=False,
+    ),
+]
+
 OptionalFloat64Array1D: TypeAlias = Annotated[
     data_alloc.NDArray | None,
     py2fgen.ArrayParamDescriptor(
diff --git a/tools/src/icon4py/tools/py2fgen/wrappers/diffusion_wrapper.py b/tools/src/icon4py/tools/py2fgen/wrappers/diffusion_wrapper.py
index 0c8044a1c3..ecfe925b2d 100644
--- a/tools/src/icon4py/tools/py2fgen/wrappers/diffusion_wrapper.py
+++ b/tools/src/icon4py/tools/py2fgen/wrappers/diffusion_wrapper.py
@@ -63,8 +63,7 @@ def diffusion_init(
     geofac_grg_y: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
     geofac_n2s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
     nudgecoeff_e: fa.EdgeField[wpfloat],
-    rbf_coeff_1: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
-    rbf_coeff_2: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
+    rbf_vec_coeff_v: wrapper_common.Float64Array3D,
     zd_cellidx: wrapper_common.OptionalInt32Array2D,
     zd_vertidx: wrapper_common.OptionalInt32Array2D,
     zd_intcoef: wrapper_common.OptionalFloat64Array2D,
@@ -92,7 +91,8 @@ def diffusion_init(
             "Need to initialise grid using 'grid_init' before running 'diffusion_init'."
         )
 
-    on_gpu = theta_ref_mc.array_ns != np  # TODO(havogt): expose `on_gpu` from py2fgen
+    xp = theta_ref_mc.array_ns
+    on_gpu = xp != np  # TODO(havogt): expose `on_gpu` from py2fgen
     actual_backend = wrapper_common.select_backend(
         wrapper_common.BackendIntEnum(backend), on_gpu=on_gpu
     )
@@ -134,7 +134,6 @@ def diffusion_init(
             dims.KDim: nlev,
         }
     )
-    xp = wgtfac_c.array_ns
 
     if zd_cellidx is None:
         # then zdiffu_t is False or the list on that rank is empty, then all of the following are not initialized
@@ -194,6 +193,15 @@ def diffusion_init(
         zd_diffcoef=zd_diffcoef,
     )
 
+    # Create separate fields for the two components of the RBF vector coefficients and swap.
+    # TODO(havogt): we could use GT4Py's named collections.
+    rbf_coeff_1 = gtx.as_field(
+        [dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 0, :]), allocator=allocator
+    )
+    rbf_coeff_2 = gtx.as_field(
+        [dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 1, :]), allocator=allocator
+    )
+
     # Interpolation state
     interpolation_state = DiffusionInterpolationState(
         e_bln_c_s=e_bln_c_s,
diff --git a/tools/src/icon4py/tools/py2fgen/wrappers/dycore_wrapper.py b/tools/src/icon4py/tools/py2fgen/wrappers/dycore_wrapper.py
index 35cfdb4843..142d00de4c 100644
--- a/tools/src/icon4py/tools/py2fgen/wrappers/dycore_wrapper.py
+++ b/tools/src/icon4py/tools/py2fgen/wrappers/dycore_wrapper.py
@@ -29,7 +29,7 @@
 from icon4py.model.atmosphere.dycore import dycore_states, solve_nonhydro
 from icon4py.model.common import dimension as dims, model_backends, utils as common_utils
 from icon4py.model.common.states.prognostic_state import PrognosticState
-from icon4py.model.common.utils import data_allocation as data_alloc
+from icon4py.model.common.utils import data_allocation as data_alloc, field_utils
 from icon4py.tools import py2fgen
 from icon4py.tools.common.logger import setup_logger
 from icon4py.tools.py2fgen.wrappers import common as wrapper_common, grid_wrapper, icon4py_export
@@ -56,10 +56,9 @@ def solve_nh_init(
     geofac_rot: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
     pos_on_tplane_e_1: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim], gtx.float64],
     pos_on_tplane_e_2: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim], gtx.float64],
-    rbf_vec_coeff_e: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2C2EDim], gtx.float64],
+    rbf_vec_coeff_e: wrapper_common.Float64Array2D,
     e_bln_c_s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2EDim], gtx.float64],
-    rbf_coeff_1: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
-    rbf_coeff_2: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
+    rbf_vec_coeff_v: wrapper_common.Float64Array3D,
     geofac_div: gtx.Field[gtx.Dims[dims.CellDim, dims.C2EDim], gtx.float64],
     geofac_n2s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
     geofac_grg_x: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
@@ -84,7 +83,7 @@ def solve_nh_init(
     theta_ref_me: gtx.Field[gtx.Dims[dims.EdgeDim, dims.KDim], gtx.float64],
     ddxn_z_full: gtx.Field[gtx.Dims[dims.EdgeDim, dims.KDim], gtx.float64],
     zdiff_gradp: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim, dims.KDim], gtx.float64],
-    vertoffset_gradp: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim, dims.KDim], gtx.int32],
+    vertidx_gradp: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim, dims.KDim], gtx.int32],
     pg_edgeidx: wrapper_common.OptionalInt32Array1D,
     pg_vertidx: wrapper_common.OptionalInt32Array1D,
     pg_exdist: wrapper_common.OptionalFloat64Array1D,
@@ -129,7 +128,8 @@ def solve_nh_init(
     if grid_wrapper.grid_state is None:
         raise Exception("Need to initialise grid using 'grid_init' before running 'solve_nh_init'.")
 
-    on_gpu = c_lin_e.array_ns != np  # TODO(havogt): expose `on_gpu` from py2fgen
+    xp = c_lin_e.array_ns
+    on_gpu = xp != np  # TODO(havogt): expose `on_gpu` from py2fgen
     actual_backend = wrapper_common.select_backend(
         wrapper_common.BackendIntEnum(backend), on_gpu=on_gpu
     )
@@ -181,6 +181,19 @@ def solve_nh_init(
     )
     nonhydro_params = solve_nonhydro.NonHydrostaticParams(config)
 
+    # Create separate fields for the two components of the RBF vector coefficients and swap.
+    # TODO(havogt): we could use GT4Py's named collections.
+    rbf_coeff_1 = gtx.as_field(
+        [dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 0, :]), allocator=allocator
+    )
+    rbf_coeff_2 = gtx.as_field(
+        [dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 1, :]), allocator=allocator
+    )
+
+    # Swap indices in rbf_vec_coeff_e. TODO(havogt): Should eventually be done on the Fortran side.
+    rbf_vec_coeff_e_transposed = gtx.as_field(
+        [dims.EdgeDim, dims.E2C2EDim], xp.transpose(rbf_vec_coeff_e), allocator=allocator
+    )
     interpolation_state = dycore_states.InterpolationState(
         c_lin_e=c_lin_e,
         c_intp=c_intp,
@@ -189,7 +202,7 @@ def solve_nh_init(
         geofac_rot=geofac_rot,
         pos_on_tplane_e_1=pos_on_tplane_e_1[:, 0:2],
         pos_on_tplane_e_2=pos_on_tplane_e_2[:, 0:2],
-        rbf_vec_coeff_e=rbf_vec_coeff_e,
+        rbf_vec_coeff_e=rbf_vec_coeff_e_transposed,
         e_bln_c_s=e_bln_c_s,
         rbf_coeff_1=rbf_coeff_1,
         rbf_coeff_2=rbf_coeff_2,
@@ -201,25 +214,25 @@ def solve_nh_init(
     )
 
     nlev = wgtfac_c.domain[dims.KDim].unit_range.stop - 1
-    k = wgtfacq_c.ndarray.shape[1]
-    cell_kflip_domain = gtx.domain(
-        {
-            dims.CellDim: wgtfac_c.domain[dims.CellDim].unit_range,
-            dims.KDim: (nlev - k, nlev),
-        }
-    )
-    k = wgtfacq_e.ndarray.shape[1]
-    edge_kflip_domain = gtx.domain(
-        {
-            dims.EdgeDim: wgtfac_e.domain[dims.EdgeDim].unit_range,
-            dims.KDim: (nlev - k, nlev),
-        }
-    )
-    wgtfacq_c = data_alloc.kflip_wgtfacq(
-        arr=wgtfacq_c.ndarray, domain=cell_kflip_domain, allocator=allocator
-    )
-    wgtfacq_e = data_alloc.kflip_wgtfacq(
-        arr=wgtfacq_e.ndarray, domain=edge_kflip_domain, allocator=allocator
+    if len(wgtfacq_c.domain[dims.KDim].unit_range) != 3:
+        raise ValueError(
+            f"Expected wgtfacq_c to have a vertical dimension of size 3, but got {len(wgtfacq_c.domain[dims.KDim].unit_range)}."
+        )
+    # uses GT4Py's embedded shift to move the domain to surface levels
+    wgtfacq_c = field_utils.flip(wgtfacq_c(dims.KDim - (nlev - 3)), dims.KDim, allocator=allocator)
+
+    if len(wgtfacq_e.domain[dims.KDim].unit_range) != 3:
+        raise ValueError(
+            f"Expected wgtfacq_e to have a vertical dimension of size 3, but got {len(wgtfacq_e.domain[dims.KDim].unit_range)}."
+        )
+    # uses GT4Py's embedded shift to move the domain to surface levels
+    wgtfacq_e = field_utils.flip(wgtfacq_e(dims.KDim - (nlev - 3)), dims.KDim, allocator=allocator)
+
+    # In Fortran `vertidx_gradp` contains `0`s in areas where the array is not used.
+    # When we translate to offsets we just subtract the current index, therefore these values will be negative.
+    # Since in Fortran accessing index `0` would be out-of-bounds, we should be safe.
+    vertoffset_gradp = field_utils.index2offset(
+        data_alloc.adjust_fortran_indices(vertidx_gradp), dims.KDim, allocator
     )
 
     metric_state_nonhydro = dycore_states.MetricStateNonHydro(
diff --git a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.f90 b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.f90
index 9f306faa0c..e7e4dde5ce 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.f90
+++ b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.f90
@@ -126,12 +126,10 @@ function diffusion_init_wrapper(theta_ref_mc, &
                                       geofac_n2s_size_1, &
                                       nudgecoeff_e, &
                                       nudgecoeff_e_size_0, &
-                                      rbf_coeff_1, &
-                                      rbf_coeff_1_size_0, &
-                                      rbf_coeff_1_size_1, &
-                                      rbf_coeff_2, &
-                                      rbf_coeff_2_size_0, &
-                                      rbf_coeff_2_size_1, &
+                                      rbf_vec_coeff_v, &
+                                      rbf_vec_coeff_v_size_0, &
+                                      rbf_vec_coeff_v_size_1, &
+                                      rbf_vec_coeff_v_size_2, &
                                       zd_cellidx, &
                                       zd_cellidx_size_0, &
                                       zd_cellidx_size_1, &
@@ -210,17 +208,13 @@ function diffusion_init_wrapper(theta_ref_mc, &
 
          integer(c_int), value :: nudgecoeff_e_size_0
 
-         type(c_ptr), value, target :: rbf_coeff_1
+         type(c_ptr), value, target :: rbf_vec_coeff_v
 
-         integer(c_int), value :: rbf_coeff_1_size_0
+         integer(c_int), value :: rbf_vec_coeff_v_size_0
 
-         integer(c_int), value :: rbf_coeff_1_size_1
+         integer(c_int), value :: rbf_vec_coeff_v_size_1
 
-         type(c_ptr), value, target :: rbf_coeff_2
-
-         integer(c_int), value :: rbf_coeff_2_size_0
-
-         integer(c_int), value :: rbf_coeff_2_size_1
+         integer(c_int), value :: rbf_vec_coeff_v_size_2
 
          type(c_ptr), value, target :: zd_cellidx
 
@@ -483,8 +477,7 @@ subroutine diffusion_init(theta_ref_mc, &
                              geofac_grg_y, &
                              geofac_n2s, &
                              nudgecoeff_e, &
-                             rbf_coeff_1, &
-                             rbf_coeff_2, &
+                             rbf_vec_coeff_v, &
                              zd_cellidx, &
                              zd_vertidx, &
                              zd_intcoef, &
@@ -525,9 +518,7 @@ subroutine diffusion_init(theta_ref_mc, &
 
       real(c_double), dimension(:), target :: nudgecoeff_e
 
-      real(c_double), dimension(:, :), target :: rbf_coeff_1
-
-      real(c_double), dimension(:, :), target :: rbf_coeff_2
+      real(c_double), dimension(:, :, :), target :: rbf_vec_coeff_v
 
       integer(c_int), dimension(:, :), pointer :: zd_cellidx
 
@@ -603,13 +594,11 @@ subroutine diffusion_init(theta_ref_mc, &
 
       integer(c_int) :: nudgecoeff_e_size_0
 
-      integer(c_int) :: rbf_coeff_1_size_0
+      integer(c_int) :: rbf_vec_coeff_v_size_0
 
-      integer(c_int) :: rbf_coeff_1_size_1
+      integer(c_int) :: rbf_vec_coeff_v_size_1
 
-      integer(c_int) :: rbf_coeff_2_size_0
-
-      integer(c_int) :: rbf_coeff_2_size_1
+      integer(c_int) :: rbf_vec_coeff_v_size_2
 
       integer(c_int) :: zd_cellidx_size_0
 
@@ -652,8 +641,7 @@ subroutine diffusion_init(theta_ref_mc, &
       !$acc host_data use_device(geofac_grg_y)
       !$acc host_data use_device(geofac_n2s)
       !$acc host_data use_device(nudgecoeff_e)
-      !$acc host_data use_device(rbf_coeff_1)
-      !$acc host_data use_device(rbf_coeff_2)
+      !$acc host_data use_device(rbf_vec_coeff_v)
       !$acc host_data use_device(zd_cellidx) if(associated(zd_cellidx))
       !$acc host_data use_device(zd_vertidx) if(associated(zd_vertidx))
       !$acc host_data use_device(zd_intcoef) if(associated(zd_intcoef))
@@ -688,11 +676,9 @@ subroutine diffusion_init(theta_ref_mc, &
 
       nudgecoeff_e_size_0 = SIZE(nudgecoeff_e, 1)
 
-      rbf_coeff_1_size_0 = SIZE(rbf_coeff_1, 1)
-      rbf_coeff_1_size_1 = SIZE(rbf_coeff_1, 2)
-
-      rbf_coeff_2_size_0 = SIZE(rbf_coeff_2, 1)
-      rbf_coeff_2_size_1 = SIZE(rbf_coeff_2, 2)
+      rbf_vec_coeff_v_size_0 = SIZE(rbf_vec_coeff_v, 1)
+      rbf_vec_coeff_v_size_1 = SIZE(rbf_vec_coeff_v, 2)
+      rbf_vec_coeff_v_size_2 = SIZE(rbf_vec_coeff_v, 3)
 
       if (associated(zd_cellidx)) then
          zd_cellidx_ptr = c_loc(zd_cellidx)
@@ -740,12 +726,10 @@ subroutine diffusion_init(theta_ref_mc, &
                                   geofac_n2s_size_1=geofac_n2s_size_1, &
                                   nudgecoeff_e=c_loc(nudgecoeff_e), &
                                   nudgecoeff_e_size_0=nudgecoeff_e_size_0, &
-                                  rbf_coeff_1=c_loc(rbf_coeff_1), &
-                                  rbf_coeff_1_size_0=rbf_coeff_1_size_0, &
-                                  rbf_coeff_1_size_1=rbf_coeff_1_size_1, &
-                                  rbf_coeff_2=c_loc(rbf_coeff_2), &
-                                  rbf_coeff_2_size_0=rbf_coeff_2_size_0, &
-                                  rbf_coeff_2_size_1=rbf_coeff_2_size_1, &
+                                  rbf_vec_coeff_v=c_loc(rbf_vec_coeff_v), &
+                                  rbf_vec_coeff_v_size_0=rbf_vec_coeff_v_size_0, &
+                                  rbf_vec_coeff_v_size_1=rbf_vec_coeff_v_size_1, &
+                                  rbf_vec_coeff_v_size_2=rbf_vec_coeff_v_size_2, &
                                   zd_cellidx=zd_cellidx_ptr, &
                                   zd_cellidx_size_0=zd_cellidx_size_0, &
                                   zd_cellidx_size_1=zd_cellidx_size_1, &
@@ -788,7 +772,6 @@ subroutine diffusion_init(theta_ref_mc, &
       !$acc end host_data
       !$acc end host_data
       !$acc end host_data
-      !$acc end host_data
    end subroutine diffusion_init
 
 end module
\ No newline at end of file
diff --git a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.h b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.h
index 46010ce5d3..a1344d8239 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.h
+++ b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.h
@@ -14,15 +14,15 @@ extern int diffusion_init_wrapper(
     double *geofac_grg_x, int geofac_grg_x_size_0, int geofac_grg_x_size_1,
     double *geofac_grg_y, int geofac_grg_y_size_0, int geofac_grg_y_size_1,
     double *geofac_n2s, int geofac_n2s_size_0, int geofac_n2s_size_1,
-    double *nudgecoeff_e, int nudgecoeff_e_size_0, double *rbf_coeff_1,
-    int rbf_coeff_1_size_0, int rbf_coeff_1_size_1, double *rbf_coeff_2,
-    int rbf_coeff_2_size_0, int rbf_coeff_2_size_1, int *zd_cellidx,
-    int zd_cellidx_size_0, int zd_cellidx_size_1, int *zd_vertidx,
-    int zd_vertidx_size_0, int zd_vertidx_size_1, double *zd_intcoef,
-    int zd_intcoef_size_0, int zd_intcoef_size_1, double *zd_diffcoef,
-    int zd_diffcoef_size_0, int ndyn_substeps, int diffusion_type, int hdiff_w,
-    int hdiff_vn, int hdiff_smag_w, int zdiffu_t, int type_t_diffu,
-    int type_vn_diffu, double hdiff_efdt_ratio, double hdiff_w_efdt_ratio,
+    double *nudgecoeff_e, int nudgecoeff_e_size_0, double *rbf_vec_coeff_v,
+    int rbf_vec_coeff_v_size_0, int rbf_vec_coeff_v_size_1,
+    int rbf_vec_coeff_v_size_2, int *zd_cellidx, int zd_cellidx_size_0,
+    int zd_cellidx_size_1, int *zd_vertidx, int zd_vertidx_size_0,
+    int zd_vertidx_size_1, double *zd_intcoef, int zd_intcoef_size_0,
+    int zd_intcoef_size_1, double *zd_diffcoef, int zd_diffcoef_size_0,
+    int ndyn_substeps, int diffusion_type, int hdiff_w, int hdiff_vn,
+    int hdiff_smag_w, int zdiffu_t, int type_t_diffu, int type_vn_diffu,
+    double hdiff_efdt_ratio, double hdiff_w_efdt_ratio,
     double smagorinski_scaling_factor, int hdiff_temp, double denom_diffu_v,
     double nudge_max_coeff, int itype_sher, int ltkeshs, int backend,
     int on_gpu);
\ No newline at end of file
diff --git a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.py b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.py
index dd26c83bfb..6c64ea4b15 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.py
+++ b/tools/tests/tools/py2fgen/wrappers/references/diffusion/diffusion.py
@@ -365,12 +365,10 @@ def diffusion_init_wrapper(
     geofac_n2s_size_1,
     nudgecoeff_e,
     nudgecoeff_e_size_0,
-    rbf_coeff_1,
-    rbf_coeff_1_size_0,
-    rbf_coeff_1_size_1,
-    rbf_coeff_2,
-    rbf_coeff_2_size_0,
-    rbf_coeff_2_size_1,
+    rbf_vec_coeff_v,
+    rbf_vec_coeff_v_size_0,
+    rbf_vec_coeff_v_size_1,
+    rbf_vec_coeff_v_size_2,
     zd_cellidx,
     zd_cellidx_size_0,
     zd_cellidx_size_1,
@@ -484,21 +482,12 @@ def diffusion_init_wrapper(
 
             nudgecoeff_e = (nudgecoeff_e, (nudgecoeff_e_size_0,), on_gpu, False)
 
-            rbf_coeff_1 = (
-                rbf_coeff_1,
+            rbf_vec_coeff_v = (
+                rbf_vec_coeff_v,
                 (
-                    rbf_coeff_1_size_0,
-                    rbf_coeff_1_size_1,
-                ),
-                on_gpu,
-                False,
-            )
-
-            rbf_coeff_2 = (
-                rbf_coeff_2,
-                (
-                    rbf_coeff_2_size_0,
-                    rbf_coeff_2_size_1,
+                    rbf_vec_coeff_v_size_0,
+                    rbf_vec_coeff_v_size_1,
+                    rbf_vec_coeff_v_size_2,
                 ),
                 on_gpu,
                 False,
@@ -561,8 +550,7 @@ def diffusion_init_wrapper(
                 geofac_grg_y=geofac_grg_y,
                 geofac_n2s=geofac_n2s,
                 nudgecoeff_e=nudgecoeff_e,
-                rbf_coeff_1=rbf_coeff_1,
-                rbf_coeff_2=rbf_coeff_2,
+                rbf_vec_coeff_v=rbf_vec_coeff_v,
                 zd_cellidx=zd_cellidx,
                 zd_vertidx=zd_vertidx,
                 zd_intcoef=zd_intcoef,
@@ -730,34 +718,18 @@ def diffusion_init_wrapper(
                     )
                     logger.debug(msg)
 
-                    rbf_coeff_1_arr = (
-                        _conversion.as_array(ffi, rbf_coeff_1, _definitions.FLOAT64)
-                        if rbf_coeff_1 is not None
-                        else None
-                    )
-                    msg = "shape of rbf_coeff_1 after computation = %s" % str(
-                        rbf_coeff_1_arr.shape if rbf_coeff_1 is not None else "None"
-                    )
-                    logger.debug(msg)
-                    msg = (
-                        "rbf_coeff_1 after computation: %s" % str(rbf_coeff_1_arr)
-                        if rbf_coeff_1 is not None
-                        else "None"
-                    )
-                    logger.debug(msg)
-
-                    rbf_coeff_2_arr = (
-                        _conversion.as_array(ffi, rbf_coeff_2, _definitions.FLOAT64)
-                        if rbf_coeff_2 is not None
+                    rbf_vec_coeff_v_arr = (
+                        _conversion.as_array(ffi, rbf_vec_coeff_v, _definitions.FLOAT64)
+                        if rbf_vec_coeff_v is not None
                         else None
                     )
-                    msg = "shape of rbf_coeff_2 after computation = %s" % str(
-                        rbf_coeff_2_arr.shape if rbf_coeff_2 is not None else "None"
+                    msg = "shape of rbf_vec_coeff_v after computation = %s" % str(
+                        rbf_vec_coeff_v_arr.shape if rbf_vec_coeff_v is not None else "None"
                     )
                     logger.debug(msg)
                     msg = (
-                        "rbf_coeff_2 after computation: %s" % str(rbf_coeff_2_arr)
-                        if rbf_coeff_2 is not None
+                        "rbf_vec_coeff_v after computation: %s" % str(rbf_vec_coeff_v_arr)
+                        if rbf_vec_coeff_v is not None
                         else "None"
                     )
                     logger.debug(msg)
diff --git a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.f90 b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.f90
index 6b29e00d95..322658ce45 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.f90
+++ b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.f90
@@ -382,12 +382,10 @@ function solve_nh_init_wrapper(c_lin_e, &
                                      e_bln_c_s, &
                                      e_bln_c_s_size_0, &
                                      e_bln_c_s_size_1, &
-                                     rbf_coeff_1, &
-                                     rbf_coeff_1_size_0, &
-                                     rbf_coeff_1_size_1, &
-                                     rbf_coeff_2, &
-                                     rbf_coeff_2_size_0, &
-                                     rbf_coeff_2_size_1, &
+                                     rbf_vec_coeff_v, &
+                                     rbf_vec_coeff_v_size_0, &
+                                     rbf_vec_coeff_v_size_1, &
+                                     rbf_vec_coeff_v_size_2, &
                                      geofac_div, &
                                      geofac_div_size_0, &
                                      geofac_div_size_1, &
@@ -457,10 +455,10 @@ function solve_nh_init_wrapper(c_lin_e, &
                                      zdiff_gradp_size_0, &
                                      zdiff_gradp_size_1, &
                                      zdiff_gradp_size_2, &
-                                     vertoffset_gradp, &
-                                     vertoffset_gradp_size_0, &
-                                     vertoffset_gradp_size_1, &
-                                     vertoffset_gradp_size_2, &
+                                     vertidx_gradp, &
+                                     vertidx_gradp_size_0, &
+                                     vertidx_gradp_size_1, &
+                                     vertidx_gradp_size_2, &
                                      pg_edgeidx, &
                                      pg_edgeidx_size_0, &
                                      pg_vertidx, &
@@ -580,17 +578,13 @@ function solve_nh_init_wrapper(c_lin_e, &
 
          integer(c_int), value :: e_bln_c_s_size_1
 
-         type(c_ptr), value, target :: rbf_coeff_1
+         type(c_ptr), value, target :: rbf_vec_coeff_v
 
-         integer(c_int), value :: rbf_coeff_1_size_0
+         integer(c_int), value :: rbf_vec_coeff_v_size_0
 
-         integer(c_int), value :: rbf_coeff_1_size_1
+         integer(c_int), value :: rbf_vec_coeff_v_size_1
 
-         type(c_ptr), value, target :: rbf_coeff_2
-
-         integer(c_int), value :: rbf_coeff_2_size_0
-
-         integer(c_int), value :: rbf_coeff_2_size_1
+         integer(c_int), value :: rbf_vec_coeff_v_size_2
 
          type(c_ptr), value, target :: geofac_div
 
@@ -730,13 +724,13 @@ function solve_nh_init_wrapper(c_lin_e, &
 
          integer(c_int), value :: zdiff_gradp_size_2
 
-         type(c_ptr), value, target :: vertoffset_gradp
+         type(c_ptr), value, target :: vertidx_gradp
 
-         integer(c_int), value :: vertoffset_gradp_size_0
+         integer(c_int), value :: vertidx_gradp_size_0
 
-         integer(c_int), value :: vertoffset_gradp_size_1
+         integer(c_int), value :: vertidx_gradp_size_1
 
-         integer(c_int), value :: vertoffset_gradp_size_2
+         integer(c_int), value :: vertidx_gradp_size_2
 
          type(c_ptr), value, target :: pg_edgeidx
 
@@ -1474,8 +1468,7 @@ subroutine solve_nh_init(c_lin_e, &
                             pos_on_tplane_e_2, &
                             rbf_vec_coeff_e, &
                             e_bln_c_s, &
-                            rbf_coeff_1, &
-                            rbf_coeff_2, &
+                            rbf_vec_coeff_v, &
                             geofac_div, &
                             geofac_n2s, &
                             geofac_grg_x, &
@@ -1500,7 +1493,7 @@ subroutine solve_nh_init(c_lin_e, &
                             theta_ref_me, &
                             ddxn_z_full, &
                             zdiff_gradp, &
-                            vertoffset_gradp, &
+                            vertidx_gradp, &
                             pg_edgeidx, &
                             pg_vertidx, &
                             pg_exdist, &
@@ -1562,9 +1555,7 @@ subroutine solve_nh_init(c_lin_e, &
 
       real(c_double), dimension(:, :), target :: e_bln_c_s
 
-      real(c_double), dimension(:, :), target :: rbf_coeff_1
-
-      real(c_double), dimension(:, :), target :: rbf_coeff_2
+      real(c_double), dimension(:, :, :), target :: rbf_vec_coeff_v
 
       real(c_double), dimension(:, :), target :: geofac_div
 
@@ -1614,7 +1605,7 @@ subroutine solve_nh_init(c_lin_e, &
 
       real(c_double), dimension(:, :, :), target :: zdiff_gradp
 
-      integer(c_int), dimension(:, :, :), target :: vertoffset_gradp
+      integer(c_int), dimension(:, :, :), target :: vertidx_gradp
 
       integer(c_int), dimension(:), pointer :: pg_edgeidx
 
@@ -1734,13 +1725,11 @@ subroutine solve_nh_init(c_lin_e, &
 
       integer(c_int) :: e_bln_c_s_size_1
 
-      integer(c_int) :: rbf_coeff_1_size_0
+      integer(c_int) :: rbf_vec_coeff_v_size_0
 
-      integer(c_int) :: rbf_coeff_1_size_1
+      integer(c_int) :: rbf_vec_coeff_v_size_1
 
-      integer(c_int) :: rbf_coeff_2_size_0
-
-      integer(c_int) :: rbf_coeff_2_size_1
+      integer(c_int) :: rbf_vec_coeff_v_size_2
 
       integer(c_int) :: geofac_div_size_0
 
@@ -1832,11 +1821,11 @@ subroutine solve_nh_init(c_lin_e, &
 
       integer(c_int) :: zdiff_gradp_size_2
 
-      integer(c_int) :: vertoffset_gradp_size_0
+      integer(c_int) :: vertidx_gradp_size_0
 
-      integer(c_int) :: vertoffset_gradp_size_1
+      integer(c_int) :: vertidx_gradp_size_1
 
-      integer(c_int) :: vertoffset_gradp_size_2
+      integer(c_int) :: vertidx_gradp_size_2
 
       integer(c_int) :: pg_edgeidx_size_0
 
@@ -1904,8 +1893,7 @@ subroutine solve_nh_init(c_lin_e, &
       !$acc host_data use_device(pos_on_tplane_e_2)
       !$acc host_data use_device(rbf_vec_coeff_e)
       !$acc host_data use_device(e_bln_c_s)
-      !$acc host_data use_device(rbf_coeff_1)
-      !$acc host_data use_device(rbf_coeff_2)
+      !$acc host_data use_device(rbf_vec_coeff_v)
       !$acc host_data use_device(geofac_div)
       !$acc host_data use_device(geofac_n2s)
       !$acc host_data use_device(geofac_grg_x)
@@ -1930,7 +1918,7 @@ subroutine solve_nh_init(c_lin_e, &
       !$acc host_data use_device(theta_ref_me)
       !$acc host_data use_device(ddxn_z_full)
       !$acc host_data use_device(zdiff_gradp)
-      !$acc host_data use_device(vertoffset_gradp)
+      !$acc host_data use_device(vertidx_gradp)
       !$acc host_data use_device(ddqz_z_full_e)
       !$acc host_data use_device(ddxt_z_full)
       !$acc host_data use_device(wgtfac_e)
@@ -1979,11 +1967,9 @@ subroutine solve_nh_init(c_lin_e, &
       e_bln_c_s_size_0 = SIZE(e_bln_c_s, 1)
       e_bln_c_s_size_1 = SIZE(e_bln_c_s, 2)
 
-      rbf_coeff_1_size_0 = SIZE(rbf_coeff_1, 1)
-      rbf_coeff_1_size_1 = SIZE(rbf_coeff_1, 2)
-
-      rbf_coeff_2_size_0 = SIZE(rbf_coeff_2, 1)
-      rbf_coeff_2_size_1 = SIZE(rbf_coeff_2, 2)
+      rbf_vec_coeff_v_size_0 = SIZE(rbf_vec_coeff_v, 1)
+      rbf_vec_coeff_v_size_1 = SIZE(rbf_vec_coeff_v, 2)
+      rbf_vec_coeff_v_size_2 = SIZE(rbf_vec_coeff_v, 3)
 
       geofac_div_size_0 = SIZE(geofac_div, 1)
       geofac_div_size_1 = SIZE(geofac_div, 2)
@@ -2054,9 +2040,9 @@ subroutine solve_nh_init(c_lin_e, &
       zdiff_gradp_size_1 = SIZE(zdiff_gradp, 2)
       zdiff_gradp_size_2 = SIZE(zdiff_gradp, 3)
 
-      vertoffset_gradp_size_0 = SIZE(vertoffset_gradp, 1)
-      vertoffset_gradp_size_1 = SIZE(vertoffset_gradp, 2)
-      vertoffset_gradp_size_2 = SIZE(vertoffset_gradp, 3)
+      vertidx_gradp_size_0 = SIZE(vertidx_gradp, 1)
+      vertidx_gradp_size_1 = SIZE(vertidx_gradp, 2)
+      vertidx_gradp_size_2 = SIZE(vertidx_gradp, 3)
 
       ddqz_z_full_e_size_0 = SIZE(ddqz_z_full_e, 1)
       ddqz_z_full_e_size_1 = SIZE(ddqz_z_full_e, 2)
@@ -2129,12 +2115,10 @@ subroutine solve_nh_init(c_lin_e, &
                                  e_bln_c_s=c_loc(e_bln_c_s), &
                                  e_bln_c_s_size_0=e_bln_c_s_size_0, &
                                  e_bln_c_s_size_1=e_bln_c_s_size_1, &
-                                 rbf_coeff_1=c_loc(rbf_coeff_1), &
-                                 rbf_coeff_1_size_0=rbf_coeff_1_size_0, &
-                                 rbf_coeff_1_size_1=rbf_coeff_1_size_1, &
-                                 rbf_coeff_2=c_loc(rbf_coeff_2), &
-                                 rbf_coeff_2_size_0=rbf_coeff_2_size_0, &
-                                 rbf_coeff_2_size_1=rbf_coeff_2_size_1, &
+                                 rbf_vec_coeff_v=c_loc(rbf_vec_coeff_v), &
+                                 rbf_vec_coeff_v_size_0=rbf_vec_coeff_v_size_0, &
+                                 rbf_vec_coeff_v_size_1=rbf_vec_coeff_v_size_1, &
+                                 rbf_vec_coeff_v_size_2=rbf_vec_coeff_v_size_2, &
                                  geofac_div=c_loc(geofac_div), &
                                  geofac_div_size_0=geofac_div_size_0, &
                                  geofac_div_size_1=geofac_div_size_1, &
@@ -2204,10 +2188,10 @@ subroutine solve_nh_init(c_lin_e, &
                                  zdiff_gradp_size_0=zdiff_gradp_size_0, &
                                  zdiff_gradp_size_1=zdiff_gradp_size_1, &
                                  zdiff_gradp_size_2=zdiff_gradp_size_2, &
-                                 vertoffset_gradp=c_loc(vertoffset_gradp), &
-                                 vertoffset_gradp_size_0=vertoffset_gradp_size_0, &
-                                 vertoffset_gradp_size_1=vertoffset_gradp_size_1, &
-                                 vertoffset_gradp_size_2=vertoffset_gradp_size_2, &
+                                 vertidx_gradp=c_loc(vertidx_gradp), &
+                                 vertidx_gradp_size_0=vertidx_gradp_size_0, &
+                                 vertidx_gradp_size_1=vertidx_gradp_size_1, &
+                                 vertidx_gradp_size_2=vertidx_gradp_size_2, &
                                  pg_edgeidx=pg_edgeidx_ptr, &
                                  pg_edgeidx_size_0=pg_edgeidx_size_0, &
                                  pg_vertidx=pg_vertidx_ptr, &
@@ -2319,7 +2303,6 @@ subroutine solve_nh_init(c_lin_e, &
       !$acc end host_data
       !$acc end host_data
       !$acc end host_data
-      !$acc end host_data
    end subroutine solve_nh_init
 
 end module
\ No newline at end of file
diff --git a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.h b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.h
index 660a2fc279..e12a36d50e 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.h
+++ b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.h
@@ -47,52 +47,51 @@ extern int solve_nh_init_wrapper(
     double *pos_on_tplane_e_2, int pos_on_tplane_e_2_size_0,
     int pos_on_tplane_e_2_size_1, double *rbf_vec_coeff_e,
     int rbf_vec_coeff_e_size_0, int rbf_vec_coeff_e_size_1, double *e_bln_c_s,
-    int e_bln_c_s_size_0, int e_bln_c_s_size_1, double *rbf_coeff_1,
-    int rbf_coeff_1_size_0, int rbf_coeff_1_size_1, double *rbf_coeff_2,
-    int rbf_coeff_2_size_0, int rbf_coeff_2_size_1, double *geofac_div,
-    int geofac_div_size_0, int geofac_div_size_1, double *geofac_n2s,
-    int geofac_n2s_size_0, int geofac_n2s_size_1, double *geofac_grg_x,
-    int geofac_grg_x_size_0, int geofac_grg_x_size_1, double *geofac_grg_y,
-    int geofac_grg_y_size_0, int geofac_grg_y_size_1, double *nudgecoeff_e,
-    int nudgecoeff_e_size_0, int *mask_prog_halo_c, int mask_prog_halo_c_size_0,
-    double *rayleigh_w, int rayleigh_w_size_0, double *exner_exfac,
-    int exner_exfac_size_0, int exner_exfac_size_1, double *exner_ref_mc,
-    int exner_ref_mc_size_0, int exner_ref_mc_size_1, double *wgtfac_c,
-    int wgtfac_c_size_0, int wgtfac_c_size_1, double *wgtfacq_c,
-    int wgtfacq_c_size_0, int wgtfacq_c_size_1, double *inv_ddqz_z_full,
-    int inv_ddqz_z_full_size_0, int inv_ddqz_z_full_size_1, double *rho_ref_mc,
-    int rho_ref_mc_size_0, int rho_ref_mc_size_1, double *theta_ref_mc,
-    int theta_ref_mc_size_0, int theta_ref_mc_size_1, double *vwind_expl_wgt,
-    int vwind_expl_wgt_size_0, double *d_exner_dz_ref_ic,
-    int d_exner_dz_ref_ic_size_0, int d_exner_dz_ref_ic_size_1,
-    double *ddqz_z_half, int ddqz_z_half_size_0, int ddqz_z_half_size_1,
-    double *theta_ref_ic, int theta_ref_ic_size_0, int theta_ref_ic_size_1,
-    double *d2dexdz2_fac1_mc, int d2dexdz2_fac1_mc_size_0,
-    int d2dexdz2_fac1_mc_size_1, double *d2dexdz2_fac2_mc,
-    int d2dexdz2_fac2_mc_size_0, int d2dexdz2_fac2_mc_size_1,
-    double *rho_ref_me, int rho_ref_me_size_0, int rho_ref_me_size_1,
-    double *theta_ref_me, int theta_ref_me_size_0, int theta_ref_me_size_1,
-    double *ddxn_z_full, int ddxn_z_full_size_0, int ddxn_z_full_size_1,
-    double *zdiff_gradp, int zdiff_gradp_size_0, int zdiff_gradp_size_1,
-    int zdiff_gradp_size_2, int *vertoffset_gradp, int vertoffset_gradp_size_0,
-    int vertoffset_gradp_size_1, int vertoffset_gradp_size_2, int *pg_edgeidx,
-    int pg_edgeidx_size_0, int *pg_vertidx, int pg_vertidx_size_0,
-    double *pg_exdist, int pg_exdist_size_0, double *ddqz_z_full_e,
-    int ddqz_z_full_e_size_0, int ddqz_z_full_e_size_1, double *ddxt_z_full,
-    int ddxt_z_full_size_0, int ddxt_z_full_size_1, double *wgtfac_e,
-    int wgtfac_e_size_0, int wgtfac_e_size_1, double *wgtfacq_e,
-    int wgtfacq_e_size_0, int wgtfacq_e_size_1, double *vwind_impl_wgt,
-    int vwind_impl_wgt_size_0, double *hmask_dd3d, int hmask_dd3d_size_0,
-    double *scalfac_dd3d, int scalfac_dd3d_size_0, double *coeff1_dwdz,
-    int coeff1_dwdz_size_0, int coeff1_dwdz_size_1, double *coeff2_dwdz,
-    int coeff2_dwdz_size_0, int coeff2_dwdz_size_1, double *coeff_gradekin,
-    int coeff_gradekin_size_0, int coeff_gradekin_size_1, int *c_owner_mask,
-    int c_owner_mask_size_0, int itime_scheme, int iadv_rhotheta,
-    int igradp_method, int rayleigh_type, double rayleigh_coeff,
-    int divdamp_order, int is_iau_active, double iau_wgt_dyn, int divdamp_type,
-    double divdamp_trans_start, double divdamp_trans_end, int l_vert_nested,
-    int ldeepatmo, double rhotheta_offctr, double veladv_offctr,
-    double nudge_max_coeff, double divdamp_fac, double divdamp_fac2,
-    double divdamp_fac3, double divdamp_fac4, double divdamp_z,
-    double divdamp_z2, double divdamp_z3, double divdamp_z4, int nflat_gradp,
-    int backend, int on_gpu);
\ No newline at end of file
+    int e_bln_c_s_size_0, int e_bln_c_s_size_1, double *rbf_vec_coeff_v,
+    int rbf_vec_coeff_v_size_0, int rbf_vec_coeff_v_size_1,
+    int rbf_vec_coeff_v_size_2, double *geofac_div, int geofac_div_size_0,
+    int geofac_div_size_1, double *geofac_n2s, int geofac_n2s_size_0,
+    int geofac_n2s_size_1, double *geofac_grg_x, int geofac_grg_x_size_0,
+    int geofac_grg_x_size_1, double *geofac_grg_y, int geofac_grg_y_size_0,
+    int geofac_grg_y_size_1, double *nudgecoeff_e, int nudgecoeff_e_size_0,
+    int *mask_prog_halo_c, int mask_prog_halo_c_size_0, double *rayleigh_w,
+    int rayleigh_w_size_0, double *exner_exfac, int exner_exfac_size_0,
+    int exner_exfac_size_1, double *exner_ref_mc, int exner_ref_mc_size_0,
+    int exner_ref_mc_size_1, double *wgtfac_c, int wgtfac_c_size_0,
+    int wgtfac_c_size_1, double *wgtfacq_c, int wgtfacq_c_size_0,
+    int wgtfacq_c_size_1, double *inv_ddqz_z_full, int inv_ddqz_z_full_size_0,
+    int inv_ddqz_z_full_size_1, double *rho_ref_mc, int rho_ref_mc_size_0,
+    int rho_ref_mc_size_1, double *theta_ref_mc, int theta_ref_mc_size_0,
+    int theta_ref_mc_size_1, double *vwind_expl_wgt, int vwind_expl_wgt_size_0,
+    double *d_exner_dz_ref_ic, int d_exner_dz_ref_ic_size_0,
+    int d_exner_dz_ref_ic_size_1, double *ddqz_z_half, int ddqz_z_half_size_0,
+    int ddqz_z_half_size_1, double *theta_ref_ic, int theta_ref_ic_size_0,
+    int theta_ref_ic_size_1, double *d2dexdz2_fac1_mc,
+    int d2dexdz2_fac1_mc_size_0, int d2dexdz2_fac1_mc_size_1,
+    double *d2dexdz2_fac2_mc, int d2dexdz2_fac2_mc_size_0,
+    int d2dexdz2_fac2_mc_size_1, double *rho_ref_me, int rho_ref_me_size_0,
+    int rho_ref_me_size_1, double *theta_ref_me, int theta_ref_me_size_0,
+    int theta_ref_me_size_1, double *ddxn_z_full, int ddxn_z_full_size_0,
+    int ddxn_z_full_size_1, double *zdiff_gradp, int zdiff_gradp_size_0,
+    int zdiff_gradp_size_1, int zdiff_gradp_size_2, int *vertidx_gradp,
+    int vertidx_gradp_size_0, int vertidx_gradp_size_1,
+    int vertidx_gradp_size_2, int *pg_edgeidx, int pg_edgeidx_size_0,
+    int *pg_vertidx, int pg_vertidx_size_0, double *pg_exdist,
+    int pg_exdist_size_0, double *ddqz_z_full_e, int ddqz_z_full_e_size_0,
+    int ddqz_z_full_e_size_1, double *ddxt_z_full, int ddxt_z_full_size_0,
+    int ddxt_z_full_size_1, double *wgtfac_e, int wgtfac_e_size_0,
+    int wgtfac_e_size_1, double *wgtfacq_e, int wgtfacq_e_size_0,
+    int wgtfacq_e_size_1, double *vwind_impl_wgt, int vwind_impl_wgt_size_0,
+    double *hmask_dd3d, int hmask_dd3d_size_0, double *scalfac_dd3d,
+    int scalfac_dd3d_size_0, double *coeff1_dwdz, int coeff1_dwdz_size_0,
+    int coeff1_dwdz_size_1, double *coeff2_dwdz, int coeff2_dwdz_size_0,
+    int coeff2_dwdz_size_1, double *coeff_gradekin, int coeff_gradekin_size_0,
+    int coeff_gradekin_size_1, int *c_owner_mask, int c_owner_mask_size_0,
+    int itime_scheme, int iadv_rhotheta, int igradp_method, int rayleigh_type,
+    double rayleigh_coeff, int divdamp_order, int is_iau_active,
+    double iau_wgt_dyn, int divdamp_type, double divdamp_trans_start,
+    double divdamp_trans_end, int l_vert_nested, int ldeepatmo,
+    double rhotheta_offctr, double veladv_offctr, double nudge_max_coeff,
+    double divdamp_fac, double divdamp_fac2, double divdamp_fac3,
+    double divdamp_fac4, double divdamp_z, double divdamp_z2, double divdamp_z3,
+    double divdamp_z4, int nflat_gradp, int backend, int on_gpu);
\ No newline at end of file
diff --git a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.py b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.py
index 15ead38e9a..e847aa080a 100644
--- a/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.py
+++ b/tools/tests/tools/py2fgen/wrappers/references/dycore/dycore.py
@@ -1197,12 +1197,10 @@ def solve_nh_init_wrapper(
     e_bln_c_s,
     e_bln_c_s_size_0,
     e_bln_c_s_size_1,
-    rbf_coeff_1,
-    rbf_coeff_1_size_0,
-    rbf_coeff_1_size_1,
-    rbf_coeff_2,
-    rbf_coeff_2_size_0,
-    rbf_coeff_2_size_1,
+    rbf_vec_coeff_v,
+    rbf_vec_coeff_v_size_0,
+    rbf_vec_coeff_v_size_1,
+    rbf_vec_coeff_v_size_2,
     geofac_div,
     geofac_div_size_0,
     geofac_div_size_1,
@@ -1272,10 +1270,10 @@ def solve_nh_init_wrapper(
     zdiff_gradp_size_0,
     zdiff_gradp_size_1,
     zdiff_gradp_size_2,
-    vertoffset_gradp,
-    vertoffset_gradp_size_0,
-    vertoffset_gradp_size_1,
-    vertoffset_gradp_size_2,
+    vertidx_gradp,
+    vertidx_gradp_size_0,
+    vertidx_gradp_size_1,
+    vertidx_gradp_size_2,
     pg_edgeidx,
     pg_edgeidx_size_0,
     pg_vertidx,
@@ -1440,21 +1438,12 @@ def solve_nh_init_wrapper(
                 False,
             )
 
-            rbf_coeff_1 = (
-                rbf_coeff_1,
+            rbf_vec_coeff_v = (
+                rbf_vec_coeff_v,
                 (
-                    rbf_coeff_1_size_0,
-                    rbf_coeff_1_size_1,
-                ),
-                on_gpu,
-                False,
-            )
-
-            rbf_coeff_2 = (
-                rbf_coeff_2,
-                (
-                    rbf_coeff_2_size_0,
-                    rbf_coeff_2_size_1,
+                    rbf_vec_coeff_v_size_0,
+                    rbf_vec_coeff_v_size_1,
+                    rbf_vec_coeff_v_size_2,
                 ),
                 on_gpu,
                 False,
@@ -1669,12 +1658,12 @@ def solve_nh_init_wrapper(
                 False,
             )
 
-            vertoffset_gradp = (
-                vertoffset_gradp,
+            vertidx_gradp = (
+                vertidx_gradp,
                 (
-                    vertoffset_gradp_size_0,
-                    vertoffset_gradp_size_1,
-                    vertoffset_gradp_size_2,
+                    vertidx_gradp_size_0,
+                    vertidx_gradp_size_1,
+                    vertidx_gradp_size_2,
                 ),
                 on_gpu,
                 False,
@@ -1790,8 +1779,7 @@ def solve_nh_init_wrapper(
                 pos_on_tplane_e_2=pos_on_tplane_e_2,
                 rbf_vec_coeff_e=rbf_vec_coeff_e,
                 e_bln_c_s=e_bln_c_s,
-                rbf_coeff_1=rbf_coeff_1,
-                rbf_coeff_2=rbf_coeff_2,
+                rbf_vec_coeff_v=rbf_vec_coeff_v,
                 geofac_div=geofac_div,
                 geofac_n2s=geofac_n2s,
                 geofac_grg_x=geofac_grg_x,
@@ -1816,7 +1804,7 @@ def solve_nh_init_wrapper(
                 theta_ref_me=theta_ref_me,
                 ddxn_z_full=ddxn_z_full,
                 zdiff_gradp=zdiff_gradp,
-                vertoffset_gradp=vertoffset_gradp,
+                vertidx_gradp=vertidx_gradp,
                 pg_edgeidx=pg_edgeidx,
                 pg_vertidx=pg_vertidx,
                 pg_exdist=pg_exdist,
@@ -2019,34 +2007,18 @@ def solve_nh_init_wrapper(
                     )
                     logger.debug(msg)
 
-                    rbf_coeff_1_arr = (
-                        _conversion.as_array(ffi, rbf_coeff_1, _definitions.FLOAT64)
-                        if rbf_coeff_1 is not None
-                        else None
-                    )
-                    msg = "shape of rbf_coeff_1 after computation = %s" % str(
-                        rbf_coeff_1_arr.shape if rbf_coeff_1 is not None else "None"
-                    )
-                    logger.debug(msg)
-                    msg = (
-                        "rbf_coeff_1 after computation: %s" % str(rbf_coeff_1_arr)
-                        if rbf_coeff_1 is not None
-                        else "None"
-                    )
-                    logger.debug(msg)
-
-                    rbf_coeff_2_arr = (
-                        _conversion.as_array(ffi, rbf_coeff_2, _definitions.FLOAT64)
-                        if rbf_coeff_2 is not None
+                    rbf_vec_coeff_v_arr = (
+                        _conversion.as_array(ffi, rbf_vec_coeff_v, _definitions.FLOAT64)
+                        if rbf_vec_coeff_v is not None
                         else None
                     )
-                    msg = "shape of rbf_coeff_2 after computation = %s" % str(
-                        rbf_coeff_2_arr.shape if rbf_coeff_2 is not None else "None"
+                    msg = "shape of rbf_vec_coeff_v after computation = %s" % str(
+                        rbf_vec_coeff_v_arr.shape if rbf_vec_coeff_v is not None else "None"
                     )
                     logger.debug(msg)
                     msg = (
-                        "rbf_coeff_2 after computation: %s" % str(rbf_coeff_2_arr)
-                        if rbf_coeff_2 is not None
+                        "rbf_vec_coeff_v after computation: %s" % str(rbf_vec_coeff_v_arr)
+                        if rbf_vec_coeff_v is not None
                         else "None"
                     )
                     logger.debug(msg)
@@ -2435,18 +2407,18 @@ def solve_nh_init_wrapper(
                     )
                     logger.debug(msg)
 
-                    vertoffset_gradp_arr = (
-                        _conversion.as_array(ffi, vertoffset_gradp, _definitions.INT32)
-                        if vertoffset_gradp is not None
+                    vertidx_gradp_arr = (
+                        _conversion.as_array(ffi, vertidx_gradp, _definitions.INT32)
+                        if vertidx_gradp is not None
                         else None
                     )
-                    msg = "shape of vertoffset_gradp after computation = %s" % str(
-                        vertoffset_gradp_arr.shape if vertoffset_gradp is not None else "None"
+                    msg = "shape of vertidx_gradp after computation = %s" % str(
+                        vertidx_gradp_arr.shape if vertidx_gradp is not None else "None"
                     )
                     logger.debug(msg)
                     msg = (
-                        "vertoffset_gradp after computation: %s" % str(vertoffset_gradp_arr)
-                        if vertoffset_gradp is not None
+                        "vertidx_gradp after computation: %s" % str(vertidx_gradp_arr)
+                        if vertidx_gradp is not None
                         else "None"
                     )
                     logger.debug(msg)
diff --git a/tools/tests/tools/py2fgen/wrappers/test_diffusion_wrapper.py b/tools/tests/tools/py2fgen/wrappers/test_diffusion_wrapper.py
index f85ec964b7..6f29687e41 100644
--- a/tools/tests/tools/py2fgen/wrappers/test_diffusion_wrapper.py
+++ b/tools/tests/tools/py2fgen/wrappers/test_diffusion_wrapper.py
@@ -96,8 +96,17 @@ def test_diffusion_wrapper_granule_inputs(
     geofac_grg_y = test_utils.array_to_array_info(geofac_grg_y_field.ndarray)
     geofac_n2s = test_utils.array_to_array_info(interpolation_savepoint.geofac_n2s().ndarray)
     nudgecoeff_e = test_utils.array_to_array_info(interpolation_savepoint.nudgecoeff_e().ndarray)
-    rbf_coeff_1 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v1().ndarray)
-    rbf_coeff_2 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v2().ndarray)
+
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    rbf_vec_coeff_v_array = np.squeeze(
+        interpolation_savepoint.serializer.read(
+            "rbf_vec_coeff_v", interpolation_savepoint.savepoint
+        ).astype(float)
+    )
+    rbf_vec_coeff_v_array = interpolation_savepoint._reduce_to_dim_size(
+        rbf_vec_coeff_v_array, [dims.V2EDim, dims.V2EDim, dims.VertexDim]
+    )
+    rbf_vec_coeff_v = test_utils.array_to_array_info(rbf_vec_coeff_v_array)
 
     # --- Extract Diagnostic and Prognostic State Parameters ---
     hdef_ic = test_utils.array_to_array_info(savepoint_diffusion_init.hdef_ic().ndarray)
@@ -158,8 +167,7 @@ def test_diffusion_wrapper_granule_inputs(
             geofac_grg_y=geofac_grg_y,
             geofac_n2s=geofac_n2s,
             nudgecoeff_e=nudgecoeff_e,
-            rbf_coeff_1=rbf_coeff_1,
-            rbf_coeff_2=rbf_coeff_2,
+            rbf_vec_coeff_v=rbf_vec_coeff_v,
             zd_cellidx=zd_cellidx,
             zd_vertidx=zd_vertidx,
             zd_intcoef=zd_intcoef,
@@ -325,8 +333,16 @@ def test_diffusion_wrapper_single_step(
     geofac_grg_y = test_utils.array_to_array_info(geofac_grg_y_field.ndarray)
     geofac_n2s = test_utils.array_to_array_info(interpolation_savepoint.geofac_n2s().ndarray)
     nudgecoeff_e = test_utils.array_to_array_info(interpolation_savepoint.nudgecoeff_e().ndarray)
-    rbf_coeff_1 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v1().ndarray)
-    rbf_coeff_2 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v2().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    rbf_vec_coeff_v_array = np.squeeze(
+        interpolation_savepoint.serializer.read(
+            "rbf_vec_coeff_v", interpolation_savepoint.savepoint
+        ).astype(float)
+    )
+    rbf_vec_coeff_v_array = interpolation_savepoint._reduce_to_dim_size(
+        rbf_vec_coeff_v_array, [dims.V2EDim, dims.V2EDim, dims.VertexDim]
+    )
+    rbf_vec_coeff_v = test_utils.array_to_array_info(rbf_vec_coeff_v_array)
 
     # Diagnostic state parameters
     hdef_ic = test_utils.array_to_array_info(savepoint_diffusion_init.hdef_ic().ndarray)
@@ -355,8 +371,7 @@ def test_diffusion_wrapper_single_step(
         geofac_grg_y=geofac_grg_y,
         geofac_n2s=geofac_n2s,
         nudgecoeff_e=nudgecoeff_e,
-        rbf_coeff_1=rbf_coeff_1,
-        rbf_coeff_2=rbf_coeff_2,
+        rbf_vec_coeff_v=rbf_vec_coeff_v,
         zd_cellidx=zd_cellidx,
         zd_vertidx=zd_vertidx,
         zd_intcoef=zd_intcoef,
diff --git a/tools/tests/tools/py2fgen/wrappers/test_dycore_wrapper.py b/tools/tests/tools/py2fgen/wrappers/test_dycore_wrapper.py
index 726a86e7bd..717aedb692 100644
--- a/tools/tests/tools/py2fgen/wrappers/test_dycore_wrapper.py
+++ b/tools/tests/tools/py2fgen/wrappers/test_dycore_wrapper.py
@@ -74,7 +74,10 @@ def solve_nh_init(
     exner_exfac = test_utils.array_to_array_info(metrics_savepoint.exner_exfac().ndarray)
     exner_ref_mc = test_utils.array_to_array_info(metrics_savepoint.exner_ref_mc().ndarray)
     wgtfac_c = test_utils.array_to_array_info(metrics_savepoint.wgtfac_c().ndarray)
-    wgtfacq_c = test_utils.array_to_array_info(metrics_savepoint.wgtfacq_c().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    wgtfacq_c = test_utils.array_to_array_info(
+        metrics_savepoint._get_field("wgtfacq_c", dims.CellDim, dims.KDim).ndarray
+    )
     inv_ddqz_z_full = test_utils.array_to_array_info(metrics_savepoint.inv_ddqz_z_full().ndarray)
     rho_ref_mc = test_utils.array_to_array_info(metrics_savepoint.rho_ref_mc().ndarray)
     theta_ref_mc = test_utils.array_to_array_info(metrics_savepoint.theta_ref_mc().ndarray)
@@ -90,15 +93,13 @@ def solve_nh_init(
     theta_ref_me = test_utils.array_to_array_info(metrics_savepoint.theta_ref_me().ndarray)
     ddxn_z_full = test_utils.array_to_array_info(metrics_savepoint.ddxn_z_full().ndarray)
 
-    zdiff_gradp_field = metrics_savepoint._get_field(
-        "zdiff_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim
-    )
-    zdiff_gradp = test_utils.array_to_array_info(zdiff_gradp_field.ndarray)
-
-    vertoffset_gradp_field = metrics_savepoint._get_field(
-        "vertoffset_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32
+    zdiff_gradp = test_utils.array_to_array_info(metrics_savepoint.zdiff_gradp().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    vertidx_gradp = test_utils.array_to_array_info(
+        metrics_savepoint._get_field(
+            "vertidx_gradp", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32
+        ).ndarray
     )
-    vertoffset_gradp = test_utils.array_to_array_info(vertoffset_gradp_field.ndarray)
 
     pg_edgeidx = test_utils.array_to_array_info(metrics_savepoint.pg_edgeidx())
     pg_vertidx = test_utils.array_to_array_info(metrics_savepoint.pg_vertidx())
@@ -106,7 +107,10 @@ def solve_nh_init(
     ddqz_z_full_e = test_utils.array_to_array_info(metrics_savepoint.ddqz_z_full_e().ndarray)
     ddxt_z_full = test_utils.array_to_array_info(metrics_savepoint.ddxt_z_full().ndarray)
     wgtfac_e = test_utils.array_to_array_info(metrics_savepoint.wgtfac_e().ndarray)
-    wgtfacq_e = test_utils.array_to_array_info(metrics_savepoint.wgtfacq_e().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    wgtfacq_e = test_utils.array_to_array_info(
+        metrics_savepoint._get_field("wgtfacq_e", dims.EdgeDim, dims.KDim).ndarray
+    )
     vwind_impl_wgt = test_utils.array_to_array_info(metrics_savepoint.vwind_impl_wgt().ndarray)
     hmask_dd3d = test_utils.array_to_array_info(metrics_savepoint.hmask_dd3d().ndarray)
     scalfac_dd3d = test_utils.array_to_array_info(metrics_savepoint.scalfac_dd3d().ndarray)
@@ -133,12 +137,22 @@ def solve_nh_init(
     )
     pos_on_tplane_e_2 = test_utils.array_to_array_info(pos_on_tplane_e_2_field.ndarray)
 
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
     rbf_vec_coeff_e = test_utils.array_to_array_info(
-        interpolation_savepoint.rbf_vec_coeff_e().ndarray
+        interpolation_savepoint._get_field("rbf_vec_coeff_e", dims.E2C2EDim, dims.EdgeDim).ndarray
     )
     e_bln_c_s = test_utils.array_to_array_info(interpolation_savepoint.e_bln_c_s().ndarray)
-    rbf_coeff_1 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v1().ndarray)
-    rbf_coeff_2 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v2().ndarray)
+
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    rbf_vec_coeff_v_array = np.squeeze(
+        interpolation_savepoint.serializer.read(
+            "rbf_vec_coeff_v", interpolation_savepoint.savepoint
+        ).astype(float)
+    )
+    rbf_vec_coeff_v_array = interpolation_savepoint._reduce_to_dim_size(
+        rbf_vec_coeff_v_array, [dims.V2EDim, dims.V2EDim, dims.VertexDim]
+    )
+    rbf_vec_coeff_v = test_utils.array_to_array_info(rbf_vec_coeff_v_array)
     geofac_div = test_utils.array_to_array_info(interpolation_savepoint.geofac_div().ndarray)
     geofac_n2s = test_utils.array_to_array_info(interpolation_savepoint.geofac_n2s().ndarray)
     geofac_grg_x = test_utils.array_to_array_info(interpolation_savepoint.geofac_grg()[0].ndarray)
@@ -161,8 +175,7 @@ def solve_nh_init(
         pos_on_tplane_e_2=pos_on_tplane_e_2,
         rbf_vec_coeff_e=rbf_vec_coeff_e,
         e_bln_c_s=e_bln_c_s,
-        rbf_coeff_1=rbf_coeff_1,
-        rbf_coeff_2=rbf_coeff_2,
+        rbf_vec_coeff_v=rbf_vec_coeff_v,
         geofac_div=geofac_div,
         geofac_n2s=geofac_n2s,
         geofac_grg_x=geofac_grg_x,
@@ -187,7 +200,7 @@ def solve_nh_init(
         theta_ref_me=theta_ref_me,
         ddxn_z_full=ddxn_z_full,
         zdiff_gradp=zdiff_gradp,
-        vertoffset_gradp=vertoffset_gradp,
+        vertidx_gradp=vertidx_gradp,
         pg_edgeidx=pg_edgeidx,
         pg_vertidx=pg_vertidx,
         pg_exdist=pg_exdist,
@@ -227,7 +240,7 @@ def solve_nh_init(
         divdamp_z3=divdamp_z3,
         divdamp_z4=divdamp_z4,
         nflat_gradp=nflat_gradp,
-        backend=wrapper_common.BackendIntEnum.DEFAULT,
+        backend=wrapper_common.BackendIntEnum.GTFN,
     )
 
 
@@ -314,7 +327,10 @@ def test_dycore_wrapper_granule_inputs(
     exner_exfac = test_utils.array_to_array_info(metrics_savepoint.exner_exfac().ndarray)
     exner_ref_mc = test_utils.array_to_array_info(metrics_savepoint.exner_ref_mc().ndarray)
     wgtfac_c = test_utils.array_to_array_info(metrics_savepoint.wgtfac_c().ndarray)
-    wgtfacq_c = test_utils.array_to_array_info(metrics_savepoint.wgtfacq_c().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    wgtfacq_c = test_utils.array_to_array_info(
+        metrics_savepoint._get_field("wgtfacq_c", dims.CellDim, dims.KDim).ndarray
+    )
     inv_ddqz_z_full = test_utils.array_to_array_info(metrics_savepoint.inv_ddqz_z_full().ndarray)
     rho_ref_mc = test_utils.array_to_array_info(metrics_savepoint.rho_ref_mc().ndarray)
     theta_ref_mc = test_utils.array_to_array_info(metrics_savepoint.theta_ref_mc().ndarray)
@@ -330,15 +346,13 @@ def test_dycore_wrapper_granule_inputs(
     theta_ref_me = test_utils.array_to_array_info(metrics_savepoint.theta_ref_me().ndarray)
     ddxn_z_full = test_utils.array_to_array_info(metrics_savepoint.ddxn_z_full().ndarray)
 
-    zdiff_gradp_field = metrics_savepoint._get_field(
-        "zdiff_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim
-    )
-    zdiff_gradp = test_utils.array_to_array_info(zdiff_gradp_field.ndarray)
-
-    vertoffset_gradp_field = metrics_savepoint._get_field(
-        "vertoffset_gradp_dsl", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32
+    zdiff_gradp = test_utils.array_to_array_info(metrics_savepoint.zdiff_gradp().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    vertidx_gradp = test_utils.array_to_array_info(
+        metrics_savepoint._get_field(
+            "vertidx_gradp", dims.EdgeDim, dims.E2CDim, dims.KDim, dtype=gtx.int32
+        ).ndarray
     )
-    vertoffset_gradp = test_utils.array_to_array_info(vertoffset_gradp_field.ndarray)
 
     pg_edgeidx = test_utils.array_to_array_info(metrics_savepoint.pg_edgeidx())
     pg_vertidx = test_utils.array_to_array_info(metrics_savepoint.pg_vertidx())
@@ -346,7 +360,10 @@ def test_dycore_wrapper_granule_inputs(
     ddqz_z_full_e = test_utils.array_to_array_info(metrics_savepoint.ddqz_z_full_e().ndarray)
     ddxt_z_full = test_utils.array_to_array_info(metrics_savepoint.ddxt_z_full().ndarray)
     wgtfac_e = test_utils.array_to_array_info(metrics_savepoint.wgtfac_e().ndarray)
-    wgtfacq_e = test_utils.array_to_array_info(metrics_savepoint.wgtfacq_e().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    wgtfacq_e = test_utils.array_to_array_info(
+        metrics_savepoint._get_field("wgtfacq_e", dims.EdgeDim, dims.KDim).ndarray
+    )
     vwind_impl_wgt = test_utils.array_to_array_info(metrics_savepoint.vwind_impl_wgt().ndarray)
     hmask_dd3d = test_utils.array_to_array_info(metrics_savepoint.hmask_dd3d().ndarray)
     scalfac_dd3d = test_utils.array_to_array_info(metrics_savepoint.scalfac_dd3d().ndarray)
@@ -373,12 +390,21 @@ def test_dycore_wrapper_granule_inputs(
     )
     pos_on_tplane_e_2 = test_utils.array_to_array_info(pos_on_tplane_e_2_field.ndarray)
 
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
     rbf_vec_coeff_e = test_utils.array_to_array_info(
-        interpolation_savepoint.rbf_vec_coeff_e().ndarray
+        interpolation_savepoint._get_field("rbf_vec_coeff_e", dims.E2C2EDim, dims.EdgeDim).ndarray
     )
     e_bln_c_s = test_utils.array_to_array_info(interpolation_savepoint.e_bln_c_s().ndarray)
-    rbf_coeff_1 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v1().ndarray)
-    rbf_coeff_2 = test_utils.array_to_array_info(interpolation_savepoint.rbf_vec_coeff_v2().ndarray)
+    # we need the raw Fortran data instead of the postprocessed GT4Py field, see dycore_wrapper.solve_nh_init
+    rbf_vec_coeff_v_array = np.squeeze(
+        interpolation_savepoint.serializer.read(
+            "rbf_vec_coeff_v", interpolation_savepoint.savepoint
+        ).astype(float)
+    )
+    rbf_vec_coeff_v_array = interpolation_savepoint._reduce_to_dim_size(
+        rbf_vec_coeff_v_array, [dims.V2EDim, dims.V2EDim, dims.VertexDim]
+    )
+    rbf_vec_coeff_v = test_utils.array_to_array_info(rbf_vec_coeff_v_array)
     geofac_div = test_utils.array_to_array_info(interpolation_savepoint.geofac_div().ndarray)
     geofac_n2s = test_utils.array_to_array_info(interpolation_savepoint.geofac_n2s().ndarray)
     geofac_grg_x = test_utils.array_to_array_info(interpolation_savepoint.geofac_grg()[0].ndarray)
@@ -473,7 +499,7 @@ def test_dycore_wrapper_granule_inputs(
         time_extrapolation_parameter_for_exner=metrics_savepoint.exner_exfac(),
         reference_exner_at_cells_on_model_levels=metrics_savepoint.exner_ref_mc(),
         wgtfac_c=metrics_savepoint.wgtfac_c(),
-        wgtfacq_c=metrics_savepoint.wgtfacq_c_dsl(),
+        wgtfacq_c=metrics_savepoint.wgtfacq_c(),
         inv_ddqz_z_full=metrics_savepoint.inv_ddqz_z_full(),
         reference_rho_at_cells_on_model_levels=metrics_savepoint.rho_ref_mc(),
         reference_theta_at_cells_on_model_levels=metrics_savepoint.theta_ref_mc(),
@@ -493,7 +519,7 @@ def test_dycore_wrapper_granule_inputs(
         ddqz_z_full_e=metrics_savepoint.ddqz_z_full_e(),
         ddxt_z_full=metrics_savepoint.ddxt_z_full(),
         wgtfac_e=metrics_savepoint.wgtfac_e(),
-        wgtfacq_e=metrics_savepoint.wgtfacq_e_dsl(),
+        wgtfacq_e=metrics_savepoint.wgtfacq_e(),
         exner_w_implicit_weight_parameter=metrics_savepoint.vwind_impl_wgt(),
         horizontal_mask_for_3d_divdamp=metrics_savepoint.hmask_dd3d(),
         scaling_factor_for_3d_divdamp=metrics_savepoint.scalfac_dd3d(),
@@ -582,8 +608,7 @@ def test_dycore_wrapper_granule_inputs(
             pos_on_tplane_e_2=pos_on_tplane_e_2,
             rbf_vec_coeff_e=rbf_vec_coeff_e,
             e_bln_c_s=e_bln_c_s,
-            rbf_coeff_1=rbf_coeff_1,
-            rbf_coeff_2=rbf_coeff_2,
+            rbf_vec_coeff_v=rbf_vec_coeff_v,
             geofac_div=geofac_div,
             geofac_n2s=geofac_n2s,
             geofac_grg_x=geofac_grg_x,
@@ -608,7 +633,7 @@ def test_dycore_wrapper_granule_inputs(
             theta_ref_me=theta_ref_me,
             ddxn_z_full=ddxn_z_full,
             zdiff_gradp=zdiff_gradp,
-            vertoffset_gradp=vertoffset_gradp,
+            vertidx_gradp=vertidx_gradp,
             pg_edgeidx=pg_edgeidx,
             pg_vertidx=pg_vertidx,
             pg_exdist=pg_exdist,