From 03621002c2478a2e2526c563cef84b38d92941d9 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Wed, 10 Sep 2025 15:38:05 +0200
Subject: [PATCH 1/3] Also using _assert_same_function_signature for Kernels

---
 parcels/field.py           | 29 +++++------------------------
 parcels/kernel.py          |  3 +++
 parcels/tools/_helpers.py  | 20 ++++++++++++++++++++
 tests/v4/test_advection.py |  2 +-
 4 files changed, 29 insertions(+), 25 deletions(-)

diff --git a/parcels/field.py b/parcels/field.py
index b248414a4a..ce437c4713 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -1,6 +1,5 @@
 from __future__ import annotations
 
-import inspect
 import warnings
 from collections.abc import Callable
 from datetime import datetime
@@ -19,6 +18,7 @@
     ZeroInterpolator_Vector,
 )
 from parcels.particle import KernelParticle
+from parcels.tools._helpers import _assert_same_function_signature
 from parcels.tools.converters import (
     UnitConverter,
     unitconverters_map,
@@ -57,25 +57,6 @@ def _deal_with_errors(error, key, vector_type: VectorType):
 }
 
 
-def _assert_same_function_signature(f: Callable, *, ref: Callable) -> None:
-    """Ensures a function `f` has the same signature as the reference function `ref`."""
-    sig_ref = inspect.signature(ref)
-    sig = inspect.signature(f)
-
-    if len(sig_ref.parameters) != len(sig.parameters):
-        raise ValueError(
-            f"Interpolation function must have {len(sig_ref.parameters)} parameters, got {len(sig.parameters)}"
-        )
-
-    for (_name1, param1), (_name2, param2) in zip(sig_ref.parameters.items(), sig.parameters.items(), strict=False):
-        if param1.kind != param2.kind:
-            raise ValueError(
-                f"Parameter '{_name2}' has incorrect parameter kind. Expected {param1.kind}, got {param2.kind}"
-            )
-        if param1.name != param2.name:
-            raise ValueError(f"Parameter '{_name2}' has incorrect name. Expected '{param1.name}', got '{param2.name}'")
-
-
 class Field:
     """The Field class that holds scalar field data.
     The `Field` object is a wrapper around a xarray.DataArray or uxarray.UxDataArray object.
@@ -157,7 +138,7 @@ def __init__(
         if interp_method is None:
             self._interp_method = _DEFAULT_INTERPOLATOR_MAPPING[type(self.grid)]
         else:
-            _assert_same_function_signature(interp_method, ref=ZeroInterpolator)
+            _assert_same_function_signature(interp_method, ref=ZeroInterpolator, context="Interpolation")
             self._interp_method = interp_method
 
         self.igrid = -1  # Default the grid index to -1
@@ -213,7 +194,7 @@ def interp_method(self):
 
     @interp_method.setter
     def interp_method(self, method: Callable):
-        _assert_same_function_signature(method, ref=ZeroInterpolator)
+        _assert_same_function_signature(method, ref=ZeroInterpolator, context="Interpolation")
         self._interp_method = method
 
     def _check_velocitysampling(self):
@@ -287,7 +268,7 @@ def __init__(
         if vector_interp_method is None:
             self._vector_interp_method = None
         else:
-            _assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector)
+            _assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation")
             self._vector_interp_method = vector_interp_method
 
     def __repr__(self):
@@ -303,7 +284,7 @@ def vector_interp_method(self):
 
     @vector_interp_method.setter
     def vector_interp_method(self, method: Callable):
-        _assert_same_function_signature(method, ref=ZeroInterpolator_Vector)
+        _assert_same_function_signature(method, ref=ZeroInterpolator_Vector, context="Interpolation")
         self._vector_interp_method = method
 
     def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 41409c15ab..612374fbdf 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -13,6 +13,7 @@
     AdvectionRK45,
 )
 from parcels.basegrid import GridType
+from parcels.tools._helpers import _assert_same_function_signature
 from parcels.tools.statuscodes import (
     StatusCode,
     _raise_field_interpolation_error,
@@ -23,6 +24,7 @@
     _raise_time_extrapolation_error,
 )
 from parcels.tools.warnings import KernelWarning
+from tests.common_kernels import DoNothing
 
 if TYPE_CHECKING:
     from collections.abc import Callable
@@ -67,6 +69,7 @@ def __init__(
         for f in pyfuncs:
             if not isinstance(f, types.FunctionType):
                 raise TypeError(f"Argument pyfunc should be a function or list of functions. Got {type(f)}")
+            _assert_same_function_signature(f, ref=DoNothing, context="Kernel")
 
         if len(pyfuncs) == 0:
             raise ValueError("List of `pyfuncs` should have at least one function.")
diff --git a/parcels/tools/_helpers.py b/parcels/tools/_helpers.py
index 4e3500e89e..22b1c0c143 100644
--- a/parcels/tools/_helpers.py
+++ b/parcels/tools/_helpers.py
@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import functools
+import inspect
 import warnings
 from collections.abc import Callable
 from datetime import timedelta
@@ -75,3 +76,22 @@ def timedelta_to_float(dt: float | timedelta | np.timedelta64) -> float:
 def should_calculate_next_ti(ti: int, tau: float, tdim: int):
     """Check if the time is beyond the last time in the field"""
     return np.greater(tau, 0) and ti < tdim - 1
+
+
+def _assert_same_function_signature(f: Callable, *, ref: Callable, context: str) -> None:
+    """Ensures a function `f` has the same signature as the reference function `ref`."""
+    sig_ref = inspect.signature(ref)
+    sig = inspect.signature(f)
+
+    if len(sig_ref.parameters) != len(sig.parameters):
+        raise ValueError(
+            f"{context} function must have {len(sig_ref.parameters)} parameters, got {len(sig.parameters)}"
+        )
+
+    for (_name1, param1), (_name2, param2) in zip(sig_ref.parameters.items(), sig.parameters.items(), strict=False):
+        if param1.kind != param2.kind:
+            raise ValueError(
+                f"Parameter '{_name2}' has incorrect parameter kind. Expected {param1.kind}, got {param2.kind}"
+            )
+        if param1.name != param2.name:
+            raise ValueError(f"Parameter '{_name2}' has incorrect name. Expected '{param1.name}', got '{param2.name}'")
diff --git a/tests/v4/test_advection.py b/tests/v4/test_advection.py
index a8bae96862..ef3c8cb3a6 100644
--- a/tests/v4/test_advection.py
+++ b/tests/v4/test_advection.py
@@ -462,7 +462,7 @@ def test_nemo_curvilinear_fieldset():
     latp = np.linspace(-70, 88, npart)
     runtime = np.timedelta64(12, "h")  # TODO increase to 160 days
 
-    def periodicBC(particle, fieldSet, time):  # pragma: no cover
+    def periodicBC(particle, fieldset, time):  # pragma: no cover
         particle.dlon = np.where(particle.lon > 180, particle.dlon - 360, particle.dlon)
 
     pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)

From 017d9ccc7ad61c249899d15ea715891ac87af69b Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Wed, 10 Sep 2025 16:00:52 +0200
Subject: [PATCH 2/3] Changing Kernel API throughout tests

---
 docs/community/v4-migration.md                |   3 +-
 parcels/application_kernels/advection.py      | 228 +++++++++---------
 .../application_kernels/advectiondiffusion.py |  68 +++---
 parcels/kernel.py                             |  30 +--
 tests/common_kernels.py                       |  19 +-
 tests/v4/test_advection.py                    |  50 ++--
 tests/v4/test_diffusion.py                    |   8 +-
 tests/v4/test_kernel.py                       |   4 +-
 tests/v4/test_particlefile.py                 |  46 ++--
 tests/v4/test_particleset.py                  |   4 +-
 tests/v4/test_particleset_execute.py          |  66 ++---
 11 files changed, 270 insertions(+), 256 deletions(-)

diff --git a/docs/community/v4-migration.md b/docs/community/v4-migration.md
index 3635d88a01..7a5bf03659 100644
--- a/docs/community/v4-migration.md
+++ b/docs/community/v4-migration.md
@@ -11,7 +11,8 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
 - Sharing state between kernels must be done via the particle data (as the kernels are not combined under the hood anymore).
 - `particl_dlon`, `particle_dlat` etc have been renamed to `particle.dlon` and `particle.dlat`.
 - `particle.dt` is a np.timedelta64 object; be careful when multiplying `particle.dt` with a velocity, as its value may be cast to nanoseconds.
-- The `time` argument in the Kernel signature is now standard `None` (and may be removed in the Kernel API before release of v4), so can't be used. Use `particle.time` instead.
+- The `time` argument in the Kernel signature has been removed in the Kernel API, so can't be used. Use `particle.time` instead.
+- The `particle` argument in the Kernel signature has been renamed to `particles`.
 - `math` functions should be replaced with array compatible equivalents (e.g., `math.sin` -> `np.sin`). Instead of `ParcelsRandom` you should use numpy's random functions.
 
 ## FieldSet
diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py
index 928da7d53a..1023f8d88e 100644
--- a/parcels/application_kernels/advection.py
+++ b/parcels/application_kernels/advection.py
@@ -16,96 +16,96 @@
 ]
 
 
-def AdvectionRK4(particle, fieldset, time):  # pragma: no cover
+def AdvectionRK4(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
-    (u1, v1) = fieldset.UV[particle]
-    lon1, lat1 = (particle.lon + u1 * 0.5 * dt, particle.lat + v1 * 0.5 * dt)
-    (u2, v2) = fieldset.UV[particle.time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle]
-    lon2, lat2 = (particle.lon + u2 * 0.5 * dt, particle.lat + v2 * 0.5 * dt)
-    (u3, v3) = fieldset.UV[particle.time + 0.5 * particle.dt, particle.depth, lat2, lon2, particle]
-    lon3, lat3 = (particle.lon + u3 * dt, particle.lat + v3 * dt)
-    (u4, v4) = fieldset.UV[particle.time + particle.dt, particle.depth, lat3, lon3, particle]
-    particle.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
-    particle.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
-
-
-def AdvectionRK4_3D(particle, fieldset, time):  # pragma: no cover
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    (u1, v1) = fieldset.UV[particles]
+    lon1, lat1 = (particles.lon + u1 * 0.5 * dt, particles.lat + v1 * 0.5 * dt)
+    (u2, v2) = fieldset.UV[particles.time + 0.5 * particles.dt, particles.depth, lat1, lon1, particles]
+    lon2, lat2 = (particles.lon + u2 * 0.5 * dt, particles.lat + v2 * 0.5 * dt)
+    (u3, v3) = fieldset.UV[particles.time + 0.5 * particles.dt, particles.depth, lat2, lon2, particles]
+    lon3, lat3 = (particles.lon + u3 * dt, particles.lat + v3 * dt)
+    (u4, v4) = fieldset.UV[particles.time + particles.dt, particles.depth, lat3, lon3, particles]
+    particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
+    particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
+
+
+def AdvectionRK4_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
-    dt = particle.dt / np.timedelta64(1, "s")
-    (u1, v1, w1) = fieldset.UVW[particle]
-    lon1 = particle.lon + u1 * 0.5 * dt
-    lat1 = particle.lat + v1 * 0.5 * dt
-    dep1 = particle.depth + w1 * 0.5 * dt
-    (u2, v2, w2) = fieldset.UVW[particle.time + 0.5 * particle.dt, dep1, lat1, lon1, particle]
-    lon2 = particle.lon + u2 * 0.5 * dt
-    lat2 = particle.lat + v2 * 0.5 * dt
-    dep2 = particle.depth + w2 * 0.5 * dt
-    (u3, v3, w3) = fieldset.UVW[particle.time + 0.5 * particle.dt, dep2, lat2, lon2, particle]
-    lon3 = particle.lon + u3 * dt
-    lat3 = particle.lat + v3 * dt
-    dep3 = particle.depth + w3 * dt
-    (u4, v4, w4) = fieldset.UVW[particle.time + particle.dt, dep3, lat3, lon3, particle]
-    particle.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
-    particle.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
-    particle.ddepth += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt
-
-
-def AdvectionRK4_3D_CROCO(particle, fieldset, time):  # pragma: no cover
+    dt = particles.dt / np.timedelta64(1, "s")
+    (u1, v1, w1) = fieldset.UVW[particles]
+    lon1 = particles.lon + u1 * 0.5 * dt
+    lat1 = particles.lat + v1 * 0.5 * dt
+    dep1 = particles.depth + w1 * 0.5 * dt
+    (u2, v2, w2) = fieldset.UVW[particles.time + 0.5 * particles.dt, dep1, lat1, lon1, particles]
+    lon2 = particles.lon + u2 * 0.5 * dt
+    lat2 = particles.lat + v2 * 0.5 * dt
+    dep2 = particles.depth + w2 * 0.5 * dt
+    (u3, v3, w3) = fieldset.UVW[particles.time + 0.5 * particles.dt, dep2, lat2, lon2, particles]
+    lon3 = particles.lon + u3 * dt
+    lat3 = particles.lat + v3 * dt
+    dep3 = particles.depth + w3 * dt
+    (u4, v4, w4) = fieldset.UVW[particles.time + particles.dt, dep3, lat3, lon3, particles]
+    particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
+    particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
+    particles.ddepth += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt
+
+
+def AdvectionRK4_3D_CROCO(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity.
     This kernel assumes the vertical velocity is the 'w' field from CROCO output and works on sigma-layers.
     """
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
-    sig_dep = particle.depth / fieldset.H[particle.time, 0, particle.lat, particle.lon]
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    sig_dep = particles.depth / fieldset.H[particles.time, 0, particles.lat, particles.lon]
 
-    (u1, v1, w1) = fieldset.UVW[particle.time, particle.depth, particle.lat, particle.lon, particle]
-    w1 *= sig_dep / fieldset.H[particle.time, 0, particle.lat, particle.lon]
-    lon1 = particle.lon + u1 * 0.5 * dt
-    lat1 = particle.lat + v1 * 0.5 * dt
+    (u1, v1, w1) = fieldset.UVW[particles.time, particles.depth, particles.lat, particles.lon, particles]
+    w1 *= sig_dep / fieldset.H[particles.time, 0, particles.lat, particles.lon]
+    lon1 = particles.lon + u1 * 0.5 * dt
+    lat1 = particles.lat + v1 * 0.5 * dt
     sig_dep1 = sig_dep + w1 * 0.5 * dt
-    dep1 = sig_dep1 * fieldset.H[particle.time, 0, lat1, lon1]
+    dep1 = sig_dep1 * fieldset.H[particles.time, 0, lat1, lon1]
 
-    (u2, v2, w2) = fieldset.UVW[particle.time + 0.5 * particle.dt, dep1, lat1, lon1, particle]
-    w2 *= sig_dep1 / fieldset.H[particle.time, 0, lat1, lon1]
-    lon2 = particle.lon + u2 * 0.5 * dt
-    lat2 = particle.lat + v2 * 0.5 * dt
+    (u2, v2, w2) = fieldset.UVW[particles.time + 0.5 * particles.dt, dep1, lat1, lon1, particles]
+    w2 *= sig_dep1 / fieldset.H[particles.time, 0, lat1, lon1]
+    lon2 = particles.lon + u2 * 0.5 * dt
+    lat2 = particles.lat + v2 * 0.5 * dt
     sig_dep2 = sig_dep + w2 * 0.5 * dt
-    dep2 = sig_dep2 * fieldset.H[particle.time, 0, lat2, lon2]
+    dep2 = sig_dep2 * fieldset.H[particles.time, 0, lat2, lon2]
 
-    (u3, v3, w3) = fieldset.UVW[particle.time + 0.5 * particle.dt, dep2, lat2, lon2, particle]
-    w3 *= sig_dep2 / fieldset.H[particle.time, 0, lat2, lon2]
-    lon3 = particle.lon + u3 * dt
-    lat3 = particle.lat + v3 * dt
+    (u3, v3, w3) = fieldset.UVW[particles.time + 0.5 * particles.dt, dep2, lat2, lon2, particles]
+    w3 *= sig_dep2 / fieldset.H[particles.time, 0, lat2, lon2]
+    lon3 = particles.lon + u3 * dt
+    lat3 = particles.lat + v3 * dt
     sig_dep3 = sig_dep + w3 * dt
-    dep3 = sig_dep3 * fieldset.H[particle.time, 0, lat3, lon3]
+    dep3 = sig_dep3 * fieldset.H[particles.time, 0, lat3, lon3]
 
-    (u4, v4, w4) = fieldset.UVW[particle.time + particle.dt, dep3, lat3, lon3, particle]
-    w4 *= sig_dep3 / fieldset.H[particle.time, 0, lat3, lon3]
-    lon4 = particle.lon + u4 * dt
-    lat4 = particle.lat + v4 * dt
+    (u4, v4, w4) = fieldset.UVW[particles.time + particles.dt, dep3, lat3, lon3, particles]
+    w4 *= sig_dep3 / fieldset.H[particles.time, 0, lat3, lon3]
+    lon4 = particles.lon + u4 * dt
+    lat4 = particles.lat + v4 * dt
     sig_dep4 = sig_dep + w4 * dt
-    dep4 = sig_dep4 * fieldset.H[particle.time, 0, lat4, lon4]
-
-    particle.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
-    particle.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
-    particle.ddepth += (
-        (dep1 - particle.depth) * 2
-        + 2 * (dep2 - particle.depth) * 2
-        + 2 * (dep3 - particle.depth)
+    dep4 = sig_dep4 * fieldset.H[particles.time, 0, lat4, lon4]
+
+    particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
+    particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
+    particles.ddepth += (
+        (dep1 - particles.depth) * 2
+        + 2 * (dep2 - particles.depth) * 2
+        + 2 * (dep3 - particles.depth)
         + dep4
-        - particle.depth
+        - particles.depth
     ) / 6
 
 
-def AdvectionEE(particle, fieldset, time):  # pragma: no cover
+def AdvectionEE(particles, fieldset):  # pragma: no cover
     """Advection of particles using Explicit Euler (aka Euler Forward) integration."""
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
-    (u1, v1) = fieldset.UV[particle]
-    particle.dlon += u1 * dt
-    particle.dlat += v1 * dt
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    (u1, v1) = fieldset.UV[particles]
+    particles.dlon += u1 * dt
+    particles.dlat += v1 * dt
 
 
-def AdvectionRK45(particle, fieldset, time):  # pragma: no cover
+def AdvectionRK45(particles, fieldset):  # pragma: no cover
     """Advection of particles using adaptive Runge-Kutta 4/5 integration.
 
     Note that this kernel requires a Particle Class that has an extra Variable 'next_dt'
@@ -115,7 +115,7 @@ def AdvectionRK45(particle, fieldset, time):  # pragma: no cover
     Time-step dt is halved if error is larger than fieldset.RK45_tol,
     and doubled if error is smaller than 1/10th of tolerance.
     """
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
 
     c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
     A = [
@@ -128,29 +128,29 @@ def AdvectionRK45(particle, fieldset, time):  # pragma: no cover
     b4 = [25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, -1.0 / 5.0]
     b5 = [16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0]
 
-    (u1, v1) = fieldset.UV[particle]
-    lon1, lat1 = (particle.lon + u1 * A[0][0] * dt, particle.lat + v1 * A[0][0] * dt)
-    (u2, v2) = fieldset.UV[particle.time + c[0] * particle.dt, particle.depth, lat1, lon1, particle]
+    (u1, v1) = fieldset.UV[particles]
+    lon1, lat1 = (particles.lon + u1 * A[0][0] * dt, particles.lat + v1 * A[0][0] * dt)
+    (u2, v2) = fieldset.UV[particles.time + c[0] * particles.dt, particles.depth, lat1, lon1, particles]
     lon2, lat2 = (
-        particle.lon + (u1 * A[1][0] + u2 * A[1][1]) * dt,
-        particle.lat + (v1 * A[1][0] + v2 * A[1][1]) * dt,
+        particles.lon + (u1 * A[1][0] + u2 * A[1][1]) * dt,
+        particles.lat + (v1 * A[1][0] + v2 * A[1][1]) * dt,
     )
-    (u3, v3) = fieldset.UV[particle.time + c[1] * particle.dt, particle.depth, lat2, lon2, particle]
+    (u3, v3) = fieldset.UV[particles.time + c[1] * particles.dt, particles.depth, lat2, lon2, particles]
     lon3, lat3 = (
-        particle.lon + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * dt,
-        particle.lat + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * dt,
+        particles.lon + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * dt,
+        particles.lat + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * dt,
     )
-    (u4, v4) = fieldset.UV[particle.time + c[2] * particle.dt, particle.depth, lat3, lon3, particle]
+    (u4, v4) = fieldset.UV[particles.time + c[2] * particles.dt, particles.depth, lat3, lon3, particles]
     lon4, lat4 = (
-        particle.lon + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * dt,
-        particle.lat + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * dt,
+        particles.lon + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * dt,
+        particles.lat + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * dt,
     )
-    (u5, v5) = fieldset.UV[particle.time + c[3] * particle.dt, particle.depth, lat4, lon4, particle]
+    (u5, v5) = fieldset.UV[particles.time + c[3] * particles.dt, particles.depth, lat4, lon4, particles]
     lon5, lat5 = (
-        particle.lon + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * dt,
-        particle.lat + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * dt,
+        particles.lon + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * dt,
+        particles.lat + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * dt,
     )
-    (u6, v6) = fieldset.UV[particle.time + c[4] * particle.dt, particle.depth, lat5, lon5, particle]
+    (u6, v6) = fieldset.UV[particles.time + c[4] * particles.dt, particles.depth, lat5, lon5, particles]
 
     lon_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * dt
     lat_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * dt
@@ -160,31 +160,31 @@ def AdvectionRK45(particle, fieldset, time):  # pragma: no cover
     kappa = np.sqrt(np.pow(lon_5th - lon_4th, 2) + np.pow(lat_5th - lat_4th, 2))
 
     good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(dt) <= np.fabs(fieldset.RK45_min_dt))
-    particle.dlon += np.where(good_particles, lon_5th, 0)
-    particle.dlat += np.where(good_particles, lat_5th, 0)
+    particles.dlon += np.where(good_particles, lon_5th, 0)
+    particles.dlat += np.where(good_particles, lat_5th, 0)
 
     increase_dt_particles = (
         good_particles & (kappa <= fieldset.RK45_tol / 10) & (np.fabs(dt * 2) <= np.fabs(fieldset.RK45_max_dt))
     )
-    particle.dt = np.where(increase_dt_particles, particle.dt * 2, particle.dt)
-    particle.dt = np.where(
-        particle.dt > fieldset.RK45_max_dt * np.timedelta64(1, "s"),
+    particles.dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
+    particles.dt = np.where(
+        particles.dt > fieldset.RK45_max_dt * np.timedelta64(1, "s"),
         fieldset.RK45_max_dt * np.timedelta64(1, "s"),
-        particle.dt,
+        particles.dt,
     )
-    particle.state = np.where(good_particles, StatusCode.Success, particle.state)
+    particles.state = np.where(good_particles, StatusCode.Success, particles.state)
 
     repeat_particles = np.invert(good_particles)
-    particle.dt = np.where(repeat_particles, particle.dt / 2, particle.dt)
-    particle.dt = np.where(
-        particle.dt < fieldset.RK45_min_dt * np.timedelta64(1, "s"),
+    particles.dt = np.where(repeat_particles, particles.dt / 2, particles.dt)
+    particles.dt = np.where(
+        particles.dt < fieldset.RK45_min_dt * np.timedelta64(1, "s"),
         fieldset.RK45_min_dt * np.timedelta64(1, "s"),
-        particle.dt,
+        particles.dt,
     )
-    particle.state = np.where(repeat_particles, StatusCode.Repeat, particle.state)
+    particles.state = np.where(repeat_particles, StatusCode.Repeat, particles.state)
 
 
-def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
+def AdvectionAnalytical(particles, fieldset):  # pragma: no cover
     """Advection of particles using 'analytical advection' integration.
 
     Based on Ariane/TRACMASS algorithm, as detailed in e.g. Doos et al (https://doi.org/10.5194/gmd-10-1733-2017).
@@ -197,17 +197,17 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
 
     tol = 1e-10
     I_s = 10  # number of intermediate time steps
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO improve API for converting dt to seconds
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO improve API for converting dt to seconds
     direction = 1.0 if dt > 0 else -1.0
     withW = True if "W" in [f.name for f in fieldset.fields.values()] else False
     withTime = True if len(fieldset.U.grid.time) > 1 else False
     tau, zeta, eta, xsi, ti, zi, yi, xi = fieldset.U._search_indices(
-        time, particle.depth, particle.lat, particle.lon, particle=particle
+        particles.depth, particles.lat, particles.lon, particles=particles
     )
     ds_t = dt
     if withTime:
         time_i = np.linspace(0, fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti], I_s)
-        ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
+        ds_t = min(ds_t, time_i[np.where(particles.time - fieldset.U.grid.time[ti] < time_i)[0][0]])
 
     if withW:
         if abs(xsi - 1) < tol:
@@ -232,7 +232,7 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
                 yi += 1
                 eta = 0
 
-    particle.ei[:] = fieldset.U.ravel_index(zi, yi, xi)
+    particles.ei[:] = fieldset.U.ravel_index(zi, yi, xi)
 
     grid = fieldset.U.grid
     if grid._gtype < 2:
@@ -242,8 +242,8 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
         px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
         py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
     if grid.mesh == "spherical":
-        px[0] = px[0] + 360 if px[0] < particle.lon - 225 else px[0]
-        px[0] = px[0] - 360 if px[0] > particle.lat + 225 else px[0]
+        px[0] = px[0] + 360 if px[0] < particles.lon - 225 else px[0]
+        px[0] = px[0] - 360 if px[0] > particles.lat + 225 else px[0]
         px[1:] = np.where(px[1:] - px[0] > 180, px[1:] - 360, px[1:])
         px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:])
     if withW:
@@ -258,7 +258,7 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
     c4 = i_u._geodetic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
     rad = np.pi / 180.0
     deg2m = 1852 * 60.0
-    meshJac = (deg2m * deg2m * math.cos(rad * particle.lat)) if grid.mesh == "spherical" else 1
+    meshJac = (deg2m * deg2m * math.cos(rad * particles.lat)) if grid.mesh == "spherical" else 1
     dxdy = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
 
     if withW:
@@ -333,26 +333,26 @@ def compute_rs(r, B, delta, s_min):
     rs_x = compute_rs(xsi, B_x, delta_x, s_min)
     rs_y = compute_rs(eta, B_y, delta_y, s_min)
 
-    particle.dlon += (
+    particles.dlon += (
         (1.0 - rs_x) * (1.0 - rs_y) * px[0]
         + rs_x * (1.0 - rs_y) * px[1]
         + rs_x * rs_y * px[2]
         + (1.0 - rs_x) * rs_y * px[3]
-        - particle.lon
+        - particles.lon
     )
-    particle.dlat += (
+    particles.dlat += (
         (1.0 - rs_x) * (1.0 - rs_y) * py[0]
         + rs_x * (1.0 - rs_y) * py[1]
         + rs_x * rs_y * py[2]
         + (1.0 - rs_x) * rs_y * py[3]
-        - particle.lat
+        - particles.lat
     )
 
     if withW:
         rs_z = compute_rs(zeta, B_z, delta_z, s_min)
-        particle.ddepth += (1.0 - rs_z) * pz[0] + rs_z * pz[1] - particle.depth
+        particles.ddepth += (1.0 - rs_z) * pz[0] + rs_z * pz[1] - particles.depth
 
-    if particle.dt > 0:
-        particle.dt = max(direction * s_min * (dxdy * dz), 1e-7).astype("timedelta64[s]")
+    if particles.dt > 0:
+        particles.dt = max(direction * s_min * (dxdy * dz), 1e-7).astype("timedelta64[s]")
     else:
-        particle.dt = min(direction * s_min * (dxdy * dz), -1e-7).astype("timedelta64[s]")
+        particles.dt = min(direction * s_min * (dxdy * dz), -1e-7).astype("timedelta64[s]")
diff --git a/parcels/application_kernels/advectiondiffusion.py b/parcels/application_kernels/advectiondiffusion.py
index afae2b5188..6ad4ad2588 100644
--- a/parcels/application_kernels/advectiondiffusion.py
+++ b/parcels/application_kernels/advectiondiffusion.py
@@ -8,7 +8,7 @@
 __all__ = ["AdvectionDiffusionEM", "AdvectionDiffusionM1", "DiffusionUniformKh"]
 
 
-def AdvectionDiffusionM1(particle, fieldset, time):  # pragma: no cover
+def AdvectionDiffusionM1(particles, fieldset):  # pragma: no cover
     """Kernel for 2D advection-diffusion, solved using the Milstein scheme at first order (M1).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
@@ -23,30 +23,34 @@ def AdvectionDiffusionM1(particle, fieldset, time):  # pragma: no cover
     The Wiener increment `dW` is normally distributed with zero
     mean and a standard deviation of sqrt(dt).
     """
-    dt = particle.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    dt = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
     # Wiener increment with zero mean and std of sqrt(dt)
     dWx = np.random.normal(0, np.sqrt(np.fabs(dt)))
     dWy = np.random.normal(0, np.sqrt(np.fabs(dt)))
 
-    Kxp1 = fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon + fieldset.dres, particle]
-    Kxm1 = fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon - fieldset.dres, particle]
+    Kxp1 = fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon + fieldset.dres, particles]
+    Kxm1 = fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon - fieldset.dres, particles]
     dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres)
 
-    u, v = fieldset.UV[particle.time, particle.depth, particle.lat, particle.lon, particle]
-    bx = np.sqrt(2 * fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon, particle])
+    u, v = fieldset.UV[particles.time, particles.depth, particles.lat, particles.lon, particles]
+    bx = np.sqrt(2 * fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon, particles])
 
-    Kyp1 = fieldset.Kh_meridional[particle.time, particle.depth, particle.lat + fieldset.dres, particle.lon, particle]
-    Kym1 = fieldset.Kh_meridional[particle.time, particle.depth, particle.lat - fieldset.dres, particle.lon, particle]
+    Kyp1 = fieldset.Kh_meridional[
+        particles.time, particles.depth, particles.lat + fieldset.dres, particles.lon, particles
+    ]
+    Kym1 = fieldset.Kh_meridional[
+        particles.time, particles.depth, particles.lat - fieldset.dres, particles.lon, particles
+    ]
     dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres)
 
-    by = np.sqrt(2 * fieldset.Kh_meridional[particle.time, particle.depth, particle.lat, particle.lon, particle])
+    by = np.sqrt(2 * fieldset.Kh_meridional[particles.time, particles.depth, particles.lat, particles.lon, particles])
 
     # Particle positions are updated only after evaluating all terms.
-    particle.dlon += u * dt + 0.5 * dKdx * (dWx**2 + dt) + bx * dWx
-    particle.dlat += v * dt + 0.5 * dKdy * (dWy**2 + dt) + by * dWy
+    particles.dlon += u * dt + 0.5 * dKdx * (dWx**2 + dt) + bx * dWx
+    particles.dlat += v * dt + 0.5 * dKdy * (dWy**2 + dt) + by * dWy
 
 
-def AdvectionDiffusionEM(particle, fieldset, time):  # pragma: no cover
+def AdvectionDiffusionEM(particles, fieldset):  # pragma: no cover
     """Kernel for 2D advection-diffusion, solved using the Euler-Maruyama scheme (EM).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
@@ -59,31 +63,35 @@ def AdvectionDiffusionEM(particle, fieldset, time):  # pragma: no cover
     The Wiener increment `dW` is normally distributed with zero
     mean and a standard deviation of sqrt(dt).
     """
-    dt = particle.dt / np.timedelta64(1, "s")
+    dt = particles.dt / np.timedelta64(1, "s")
     # Wiener increment with zero mean and std of sqrt(dt)
     dWx = np.random.normal(0, np.sqrt(np.fabs(dt)))
     dWy = np.random.normal(0, np.sqrt(np.fabs(dt)))
 
-    u, v = fieldset.UV[particle.time, particle.depth, particle.lat, particle.lon, particle]
+    u, v = fieldset.UV[particles.time, particles.depth, particles.lat, particles.lon, particles]
 
-    Kxp1 = fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon + fieldset.dres, particle]
-    Kxm1 = fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon - fieldset.dres, particle]
+    Kxp1 = fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon + fieldset.dres, particles]
+    Kxm1 = fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon - fieldset.dres, particles]
     dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres)
     ax = u + dKdx
-    bx = np.sqrt(2 * fieldset.Kh_zonal[particle.time, particle.depth, particle.lat, particle.lon, particle])
-
-    Kyp1 = fieldset.Kh_meridional[particle.time, particle.depth, particle.lat + fieldset.dres, particle.lon, particle]
-    Kym1 = fieldset.Kh_meridional[particle.time, particle.depth, particle.lat - fieldset.dres, particle.lon, particle]
+    bx = np.sqrt(2 * fieldset.Kh_zonal[particles.time, particles.depth, particles.lat, particles.lon, particles])
+
+    Kyp1 = fieldset.Kh_meridional[
+        particles.time, particles.depth, particles.lat + fieldset.dres, particles.lon, particles
+    ]
+    Kym1 = fieldset.Kh_meridional[
+        particles.time, particles.depth, particles.lat - fieldset.dres, particles.lon, particles
+    ]
     dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres)
     ay = v + dKdy
-    by = np.sqrt(2 * fieldset.Kh_meridional[particle.time, particle.depth, particle.lat, particle.lon, particle])
+    by = np.sqrt(2 * fieldset.Kh_meridional[particles.time, particles.depth, particles.lat, particles.lon, particles])
 
     # Particle positions are updated only after evaluating all terms.
-    particle.dlon += ax * dt + bx * dWx
-    particle.dlat += ay * dt + by * dWy
+    particles.dlon += ax * dt + bx * dWx
+    particles.dlat += ay * dt + by * dWy
 
 
-def DiffusionUniformKh(particle, fieldset, time):  # pragma: no cover
+def DiffusionUniformKh(particles, fieldset):  # pragma: no cover
     """Kernel for simple 2D diffusion where diffusivity (Kh) is assumed uniform.
 
     Assumes that fieldset has constant fields `Kh_zonal` and `Kh_meridional`.
@@ -101,15 +109,13 @@ def DiffusionUniformKh(particle, fieldset, time):  # pragma: no cover
     The Wiener increment `dW` is normally distributed with zero
     mean and a standard deviation of sqrt(dt).
     """
-    dt = particle.dt / np.timedelta64(1, "s")
+    dt = particles.dt / np.timedelta64(1, "s")
     # Wiener increment with zero mean and std of sqrt(dt)
     dWx = np.random.normal(0, np.sqrt(np.fabs(dt)))
     dWy = np.random.normal(0, np.sqrt(np.fabs(dt)))
 
-    print(particle)
-
-    bx = np.sqrt(2 * fieldset.Kh_zonal[particle])
-    by = np.sqrt(2 * fieldset.Kh_meridional[particle])
+    bx = np.sqrt(2 * fieldset.Kh_zonal[particles])
+    by = np.sqrt(2 * fieldset.Kh_meridional[particles])
 
-    particle.dlon += bx * dWx
-    particle.dlat += by * dWy
+    particles.dlon += bx * dWx
+    particles.dlat += by * dWy
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 612374fbdf..9b3bc9978f 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -115,22 +115,22 @@ def remove_deleted(self, pset):
 
     def add_positionupdate_kernels(self):
         # Adding kernels that set and update the coordinate changes
-        def Setcoords(particle, fieldset, time):  # pragma: no cover
+        def Setcoords(particles, fieldset):  # pragma: no cover
             import numpy as np  # noqa
 
-            particle.dlon = 0
-            particle.dlat = 0
-            particle.ddepth = 0
-            particle.lon = particle.lon_nextloop
-            particle.lat = particle.lat_nextloop
-            particle.depth = particle.depth_nextloop
-            particle.time = particle.time_nextloop
+            particles.dlon = 0
+            particles.dlat = 0
+            particles.ddepth = 0
+            particles.lon = particles.lon_nextloop
+            particles.lat = particles.lat_nextloop
+            particles.depth = particles.depth_nextloop
+            particles.time = particles.time_nextloop
 
-        def Updatecoords(particle, fieldset, time):  # pragma: no cover
-            particle.lon_nextloop = particle.lon + particle.dlon
-            particle.lat_nextloop = particle.lat + particle.dlat
-            particle.depth_nextloop = particle.depth + particle.ddepth
-            particle.time_nextloop = particle.time + particle.dt
+        def Updatecoords(particles, fieldset):  # pragma: no cover
+            particles.lon_nextloop = particles.lon + particles.dlon
+            particles.lat_nextloop = particles.lat + particles.dlat
+            particles.depth_nextloop = particles.depth + particles.ddepth
+            particles.time_nextloop = particles.time + particles.dt
 
         self._pyfuncs = (Setcoords + self + Updatecoords)._pyfuncs
 
@@ -258,12 +258,12 @@ def execute(self, pset, endtime, dt):
             # run kernels for all particles that need to be evaluated
             evaluate_particles = (pset.state == StatusCode.Evaluate) & (pset.dt != 0)
             for f in self._pyfuncs:
-                f(pset[evaluate_particles], self._fieldset, None)
+                f(pset[evaluate_particles], self._fieldset)
 
                 # check for particles that have to be repeated
                 repeat_particles = pset.state == StatusCode.Repeat
                 while np.any(repeat_particles):
-                    f(pset[repeat_particles], self._fieldset, None)
+                    f(pset[repeat_particles], self._fieldset)
                     repeat_particles = pset.state == StatusCode.Repeat
 
             # revert to original dt (unless in RK45 mode)
diff --git a/tests/common_kernels.py b/tests/common_kernels.py
index 3ff946ac02..58f674cf82 100644
--- a/tests/common_kernels.py
+++ b/tests/common_kernels.py
@@ -1,18 +1,21 @@
 """Shared kernels between tests."""
 
+import numpy as np
 
-def DoNothing(particle, fieldset, time):  # pragma: no cover
+from parcels.tools.statuscodes import StatusCode
+
+
+def DoNothing(particles, fieldset):  # pragma: no cover
     pass
 
 
-def DeleteParticle(particle, fieldset, time):  # pragma: no cover
-    if particle.state >= 50:  # This captures all Errors
-        particle.delete()
+def DeleteParticle(particles, fieldset):  # pragma: no cover
+    particles.state = np.where(particles.state >= 50, StatusCode.Delete, particles.state)
 
 
-def MoveEast(particle, fieldset, time):  # pragma: no cover
-    particle.dlon += 0.1
+def MoveEast(particles, fieldset):  # pragma: no cover
+    particles.dlon += 0.1
 
 
-def MoveNorth(particle, fieldset, time):  # pragma: no cover
-    particle.dlat += 0.1
+def MoveNorth(particles, fieldset):  # pragma: no cover
+    particles.dlat += 0.1
diff --git a/tests/v4/test_advection.py b/tests/v4/test_advection.py
index ef3c8cb3a6..5308fa7ed6 100644
--- a/tests/v4/test_advection.py
+++ b/tests/v4/test_advection.py
@@ -74,9 +74,9 @@ def test_advection_zonal_with_particlefile(tmp_store):
     np.testing.assert_allclose(ds.isel(obs=-1).lon.values, pset.lon)
 
 
-def periodicBC(particle, fieldset, time):
-    particle.total_dlon += particle.dlon
-    particle.lon = np.fmod(particle.lon, 2)
+def periodicBC(particles, fieldset):
+    particles.total_dlon += particles.dlon
+    particles.lon = np.fmod(particles.lon, 2)
 
 
 def test_advection_zonal_periodic():
@@ -138,23 +138,25 @@ def test_advection_3D_outofbounds(direction, wErrorThroughSurface):
     UV = VectorField("UV", U, V)
     fieldset = FieldSet([U, V, W, UVW, UV])
 
-    def DeleteParticle(particle, fieldset, time):  # pragma: no cover
-        particle.state = np.where(particle.state == StatusCode.ErrorOutOfBounds, StatusCode.Delete, particle.state)
-        particle.state = np.where(particle.state == StatusCode.ErrorThroughSurface, StatusCode.Delete, particle.state)
+    def DeleteParticle(particles, fieldset):  # pragma: no cover
+        particles.state = np.where(particles.state == StatusCode.ErrorOutOfBounds, StatusCode.Delete, particles.state)
+        particles.state = np.where(
+            particles.state == StatusCode.ErrorThroughSurface, StatusCode.Delete, particles.state
+        )
 
-    def SubmergeParticle(particle, fieldset, time):  # pragma: no cover
-        if len(particle.state) == 0:
+    def SubmergeParticle(particles, fieldset):  # pragma: no cover
+        if len(particles.state) == 0:
             return
-        inds = np.argwhere(particle.state == StatusCode.ErrorThroughSurface).flatten()
+        inds = np.argwhere(particles.state == StatusCode.ErrorThroughSurface).flatten()
         if len(inds) == 0:
             return
-        dt = particle.dt / np.timedelta64(1, "s")
-        (u, v) = fieldset.UV[particle[inds]]
-        particle[inds].dlon = u * dt
-        particle[inds].dlat = v * dt
-        particle[inds].ddepth = 0.0
-        particle[inds].depth = 0
-        particle[inds].state = StatusCode.Evaluate
+        dt = particles.dt / np.timedelta64(1, "s")
+        (u, v) = fieldset.UV[particles[inds]]
+        particles[inds].dlon = u * dt
+        particles[inds].dlat = v * dt
+        particles[inds].ddepth = 0.0
+        particles[inds].depth = 0
+        particles[inds].state = StatusCode.Evaluate
 
     kernels = [AdvectionRK4_3D]
     if wErrorThroughSurface:
@@ -381,9 +383,9 @@ def test_stommelgyre_fieldset(method, rtol, grid_type):
         [Variable("p", initial=0.0, dtype=np.float32), Variable("p_start", initial=0.0, dtype=np.float32)]
     )
 
-    def UpdateP(particle, fieldset, time):  # pragma: no cover
-        particle.p = fieldset.P[particle.time, particle.depth, particle.lat, particle.lon]
-        particle.p_start = np.where(particle.time == 0, particle.p, particle.p_start)
+    def UpdateP(particles, fieldset):  # pragma: no cover
+        particles.p = fieldset.P[particles.time, particles.depth, particles.lat, particles.lon]
+        particles.p_start = np.where(particles.time == 0, particles.p, particles.p_start)
 
     pset = ParticleSet(fieldset, pclass=SampleParticle, lon=start_lon, lat=start_lat, time=np.timedelta64(0, "s"))
     pset.execute([kernel[method], UpdateP], dt=dt, runtime=runtime)
@@ -420,9 +422,9 @@ def test_peninsula_fieldset(method, rtol, grid_type):
         [Variable("p", initial=0.0, dtype=np.float32), Variable("p_start", initial=0.0, dtype=np.float32)]
     )
 
-    def UpdateP(particle, fieldset, time):  # pragma: no cover
-        particle.p = fieldset.P[particle.time, particle.depth, particle.lat, particle.lon]
-        particle.p_start = np.where(particle.time == 0, particle.p, particle.p_start)
+    def UpdateP(particles, fieldset):  # pragma: no cover
+        particles.p = fieldset.P[particles.time, particles.depth, particles.lat, particles.lon]
+        particles.p_start = np.where(particles.time == 0, particles.p, particles.p_start)
 
     pset = ParticleSet(fieldset, pclass=SampleParticle, lon=start_lon, lat=start_lat, time=np.timedelta64(0, "s"))
     pset.execute([kernel[method], UpdateP], dt=dt, runtime=runtime)
@@ -462,8 +464,8 @@ def test_nemo_curvilinear_fieldset():
     latp = np.linspace(-70, 88, npart)
     runtime = np.timedelta64(12, "h")  # TODO increase to 160 days
 
-    def periodicBC(particle, fieldset, time):  # pragma: no cover
-        particle.dlon = np.where(particle.lon > 180, particle.dlon - 360, particle.dlon)
+    def periodicBC(particles, fieldset):  # pragma: no cover
+        particles.dlon = np.where(particles.lon > 180, particles.dlon - 360, particles.dlon)
 
     pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)
     pset.execute([AdvectionEE, periodicBC], runtime=runtime, dt=np.timedelta64(6, "h"))
diff --git a/tests/v4/test_diffusion.py b/tests/v4/test_diffusion.py
index 1a47693023..65c9aad008 100644
--- a/tests/v4/test_diffusion.py
+++ b/tests/v4/test_diffusion.py
@@ -102,9 +102,9 @@ def test_randomexponential(lambd):
 
     pset = ParticleSet(fieldset=fieldset, lon=np.zeros(npart), lat=np.zeros(npart), depth=np.zeros(npart))
 
-    def vertical_randomexponential(particle, fieldset, time):  # pragma: no cover
+    def vertical_randomexponential(particles, fieldset):  # pragma: no cover
         # Kernel for random exponential variable in depth direction
-        particle.depth = np.random.exponential(scale=1 / fieldset.lambd, size=len(particle))
+        particles.depth = np.random.exponential(scale=1 / fieldset.lambd, size=len(particles))
 
     pset.execute(vertical_randomexponential, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"))
 
@@ -130,8 +130,8 @@ def test_randomvonmises(mu, kappa):
         fieldset=fieldset, pclass=AngleParticle, lon=np.zeros(npart), lat=np.zeros(npart), depth=np.zeros(npart)
     )
 
-    def vonmises(particle, fieldset, time):  # pragma: no cover
-        particle.angle = np.array([random.vonmisesvariate(fieldset.mu, fieldset.kappa) for _ in range(len(particle))])
+    def vonmises(particles, fieldset):  # pragma: no cover
+        particles.angle = np.array([random.vonmisesvariate(fieldset.mu, fieldset.kappa) for _ in range(len(particles))])
 
     pset.execute(vonmises, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"))
 
diff --git a/tests/v4/test_kernel.py b/tests/v4/test_kernel.py
index acceb8dff0..c2da6259ea 100644
--- a/tests/v4/test_kernel.py
+++ b/tests/v4/test_kernel.py
@@ -26,8 +26,8 @@ def fieldset() -> FieldSet:
 def test_unknown_var_in_kernel(fieldset):
     pset = ParticleSet(fieldset, lon=[0.5], lat=[0.5])
 
-    def ErrorKernel(particle, fieldset, time):  # pragma: no cover
-        particle.unknown_varname += 0.2
+    def ErrorKernel(particles, fieldset):  # pragma: no cover
+        particles.unknown_varname += 0.2
 
     with pytest.raises(KeyError, match="'unknown_varname'"):
         pset.execute(ErrorKernel, runtime=np.timedelta64(2, "s"))
diff --git a/tests/v4/test_particlefile.py b/tests/v4/test_particlefile.py
index 567793b4e6..5937f59f7f 100755
--- a/tests/v4/test_particlefile.py
+++ b/tests/v4/test_particlefile.py
@@ -118,8 +118,8 @@ def test_pfile_array_remove_all_particles(fieldset, chunks_obs, tmp_zarrfile):
 
 @pytest.mark.skip(reason="TODO v4: stuck in infinite loop")
 def test_variable_write_double(fieldset, tmp_zarrfile):
-    def Update_lon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += 0.1
+    def Update_lon(particles, fieldset):  # pragma: no cover
+        particles.dlon += 0.1
 
     particle = get_default_particle(np.float64)
     pset = ParticleSet(fieldset, pclass=particle, lon=[0], lat=[0])
@@ -198,12 +198,12 @@ def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_zarrfile, d
     )
     pfile = ParticleFile(tmp_zarrfile, outputdt=abs(dt), chunks=(1, 1))
 
-    def IncrLon(particle, fieldset, time):  # pragma: no cover
-        particle.sample_var += 1.0
-        particle.state = np.where(
-            particle.sample_var > fieldset.maxvar,
+    def IncrLon(particles, fieldset):  # pragma: no cover
+        particles.sample_var += 1.0
+        particles.state = np.where(
+            particles.sample_var > fieldset.maxvar,
             StatusCode.Delete,
-            particle.state,
+            particles.state,
         )
 
     for _ in range(npart):
@@ -224,9 +224,9 @@ def IncrLon(particle, fieldset, time):  # pragma: no cover
 
 @pytest.mark.xfail(reason="need to fix backwards in time")
 def test_write_timebackward(fieldset, tmp_zarrfile):
-    def Update_lon(particle, fieldset, time):  # pragma: no cover
-        dt = particle.dt / np.timedelta64(1, "s")
-        particle.dlon -= 0.1 * dt
+    def Update_lon(particles, fieldset):  # pragma: no cover
+        dt = particles.dt / np.timedelta64(1, "s")
+        particles.dlon -= 0.1 * dt
 
     pset = ParticleSet(
         fieldset,
@@ -259,19 +259,19 @@ def test_write_xiyi(fieldset, tmp_zarrfile):
         ]
     )
 
-    def Get_XiYi(particle, fieldset, time):  # pragma: no cover
+    def Get_XiYi(particles, fieldset):  # pragma: no cover
         """Kernel to sample the grid indices of the particle.
         Note that this sampling should be done _before_ the advection kernel
         and that the first outputted value is zero.
         Be careful when using multiple grids, as the index may be different for the grids.
         """
-        particle.pxi0 = fieldset.U.unravel_index(particle.ei)[2]
-        particle.pxi1 = fieldset.P.unravel_index(particle.ei)[2]
-        particle.pyi = fieldset.U.unravel_index(particle.ei)[1]
+        particles.pxi0 = fieldset.U.unravel_index(particles.ei)[2]
+        particles.pxi1 = fieldset.P.unravel_index(particles.ei)[2]
+        particles.pyi = fieldset.U.unravel_index(particles.ei)[1]
 
-    def SampleP(particle, fieldset, time):  # pragma: no cover
-        if time > 5 * 3600:
-            _ = fieldset.P[particle]  # To trigger sampling of the P field
+    def SampleP(particles, fieldset):  # pragma: no cover
+        if np.any(particles.time > 5 * 3600):
+            _ = fieldset.P[particles]  # To trigger sampling of the P field
 
     pset = ParticleSet(fieldset, pclass=XiYiParticle, lon=[0, 0.2], lat=[0.2, 1])
     pfile = ParticleFile(tmp_zarrfile, outputdt=dt)
@@ -302,8 +302,8 @@ def test_reset_dt(fieldset, tmp_zarrfile):
     # Assert that p.dt gets reset when a write_time is not a multiple of dt
     # for p.dt=0.02 to reach outputdt=0.05 and endtime=0.1, the steps should be [0.2, 0.2, 0.1, 0.2, 0.2, 0.1], resulting in 6 kernel executions
 
-    def Update_lon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += 0.1
+    def Update_lon(particles, fieldset):  # pragma: no cover
+        particles.dlon += 0.1
 
     particle = get_default_particle(np.float64)
     pset = ParticleSet(fieldset, pclass=particle, lon=[0], lat=[0])
@@ -319,8 +319,8 @@ def Update_lon(particle, fieldset, time):  # pragma: no cover
 def test_correct_misaligned_outputdt_dt(fieldset, tmp_zarrfile):
     """Testing that outputdt does not need to be a multiple of dt."""
 
-    def Update_lon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += particle.dt / np.timedelta64(1, "s")
+    def Update_lon(particles, fieldset):  # pragma: no cover
+        particles.dlon += particles.dt / np.timedelta64(1, "s")
 
     particle = get_default_particle(np.float64)
     pset = ParticleSet(fieldset, pclass=particle, lon=[0], lat=[0])
@@ -470,8 +470,8 @@ def test_pfile_set_towrite_False(fieldset, tmp_zarrfile):
     pset.set_variable_write_status("lat", False)
     pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1)
 
-    def Update_lon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += 0.1
+    def Update_lon(particles, fieldset):  # pragma: no cover
+        particles.dlon += 0.1
 
     pset.execute(Update_lon, runtime=10, output_file=pfile)
 
diff --git a/tests/v4/test_particleset.py b/tests/v4/test_particleset.py
index f3f2c48fad..e7a29e6ca4 100644
--- a/tests/v4/test_particleset.py
+++ b/tests/v4/test_particleset.py
@@ -119,8 +119,8 @@ def test_pset_starttime_not_multiple_dt(fieldset):
     datetimes = [fieldset.time_interval.left + np.timedelta64(t, "s") for t in times]
     pset = ParticleSet(fieldset, lon=[0] * len(times), lat=[0] * len(times), pclass=Particle, time=datetimes)
 
-    def Addlon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += particle.dt / np.timedelta64(1, "s")
+    def Addlon(particles, fieldset):  # pragma: no cover
+        particles.dlon += particles.dt / np.timedelta64(1, "s")
 
     pset.execute(Addlon, dt=np.timedelta64(2, "s"), runtime=np.timedelta64(8, "s"), verbose_progress=False)
     assert np.allclose([p.lon_nextloop for p in pset], [8 - t for t in times])
diff --git a/tests/v4/test_particleset_execute.py b/tests/v4/test_particleset_execute.py
index 80087d1207..c5aef6905c 100644
--- a/tests/v4/test_particleset_execute.py
+++ b/tests/v4/test_particleset_execute.py
@@ -154,8 +154,8 @@ def test_pset_remove_particle_in_kernel(fieldset):
     npart = 100
     pset = ParticleSet(fieldset, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart))
 
-    def DeleteKernel(particle, fieldset, time):  # pragma: no cover
-        particle.state = np.where((particle.lon >= 0.4) & (particle.lon <= 0.6), StatusCode.Delete, particle.state)
+    def DeleteKernel(particles, fieldset):  # pragma: no cover
+        particles.state = np.where((particles.lon >= 0.4) & (particles.lon <= 0.6), StatusCode.Delete, particles.state)
 
     pset.execute(pset.Kernel(DeleteKernel), runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"))
     indices = [i for i in range(npart) if not (40 <= i < 60)]
@@ -169,8 +169,10 @@ def DeleteKernel(particle, fieldset, time):  # pragma: no cover
 def test_pset_stop_simulation(fieldset, npart):
     pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), pclass=Particle)
 
-    def Delete(particle, fieldset, time):  # pragma: no cover
-        particle[particle.time >= fieldset.time_interval.left + np.timedelta64(4, "s")].state = StatusCode.StopExecution
+    def Delete(particles, fieldset):  # pragma: no cover
+        particles[
+            particles.time >= fieldset.time_interval.left + np.timedelta64(4, "s")
+        ].state = StatusCode.StopExecution
 
     pset.execute(Delete, dt=np.timedelta64(1, "s"), runtime=np.timedelta64(21, "s"))
     assert pset[0].time == fieldset.time_interval.left + np.timedelta64(4, "s")
@@ -180,8 +182,8 @@ def Delete(particle, fieldset, time):  # pragma: no cover
 def test_pset_multi_execute(fieldset, with_delete, npart=10, n=5):
     pset = ParticleSet(fieldset, lon=np.linspace(0, 1, npart), lat=np.zeros(npart))
 
-    def AddLat(particle, fieldset, time):  # pragma: no cover
-        particle.dlat += 0.1
+    def AddLat(particles, fieldset):  # pragma: no cover
+        particles.dlat += 0.1
 
     k_add = pset.Kernel(AddLat)
     for _ in range(n + 1):
@@ -213,8 +215,8 @@ def test_dont_run_particles_outside_starttime(fieldset):
     start_times = [fieldset.time_interval.left + np.timedelta64(t, "s") for t in [0, 2, 10]]
     endtime = fieldset.time_interval.left + np.timedelta64(8, "s")
 
-    def AddLon(particle, fieldset, time):  # pragma: no cover
-        particle.lon += 1
+    def AddLon(particles, fieldset):  # pragma: no cover
+        particles.lon += 1
 
     pset = ParticleSet(fieldset, lon=np.zeros(len(start_times)), lat=np.zeros(len(start_times)), time=start_times)
     pset.execute(AddLon, dt=np.timedelta64(1, "s"), endtime=endtime)
@@ -245,12 +247,12 @@ def test_some_particles_throw_outofbounds(zonal_flow_fieldset):
 
 
 def test_delete_on_all_errors(fieldset):
-    def MoveRight(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += 1
-        fieldset.U[particle.time, particle.depth, particle.lat, particle.lon, particle]
+    def MoveRight(particles, fieldset):  # pragma: no cover
+        particles.dlon += 1
+        fieldset.U[particles.time, particles.depth, particles.lat, particles.lon, particles]
 
-    def DeleteAllErrorParticles(particle, fieldset, time):  # pragma: no cover
-        particle[particle.state > 20].state = StatusCode.Delete
+    def DeleteAllErrorParticles(particles, fieldset):  # pragma: no cover
+        particles[particles.state > 20].state = StatusCode.Delete
 
     pset = ParticleSet(fieldset, lon=[1e5, 2], lat=[0, 0])
     pset.execute([MoveRight, DeleteAllErrorParticles], runtime=np.timedelta64(10, "s"), dt=np.timedelta64(1, "s"))
@@ -261,8 +263,8 @@ def test_some_particles_throw_outoftime(fieldset):
     time = [fieldset.time_interval.left + np.timedelta64(t, "D") for t in [0, 350]]
     pset = ParticleSet(fieldset, lon=np.zeros_like(time), lat=np.zeros_like(time), time=time)
 
-    def FieldAccessOutsideTime(particle, fieldset, time):  # pragma: no cover
-        fieldset.U[particle.time + np.timedelta64(400, "D"), particle.depth, particle.lat, particle.lon, particle]
+    def FieldAccessOutsideTime(particles, fieldset):  # pragma: no cover
+        fieldset.U[particles.time + np.timedelta64(400, "D"), particles.depth, particles.lat, particles.lon, particles]
 
     with pytest.raises(TimeExtrapolationError):
         pset.execute(FieldAccessOutsideTime, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(10, "D"))
@@ -278,8 +280,8 @@ def test_errorinterpolation(fieldset):
     def NaNInterpolator(field, ti, position, tau, t, z, y, x):  # pragma: no cover
         return np.nan * np.zeros_like(x)
 
-    def SampleU(particle, fieldset, time):  # pragma: no cover
-        fieldset.U[particle.time, particle.depth, particle.lat, particle.lon, particle]
+    def SampleU(particles, fieldset):  # pragma: no cover
+        fieldset.U[particles.time, particles.depth, particles.lat, particles.lon, particles]
 
     fieldset.U.interp_method = NaNInterpolator
     pset = ParticleSet(fieldset, lon=[0, 2], lat=[0, 0])
@@ -288,9 +290,9 @@ def SampleU(particle, fieldset, time):  # pragma: no cover
 
 
 def test_execution_check_stopallexecution(fieldset):
-    def addoneLon(particle, fieldset, time):  # pragma: no cover
-        particle.dlon += 1
-        particle[particle.lon + particle.dlon >= 10].state = StatusCode.StopAllExecution
+    def addoneLon(particles, fieldset):  # pragma: no cover
+        particles.dlon += 1
+        particles[particles.lon + particles.dlon >= 10].state = StatusCode.StopAllExecution
 
     pset = ParticleSet(fieldset, lon=[0, 0], lat=[0, 0])
     pset.execute(addoneLon, runtime=np.timedelta64(20, "s"), dt=np.timedelta64(1, "s"))
@@ -301,15 +303,15 @@ def addoneLon(particle, fieldset, time):  # pragma: no cover
 def test_execution_recover_out_of_bounds(fieldset):
     npart = 2
 
-    def MoveRight(particle, fieldset, time):  # pragma: no cover
-        fieldset.U[particle.time, particle.depth, particle.lat, particle.lon + 0.1, particle]
-        particle.dlon += 0.1
+    def MoveRight(particles, fieldset):  # pragma: no cover
+        fieldset.U[particles.time, particles.depth, particles.lat, particles.lon + 0.1, particles]
+        particles.dlon += 0.1
 
-    def MoveLeft(particle, fieldset, time):  # pragma: no cover
-        inds = np.where(particle.state == StatusCode.ErrorOutOfBounds)
-        print(inds, particle.state)
-        particle[inds].dlon -= 1.0
-        particle[inds].state = StatusCode.Success
+    def MoveLeft(particles, fieldset):  # pragma: no cover
+        inds = np.where(particles.state == StatusCode.ErrorOutOfBounds)
+        print(inds, particles.state)
+        particles[inds].dlon -= 1.0
+        particles[inds].state = StatusCode.Success
 
     lon = np.linspace(0.05, 6.95, npart)
     lat = np.linspace(1, 0, npart)
@@ -337,8 +339,8 @@ def test_execution_runtime(fieldset, starttime, runtime, dt, npart):
 
 
 def test_changing_dt_in_kernel(fieldset):
-    def KernelCounter(particle, fieldset, time):  # pragma: no cover
-        particle.lon += 1
+    def KernelCounter(particles, fieldset):  # pragma: no cover
+        particles.lon += 1
 
     pset = ParticleSet(fieldset, lon=np.zeros(1), lat=np.zeros(1))
     pset.execute(KernelCounter, dt=np.timedelta64(2, "s"), runtime=np.timedelta64(5, "s"))
@@ -351,8 +353,8 @@ def KernelCounter(particle, fieldset, time):  # pragma: no cover
 def test_execution_fail_python_exception(fieldset, npart):
     pset = ParticleSet(fieldset, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart))
 
-    def PythonFail(particle, fieldset, time):  # pragma: no cover
-        inds = np.argwhere(particle.time >= fieldset.time_interval.left + np.timedelta64(10, "s"))
+    def PythonFail(particles, fieldset):  # pragma: no cover
+        inds = np.argwhere(particles.time >= fieldset.time_interval.left + np.timedelta64(10, "s"))
         if inds.size > 0:
             raise RuntimeError("Enough is enough!")
 

From f411b45468d6281a2bfc866344b7f81ec087b06a Mon Sep 17 00:00:00 2001
From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com>
Date: Fri, 12 Sep 2025 10:26:00 +0200
Subject: [PATCH 3/3] Add test test_kernel_signature

---
 parcels/tools/_helpers.py |  8 +++++---
 tests/v4/test_kernel.py   | 40 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+), 3 deletions(-)

diff --git a/parcels/tools/_helpers.py b/parcels/tools/_helpers.py
index 22b1c0c143..7432572767 100644
--- a/parcels/tools/_helpers.py
+++ b/parcels/tools/_helpers.py
@@ -88,10 +88,12 @@ def _assert_same_function_signature(f: Callable, *, ref: Callable, context: str)
             f"{context} function must have {len(sig_ref.parameters)} parameters, got {len(sig.parameters)}"
         )
 
-    for (_name1, param1), (_name2, param2) in zip(sig_ref.parameters.items(), sig.parameters.items(), strict=False):
+    for param1, param2 in zip(sig_ref.parameters.values(), sig.parameters.values(), strict=False):
         if param1.kind != param2.kind:
             raise ValueError(
-                f"Parameter '{_name2}' has incorrect parameter kind. Expected {param1.kind}, got {param2.kind}"
+                f"Parameter '{param2.name}' has incorrect parameter kind. Expected {param1.kind}, got {param2.kind}"
             )
         if param1.name != param2.name:
-            raise ValueError(f"Parameter '{_name2}' has incorrect name. Expected '{param1.name}', got '{param2.name}'")
+            raise ValueError(
+                f"Parameter '{param2.name}' has incorrect name. Expected '{param1.name}', got '{param2.name}'"
+            )
diff --git a/tests/v4/test_kernel.py b/tests/v4/test_kernel.py
index c2da6259ea..3ece81a744 100644
--- a/tests/v4/test_kernel.py
+++ b/tests/v4/test_kernel.py
@@ -85,3 +85,43 @@ def test_kernel_from_list_error_checking(fieldset):
     with pytest.raises(ValueError, match="Argument `pyfunc_list` should be a list of functions."):
         kernels_mixed = pset.Kernel([pset.Kernel(AdvectionRK4), MoveEast, MoveNorth])
         assert kernels_mixed.funcname == "AdvectionRK4MoveEastMoveNorth"
+
+
+def test_kernel_signature(fieldset):
+    pset = ParticleSet(fieldset, lon=[0.5], lat=[0.5])
+
+    def good_kernel(particles, fieldset):
+        pass
+
+    def version_3_kernel(particle, fieldset, time):
+        pass
+
+    def version_3_kernel_without_time(particle, fieldset):
+        pass
+
+    def kernel_switched_args(fieldset, particle):
+        pass
+
+    def kernel_with_forced_kwarg(particles, *, fieldset=0):
+        pass
+
+    pset.Kernel(good_kernel)
+
+    with pytest.raises(ValueError, match="Kernel function must have 2 parameters, got 3"):
+        pset.Kernel(version_3_kernel)
+
+    with pytest.raises(
+        ValueError, match="Parameter 'particle' has incorrect name. Expected 'particles', got 'particle'"
+    ):
+        pset.Kernel(version_3_kernel_without_time)
+
+    with pytest.raises(
+        ValueError, match="Parameter 'fieldset' has incorrect name. Expected 'particles', got 'fieldset'"
+    ):
+        pset.Kernel(kernel_switched_args)
+
+    with pytest.raises(
+        ValueError,
+        match="Parameter 'fieldset' has incorrect parameter kind. Expected POSITIONAL_OR_KEYWORD, got KEYWORD_ONLY",
+    ):
+        pset.Kernel(kernel_with_forced_kwarg)