diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_nonhydro.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_nonhydro.py
index 8b939ed6d6..5871ca0a82 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_nonhydro.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_nonhydro.py
@@ -726,7 +726,6 @@ def __init__(
                 "start_cell_halo_level_2": self._start_cell_halo_level_2,
                 "end_cell_halo": self._end_cell_halo,
                 "end_cell_halo_level_2": self._end_cell_halo_level_2,
-                "start_cell_lateral_boundary": self._start_cell_lateral_boundary,
             },
             vertical_sizes={
                 "nflatlev": self._vertical_params.nflatlev,
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/compute_cell_diagnostics_for_dycore.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/compute_cell_diagnostics_for_dycore.py
index e976a24809..28dc64a64a 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/compute_cell_diagnostics_for_dycore.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/compute_cell_diagnostics_for_dycore.py
@@ -6,33 +6,22 @@
 # Please, refer to the LICENSE file in the root directory.
 # SPDX-License-Identifier: BSD-3-Clause
 
-# 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
-
 from typing import Final
 
 import gt4py.next as gtx
-from gt4py.next import astype, broadcast, maximum
+from gt4py.next import astype, maximum
 from gt4py.next.experimental import concat_where
 
 from icon4py.model.atmosphere.dycore.dycore_states import HorizontalPressureDiscretizationType
 from icon4py.model.atmosphere.dycore.stencils.compute_perturbation_of_rho_and_theta import (
     _compute_perturbation_of_rho_and_theta,
 )
+from icon4py.model.atmosphere.dycore.stencils.extrapolate_quadratically_to_surface import (
+    _extrapolate_quadratically_to_surface,
+)
 from icon4py.model.atmosphere.dycore.stencils.extrapolate_temporally_exner_pressure import (
     _extrapolate_temporally_exner_pressure,
 )
-from icon4py.model.atmosphere.dycore.stencils.interpolate_to_surface import _interpolate_to_surface
 from icon4py.model.common import dimension as dims, field_type_aliases as fa, type_alias as ta
 from icon4py.model.common.dimension import Koff
 from icon4py.model.common.interpolation.stencils.interpolate_cell_field_to_half_levels_vp import (
@@ -68,6 +57,9 @@ def _calculate_pressure_buoyancy_acceleration_at_cells_on_half_levels(
 
 @gtx.field_operator
 def _compute_perturbed_quantities_and_interpolation(
+    time_extrapolation_parameter_for_exner: fa.CellKField[ta.vpfloat],
+    current_exner: fa.CellKField[ta.wpfloat],
+    reference_exner_at_cells_on_model_levels: fa.CellKField[ta.vpfloat],
     current_rho: fa.CellKField[ta.wpfloat],
     reference_rho_at_cells_on_model_levels: fa.CellKField[ta.wpfloat],
     current_theta_v: fa.CellKField[ta.wpfloat],
@@ -77,30 +69,50 @@ def _compute_perturbed_quantities_and_interpolation(
     perturbed_exner_at_cells_on_model_levels: fa.CellKField[ta.wpfloat],
     ddz_of_reference_exner_at_cells_on_half_levels: fa.CellKField[ta.wpfloat],
     ddqz_z_half: fa.CellKField[ta.vpfloat],
-    pressure_buoyancy_acceleration_at_cells_on_half_levels: fa.CellKField[ta.vpfloat],
-    rho_at_cells_on_half_levels: fa.CellKField[ta.wpfloat],
     exner_at_cells_on_half_levels: fa.CellKField[ta.vpfloat],
-    temporal_extrapolation_of_perturbed_exner: fa.CellKField[ta.vpfloat],
-    theta_v_at_cells_on_half_levels: fa.CellKField[ta.wpfloat],
+    wgtfacq_c: fa.CellKField[ta.vpfloat],
+    reference_theta_at_cells_on_half_levels: fa.CellKField[ta.vpfloat],
+    inv_ddqz_z_full: fa.CellKField[ta.vpfloat],
+    d2dexdz2_fac1_mc: fa.CellKField[ta.vpfloat],
+    d2dexdz2_fac2_mc: fa.CellKField[ta.vpfloat],
+    ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels: fa.CellKField[ta.vpfloat],
+    d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels: fa.CellKField[ta.vpfloat],
     igradp_method: gtx.int32,
-    nflatlev: gtx.int32,
+    surface_level: gtx.int32,
+    start_cell_halo_level_2: gtx.int32,
+    end_cell_halo_level_2: gtx.int32,
 ) -> tuple[
-    fa.CellKField[ta.vpfloat],
+    fa.CellKField[ta.wpfloat],
+    fa.CellKField[ta.wpfloat],
     fa.CellKField[ta.vpfloat],
     fa.CellKField[ta.vpfloat],
     fa.CellKField[ta.wpfloat],
     fa.CellKField[ta.wpfloat],
+    fa.CellKField[ta.wpfloat],
+    fa.CellKField[ta.wpfloat],
+    fa.CellKField[ta.wpfloat],
     fa.CellKField[ta.vpfloat],
     fa.CellKField[ta.vpfloat],
-    fa.CellKField[ta.wpfloat],
 ]:
+    (
+        temporal_extrapolation_of_perturbed_exner,
+        perturbed_exner_at_cells_on_model_levels,
+    ) = _extrapolate_temporally_exner_pressure(
+        exner_exfac=time_extrapolation_parameter_for_exner,
+        exner=current_exner,
+        exner_ref_mc=reference_exner_at_cells_on_model_levels,
+        exner_pr=perturbed_exner_at_cells_on_model_levels,
+    )
+
     exner_at_cells_on_half_levels = (
         concat_where(
-            (maximum(1, nflatlev) <= dims.KDim),
+            dims.KDim < surface_level - 1,
             _interpolate_cell_field_to_half_levels_vp(
                 wgtfac_c=wgtfac_c, interpolant=temporal_extrapolation_of_perturbed_exner
             ),
-            exner_at_cells_on_half_levels,
+            _extrapolate_quadratically_to_surface(
+                wgtfacq_c=wgtfacq_c, interpolant=temporal_extrapolation_of_perturbed_exner
+            ),
         )
         if igradp_method == horzpres_discr_type.TAYLOR_HYDRO
         else exner_at_cells_on_half_levels
@@ -116,34 +128,31 @@ def _compute_perturbed_quantities_and_interpolation(
         reference_theta_at_cells_on_model_levels,
     )
 
-    rho_at_cells_on_half_levels = concat_where(
-        dims.KDim >= 1,
-        _interpolate_cell_field_to_half_levels_wp(wgtfac_c, current_rho),
-        rho_at_cells_on_half_levels,
-    )
+    rho_at_cells_on_half_levels = _interpolate_cell_field_to_half_levels_wp(wgtfac_c, current_rho)
 
     wgtfac_c_wp = astype(wgtfac_c, wpfloat)
 
     perturbed_theta_v_at_cells_on_half_levels = concat_where(
-        dims.KDim >= 1,
+        dims.KDim < surface_level - 1,
         _interpolate_cell_field_to_half_levels_vp(
             wgtfac_c=wgtfac_c, interpolant=perturbed_theta_v_at_cells_on_model_levels
         ),
-        broadcast(0.0, (dims.CellDim, dims.KDim)),
+        _extrapolate_quadratically_to_surface(
+            wgtfacq_c=wgtfacq_c, interpolant=perturbed_theta_v_at_cells_on_model_levels
+        ),
     )
 
     theta_v_at_cells_on_half_levels = concat_where(
-        dims.KDim >= 1,
+        dims.KDim < surface_level - 1,
         _interpolate_cell_field_to_half_levels_wp(
             wgtfac_c=wgtfac_c_wp, interpolant=current_theta_v
         ),
-        theta_v_at_cells_on_half_levels,
+        reference_theta_at_cells_on_half_levels + perturbed_theta_v_at_cells_on_half_levels,
     )
 
     ddqz_z_half_wp = astype(ddqz_z_half, wpfloat)
 
-    pressure_buoyancy_acceleration_at_cells_on_half_levels = concat_where(
-        dims.KDim >= 1,
+    pressure_buoyancy_acceleration_at_cells_on_half_levels = (
         _calculate_pressure_buoyancy_acceleration_at_cells_on_half_levels(
             exner_w_explicit_weight_parameter,
             theta_v_at_cells_on_half_levels,
@@ -151,115 +160,62 @@ def _compute_perturbed_quantities_and_interpolation(
             ddqz_z_half_wp,
             perturbed_theta_v_at_cells_on_half_levels,
             ddz_of_reference_exner_at_cells_on_half_levels,
-        ),
-        pressure_buoyancy_acceleration_at_cells_on_half_levels,
-    )
-
-    return (
-        perturbed_rho_at_cells_on_model_levels,
-        perturbed_theta_v_at_cells_on_model_levels,
-        perturbed_exner_at_cells_on_model_levels,
-        rho_at_cells_on_half_levels,
-        exner_at_cells_on_half_levels,
-        perturbed_theta_v_at_cells_on_half_levels,
-        theta_v_at_cells_on_half_levels,
-        pressure_buoyancy_acceleration_at_cells_on_half_levels,
-    )
-
-
-@gtx.field_operator
-def _surface_computations(
-    wgtfacq_c: fa.CellKField[ta.wpfloat],
-    exner_at_cells_on_half_levels: fa.CellKField[ta.vpfloat],
-    temporal_extrapolation_of_perturbed_exner: fa.CellKField[vpfloat],
-    igradp_method: gtx.int32,
-) -> fa.CellKField[ta.vpfloat]:
-    exner_at_cells_on_half_levels = (
-        _interpolate_to_surface(
-            wgtfacq_c=wgtfacq_c, interpolant=temporal_extrapolation_of_perturbed_exner
         )
-        if igradp_method == horzpres_discr_type.TAYLOR_HYDRO
-        else exner_at_cells_on_half_levels
     )
 
-    return exner_at_cells_on_half_levels
-
-
-@gtx.field_operator
-def _compute_first_and_second_vertical_derivative_of_exner(
-    exner_at_cells_on_half_levels: fa.CellKField[vpfloat],
-    inv_ddqz_z_full: fa.CellKField[vpfloat],
-    ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels: fa.CellKField[vpfloat],
-    d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels: fa.CellKField[vpfloat],
-    perturbed_theta_v_at_cells_on_half_levels: fa.CellKField[vpfloat],
-    d2dexdz2_fac1_mc: fa.CellKField[vpfloat],
-    d2dexdz2_fac2_mc: fa.CellKField[vpfloat],
-    perturbed_theta_v_at_cells_on_model_levels: fa.CellKField[vpfloat],
-    igradp_method: gtx.int32,
-    nflatlev: gtx.int32,
-    nflat_gradp: gtx.int32,
-) -> tuple[
-    fa.CellKField[vpfloat],
-    fa.CellKField[vpfloat],
-]:
     ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels = (
-        concat_where(
-            (nflatlev <= dims.KDim),
+        (
             _compute_first_vertical_derivative_at_cells(
                 exner_at_cells_on_half_levels, inv_ddqz_z_full
-            ),
-            ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
+            )
         )
         if igradp_method == horzpres_discr_type.TAYLOR_HYDRO
         else ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels
     )
 
     d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels = (
-        concat_where(
-            (nflat_gradp <= dims.KDim),
-            -vpfloat("0.5")
-            * (
-                (
-                    perturbed_theta_v_at_cells_on_half_levels
-                    - perturbed_theta_v_at_cells_on_half_levels(Koff[1])
-                )
-                * d2dexdz2_fac1_mc
-                + perturbed_theta_v_at_cells_on_model_levels * d2dexdz2_fac2_mc
-            ),
-            d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
+        -vpfloat("0.5")
+        * (
+            (
+                perturbed_theta_v_at_cells_on_half_levels
+                - perturbed_theta_v_at_cells_on_half_levels(Koff[1])
+            )
+            * d2dexdz2_fac1_mc
+            + perturbed_theta_v_at_cells_on_model_levels * d2dexdz2_fac2_mc
         )
         if igradp_method == horzpres_discr_type.TAYLOR_HYDRO
         else d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels
     )
 
-    return (
-        ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
-        d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
-    )
-
-
-@gtx.field_operator
-def _set_theta_v_and_exner_on_surface_level(
-    temporal_extrapolation_of_perturbed_exner: fa.CellKField[vpfloat],
-    wgtfacq_c: fa.CellKField[vpfloat],
-    perturbed_theta_v_at_cells_on_model_levels: fa.CellKField[vpfloat],
-    reference_theta_at_cells_on_half_levels: fa.CellKField[vpfloat],
-) -> tuple[fa.CellKField[vpfloat], fa.CellKField[wpfloat], fa.CellKField[vpfloat]]:
-    perturbed_theta_v_at_cells_on_half_levels = _interpolate_to_surface(
-        wgtfacq_c=wgtfacq_c, interpolant=perturbed_theta_v_at_cells_on_model_levels
-    )
-    theta_v_at_cells_on_half_levels = (
-        reference_theta_at_cells_on_half_levels + perturbed_theta_v_at_cells_on_half_levels
-    )
-
-    exner_at_cells_on_half_levels = _interpolate_to_surface(
-        wgtfacq_c=wgtfacq_c, interpolant=temporal_extrapolation_of_perturbed_exner
+    (
+        perturbed_rho_at_cells_on_model_levels,
+        perturbed_theta_v_at_cells_on_model_levels,
+    ) = concat_where(
+        (dims.CellDim >= start_cell_halo_level_2) & (dims.CellDim < end_cell_halo_level_2),
+        _compute_perturbation_of_rho_and_theta(
+            rho=current_rho,
+            rho_ref_mc=reference_rho_at_cells_on_model_levels,
+            theta_v=current_theta_v,
+            theta_ref_mc=reference_theta_at_cells_on_model_levels,
+        ),
+        (
+            perturbed_rho_at_cells_on_model_levels,
+            perturbed_theta_v_at_cells_on_model_levels,
+        ),
     )
 
     return (
-        perturbed_theta_v_at_cells_on_half_levels,
-        astype(theta_v_at_cells_on_half_levels, wpfloat),
+        perturbed_rho_at_cells_on_model_levels,
+        perturbed_theta_v_at_cells_on_model_levels,
+        perturbed_exner_at_cells_on_model_levels,
+        rho_at_cells_on_half_levels,
         exner_at_cells_on_half_levels,
+        perturbed_theta_v_at_cells_on_half_levels,
+        theta_v_at_cells_on_half_levels,
+        pressure_buoyancy_acceleration_at_cells_on_half_levels,
+        temporal_extrapolation_of_perturbed_exner,
+        ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
+        d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
     )
 
 
@@ -295,7 +251,6 @@ def compute_perturbed_quantities_and_interpolation(
     igradp_method: gtx.int32,
     nflatlev: gtx.int32,
     nflat_gradp: gtx.int32,
-    start_cell_lateral_boundary: gtx.int32,
     start_cell_lateral_boundary_level_3: gtx.int32,
     start_cell_halo_level_2: gtx.int32,
     end_cell_halo: gtx.int32,
@@ -362,31 +317,10 @@ def compute_perturbed_quantities_and_interpolation(
         - d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels
     """
 
-    _extrapolate_temporally_exner_pressure(
-        exner_exfac=time_extrapolation_parameter_for_exner,
-        exner=current_exner,
-        exner_ref_mc=reference_exner_at_cells_on_model_levels,
-        exner_pr=perturbed_exner_at_cells_on_model_levels,
-        out=(temporal_extrapolation_of_perturbed_exner, perturbed_exner_at_cells_on_model_levels),
-        domain={
-            dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
-            dims.KDim: (model_top, surface_level - 1),
-        },
-    )
-
-    _surface_computations(
-        wgtfacq_c=wgtfacq_c,
-        exner_at_cells_on_half_levels=exner_at_cells_on_half_levels,
-        temporal_extrapolation_of_perturbed_exner=temporal_extrapolation_of_perturbed_exner,
-        igradp_method=igradp_method,
-        out=exner_at_cells_on_half_levels,
-        domain={
-            dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
-            dims.KDim: (surface_level - 1, surface_level),
-        },
-    )
-
     _compute_perturbed_quantities_and_interpolation(
+        time_extrapolation_parameter_for_exner=time_extrapolation_parameter_for_exner,
+        current_exner=current_exner,
+        reference_exner_at_cells_on_model_levels=reference_exner_at_cells_on_model_levels,
         current_rho=current_rho,
         reference_rho_at_cells_on_model_levels=reference_rho_at_cells_on_model_levels,
         current_theta_v=current_theta_v,
@@ -396,13 +330,18 @@ def compute_perturbed_quantities_and_interpolation(
         perturbed_exner_at_cells_on_model_levels=perturbed_exner_at_cells_on_model_levels,
         ddz_of_reference_exner_at_cells_on_half_levels=ddz_of_reference_exner_at_cells_on_half_levels,
         ddqz_z_half=ddqz_z_half,
-        pressure_buoyancy_acceleration_at_cells_on_half_levels=pressure_buoyancy_acceleration_at_cells_on_half_levels,
-        rho_at_cells_on_half_levels=rho_at_cells_on_half_levels,
         exner_at_cells_on_half_levels=exner_at_cells_on_half_levels,
-        temporal_extrapolation_of_perturbed_exner=temporal_extrapolation_of_perturbed_exner,
-        theta_v_at_cells_on_half_levels=theta_v_at_cells_on_half_levels,
+        wgtfacq_c=wgtfacq_c,
+        reference_theta_at_cells_on_half_levels=reference_theta_at_cells_on_half_levels,
+        inv_ddqz_z_full=inv_ddqz_z_full,
+        d2dexdz2_fac1_mc=d2dexdz2_fac1_mc,
+        d2dexdz2_fac2_mc=d2dexdz2_fac2_mc,
+        ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels=ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
+        d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels=d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
         igradp_method=igradp_method,
-        nflatlev=nflatlev,
+        surface_level=surface_level,
+        start_cell_halo_level_2=start_cell_halo_level_2,
+        end_cell_halo_level_2=end_cell_halo_level_2,
         out=(
             perturbed_rho_at_cells_on_model_levels,
             perturbed_theta_v_at_cells_on_model_levels,
@@ -412,64 +351,56 @@ def compute_perturbed_quantities_and_interpolation(
             perturbed_theta_v_at_cells_on_half_levels,
             theta_v_at_cells_on_half_levels,
             pressure_buoyancy_acceleration_at_cells_on_half_levels,
-        ),
-        domain={
-            dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
-            dims.KDim: (model_top, surface_level - 1),
-        },
-    )
-
-    _set_theta_v_and_exner_on_surface_level(
-        temporal_extrapolation_of_perturbed_exner=temporal_extrapolation_of_perturbed_exner,
-        wgtfacq_c=wgtfacq_c,
-        perturbed_theta_v_at_cells_on_model_levels=perturbed_theta_v_at_cells_on_model_levels,
-        reference_theta_at_cells_on_half_levels=reference_theta_at_cells_on_half_levels,
-        out=(
-            perturbed_theta_v_at_cells_on_half_levels,
-            theta_v_at_cells_on_half_levels,
-            exner_at_cells_on_half_levels,
-        ),
-        domain={
-            dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
-            dims.KDim: (surface_level - 1, surface_level),
-        },
-    )
-
-    _compute_first_and_second_vertical_derivative_of_exner(
-        exner_at_cells_on_half_levels=exner_at_cells_on_half_levels,
-        inv_ddqz_z_full=inv_ddqz_z_full,
-        ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels=ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
-        d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels=d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
-        perturbed_theta_v_at_cells_on_half_levels=perturbed_theta_v_at_cells_on_half_levels,
-        d2dexdz2_fac1_mc=d2dexdz2_fac1_mc,
-        d2dexdz2_fac2_mc=d2dexdz2_fac2_mc,
-        perturbed_theta_v_at_cells_on_model_levels=perturbed_theta_v_at_cells_on_model_levels,
-        igradp_method=igradp_method,
-        nflatlev=nflatlev,
-        nflat_gradp=nflat_gradp,
-        out=(
+            temporal_extrapolation_of_perturbed_exner,
             ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
             d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels,
         ),
-        domain={
-            dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
-            dims.KDim: (model_top, surface_level - 1),
-        },
-    )
-
-    _compute_perturbation_of_rho_and_theta(
-        rho=current_rho,
-        rho_ref_mc=reference_rho_at_cells_on_model_levels,
-        theta_v=current_theta_v,
-        theta_ref_mc=reference_theta_at_cells_on_model_levels,
-        out=(
-            perturbed_rho_at_cells_on_model_levels,
-            perturbed_theta_v_at_cells_on_model_levels,
+        domain=(
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (model_top, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (model_top, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (model_top, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (1, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (maximum(1, nflatlev), surface_level),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (1, surface_level),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (1, surface_level),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (1, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (model_top, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (nflatlev, surface_level - 1),
+            },
+            {
+                dims.CellDim: (start_cell_lateral_boundary_level_3, end_cell_halo),
+                dims.KDim: (nflat_gradp, surface_level - 1),
+            },
         ),
-        domain={
-            dims.CellDim: (start_cell_halo_level_2, end_cell_halo_level_2),
-            dims.KDim: (model_top, surface_level - 1),
-        },
     )
 
 
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/interpolate_to_surface.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/extrapolate_quadratically_to_surface.py
similarity index 83%
rename from model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/interpolate_to_surface.py
rename to model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/extrapolate_quadratically_to_surface.py
index 9395878611..7670c5a48b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/interpolate_to_surface.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/extrapolate_quadratically_to_surface.py
@@ -13,21 +13,21 @@
 
 
 @gtx.field_operator
-def _interpolate_to_surface(
+def _extrapolate_quadratically_to_surface(
     wgtfacq_c: fa.CellKField[vpfloat],
     interpolant: fa.CellKField[vpfloat],
 ) -> fa.CellKField[vpfloat]:
     """Formerly known as _mo_solve_nonhydro_stencil_04."""
-    interpolation_to_surface = (
+    extrapolate_quadratically_to_surface = (  # rename extrapolation_to_surface
         wgtfacq_c(Koff[-1]) * interpolant(Koff[-1])
         + wgtfacq_c(Koff[-2]) * interpolant(Koff[-2])
         + wgtfacq_c(Koff[-3]) * interpolant(Koff[-3])
     )
-    return interpolation_to_surface
+    return extrapolate_quadratically_to_surface
 
 
 @gtx.program(grid_type=gtx.GridType.UNSTRUCTURED)
-def interpolate_to_surface(
+def extrapolate_quadratically_to_surface(
     wgtfacq_c: fa.CellKField[vpfloat],
     interpolant: fa.CellKField[vpfloat],
     interpolation_to_surface: fa.CellKField[vpfloat],
@@ -36,7 +36,7 @@ def interpolate_to_surface(
     vertical_start: gtx.int32,
     vertical_end: gtx.int32,
 ):
-    _interpolate_to_surface(
+    _extrapolate_quadratically_to_surface(
         wgtfacq_c,
         interpolant,
         out=interpolation_to_surface,
diff --git a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/set_theta_v_prime_ic_at_lower_boundary.py b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/set_theta_v_prime_ic_at_lower_boundary.py
index 8f48be0c81..58aa6ba91b 100644
--- a/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/set_theta_v_prime_ic_at_lower_boundary.py
+++ b/model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/set_theta_v_prime_ic_at_lower_boundary.py
@@ -8,7 +8,9 @@
 import gt4py.next as gtx
 from gt4py.next import astype
 
-from icon4py.model.atmosphere.dycore.stencils.interpolate_to_surface import _interpolate_to_surface
+from icon4py.model.atmosphere.dycore.stencils.extrapolate_quadratically_to_surface import (
+    _extrapolate_quadratically_to_surface,
+)
 from icon4py.model.common import dimension as dims, field_type_aliases as fa
 from icon4py.model.common.type_alias import vpfloat, wpfloat
 
@@ -20,7 +22,9 @@ def _set_theta_v_prime_ic_at_lower_boundary(
     theta_ref_ic: fa.CellKField[vpfloat],
 ) -> tuple[fa.CellKField[vpfloat], fa.CellKField[wpfloat]]:
     """Formerly known as _mo_solve_nonhydro_stencil_11_upper."""
-    z_theta_v_pr_ic_vp = _interpolate_to_surface(wgtfacq_c=wgtfacq_c, interpolant=z_rth_pr)
+    z_theta_v_pr_ic_vp = _extrapolate_quadratically_to_surface(
+        wgtfacq_c=wgtfacq_c, interpolant=z_rth_pr
+    )
     theta_v_ic_vp = theta_ref_ic + z_theta_v_pr_ic_vp
     return z_theta_v_pr_ic_vp, astype(theta_v_ic_vp, wpfloat)
 
diff --git a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
index e2711f3c82..0afefd5a53 100644
--- a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
+++ b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py
@@ -156,8 +156,6 @@ def test_time_step_flags(
     ],
 )
 def test_nonhydro_predictor_step(
-    istep_init,
-    istep_exit,
     substep_init,
     step_date_init,
     step_date_exit,
@@ -172,7 +170,6 @@ def test_nonhydro_predictor_step(
     interpolation_savepoint,
     savepoint_nonhydro_exit,
     experiment,
-    ndyn_substeps,
     at_initial_timestep,
     caplog,
     backend,
@@ -234,7 +231,7 @@ def test_nonhydro_predictor_step(
     cell_domain = h_grid.domain(dims.CellDim)
     edge_domain = h_grid.domain(dims.EdgeDim)
 
-    cell_start_lateral_boundary_level_2 = icon_grid.start_index(
+    cell_start_lateral_boundary_level_3 = icon_grid.start_index(
         cell_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_3)
     )
 
@@ -250,90 +247,90 @@ def test_nonhydro_predictor_step(
     # stencils 2, 3
     assert test_utils.dallclose(
         diagnostic_state_nh.perturbed_exner_at_cells_on_model_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.exner_pr().asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.exner_pr().asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
     assert test_utils.dallclose(
         solve_nonhydro.temporal_extrapolation_of_perturbed_exner.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.z_exner_ex_pr().asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.z_exner_ex_pr().asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
 
     # stencils 4,5
     assert test_utils.dallclose(
         solve_nonhydro.exner_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, nlev - 1
+            cell_start_lateral_boundary_level_3:, nlev - 1
         ],
-        sp_exit.z_exner_ic().asnumpy()[cell_start_lateral_boundary_level_2:, nlev - 1],
+        sp_exit.z_exner_ic().asnumpy()[cell_start_lateral_boundary_level_3:, nlev - 1],
     )
     nflatlev = vertical_params.nflatlev
     assert test_utils.dallclose(
         solve_nonhydro.exner_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, nflatlev : nlev - 1
+            cell_start_lateral_boundary_level_3:, nflatlev : nlev - 1
         ],
-        sp_exit.z_exner_ic().asnumpy()[cell_start_lateral_boundary_level_2:, nflatlev : nlev - 1],
+        sp_exit.z_exner_ic().asnumpy()[cell_start_lateral_boundary_level_3:, nflatlev : nlev - 1],
         rtol=1.0e-9,
     )
     # stencil 6
     assert test_utils.dallclose(
         solve_nonhydro.ddz_of_temporal_extrapolation_of_perturbed_exner_on_model_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, nflatlev:
+            cell_start_lateral_boundary_level_3:, nflatlev:
         ],
-        sp_exit.z_dexner_dz_c(0).asnumpy()[cell_start_lateral_boundary_level_2:, nflatlev:],
+        sp_exit.z_dexner_dz_c(0).asnumpy()[cell_start_lateral_boundary_level_3:, nflatlev:],
         atol=5e-18,
     )
 
     # stencils 7,8,9
     assert test_utils.dallclose(
         diagnostic_state_nh.rho_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.rho_ic().asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.rho_ic().asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
     assert test_utils.dallclose(
         solve_nonhydro.pressure_buoyancy_acceleration_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, 1:
+            cell_start_lateral_boundary_level_3:, 1:
         ],
-        sp_exit.z_th_ddz_exner_c().asnumpy()[cell_start_lateral_boundary_level_2:, 1:],
+        sp_exit.z_th_ddz_exner_c().asnumpy()[cell_start_lateral_boundary_level_3:, 1:],
         rtol=2.0e-12,
     )
 
     # stencils 7,8,9, 11
     assert test_utils.dallclose(
         solve_nonhydro.perturbed_theta_v_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.z_theta_v_pr_ic().asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.z_theta_v_pr_ic().asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
     assert test_utils.dallclose(
         diagnostic_state_nh.theta_v_at_cells_on_half_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.theta_v_ic().asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.theta_v_ic().asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
     # stencils 7,8,9, 13
     assert test_utils.dallclose(
         solve_nonhydro.perturbed_rho_at_cells_on_model_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.z_rth_pr(0).asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.z_rth_pr(0).asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
     assert test_utils.dallclose(
         solve_nonhydro.perturbed_theta_v_at_cells_on_model_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, :
+            cell_start_lateral_boundary_level_3:, :
         ],
-        sp_exit.z_rth_pr(1).asnumpy()[cell_start_lateral_boundary_level_2:, :],
+        sp_exit.z_rth_pr(1).asnumpy()[cell_start_lateral_boundary_level_3:, :],
     )
 
     # stencils 12
     nflat_gradp = grid_savepoint.nflat_gradp()
     assert test_utils.dallclose(
         solve_nonhydro.d2dz2_of_temporal_extrapolation_of_perturbed_exner_on_model_levels.asnumpy()[
-            cell_start_lateral_boundary_level_2:, nflat_gradp:
+            cell_start_lateral_boundary_level_3:, nflat_gradp:
         ],
-        sp_exit.z_dexner_dz_c(1).asnumpy()[cell_start_lateral_boundary_level_2:, nflat_gradp:],
+        sp_exit.z_dexner_dz_c(1).asnumpy()[cell_start_lateral_boundary_level_3:, nflat_gradp:],
         atol=1e-22,
     )
 
@@ -810,7 +807,6 @@ def test_run_solve_nonhydro_single_step(
     )
 
 
-# why is this not run for APE?
 @pytest.mark.embedded_remap_error
 @pytest.mark.datatest
 @pytest.mark.parametrize("experiment", [definitions.Experiments.MCH_CH_R04B09])
@@ -1023,7 +1019,6 @@ def test_compute_perturbed_quantities_and_interpolation(
     step_date_init,
     step_date_exit,
     *,
-    ndyn_substeps,
     icon_grid,
     lowest_layer_thickness,
     model_top_height,
@@ -1031,9 +1026,6 @@ def test_compute_perturbed_quantities_and_interpolation(
     damping_height,
     grid_savepoint,
     metrics_savepoint,
-    interpolation_savepoint,
-    substep_init,
-    substep_exit,
     savepoint_nonhydro_init,
     savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_init,
     savepoint_nonhydro_exit,
@@ -1091,13 +1083,16 @@ def test_compute_perturbed_quantities_and_interpolation(
     nflat_gradp = grid_savepoint.nflat_gradp()
 
     cell_domain = h_grid.domain(dims.CellDim)
-    start_cell_lateral_boundary = icon_grid.start_index(cell_domain(h_grid.Zone.LATERAL_BOUNDARY))
     start_cell_lateral_boundary_level_3 = icon_grid.start_index(
         cell_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_3)
     )
     start_cell_halo_level_2 = icon_grid.start_index(cell_domain(h_grid.Zone.HALO_LEVEL_2))
     end_cell_halo = icon_grid.end_index(cell_domain(h_grid.Zone.HALO))
     end_cell_halo_level_2 = icon_grid.end_index(cell_domain(h_grid.Zone.HALO_LEVEL_2))
+    print("start_cell_lateral_boundary_level_3", start_cell_lateral_boundary_level_3)
+    print("end_cell_halo", end_cell_halo)
+    print("start_cell_halo_level_2", start_cell_halo_level_2)
+    print("end_cell_halo_level_2", end_cell_halo_level_2)
 
     reference_rho_at_cells_on_model_levels = metrics_savepoint.rho_ref_mc()
     reference_theta_at_cells_on_model_levels = metrics_savepoint.theta_ref_mc()
@@ -1157,7 +1152,6 @@ def test_compute_perturbed_quantities_and_interpolation(
         igradp_method=igradp_method,
         nflatlev=nflatlev,
         nflat_gradp=nflat_gradp,
-        start_cell_lateral_boundary=start_cell_lateral_boundary,
         start_cell_lateral_boundary_level_3=start_cell_lateral_boundary_level_3,
         start_cell_halo_level_2=start_cell_halo_level_2,
         end_cell_halo=end_cell_halo,
@@ -1182,8 +1176,8 @@ def test_compute_perturbed_quantities_and_interpolation(
     assert test_utils.dallclose(
         perturbed_exner_at_cells_on_model_levels.asnumpy(), exner_pr_ref.asnumpy()
     )
-    assert test_utils.dallclose(rho_at_cells_on_half_levels.asnumpy(), rho_ic_ref.asnumpy())
 
+    assert test_utils.dallclose(rho_at_cells_on_half_levels.asnumpy(), rho_ic_ref.asnumpy())
     assert test_utils.dallclose(
         exner_at_cells_on_half_levels.asnumpy()[:, nflatlev:],
         z_exner_ic_ref.asnumpy()[:, nflatlev:],
@@ -1238,17 +1232,8 @@ def test_interpolate_rho_theta_v_to_half_levels_and_compute_pressure_buoyancy_ac
     step_date_init,
     step_date_exit,
     *,
-    ndyn_substeps,
     icon_grid,
-    lowest_layer_thickness,
-    model_top_height,
-    stretch_factor,
-    damping_height,
-    grid_savepoint,
     metrics_savepoint,
-    interpolation_savepoint,
-    substep_init,
-    substep_exit,
     savepoint_nonhydro_init,
     savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_init,
     savepoint_nonhydro_exit,
@@ -1380,7 +1365,6 @@ def test_compute_theta_rho_face_values_and_pressure_gradient_and_update_vn(
     step_date_init,
     step_date_exit,
     *,
-    ndyn_substeps,
     icon_grid,
     savepoint_nonhydro_init,
     lowest_layer_thickness,
@@ -1391,9 +1375,6 @@ def test_compute_theta_rho_face_values_and_pressure_gradient_and_update_vn(
     metrics_savepoint,
     interpolation_savepoint,
     savepoint_nonhydro_exit,
-    istep_init,
-    substep_init,
-    substep_exit,
     savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_init,
     savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_exit,
     backend,
@@ -1596,19 +1577,13 @@ def test_apply_divergence_damping_and_update_vn(
     step_date_init,
     step_date_exit,
     *,
-    ndyn_substeps,
     icon_grid,
     savepoint_nonhydro_init,
-    lowest_layer_thickness,
-    model_top_height,
-    stretch_factor,
-    damping_height,
     grid_savepoint,
     metrics_savepoint,
     interpolation_savepoint,
     savepoint_nonhydro_exit,
     savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_init,
-    savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_exit,
     backend,
 ):
     sp_nh_init = savepoint_nonhydro_init
@@ -1721,15 +1696,10 @@ def test_apply_divergence_damping_and_update_vn(
     ],
 )
 def test_compute_horizontal_velocity_quantities_and_fluxes(
-    istep_init,
-    istep_exit,
-    substep_init,
-    substep_exit,
     step_date_init,
     step_date_exit,
     experiment,
     icon_grid,
-    ndyn_substeps,
     grid_savepoint,
     lowest_layer_thickness,
     model_top_height,
@@ -1739,8 +1709,6 @@ def test_compute_horizontal_velocity_quantities_and_fluxes(
     savepoint_dycore_30_to_38_exit,
     interpolation_savepoint,
     metrics_savepoint,
-    savepoint_nonhydro_init,
-    savepoint_nonhydro_exit,
     backend,
 ):
     edge_domain = h_grid.domain(dims.EdgeDim)
@@ -1903,21 +1871,17 @@ def test_compute_horizontal_velocity_quantities_and_fluxes(
 def test_compute_averaged_vn_and_fluxes_and_prepare_tracer_advection(
     istep_init,
     istep_exit,
-    substep_init,
-    substep_exit,
     step_date_init,
     step_date_exit,
     experiment,
     icon_grid,
     ndyn_substeps,
     at_first_substep,
-    grid_savepoint,
     savepoint_dycore_30_to_38_init,
     savepoint_dycore_30_to_38_exit,
     interpolation_savepoint,
     metrics_savepoint,
     savepoint_nonhydro_init,
-    savepoint_nonhydro_exit,
     backend,
 ):
     edge_domain = h_grid.domain(dims.EdgeDim)
@@ -2026,7 +1990,6 @@ def test_vertically_implicit_solver_at_predictor_step(
     step_date_init,
     step_date_exit,
     *,
-    ndyn_substeps,
     icon_grid,
     savepoint_nonhydro_init,
     lowest_layer_thickness,
@@ -2037,9 +2000,6 @@ def test_vertically_implicit_solver_at_predictor_step(
     metrics_savepoint,
     interpolation_savepoint,
     savepoint_nonhydro_exit,
-    istep_init,
-    istep_exit,
-    substep_exit,
     savepoint_vertically_implicit_dycore_solver_init,
     backend,
 ):
diff --git a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_compute_perturbed_quantities_and_interpolation.py b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_compute_perturbed_quantities_and_interpolation.py
index a173e626f1..c6c9d10921 100644
--- a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_compute_perturbed_quantities_and_interpolation.py
+++ b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_compute_perturbed_quantities_and_interpolation.py
@@ -6,18 +6,6 @@
 # Please, refer to the LICENSE file in the root directory.
 # SPDX-License-Identifier: BSD-3-Clause
 
-# 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
 from typing import Any
 
 import gt4py.next as gtx
@@ -44,11 +32,11 @@
 from .test_compute_virtual_potential_temperatures_and_pressure_gradient import (
     compute_virtual_potential_temperatures_and_pressure_gradient_numpy,
 )
+from .test_extrapolate_quadratically_to_surface import extrapolate_quadratically_to_surface_numpy
 from .test_extrapolate_temporally_exner_pressure import extrapolate_temporally_exner_pressure_numpy
 from .test_interpolate_cell_field_to_half_levels_vp import (
     interpolate_cell_field_to_half_levels_vp_numpy,
 )
-from .test_interpolate_to_surface import interpolate_to_surface_numpy
 from .test_set_theta_v_prime_ic_at_lower_boundary import (
     set_theta_v_prime_ic_at_lower_boundary_numpy,
 )
@@ -87,7 +75,6 @@ class TestComputePerturbedQuantitiesAndInterpolation(stencil_tests.StencilTest):
             "igradp_method",
             "nflatlev",
             "nflat_gradp",
-            "start_cell_lateral_boundary",
             "start_cell_lateral_boundary_level_3",
             "start_cell_halo_level_2",
             "end_cell_halo",
@@ -137,7 +124,6 @@ def reference(
         igradp_method: gtx.int32,
         nflatlev: gtx.int32,
         nflat_gradp: gtx.int32,
-        start_cell_lateral_boundary: gtx.int32,
         start_cell_lateral_boundary_level_3: gtx.int32,
         start_cell_halo_level_2: gtx.int32,
         end_cell_halo: gtx.int32,
@@ -153,8 +139,7 @@ def reference(
             perturbed_rho_at_cells_on_model_levels,
             perturbed_theta_v_at_cells_on_model_levels[:, : surface_level - 1],
         ) = np.where(
-            (start_cell_lateral_boundary <= horz_idx)
-            & (horz_idx < start_cell_lateral_boundary_level_3),
+            horz_idx < start_cell_lateral_boundary_level_3,
             (
                 np.zeros_like(perturbed_rho_at_cells_on_model_levels),
                 np.zeros_like(perturbed_theta_v_at_cells_on_model_levels[:, : surface_level - 1]),
@@ -195,7 +180,7 @@ def reference(
                 (start_cell_lateral_boundary_level_3 <= horz_idx)
                 & (horz_idx < end_cell_halo)
                 & (vert_idx == surface_level - 1),
-                interpolate_to_surface_numpy(
+                extrapolate_quadratically_to_surface_numpy(
                     interpolant=temporal_extrapolation_of_perturbed_exner,
                     wgtfacq_c=wgtfacq_c,
                     interpolation_to_surface=exner_at_cells_on_half_levels,
@@ -412,7 +397,6 @@ def input_data(self, grid: base.Grid) -> dict[str, gtx.Field | state_utils.Scala
         igradp_method = horzpres_discr_type.TAYLOR_HYDRO
 
         cell_domain = h_grid.domain(dims.CellDim)
-        start_cell_lateral_boundary = grid.start_index(cell_domain(h_grid.Zone.LATERAL_BOUNDARY))
         start_cell_lateral_boundary_level_3 = grid.start_index(
             cell_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_3)
         )
@@ -454,7 +438,6 @@ def input_data(self, grid: base.Grid) -> dict[str, gtx.Field | state_utils.Scala
             igradp_method=igradp_method,
             nflatlev=nflatlev,
             nflat_gradp=nflat_gradp,
-            start_cell_lateral_boundary=start_cell_lateral_boundary,
             start_cell_lateral_boundary_level_3=start_cell_lateral_boundary_level_3,
             start_cell_halo_level_2=start_cell_halo_level_2,
             end_cell_halo=end_cell_halo,
diff --git a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_interpolate_to_surface.py b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_extrapolate_quadratically_to_surface.py
similarity index 88%
rename from model/atmosphere/dycore/tests/dycore/stencil_tests/test_interpolate_to_surface.py
rename to model/atmosphere/dycore/tests/dycore/stencil_tests/test_extrapolate_quadratically_to_surface.py
index 2dc775909d..d85b91728f 100644
--- a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_interpolate_to_surface.py
+++ b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_extrapolate_quadratically_to_surface.py
@@ -11,7 +11,9 @@
 import numpy as np
 import pytest
 
-from icon4py.model.atmosphere.dycore.stencils.interpolate_to_surface import interpolate_to_surface
+from icon4py.model.atmosphere.dycore.stencils.extrapolate_quadratically_to_surface import (
+    extrapolate_quadratically_to_surface,
+)
 from icon4py.model.common import dimension as dims
 from icon4py.model.common.grid import base
 from icon4py.model.common.states import utils as state_utils
@@ -20,7 +22,7 @@
 from icon4py.model.testing.stencil_tests import StencilTest
 
 
-def interpolate_to_surface_numpy(
+def extrapolate_quadratically_to_surface_numpy(
     interpolant: np.ndarray, wgtfacq_c: np.ndarray, interpolation_to_surface: np.ndarray
 ) -> np.ndarray:
     interpolation_to_surface = np.copy(interpolation_to_surface)
@@ -33,7 +35,7 @@ def interpolate_to_surface_numpy(
 
 
 class TestInterpolateToSurface(StencilTest):
-    PROGRAM = interpolate_to_surface
+    PROGRAM = extrapolate_quadratically_to_surface
     OUTPUTS = ("interpolation_to_surface",)
 
     @staticmethod
@@ -44,7 +46,7 @@ def reference(
         interpolation_to_surface: np.ndarray,
         **kwargs: Any,
     ) -> dict:
-        interpolation_to_surface = interpolate_to_surface_numpy(
+        interpolation_to_surface = extrapolate_quadratically_to_surface_numpy(
             wgtfacq_c=wgtfacq_c,
             interpolant=interpolant,
             interpolation_to_surface=interpolation_to_surface,
diff --git a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_set_theta_v_prime_ic_at_lower_boundary.py b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_set_theta_v_prime_ic_at_lower_boundary.py
index 122f1913ef..8198cbbf6d 100644
--- a/model/atmosphere/dycore/tests/dycore/stencil_tests/test_set_theta_v_prime_ic_at_lower_boundary.py
+++ b/model/atmosphere/dycore/tests/dycore/stencil_tests/test_set_theta_v_prime_ic_at_lower_boundary.py
@@ -21,7 +21,7 @@
 from icon4py.model.common.utils.data_allocation import random_field, zero_field
 from icon4py.model.testing.stencil_tests import StencilTest
 
-from .test_interpolate_to_surface import interpolate_to_surface_numpy
+from .test_extrapolate_quadratically_to_surface import extrapolate_quadratically_to_surface_numpy
 
 
 def set_theta_v_prime_ic_at_lower_boundary_numpy(
@@ -31,7 +31,7 @@ def set_theta_v_prime_ic_at_lower_boundary_numpy(
     z_theta_v_pr_ic: np.ndarray,
     theta_v_ic: np.ndarray,
 ) -> tuple[np.ndarray, ...]:
-    z_theta_v_pr_ic = interpolate_to_surface_numpy(
+    z_theta_v_pr_ic = extrapolate_quadratically_to_surface_numpy(
         wgtfacq_c=wgtfacq_c,
         interpolant=z_rth_pr,
         interpolation_to_surface=z_theta_v_pr_ic,
diff --git a/model/common/tests/common/interpolation/stencil_tests/test_mo_intp_rbf_rbf_vec_interpol_vertex.py b/model/common/tests/common/interpolation/stencil_tests/test_mo_intp_rbf_rbf_vec_interpol_vertex.py
index 07104cbd89..3093154f5b 100644
--- a/model/common/tests/common/interpolation/stencil_tests/test_mo_intp_rbf_rbf_vec_interpol_vertex.py
+++ b/model/common/tests/common/interpolation/stencil_tests/test_mo_intp_rbf_rbf_vec_interpol_vertex.py
@@ -48,7 +48,7 @@ def reference(
         horizontal_start: int,
         horizontal_end: int,
         **kwargs: Any,
-    ) -> dict[str, np.ndarray]:
+    ) -> dict:
         v2e = connectivities[dims.V2EDim]
         ptr_coeff_1 = np.expand_dims(ptr_coeff_1, axis=-1)
         p_u_out = np.sum(
diff --git a/model/common/tests/common/math/unit_tests/test_operators.py b/model/common/tests/common/math/unit_tests/test_operators.py
index b487d1a50b..3c75a92f9b 100644
--- a/model/common/tests/common/math/unit_tests/test_operators.py
+++ b/model/common/tests/common/math/unit_tests/test_operators.py
@@ -33,7 +33,7 @@ def reference(
         psi_c: np.ndarray,
         geofac_n2s: np.ndarray,
         **kwargs: Any,
-    ) -> dict[str, np.ndarray]:
+    ) -> dict:
         nabla2_psi_c_np = reference_funcs.nabla2_on_cell_numpy(connectivities, psi_c, geofac_n2s)
         return dict(nabla2_psi_c=nabla2_psi_c_np)