diff --git a/.github/workflows/icon4py-qa.yml b/.github/workflows/icon4py-qa.yml
index e40d32042b..4dbcdd3541 100644
--- a/.github/workflows/icon4py-qa.yml
+++ b/.github/workflows/icon4py-qa.yml
@@ -29,8 +29,7 @@ jobs:
                     **/base-requirements-dev.txt
                     **/requirements.txt
                     **/requirements-dev.txt
-            - name: Install icon4py-model packages
-              working-directory: model
+            - name: Install all icon4py namespace packages
               run: |
                   python -m pip install --upgrade pip setuptools wheel
                   python -m pip install -r ./requirements-dev.txt
diff --git a/.gitignore b/.gitignore
index 3c5a2f48dc..9cf32d256d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,7 @@
+# nvidia
+*.ncu-rep
+
+
 ### Custom ####
 _build
 _local
@@ -7,6 +11,7 @@ tmp
 testdata
 simple_mesh*.nc
 /prof
+.profile
 
 ### GT4Py ####
 .gt_cache/
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion.py
index 85c6b8c30a..ac9a332aaa 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion.py
@@ -18,7 +18,6 @@
 from enum import Enum
 from typing import Final, Optional
 
-import numpy as np
 from gt4py.next import as_field
 from gt4py.next.common import Dimension
 from gt4py.next.ffront.fbuiltins import Field, int32
@@ -73,8 +72,7 @@
     mo_intp_rbf_rbf_vec_interpol_vertex,
 )
 from icon4py.model.common.states.prognostic_state import PrognosticState
-from icon4py.model.common.model_backend import backend
-
+from icon4py.model.common.settings import xp
 
 """
 Diffusion module ported from ICON mo_nh_diffusion.f90.
@@ -257,9 +255,6 @@ def _validate(self):
             self.apply_to_temperature = False
             self.apply_to_horizontal_wind = False
             self.apply_to_vertical_wind = False
-        else:
-            self.apply_to_temperature = True
-            self.apply_to_horizontal_wind = True
 
         if self.shear_type not in (
             TurbulenceShearForcingType.VERTICAL_OF_HORIZONTAL_WIND,
@@ -474,7 +469,7 @@ def _allocate(*dims: Dimension):
 
         def _index_field(dim: Dimension, size=None):
             size = size if size else self.grid.size[dim]
-            return as_field((dim,), np.arange(size, dtype=int32))
+            return as_field((dim,), xp.arange(size, dtype=int32))
 
         self.diff_multfac_vn = _allocate(KDim)
 
@@ -492,7 +487,7 @@ def _index_field(dim: Dimension, size=None):
         self.horizontal_cell_index = _index_field(CellDim)
         self.horizontal_edge_index = _index_field(EdgeDim)
         self.w_tmp = as_field(
-            (CellDim, KDim), np.zeros((self.grid.num_cells, self.grid.num_levels + 1), dtype=float)
+            (CellDim, KDim), xp.zeros((self.grid.num_cells, self.grid.num_levels + 1), dtype=float)
         )
 
     def initial_run(
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_states.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_states.py
index 96f3be04f9..87db9252fc 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_states.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_states.py
@@ -88,11 +88,11 @@ class DiffusionInterpolationState:
 
     @functools.cached_property
     def geofac_n2s_c(self) -> Field[[CellDim], float]:
-        return as_field((CellDim,), data=self.geofac_n2s.asnumpy()[:, 0])
+        return as_field((CellDim,), data=self.geofac_n2s.ndarray[:, 0])
 
     @functools.cached_property
     def geofac_n2s_nbh(self) -> Field[[CECDim], float]:
-        geofac_nbh_ar = self.geofac_n2s.asnumpy()[:, 1:]
+        geofac_nbh_ar = self.geofac_n2s.ndarray[:, 1:]
         old_shape = geofac_nbh_ar.shape
         return as_field(
             (CECDim,),
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_utils.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_utils.py
index 523b8e951c..94e23cc0aa 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_utils.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_utils.py
@@ -12,7 +12,6 @@
 # SPDX-License-Identifier: GPL-3.0-or-later
 from typing import Tuple
 
-import numpy as np
 from gt4py.next import as_field
 from gt4py.next.common import Dimension, Field, GridType
 from gt4py.next.ffront.decorator import field_operator, program
@@ -20,13 +19,13 @@
 
 from icon4py.model.common.dimension import CellDim, EdgeDim, KDim, VertexDim
 from icon4py.model.common.math.smagorinsky import _en_smag_fac_for_zero_nshift
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend, xp
 
 
 # TODO(Magdalena): fix duplication: duplicated from test testutils/utils.py
 def zero_field(grid, *dims: Dimension, dtype=float):
     shapex = tuple(map(lambda x: grid.size[x], dims))
-    return as_field(dims, np.zeros(shapex, dtype=dtype))
+    return as_field(dims, xp.zeros(shapex, dtype=dtype))
 
 
 @field_operator
@@ -186,8 +185,8 @@ def init_nabla2_factor_in_upper_damping_zone(
         physcial_heights: vector of physical heights [m] of the height levels
     """
     # TODO(Magdalena): fix with as_offset in gt4py
-    heights = physical_heights.asnumpy()
-    buffer = np.zeros(k_size)
+    heights = physical_heights.ndarray
+    buffer = xp.zeros(k_size)
     buffer[1 : nrdmax + 1] = (
         1.0
         / 12.0
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_theta_and_exner.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_theta_and_exner.py
index 541f3a50ef..08bfcc8af1 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_theta_and_exner.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_theta_and_exner.py
@@ -28,7 +28,7 @@
     _update_theta_and_exner,
 )
 from icon4py.model.common.dimension import CECDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_vn.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_vn.py
index 2eda98d777..591d1605d2 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_vn.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_vn.py
@@ -25,7 +25,7 @@
 )
 from icon4py.model.atmosphere.diffusion.stencils.calculate_nabla4 import _calculate_nabla4
 from icon4py.model.common.dimension import ECVDim, EdgeDim, KDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_w_and_compute_horizontal_gradients_for_turbulence.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_w_and_compute_horizontal_gradients_for_turbulence.py
index d4160e460a..8dd62df576 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_w_and_compute_horizontal_gradients_for_turbulence.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_diffusion_to_w_and_compute_horizontal_gradients_for_turbulence.py
@@ -26,7 +26,7 @@
     _calculate_nabla2_for_w,
 )
 from icon4py.model.common.dimension import C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_global_to_vn.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_global_to_vn.py
index 6da0a1fee8..e6e23155e0 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_global_to_vn.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_global_to_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_to_vn.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_to_vn.py
index 43c65ca23a..f3c76a6713 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_to_vn.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_and_nabla4_to_vn.py
@@ -15,7 +15,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, broadcast, int32, maximum
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_vn_in_lateral_boundary.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_vn_in_lateral_boundary.py
index f8750138c5..c46da23269 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_vn_in_lateral_boundary.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_vn_in_lateral_boundary.py
@@ -15,7 +15,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w.py
index 0503f3cbcc..1361c914d2 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2E2CO, C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w_in_upper_damping_layer.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w_in_upper_damping_layer.py
index 1efa53ff05..9c2c854925 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w_in_upper_damping_layer.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/apply_nabla2_to_w_in_upper_damping_layer.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, broadcast, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostic_quantities_for_turbulence.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostic_quantities_for_turbulence.py
index 0aa4a4feb2..000764f073 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostic_quantities_for_turbulence.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostic_quantities_for_turbulence.py
@@ -22,7 +22,7 @@
     _temporary_fields_for_turbulence_diagnostics,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostics_for_turbulence.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostics_for_turbulence.py
index 4a0721e5d7..9dd4f5d307 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostics_for_turbulence.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_diagnostics_for_turbulence.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_enhanced_diffusion_coefficients_for_grid_point_cold_pools.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_enhanced_diffusion_coefficients_for_grid_point_cold_pools.py
index 4aea5bf44e..afb3b45cd3 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_enhanced_diffusion_coefficients_for_grid_point_cold_pools.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_enhanced_diffusion_coefficients_for_grid_point_cold_pools.py
@@ -22,7 +22,7 @@
     _temporary_field_for_grid_point_cold_pools_enhancement,
 )
 from icon4py.model.common.dimension import CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_horizontal_gradients_for_turbulence.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_horizontal_gradients_for_turbulence.py
index 7e9bacab61..84965b3f61 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_horizontal_gradients_for_turbulence.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_horizontal_gradients_for_turbulence.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2E2CO, C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_and_smag_coefficients_for_vn.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_and_smag_coefficients_for_vn.py
index 7435b98555..779027d81f 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_and_smag_coefficients_for_vn.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_and_smag_coefficients_for_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, maximum, minimum, sqrt
 
 from icon4py.model.common.dimension import E2C2V, E2ECV, ECVDim, EdgeDim, KDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_theta.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_theta.py
index 69b7b1779e..00e93d3b6d 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_theta.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_theta.py
@@ -22,7 +22,7 @@
     _calculate_nabla2_of_theta,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_w.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_w.py
index df71c48848..2c04213595 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_w.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_w.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2E2CO, C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_z.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_z.py
index fbe32388e9..2dfdcc64bf 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_z.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_for_z.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_of_theta.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_of_theta.py
index d6ecb4f742..6d1e007041 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_of_theta.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla2_of_theta.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla4.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla4.py
index b1fcd4526c..23c75d63e6 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla4.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/calculate_nabla4.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C2V, E2ECV, ECVDim, EdgeDim, KDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/enhance_diffusion_coefficient_for_grid_point_cold_pools.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/enhance_diffusion_coefficient_for_grid_point_cold_pools.py
index cb1a98a555..39c0431e22 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/enhance_diffusion_coefficient_for_grid_point_cold_pools.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/enhance_diffusion_coefficient_for_grid_point_cold_pools.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32, max_over, maximum
 
 from icon4py.model.common.dimension import E2C, CellDim, E2CDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_field_for_grid_point_cold_pools_enhancement.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_field_for_grid_point_cold_pools_enhancement.py
index 9208fdf09f..05b1dddcf9 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_field_for_grid_point_cold_pools_enhancement.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_field_for_grid_point_cold_pools_enhancement.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum, where
 
 from icon4py.model.common.dimension import C2E2C, C2E2CDim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_fields_for_turbulence_diagnostics.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_fields_for_turbulence_diagnostics.py
index 941ca7ecdf..cbcaaad853 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_fields_for_turbulence_diagnostics.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/temporary_fields_for_turbulence_diagnostics.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/truly_horizontal_diffusion_nabla_of_theta_over_steep_points.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/truly_horizontal_diffusion_nabla_of_theta_over_steep_points.py
index 8699498f4b..431758cf86 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/truly_horizontal_diffusion_nabla_of_theta_over_steep_points.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/truly_horizontal_diffusion_nabla_of_theta_over_steep_points.py
@@ -17,7 +17,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, where
 
 from icon4py.model.common.dimension import C2CEC, C2E2C, CECDim, CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/update_theta_and_exner.py b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/update_theta_and_exner.py
index 0d5e1933a1..85d608f0b0 100644
--- a/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/update_theta_and_exner.py
+++ b/model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/stencils/update_theta_and_exner.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion.py b/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion.py
index 4eb3d39efd..b4e0748473 100644
--- a/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion.py
+++ b/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion.py
@@ -18,7 +18,7 @@
 from icon4py.model.atmosphere.diffusion.diffusion_utils import scale_k
 from icon4py.model.common.grid.horizontal import CellParams, EdgeParams
 from icon4py.model.common.grid.vertical import VerticalModelParams
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.test_utils.datatest_utils import GLOBAL_EXPERIMENT, REGIONAL_EXPERIMENT
 from icon4py.model.common.test_utils.helpers import dallclose
 from icon4py.model.common.test_utils.reference_funcs import enhanced_smagorinski_factor_numpy
diff --git a/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion_utils.py b/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion_utils.py
index 3d56479198..e060efd695 100644
--- a/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion_utils.py
+++ b/model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion_utils.py
@@ -24,7 +24,7 @@
 )
 from icon4py.model.common.dimension import KDim, VertexDim
 from icon4py.model.common.grid.simple import SimpleGrid
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.test_utils.helpers import random_field, zero_field
 
 from .utils import construct_config, diff_multfac_vn_numpy, smag_limit_numpy
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/accumulate_prep_adv_fields.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/accumulate_prep_adv_fields.py
index 1c6baa698d..d278725eb9 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/accumulate_prep_adv_fields.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/accumulate_prep_adv_fields.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_from_data_assimilation.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_from_data_assimilation.py
index 7d7592fa90..e54837a65b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_from_data_assimilation.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_from_data_assimilation.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_to_vn.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_to_vn.py
index fcb381410f..1e10ba0fab 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_to_vn.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_analysis_increments_to_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_w_con_approaching_cfl.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_w_con_approaching_cfl.py
index 6694c74bea..2e2e0875d1 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_w_con_approaching_cfl.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_w_con_approaching_cfl.py
@@ -24,7 +24,7 @@
 )
 
 from icon4py.model.common.dimension import C2E2CO, C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_wn_approaching_cfl.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_wn_approaching_cfl.py
index 38c1948d82..c86137e4ba 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_wn_approaching_cfl.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_extra_diffusion_for_wn_approaching_cfl.py
@@ -38,7 +38,7 @@
     Koff,
     VertexDim,
 )
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_interpolated_horizontal_advection_of_w.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_interpolated_horizontal_advection_of_w.py
index b35f23a745..c49474c7ce 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_interpolated_horizontal_advection_of_w.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_interpolated_horizontal_advection_of_w.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn.py
index fa07e89254..8e2b496de3 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn_by_interpolating_between_time_levels.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn_by_interpolating_between_time_levels.py
index d1d1bf5b79..3d9949f2a8 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn_by_interpolating_between_time_levels.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_temporal_tendencies_to_vn_by_interpolating_between_time_levels.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_vertical_wind_derivative_to_divergence_damping.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_vertical_wind_derivative_to_divergence_damping.py
index 12668ac273..753903513c 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_vertical_wind_derivative_to_divergence_damping.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/add_vertical_wind_derivative_to_divergence_damping.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, broadcast, int32
 
 from icon4py.model.common.dimension import E2C, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_2nd_order_divergence_damping.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_2nd_order_divergence_damping.py
index c5343b5ef2..bdf3b23007 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_2nd_order_divergence_damping.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_2nd_order_divergence_damping.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_4th_order_divergence_damping.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_4th_order_divergence_damping.py
index 5660973b79..ff6edabdee 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_4th_order_divergence_damping.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_4th_order_divergence_damping.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, broadcast, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_hydrostatic_correction_to_horizontal_gradient_of_exner_pressure.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_hydrostatic_correction_to_horizontal_gradient_of_exner_pressure.py
index d8341df31e..f9e51a6ad1 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_hydrostatic_correction_to_horizontal_gradient_of_exner_pressure.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_hydrostatic_correction_to_horizontal_gradient_of_exner_pressure.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32, where
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_rayleigh_damping_mechanism.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_rayleigh_damping_mechanism.py
index 2d83e534e4..30dbb06d0c 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_rayleigh_damping_mechanism.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_rayleigh_damping_mechanism.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, broadcast, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_weighted_2nd_and_4th_order_divergence_damping.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_weighted_2nd_and_4th_order_divergence_damping.py
index 7fe36d9aee..6f4a2d1b35 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_weighted_2nd_and_4th_order_divergence_damping.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/apply_weighted_2nd_and_4th_order_divergence_damping.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, broadcast, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_normal_wind_tendency.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_normal_wind_tendency.py
index 2d04c9ad98..f0d89744e2 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_normal_wind_tendency.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_normal_wind_tendency.py
@@ -28,7 +28,7 @@
     Koff,
     VertexDim,
 )
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_vertical_wind_tendency.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_vertical_wind_tendency.py
index 743cbbec0b..dc901a012d 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_vertical_wind_tendency.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_advective_vertical_wind_tendency.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_airmass.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_airmass.py
index c468bb88d7..fd42c8911a 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_airmass.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_airmass.py
@@ -15,7 +15,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_approx_of_2nd_vertical_derivative_of_exner.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_approx_of_2nd_vertical_derivative_of_exner.py
index 3dfa426a99..c5fbd1ee58 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_approx_of_2nd_vertical_derivative_of_exner.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_approx_of_2nd_vertical_derivative_of_exner.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn.py
index 44ebd75b8f..fbdadd4fd6 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32, neighbor_sum
 
 from icon4py.model.common.dimension import E2C2EO, E2C2EODim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn_and_graddiv_vn_and_vt.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn_and_graddiv_vn_and_vt.py
index 67c8541b39..97db77ca8b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn_and_graddiv_vn_and_vt.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_avg_vn_and_graddiv_vn_and_vt.py
@@ -18,7 +18,7 @@
 from icon4py.model.atmosphere.dycore.compute_avg_vn import _compute_avg_vn
 from icon4py.model.atmosphere.dycore.compute_tangential_wind import _compute_tangential_wind
 from icon4py.model.common.dimension import E2C2EO, E2C2EDim, E2C2EODim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction.py
index 512f9e1755..9ac6fd8271 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w.py
index 2979307c73..5c9b2c58a6 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w_for_lower_boundary.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w_for_lower_boundary.py
index 008984e19d..18fec6aba8 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w_for_lower_boundary.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_contravariant_correction_of_w_for_lower_boundary.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_divergence_of_fluxes_of_rho_and_theta.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_divergence_of_fluxes_of_rho_and_theta.py
index 8d82296c7a..ce1b0d60dd 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_divergence_of_fluxes_of_rho_and_theta.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_divergence_of_fluxes_of_rho_and_theta.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_dwdz_for_divergence_damping.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_dwdz_for_divergence_damping.py
index be487494eb..906c044977 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_dwdz_for_divergence_damping.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_dwdz_for_divergence_damping.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_exner_from_rhotheta.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_exner_from_rhotheta.py
index df3aa474f4..68470e7ac8 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_exner_from_rhotheta.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_exner_from_rhotheta.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, exp, int32, log
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_part_for_rho_and_exner.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_part_for_rho_and_exner.py
index c170d1d899..63992e1995 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_part_for_rho_and_exner.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_part_for_rho_and_exner.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_from_advection_and_vertical_wind_density.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_from_advection_and_vertical_wind_density.py
index 7398bacbda..0f00b88c74 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_from_advection_and_vertical_wind_density.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_from_advection_and_vertical_wind_density.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_speed_and_vertical_wind_times_density.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_speed_and_vertical_wind_times_density.py
index 1f42c1f817..879643a31c 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_speed_and_vertical_wind_times_density.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_explicit_vertical_wind_speed_and_vertical_wind_times_density.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_first_vertical_derivative.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_first_vertical_derivative.py
index c3173c8444..75406eec02 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_first_vertical_derivative.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_first_vertical_derivative.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_graddiv2_of_vn.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_graddiv2_of_vn.py
index 0995c3eeb3..d136ebd8b7 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_graddiv2_of_vn.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_graddiv2_of_vn.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import E2C2EO, E2C2EODim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_of_rho_and_theta.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_of_rho_and_theta.py
index 6bb98fc047..65a2c2029e 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_of_rho_and_theta.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_of_rho_and_theta.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, where
 
 from icon4py.model.common.dimension import E2C, E2EC, CellDim, ECDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_term_for_vertical_velocity.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_term_for_vertical_velocity.py
index d034dd8f78..36231bb463 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_term_for_vertical_velocity.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_advection_term_for_vertical_velocity.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C, E2V, CellDim, EdgeDim, KDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_flat_coordinates.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_flat_coordinates.py
index 2ea32a5685..47659e2d13 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_flat_coordinates.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_flat_coordinates.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_nonflat_coordinates.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_nonflat_coordinates.py
index 32f2837eb1..dc1a41750f 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_nonflat_coordinates.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_exner_pressure_for_nonflat_coordinates.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import E2C, CellDim, E2CDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_extner_pressure_for_multiple_levels.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_extner_pressure_for_multiple_levels.py
index eed15e1672..9ed0c65255 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_extner_pressure_for_multiple_levels.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_gradient_of_extner_pressure_for_multiple_levels.py
@@ -17,7 +17,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C, E2EC, CellDim, ECDim, EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_kinetic_energy.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_kinetic_energy.py
index 1d15f87a23..478acd6d99 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_kinetic_energy.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_horizontal_kinetic_energy.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_hydrostatic_correction_term.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_hydrostatic_correction_term.py
index abe05ee4e7..23789d7d2b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_hydrostatic_correction_term.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_hydrostatic_correction_term.py
@@ -17,7 +17,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import E2C, E2EC, CellDim, ECDim, EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_mass_flux.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_mass_flux.py
index 0237074656..a097675a95 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_mass_flux.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_mass_flux.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_maximum_cfl_and_clip_contravariant_vertical_velocity.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_maximum_cfl_and_clip_contravariant_vertical_velocity.py
index f08ffd62fc..31d21ae006 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_maximum_cfl_and_clip_contravariant_vertical_velocity.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_maximum_cfl_and_clip_contravariant_vertical_velocity.py
@@ -23,7 +23,7 @@
 )
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta.py
index 8f650634dd..ba1602c8be 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta_and_rho_at_ic.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta_and_rho_at_ic.py
index c36f261d75..93b46ce265 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta_and_rho_at_ic.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_pertubation_of_rho_and_theta_and_rho_at_ic.py
@@ -19,7 +19,7 @@
     _compute_pertubation_of_rho_and_theta,
 )
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_results_for_thermodynamic_variables.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_results_for_thermodynamic_variables.py
index 3ab0ae431b..31147db786 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_results_for_thermodynamic_variables.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_results_for_thermodynamic_variables.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_rho_virtual_potential_temperatures_and_pressure_gradient.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_rho_virtual_potential_temperatures_and_pressure_gradient.py
index a44872ad7b..c9322ab14b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_rho_virtual_potential_temperatures_and_pressure_gradient.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_rho_virtual_potential_temperatures_and_pressure_gradient.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_solver_coefficients_matrix.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_solver_coefficients_matrix.py
index 009815fad9..987fc68579 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_solver_coefficients_matrix.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_solver_coefficients_matrix.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_tangential_wind.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_tangential_wind.py
index 0d376c5796..03c67c522e 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_tangential_wind.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_tangential_wind.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import E2C2E, E2C2EDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_theta_and_exner.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_theta_and_exner.py
index f55c48b776..da5df4473c 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_theta_and_exner.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_theta_and_exner.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, exp, int32, log, where
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_virtual_potential_temperatures_and_pressure_gradient.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_virtual_potential_temperatures_and_pressure_gradient.py
index b4f6b9da42..5ff7adee0a 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_virtual_potential_temperatures_and_pressure_gradient.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_virtual_potential_temperatures_and_pressure_gradient.py
@@ -19,7 +19,7 @@
     _interpolate_to_half_levels_vp,
 )
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_vn_on_lateral_boundary.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_vn_on_lateral_boundary.py
index ea5f37dd4a..b798a36e04 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_vn_on_lateral_boundary.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/compute_vn_on_lateral_boundary.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/copy_cell_kdim_field_to_vp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/copy_cell_kdim_field_to_vp.py
index e82b324cec..8af01d5be1 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/copy_cell_kdim_field_to_vp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/copy_cell_kdim_field_to_vp.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/correct_contravariant_vertical_velocity.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/correct_contravariant_vertical_velocity.py
index cde2fc3242..38366f5093 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/correct_contravariant_vertical_velocity.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/correct_contravariant_vertical_velocity.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_at_top.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_at_top.py
index e7da5e5081..8ce8b4818a 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_at_top.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_at_top.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_temporally_exner_pressure.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_temporally_exner_pressure.py
index 6d42ad17c1..881092e38c 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_temporally_exner_pressure.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/extrapolate_temporally_exner_pressure.py
@@ -15,7 +15,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_solve_nonhydro_stencil_39_40.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_solve_nonhydro_stencil_39_40.py
index efad6ee216..a55b426725 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_solve_nonhydro_stencil_39_40.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_solve_nonhydro_stencil_39_40.py
@@ -21,7 +21,7 @@
     _compute_contravariant_correction_of_w_for_lower_boundary,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_15_to_18.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_15_to_18.py
index 8eb521a2a2..deaa23b1b8 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_15_to_18.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_15_to_18.py
@@ -27,7 +27,7 @@
     _interpolate_contravatiant_vertical_verlocity_to_full_levels,
 )
 from icon4py.model.common.dimension import C2E2CODim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_19_to_20.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_19_to_20.py
index cfca1626f2..e7e450e37b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_19_to_20.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_19_to_20.py
@@ -33,7 +33,7 @@
     V2EDim,
     VertexDim,
 )
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_1_to_7.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_1_to_7.py
index 9082d298f6..87303c38ed 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_1_to_7.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_1_to_7.py
@@ -33,7 +33,7 @@
     _mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl,
 )
 from icon4py.model.common.dimension import CellDim, E2C2EDim, EdgeDim, KDim, V2CDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_13.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_13.py
index 3167a6c8e9..bb8da96ad8 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_13.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_13.py
@@ -26,7 +26,7 @@
     _set_cell_kdim_field_to_zero_vp,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_14.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_14.py
index 695b77e57a..ea0074ed2f 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_14.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/fused_velocity_advection_stencil_8_to_14.py
@@ -29,7 +29,7 @@
     _set_cell_kdim_field_to_zero_vp,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_contravatiant_vertical_verlocity_to_full_levels.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_contravatiant_vertical_verlocity_to_full_levels.py
index a4d496998f..29a8527a4b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_contravatiant_vertical_verlocity_to_full_levels.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_contravatiant_vertical_verlocity_to_full_levels.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_cell_center.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_cell_center.py
index 17b1649386..1bf776d48b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_cell_center.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_cell_center.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2CE, C2E, C2EDim, CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_half_levels_vp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_half_levels_vp.py
index 4c1332af94..350f006a0f 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_half_levels_vp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_half_levels_vp.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_surface.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_surface.py
index 5467242f68..8cfabb61ed 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_surface.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_to_surface.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_and_vt_to_ie_and_compute_ekin_on_edges.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_and_vt_to_ie_and_compute_ekin_on_edges.py
index 895a7d9eb1..b33ed00b50 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_and_vt_to_ie_and_compute_ekin_on_edges.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_and_vt_to_ie_and_compute_ekin_on_edges.py
@@ -20,7 +20,7 @@
 )
 from icon4py.model.atmosphere.dycore.interpolate_vt_to_ie import _interpolate_vt_to_ie
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_to_ie_and_compute_ekin_on_edges.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_to_ie_and_compute_ekin_on_edges.py
index a1ff78abec..47af3271f7 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_to_ie_and_compute_ekin_on_edges.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vn_to_ie_and_compute_ekin_on_edges.py
@@ -19,7 +19,7 @@
     _compute_horizontal_kinetic_energy,
 )
 from icon4py.model.common.dimension import EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vt_to_ie.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vt_to_ie.py
index 84737c8f66..3a930515cc 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vt_to_ie.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/interpolate_vt_to_ie.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl.py
index daf11cef73..d49f92afb0 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import V2C, CellDim, KDim, V2CDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_divrot_rot_vertex_ri_dsl.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_divrot_rot_vertex_ri_dsl.py
index 617ff9f29e..5ea61aa6ab 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_divrot_rot_vertex_ri_dsl.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_divrot_rot_vertex_ri_dsl.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import V2E, EdgeDim, KDim, V2EDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_gradients_grad_green_gauss_cell_dsl.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_gradients_grad_green_gauss_cell_dsl.py
index 47724a039a..dafe3e5318 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_gradients_grad_green_gauss_cell_dsl.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_math_gradients_grad_green_gauss_cell_dsl.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32, neighbor_sum
 
 from icon4py.model.common.dimension import C2E2CO, C2E2CODim, CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_solve_nonhydro_stencil_51.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_solve_nonhydro_stencil_51.py
index 6f921b7601..3b92096990 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_solve_nonhydro_stencil_51.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/mo_solve_nonhydro_stencil_51.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 
 
 @field_operator
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro.py
index db3787c9a5..4b802d3e85 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro.py
@@ -168,7 +168,7 @@
 from icon4py.model.common.grid.vertical import VerticalModelParams
 from icon4py.model.common.math.smagorinsky import en_smag_fac_for_zero_nshift
 from icon4py.model.common.states.prognostic_state import PrognosticState
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 
 
 # flake8: noqa
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro_program.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro_program.py
index c70244c11b..c51edccf87 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro_program.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/nh_solve/solve_nonhydro_program.py
@@ -82,7 +82,7 @@
 from icon4py.model.atmosphere.dycore.update_densety_exener_wind import _update_densety_exener_wind
 from icon4py.model.atmosphere.dycore.update_wind import _update_wind
 from icon4py.model.common.dimension import CEDim, CellDim, ECDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 
 
 # TODO: abishekg7 move this to tests
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_vp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_vp.py
index fca20c21e4..ab45b2e8ec 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_vp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_vp.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, broadcast, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_wp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_wp.py
index 8837caa0da..fb053afad2 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_wp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_cell_kdim_field_to_zero_wp.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, broadcast, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_lower_boundary_condition_for_w_and_contravariant_correction.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_lower_boundary_condition_for_w_and_contravariant_correction.py
index 9c89fe5fcb..d4bf7addfc 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_lower_boundary_condition_for_w_and_contravariant_correction.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_lower_boundary_condition_for_w_and_contravariant_correction.py
@@ -19,7 +19,7 @@
     _set_cell_kdim_field_to_zero_wp,
 )
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_theta_v_prime_ic_at_lower_boundary.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_theta_v_prime_ic_at_lower_boundary.py
index 8eabaf9de6..89bd442968 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_theta_v_prime_ic_at_lower_boundary.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_theta_v_prime_ic_at_lower_boundary.py
@@ -17,7 +17,7 @@
 
 from icon4py.model.atmosphere.dycore.interpolate_to_surface import _interpolate_to_surface
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_index_to_zero_vp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_index_to_zero_vp.py
index dbdd114ba0..d4c8f0b658 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_index_to_zero_vp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_index_to_zero_vp.py
@@ -19,7 +19,7 @@
     _set_cell_kdim_field_to_zero_vp,
 )
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_vp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_vp.py
index 6942d1f006..a315f4bac2 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_vp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_vp.py
@@ -19,7 +19,7 @@
     _set_cell_kdim_field_to_zero_vp,
 )
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_wp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_wp.py
index 29b3635df6..990cc677d1 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_wp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_cell_kdim_fields_to_zero_wp.py
@@ -19,7 +19,7 @@
     _set_cell_kdim_field_to_zero_wp,
 )
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_edge_kdim_fields_to_zero_wp.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_edge_kdim_fields_to_zero_wp.py
index 888de75f9b..0094f5a259 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_edge_kdim_fields_to_zero_wp.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/set_two_edge_kdim_fields_to_zero_wp.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, broadcast, int32
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_back_substitution.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_back_substitution.py
index 04d5f40265..bae3536819 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_back_substitution.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_back_substitution.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_forward_sweep.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_forward_sweep.py
index fe8e6e14f3..e43b221783 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_forward_sweep.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_tridiagonal_matrix_for_w_forward_sweep.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/state_utils/utils.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/state_utils/utils.py
index ee5acd308a..5aa3b973f5 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/state_utils/utils.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/state_utils/utils.py
@@ -23,7 +23,7 @@
 )
 
 from icon4py.model.common.dimension import EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_densety_exener_wind.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_densety_exener_wind.py
index 4579dd3df9..ab9c46d617 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_densety_exener_wind.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_densety_exener_wind.py
@@ -17,7 +17,7 @@
 
 from icon4py.model.atmosphere.dycore.update_wind import _update_wind
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_dynamical_exner_time_increment.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_dynamical_exner_time_increment.py
index 0860757117..02d61d27ba 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_dynamical_exner_time_increment.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_dynamical_exner_time_increment.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_flux_weighted.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_flux_weighted.py
index c8bcfc7df3..b5dee86658 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_flux_weighted.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_flux_weighted.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, astype, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_volume_flux.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_volume_flux.py
index 4090a78e97..4a67a0f867 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_volume_flux.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_mass_volume_flux.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_theta_v.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_theta_v.py
index f121943310..fe661422db 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_theta_v.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_theta_v.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32, where
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_wind.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_wind.py
index fcd1508753..e3c64294fe 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_wind.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/update_wind.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32
 
 from icon4py.model.common.dimension import CellDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/velocity/velocity_advection_program.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/velocity/velocity_advection_program.py
index 862acfd611..386a365066 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/velocity/velocity_advection_program.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/velocity/velocity_advection_program.py
@@ -42,7 +42,7 @@
     _set_cell_kdim_field_to_zero_vp,
 )
 from icon4py.model.common.dimension import CEDim, CellDim, EdgeDim, KDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 
 
 @field_operator
diff --git a/model/atmosphere/dycore/tests/dycore_tests/test_solve_nonhydro.py b/model/atmosphere/dycore/tests/dycore_tests/test_solve_nonhydro.py
index 4cc4454574..a73fca0ff6 100644
--- a/model/atmosphere/dycore/tests/dycore_tests/test_solve_nonhydro.py
+++ b/model/atmosphere/dycore/tests/dycore_tests/test_solve_nonhydro.py
@@ -39,7 +39,7 @@
 )
 from icon4py.model.common.grid.vertical import VerticalModelParams
 from icon4py.model.common.math.smagorinsky import en_smag_fac_for_zero_nshift
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.states.prognostic_state import PrognosticState
 from icon4py.model.common.test_utils.datatest_utils import (
     GLOBAL_EXPERIMENT,
diff --git a/model/atmosphere/dycore/tests/dycore_tests/test_utils.py b/model/atmosphere/dycore/tests/dycore_tests/test_utils.py
index 276e3a0adf..b089faaed0 100644
--- a/model/atmosphere/dycore/tests/dycore_tests/test_utils.py
+++ b/model/atmosphere/dycore/tests/dycore_tests/test_utils.py
@@ -22,7 +22,7 @@
 from icon4py.model.common import constants
 from icon4py.model.common.dimension import KDim
 from icon4py.model.common.grid.simple import SimpleGrid
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.test_utils.helpers import dallclose, random_field, zero_field
 
 
diff --git a/model/common/src/icon4py/model/common/config.py b/model/common/src/icon4py/model/common/config.py
new file mode 100644
index 0000000000..8c9f0eb95c
--- /dev/null
+++ b/model/common/src/icon4py/model/common/config.py
@@ -0,0 +1,71 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# This file is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or any later
+# version. See the LICENSE.txt file at the top-level directory of this
+# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+# Assuming this code is in a module called icon4py_config.py
+import dataclasses
+import os
+from enum import Enum
+from functools import cached_property
+
+import numpy as np
+from gt4py.next.program_processors.runners.gtfn import (
+    run_gtfn_cached,
+    run_gtfn_gpu_cached,
+)
+from gt4py.next.program_processors.runners.roundtrip import backend as run_roundtrip
+
+
+class Device(Enum):
+    CPU = "CPU"
+    GPU = "GPU"
+
+
+class GT4PyBackend(Enum):
+    CPU = "run_gtfn_cached"
+    GPU = "run_gtfn_gpu_cached"
+    ROUNDTRIP = "run_roundtrip"
+
+
+@dataclasses.dataclass
+class Icon4PyConfig:
+    @cached_property
+    def icon4py_backend(self):
+        return os.environ.get("ICON4PY_BACKEND", "CPU")
+
+    @cached_property
+    def array_ns(self):
+        if self.device == Device.GPU:
+            import cupy as cp  # type: ignore[import-untyped]
+
+            return cp
+        else:
+            return np
+
+    @cached_property
+    def gt4py_runner(self):
+        backend_map = {
+            GT4PyBackend.CPU.name: run_gtfn_cached,
+            GT4PyBackend.GPU.name: run_gtfn_gpu_cached,
+            GT4PyBackend.ROUNDTRIP.name: run_roundtrip,
+        }
+        return backend_map[self.icon4py_backend]
+
+    @cached_property
+    def device(self):
+        device_map = {
+            GT4PyBackend.CPU.name: Device.CPU,
+            GT4PyBackend.GPU.name: Device.GPU,
+            GT4PyBackend.ROUNDTRIP.name: Device.CPU,
+        }
+        device = device_map[self.icon4py_backend]
+        return device
diff --git a/model/common/src/icon4py/model/common/grid/vertical.py b/model/common/src/icon4py/model/common/grid/vertical.py
index f8434dd796..8609f282e0 100644
--- a/model/common/src/icon4py/model/common/grid/vertical.py
+++ b/model/common/src/icon4py/model/common/grid/vertical.py
@@ -19,6 +19,7 @@
 from gt4py.next.ffront.fbuiltins import int32
 
 from icon4py.model.common.dimension import KDim
+from icon4py.model.common.settings import xp
 
 
 log = logging.getLogger(__name__)
@@ -49,7 +50,7 @@ class VerticalModelParams:
     nflat_gradp: Final[int32] = None
 
     def __post_init__(self):
-        vct_a_array = self.vct_a.asnumpy()
+        vct_a_array = self.vct_a.ndarray
         object.__setattr__(
             self,
             "index_of_damping_layer",
@@ -79,11 +80,11 @@ def _determine_kstart_moist(
     ) -> int32:
         n_levels = vct_a.shape[0]
         interface_height = 0.5 * (vct_a[: n_levels - 1 - nshift_total] + vct_a[1 + nshift_total :])
-        return int32(np.min(np.where(interface_height < top_moist_threshold)))
+        return int32(xp.min(xp.where(interface_height < top_moist_threshold)[0]).item())
 
     @classmethod
     def _determine_damping_height_index(cls, vct_a: np.ndarray, damping_height: float):
-        return int32(np.argmax(np.where(vct_a >= damping_height)))
+        return int32(xp.argmax(xp.where(vct_a >= damping_height)[0]).item())
 
     @property
     def physical_heights(self) -> Field[[KDim], float]:
diff --git a/model/common/src/icon4py/model/common/interpolation/stencils/mo_intp_rbf_rbf_vec_interpol_vertex.py b/model/common/src/icon4py/model/common/interpolation/stencils/mo_intp_rbf_rbf_vec_interpol_vertex.py
index 979d96c70e..2d867d0a93 100644
--- a/model/common/src/icon4py/model/common/interpolation/stencils/mo_intp_rbf_rbf_vec_interpol_vertex.py
+++ b/model/common/src/icon4py/model/common/interpolation/stencils/mo_intp_rbf_rbf_vec_interpol_vertex.py
@@ -16,7 +16,7 @@
 from gt4py.next.ffront.fbuiltins import Field, int32, neighbor_sum
 
 from icon4py.model.common.dimension import V2E, EdgeDim, KDim, V2EDim, VertexDim
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 from icon4py.model.common.type_alias import wpfloat
 
 
diff --git a/model/common/src/icon4py/model/common/math/smagorinsky.py b/model/common/src/icon4py/model/common/math/smagorinsky.py
index b82366302c..9e6765b1f2 100644
--- a/model/common/src/icon4py/model/common/math/smagorinsky.py
+++ b/model/common/src/icon4py/model/common/math/smagorinsky.py
@@ -15,7 +15,7 @@
 from gt4py.next.ffront.fbuiltins import broadcast, maximum, minimum
 
 from icon4py.model.common.dimension import KDim, Koff
-from icon4py.model.common.model_backend import backend
+from icon4py.model.common.settings import backend
 
 
 @field_operator
diff --git a/model/common/src/icon4py/model/common/model_backend.py b/model/common/src/icon4py/model/common/settings.py
similarity index 65%
rename from model/common/src/icon4py/model/common/model_backend.py
rename to model/common/src/icon4py/model/common/settings.py
index 3f2b518a9e..0f5d1a82c0 100644
--- a/model/common/src/icon4py/model/common/model_backend.py
+++ b/model/common/src/icon4py/model/common/settings.py
@@ -10,16 +10,10 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
-from gt4py.next.program_processors.runners.gtfn import (
-    run_gtfn,
-    run_gtfn_cached,
-    run_gtfn_gpu,
-    run_gtfn_imperative,
-)
+from icon4py.model.common.config import Icon4PyConfig
 
 
-cached_backend = run_gtfn_cached
-compiled_backend = run_gtfn
-imperative_backend = run_gtfn_imperative
-gpu_backend = run_gtfn_gpu
-backend = cached_backend
+config = Icon4PyConfig()
+backend = config.gt4py_runner
+xp = config.array_ns
+device = config.device
diff --git a/model/common/src/icon4py/model/common/states/prognostic_state.py b/model/common/src/icon4py/model/common/states/prognostic_state.py
index 9a0719faf5..a248211314 100644
--- a/model/common/src/icon4py/model/common/states/prognostic_state.py
+++ b/model/common/src/icon4py/model/common/states/prognostic_state.py
@@ -36,4 +36,4 @@ class PrognosticState:
 
     @property
     def w_1(self) -> Field[[CellDim], float]:
-        return as_field((CellDim,), self.w.asnumpy()[:, 0])
+        return as_field((CellDim,), self.w.ndarray[:, 0])
diff --git a/model/common/src/icon4py/model/common/test_utils/helpers.py b/model/common/src/icon4py/model/common/test_utils/helpers.py
index 25cf8357da..7de2eea602 100644
--- a/model/common/src/icon4py/model/common/test_utils/helpers.py
+++ b/model/common/src/icon4py/model/common/test_utils/helpers.py
@@ -117,7 +117,7 @@ def constant_field(
 
 def as_1D_sparse_field(field: gt_common.Field, target_dim: gt_common.Dimension) -> gt_common.Field:
     """Convert a 2D sparse field to a 1D flattened (Felix-style) sparse field."""
-    buffer = field.asnumpy()
+    buffer = field.ndarray
     return numpy_to_1D_sparse_field(buffer, target_dim)
 
 
@@ -131,7 +131,7 @@ def numpy_to_1D_sparse_field(field: np.ndarray, dim: gt_common.Dimension) -> gt_
 
 def flatten_first_two_dims(*dims: gt_common.Dimension, field: gt_common.Field) -> gt_common.Field:
     """Convert a n-D sparse field to a (n-1)-D flattened (Felix-style) sparse field."""
-    buffer = field.asnumpy()
+    buffer = field.ndarray
     old_shape = buffer.shape
     assert len(old_shape) >= 2
     flattened_size = old_shape[0] * old_shape[1]
diff --git a/tools/README.md b/tools/README.md
index 9a311082eb..d6288ad2f7 100644
--- a/tools/README.md
+++ b/tools/README.md
@@ -314,7 +314,7 @@ OUTPUT_FILEPATH   A path to the output Fortran source file to be generated.
 
 `py2fgen` is a command-line interface (CLI) tool designed to generate C and Fortran 90 (F90) wrappers, as well as a C library, for embedding a Python module into C and Fortran applications. This tool facilitates the embedding of Python code into Fortran programs by utilizing the [`CFFI`](https://cffi.readthedocs.io/en/latest/embedding) library. `CFFI` instantiates a Python interpreter to execute Python code which is "frozen" into the dynamic library generated by `CFFI`.
 
-**Note:** `py2fgen` is currently in an **experimental** stage. Simple examples have been successfully tested within Fortran code, but it is not yet production-ready. Further testing with complex Python code is necessary. The performance implications of invoking a Python interpreter from within Fortran are also yet to be fully understood.
+**Note:** `py2fgen` is currently in an **proof of concept** stage. Although embedding of more complex functions in Fortran such as the diffusion granule has been done in ICON, the performance implications of invoking a Python interpreter from within Fortran yet to be fully understood.
 
 ## Usage
 
@@ -324,16 +324,26 @@ OUTPUT_FILEPATH   A path to the output Fortran source file to be generated.
 py2fgen [OPTIONS] MODULE_IMPORT_PATH FUNCTION_NAME
 
 Arguments:
-  MODULE_IMPORT_PATH  The Python module containing the function to embed.
-  FUNCTION_NAME       The function within the module to embed.
-
+  MODULE_IMPORT_PATH  The Python module import path to the module where the functions to be embedded are defined.
+  FUNCTIONS           A comma-separated list of functions to be embedded in the case of multiple, otherwise just the function name.
+  PLUGIN_NAME         The name of the plugin used for creating the shared library and bindings.
 Options:
-  -b, --build-path PATH     Specify the directory for generated code and compiled libraries.
-  -d, --debug-mode          Enable debug mode to print additional runtime information.
-  -g, --gt4py-backend TYPE  Set the gt4py backend to use (options: CPU, GPU, ROUNDTRIP).
-  --help                    Display the help message and exit.
+  -o, --output-path PATH          Specify the directory for generated code and
+                                  compiled libraries.
+  -b, --backend [CPU|GPU|ROUNDTRIP]
+                                  Set the backend to use, thereby unpacking
+                                  Fortran pointers into NumPy or CuPy arrays
+                                  respectively.
+  -d, --debug-mode                Enable debug mode to log additional Python
+                                  runtime information.
+  --help                          Show this message and exit.
 ```
 
+## Important Environment Variables
+
+When embedding granules it may be necessary to use an ICON grid file, as is the case in the diffusion granule wrapper.
+The granule expects an `ICON_GRID_LOC` environment variable with the path to the folder holding the grid netcdf file.
+
 ### Example
 
 To create a Fortran interface along with the dynamic library for a Python function named `square` within the module `example.functions`, execute:
@@ -342,6 +352,12 @@ To create a Fortran interface along with the dynamic library for a Python functi
 py2fgen example.functions square
 ```
 
+It is also possible to generate bindings for more than one function at a time by using a comma-separated list of function names:
+
+```bash
+py2fgen example.functions square,square2
+```
+
 `py2fgen` can accept two types of functions:
 
 - **Simple Function:** Any Python function can be exposed.
@@ -368,15 +384,9 @@ Examples can be found under `tools/tests/py2fgen/fortran_samples`.
 
 ## Compilation
 
-Compiling your Fortran driver code requires a Fortran compiler, such as `gfortran`. Follow these steps:
+Compiling your Fortran driver code requires a Fortran compiler, such as `gfortran` or `nvfortran`. Follow these steps:
 
-1. Compile (without linking) .f90 interface:
-
-```bash
-gfortran -c <function_name>_plugin.f90
-```
-
-2. Compile the Fortran driver code along with the Fortran interface and dynamic library:
+1. Compile and link the Fortran driver code along with the Fortran interface and dynamic library:
 
 ```bash
 gfortran -I. -Wl,-rpath=. -L. <function_name>_plugin.f90 <fortran_driver>.f90 -l<function_name>_plugin -o <executable_name>
@@ -385,3 +395,27 @@ gfortran -I. -Wl,-rpath=. -L. <function_name>_plugin.f90 <fortran_driver>.f90 -l
 Replace `<function_name>`, `<fortran_driver>`, and `<executable_name>` with the appropriate names for your project.
 
 **Note:** When executing the compiled binary make sure that you have sourced a Python virtual environment where all required dependencies to run the embedded Python code are present.
+
+## Error handling
+
+All generated Python wrapper code is wrapped in a `try: ... except: ...` block, with the wrapper function returning an error code `1` if an Exception ocurred
+and `0` otherwise. In case of an exception the error message is written to the `py2f_cffi.log` file, which is located in the same directory as the generated bindings.
+This means that on the Fortran side we can handle errors gracefully as follows:
+
+```Fortran
+integer(c_int) :: rc
+real(c_double), dimension(:, :), allocatable :: input, result
+
+call square(input, result, rc)
+
+! handle the Python error here
+ if (rc /= 0) then
+     print *, "Python failed with exit code = ", rc
+     call exit(1)
+ end if
+```
+
+## Other requirements
+
+- Embedded Python functions must have type hints for all function parameters, as these are used to derive the corresponding C and Fortran types.
+- Embedded Python functions are assumed to modify function parameters in-place. Explicitly returning anything is currently not supported.
diff --git a/tools/pyproject.toml b/tools/pyproject.toml
index 2f453c0b2c..8bb7d3df53 100644
--- a/tools/pyproject.toml
+++ b/tools/pyproject.toml
@@ -36,8 +36,8 @@ readme = 'README.md'
 requires-python = '>=3.10'
 
 [project.optional-dependencies]
-all = ["icon4pytools[py2f]"]
-py2f = ['cffi>=1.5']
+cupy11 = ['cupy-cuda11x']
+cupy12 = ['cupy-cuda12x']
 
 [project.scripts]
 f2ser = 'icon4pytools.f2ser.cli:main'
diff --git a/tools/src/icon4pytools/py2fgen/cffi.py b/tools/src/icon4pytools/py2fgen/cffi.py
deleted file mode 100644
index cb364e31ad..0000000000
--- a/tools/src/icon4pytools/py2fgen/cffi.py
+++ /dev/null
@@ -1,116 +0,0 @@
-# ICON4Py - ICON inspired code in Python and GT4Py
-#
-# Copyright (c) 2022, ETH Zurich and MeteoSwiss
-# All rights reserved.
-#
-# This file is free software: you can redistribute it and/or modify it under
-# the terms of the GNU General Public License as published by the
-# Free Software Foundation, either version 3 of the License, or any later
-# version. See the LICENSE.txt file at the top-level directory of this
-# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
-#
-# SPDX-License-Identifier: GPL-3.0-or-later
-import logging
-from pathlib import Path
-
-import cffi
-import numpy as np
-from cffi import FFI
-
-from icon4pytools.common.logger import setup_logger
-
-
-ffi = FFI()
-
-logger = setup_logger(__name__)
-
-
-def unpack(ptr, *sizes: int) -> np.ndarray:
-    """
-    Unpacks an n-dimensional Fortran (column-major) array into a NumPy array (row-major).
-
-    This function is particularly useful when interfacing with C code that returns arrays
-    in column-major order (Fortran order), but you want to work with the array in row-major
-    order (NumPy's default order).
-
-    Args:
-        ptr: A C pointer to the field.
-        *sizes: Size arguments representing the length of each dimension of the array in
-            row-major order. The length of this argument must match the number of dimensions
-            of the array.
-
-    Returns:
-        A NumPy array with shape specified by the sizes and dtype determined by the ctype
-        (C data type) of the pointer.
-    """
-    length = np.prod(sizes)
-    c_type = ffi.getctype(ffi.typeof(ptr).item)
-
-    # special casing different types
-    if c_type == "int":
-        dtype = np.int32
-    else:
-        dtype = np.dtype(c_type)  # type: ignore
-
-    arr = np.frombuffer(  # type: ignore
-        ffi.buffer(ptr, length * ffi.sizeof(c_type)),
-        dtype=dtype,
-        count=-1,
-        offset=0,
-    ).reshape(sizes)
-    return arr
-
-
-def generate_and_compile_cffi_plugin(
-    plugin_name: str, c_header: str, python_wrapper: str, build_path: Path
-) -> None:
-    """
-    Create and compile a CFFI plugin.
-
-    This function generates a C shared library and Fortran interface for Python functions
-    to be exposed in the {plugin_name} module. It creates a linkable C library named
-    'lib{plugin_name}.so' in the specified build directory.
-
-    Args:
-        plugin_name: Name of the plugin.
-        c_header: C header signatures for the Python functions.
-        python_wrapper: Python code wrapping the original function to be exposed.
-        build_path: Path to the build directory.
-    """
-    try:
-        header_file_path = write_c_header(build_path, plugin_name, c_header)
-        compile_cffi_plugin(
-            builder=configure_cffi_builder(c_header, plugin_name, header_file_path),
-            python_wrapper=python_wrapper,
-            build_path=str(build_path),
-            plugin_name=plugin_name,
-        )
-    except Exception as e:
-        logging.error(f"Error generating and compiling CFFI plugin: {e}")
-        raise
-
-
-def write_c_header(build_path: Path, plugin_name: str, c_header: str) -> Path:
-    """Write the C header file to the specified path."""
-    c_header_file = plugin_name + ".h"
-    header_file_path = build_path / c_header_file
-    with open(header_file_path, "w") as f:
-        f.write(c_header)
-    return header_file_path
-
-
-def configure_cffi_builder(c_header: str, plugin_name: str, header_file_path: Path) -> cffi.FFI:
-    """Configure and returns a CFFI FFI builder instance."""
-    builder = cffi.FFI()
-    builder.embedding_api(c_header)
-    builder.set_source(plugin_name, f'#include "{header_file_path.name}"')
-    return builder
-
-
-def compile_cffi_plugin(
-    builder: cffi.FFI, python_wrapper: str, build_path: str, plugin_name: str
-) -> None:
-    """Compile the CFFI plugin with the given configuration."""
-    logger.info("Compiling CFFI dynamic library...")
-    builder.embedding_init_code(python_wrapper)
-    builder.compile(tmpdir=build_path, target=f"lib{plugin_name}.*", verbose=True)
diff --git a/tools/src/icon4pytools/py2fgen/cli.py b/tools/src/icon4pytools/py2fgen/cli.py
index 01e15520b4..99db23a5c8 100644
--- a/tools/src/icon4pytools/py2fgen/cli.py
+++ b/tools/src/icon4pytools/py2fgen/cli.py
@@ -10,20 +10,24 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
-
 import pathlib
 
 import click
+from icon4py.model.common.config import GT4PyBackend
 
 from icon4pytools.icon4pygen.bindings.utils import write_string
-from icon4pytools.py2fgen.cffi import generate_and_compile_cffi_plugin
 from icon4pytools.py2fgen.generate import (
     generate_c_header,
     generate_f90_interface,
     generate_python_wrapper,
 )
 from icon4pytools.py2fgen.parsing import parse
-from icon4pytools.py2fgen.utils import Backend
+from icon4pytools.py2fgen.plugin import generate_and_compile_cffi_plugin
+
+
+def parse_comma_separated_list(ctx, param, value) -> list[str]:
+    # Splits the input string by commas and strips any leading/trailing whitespace from the strings
+    return [item.strip() for item in value.split(",")]
 
 
 @click.command("py2fgen")
@@ -31,49 +35,48 @@
     "module_import_path",
     type=str,
 )
-@click.argument("function_name", type=str)
+@click.argument("functions", type=str, callback=parse_comma_separated_list)
+@click.argument("plugin_name", type=str)
 @click.option(
-    "--build-path",
-    "-b",
+    "--output-path",
+    "-o",
     type=click.Path(dir_okay=True, resolve_path=True, path_type=pathlib.Path),
     default=".",
     help="Specify the directory for generated code and compiled libraries.",
 )
+@click.option(
+    "--backend",
+    "-b",
+    type=click.Choice([e.name for e in GT4PyBackend], case_sensitive=False),
+    default="CPU",
+    help="Set the backend to use, thereby unpacking Fortran pointers into NumPy or CuPy arrays respectively.",
+)
 @click.option(
     "--debug-mode",
     "-d",
     is_flag=True,
-    help="Enable debug mode to print additional runtime information.",
-)
-@click.option(
-    "--gt4py-backend",
-    "-g",
-    type=click.Choice([e.name for e in Backend], case_sensitive=False),
-    default="ROUNDTRIP",
-    help="Set the gt4py backend to use.",
+    help="Enable debug mode to log additional Python runtime information.",
 )
 def main(
     module_import_path: str,
-    function_name: str,
-    build_path: pathlib.Path,
+    functions: list[str],
+    plugin_name: str,
+    output_path: pathlib.Path,
     debug_mode: bool,
-    gt4py_backend: str,
+    backend: str,
 ) -> None:
     """Generate C and F90 wrappers and C library for embedding a Python module in C and Fortran."""
-    backend = Backend[gt4py_backend]
-    build_path.mkdir(exist_ok=True, parents=True)
+    output_path.mkdir(exist_ok=True, parents=True)
 
-    plugin = parse(module_import_path, function_name)
+    plugin = parse(module_import_path, functions, plugin_name)
 
     c_header = generate_c_header(plugin)
-    python_wrapper = generate_python_wrapper(plugin, backend.value, debug_mode)
+    python_wrapper = generate_python_wrapper(plugin, backend, debug_mode)
     f90_interface = generate_f90_interface(plugin)
 
-    generate_and_compile_cffi_plugin(plugin.plugin_name, c_header, python_wrapper, build_path)
-    write_string(f90_interface, build_path, f"{plugin.plugin_name}.f90")
-
-    if debug_mode:
-        write_string(python_wrapper, build_path, f"{plugin.plugin_name}.py")
+    generate_and_compile_cffi_plugin(plugin.plugin_name, c_header, python_wrapper, output_path)
+    write_string(f90_interface, output_path, f"{plugin.plugin_name}.f90")
+    write_string(python_wrapper, output_path, f"{plugin.plugin_name}.py")
 
 
 if __name__ == "__main__":
diff --git a/tools/src/icon4pytools/py2fgen/generate.py b/tools/src/icon4pytools/py2fgen/generate.py
index 143a4f8c9a..9812d8a3f8 100644
--- a/tools/src/icon4pytools/py2fgen/generate.py
+++ b/tools/src/icon4pytools/py2fgen/generate.py
@@ -20,6 +20,7 @@
 from icon4pytools.py2fgen.template import (
     CffiPlugin,
     CHeaderGenerator,
+    F90Interface,
     F90InterfaceGenerator,
     PythonWrapper,
     PythonWrapperGenerator,
@@ -44,15 +45,13 @@ def generate_c_header(plugin: CffiPlugin) -> str:
     return codegen.format_source("cpp", generated_code, style="LLVM")
 
 
-def generate_python_wrapper(
-    plugin: CffiPlugin, gt4py_backend: Optional[str], debug_mode: bool
-) -> str:
+def generate_python_wrapper(plugin: CffiPlugin, backend: Optional[str], debug_mode: bool) -> str:
     """
     Generate Python wrapper code.
 
     Args:
         plugin: The CffiPlugin instance containing information for code generation.
-        gt4py_backend: Optional gt4py backend specification.
+        backend: Optional gt4py backend specification.
         debug_mode: Flag indicating if debug mode is enabled.
 
     Returns:
@@ -62,9 +61,9 @@ def generate_python_wrapper(
     python_wrapper = PythonWrapper(
         module_name=plugin.module_name,
         plugin_name=plugin.plugin_name,
-        function=plugin.function,
+        functions=plugin.functions,
         imports=plugin.imports,
-        gt4py_backend=gt4py_backend,
+        backend=backend,
         debug_mode=debug_mode,
     )
 
@@ -80,5 +79,5 @@ def generate_f90_interface(plugin: CffiPlugin) -> str:
         plugin: The CffiPlugin instance containing information for code generation.
     """
     logger.info("Generating Fortran interface...")
-    generated_code = F90InterfaceGenerator.apply(plugin)
+    generated_code = F90InterfaceGenerator.apply(F90Interface(cffi_plugin=plugin))
     return format_fortran_code(generated_code)
diff --git a/tools/src/icon4pytools/py2fgen/parsing.py b/tools/src/icon4pytools/py2fgen/parsing.py
index eb22becc8d..1cc8ede610 100644
--- a/tools/src/icon4pytools/py2fgen/parsing.py
+++ b/tools/src/icon4pytools/py2fgen/parsing.py
@@ -27,7 +27,7 @@
 from icon4pytools.py2fgen.utils import parse_type_spec
 
 
-class ImportExtractor(ast.NodeVisitor):
+class ImportStmtVisitor(ast.NodeVisitor):
     """AST Visitor to extract import statements."""
 
     def __init__(self):
@@ -48,8 +48,8 @@ def visit_ImportFrom(self, node):
             self.import_statements.append(import_statement)
 
 
-class FuncDefVisitor(ast.NodeVisitor):
-    """AST Visitor to extract function definitions and type hints."""
+class TypeHintVisitor(ast.NodeVisitor):
+    """AST Visitor to extract function parameter type hints."""
 
     def __init__(self):
         self.type_hints: dict[str, str] = {}
@@ -59,20 +59,36 @@ def visit_FunctionDef(self, node):
             if arg.annotation:
                 annotation = ast.unparse(arg.annotation)
                 self.type_hints[arg.arg] = annotation
+            else:
+                raise TypeError(
+                    f"Missing type hint for parameter '{arg.arg}' in function '{node.name}'"
+                )
 
 
-def parse(module_name: str, function_name: str) -> CffiPlugin:
+def parse(module_name: str, functions: list[str], plugin_name: str) -> CffiPlugin:
     module = importlib.import_module(module_name)
-    parsed_function = _parse_function(module, function_name)
     parsed_imports = _extract_import_statements(module)
+
+    parsed_functions: list[Func] = []
+    for f in functions:
+        parsed_functions.append(_parse_function(module, f))
+
     return CffiPlugin(
         module_name=module_name,
-        plugin_name=f"{function_name}_plugin",
-        function=parsed_function,
+        plugin_name=plugin_name,
+        functions=parsed_functions,
         imports=parsed_imports,
     )
 
 
+def _extract_import_statements(module: ModuleType) -> list[str]:
+    src = inspect.getsource(module)
+    tree = ast.parse(src)
+    visitor = ImportStmtVisitor()
+    visitor.visit(tree)
+    return visitor.import_statements
+
+
 def _parse_function(module: ModuleType, function_name: str) -> Func:
     func = unwrap(getattr(module, function_name))
     is_gt4py_program = isinstance(func, Program)
@@ -87,6 +103,32 @@ def _parse_function(module: ModuleType, function_name: str) -> Func:
     return Func(name=function_name, args=params, is_gt4py_program=is_gt4py_program)
 
 
+def _extract_type_hint_strings(
+    module: ModuleType, func: Callable, is_gt4py_program: bool, function_name: str
+):
+    src = extract_function_signature(
+        inspect.getsource(module) if is_gt4py_program else inspect.getsource(func), function_name
+    )
+    tree = ast.parse(src)
+    visitor = TypeHintVisitor()
+    visitor.visit(tree)
+    return visitor.type_hints
+
+
+def extract_function_signature(code: str, function_name: str) -> str:
+    # This pattern attempts to match function definitions
+    pattern = rf"\bdef\s+{re.escape(function_name)}\s*\(([\s\S]*?)\)\s*:"
+
+    match = re.search(pattern, code)
+
+    if match:
+        # Constructing the full signature with empty return for ease of parsing by AST visitor
+        signature = match.group()
+        return signature.strip() + "\n  return None"
+    else:
+        raise Exception(f"Could not parse function signature from the following code:\n {code}")
+
+
 def _get_gt4py_func_params(func: Program, type_hints: dict[str, str]) -> List[FuncParameter]:
     return [
         FuncParameter(
@@ -113,37 +155,3 @@ def _get_simple_func_params(func: Callable, type_hints: dict[str, str]) -> List[
         )
         for s, param in sig_params.items()
     ]
-
-
-def _extract_type_hint_strings(
-    module: ModuleType, func: Callable, is_gt4py_program: bool, function_name: str
-):
-    src = extract_function_signature(
-        inspect.getsource(module) if is_gt4py_program else inspect.getsource(func), function_name
-    )
-    tree = ast.parse(src)
-    visitor = FuncDefVisitor()
-    visitor.visit(tree)
-    return visitor.type_hints
-
-
-def extract_function_signature(code: str, function_name: str) -> str:
-    # This pattern attempts to match function definitions
-    pattern = rf"\bdef\s+{re.escape(function_name)}\s*\(([\s\S]*?)\)\s*:"
-
-    match = re.search(pattern, code)
-
-    if match:
-        # Constructing the full signature with empty return for ease of parsing by AST visitor
-        signature = match.group()
-        return signature.strip() + "\n  return None"
-    else:
-        raise Exception(f"Could not parse function signature from the following code:\n {code}")
-
-
-def _extract_import_statements(module: ModuleType) -> list[str]:
-    src = inspect.getsource(module)
-    tree = ast.parse(src)
-    visitor = ImportExtractor()
-    visitor.visit(tree)
-    return visitor.import_statements
diff --git a/tools/src/icon4pytools/py2fgen/plugin.py b/tools/src/icon4pytools/py2fgen/plugin.py
new file mode 100644
index 0000000000..8a85c90594
--- /dev/null
+++ b/tools/src/icon4pytools/py2fgen/plugin.py
@@ -0,0 +1,182 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# This file is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or any later
+# version. See the LICENSE.txt file at the top-level directory of this
+# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+import logging
+import typing
+from pathlib import Path
+
+import cffi
+import numpy as np
+from cffi import FFI
+from numpy.typing import NDArray
+
+from icon4pytools.common.logger import setup_logger
+
+
+if typing.TYPE_CHECKING:
+    import cupy as cp  # type: ignore
+
+ffi = FFI()  # needed for unpack and unpack_gpu functions
+
+logger = setup_logger(__name__)
+
+
+def unpack(ptr, *sizes: int) -> NDArray:
+    """
+    Converts a C pointer into a NumPy array to directly manipulate memory allocated in Fortran.
+    This function is needed for operations requiring in-place modification of CPU data, enabling
+    changes made in Python to reflect immediately in the original Fortran memory space.
+
+    Args:
+        ptr (CData): A CFFI pointer to the beginning of the data array in CPU memory. This pointer
+                     should reference a contiguous block of memory whose total size matches the product
+                     of the specified dimensions.
+        *sizes (int): Variable length argument list specifying the dimensions of the array.
+                      These sizes determine the shape of the resulting NumPy array.
+
+    Returns:
+        np.ndarray: A NumPy array that provides a direct view of the data pointed to by `ptr`.
+                    This array shares the underlying data with the original Fortran code, allowing
+                    modifications made through the array to affect the original data.
+    """
+    length = np.prod(sizes)
+    c_type = ffi.getctype(ffi.typeof(ptr).item)
+
+    # Map C data types to NumPy dtypes
+    dtype_map: dict[str, np.dtype] = {
+        "int": np.dtype(np.int32),
+        "double": np.dtype(np.float64),
+    }
+    dtype = dtype_map.get(c_type, np.dtype(c_type))
+
+    # Create a NumPy array from the buffer, specifying the Fortran order
+    arr = np.frombuffer(ffi.buffer(ptr, length * ffi.sizeof(c_type)), dtype=dtype).reshape(  # type: ignore
+        sizes, order="F"
+    )
+    return arr
+
+
+def unpack_gpu(ptr, *sizes: int):
+    """
+    Converts a C pointer into a CuPy array to directly manipulate memory allocated in Fortran.
+    This function is needed for operations that require in-place modification of GPU data,
+    enabling changes made in Python to reflect immediately in the original Fortran memory space.
+
+    Args:
+        ptr (cffi.CData): A CFFI pointer to GPU memory allocated by OpenACC, representing
+                          the starting address of the data. This pointer must correspond to
+                          a contiguous block of memory whose total size matches the product
+                          of the specified dimensions.
+        *sizes (int): Variable length argument list specifying the dimensions of the array.
+                      These sizes determine the shape of the resulting CuPy array.
+
+    Returns:
+        cp.ndarray: A CuPy array that provides a direct view of the data pointed to by `ptr`.
+                    This array shares the underlying data with the original Fortran code, allowing
+                    modifications made through the array to affect the original data.
+    """
+
+    if not sizes:
+        raise ValueError("Sizes must be provided to determine the array shape.")
+
+    length = np.prod(sizes)
+    c_type = ffi.getctype(ffi.typeof(ptr).item)
+
+    dtype_map = {
+        "int": cp.int32,
+        "double": cp.float64,
+    }
+    dtype = dtype_map.get(c_type, None)
+    if dtype is None:
+        raise ValueError(f"Unsupported C data type: {c_type}")
+
+    itemsize = ffi.sizeof(c_type)
+    total_size = length * itemsize
+
+    # cupy array from OpenACC device pointer
+    current_device = cp.cuda.Device()
+    ptr_val = int(ffi.cast("uintptr_t", ptr))
+    mem = cp.cuda.UnownedMemory(ptr_val, total_size, owner=ptr, device_id=current_device.id)
+    memptr = cp.cuda.MemoryPointer(mem, 0)
+    arr = cp.ndarray(shape=sizes, dtype=dtype, memptr=memptr, order="F")
+
+    return arr
+
+
+def int_array_to_bool_array(int_array: NDArray) -> NDArray:
+    """
+    Converts a NumPy array of integers to a boolean array.
+    In the input array, 0 represents False, and any non-zero value (1 or -1) represents True.
+
+    Args:
+        int_array: A NumPy array of integers.
+
+    Returns:
+        A NumPy array of booleans.
+    """
+    bool_array = int_array != 0
+    return bool_array
+
+
+def generate_and_compile_cffi_plugin(
+    plugin_name: str, c_header: str, python_wrapper: str, build_path: Path
+) -> None:
+    """
+    Create and compile a CFFI plugin.
+
+    This function generates a C shared library and Fortran interface for Python functions
+    to be exposed in the {plugin_name} module. It creates a linkable C library named
+    'lib{plugin_name}.so' in the specified build directory.
+
+    Args:
+        plugin_name: Name of the plugin.
+        c_header: C header signatures for the Python functions.
+        python_wrapper: Python code wrapping the original function to be exposed.
+        build_path: Path to the build directory.
+    """
+    try:
+        header_file_path = write_c_header(build_path, plugin_name, c_header)
+        compile_cffi_plugin(
+            builder=configure_cffi_builder(c_header, plugin_name, header_file_path),
+            python_wrapper=python_wrapper,
+            build_path=str(build_path),
+            plugin_name=plugin_name,
+        )
+    except Exception as e:
+        logging.error(f"Error generating and compiling CFFI plugin: {e}")
+        raise
+
+
+def write_c_header(build_path: Path, plugin_name: str, c_header: str) -> Path:
+    """Write the C header file to the specified path."""
+    c_header_file = plugin_name + ".h"
+    header_file_path = build_path / c_header_file
+    with open(header_file_path, "w") as f:
+        f.write(c_header)
+    return header_file_path
+
+
+def configure_cffi_builder(c_header: str, plugin_name: str, header_file_path: Path) -> cffi.FFI:
+    """Configure and returns a CFFI FFI builder instance."""
+    builder = cffi.FFI()
+    builder.embedding_api(c_header)
+    builder.set_source(plugin_name, f'#include "{header_file_path.name}"')
+    return builder
+
+
+def compile_cffi_plugin(
+    builder: cffi.FFI, python_wrapper: str, build_path: str, plugin_name: str
+) -> None:
+    """Compile the CFFI plugin with the given configuration."""
+    logger.info("Compiling CFFI dynamic library...")
+    builder.embedding_init_code(python_wrapper)
+    builder.compile(tmpdir=build_path, target=f"lib{plugin_name}.*", verbose=True)
diff --git a/tools/src/icon4pytools/py2fgen/template.py b/tools/src/icon4pytools/py2fgen/template.py
index 81652efcb8..b2be3fbbda 100644
--- a/tools/src/icon4pytools/py2fgen/template.py
+++ b/tools/src/icon4pytools/py2fgen/template.py
@@ -11,19 +11,20 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 import inspect
-from typing import Any, Optional, Sequence
+from typing import Any, Sequence
 
 from gt4py.eve import Node, datamodels
 from gt4py.eve.codegen import JinjaTemplate as as_jinja, TemplatedGenerator
 from gt4py.next import Dimension
 from gt4py.next.type_system.type_specifications import ScalarKind
+from icon4py.model.common.config import GT4PyBackend
 
 from icon4pytools.icon4pygen.bindings.codegen.type_conversion import (
     BUILTIN_TO_CPP_TYPE,
     BUILTIN_TO_ISO_C_TYPE,
 )
-from icon4pytools.py2fgen.cffi import unpack
-from icon4pytools.py2fgen.utils import build_array_size_args, flatten_and_get_unique_elts
+from icon4pytools.py2fgen.plugin import int_array_to_bool_array, unpack, unpack_gpu
+from icon4pytools.py2fgen.utils import flatten_and_get_unique_elts
 
 
 CFFI_DECORATOR = "@ffi.def_extern()"
@@ -47,7 +48,12 @@ class FuncParameter(Node):
     def __post_init__(self):
         self.size_args = dims_to_size_strings(self.dimensions)
         self.is_array = True if len(self.dimensions) >= 1 else False
-        self.gtdims = [dimension.value + "Dim" for dimension in self.dimensions]
+        # We need some fields to have nlevp1 levels on the fortran wrapper side, which we make
+        # happen by using KHalfDim as a type hint. However, this is not yet supported on the icon4py
+        # side. So before generating the python wrapper code, we replace occurences of KHalfDim with KDim
+        self.gtdims = [
+            dimension.value.replace("KHalf", "K") + "Dim" for dimension in self.dimensions
+        ]
 
 
 class Func(Node):
@@ -66,14 +72,36 @@ class CffiPlugin(Node):
     module_name: str
     plugin_name: str
     imports: list[str]
-    function: Func
+    functions: list[Func]
 
 
 class PythonWrapper(CffiPlugin):
-    gt4py_backend: Optional[str]
+    backend: str
     debug_mode: bool
     cffi_decorator: str = CFFI_DECORATOR
     cffi_unpack: str = inspect.getsource(unpack)
+    cffi_unpack_gpu: str = inspect.getsource(unpack_gpu)
+    int_to_bool: str = inspect.getsource(int_array_to_bool_array)
+    gt4py_backend: str = datamodels.field(init=False)
+    is_gt4py_program_present: bool = datamodels.field(init=False)
+
+    def __post_init__(self, *args: Any, **kwargs: Any) -> None:
+        self.gt4py_backend = GT4PyBackend[self.backend].value
+        self.is_gt4py_program_present = any(func.is_gt4py_program for func in self.functions)
+
+
+def build_array_size_args() -> dict[str, str]:
+    array_size_args = {}
+    from icon4py.model.common import dimension
+
+    for var_name, var in vars(dimension).items():
+        if isinstance(var, Dimension):
+            dim_name = var_name.replace(
+                "Dim", ""
+            )  # Assumes we keep suffixing each Dimension with Dim in icon4py.common.dimension module
+            size_name = f"n_{dim_name}"
+            array_size_args[dim_name] = size_name
+    return array_size_args
 
 
 def to_c_type(scalar_type: ScalarKind) -> str:
@@ -81,7 +109,7 @@ def to_c_type(scalar_type: ScalarKind) -> str:
     return BUILTIN_TO_CPP_TYPE[scalar_type]
 
 
-def to_f_type(scalar_type: ScalarKind) -> str:
+def to_iso_c_type(scalar_type: ScalarKind) -> str:
     """Convert a scalar type to its corresponding ISO C type."""
     return BUILTIN_TO_ISO_C_TYPE[scalar_type]
 
@@ -104,7 +132,7 @@ def render_c_pointer(param: FuncParameter) -> str:
     return "*" if len(param.dimensions) > 0 else ""
 
 
-def render_fortran_array_dimensions(param: FuncParameter) -> str:
+def render_fortran_array_dimensions(param: FuncParameter, assumed_size_array: bool) -> str:
     """
     Render Fortran array dimensions for array types.
 
@@ -114,9 +142,13 @@ def render_fortran_array_dimensions(param: FuncParameter) -> str:
     Returns:
         A string representing Fortran array dimensions.
     """
-    if len(param.dimensions) > 0:
+    if len(param.dimensions) > 0 and assumed_size_array:
+        return "dimension(*),"
+
+    if len(param.dimensions) > 0 and not assumed_size_array:
         dims = ",".join(":" for _ in param.dimensions)
         return f"dimension({dims}),"
+
     return ""
 
 
@@ -162,77 +194,144 @@ def render_fortran_array_sizes(param: FuncParameter) -> str:
 class PythonWrapperGenerator(TemplatedGenerator):
     PythonWrapper = as_jinja(
         """\
-# necessary imports for generated code to work
+# imports for generated wrapper code
+import logging
 from {{ plugin_name }} import ffi
 import numpy as np
-from gt4py.next.ffront.fbuiltins import int32
+{% if _this_node.backend == 'GPU' %}import cupy as cp {% endif %}
+from numpy.typing import NDArray
 from gt4py.next.iterator.embedded import np_as_located_field
-from gt4py.next.program_processors.runners.gtfn import run_gtfn, run_gtfn_gpu
+
+{% if _this_node.is_gt4py_program_present %}
+# necessary imports when embedding a gt4py program directly
+from gt4py.next.program_processors.runners.gtfn import run_gtfn_cached, run_gtfn_gpu_cached
 from gt4py.next.program_processors.runners.roundtrip import backend as run_roundtrip
 from icon4py.model.common.grid.simple import SimpleGrid
 
-# all other imports from the module from which the function is being wrapped
+# We need a grid to pass offset providers to the embedded gt4py program (granules load their own grid at runtime)
+grid = SimpleGrid()
+{% endif %}
+
+# logger setup
+log_format = '%(asctime)s.%(msecs)03d - %(levelname)s - %(message)s'
+logging.basicConfig(level=logging.{%- if _this_node.debug_mode -%}DEBUG{%- else -%}ERROR{%- endif -%},
+                    format=log_format,
+                    datefmt='%Y-%m-%d %H:%M:%S')
+{% if _this_node.backend == 'GPU' %}logging.info(cp.show_config()) {% endif %}
+
+# embedded module imports
 {% for stmt in imports -%}
 {{ stmt }}
 {% endfor %}
 
-# We need a grid to pass offset providers
-grid = SimpleGrid()
-
-from {{ module_name }} import {{ _this_node.function.name }}
+# embedded function imports
+{% for func in _this_node.functions -%}
+from {{ module_name }} import {{ func.name }}
+{% endfor %}
 
+{% if _this_node.backend == 'GPU' %}
+{{ cffi_unpack_gpu }}
+{% else %}
 {{ cffi_unpack }}
+{% endif %}
+
+{{ int_to_bool }}
+
+{% for func in _this_node.functions %}
 
 {{ cffi_decorator }}
-def {{ _this_node.function.name }}_wrapper(
-{%- for arg in _this_node.function.args -%}
-{{ arg.name }}: {{ arg.py_type_hint }}{% if not loop.last or _this_node.function.global_size_args %}, {% endif %}
+def {{ func.name }}_wrapper(
+{%- for arg in func.args -%}
+{{ arg.name }}: {{ arg.py_type_hint | replace("KHalfDim","KDim") }}{% if not loop.last or func.global_size_args %}, {% endif %}
 {%- endfor %}
-{%- for arg in _this_node.function.global_size_args -%}
+{%- for arg in func.global_size_args -%}
 {{ arg }}: int32{{ ", " if not loop.last else "" }}
 {%- endfor -%}
 ):
+    try:
+        {%- if _this_node.debug_mode %}
+        logging.info("Python Execution Context Start")
+        {% endif %}
+
+        # Unpack pointers into Ndarrays
+        {% for arg in func.args %}
+        {% if arg.is_array %}
+        {%- if _this_node.debug_mode %}
+        msg = '{{ arg.name }} before unpacking: %s' % str({{ arg.name}})
+        logging.debug(msg)
+        {% endif %}
+        {{ arg.name }} = unpack{%- if _this_node.backend == 'GPU' -%}_gpu{%- endif -%}({{ arg.name }}, {{ ", ".join(arg.size_args) }})
+
+        {%- if arg.d_type.name == "BOOL" %}
+        {{ arg.name }} = int_array_to_bool_array({{ arg.name }})
+        {%- endif %}
+
+        {%- if _this_node.debug_mode %}
+        msg = '{{ arg.name }} after unpacking: %s' % str({{ arg.name}})
+        logging.debug(msg)
+        msg = 'shape of {{ arg.name }} after unpacking = %s' % str({{ arg.name}}.shape)
+        logging.debug(msg)
+        {% endif %}
+        {% endif %}
+        {% endfor %}
+
+        # Allocate GT4Py Fields
+        {% for arg in func.args %}
+        {% if arg.is_array %}
+        {{ arg.name }} = np_as_located_field({{ ", ".join(arg.gtdims) }})({{ arg.name }})
+        {%- if _this_node.debug_mode %}
+        msg = 'shape of {{ arg.name }} after allocating as field = %s' % str({{ arg.name}}.shape)
+        logging.debug(msg)
+        msg = '{{ arg.name }} after allocating as field: %s' % str({{ arg.name }}.ndarray)
+        logging.debug(msg)
+        {% endif %}
+        {% endif %}
+        {% endfor %}
+
+        {{ func.name }}
+        {%- if func.is_gt4py_program -%}.with_backend({{ _this_node.gt4py_backend }}){%- endif -%}(
+        {%- for arg in func.args -%}
+        {{ arg.name }}{{ ", " if not loop.last or func.is_gt4py_program else "" }}
+        {%- endfor -%}
+        {%- if func.is_gt4py_program -%}
+        offset_provider=grid.offset_providers
+        {%- endif -%}
+        )
 
-    # Unpack pointers into Ndarrays
-    {% for arg in _this_node.function.args %}
-    {% if arg.is_array %}
-    {{ arg.name }} = unpack({{ arg.name }}, {{ ", ".join(arg.size_args) }})
-    {%- if _this_node.debug_mode %}
-    msg = 'shape of {{ arg.name }} = %s' % str({{ arg.name}}.shape)
-    print(msg)
-    {% endif %}
-    {% endif %}
-    {% endfor %}
+        {% if _this_node.debug_mode %}
+        # debug info
+        {% for arg in func.args %}
+        {% if arg.is_array %}
+        msg = 'shape of {{ arg.name }} after computation = %s' % str({{ arg.name}}.shape)
+        logging.debug(msg)
+        msg = '{{ arg.name }} after computation: %s' % str({{ arg.name }}.ndarray)
+        logging.debug(msg)
+        {% endif %}
+        {% endfor %}
+        {% endif %}
+
+        {%- if _this_node.debug_mode %}
+        logging.info("Python Execution Context End")
+        {% endif %}
+
+    except Exception as e:
+        logging.exception(f"A Python error occurred: {e}")
+        return 1
+
+    return 0
+
+        {% endfor %}
 
-    # Allocate GT4Py Fields
-    {% for arg in _this_node.function.args %}
-    {% if arg.is_array %}
-    {{ arg.name }} = np_as_located_field({{ ", ".join(arg.gtdims) }})({{ arg.name }})
-    {%- if _this_node.debug_mode %}
-    print({{ arg.name }})
-    print({{ arg.name }}.shape)
-    {% endif %}
-    {% endif %}
-    {% endfor %}
 
-    {{ _this_node.function.name }}
-    {%- if _this_node.function.is_gt4py_program -%}.with_backend({{ _this_node.gt4py_backend }}){%- endif -%}(
-    {%- for arg in _this_node.function.args -%}
-    {{ arg.name }}{{ ", " if not loop.last or _this_node.function.is_gt4py_program else "" }}
-    {%- endfor -%}
-    {%- if _this_node.function.is_gt4py_program -%}
-    offset_provider=grid.offset_providers
-    {%- endif -%}
-)
 """
     )
 
 
 class CHeaderGenerator(TemplatedGenerator):
-    CffiPlugin = as_jinja("""extern void {{_this_node.function.name}}_wrapper({{function}});""")
+    CffiPlugin = as_jinja("""{{'\n'.join(functions)}}""")
 
     Func = as_jinja(
-        "{%- for arg in args -%}{{ arg }}{% if not loop.last or global_size_args|length > 0 %}, {% endif %}{% endfor -%}{%- for sarg in global_size_args -%} int {{ sarg }}{% if not loop.last %}, {% endif %}{% endfor -%}"
+        "extern int {{ name }}_wrapper({%- for arg in args -%}{{ arg }}{% if not loop.last or global_size_args|length > 0 %}, {% endif %}{% endfor -%}{%- for sarg in global_size_args -%} int {{ sarg }}{% if not loop.last %}, {% endif %}{% endfor -%});"
     )
 
     def visit_FuncParameter(self, param: FuncParameter):
@@ -243,38 +342,134 @@ def visit_FuncParameter(self, param: FuncParameter):
     FuncParameter = as_jinja("""{{rendered_type}}{{pointer}} {{name}}""")
 
 
+class F90FunctionDeclaration(Func):
+    def __post_init__(self, *args: Any, **kwargs: Any) -> None:
+        super().__post_init__()  # call Func __post_init__
+
+
+class DimensionPosition(Node):
+    variable: str
+    size_arg: str
+    index: int
+
+
+class F90FunctionDefinition(Func):
+    dimension_size_declarations: Sequence[DimensionPosition] = datamodels.field(init=False)
+
+    def __post_init__(self, *args: Any, **kwargs: Any) -> None:
+        super().__post_init__()  # call Func __post_init__
+
+        dim_positions = []
+        for arg in self.args:
+            for index, size_arg in enumerate(arg.size_args):
+                dim_positions.append(
+                    DimensionPosition(variable=str(arg.name), size_arg=size_arg, index=index + 1)
+                )  # Use Fortran indexing
+
+        self.dimension_size_declarations = dim_positions
+
+
+class F90Interface(Node):
+    cffi_plugin: CffiPlugin
+    function_declaration: list[F90FunctionDeclaration] = datamodels.field(init=False)
+    function_definition: list[F90FunctionDefinition] = datamodels.field(init=False)
+
+    def __post_init__(self, *args: Any, **kwargs: Any) -> None:
+        functions = self.cffi_plugin.functions
+        self.function_declaration = [
+            F90FunctionDeclaration(name=f.name, args=f.args, is_gt4py_program=f.is_gt4py_program)
+            for f in functions
+        ]
+        self.function_definition = [
+            F90FunctionDefinition(name=f.name, args=f.args, is_gt4py_program=f.is_gt4py_program)
+            for f in functions
+        ]
+
+
 class F90InterfaceGenerator(TemplatedGenerator):
-    CffiPlugin = as_jinja(
+    F90Interface = as_jinja(
         """\
-    module {{ plugin_name }}
-    use, intrinsic:: iso_c_binding
+module {{ _this_node.cffi_plugin.plugin_name }}
+    use, intrinsic :: iso_c_binding
     implicit none
 
-    public
-    interface
-    {{ function }}
-    end interface
-    end module
-    """
+    {% for func in _this_node.cffi_plugin.functions %}
+    public :: {{ func.name }}
+    {% endfor %}
+
+interface
+    {{ '\n'.join(function_declaration) }}
+end interface
+
+contains
+    {{ '\n'.join(function_definition) }}
+end module
+"""
     )
 
-    def visit_Func(self, func: Func, **kwargs):
+    def visit_F90FunctionDeclaration(self, func: F90FunctionDeclaration, **kwargs):
         arg_names = ", &\n ".join(map(lambda x: x.name, func.args))
         if func.global_size_args:
             arg_names += ",&\n" + ", &\n".join(func.global_size_args)
+        return self.generic_visit(func, assumed_size_array=True, param_names=arg_names)
+
+    F90FunctionDeclaration = as_jinja(
+        """
+function {{name}}_wrapper({{param_names}}) bind(c, name="{{name}}_wrapper") result(rc)
+   import :: c_int, c_double, c_bool, c_ptr
+   {% for size_arg in global_size_args %}
+   integer(c_int), value :: {{ size_arg }}
+   {% endfor %}
+   integer(c_int) :: rc  ! Stores the return code
+   {% for arg in args %}
+   {{ arg }}
+   {% endfor %}
+end function {{name}}_wrapper
+    """
+    )
 
-        return self.generic_visit(func, param_names=arg_names)
+    def visit_F90FunctionDefinition(self, func: F90FunctionDefinition, **kwargs):
+        if len(func.args) < 1:
+            arg_names, param_names_with_size_args = "", ""
+        else:
+            arg_names = ", &\n ".join(map(lambda x: x.name, func.args))
+            param_names_with_size_args = arg_names + ",&\n" + ", &\n".join(func.global_size_args)
 
-    Func = as_jinja(
-        """subroutine {{name}}_wrapper({{param_names}}) bind(c, name='{{name}}_wrapper')
-       use, intrinsic :: iso_c_binding
-       {% for size_arg in global_size_args -%}
-       integer(c_int), value, target :: {{ size_arg }}
-       {% endfor %}
-       {% for arg in args: %}\
-       {{arg}}\
-       {% endfor %}\
-    end subroutine {{name}}_wrapper
+        return self.generic_visit(
+            func,
+            assumed_size_array=False,
+            param_names=arg_names,
+            param_names_with_size_args=param_names_with_size_args,
+            arrays=[arg for arg in func.args if arg.is_array],
+        )
+
+    # todo(samkellerhals): Consider using unique SIZE args
+    F90FunctionDefinition = as_jinja(
+        """
+subroutine {{name}}({{param_names}}, &\nrc)
+   use, intrinsic :: iso_c_binding
+   {% for size_arg in global_size_args %}
+   integer(c_int) :: {{ size_arg }}
+   {% endfor %}
+   {% for arg in args %}
+   {{ arg }}
+   {% endfor %}
+   integer(c_int) :: rc  ! Stores the return code
+
+   !$ACC host_data use_device( &
+   {%- for arr in arrays %}
+       !$ACC {{ arr.name }}{% if not loop.last %}, &{% else %} &{% endif %}
+   {%- endfor %}
+   !$ACC )
+
+   {% for d in _this_node.dimension_size_declarations %}
+   {{ d.size_arg }} = SIZE({{ d.variable }}, {{ d.index }})
+   {% endfor %}
+
+   rc = {{ name }}_wrapper({{ param_names_with_size_args }})
+
+   !$acc end host_data
+end subroutine {{name}}
     """
     )
 
@@ -282,12 +477,9 @@ def visit_FuncParameter(self, param: FuncParameter, **kwargs):
         return self.generic_visit(
             param,
             value=as_f90_value(param),
-            rendered_type=to_f_type(param.d_type),
-            dim=render_fortran_array_dimensions(param),
+            iso_c_type=to_iso_c_type(param.d_type),
+            dim=render_fortran_array_dimensions(param, kwargs["assumed_size_array"]),
             explicit_size=render_fortran_array_sizes(param),
         )
 
-    FuncParameter = as_jinja(
-        """{{rendered_type}}, {{dim}} {{value}} target :: {{name}}{{ explicit_size }}
-    """
-    )
+    FuncParameter = as_jinja("""{{iso_c_type}}, {{dim}} {{value}} target :: {{name}}""")
diff --git a/tools/src/icon4pytools/py2fgen/utils.py b/tools/src/icon4pytools/py2fgen/utils.py
index 28923c776f..0313c04aa6 100644
--- a/tools/src/icon4pytools/py2fgen/utils.py
+++ b/tools/src/icon4pytools/py2fgen/utils.py
@@ -10,33 +10,14 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
-
-from enum import Enum
+import importlib
+import os
+from pathlib import Path
 
 from gt4py.next.common import Dimension
 from gt4py.next.type_system.type_specifications import FieldType, ScalarKind, ScalarType, TypeSpec
 
 
-def build_array_size_args() -> dict[str, str]:
-    array_size_args = {}
-    from icon4py.model.common import dimension
-
-    for var_name, var in vars(dimension).items():
-        if isinstance(var, Dimension):
-            dim_name = var_name.replace(
-                "Dim", ""
-            )  # Assumes we keep suffixing each Dimension with Dim
-            size_name = f"n_{dim_name}"
-            array_size_args[dim_name] = size_name
-    return array_size_args
-
-
-class Backend(Enum):
-    CPU = "run_gtfn"
-    GPU = "run_gtfn_gpu"
-    ROUNDTRIP = "run_roundtrip"
-
-
 def parse_type_spec(type_spec: TypeSpec) -> tuple[list[Dimension], ScalarKind]:
     if isinstance(type_spec, ScalarType):
         return [], type_spec.kind
@@ -48,3 +29,34 @@ def parse_type_spec(type_spec: TypeSpec) -> tuple[list[Dimension], ScalarKind]:
 
 def flatten_and_get_unique_elts(list_of_lists: list[list[str]]) -> list[str]:
     return sorted(set(item for sublist in list_of_lists for item in sublist))
+
+
+def get_local_test_grid():
+    test_folder = "testdata"
+    module_spec = importlib.util.find_spec("icon4pytools")
+
+    if module_spec and module_spec.origin:
+        # following namespace package conventions the root is three levels down
+        repo_root = Path(module_spec.origin).parents[3]
+        return os.path.join(repo_root, test_folder)
+    else:
+        raise FileNotFoundError(
+            "The `icon4pytools` package could not be found. Ensure the package is installed "
+            "and accessible. Alternatively, set the 'ICON_GRID_LOC' environment variable "
+            "explicitly to specify the location."
+        )
+
+
+def get_icon_grid_loc():
+    env_path = os.environ.get("ICON_GRID_LOC")
+    if env_path is not None:
+        return env_path
+    else:
+        return get_local_test_grid()
+
+
+def get_grid_filename():
+    env_path = os.environ.get("ICON_GRID_NAME")
+    if env_path is not None:
+        return env_path
+    return "grid.nc"
diff --git a/tools/src/icon4pytools/py2fgen/wrappers/diffusion.py b/tools/src/icon4pytools/py2fgen/wrappers/diffusion.py
index 8a8cb6ed71..855236c890 100644
--- a/tools/src/icon4pytools/py2fgen/wrappers/diffusion.py
+++ b/tools/src/icon4pytools/py2fgen/wrappers/diffusion.py
@@ -10,10 +10,7 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
-# flake8: noqa
-
-# We use gt4py type annotations and thus need to ignore this in MyPy
-# mypy: disable-error-code="valid-type"
+# type: ignore
 
 """
 Wrapper module for diffusion granule.
@@ -24,202 +21,189 @@
 - passing of scalar types or fields of simple types
 """
 
-import numpy as np
 from gt4py.next.common import Field
-from gt4py.next.ffront.fbuiltins import int32
+from gt4py.next.ffront.fbuiltins import float64, int32
 from icon4py.model.atmosphere.diffusion.diffusion import (
     Diffusion,
     DiffusionConfig,
+    DiffusionParams,
+    TurbulenceShearForcingType,
+)
+from icon4py.model.atmosphere.diffusion.diffusion_states import (
     DiffusionDiagnosticState,
     DiffusionInterpolationState,
     DiffusionMetricState,
-    DiffusionParams,
-    DiffusionType,
 )
 from icon4py.model.common.dimension import (
     C2E2CDim,
     C2E2CODim,
     C2EDim,
     CECDim,
+    CEDim,
     CellDim,
     E2C2VDim,
     E2CDim,
-    E2VDim,
     ECVDim,
     EdgeDim,
     KDim,
+    KHalfDim,
     V2EDim,
     VertexDim,
 )
-from icon4py.model.common.grid.base import GridConfig
-from icon4py.model.common.grid.horizontal import CellParams, EdgeParams, HorizontalGridSize
-from icon4py.model.common.grid.icon import IconGrid
-from icon4py.model.common.grid.vertical import VerticalGridSize, VerticalModelParams
+from icon4py.model.common.grid.horizontal import CellParams, EdgeParams
+from icon4py.model.common.grid.vertical import VerticalModelParams
+from icon4py.model.common.settings import device
 from icon4py.model.common.states.prognostic_state import PrognosticState
+from icon4py.model.common.test_utils.grid_utils import _load_from_gridfile
+from icon4py.model.common.test_utils.helpers import as_1D_sparse_field, flatten_first_two_dims
 
+from icon4pytools.common.logger import setup_logger
+from icon4pytools.py2fgen.utils import get_grid_filename, get_icon_grid_loc
 
-# TODO (magdalena) Revise interface architecture with Fortran granules:
-# The module variable to match the Fortran interface: where only fields are passed.
-# We should rather instantiate the object init and return it.
-diffusion: Diffusion()
 
-nproma = 50000
-field_sizes = {EdgeDim: nproma, CellDim: nproma, VertexDim: nproma}
+logger = setup_logger(__name__)
+
+# global diffusion object
+diffusion_granule: Diffusion = Diffusion()
 
 
 def diffusion_init(
-    nproma: int,
-    nlev: int,
-    n_shift_total: int,
-    nrdmax: float,
-    n_dyn_substeps: int,
-    nudge_max_coeff: float,
-    denom_diffu_v: float,
-    hdiff_smag_fac: float,
-    hdiff_order: int,
-    hdiff_efdt_ratio: float,
-    itype_sher: int,
-    itype_vn_diffu: int,
-    itype_t_diffu: int,
-    lhdiff_w: bool,
-    lhdiff_temp: bool,
-    l_limited_area: bool,
-    l_zdiffu_t: bool,
-    lsmag_3d: bool,
-    vct_a: Field[[KDim], float],
-    e_bln_c_s: Field[[CellDim, C2EDim], float],
-    geofac_div: Field[[CellDim, C2EDim], float],
-    geofac_n2s: Field[[CellDim, C2E2CODim], float],
-    geofac_grg_x: Field[[CellDim, C2E2CODim], float],
-    geofac_grg_y: Field[[CellDim, C2E2CODim], float],
-    nudgecoeff_e: Field[[EdgeDim], float],
-    zd_intcoef: Field[
-        [CellDim, C2E2CDim, KDim], float
-    ],  # special DSL field:   zd_intcoef_dsl in mo_vertical_grid.f90
-    zd_diffcoef: Field[
-        [CellDim, KDim], float
-    ],  # special DSL field mask instead of list: zd_diffcoef_dsl in mo_vertical_grid.f90
-    wgtfac_c: Field[[CellDim, KDim], float],
-    theta_ref_mc: Field[[CellDim, KDim], float],
-    edges_vertex_idx: Field[[EdgeDim, E2VDim], int32],
-    edges_cell_idx: Field[[EdgeDim, E2CDim], int32],
-    edges_tangent_orientation: Field[[EdgeDim], float],
-    edges_primal_normal_vert_1: Field[[ECVDim], float],  # shallow derived type in Fortran
-    edges_primal_normal_vert_2: Field[[ECVDim], float],  # shallow derived type in Fortran
-    edges_dual_normal_vert_1: Field[[ECVDim], float],  # shallow derived type in Fortran
-    edges_dual_normal_vert_2: Field[[ECVDim], float],  # shallow derived type in Fortran
-    edges_inv_vert_vert_lengths: Field[[EdgeDim], float],
-    edges_inv_primal_edge_length: Field[[EdgeDim], float],
-    edges_inv_dual_edge_length: Field[[EdgeDim], float],
-    edges_area_edge: Field[[EdgeDim], float],
-    cells_neighbor_idx: Field[[CellDim, C2E2CDim], int32],
-    cells_edge_idx: Field[[CellDim, C2EDim], int32],
-    cells_area: Field[[CellDim], float],
-    # dsl specific additional args
-    mean_cell_area: float,
+    vct_a: Field[[KHalfDim], float64],
+    theta_ref_mc: Field[[CellDim, KDim], float64],
+    wgtfac_c: Field[[CellDim, KHalfDim], float64],
+    e_bln_c_s: Field[[CellDim, C2EDim], float64],
+    geofac_div: Field[[CellDim, C2EDim], float64],
+    geofac_grg_x: Field[[CellDim, C2E2CODim], float64],
+    geofac_grg_y: Field[[CellDim, C2E2CODim], float64],
+    geofac_n2s: Field[[CellDim, C2E2CODim], float64],
+    nudgecoeff_e: Field[[EdgeDim], float64],
+    rbf_coeff_1: Field[[VertexDim, V2EDim], float64],
+    rbf_coeff_2: Field[[VertexDim, V2EDim], float64],
     mask_hdiff: Field[[CellDim, KDim], bool],
-    zd_vertoffset: Field[
-        [CECDim], int32
-    ],  # vertical offsets used in DSL for absolute indices zd_vertidx in mo_vertical_grid.f90
-    rbf_coeff_1: Field[[VertexDim, V2EDim], float],  # -> used in rbf_vec_interpol_vertex
-    rbf_coeff_2: Field[[VertexDim, V2EDim], float],  # -> used in rbf_vec_interpol_vertex
-    verts_edge_idx: Field[[VertexDim, V2EDim], int32],  # -> mo_intp_rbf_rbf_vec_interpol_vertex
+    zd_diffcoef: Field[[CellDim, KDim], float64],
+    zd_vertoffset: Field[[CellDim, C2E2CDim, KDim], int32],
+    zd_intcoef: Field[[CellDim, C2E2CDim, KDim], float64],
+    num_levels: int32,
+    mean_cell_area: float64,
+    ndyn_substeps: int32,
+    rayleigh_damping_height: float64,
+    nflatlev: int32,
+    nflat_gradp: int32,
+    diffusion_type: int32,
+    hdiff_w: bool,
+    hdiff_vn: bool,
+    zdiffu_t: bool,
+    type_t_diffu: int32,
+    type_vn_diffu: int32,
+    hdiff_efdt_ratio: float64,
+    smagorinski_scaling_factor: float64,
+    hdiff_temp: bool,
+    tangent_orientation: Field[[EdgeDim], float64],
+    inverse_primal_edge_lengths: Field[[EdgeDim], float64],
+    inv_dual_edge_length: Field[[EdgeDim], float64],
+    inv_vert_vert_length: Field[[EdgeDim], float64],
+    edge_areas: Field[[EdgeDim], float64],
+    f_e: Field[[EdgeDim], float64],
+    cell_areas: Field[[CellDim], float64],
+    primal_normal_vert_x: Field[[EdgeDim, E2C2VDim], float64],
+    primal_normal_vert_y: Field[[EdgeDim, E2C2VDim], float64],
+    dual_normal_vert_x: Field[[EdgeDim, E2C2VDim], float64],
+    dual_normal_vert_y: Field[[EdgeDim, E2C2VDim], float64],
+    primal_normal_cell_x: Field[[EdgeDim, E2CDim], float64],
+    primal_normal_cell_y: Field[[EdgeDim, E2CDim], float64],
+    dual_normal_cell_x: Field[[EdgeDim, E2CDim], float64],
+    dual_normal_cell_y: Field[[EdgeDim, E2CDim], float64],
 ):
-    """
-    Instantiate and Initialize the diffusion object.
-
-    should only accept simple fields as arguments for compatibility with the standalone
-    Fortran ICON Diffusion component (aka Diffusion granule)
+    logger.info(f"Using Device = {device}")
 
-    """
-    if diffusion.initialized:
-        raise DuplicateInitializationException("Diffusion has already been already initialized")
+    # ICON grid
+    if device.name == "GPU":
+        on_gpu = True
+    else:
+        on_gpu = False
 
-    horizontal_size = HorizontalGridSize(num_cells=nproma, num_vertices=nproma, num_edges=nproma)
-    vertical_size = VerticalGridSize(num_lev=nlev)
-    mesh_config = GridConfig(
-        horizontal_config=horizontal_size,
-        vertical_config=vertical_size,
-        limited_area=l_limited_area,
-        n_shift_total=n_shift_total,
+    icon_grid = _load_from_gridfile(
+        file_path=get_icon_grid_loc(),
+        filename=get_grid_filename(),
+        num_levels=num_levels,
+        on_gpu=on_gpu,
     )
-    # we need the start, end indices in order for the grid to be functional those are not passed
-    # to init in the Fortran diffusion granule since they are hidden away in the get_indices_[c,e,v]
-    # diffusion_run does take p_patch in order to pass it on to other subroutines (interpolation, get_indices...
 
-    c2e2c0 = np.column_stack(
-        (
-            (np.asarray(cells_neighbor_idx)),
-            (np.asarray(range(np.asarray(cells_neighbor_idx).shape[0]))),
-        )
-    )
-    e2c2v = np.asarray(edges_vertex_idx)
-    e2v = e2c2v[:, 0:2]
-    connectivities = {
-        E2CDim: np.asarray(edges_cell_idx),
-        E2C2VDim: e2c2v,
-        E2VDim: e2v,
-        C2EDim: np.asarray(cells_edge_idx),
-        C2E2CDim: (np.asarray(cells_neighbor_idx)),
-        C2E2CODim: c2e2c0,
-        V2EDim: np.asarray(verts_edge_idx),
-    }
-    # TODO (Magdalena) we need start_index, end_index: those are not passed in the fortran granules,
-    #  because they are used through get_indices only
-    grid = IconGrid().with_config(mesh_config).with_connectivities(connectivities)
+    # Edge geometry
     edge_params = EdgeParams(
-        tangent_orientation=edges_tangent_orientation,
-        inverse_primal_edge_lengths=edges_inv_primal_edge_length,
-        dual_normal_vert_x=edges_dual_normal_vert_1,
-        dual_normal_vert_y=edges_dual_normal_vert_2,
-        inverse_dual_edge_lengths=edges_inv_dual_edge_length,
-        inverse_vertex_vertex_lengths=edges_inv_vert_vert_lengths,
-        primal_normal_vert_x=edges_primal_normal_vert_1,
-        primal_normal_vert_y=edges_primal_normal_vert_2,
-        edge_areas=edges_area_edge,
+        tangent_orientation=tangent_orientation,
+        inverse_primal_edge_lengths=inverse_primal_edge_lengths,
+        inverse_dual_edge_lengths=inv_dual_edge_length,
+        inverse_vertex_vertex_lengths=inv_vert_vert_length,
+        primal_normal_vert_x=as_1D_sparse_field(primal_normal_vert_x, ECVDim),
+        primal_normal_vert_y=as_1D_sparse_field(primal_normal_vert_y, ECVDim),
+        dual_normal_vert_x=as_1D_sparse_field(dual_normal_vert_x, ECVDim),
+        dual_normal_vert_y=as_1D_sparse_field(dual_normal_vert_y, ECVDim),
+        primal_normal_cell_x=as_1D_sparse_field(primal_normal_cell_x, ECVDim),
+        primal_normal_cell_y=as_1D_sparse_field(primal_normal_cell_y, ECVDim),
+        dual_normal_cell_x=as_1D_sparse_field(dual_normal_cell_x, ECVDim),
+        dual_normal_cell_y=as_1D_sparse_field(dual_normal_cell_y, ECVDim),
+        edge_areas=edge_areas,
+        f_e=f_e,
     )
-    cell_params = CellParams(area=cells_area, mean_cell_area=mean_cell_area)
-    config: DiffusionConfig = DiffusionConfig(
-        diffusion_type=DiffusionType(hdiff_order),
-        hdiff_w=lhdiff_w,
-        hdiff_temp=lhdiff_temp,
-        type_vn_diffu=itype_vn_diffu,
-        smag_3d=lsmag_3d,
-        type_t_diffu=itype_t_diffu,
+
+    # cell geometry
+    cell_params = CellParams(area=cell_areas, mean_cell_area=mean_cell_area)
+
+    # diffusion parameters
+    config = DiffusionConfig(
+        diffusion_type=diffusion_type,
+        hdiff_w=hdiff_w,
+        hdiff_vn=hdiff_vn,
+        zdiffu_t=zdiffu_t,
+        type_t_diffu=type_t_diffu,
+        type_vn_diffu=type_vn_diffu,
         hdiff_efdt_ratio=hdiff_efdt_ratio,
-        hdiff_w_efdt_ratio=hdiff_efdt_ratio,
-        smagorinski_scaling_factor=hdiff_smag_fac,
-        n_substeps=n_dyn_substeps,
-        zdiffu_t=l_zdiffu_t,
-        velocity_boundary_diffusion_denom=denom_diffu_v,
-        max_nudging_coeff=nudge_max_coeff,
+        smagorinski_scaling_factor=smagorinski_scaling_factor,
+        hdiff_temp=hdiff_temp,
+        n_substeps=ndyn_substeps,
+        thslp_zdiffu=0.02,
+        thhgtd_zdiffu=125.0,
+        velocity_boundary_diffusion_denom=150.0,
+        max_nudging_coeff=0.075,
+        shear_type=TurbulenceShearForcingType.VERTICAL_HORIZONTAL_OF_HORIZONTAL_VERTICAL_WIND,
+    )
+
+    diffusion_params = DiffusionParams(config)
+
+    # vertical parameters
+    vertical_params = VerticalModelParams(
+        vct_a=vct_a,
+        rayleigh_damping_height=rayleigh_damping_height,
+        nflatlev=nflatlev,
+        nflat_gradp=nflat_gradp,
     )
-    vertical_params = VerticalModelParams(vct_a=vct_a, rayleigh_damping_height=nrdmax)
 
-    derived_diffusion_params = DiffusionParams(config)
+    # metric state
     metric_state = DiffusionMetricState(
+        mask_hdiff=mask_hdiff,
         theta_ref_mc=theta_ref_mc,
         wgtfac_c=wgtfac_c,
-        mask_hdiff=mask_hdiff,
-        zd_vertoffset=zd_vertoffset,
+        zd_intcoef=flatten_first_two_dims(CECDim, KDim, field=zd_intcoef),
+        zd_vertoffset=flatten_first_two_dims(CECDim, KDim, field=zd_vertoffset),
         zd_diffcoef=zd_diffcoef,
-        zd_intcoef=zd_intcoef,
     )
+
+    # interpolation state
     interpolation_state = DiffusionInterpolationState(
-        e_bln_c_s,
-        rbf_coeff_1,
-        rbf_coeff_2,
-        geofac_div,
-        geofac_n2s,
-        geofac_grg_x,
-        geofac_grg_y,
-        nudgecoeff_e,
+        e_bln_c_s=as_1D_sparse_field(e_bln_c_s, CEDim),
+        rbf_coeff_1=rbf_coeff_1,
+        rbf_coeff_2=rbf_coeff_2,
+        geofac_div=as_1D_sparse_field(geofac_div, CEDim),
+        geofac_n2s=geofac_n2s,
+        geofac_grg_x=geofac_grg_x,
+        geofac_grg_y=geofac_grg_y,
+        nudgecoeff_e=nudgecoeff_e,
     )
-
-    diffusion.init(
-        grid=grid,
+    diffusion_granule.init(
+        grid=icon_grid,
         config=config,
-        params=derived_diffusion_params,
+        params=diffusion_params,
         vertical_params=vertical_params,
         metric_state=metric_state,
         interpolation_state=interpolation_state,
@@ -229,37 +213,41 @@ def diffusion_init(
 
 
 def diffusion_run(
-    dtime: float,
+    w: Field[[CellDim, KHalfDim], float64],
+    vn: Field[[EdgeDim, KDim], float64],
+    exner: Field[[CellDim, KDim], float64],
+    theta_v: Field[[CellDim, KDim], float64],
+    rho: Field[[CellDim, KDim], float64],
+    hdef_ic: Field[[CellDim, KHalfDim], float64],
+    div_ic: Field[[CellDim, KHalfDim], float64],
+    dwdx: Field[[CellDim, KHalfDim], float64],
+    dwdy: Field[[CellDim, KHalfDim], float64],
+    dtime: float64,
     linit: bool,
-    vn: Field[[EdgeDim, KDim], float],
-    w: Field[[CellDim, KDim], float],
-    theta_v: Field[[CellDim, KDim], float],
-    exner: Field[[CellDim, KDim], float],
-    div_ic: Field[[CellDim, KDim], float],
-    hdef_ic: Field[[CellDim, KDim], float],
-    dwdx: Field[[CellDim, KDim], float],
-    dwdy: Field[[CellDim, KDim], float],
-    rho: Field[[CellDim, KDim], float],
 ):
-    diagnostic_state = DiffusionDiagnosticState(hdef_ic, div_ic, dwdx, dwdy)
+    # prognostic and diagnostic variables
     prognostic_state = PrognosticState(
-        rho=rho,
         w=w,
         vn=vn,
         exner=exner,
         theta_v=theta_v,
+        rho=rho,
     )
+
+    diagnostic_state = DiffusionDiagnosticState(
+        hdef_ic=hdef_ic,
+        div_ic=div_ic,
+        dwdx=dwdx,
+        dwdy=dwdy,
+    )
+
     if linit:
-        diffusion.initial_run(
+        diffusion_granule.initial_run(
             diagnostic_state,
             prognostic_state,
             dtime,
         )
     else:
-        diffusion.run(diagnostic_state, prognostic_state, dtime)
-
-
-class DuplicateInitializationException(Exception):
-    """Raised if the component is already initilalized."""
-
-    pass
+        diffusion_granule.run(
+            prognostic_state=prognostic_state, diagnostic_state=diagnostic_state, dtime=dtime
+        )
diff --git a/tools/src/icon4pytools/py2fgen/wrappers/simple.py b/tools/src/icon4pytools/py2fgen/wrappers/simple.py
index 0b66da34b5..f7e072ff02 100644
--- a/tools/src/icon4pytools/py2fgen/wrappers/simple.py
+++ b/tools/src/icon4pytools/py2fgen/wrappers/simple.py
@@ -78,3 +78,10 @@ def multi_return(
             KDim: (vertical_start, vertical_end),
         },
     )
+
+
+def square_error(
+    inp: Field[[CellDim, KDim], float64],
+    result: Field[[CellDim, KDim], float64],
+):
+    raise Exception("Exception foo occurred")
diff --git a/tools/tests/py2fgen/fortran_samples/test_diffusion.f90 b/tools/tests/py2fgen/fortran_samples/test_diffusion.f90
new file mode 100644
index 0000000000..460de43512
--- /dev/null
+++ b/tools/tests/py2fgen/fortran_samples/test_diffusion.f90
@@ -0,0 +1,331 @@
+module random_utils
+contains
+   subroutine fill_random_1d(array, low, high)
+      use, intrinsic :: iso_c_binding, only: c_double
+      implicit none
+
+      real(c_double), intent(inout) :: array(:)
+      real(c_double), intent(in) :: low, high
+      integer :: i
+      real(c_double) :: rnd
+
+      do i = 1, size(array)
+         call random_number(rnd)
+         array(i) = low + rnd*(high - low)
+      end do
+   end subroutine fill_random_1d
+
+   subroutine fill_random_2d(array, low, high)
+      use, intrinsic :: iso_c_binding, only: c_double
+      implicit none
+
+      real(c_double), intent(inout) :: array(:, :)
+      real(c_double), intent(in) :: low, high
+      integer :: i, j
+      real(c_double) :: rnd
+
+      do i = 1, size(array, 1)
+         do j = 1, size(array, 2)
+            call random_number(rnd)
+            array(i, j) = low + rnd*(high - low)
+         end do
+      end do
+   end subroutine fill_random_2d
+
+    subroutine fill_random_3d_int(array, low, high)
+    use, intrinsic :: iso_c_binding, only: c_int
+    implicit none
+
+    integer(c_int), intent(inout) :: array(:, :, :)
+    integer(c_int), intent(in) :: low, high
+    integer :: i, j, k
+    real :: rnd  ! real number between 0 and 1
+
+    do i = 1, size(array, 1)
+        do j = 1, size(array, 2)
+            do k = 1, size(array, 3)
+                call random_number(rnd)
+                array(i, j, k) = int(low + rnd * (high - low))
+            end do
+        end do
+    end do
+    end subroutine fill_random_3d_int
+
+    subroutine fill_random_3d(array, low, high)
+    use, intrinsic :: iso_c_binding, only: c_double
+    implicit none
+
+    real(c_double), intent(inout) :: array(:, :, :)
+    real(c_double), intent(in) :: low, high
+    integer :: i, j, k
+    real :: rnd  ! real number between 0 and 1
+
+    do i = 1, size(array, 1)
+        do j = 1, size(array, 2)
+            do k = 1, size(array, 3)
+                call random_number(rnd)
+                array(i, j, k) = low + rnd * (high - low)
+            end do
+        end do
+    end do
+    end subroutine fill_random_3d
+
+    subroutine fill_random_2d_bool(array)
+    use, intrinsic :: iso_c_binding, only: c_int
+    implicit none
+
+    logical(c_int), intent(inout) :: array(:, :)
+    integer :: i, j
+    real :: rnd  ! real number between 0 and 1
+
+    do i = 1, size(array, 1)
+        do j = 1, size(array, 2)
+            call random_number(rnd)
+            if (rnd < 0.5) then
+                array(i, j) = .false.
+            else
+                array(i, j) = .true.
+            endif
+        end do
+    end do
+end subroutine fill_random_2d_bool
+end module random_utils
+
+program diffusion_simulation
+   use, intrinsic :: iso_c_binding
+   use random_utils, only: fill_random_1d, fill_random_2d, fill_random_2d_bool, fill_random_3d_int, fill_random_3d
+   use diffusion_plugin
+   implicit none
+
+   integer(c_int) :: rc
+
+   ! Constants and types
+   integer(c_int), parameter :: num_cells = 20480
+   integer(c_int), parameter :: num_edges = 30720
+   integer(c_int), parameter :: num_vertices = 10242
+   integer(c_int), parameter :: num_levels = 60
+   integer(c_int), parameter :: num_c2ec2o = 4
+   integer(c_int), parameter :: num_v2e = 6
+   integer(c_int), parameter :: num_c2e = 2
+   integer(c_int), parameter :: num_e2c2v = 4
+   integer(c_int), parameter :: num_c2e2c = 3
+   integer(c_int), parameter :: num_e2c = 3
+   real(c_double), parameter :: mean_cell_area = 24907282236.708576
+   integer(c_int), parameter :: ndyn_substeps = 2
+   real(c_double), parameter :: dtime = 2.0
+   real(c_double), parameter :: rayleigh_damping_height = 50000
+   integer(c_int), parameter :: nflatlev = 30
+   integer(c_int), parameter :: nflat_gradp = 59
+   integer(c_int), parameter :: diffusion_type = 5 ! Assuming DiffusionType.SMAGORINSKY_4TH_ORDER is represented by 5
+   logical(c_int), parameter :: hdiff_w = .true.
+   logical(c_int), parameter :: hdiff_vn = .true.
+   logical(c_int), parameter :: zdiffu_t = .true. ! this runs stencil 15 which uses the boolean mask
+   integer(c_int), parameter :: type_t_diffu = 2
+   integer(c_int), parameter :: type_vn_diffu = 1
+   real(c_double), parameter :: hdiff_efdt_ratio = 24.0
+   real(c_double), parameter :: smagorinski_scaling_factor = 0.025
+   logical(c_int), parameter :: hdiff_temp = .true.
+   logical(c_int), parameter :: linit = .true.
+
+   ! Declaring arrays for diffusion_init and diffusion_run
+   real(c_double), dimension(:), allocatable :: vct_a
+   real(c_double), dimension(:), allocatable :: nudgecoeff_e
+   real(c_double), dimension(:), allocatable :: tangent_orientation
+   real(c_double), dimension(:), allocatable :: inverse_primal_edge_lengths
+   real(c_double), dimension(:), allocatable :: inv_dual_edge_length
+   real(c_double), dimension(:), allocatable :: inv_vert_vert_length
+   real(c_double), dimension(:), allocatable :: edge_areas
+   real(c_double), dimension(:), allocatable :: f_e
+   real(c_double), dimension(:), allocatable :: cell_areas
+
+   real(c_double), dimension(:, :), allocatable :: theta_ref_mc
+   real(c_double), dimension(:, :), allocatable :: wgtfac_c
+   real(c_double), dimension(:, :), allocatable :: e_bln_c_s
+   real(c_double), dimension(:, :), allocatable :: geofac_div
+   real(c_double), dimension(:, :), allocatable :: geofac_grg_x
+   real(c_double), dimension(:, :), allocatable :: geofac_grg_y
+   real(c_double), dimension(:, :), allocatable :: geofac_n2s
+   real(c_double), dimension(:, :), allocatable :: rbf_coeff_1
+   real(c_double), dimension(:, :), allocatable :: rbf_coeff_2
+   real(c_double), dimension(:, :), allocatable :: dwdx
+   real(c_double), dimension(:, :), allocatable :: dwdy
+   real(c_double), dimension(:, :), allocatable :: hdef_ic
+   real(c_double), dimension(:, :), allocatable :: div_ic
+   real(c_double), dimension(:, :), allocatable :: w
+   real(c_double), dimension(:, :), allocatable :: vn
+   real(c_double), dimension(:, :), allocatable :: exner
+   real(c_double), dimension(:, :), allocatable :: theta_v
+   real(c_double), dimension(:, :), allocatable :: rho
+   real(c_double), dimension(:, :), allocatable :: dual_normal_cell_x
+   real(c_double), dimension(:, :), allocatable :: dual_normal_cell_y
+   real(c_double), dimension(:, :), allocatable :: dual_normal_vert_x
+   real(c_double), dimension(:, :), allocatable :: dual_normal_vert_y
+   real(c_double), dimension(:, :), allocatable :: primal_normal_cell_x
+   real(c_double), dimension(:, :), allocatable :: primal_normal_cell_y
+   real(c_double), dimension(:, :), allocatable :: primal_normal_vert_x
+   real(c_double), dimension(:, :), allocatable :: primal_normal_vert_y
+   real(c_double), dimension(:, :), allocatable :: zd_diffcoef
+   logical(c_int), dimension(:, :), allocatable :: mask_hdiff
+
+   integer(c_int), dimension(:, :, :), allocatable :: zd_vertoffset
+   real(c_double), dimension(:, :, :), allocatable :: zd_intcoef
+
+    !$acc enter data create(vct_a, theta_ref_mc, wgtfac_c, e_bln_c_s, geofac_div, &
+    !$acc geofac_grg_x, geofac_grg_y, geofac_n2s, nudgecoeff_e, rbf_coeff_1, &
+    !$acc rbf_coeff_2, dwdx, dwdy, hdef_ic, div_ic, w, vn, exner, theta_v, rho, &
+    !$acc dual_normal_cell_x, dual_normal_cell_y, dual_normal_vert_x, &
+    !$acc dual_normal_vert_y, primal_normal_cell_x, primal_normal_cell_y, &
+    !$acc primal_normal_vert_x, primal_normal_vert_y, tangent_orientation, &
+    !$acc inverse_primal_edge_lengths, inv_dual_edge_length, inv_vert_vert_length, &
+    !$acc edge_areas, cell_areas, f_e, zd_diffcoef, zd_vertoffset, zd_intcoef, &
+    !$acc mask_hdiff)
+
+   ! allocating arrays
+   allocate(zd_diffcoef(num_cells, num_levels))
+   allocate(zd_vertoffset(num_cells, num_c2e2c, num_levels))
+   allocate(zd_intcoef(num_cells, num_c2e2c, num_levels))
+   allocate(mask_hdiff(num_cells, num_levels))
+   allocate (vct_a(num_levels))
+   allocate (theta_ref_mc(num_cells, num_levels))
+   allocate (wgtfac_c(num_cells, num_levels + 1))
+   allocate (e_bln_c_s(num_cells, num_c2e))
+   allocate (geofac_div(num_cells, num_c2e))
+   allocate (geofac_grg_x(num_cells, num_c2ec2o))
+   allocate (geofac_grg_y(num_cells, num_c2ec2o))
+   allocate (geofac_n2s(num_cells, num_c2ec2o))
+   allocate (nudgecoeff_e(num_edges))
+   allocate (rbf_coeff_1(num_vertices, num_v2e))
+   allocate (rbf_coeff_2(num_vertices, num_v2e))
+   allocate (dwdx(num_cells, num_levels))
+   allocate (dwdy(num_cells, num_levels))
+   allocate (hdef_ic(num_cells, num_levels + 1))
+   allocate (div_ic(num_cells, num_levels + 1))
+   allocate (w(num_cells, num_levels + 1))
+   allocate (vn(num_edges, num_levels))
+   allocate (exner(num_cells, num_levels))
+   allocate (theta_v(num_cells, num_levels))
+   allocate (rho(num_cells, num_levels))
+   allocate (dual_normal_cell_x(num_edges, num_e2c))
+   allocate (dual_normal_cell_y(num_edges, num_e2c))
+   allocate (dual_normal_vert_x(num_edges, num_e2c2v))
+   allocate (dual_normal_vert_y(num_edges, num_e2c2v))
+   allocate (primal_normal_cell_x(num_edges, num_e2c))
+   allocate (primal_normal_cell_y(num_edges, num_e2c))
+   allocate (primal_normal_vert_x(num_edges, num_e2c))
+   allocate (primal_normal_vert_y(num_edges, num_e2c))
+   allocate (tangent_orientation(num_edges))
+   allocate (inverse_primal_edge_lengths(num_edges))
+   allocate (inv_dual_edge_length(num_edges))
+   allocate (inv_vert_vert_length(num_edges))
+   allocate (edge_areas(num_edges))
+   allocate (f_e(num_edges))
+   allocate (cell_areas(num_cells))
+
+   ! Fill arrays with random numbers
+   ! For 1D arrays
+   call fill_random_1d(vct_a, 0.0_c_double, 75000.0_c_double) ! needs to be above 50000 damping height restriction
+   call fill_random_1d(nudgecoeff_e, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(tangent_orientation, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(inverse_primal_edge_lengths, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(inv_dual_edge_length, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(inv_vert_vert_length, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(edge_areas, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(f_e, 0.0_c_double, 1.0_c_double)
+   call fill_random_1d(cell_areas, 0.0_c_double, 1.0_c_double)
+
+   ! For 2D arrays
+   call fill_random_2d(theta_ref_mc, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(wgtfac_c, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(geofac_grg_x, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(geofac_grg_y, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(geofac_n2s, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(rbf_coeff_1, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(rbf_coeff_2, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dwdx, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dwdy, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(hdef_ic, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(div_ic, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(w, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(vn, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(exner, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(theta_v, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(rho, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(zd_diffcoef, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d_bool(mask_hdiff)
+
+   call fill_random_2d(e_bln_c_s, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(geofac_div, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dual_normal_vert_x, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dual_normal_vert_y, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(primal_normal_vert_x, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(primal_normal_vert_y, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dual_normal_cell_x, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(dual_normal_cell_y, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(primal_normal_cell_x, 0.0_c_double, 1.0_c_double)
+   call fill_random_2d(primal_normal_cell_y, 0.0_c_double, 1.0_c_double)
+
+   ! For 3D arrays
+   call fill_random_3d(zd_intcoef, 0.0_c_double, 1.0_c_double)
+   call fill_random_3d_int(zd_vertoffset, 0, 1)
+
+
+    !$acc data copyin(vct_a, theta_ref_mc, wgtfac_c, e_bln_c_s, geofac_div, &
+    !$acc geofac_grg_x, geofac_grg_y, geofac_n2s, nudgecoeff_e, rbf_coeff_1, &
+    !$acc rbf_coeff_2, dwdx, dwdy, hdef_ic, div_ic, w, vn, exner, theta_v, rho, &
+    !$acc dual_normal_cell_x, dual_normal_cell_y, dual_normal_vert_x, &
+    !$acc dual_normal_vert_y, primal_normal_cell_x, primal_normal_cell_y, &
+    !$acc primal_normal_vert_x, primal_normal_vert_y, tangent_orientation, &
+    !$acc inverse_primal_edge_lengths, inv_dual_edge_length, inv_vert_vert_length, &
+    !$acc edge_areas, cell_areas, f_e, zd_diffcoef, zd_vertoffset, zd_intcoef, &
+    !$acc mask_hdiff)
+
+   ! Call diffusion_init
+   call diffusion_init(vct_a, theta_ref_mc, wgtfac_c, e_bln_c_s, geofac_div, &
+                       geofac_grg_x, geofac_grg_y, geofac_n2s, nudgecoeff_e, rbf_coeff_1, &
+                       rbf_coeff_2, mask_hdiff, zd_diffcoef, zd_vertoffset, zd_intcoef, &
+                       num_levels, mean_cell_area, ndyn_substeps, rayleigh_damping_height, &
+                       nflatlev, nflat_gradp, diffusion_type, &
+                       hdiff_w, hdiff_vn, zdiffu_t, type_t_diffu, type_vn_diffu, &
+                       hdiff_efdt_ratio, smagorinski_scaling_factor, hdiff_temp, &
+                       tangent_orientation, inverse_primal_edge_lengths, inv_dual_edge_length, &
+                       inv_vert_vert_length, edge_areas, f_e, cell_areas, primal_normal_vert_x, &
+                       primal_normal_vert_y, dual_normal_vert_x, dual_normal_vert_y, &
+                       primal_normal_cell_x, primal_normal_cell_y, dual_normal_cell_x, &
+                       dual_normal_cell_y, rc)
+
+   print *, "Python exit code = ", rc
+   if (rc /= 0) then
+       call exit(1)
+   end if
+
+   ! Call diffusion_run
+   call diffusion_run(w, vn, exner, theta_v, rho, hdef_ic, div_ic, dwdx, dwdy, dtime, linit, rc)
+
+   print *, "Python exit code = ", rc
+   if (rc /= 0) then
+       call exit(1)
+   end if
+
+    !$acc update host(vct_a, theta_ref_mc, wgtfac_c, e_bln_c_s, geofac_div, &
+    !$acc geofac_grg_x, geofac_grg_y, geofac_n2s, nudgecoeff_e, rbf_coeff_1, &
+    !$acc rbf_coeff_2, dwdx, dwdy, hdef_ic, div_ic, w, vn, exner, theta_v, rho, &
+    !$acc dual_normal_cell_x, dual_normal_cell_y, dual_normal_vert_x, &
+    !$acc dual_normal_vert_y, primal_normal_cell_x, primal_normal_cell_y, &
+    !$acc primal_normal_vert_x, primal_normal_vert_y, tangent_orientation, &
+    !$acc inverse_primal_edge_lengths, inv_dual_edge_length, inv_vert_vert_length, &
+    !$acc edge_areas, cell_areas, f_e, zd_diffcoef, zd_vertoffset, zd_intcoef, &
+    !$acc mask_hdiff)
+
+   print *, "passed: could run diffusion"
+
+    !$acc end data
+    !$acc exit data delete(vct_a, theta_ref_mc, wgtfac_c, e_bln_c_s, geofac_div, &
+    !$acc geofac_grg_x, geofac_grg_y, geofac_n2s, nudgecoeff_e, rbf_coeff_1, &
+    !$acc rbf_coeff_2, dwdx, dwdy, hdef_ic, div_ic, w, vn, exner, theta_v, rho, &
+    !$acc dual_normal_cell_x, dual_normal_cell_y, dual_normal_vert_x, &
+    !$acc dual_normal_vert_y, primal_normal_cell_x, primal_normal_cell_y, &
+    !$acc primal_normal_vert_x, primal_normal_vert_y, tangent_orientation, &
+    !$acc inverse_primal_edge_lengths, inv_dual_edge_length, inv_vert_vert_length, &
+    !$acc edge_areas, cell_areas, f_e, zd_diffcoef, zd_vertoffset, zd_intcoef, &
+    !$acc mask_hdiff)
+end program diffusion_simulation
diff --git a/tools/tests/py2fgen/fortran_samples/test_multi_return.f90 b/tools/tests/py2fgen/fortran_samples/test_multi_return.f90
index c0b324dd4c..86c13b7495 100644
--- a/tools/tests/py2fgen/fortran_samples/test_multi_return.f90
+++ b/tools/tests/py2fgen/fortran_samples/test_multi_return.f90
@@ -1,83 +1,96 @@
 program call_multi_return_cffi_plugin
-    use, intrinsic :: iso_c_binding
-    use multi_return_plugin
-    implicit none
+   use, intrinsic :: iso_c_binding
+   use multi_return_plugin
+   implicit none
 
-    integer(c_int) :: edim, kdim, i, j,  horizontal_start, horizontal_end, vertical_start, vertical_end
-    logical :: computation_correct
-    character(len=100) :: str_buffer
-    real(c_double) :: r_nsubsteps
-    real(c_double), dimension(:,:), allocatable :: z_vn_avg, mass_fl_e, vn_traj, mass_flx_me
+   integer(c_int) :: edim, kdim, i, j, horizontal_start, horizontal_end, vertical_start, vertical_end, rc
+   logical :: computation_correct
+   character(len=100) :: str_buffer
+   real(c_double) :: r_nsubsteps
+   real(c_double), dimension(:, :), allocatable :: z_vn_avg, mass_fl_e, vn_traj, mass_flx_me
 
-    ! array dimensions
-    edim = 27
-    kdim = 10
+   ! array dimensions
+   edim = 27
+   kdim = 10
 
-    ! allocate arrays (allocate in column-major order)
-    allocate(z_vn_avg(kdim, edim))
-    allocate(mass_fl_e(kdim, edim))
-    allocate(vn_traj(kdim, edim))
-    allocate(mass_flx_me(kdim, edim))
+   !$ACC enter data create(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me)
 
-    ! initialize arrays and variables
-    z_vn_avg = 1.0d0
-    mass_fl_e = 2.0d0
-    vn_traj = 3.0d0
-    mass_flx_me = 4.0d0
-    r_nsubsteps = 9.0d0
-    horizontal_start = 0
-    horizontal_end = edim
-    vertical_start = 0
-    vertical_end = kdim
+   ! allocate arrays (allocate in column-major order)
+   allocate (z_vn_avg(edim, kdim))
+   allocate (mass_fl_e(edim, kdim))
+   allocate (vn_traj(edim, kdim))
+   allocate (mass_flx_me(edim, kdim))
 
-    ! print array shapes and values before computation
-    print *, "Arrays before computation:"
-    write(str_buffer, '("Shape of z_vn_avg = ", I2, ",", I2)') size(z_vn_avg, 1), size(z_vn_avg, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of mass_fl_e = ", I2, ",", I2)') size(mass_fl_e, 1), size(mass_fl_e, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of vn_traj = ", I2, ",", I2)') size(vn_traj, 1), size(vn_traj, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of mass_flx_me = ", I2, ",", I2)') size(mass_flx_me, 1), size(mass_flx_me, 2)
-    print *, trim(str_buffer)
-    print *
+   ! initialize arrays and variables
+   z_vn_avg = 1.0d0
+   mass_fl_e = 2.0d0
+   vn_traj = 3.0d0
+   mass_flx_me = 4.0d0
+   r_nsubsteps = 9.0d0
+   horizontal_start = 0
+   horizontal_end = edim
+   vertical_start = 0
+   vertical_end = kdim
 
-    ! call the cffi plugin
-    call multi_return_wrapper(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me, r_nsubsteps, &
-            horizontal_start, horizontal_end, vertical_start, vertical_end, edim, kdim)
+   !$ACC data copyin(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me, r_nsubsteps)
 
-    ! print array shapes and values before computation
-    print *, "Arrays after computation:"
-    write(str_buffer, '("Shape of z_vn_avg = ", I2, ",", I2)') size(z_vn_avg, 1), size(z_vn_avg, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of mass_fl_e = ", I2, ",", I2)') size(mass_fl_e, 1), size(mass_fl_e, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of vn_traj = ", I2, ",", I2)') size(vn_traj, 1), size(vn_traj, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of mass_flx_me = ", I2, ",", I2)') size(mass_flx_me, 1), size(mass_flx_me, 2)
-    print *, trim(str_buffer)
-    print *
+   ! print array shapes and values before computation
+   print *, "Arrays before computation:"
+   write (str_buffer, '("Shape of z_vn_avg = ", I2, ",", I2)') size(z_vn_avg, 1), size(z_vn_avg, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of mass_fl_e = ", I2, ",", I2)') size(mass_fl_e, 1), size(mass_fl_e, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of vn_traj = ", I2, ",", I2)') size(vn_traj, 1), size(vn_traj, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of mass_flx_me = ", I2, ",", I2)') size(mass_flx_me, 1), size(mass_flx_me, 2)
+   print *, trim(str_buffer)
+   print *
 
-    ! Assert vn_traj == 12 and mass_flx_me == 22
-    computation_correct = .true.
-    do i = 1, edim
-        do j = 1, kdim
-            if (vn_traj(i, j) /= 12.0d0 .or. mass_flx_me(i, j) /= 22.0d0) then
-                computation_correct = .false.
-                exit
-            end if
-        end do
-        if (.not. computation_correct) exit
-    end do
+   ! call the cffi plugin
+   call multi_return(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me, r_nsubsteps, &
+                     horizontal_start, horizontal_end, vertical_start, vertical_end, rc)
+   print *, "Python exit code = ", rc
+   if (rc /= 0) then
+       call exit(1)
+   end if
 
-    ! deallocate arrays
-    deallocate(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me)
+   !$ACC update host(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me)
 
-    ! Check and print the result of the assertion
-    if (computation_correct) then
-        print *, "passed: vn_traj and mass_flx_me have expected values."
-    else
-        print *, "failed: vn_traj or mass_flx_me does not have the expected values."
-        stop 1
-    end if
+   ! print array shapes and values before computation
+   print *, "Arrays after computation:"
+   write (str_buffer, '("Shape of z_vn_avg = ", I2, ",", I2)') size(z_vn_avg, 1), size(z_vn_avg, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of mass_fl_e = ", I2, ",", I2)') size(mass_fl_e, 1), size(mass_fl_e, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of vn_traj = ", I2, ",", I2)') size(vn_traj, 1), size(vn_traj, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of mass_flx_me = ", I2, ",", I2)') size(mass_flx_me, 1), size(mass_flx_me, 2)
+   print *, trim(str_buffer)
+   print *
+
+   ! Assert vn_traj == 12 and mass_flx_me == 22
+   computation_correct = .true.
+   do i = 1, edim
+      do j = 1, kdim
+         if (vn_traj(i, j) /= 12.0d0 .or. mass_flx_me(i, j) /= 22.0d0) then
+            computation_correct = .false.
+            exit
+         end if
+      end do
+      if (.not. computation_correct) exit
+   end do
+
+   !$ACC end data
+   !$ACC exit data delete(z_vn_avg, mass_fl_e, vn_traj, mass_flx_me)
+
+   ! deallocate arrays
+   deallocate (z_vn_avg, mass_fl_e, vn_traj, mass_flx_me)
+
+   ! Check and print the result of the assertion
+   if (computation_correct) then
+      print *, "passed: vn_traj and mass_flx_me have expected values."
+   else
+      print *, "failed: vn_traj or mass_flx_me does not have the expected values."
+      stop 1
+   end if
 end program call_multi_return_cffi_plugin
diff --git a/tools/tests/py2fgen/fortran_samples/test_square.f90 b/tools/tests/py2fgen/fortran_samples/test_square.f90
index faf3891262..24bebd0ea0 100644
--- a/tools/tests/py2fgen/fortran_samples/test_square.f90
+++ b/tools/tests/py2fgen/fortran_samples/test_square.f90
@@ -1,80 +1,90 @@
 program call_square_wrapper_cffi_plugin
-    use, intrinsic :: iso_c_binding
-#ifdef USE_SQUARE_FROM_FUNCTION
-    use square_from_function_plugin
-#else
-    use square_plugin
-#endif
-    implicit none
-    character(len=100) :: str_buffer
-    integer(c_int) :: cdim, kdim, i, j
-    logical :: computation_correct
-    real(c_double), dimension(:, :), allocatable :: input, result
+   use, intrinsic :: iso_c_binding
+   use square_plugin
+   implicit none
+   character(len=100) :: str_buffer
+   integer(c_int) :: cdim, kdim, i, j, rc
+   logical :: computation_correct
+   real(c_double), dimension(:, :), allocatable :: input, result
+
+   ! array dimensions
+   cdim = 18
+   kdim = 10
 
-    ! array dimensions
-    cdim = 18
-    kdim = 10
+   !$ACC enter data create(input, result)
 
-    ! allocate arrays (allocate in column-major order)
-    allocate(input(kdim, cdim))
-    allocate(result(kdim, cdim))
+   ! allocate arrays (allocate in column-major order)
+   allocate (input(cdim, kdim))
+   allocate (result(cdim, kdim))
 
-    ! initialise arrays
-    input = 5.0d0
-    result = 0.0d0
+   ! initialise arrays
+   input = 5.0d0
+   result = 0.0d0
 
-    ! print array shapes and values before computation
-    print *, "Arrays before computation:"
-    write(str_buffer, '("Shape of input = ", I2, ",", I2)') size(input, 1), size(input, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of result = ", I2, ",", I2)') size(result, 1), size(result, 2)
-    print *, trim(str_buffer)
-    print *
-    print *, "input = ", input
-    print *
-    print *, "result = ", result
-    print *
+   !$ACC data copyin(input, result)
+
+   ! print array shapes and values before computation
+   print *, "Fortran Arrays before calling Python:"
+   write (str_buffer, '("Shape of input = ", I2, ",", I2)') size(input, 1), size(input, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of result = ", I2, ",", I2)') size(result, 1), size(result, 2)
+   print *, trim(str_buffer)
+   print *
+   print *, "input = ", input
+   print *
+   print *, "result = ", result
+   print *
 
-    ! Call the appropriate cffi plugin
 #ifdef USE_SQUARE_FROM_FUNCTION
-    call square_from_function_wrapper(input, result, cdim, kdim)
+   call square_from_function(input, result, rc)
+#elif USE_SQUARE_ERROR
+   call square_error(input, result, rc)
 #else
-    call square_wrapper(input, result, cdim, kdim)
+   call square(input, result, rc)
 #endif
+   if (rc /= 0) then
+       print *, "Python failed with exit code = ", rc
+       call exit(1)
+   end if
+
+   !$ACC update host(input, result)
+
+   ! print array shapes and values before computation
+   print *, "Fortran arrays after calling Python:"
+   write (str_buffer, '("Shape of input = ", I2, ",", I2)') size(input, 1), size(input, 2)
+   print *, trim(str_buffer)
+   write (str_buffer, '("Shape of result = ", I2, ",", I2)') size(result, 1), size(result, 2)
+   print *, trim(str_buffer)
+   print *
+   print *, "input = ", input
+   print *
+   print *, "result = ", result
+   print *
 
-     ! print array shapes and values before computation
-    print *, "Arrays after computation:"
-    write(str_buffer, '("Shape of input = ", I2, ",", I2)') size(input, 1), size(input, 2)
-    print *, trim(str_buffer)
-    write(str_buffer, '("Shape of result = ", I2, ",", I2)') size(result, 1), size(result, 2)
-    print *, trim(str_buffer)
-    print *
-    print *, "input = ", input
-    print *
-    print *, "result = ", result
-    print *
+   ! Assert each element of result is the square of the corresponding element in input
+   computation_correct = .true.
+   do i = 1, cdim
+      do j = 1, kdim
+         if (result(i, j) /= input(i, j)**2) then
+            print *, "Error: result(", i, ",", j, ") =", result(i, j), &
+               "is not the square of input(", i, ",", j, ") =", input(i, j)
+            computation_correct = .false.
+            exit
+         end if
+      end do
+      if (.not. computation_correct) exit
+   end do
 
-    ! Assert each element of result is the square of the corresponding element in input
-    computation_correct = .true.
-    do i = 1, cdim
-        do j = 1, kdim
-            if (result(i, j) /= input(i, j)**2) then
-                print *, "Error: result(", i, ",", j, ") =", result(i, j), &
-                        "is not the square of input(", i, ",", j, ") =", input(i, j)
-                computation_correct = .false.
-                exit
-            endif
-        enddo
-        if (.not. computation_correct) exit
-    enddo
+   !$ACC end data
+   !$ACC exit data delete(input, result)
 
-    ! deallocate arrays
-    deallocate(input, result)
+   ! deallocate arrays
+   deallocate (input, result)
 
-    ! Check and print the result of the assertion
-    if (computation_correct) then
-        print *, "passed: result has expected values."
-    else
-        print *, "failed: result does not have the expected values."
-    end if
+   ! Check and print the result of the assertion
+   if (computation_correct) then
+      print *, "passed: result has expected values."
+   else
+      print *, "failed: result does not have the expected values."
+   end if
 end program call_square_wrapper_cffi_plugin
diff --git a/tools/tests/py2fgen/test_cffi.py b/tools/tests/py2fgen/test_cffi.py
index fdd16cba74..492c00cd66 100644
--- a/tools/tests/py2fgen/test_cffi.py
+++ b/tools/tests/py2fgen/test_cffi.py
@@ -10,12 +10,15 @@
 # distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
+import subprocess
+import tempfile
+from pathlib import Path
 
 import numpy as np
 import pytest
 from cffi import FFI
 
-from icon4pytools.py2fgen.cffi import unpack
+from icon4pytools.py2fgen.plugin import generate_and_compile_cffi_plugin, unpack
 
 
 @pytest.fixture
@@ -26,11 +29,11 @@ def ffi():
 @pytest.mark.parametrize(
     "data, expected_result",
     [
-        ([1.0, 2.0, 3.0, 4.0], np.array([[1.0, 2.0], [3.0, 4.0]])),
-        ([1, 2, 3, 4], np.array([[1, 2], [3, 4]])),
+        ([1.0, 2.0, 3.0, 4.0], np.array([[1.0, 3.0], [2.0, 4.0]])),
+        ([1, 2, 3, 4], np.array([[1, 3], [2, 4]])),
     ],
 )
-def test_unpack(data, expected_result, ffi):
+def test_unpack_column_major(data, expected_result, ffi):
     ptr = ffi.new("double[]", data) if isinstance(data[0], float) else ffi.new("int[]", data)
 
     rows, cols = expected_result.shape
@@ -38,3 +41,72 @@ def test_unpack(data, expected_result, ffi):
     result = unpack(ptr, rows, cols)
 
     assert np.array_equal(result, expected_result)
+
+
+def test_compile_and_run_cffi_plugin_from_C():
+    plugin_name = "test_plugin"
+    c_header = "int test_function();"
+    c_source_code = f"""
+    #include <stdio.h>
+    #include "{plugin_name}.h"
+
+    int main() {{
+        printf("%d\\n", test_function());
+        return 0;
+    }}
+    """
+
+    python_wrapper = """
+    from test_plugin import ffi
+
+    @ffi.def_extern()
+    def test_function():
+        return 42
+    """
+
+    with tempfile.TemporaryDirectory() as tmpdirname:
+        build_path = Path(tmpdirname)
+
+        try:
+            # Generate and compile the CFFI plugin, which creates lib{plugin_name}.so
+            generate_and_compile_cffi_plugin(plugin_name, c_header, python_wrapper, build_path)
+            compiled_library_path = build_path / f"lib{plugin_name}.so"
+
+            # Verify the shared library was created
+            assert (
+                compiled_library_path.exists()
+            ), f"Compiled library {compiled_library_path} does not exist."
+            assert compiled_library_path.stat().st_size > 0, "Compiled library is empty."
+
+            # Write the main C program to a file
+            main_program_path = build_path / "main.c"
+            with open(main_program_path, "w") as main_program_file:
+                main_program_file.write(c_source_code)
+
+            # Compile the main C program against the shared library
+            subprocess.run(
+                [
+                    "gcc",
+                    "-o",
+                    build_path / "test_program",
+                    str(main_program_path),
+                    "-L" + str(build_path),
+                    "-l" + plugin_name,
+                    "-Wl,-rpath=" + str(build_path),
+                ],
+                check=True,
+            )
+
+            # Execute the compiled program and capture its output
+            result = subprocess.run(
+                [build_path / "test_program"], stdout=subprocess.PIPE, check=True
+            )
+            output = result.stdout.decode().strip()
+
+            # Assert the output of test_function called within the C program
+            assert output == "42", f"Expected '42', got '{output}'"
+
+        except Exception as e:
+            pytest.fail(
+                f"Unexpected error during plugin generation, compilation, or execution: {e}"
+            )
diff --git a/tools/tests/py2fgen/test_cli.py b/tools/tests/py2fgen/test_cli.py
index d40d8c47a7..ce51e8eb30 100644
--- a/tools/tests/py2fgen/test_cli.py
+++ b/tools/tests/py2fgen/test_cli.py
@@ -29,83 +29,117 @@ def wrapper_module():
     return "icon4pytools.py2fgen.wrappers.simple"
 
 
+@pytest.fixture
+def diffusion_module():
+    return "icon4pytools.py2fgen.wrappers.diffusion_test_case"
+
+
 def run_test_case(
     cli,
     module: str,
     function: str,
+    plugin_name: str,
     backend: str,
     samples_path: Path,
     fortran_driver: str,
+    compiler: str = "gfortran",
     extra_compiler_flags: tuple[str, ...] = (),
+    expected_error_code: int = 0,
 ):
     with cli.isolated_filesystem():
-        result = cli.invoke(main, [module, function, "--gt4py-backend", backend, "-d"])
+        result = cli.invoke(main, [module, function, plugin_name, "-b", backend, "-d"])
         assert result.exit_code == 0, "CLI execution failed"
 
         try:
-            compile_fortran_code(function, samples_path, fortran_driver, extra_compiler_flags)
+            compile_fortran_code(
+                plugin_name, samples_path, fortran_driver, compiler, extra_compiler_flags
+            )
         except subprocess.CalledProcessError as e:
             pytest.fail(f"Compilation failed: {e}")
 
         try:
-            fortran_result = run_fortran_executable(function)
-            assert "passed" in fortran_result.stdout
+            fortran_result = run_fortran_executable(plugin_name)
+            if expected_error_code == 0:
+                assert "passed" in fortran_result.stdout
+            else:
+                assert "failed" in fortran_result.stdout
         except subprocess.CalledProcessError as e:
             pytest.fail(f"Execution of compiled Fortran code failed: {e}\nOutput:\n{e.stdout}")
 
 
 def compile_fortran_code(
-    function: str, samples_path: Path, fortran_driver: str, extra_compiler_flags: tuple[str, ...]
+    plugin_name: str,
+    samples_path: Path,
+    fortran_driver: str,
+    compiler: str,
+    extra_compiler_flags: tuple[str, ...],
 ):
-    subprocess.run(["gfortran", "-c", f"{function}_plugin.f90", "."], check=True)
+    command = [
+        f"{compiler}",
+        "-cpp",
+        "-I.",
+        "-Wl,-rpath=.",
+        "-L.",
+        f"{plugin_name}.f90",
+        str(samples_path / f"{fortran_driver}.f90"),
+        f"-l{plugin_name}",
+        "-o",
+        plugin_name,
+    ] + [f for f in extra_compiler_flags]
     subprocess.run(
-        [
-            "gfortran",
-            "-cpp",
-            "-I.",
-            "-Wl,-rpath=.",
-            "-L.",
-            f"{function}_plugin.f90",
-            str(samples_path / f"{fortran_driver}.f90"),
-            f"-l{function}_plugin",
-            "-o",
-            function,
-        ]
-        + [f for f in extra_compiler_flags],
+        command,
         check=True,
     )
 
 
-def run_fortran_executable(function: str):
-    return subprocess.run([f"./{function}"], capture_output=True, text=True, check=True)
-
-
-@pytest.mark.parametrize("backend", ("CPU", "ROUNDTRIP"))
-def test_py2fgen_compilation_and_execution_square(
-    cli_runner, backend, samples_path, wrapper_module
+def run_fortran_executable(plugin_name: str):
+    try:
+        result = subprocess.run([f"./{plugin_name}"], capture_output=True, text=True, check=True)
+    except subprocess.CalledProcessError as e:
+        # If an error occurs, use the exception's `stdout` and `stderr`.
+        result = e
+    return result
+
+
+@pytest.mark.parametrize(
+    "backend, extra_flags",
+    [
+        ("CPU", ("-DUSE_SQUARE_FROM_FUNCTION",)),
+        ("ROUNDTRIP", ""),
+        ("CPU", ("-DUSE_SQUARE_FROM_FUNCTION",)),
+        ("ROUNDTRIP", ""),
+    ],
+)
+def test_py2fgen_compilation_and_execution_square_cpu(
+    cli_runner, backend, samples_path, wrapper_module, extra_flags
 ):
+    """Tests embedding Python functions, and GT4Py program directly.
+    Also tests embedding multiple functions in one shared library.
+    """
     run_test_case(
         cli_runner,
         wrapper_module,
-        "square",
+        "square,square_from_function",
+        "square_plugin",
         backend,
         samples_path,
         "test_square",
+        extra_compiler_flags=extra_flags,
     )
 
 
-@pytest.mark.parametrize("backend", ("CPU", "ROUNDTRIP"))
-def test_py2fgen_compilation_and_execution_square_from_function(
-    cli_runner, backend, samples_path, wrapper_module
-):
+def test_py2fgen_python_error_propagation_to_fortran(cli_runner, samples_path, wrapper_module):
+    """Tests that Exceptions triggered in Python propagate an error code (1) up to Fortran."""
     run_test_case(
         cli_runner,
         wrapper_module,
-        "square_from_function",
-        backend,
+        "square_error",
+        "square_plugin",
+        "ROUNDTRIP",
         samples_path,
         "test_square",
-        ("-DUSE_SQUARE_FROM_FUNCTION",),
+        extra_compiler_flags=("-DUSE_SQUARE_ERROR",),
+        expected_error_code=1,
     )
 
 
@@ -113,11 +147,83 @@ def test_py2fgen_compilation_and_execution_square_from_function(
 def test_py2fgen_compilation_and_execution_multi_return(
     cli_runner, backend, samples_path, wrapper_module
 ):
+    """Tests embedding multi return gt4py program."""
     run_test_case(
         cli_runner,
         wrapper_module,
         "multi_return",
+        "multi_return_plugin",
         backend,
         samples_path,
         "test_multi_return",
     )
+
+
+# todo: enable on CI
+@pytest.mark.skip("Requires setting ICON_GRID_LOC environment variable.")
+def test_py2fgen_compilation_and_execution_diffusion(cli_runner, samples_path):
+    run_test_case(
+        cli_runner,
+        "icon4pytools.py2fgen.wrappers.diffusion",
+        "diffusion_init,diffusion_run",
+        "diffusion_plugin",
+        "CPU",
+        samples_path,
+        "test_diffusion",
+    )
+
+
+# todo: enable on CI
+@pytest.mark.skip("Requires setting various environment variables.")
+@pytest.mark.parametrize(
+    "function_name, plugin_name, test_name, backend, extra_flags",
+    [
+        ("square", "square_plugin", "test_square", "GPU", ("-acc", "-Minfo=acc")),
+        ("multi_return", "multi_return_plugin", "test_multi_return", "GPU", ("-acc", "-Minfo=acc")),
+    ],
+)
+def test_py2fgen_compilation_and_execution_gpu(
+    cli_runner,
+    function_name,
+    plugin_name,
+    test_name,
+    backend,
+    samples_path,
+    wrapper_module,
+    extra_flags,
+):
+    run_test_case(
+        cli_runner,
+        wrapper_module,
+        function_name,
+        plugin_name,
+        backend,
+        samples_path,
+        test_name,
+        "/opt/nvidia/hpc_sdk/Linux_x86_64/2024/compilers/bin/nvfortran",  # Ensure NVFORTRAN_COMPILER is set in your environment variables
+        extra_flags,
+    )
+
+
+# todo: enable on CI
+# Need to compile using nvfortran, and set CUDACXX path to nvcc cuda compiler. Also need to set ICON_GRID_LOC for path to gridfile, and ICON4PY_BACKEND to determine device at runtime.
+@pytest.mark.skip("Requires setting various environment variables.")
+@pytest.mark.parametrize(
+    "backend, extra_flags",
+    [("GPU", ("-acc", "-Minfo=acc"))],
+)
+def test_py2fgen_compilation_and_execution_diffusion_gpu(
+    cli_runner, samples_path, backend, extra_flags
+):
+    # todo: requires setting ICON_GRID_LOC
+    run_test_case(
+        cli_runner,
+        "icon4pytools.py2fgen.wrappers.diffusion",
+        "diffusion_init,diffusion_run",
+        "diffusion_plugin",
+        backend,
+        samples_path,
+        "test_diffusion",
+        "/opt/nvidia/hpc_sdk/Linux_x86_64/2024/compilers/bin/nvfortran",  # todo: set nvfortran location in base.yml file.
+        extra_flags,
+    )
diff --git a/tools/tests/py2fgen/test_codegen.py b/tools/tests/py2fgen/test_codegen.py
index 03bae08c30..1ff57f9410 100644
--- a/tools/tests/py2fgen/test_codegen.py
+++ b/tools/tests/py2fgen/test_codegen.py
@@ -16,10 +16,14 @@
 from gt4py.next.type_system.type_specifications import ScalarKind
 from icon4py.model.common.dimension import CellDim, KDim
 
+from icon4pytools.py2fgen.generate import (
+    generate_c_header,
+    generate_f90_interface,
+    generate_python_wrapper,
+)
 from icon4pytools.py2fgen.template import (
     CffiPlugin,
     CHeaderGenerator,
-    F90InterfaceGenerator,
     Func,
     FuncParameter,
     as_f90_value,
@@ -82,20 +86,20 @@ def test_as_target(param, expected):
 
 def test_cheader_generation_for_single_function():
     plugin = CffiPlugin(
-        module_name="libtest", plugin_name="libtest_plugin", function=foo, imports=["import foo"]
+        module_name="libtest", plugin_name="libtest_plugin", functions=[foo], imports=["import foo"]
     )
 
     header = CHeaderGenerator.apply(plugin)
-    assert header == "extern void foo_wrapper(int one, double* two, int n_Cell, int n_K);"
+    assert header == "extern int foo_wrapper(int one, double* two, int n_Cell, int n_K);"
 
 
 def test_cheader_for_pointer_args():
     plugin = CffiPlugin(
-        module_name="libtest", plugin_name="libtest_plugin", function=bar, imports=["import bar"]
+        module_name="libtest", plugin_name="libtest_plugin", functions=[bar], imports=["import bar"]
     )
 
     header = CHeaderGenerator.apply(plugin)
-    assert header == "extern void bar_wrapper(float* one, int two, int n_Cell, int n_K);"
+    assert header == "extern int bar_wrapper(float* one, int two, int n_Cell, int n_K);"
 
 
 def compare_ignore_whitespace(s1: str, s2: str):
@@ -103,29 +107,263 @@ def compare_ignore_whitespace(s1: str, s2: str):
     return s1.translate(no_whitespace) == s2.translate(no_whitespace)
 
 
-def test_fortran_interface():
-    plugin = CffiPlugin(
-        module_name="libtest", plugin_name="libtest_plugin", function=foo, imports=["import foo"]
+@pytest.fixture
+def dummy_plugin():
+    return CffiPlugin(
+        module_name="libtest",
+        plugin_name="libtest_plugin",
+        functions=[foo, bar],
+        imports=["import foo_module_x\nimport bar_module_y"],
     )
-    interface = F90InterfaceGenerator.apply(plugin)
+
+
+def test_fortran_interface(dummy_plugin):
+    interface = generate_f90_interface(dummy_plugin)
     expected = """
     module libtest_plugin
-    use, intrinsic:: iso_c_binding
-    implicit none
+   use, intrinsic :: iso_c_binding
+   implicit none
+
+   public :: foo
+
+   public :: bar
+
+   interface
+
+      function foo_wrapper(one, &
+                           two, &
+                           n_Cell, &
+                           n_K) bind(c, name="foo_wrapper") result(rc)
+         import :: c_int, c_double, c_bool, c_ptr
+
+         integer(c_int), value :: n_Cell
+
+         integer(c_int), value :: n_K
+
+         integer(c_int) :: rc  ! Stores the return code
+
+         integer(c_int), value, target :: one
+
+         real(c_double), dimension(*), target :: two
+
+      end function foo_wrapper
+
+      function bar_wrapper(one, &
+                           two, &
+                           n_Cell, &
+                           n_K) bind(c, name="bar_wrapper") result(rc)
+         import :: c_int, c_double, c_bool, c_ptr
+
+         integer(c_int), value :: n_Cell
+
+         integer(c_int), value :: n_K
+
+         integer(c_int) :: rc  ! Stores the return code
+
+         real(c_float), dimension(*), target :: one
+
+         integer(c_int), value, target :: two
+
+      end function bar_wrapper
 
-    public
-    interface
-        subroutine foo_wrapper(one, &
+   end interface
+
+contains
+
+   subroutine foo(one, &
+                  two, &
+                  rc)
+      use, intrinsic :: iso_c_binding
+
+      integer(c_int) :: n_Cell
+
+      integer(c_int) :: n_K
+
+      integer(c_int), value, target :: one
+
+      real(c_double), dimension(:, :), target :: two
+
+      integer(c_int) :: rc  ! Stores the return code
+
+      !$ACC host_data use_device( &
+      !$ACC two &
+      !$ACC )
+
+      n_Cell = SIZE(two, 1)
+
+      n_K = SIZE(two, 2)
+
+      rc = foo_wrapper(one, &
                        two, &
                        n_Cell, &
-                       n_K) bind(c, name='foo_wrapper')
-            use, intrinsic :: iso_c_binding
-            integer(c_int), value, target :: n_Cell
-            integer(c_int), value, target :: n_K
-            integer(c_int), value, target :: one
-            real(c_double), dimension(:, :), target :: two(n_Cell, n_K)
-        end subroutine foo_wrapper
-    end interface
-    end module
+                       n_K)
+
+      !$acc end host_data
+   end subroutine foo
+
+   subroutine bar(one, &
+                  two, &
+                  rc)
+      use, intrinsic :: iso_c_binding
+
+      integer(c_int) :: n_Cell
+
+      integer(c_int) :: n_K
+
+      real(c_float), dimension(:, :), target :: one
+
+      integer(c_int), value, target :: two
+
+      integer(c_int) :: rc  ! Stores the return code
+
+      !$ACC host_data use_device( &
+      !$ACC one &
+      !$ACC )
+
+      n_Cell = SIZE(one, 1)
+
+      n_K = SIZE(one, 2)
+
+      rc = bar_wrapper(one, &
+                       two, &
+                       n_Cell, &
+                       n_K)
+
+      !$acc end host_data
+   end subroutine bar
+
+end module
+"""
+    assert compare_ignore_whitespace(interface, expected)
+
+
+def test_python_wrapper(dummy_plugin):
+    interface = generate_python_wrapper(dummy_plugin, "GPU", False)
+    expected = '''
+# imports for generated wrapper code
+import logging
+from libtest_plugin import ffi
+import numpy as np
+import cupy as cp
+from numpy.typing import NDArray
+from gt4py.next.iterator.embedded import np_as_located_field
+
+# logger setup
+log_format = '%(asctime)s.%(msecs)03d - %(levelname)s - %(message)s'
+logging.basicConfig(level=logging.ERROR,
+                    format=log_format,
+                    datefmt='%Y-%m-%d %H:%M:%S')
+logging.info(cp.show_config())
+
+# embedded module imports
+import foo_module_x
+import bar_module_y
+
+# embedded function imports
+from libtest import foo
+from libtest import bar
+
+def unpack_gpu(ptr, *sizes: int):
+    """
+    Converts a C pointer into a CuPy array to directly manipulate memory allocated in Fortran.
+    This function is needed for operations that require in-place modification of GPU data,
+    enabling changes made in Python to reflect immediately in the original Fortran memory space.
+
+    Args:
+        ptr (cffi.CData): A CFFI pointer to GPU memory allocated by OpenACC, representing
+                          the starting address of the data. This pointer must correspond to
+                          a contiguous block of memory whose total size matches the product
+                          of the specified dimensions.
+        *sizes (int): Variable length argument list specifying the dimensions of the array.
+                      These sizes determine the shape of the resulting CuPy array.
+
+    Returns:
+        cp.ndarray: A CuPy array that provides a direct view of the data pointed to by `ptr`.
+                    This array shares the underlying data with the original Fortran code, allowing
+                    modifications made through the array to affect the original data.
+    """
+
+    if not sizes:
+        raise ValueError("Sizes must be provided to determine the array shape.")
+
+    length = np.prod(sizes)
+    c_type = ffi.getctype(ffi.typeof(ptr).item)
+
+    dtype_map = {
+        "int": cp.int32,
+        "double": cp.float64,
+    }
+    dtype = dtype_map.get(c_type, None)
+    if dtype is None:
+        raise ValueError(f"Unsupported C data type: {c_type}")
+
+    itemsize = ffi.sizeof(c_type)
+    total_size = length * itemsize
+
+    # cupy array from OpenACC device pointer
+    current_device = cp.cuda.Device()
+    ptr_val = int(ffi.cast("uintptr_t", ptr))
+    mem = cp.cuda.UnownedMemory(ptr_val, total_size, owner=ptr, device_id=current_device.id)
+    memptr = cp.cuda.MemoryPointer(mem, 0)
+    arr = cp.ndarray(shape=sizes, dtype=dtype, memptr=memptr, order="F")
+
+    return arr
+
+def int_array_to_bool_array(int_array: NDArray) -> NDArray:
+    """
+    Converts a NumPy array of integers to a boolean array.
+    In the input array, 0 represents False, and any non-zero value (1 or -1) represents True.
+
+    Args:
+        int_array: A NumPy array of integers.
+
+    Returns:
+        A NumPy array of booleans.
+    """
+    bool_array = int_array != 0
+    return bool_array
+
+@ffi.def_extern()
+def foo_wrapper(one: int32, two: Field[CellDim, KDim], float64], n_Cell: int32, n_K: int32):
+    try:
+        # Unpack pointers into Ndarrays
+        two = unpack_gpu(two, n_Cell, n_K)
+
+        # Allocate GT4Py Fields
+        two = np_as_located_field(CellDim, KDim)(two)
+
+        foo(one, two)
+
+    except Exception as e:
+        logging.exception(f"A Python error occurred: {e}")
+        return 1
+
+    return 0
+
+@ffi.def_extern()
+def bar_wrapper(one: Field[CellDim, KDim], float64], two: int32, n_Cell: int32, n_K: int32):
+    try:
+        # Unpack pointers into Ndarrays
+        one = unpack_gpu(one, n_Cell, n_K)
+
+        # Allocate GT4Py Fields
+        one = np_as_located_field(CellDim, KDim)(one)
+
+        bar(one, two)
+
+    except Exception as e:
+        logging.exception(f"A Python error occurred: {e}")
+        return 1
+
+    return 0
+    '''
+    assert compare_ignore_whitespace(interface, expected)
+
+
+def test_c_header(dummy_plugin):
+    interface = generate_c_header(dummy_plugin)
+    expected = """
+    extern int foo_wrapper(int one, double *two, int n_Cell, int n_K);
+    extern int bar_wrapper(float *one, int two, int n_Cell, int n_K);
     """
     assert compare_ignore_whitespace(interface, expected)
diff --git a/tools/tests/py2fgen/test_diffusion_wrapper.py b/tools/tests/py2fgen/test_diffusion_wrapper.py
new file mode 100644
index 0000000000..cb49e594c1
--- /dev/null
+++ b/tools/tests/py2fgen/test_diffusion_wrapper.py
@@ -0,0 +1,219 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# This file is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or any later
+# version. See the LICENSE.txt file at the top-level directory of this
+# distribution for a copy of the license or check <https://www.gnu.org/licenses/>.
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+# type: ignore
+import pytest
+from gt4py.next import as_field
+from icon4py.model.atmosphere.diffusion.diffusion import DiffusionType
+from icon4py.model.common.dimension import (
+    C2E2CDim,
+    C2E2CODim,
+    C2EDim,
+    CellDim,
+    E2C2VDim,
+    E2CDim,
+    EdgeDim,
+    KDim,
+    V2EDim,
+    VertexDim,
+)
+from icon4py.model.common.settings import xp
+
+from icon4pytools.py2fgen.wrappers.diffusion import diffusion_init, diffusion_run
+
+
+# todo(samkellerhals): turn on and off using a marker/option (required ICON_GRID_LOC)
+@pytest.mark.skip
+def test_diffusion_wrapper_py():
+    # grid parameters
+    num_cells = 20480
+    num_edges = 30720
+    num_vertices = 10242
+    num_levels = 60
+    num_c2ec2o = 4
+    num_v2e = 6
+    num_c2e = 2
+    num_e2c2v = 4
+    num_c2e2c = 3
+    num_e2c = 3
+    mean_cell_area = 24907282236.708576
+
+    # other configuration parameters
+    ndyn_substeps = 2
+    dtime = 2.0
+    rayleigh_damping_height = 50000
+    nflatlev = 30
+    nflat_gradp = 59
+
+    # diffusion configuration
+    diffusion_type = DiffusionType.SMAGORINSKY_4TH_ORDER  # 5
+    hdiff_w = True
+    hdiff_vn = True
+    zdiffu_t = True  # Setting this to true triggers running stencil 15 with the boolean masks.
+    type_t_diffu = 2
+    type_vn_diffu = 1
+    hdiff_efdt_ratio = 24.0
+    smagorinski_scaling_factor = 0.025
+    hdiff_temp = True
+
+    # input data - numpy
+    rng = xp.random.default_rng()
+
+    vct_a = rng.uniform(
+        low=0, high=75000, size=(num_levels,)
+    )  # has to be from 0 to 75000, must have larger values than rayleigh damping height
+    theta_ref_mc = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    wgtfac_c = rng.uniform(low=0, high=1, size=(num_cells, num_levels + 1))
+    e_bln_c_s = rng.uniform(low=0, high=1, size=(num_cells, num_c2e))
+    geofac_div = rng.uniform(low=0, high=1, size=(num_cells, num_c2e))
+    geofac_grg_x = rng.uniform(low=0, high=1, size=(num_cells, num_c2ec2o))
+    geofac_grg_y = rng.uniform(low=0, high=1, size=(num_cells, num_c2ec2o))
+    geofac_n2s = rng.uniform(low=0, high=1, size=(num_cells, num_c2ec2o))
+    nudgecoeff_e = xp.zeros((num_edges,))
+    rbf_coeff_1 = rng.uniform(low=0, high=1, size=(num_vertices, num_v2e))
+    rbf_coeff_2 = rng.uniform(low=0, high=1, size=(num_vertices, num_v2e))
+    dwdx = xp.zeros((num_cells, num_levels))
+    dwdy = xp.zeros((num_cells, num_levels))
+    hdef_ic = xp.zeros((num_cells, num_levels + 1))
+    div_ic = xp.zeros((num_cells, num_levels + 1))
+    w = rng.uniform(low=0, high=1, size=(num_cells, num_levels + 1))
+    vn = rng.uniform(low=0, high=1, size=(num_edges, num_levels))
+    exner = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    theta_v = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    rho = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    dual_normal_cell_x = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    dual_normal_cell_y = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    dual_normal_vert_x = rng.uniform(low=0, high=1, size=(num_edges, num_e2c2v))
+    dual_normal_vert_y = rng.uniform(low=0, high=1, size=(num_edges, num_e2c2v))
+    primal_normal_cell_x = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    primal_normal_cell_y = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    primal_normal_vert_x = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    primal_normal_vert_y = rng.uniform(low=0, high=1, size=(num_edges, num_e2c))
+    tangent_orientation = rng.uniform(low=0, high=1, size=(num_edges))
+    inverse_primal_edge_lengths = rng.uniform(low=0, high=1, size=(num_edges))
+    inv_dual_edge_length = rng.uniform(low=0, high=1, size=(num_edges))
+    inv_vert_vert_length = rng.uniform(low=0, high=1, size=(num_edges))
+    edge_areas = rng.uniform(low=0, high=1, size=(num_edges))
+    f_e = rng.uniform(low=0, high=1, size=(num_edges))
+    cell_areas = rng.uniform(low=0, high=1, size=(num_cells))
+    zd_diffcoef = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    zd_vertoffset = xp.round(
+        rng.uniform(low=0, high=1, size=(num_cells, num_c2e2c, num_levels))
+    ).astype(xp.int32)
+    zd_intcoef = rng.uniform(low=0, high=1, size=(num_cells, num_c2e2c, num_levels))
+
+    # Create a boolean array based on a condition
+    mask_hdiff = rng.uniform(low=0, high=1, size=(num_cells, num_levels))
+    mask_hdiff = mask_hdiff < 0.5
+
+    # input data - gt4py fields
+    theta_ref_mc = as_field((CellDim, KDim), theta_ref_mc)
+    wgtfac_c = as_field((CellDim, KDim), wgtfac_c)
+    vct_a = as_field((KDim,), vct_a)
+    e_bln_c_s = as_field((CellDim, C2EDim), e_bln_c_s)
+    geofac_div = as_field((CellDim, C2EDim), geofac_div)
+    geofac_grg_x = as_field((CellDim, C2E2CODim), geofac_grg_x)
+    geofac_grg_y = as_field((CellDim, C2E2CODim), geofac_grg_y)
+    geofac_n2s = as_field((CellDim, C2E2CODim), geofac_n2s)
+    nudgecoeff_e = as_field((EdgeDim,), nudgecoeff_e)
+    rbf_coeff_1 = as_field((VertexDim, V2EDim), rbf_coeff_1)
+    rbf_coeff_2 = as_field((VertexDim, V2EDim), rbf_coeff_2)
+    dwdx = as_field((CellDim, KDim), dwdx)
+    dwdy = as_field((CellDim, KDim), dwdy)
+    hdef_ic = as_field((CellDim, KDim), hdef_ic)
+    div_ic = as_field((CellDim, KDim), div_ic)
+    mask_hdiff = as_field((CellDim, KDim), mask_hdiff)
+    zd_diffcoef = as_field((CellDim, KDim), zd_diffcoef)
+    zd_vertoffset = as_field((CellDim, C2E2CDim, KDim), zd_vertoffset)
+    zd_intcoef = as_field((CellDim, C2E2CDim, KDim), zd_intcoef)
+    w = as_field((CellDim, KDim), w)
+    vn = as_field((EdgeDim, KDim), vn)
+    exner = as_field((CellDim, KDim), exner)
+    theta_v = as_field((CellDim, KDim), theta_v)
+    rho = as_field((CellDim, KDim), rho)
+    dual_normal_cell_x = as_field((EdgeDim, E2CDim), dual_normal_cell_x)
+    dual_normal_cell_y = as_field((EdgeDim, E2CDim), dual_normal_cell_y)
+    dual_normal_vert_x = as_field((EdgeDim, E2C2VDim), dual_normal_vert_x)
+    dual_normal_vert_y = as_field((EdgeDim, E2C2VDim), dual_normal_vert_y)
+    primal_normal_cell_x = as_field((EdgeDim, E2CDim), primal_normal_cell_x)
+    primal_normal_cell_y = as_field((EdgeDim, E2CDim), primal_normal_cell_y)
+    primal_normal_vert_x = as_field((EdgeDim, E2C2VDim), primal_normal_vert_x)
+    primal_normal_vert_y = as_field((EdgeDim, E2C2VDim), primal_normal_vert_y)
+    tangent_orientation = as_field((EdgeDim,), tangent_orientation)
+    inverse_primal_edge_lengths = as_field((EdgeDim,), inverse_primal_edge_lengths)
+    inv_dual_edge_length = as_field((EdgeDim,), inv_dual_edge_length)
+    inv_vert_vert_length = as_field((EdgeDim,), inv_vert_vert_length)
+    edge_areas = as_field((EdgeDim,), edge_areas)
+    f_e = as_field((EdgeDim,), f_e)
+    cell_areas = as_field((CellDim,), cell_areas)
+
+    diffusion_init(
+        vct_a=vct_a,
+        theta_ref_mc=theta_ref_mc,
+        wgtfac_c=wgtfac_c,
+        e_bln_c_s=e_bln_c_s,
+        geofac_div=geofac_div,
+        geofac_grg_x=geofac_grg_x,
+        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,
+        num_levels=num_levels,
+        mean_cell_area=mean_cell_area,
+        ndyn_substeps=ndyn_substeps,
+        rayleigh_damping_height=rayleigh_damping_height,
+        nflatlev=nflatlev,
+        nflat_gradp=nflat_gradp,
+        diffusion_type=diffusion_type,
+        hdiff_w=hdiff_w,
+        hdiff_vn=hdiff_vn,
+        zdiffu_t=zdiffu_t,
+        type_t_diffu=type_t_diffu,
+        type_vn_diffu=type_vn_diffu,
+        hdiff_efdt_ratio=hdiff_efdt_ratio,
+        smagorinski_scaling_factor=smagorinski_scaling_factor,
+        hdiff_temp=hdiff_temp,
+        tangent_orientation=tangent_orientation,
+        inverse_primal_edge_lengths=inverse_primal_edge_lengths,
+        inv_dual_edge_length=inv_dual_edge_length,
+        inv_vert_vert_length=inv_vert_vert_length,
+        edge_areas=edge_areas,
+        f_e=f_e,
+        cell_areas=cell_areas,
+        primal_normal_vert_x=primal_normal_vert_x,
+        primal_normal_vert_y=primal_normal_vert_y,
+        dual_normal_vert_x=dual_normal_vert_x,
+        dual_normal_vert_y=dual_normal_vert_y,
+        primal_normal_cell_x=primal_normal_cell_x,
+        primal_normal_cell_y=primal_normal_cell_y,
+        dual_normal_cell_x=dual_normal_cell_x,
+        dual_normal_cell_y=dual_normal_cell_y,
+        mask_hdiff=mask_hdiff,
+        zd_diffcoef=zd_diffcoef,
+        zd_vertoffset=zd_vertoffset,
+        zd_intcoef=zd_intcoef,
+    )
+
+    diffusion_run(
+        w=w,
+        vn=vn,
+        exner=exner,
+        theta_v=theta_v,
+        rho=rho,
+        hdef_ic=hdef_ic,
+        div_ic=div_ic,
+        dwdx=dwdx,
+        dwdy=dwdy,
+        dtime=dtime,
+        linit=False,
+    )
diff --git a/tools/tests/py2fgen/test_parsing.py b/tools/tests/py2fgen/test_parsing.py
index 3124d0bac0..71c59fa1f5 100644
--- a/tools/tests/py2fgen/test_parsing.py
+++ b/tools/tests/py2fgen/test_parsing.py
@@ -11,13 +11,48 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-from icon4pytools.py2fgen.parsing import parse
+import ast
+
+import pytest
+
+from icon4pytools.py2fgen.parsing import ImportStmtVisitor, TypeHintVisitor, parse
 from icon4pytools.py2fgen.template import CffiPlugin
 
 
+source = """
+import foo
+import bar
+
+def test_function(x: Field[[EdgeDim, KDim], float64], y: int):
+    return x * y
+"""
+
+
 def test_parse_functions_on_wrapper():
     module_path = "icon4pytools.py2fgen.wrappers.diffusion"
-    function_name = "diffusion_init"
-    plugin = parse(module_path, function_name)
+    functions = ["diffusion_init", "diffusion_run"]
+    plugin = parse(module_path, functions, "diffusion_plugin")
     assert isinstance(plugin, CffiPlugin)
-    assert plugin.plugin_name == "diffusion_init_plugin"
+
+
+def test_import_visitor():
+    tree = ast.parse(source)
+    extractor = ImportStmtVisitor()
+    extractor.visit(tree)
+    expected_imports = ["import foo", "import bar"]
+    assert extractor.import_statements == expected_imports
+
+
+def test_type_hint_visitor():
+    tree = ast.parse(source)
+    visitor = TypeHintVisitor()
+    visitor.visit(tree)
+    expected_type_hints = {"x": "Field[[EdgeDim, KDim], float64]", "y": "int"}
+    assert visitor.type_hints == expected_type_hints
+
+
+def test_function_missing_type_hints():
+    tree = ast.parse(source.replace(": int", ""))
+    visitor = TypeHintVisitor()
+    with pytest.raises(TypeError):
+        visitor.visit(tree)
diff --git a/tools/tox.ini b/tools/tox.ini
index be68d08e88..0c8b7eb934 100644
--- a/tools/tox.ini
+++ b/tools/tox.ini
@@ -15,7 +15,7 @@ deps =
     -r {toxinidir}/requirements-dev.txt
 commands =
     -pytest -v -s -n auto -cache-clear --cov --cov-config=pyproject.toml --cov-reset --doctest-modules src/icon4pytools/
-    pytest -vvv -l -s --full-trace -n auto --cov-config=pyproject.toml
+    pytest -vvv -s -n auto --cov-config=pyproject.toml
 commands_post =
     rm -rf tests/_reports/coverage_html
     -coverage html