From 7ba5168ded19e2a08f8048082c5de74417b95f72 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Mon, 15 Jun 2026 14:20:27 +0200
Subject: [PATCH 01/11] Moving RK2 and RK4 methods to separate integrators.py
 functions

---
 src/parcels/kernels/_advection.py   | 71 ++++++++++++++---------------
 src/parcels/kernels/_integration.py | 51 +++++++++++++++++++++
 2 files changed, 84 insertions(+), 38 deletions(-)
 create mode 100644 src/parcels/kernels/_integration.py

diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py
index 792c8ee17f..639b6f5890 100644
--- a/src/parcels/kernels/_advection.py
+++ b/src/parcels/kernels/_advection.py
@@ -5,6 +5,7 @@
 import numpy as np
 
 from parcels._core.statuscodes import StatusCode
+from parcels.kernels._integration import RK2, RK4
 
 __all__ = [
     "AdvectionAnalytical",
@@ -34,59 +35,53 @@ def _constrain_dt_to_within_time_interval(time_interval, time, dt):
 def AdvectionRK2(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
-    (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 * dt, particles.z, lat1, lon1, particles]
-    particles.dlon += u2 * dt
-    particles.dlat += v2 * dt
+
+    def rhs(fieldset, time, depth, lat, lon, particles):
+        return fieldset.UV[time, depth, lat, lon, particles]
+
+    u, v = RK2(fieldset, particles, rhs)
+
+    particles.dlon += u * dt
+    particles.dlat += v * dt
 
 
 def AdvectionRK2_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
-    (u1, v1, w1) = fieldset.UVW[particles]
-    lon1 = particles.lon + u1 * 0.5 * dt
-    lat1 = particles.lat + v1 * 0.5 * dt
-    z1 = particles.z + w1 * 0.5 * dt
-    (u2, v2, w2) = fieldset.UVW[particles.time + 0.5 * dt, z1, lat1, lon1, particles]
-    particles.dlon += u2 * dt
-    particles.dlat += v2 * dt
-    particles.dz += w2 * dt
+
+    def rhs(fieldset, time, depth, lat, lon, particles):
+        return fieldset.UVW[time, depth, lat, lon, particles]
+
+    u, v, w = RK2(fieldset, particles, rhs)
+
+    particles.dlon += u * dt
+    particles.dlat += v * dt
+    particles.dz += w * dt
 
 
 def AdvectionRK4(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
-    (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 * dt, particles.z, 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 * dt, particles.z, lat2, lon2, particles]
-    lon3, lat3 = (particles.lon + u3 * dt, particles.lat + v3 * dt)
-    (u4, v4) = fieldset.UV[particles.time + dt, particles.z, 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 rhs(fieldset, time, depth, lat, lon, particles):
+        return fieldset.UV[time, depth, lat, lon, particles]
+
+    u, v = RK4(fieldset, particles, rhs)
+    particles.dlon += u * dt
+    particles.dlat += v * dt
 
 
 def AdvectionRK4_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
-    (u1, v1, w1) = fieldset.UVW[particles]
-    lon1 = particles.lon + u1 * 0.5 * dt
-    lat1 = particles.lat + v1 * 0.5 * dt
-    z1 = particles.z + w1 * 0.5 * dt
-    (u2, v2, w2) = fieldset.UVW[particles.time + 0.5 * dt, z1, lat1, lon1, particles]
-    lon2 = particles.lon + u2 * 0.5 * dt
-    lat2 = particles.lat + v2 * 0.5 * dt
-    z2 = particles.z + w2 * 0.5 * dt
-    (u3, v3, w3) = fieldset.UVW[particles.time + 0.5 * dt, z2, lat2, lon2, particles]
-    lon3 = particles.lon + u3 * dt
-    lat3 = particles.lat + v3 * dt
-    z3 = particles.z + w3 * dt
-    (u4, v4, w4) = fieldset.UVW[particles.time + dt, z3, lat3, lon3, particles]
-    particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
-    particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
-    particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt
+
+    def rhs(fieldset, time, depth, lat, lon, particles):
+        return fieldset.UVW[time, depth, lat, lon, particles]
+
+    u, v, w = RK4(fieldset, particles, rhs)
+    particles.dlon += u * dt
+    particles.dlat += v * dt
+    particles.dz += w * dt
 
 
 def AdvectionEE(particles, fieldset):  # pragma: no cover
diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
new file mode 100644
index 0000000000..228344977a
--- /dev/null
+++ b/src/parcels/kernels/_integration.py
@@ -0,0 +1,51 @@
+"""Collection of time integrators for use in Parcels Kernels"""
+
+
+def RK2(fieldset, particles, rhs):
+    z, lat, lon = particles.z, particles.lat, particles.lon
+    fields = rhs(fieldset, particles.time, z, lat, lon, particles)
+    if len(fields) == 1:
+        raise NotImplementedError("RK2 integration is not implemented for Fields that return only one component.")
+    if len(fields) > 1:
+        lon = particles.lon + fields[0] * 0.5 * particles.dt
+        lat = particles.lat + fields[1] * 0.5 * particles.dt
+    if len(fields) > 2:
+        z = particles.z + fields[2] * 0.5 * particles.dt
+    t = particles.time + 0.5 * particles.dt
+    return rhs(fieldset, t, z, lat, lon, particles)
+
+
+def RK4(fieldset, particles, rhs):
+    z, lat, lon = particles.z, particles.lat, particles.lon
+    k1 = rhs(fieldset, particles.time, z, lat, lon, particles)
+    if len(k1) == 1:
+        raise NotImplementedError("RK4 integration is not implemented for Fields that return only one component.")
+    if len(k1) > 1:
+        lon = particles.lon + k1[0] * 0.5 * particles.dt
+        lat = particles.lat + k1[1] * 0.5 * particles.dt
+    if len(k1) > 2:
+        z = particles.z + k1[2] * 0.5 * particles.dt
+    t = particles.time + 0.5 * particles.dt
+    k2 = rhs(fieldset, t, z, lat, lon, particles)
+    if len(k2) > 1:
+        lon = particles.lon + k2[0] * 0.5 * particles.dt
+        lat = particles.lat + k2[1] * 0.5 * particles.dt
+    if len(k2) > 2:
+        z = particles.z + k2[2] * 0.5 * particles.dt
+    t = particles.time + 0.5 * particles.dt
+    k3 = rhs(fieldset, t, z, lat, lon, particles)
+    if len(k3) > 1:
+        lon = particles.lon + k3[0] * particles.dt
+        lat = particles.lat + k3[1] * particles.dt
+    if len(k3) > 2:
+        z = particles.z + k3[2] * particles.dt
+    t = particles.time + particles.dt
+    k4 = rhs(fieldset, t, z, lat, lon, particles)
+    if len(k4) == 2:
+        return ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0, (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0)
+    if len(k4) == 3:
+        return (
+            (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0,
+            (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0,
+            (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.0,
+        )

From 343afd2d34f194910939524a6ede6a73bcf4bd23 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Mon, 15 Jun 2026 16:42:16 +0200
Subject: [PATCH 02/11] Moving rhs definitions to outside Kernels

For potential speed-up since no redefinition necessary
---
 src/parcels/kernels/_advection.py | 30 ++++++++++++------------------
 1 file changed, 12 insertions(+), 18 deletions(-)

diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py
index 639b6f5890..7b3164e83d 100644
--- a/src/parcels/kernels/_advection.py
+++ b/src/parcels/kernels/_advection.py
@@ -32,15 +32,19 @@ def _constrain_dt_to_within_time_interval(time_interval, time, dt):
     return dt
 
 
+def rhs_uv(fieldset, time, depth, lat, lon, particles):
+    return fieldset.UV[time, depth, lat, lon, particles]
+
+
+def rhs_uvw(fieldset, time, depth, lat, lon, particles):
+    return fieldset.UVW[time, depth, lat, lon, particles]
+
+
 def AdvectionRK2(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
 
-    def rhs(fieldset, time, depth, lat, lon, particles):
-        return fieldset.UV[time, depth, lat, lon, particles]
-
-    u, v = RK2(fieldset, particles, rhs)
-
+    u, v = RK2(fieldset, particles, rhs_uv)
     particles.dlon += u * dt
     particles.dlat += v * dt
 
@@ -49,11 +53,7 @@ def AdvectionRK2_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
 
-    def rhs(fieldset, time, depth, lat, lon, particles):
-        return fieldset.UVW[time, depth, lat, lon, particles]
-
-    u, v, w = RK2(fieldset, particles, rhs)
-
+    u, v, w = RK2(fieldset, particles, rhs_uvw)
     particles.dlon += u * dt
     particles.dlat += v * dt
     particles.dz += w * dt
@@ -63,10 +63,7 @@ def AdvectionRK4(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
 
-    def rhs(fieldset, time, depth, lat, lon, particles):
-        return fieldset.UV[time, depth, lat, lon, particles]
-
-    u, v = RK4(fieldset, particles, rhs)
+    u, v = RK4(fieldset, particles, rhs_uv)
     particles.dlon += u * dt
     particles.dlat += v * dt
 
@@ -75,10 +72,7 @@ def AdvectionRK4_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
 
-    def rhs(fieldset, time, depth, lat, lon, particles):
-        return fieldset.UVW[time, depth, lat, lon, particles]
-
-    u, v, w = RK4(fieldset, particles, rhs)
+    u, v, w = RK4(fieldset, particles, rhs_uvw)
     particles.dlon += u * dt
     particles.dlat += v * dt
     particles.dz += w * dt

From a9c64554abbd6ff2915a0786db1b3bad16fa5e73 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Mon, 15 Jun 2026 16:42:59 +0200
Subject: [PATCH 03/11] Adding possibility to ingest Copernicus wind data

Since it's on a different grid, has to be a different item in the dict
---
 src/parcels/_datasets/remote.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/parcels/_datasets/remote.py b/src/parcels/_datasets/remote.py
index eb2206591a..35253a8f1c 100644
--- a/src/parcels/_datasets/remote.py
+++ b/src/parcels/_datasets/remote.py
@@ -52,6 +52,7 @@ def _get_data_home() -> Path:
         "data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m_uo-vo_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc",
         "data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-so_anfc_0.083deg_P1D-m_so_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc",
         "data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-thetao_anfc_0.083deg_P1D-m_thetao_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc",
+        "data/CopernicusMarine_data_for_Argo_tutorial/cmems_obs-wind_glo_phy_my_l4_0.125deg_PT1H_multi-vars_31.06E-32.94E_32.94S-30.06S_2024-01-01-2024-02-01.nc",
     ]
     + [
         "data/DecayingMovingEddy_data/decaying_moving_eddyU.nc",
@@ -216,6 +217,7 @@ class _Purpose(enum.Enum):
     ("Peninsula_data/T", (_V3Dataset(_ODIE,"data/Peninsula_data/peninsulaT.nc"), _Purpose.TUTORIAL)),
     ("GlobCurrent_example_data/data", (_V3Dataset(_ODIE,"data/GlobCurrent_example_data/*000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc", pre_decode_cf_callable=patch_dataset_v4_compat), _Purpose.TUTORIAL)),
     ("CopernicusMarine_data_for_Argo_tutorial/data", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-*.nc"), _Purpose.TUTORIAL)),
+    ("CopernicusMarine_data_for_Argo_tutorial/data_wind", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_Argo_tutorial/cmems_obs-wind_glo_phy*.nc"), _Purpose.TUTORIAL)),
     ("DecayingMovingEddy_data/U", (_V3Dataset(_ODIE,"data/DecayingMovingEddy_data/decaying_moving_eddyU.nc"), _Purpose.TUTORIAL)),
     ("DecayingMovingEddy_data/V", (_V3Dataset(_ODIE,"data/DecayingMovingEddy_data/decaying_moving_eddyV.nc"), _Purpose.TUTORIAL)),
     ("FESOM_periodic_channel/fesom_channel", (_V3Dataset(_ODIE,"data/FESOM_periodic_channel/fesom_channel.nc"), _Purpose.TUTORIAL)),

From 0011abe0b3b88b12c552d292712abda31988606a Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Mon, 15 Jun 2026 16:44:09 +0200
Subject: [PATCH 04/11] Expend dt integrators tutorial with section on using
 RK2 function for custom Kernel

---
 .../examples/tutorial_dt_integrators.ipynb    | 145 ++++++++++++++++--
 1 file changed, 136 insertions(+), 9 deletions(-)

diff --git a/docs/user_guide/examples/tutorial_dt_integrators.ipynb b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
index 2acd3ebda5..44fc3deb5f 100644
--- a/docs/user_guide/examples/tutorial_dt_integrators.ipynb
+++ b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
@@ -523,7 +523,9 @@
    "outputs": [],
    "source": [
     "# And the second-order Runge-Kutta (RK2) method\n",
-    "print(inspect.getsource(parcels.kernels.AdvectionRK2))"
+    "print(inspect.getsource(parcels.kernels._advection.rhs_uv))\n",
+    "print(inspect.getsource(parcels.kernels.AdvectionRK2))\n",
+    "print(inspect.getsource(parcels.kernels._integration.RK2))"
    ]
   },
   {
@@ -534,7 +536,9 @@
    "outputs": [],
    "source": [
     "# And RK4:\n",
-    "print(inspect.getsource(parcels.kernels.AdvectionRK4))"
+    "print(inspect.getsource(parcels.kernels._advection.rhs_uv))\n",
+    "print(inspect.getsource(parcels.kernels.AdvectionRK4))\n",
+    "print(inspect.getsource(parcels.kernels._integration.RK4))"
    ]
   },
   {
@@ -795,6 +799,135 @@
    "cell_type": "markdown",
    "id": "36",
    "metadata": {},
+   "source": [
+    "## Writing custom Kernels with higher-order integrators\n",
+    "\n",
+    "Similarly to the built-in AdvectionRK2 Kernel, you can also write your own custom Kernels that use higher-order integrators. To do so, you can use the `RK2` and `RK4` integrators from `parcels.integrators` in your Kernel definition. \n",
+    "\n",
+    "The 'trick' is to write a right-hand-side function (`rhs`) that computes the change in position at a given time and location, and then use this function in the `RK2` or `RK4` integrator to compute the change in position over a given timestep.\n",
+    "\n",
+    "For example, let's say we want to write a custom Kernel that computes advection and a simple windage effect. We can write the right-hand-side function as follows"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def rhs(fieldset, time, depth, lat, lon, particles):\n",
+    "    uo, vo = fieldset.UV[time, depth, lat, lon, particles]\n",
+    "    uw, vw = fieldset.UVw[time, depth, lat, lon, particles]\n",
+    "    return (uo + fieldset.wind_coeff * (uw - uo), vo + fieldset.wind_coeff * (vw - vo))\n",
+    "\n",
+    "\n",
+    "def WindageRK2(particles, fieldset):\n",
+    "\n",
+    "    u, v = parcels.kernels._integration.RK2(fieldset, particles, rhs)\n",
+    "\n",
+    "    particles.dlon += u * particles.dt\n",
+    "    particles.dlat += v * particles.dt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38",
+   "metadata": {},
+   "source": [
+    "Now, if we also load in the wind data for the same region with the code below, we can execute this Kernel to compute the trajectories of particles that are advected by the ocean currents and affected by the wind. The solution will be second-order accurate in time, since we are using the RK2 integrator; but it is easy to swap in the RK4 integrator if you want to increase the accuracy even more."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_w = parcels.tutorial.open_dataset(\n",
+    "    \"CopernicusMarine_data_for_Argo_tutorial/data_wind\"\n",
+    ")\n",
+    "ds_w = ds_w.expand_dims({\"depth\": [0]})\n",
+    "ds_w.load()  # load the dataset into memory\n",
+    "\n",
+    "fields = {\"Uw\": ds_w[\"eastward_wind\"], \"Vw\": ds_w[\"northward_wind\"]}\n",
+    "ds_wind = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
+    "fieldset_w = parcels.FieldSet.from_sgrid_conventions(ds_wind)\n",
+    "\n",
+    "fieldset.add_field(fieldset_w.Uw)\n",
+    "fieldset.add_field(fieldset_w.Vw)\n",
+    "UVw = parcels.VectorField(\n",
+    "    \"UVw\",\n",
+    "    fieldset_w.Uw,\n",
+    "    fieldset_w.Vw,\n",
+    "    vector_interp_method=fieldset.UV.vector_interp_method,\n",
+    ")\n",
+    "fieldset.add_field(UVw)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "40",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fieldset.wind_coeff = 0.01  # windage coefficient (1% of the relative wind speed)\n",
+    "\n",
+    "pset = parcels.ParticleSet(\n",
+    "    fieldset=fieldset,\n",
+    "    pclass=parcels.Particle,\n",
+    "    time=initial_release_times,\n",
+    "    z=initial_release_zs,\n",
+    "    lat=initial_release_lats,\n",
+    "    lon=initial_release_lons,\n",
+    ")\n",
+    "\n",
+    "pfile = parcels.ParticleFile(\n",
+    "    path=OUTPUT_FOLDER / \"WindRK2.parquet\",\n",
+    "    outputdt=outputdt,\n",
+    "    mode=\"w\",\n",
+    ")\n",
+    "\n",
+    "pset.execute(\n",
+    "    WindageRK2,\n",
+    "    runtime=runtime,\n",
+    "    dt=np.timedelta64(1, \"h\"),\n",
+    "    output_file=pfile,\n",
+    "    verbose_progress=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "41",
+   "metadata": {},
+   "source": [
+    "As you can see, the particles move very differently when we include the windage effect!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "42",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = parcels.read_particlefile(pfile.path)\n",
+    "for traj in df.partition_by(\"particle_id\"):\n",
+    "    plt.plot(\n",
+    "        traj[\"lon\"],\n",
+    "        traj[\"lat\"],\n",
+    "        alpha=0.75,\n",
+    "    )\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "43",
+   "metadata": {},
    "source": [
     "## Summary\n",
     "If you want to test the accuracy and efficiency of your own simulation, here is a brief list of things to consider:\n",
@@ -802,17 +935,11 @@
     "- **Run the simulation for the full runtime.** Many errors grow over time so to quantify these, we must take into account the full length of the simulation.\n",
     "- **Consider which step introduces the largest error.** If we compute advection and diffusion, increasing the accuracy of one process much more than the other will not improve the results."
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "37",
-   "metadata": {},
-   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "docs",
+   "display_name": "Parcels:docs (3.14.4)",
    "language": "python",
    "name": "python3"
   },

From 1d273bce62ea97289853e2a5510a718b4b9f913f Mon Sep 17 00:00:00 2001
From: Erik van Sebille <erikvansebille@gmail.com>
Date: Mon, 15 Jun 2026 16:57:01 +0200
Subject: [PATCH 05/11] Merge conflicts

---
 docs/user_guide/examples/tutorial_dt_integrators.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/user_guide/examples/tutorial_dt_integrators.ipynb b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
index 44fc3deb5f..0823c597c9 100644
--- a/docs/user_guide/examples/tutorial_dt_integrators.ipynb
+++ b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
@@ -861,7 +861,7 @@
     "    \"UVw\",\n",
     "    fieldset_w.Uw,\n",
     "    fieldset_w.Vw,\n",
-    "    vector_interp_method=fieldset.UV.vector_interp_method,\n",
+    "    interp_method=fieldset.UV.interp_method,\n",
     ")\n",
     "fieldset.add_field(UVw)"
    ]

From 7afa8383da67c6ba8438bb6cb0301905ae48f851 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Thu, 8 Oct 2026 09:31:01 +0200
Subject: [PATCH 06/11] Updating RK-integrators to use t,z,y,x

---
 .../examples/tutorial_dt_integrators.ipynb    | 43 +++++++------------
 src/parcels/kernels/_advection.py             |  8 ++--
 src/parcels/kernels/_integration.py           | 40 ++++++++---------
 3 files changed, 40 insertions(+), 51 deletions(-)

diff --git a/docs/user_guide/examples/tutorial_dt_integrators.ipynb b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
index fd90387f78..1357ebba3e 100644
--- a/docs/user_guide/examples/tutorial_dt_integrators.ipynb
+++ b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
@@ -813,9 +813,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def rhs(fieldset, time, depth, lat, lon, particles):\n",
-    "    uo, vo = fieldset.UV[time, depth, lat, lon, particles]\n",
-    "    uw, vw = fieldset.UVw[time, depth, lat, lon, particles]\n",
+    "def rhs(fieldset, t, z, y, x, particles):\n",
+    "    uo, vo = fieldset.UV[t, z, y, x, particles]\n",
+    "    uw, vw = fieldset.UVw[t, z, y, x, particles]\n",
     "    return (uo + fieldset.wind_coeff * (uw - uo), vo + fieldset.wind_coeff * (vw - vo))\n",
     "\n",
     "\n",
@@ -823,8 +823,8 @@
     "\n",
     "    u, v = parcels.kernels._integration.RK2(fieldset, particles, rhs)\n",
     "\n",
-    "    particles.dlon += u * particles.dt\n",
-    "    particles.dlat += v * particles.dt"
+    "    particles.dx += u * particles.dt\n",
+    "    particles.dy += v * particles.dt"
    ]
   },
   {
@@ -850,17 +850,10 @@
     "\n",
     "fields = {\"Uw\": ds_w[\"eastward_wind\"], \"Vw\": ds_w[\"northward_wind\"]}\n",
     "ds_wind = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
-    "fieldset_w = parcels.FieldSet.from_sgrid_conventions(ds_wind)\n",
-    "\n",
-    "fieldset.add_field(fieldset_w.Uw)\n",
-    "fieldset.add_field(fieldset_w.Vw)\n",
-    "UVw = parcels.VectorField(\n",
-    "    \"UVw\",\n",
-    "    fieldset_w.Uw,\n",
-    "    fieldset_w.Vw,\n",
-    "    interp_method=fieldset.UV.interp_method,\n",
-    ")\n",
-    "fieldset.add_field(UVw)"
+    "fieldset_w = parcels.FieldSet.from_sgrid_conventions(ds_wind, vector_fields={\"UVw\": (\"Uw\", \"Vw\")})\n",
+    "\n",
+    "fieldset += fieldset_w\n",
+    "fieldset.describe()"
    ]
   },
   {
@@ -875,10 +868,10 @@
     "pset = parcels.ParticleSet(\n",
     "    fieldset=fieldset,\n",
     "    pclass=parcels.Particle,\n",
-    "    time=initial_release_times,\n",
-    "    z=initial_release_zs,\n",
-    "    lat=initial_release_lats,\n",
-    "    lon=initial_release_lons,\n",
+    "    t=initial_release_times,\n",
+    "    z=np.ones_like(initial_release_times),\n",
+    "    y=initial_release_lats,\n",
+    "    x=initial_release_lons,\n",
     ")\n",
     "\n",
     "pfile = parcels.ParticleFile(\n",
@@ -913,11 +906,7 @@
    "source": [
     "df = parcels.read_particlefile(pfile.path)\n",
     "for traj in df.partition_by(\"particle_id\"):\n",
-    "    plt.plot(\n",
-    "        traj[\"lon\"],\n",
-    "        traj[\"lat\"],\n",
-    "        alpha=0.75,\n",
-    "    )\n",
+    "    plt.plot(traj[\"x\"], traj[\"y\"], alpha=0.75)\n",
     "plt.show()"
    ]
   },
@@ -936,7 +925,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Parcels:docs (3.14.4)",
+   "display_name": "Parcels:default (3.14.7)",
    "language": "python",
    "name": "python3"
   },
@@ -950,7 +939,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.14.6"
+   "version": "3.14.7"
   }
  },
  "nbformat": 4,
diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py
index e55a0be4ca..1006a78566 100644
--- a/src/parcels/kernels/_advection.py
+++ b/src/parcels/kernels/_advection.py
@@ -18,12 +18,12 @@
 ]
 
 
-def rhs_uv(fieldset, time, depth, lat, lon, particles):
-    return fieldset.UV[time, depth, lat, lon, particles]
+def rhs_uv(fieldset, t, z, y, x, particles):
+    return fieldset.UV[t, z, y, x, particles]
 
 
-def rhs_uvw(fieldset, time, depth, lat, lon, particles):
-    return fieldset.UVW[time, depth, lat, lon, particles]
+def rhs_uvw(fieldset, t, z, y, x, particles):
+    return fieldset.UVW[t, z, y, x, particles]
 
 
 def AdvectionRK2(particles, fieldset):  # pragma: no cover
diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
index 228344977a..3df5700a97 100644
--- a/src/parcels/kernels/_integration.py
+++ b/src/parcels/kernels/_integration.py
@@ -2,45 +2,45 @@
 
 
 def RK2(fieldset, particles, rhs):
-    z, lat, lon = particles.z, particles.lat, particles.lon
-    fields = rhs(fieldset, particles.time, z, lat, lon, particles)
+    z, y, x = particles.z, particles.y, particles.x
+    fields = rhs(fieldset, particles.t, z, y, x, particles)
     if len(fields) == 1:
         raise NotImplementedError("RK2 integration is not implemented for Fields that return only one component.")
     if len(fields) > 1:
-        lon = particles.lon + fields[0] * 0.5 * particles.dt
-        lat = particles.lat + fields[1] * 0.5 * particles.dt
+        x = particles.x + fields[0] * 0.5 * particles.dt
+        y = particles.y + fields[1] * 0.5 * particles.dt
     if len(fields) > 2:
         z = particles.z + fields[2] * 0.5 * particles.dt
-    t = particles.time + 0.5 * particles.dt
-    return rhs(fieldset, t, z, lat, lon, particles)
+    t = particles.t + 0.5 * particles.dt
+    return rhs(fieldset, t, z, y, x, particles)
 
 
 def RK4(fieldset, particles, rhs):
-    z, lat, lon = particles.z, particles.lat, particles.lon
-    k1 = rhs(fieldset, particles.time, z, lat, lon, particles)
+    z, y, x = particles.z, particles.y, particles.x
+    k1 = rhs(fieldset, particles.t, z, y, x, particles)
     if len(k1) == 1:
         raise NotImplementedError("RK4 integration is not implemented for Fields that return only one component.")
     if len(k1) > 1:
-        lon = particles.lon + k1[0] * 0.5 * particles.dt
-        lat = particles.lat + k1[1] * 0.5 * particles.dt
+        x = particles.x + k1[0] * 0.5 * particles.dt
+        y = particles.y + k1[1] * 0.5 * particles.dt
     if len(k1) > 2:
         z = particles.z + k1[2] * 0.5 * particles.dt
-    t = particles.time + 0.5 * particles.dt
-    k2 = rhs(fieldset, t, z, lat, lon, particles)
+    t = particles.t + 0.5 * particles.dt
+    k2 = rhs(fieldset, t, z, y, x, particles)
     if len(k2) > 1:
-        lon = particles.lon + k2[0] * 0.5 * particles.dt
-        lat = particles.lat + k2[1] * 0.5 * particles.dt
+        x = particles.x + k2[0] * 0.5 * particles.dt
+        y = particles.y + k2[1] * 0.5 * particles.dt
     if len(k2) > 2:
         z = particles.z + k2[2] * 0.5 * particles.dt
-    t = particles.time + 0.5 * particles.dt
-    k3 = rhs(fieldset, t, z, lat, lon, particles)
+    t = particles.t + 0.5 * particles.dt
+    k3 = rhs(fieldset, t, z, y, x, particles)
     if len(k3) > 1:
-        lon = particles.lon + k3[0] * particles.dt
-        lat = particles.lat + k3[1] * particles.dt
+        x = particles.x + k3[0] * particles.dt
+        y = particles.y + k3[1] * particles.dt
     if len(k3) > 2:
         z = particles.z + k3[2] * particles.dt
-    t = particles.time + particles.dt
-    k4 = rhs(fieldset, t, z, lat, lon, particles)
+    t = particles.t + particles.dt
+    k4 = rhs(fieldset, t, z, y, x, particles)
     if len(k4) == 2:
         return ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0, (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0)
     if len(k4) == 3:

From 6d818122632deaaeb252b65f09d816b67f29fa46 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 8 Oct 2026 07:31:23 +0000
Subject: [PATCH 07/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 docs/user_guide/examples/tutorial_dt_integrators.ipynb | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/user_guide/examples/tutorial_dt_integrators.ipynb b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
index 1357ebba3e..b0c6ebda39 100644
--- a/docs/user_guide/examples/tutorial_dt_integrators.ipynb
+++ b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
@@ -850,7 +850,9 @@
     "\n",
     "fields = {\"Uw\": ds_w[\"eastward_wind\"], \"Vw\": ds_w[\"northward_wind\"]}\n",
     "ds_wind = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
-    "fieldset_w = parcels.FieldSet.from_sgrid_conventions(ds_wind, vector_fields={\"UVw\": (\"Uw\", \"Vw\")})\n",
+    "fieldset_w = parcels.FieldSet.from_sgrid_conventions(\n",
+    "    ds_wind, vector_fields={\"UVw\": (\"Uw\", \"Vw\")}\n",
+    ")\n",
     "\n",
     "fieldset += fieldset_w\n",
     "fieldset.describe()"

From 3fa27b2eea295304320030068edfa6c0831ba86c Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 9 Oct 2026 08:42:03 +0200
Subject: [PATCH 08/11] Implementing review suggestions

---
 src/parcels/kernels/_advection.py   | 8 ++++----
 src/parcels/kernels/_integration.py | 5 ++---
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py
index 1006a78566..03d3aa0f08 100644
--- a/src/parcels/kernels/_advection.py
+++ b/src/parcels/kernels/_advection.py
@@ -28,14 +28,14 @@ def rhs_uvw(fieldset, t, z, y, x, particles):
 
 def AdvectionRK2(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration."""
-    u, v = RK2(fieldset, particles, rhs_uv)
+    u, v = RK2(particles, fieldset, rhs_uv)
     particles.dx += u * particles.dt
     particles.dy += v * particles.dt
 
 
 def AdvectionRK2_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
-    u, v, w = RK2(fieldset, particles, rhs_uvw)
+    u, v, w = RK2(particles, fieldset, rhs_uvw)
     particles.dx += u * particles.dt
     particles.dy += v * particles.dt
     particles.dz += w * particles.dt
@@ -43,14 +43,14 @@ def AdvectionRK2_3D(particles, fieldset):  # pragma: no cover
 
 def AdvectionRK4(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
-    u, v = RK4(fieldset, particles, rhs_uv)
+    u, v = RK4(particles, fieldset, rhs_uv)
     particles.dx += u * particles.dt
     particles.dy += v * particles.dt
 
 
 def AdvectionRK4_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
-    u, v, w = RK4(fieldset, particles, rhs_uvw)
+    u, v, w = RK4(particles, fieldset, rhs_uvw)
     particles.dx += u * particles.dt
     particles.dy += v * particles.dt
     particles.dz += w * particles.dt
diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
index 3df5700a97..7d11ed0f1c 100644
--- a/src/parcels/kernels/_integration.py
+++ b/src/parcels/kernels/_integration.py
@@ -1,7 +1,7 @@
 """Collection of time integrators for use in Parcels Kernels"""
 
 
-def RK2(fieldset, particles, rhs):
+def RK2(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     fields = rhs(fieldset, particles.t, z, y, x, particles)
     if len(fields) == 1:
@@ -15,7 +15,7 @@ def RK2(fieldset, particles, rhs):
     return rhs(fieldset, t, z, y, x, particles)
 
 
-def RK4(fieldset, particles, rhs):
+def RK4(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     k1 = rhs(fieldset, particles.t, z, y, x, particles)
     if len(k1) == 1:
@@ -32,7 +32,6 @@ def RK4(fieldset, particles, rhs):
         y = particles.y + k2[1] * 0.5 * particles.dt
     if len(k2) > 2:
         z = particles.z + k2[2] * 0.5 * particles.dt
-    t = particles.t + 0.5 * particles.dt
     k3 = rhs(fieldset, t, z, y, x, particles)
     if len(k3) > 1:
         x = particles.x + k3[0] * particles.dt

From f9126c4fd39253703c61de76db2e6b3189b76c26 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 9 Oct 2026 09:05:03 +0200
Subject: [PATCH 09/11] Fix error message when calling integrators with one
 field

---
 src/parcels/kernels/__init__.py     |  4 ++++
 src/parcels/kernels/_integration.py | 17 +++++++++++++++++
 tests/test_advection.py             | 21 +++++++++++++++++++++
 3 files changed, 42 insertions(+)

diff --git a/src/parcels/kernels/__init__.py b/src/parcels/kernels/__init__.py
index 710fb633da..8ea1d0686a 100644
--- a/src/parcels/kernels/__init__.py
+++ b/src/parcels/kernels/__init__.py
@@ -12,6 +12,7 @@
     AdvectionDiffusionM1,
     DiffusionUniformKh,
 )
+from ._integration import RK2, RK4
 from ._sigmagrids import (
     AdvectionRK2_3D_CROCO,
     SampleOmegaCroco,
@@ -35,4 +36,7 @@
     "AdvectionRK2_3D_CROCO",
     "SampleOmegaCroco",
     "convert_z_to_sigma_croco",
+    # integration
+    "RK2",
+    "RK4",
 ]
diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
index 7d11ed0f1c..7456c66257 100644
--- a/src/parcels/kernels/_integration.py
+++ b/src/parcels/kernels/_integration.py
@@ -1,9 +1,22 @@
 """Collection of time integrators for use in Parcels Kernels"""
 
 
+def _validate_rhs_output(output, integrator_name):
+    """Validate that an rhs function returned a 2- or 3-component tuple (u, v) or (u, v, w)."""
+    if not isinstance(output, tuple):
+        raise TypeError(
+            f"{integrator_name}: rhs must return a tuple of (u, v) or (u, v, w), got {type(output).__name__}."
+        )
+    if len(output) not in (2, 3):
+        raise ValueError(
+            f"{integrator_name}: rhs must return a tuple of 2 (u, v) or 3 (u, v, w) components, got {len(output)}."
+        )
+
+
 def RK2(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     fields = rhs(fieldset, particles.t, z, y, x, particles)
+    _validate_rhs_output(fields, "RK2")
     if len(fields) == 1:
         raise NotImplementedError("RK2 integration is not implemented for Fields that return only one component.")
     if len(fields) > 1:
@@ -18,6 +31,7 @@ def RK2(particles, fieldset, rhs):
 def RK4(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     k1 = rhs(fieldset, particles.t, z, y, x, particles)
+    _validate_rhs_output(k1, "RK4")
     if len(k1) == 1:
         raise NotImplementedError("RK4 integration is not implemented for Fields that return only one component.")
     if len(k1) > 1:
@@ -27,12 +41,14 @@ def RK4(particles, fieldset, rhs):
         z = particles.z + k1[2] * 0.5 * particles.dt
     t = particles.t + 0.5 * particles.dt
     k2 = rhs(fieldset, t, z, y, x, particles)
+    _validate_rhs_output(k2, "RK4")
     if len(k2) > 1:
         x = particles.x + k2[0] * 0.5 * particles.dt
         y = particles.y + k2[1] * 0.5 * particles.dt
     if len(k2) > 2:
         z = particles.z + k2[2] * 0.5 * particles.dt
     k3 = rhs(fieldset, t, z, y, x, particles)
+    _validate_rhs_output(k3, "RK4")
     if len(k3) > 1:
         x = particles.x + k3[0] * particles.dt
         y = particles.y + k3[1] * particles.dt
@@ -40,6 +56,7 @@ def RK4(particles, fieldset, rhs):
         z = particles.z + k3[2] * particles.dt
     t = particles.t + particles.dt
     k4 = rhs(fieldset, t, z, y, x, particles)
+    _validate_rhs_output(k4, "RK4")
     if len(k4) == 2:
         return ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0, (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0)
     if len(k4) == 3:
diff --git a/tests/test_advection.py b/tests/test_advection.py
index 451e5e485c..cb85c7a993 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -28,6 +28,8 @@
 )
 from parcels._datasets.structured.generic import datasets_sgrid
 from parcels.kernels import (
+    RK2,
+    RK4,
     AdvectionDiffusionEM,
     AdvectionDiffusionM1,
     AdvectionEE,
@@ -232,6 +234,25 @@ def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice)
         np.testing.assert_allclose(np.array([p.z - z0 for p in pset]), 4 * w_value, atol=1e-5)
 
 
+@pytest.mark.parametrize("npart", [1, 2, 3, 10])
+@pytest.mark.parametrize("integrator", [RK2, RK4])
+def test_advection_with_integrators_singlefield(npart, integrator):
+    """Test that the RK2 and RK4 integrators don't work with a single field (e.g. U)"""
+    ds = simple_UV_dataset(mesh="flat").rename({"U": "T"})
+    fset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="flat")
+
+    pset = parcels.ParticleSet(fset, x=np.zeros(npart), y=np.arange(npart))
+
+    def rhs(fieldset, t, z, y, x, particles):
+        return fieldset.T[t, z, y, x, particles]
+
+    def SingleField(particles, fieldset):
+        integrator(particles, fieldset, rhs)
+
+    with pytest.raises(TypeError):
+        pset.execute(SingleField, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"))
+
+
 def test_radialrotation(npart=10):
     ds = radial_rotation_dataset()
     fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="flat")

From 80b2d755b25043b2f91bdd94e0d96798effb43c8 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 9 Oct 2026 09:08:26 +0200
Subject: [PATCH 10/11] Cleaning up duplicate error messages

Following on from f9126c4fd39253703c61de76db2e6b3189b76c26
---
 src/parcels/kernels/_integration.py | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
index 7456c66257..cb56770884 100644
--- a/src/parcels/kernels/_integration.py
+++ b/src/parcels/kernels/_integration.py
@@ -17,8 +17,6 @@ def RK2(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     fields = rhs(fieldset, particles.t, z, y, x, particles)
     _validate_rhs_output(fields, "RK2")
-    if len(fields) == 1:
-        raise NotImplementedError("RK2 integration is not implemented for Fields that return only one component.")
     if len(fields) > 1:
         x = particles.x + fields[0] * 0.5 * particles.dt
         y = particles.y + fields[1] * 0.5 * particles.dt
@@ -32,8 +30,6 @@ def RK4(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     k1 = rhs(fieldset, particles.t, z, y, x, particles)
     _validate_rhs_output(k1, "RK4")
-    if len(k1) == 1:
-        raise NotImplementedError("RK4 integration is not implemented for Fields that return only one component.")
     if len(k1) > 1:
         x = particles.x + k1[0] * 0.5 * particles.dt
         y = particles.y + k1[1] * 0.5 * particles.dt
@@ -41,14 +37,12 @@ def RK4(particles, fieldset, rhs):
         z = particles.z + k1[2] * 0.5 * particles.dt
     t = particles.t + 0.5 * particles.dt
     k2 = rhs(fieldset, t, z, y, x, particles)
-    _validate_rhs_output(k2, "RK4")
     if len(k2) > 1:
         x = particles.x + k2[0] * 0.5 * particles.dt
         y = particles.y + k2[1] * 0.5 * particles.dt
     if len(k2) > 2:
         z = particles.z + k2[2] * 0.5 * particles.dt
     k3 = rhs(fieldset, t, z, y, x, particles)
-    _validate_rhs_output(k3, "RK4")
     if len(k3) > 1:
         x = particles.x + k3[0] * particles.dt
         y = particles.y + k3[1] * particles.dt
@@ -56,7 +50,6 @@ def RK4(particles, fieldset, rhs):
         z = particles.z + k3[2] * particles.dt
     t = particles.t + particles.dt
     k4 = rhs(fieldset, t, z, y, x, particles)
-    _validate_rhs_output(k4, "RK4")
     if len(k4) == 2:
         return ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0, (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0)
     if len(k4) == 3:

From 85af999eb0e6c36f70ed3559416fa9a0ae798e93 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 9 Oct 2026 15:46:43 +0200
Subject: [PATCH 11/11] Reviewer feedback

---
 src/parcels/kernels/_integration.py | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/src/parcels/kernels/_integration.py b/src/parcels/kernels/_integration.py
index cb56770884..44d1c94923 100644
--- a/src/parcels/kernels/_integration.py
+++ b/src/parcels/kernels/_integration.py
@@ -17,9 +17,8 @@ def RK2(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     fields = rhs(fieldset, particles.t, z, y, x, particles)
     _validate_rhs_output(fields, "RK2")
-    if len(fields) > 1:
-        x = particles.x + fields[0] * 0.5 * particles.dt
-        y = particles.y + fields[1] * 0.5 * particles.dt
+    x = particles.x + fields[0] * 0.5 * particles.dt
+    y = particles.y + fields[1] * 0.5 * particles.dt
     if len(fields) > 2:
         z = particles.z + fields[2] * 0.5 * particles.dt
     t = particles.t + 0.5 * particles.dt
@@ -30,22 +29,19 @@ def RK4(particles, fieldset, rhs):
     z, y, x = particles.z, particles.y, particles.x
     k1 = rhs(fieldset, particles.t, z, y, x, particles)
     _validate_rhs_output(k1, "RK4")
-    if len(k1) > 1:
-        x = particles.x + k1[0] * 0.5 * particles.dt
-        y = particles.y + k1[1] * 0.5 * particles.dt
+    x = particles.x + k1[0] * 0.5 * particles.dt
+    y = particles.y + k1[1] * 0.5 * particles.dt
     if len(k1) > 2:
         z = particles.z + k1[2] * 0.5 * particles.dt
     t = particles.t + 0.5 * particles.dt
     k2 = rhs(fieldset, t, z, y, x, particles)
-    if len(k2) > 1:
-        x = particles.x + k2[0] * 0.5 * particles.dt
-        y = particles.y + k2[1] * 0.5 * particles.dt
+    x = particles.x + k2[0] * 0.5 * particles.dt
+    y = particles.y + k2[1] * 0.5 * particles.dt
     if len(k2) > 2:
         z = particles.z + k2[2] * 0.5 * particles.dt
     k3 = rhs(fieldset, t, z, y, x, particles)
-    if len(k3) > 1:
-        x = particles.x + k3[0] * particles.dt
-        y = particles.y + k3[1] * particles.dt
+    x = particles.x + k3[0] * particles.dt
+    y = particles.y + k3[1] * particles.dt
     if len(k3) > 2:
         z = particles.z + k3[2] * particles.dt
     t = particles.t + particles.dt