diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb
index a81b4cb97..19b48519e 100644
--- a/docs/user_guide/examples/tutorial_schism.ipynb
+++ b/docs/user_guide/examples/tutorial_schism.ipynb
@@ -174,8 +174,24 @@
     "U = u.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(\n",
     "    rename\n",
     ")  # reverse to surface-first\n",
-    "V = v.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename)\n",
-    "\n",
+    "V = v.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "Now, we can construct the `ux.UxDataset` object that we will use to create the `Fieldset` in Parcels with `Field` objects for `U` and `V`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "uxds = ux.UxDataset(\n",
     "    xr.Dataset(\n",
     "        {\n",
@@ -191,7 +207,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "9",
+   "id": "11",
    "metadata": {},
    "source": [
     "## Build the `FieldSet`\n",
@@ -207,7 +223,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "10",
+   "id": "12",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -220,29 +236,30 @@
   },
   {
    "cell_type": "markdown",
-   "id": "11",
+   "id": "13",
    "metadata": {},
    "source": [
     "## Stop particles that end up below the bathymetry\n",
     "\n",
     "Because SCHISM's LSC2 grid keeps below-seabed levels as `NaN`, a particle whose (fixed) depth ends up\n",
     "deeper than the local water column will sample `NaN`. Rather than feed it a fabricated velocity, we let\n",
-    "Parcels flag it: sampling `NaN` sets the particle state to {py:obj}`parcels.StatusCode.ErrorInterpolation`, and ending up outside\n",
-    "the mesh sets {py:obj}`parcels.StatusCode.ErrorOutOfBounds`. We add a small kernel that runs after advection, records that the\n",
-    "particle went out of bounds, and sets its state to {py:obj}`parcels.StatusCode.Delete` so it stops being advected (its trajectory\n",
+    "Parcels flag it: sampling `NaN` sets the particle state to `parcels.StatusCode.ErrorInterpolation`, and ending up outside\n",
+    "the mesh sets `parcels.StatusCode.ErrorOutOfBounds` or `parcels.StatusCode.ErrorGridSearching`. We add a small kernel that runs after advection, records that the\n",
+    "particle went out of bounds, and sets its state to `parcels.StatusCode.Delete` so it stops being advected (its trajectory\n",
     "up to that point is kept).\n",
     "\n",
     "```{note}\n",
     "Most of the Lake Ontario mesh is shallow nearshore water with only a couple of valid vertical levels;\n",
-    "the deep central basin has the full set. We therefore release particles in the deep basin (nodes with\n",
-    "many valid levels) so they start within the resolved water column.\n",
+    "the deep central basin has the full set. We therefore release most particles at 2 m depth in the deep\n",
+    "basin (nodes with many valid levels), and a smaller set at 1 m depth in the finely resolved nearshore\n",
+    "regions.\n",
     "```"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "12",
+   "id": "14",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -250,6 +267,7 @@
     "    parcels.StatusCode.ErrorInterpolation,  # sampled NaN (below the local seabed)\n",
     "    parcels.StatusCode.ErrorOutOfBounds,  # left the horizontal mesh / below the grid\n",
     "    parcels.StatusCode.ErrorThroughSurface,  # above the surface\n",
+    "    parcels.StatusCode.ErrorGridSearching,  # left the mesh during a step\n",
     "]\n",
     "\n",
     "SchismParticle = parcels.Particle.add_variable(\n",
@@ -266,16 +284,104 @@
   },
   {
    "cell_type": "markdown",
-   "id": "13",
+   "id": "15",
    "metadata": {},
    "source": [
-    "## Release particles and advect"
+    "## Limit the timestep by the local CFL number\n",
+    "\n",
+    "Unstructured meshes are often strongly refined near the coast, so a single `dt` that suits the large\n",
+    "cells in the basin can carry a particle across several small cells in one step. This is measured by\n",
+    "the CFL (Courant–Friedrichs–Lewy) number\n",
+    "\n",
+    "$$\n",
+    "C = \\frac{|u|\\,\\Delta t}{\\Delta x},\n",
+    "$$\n",
+    "\n",
+    "the number of cells of width $\\Delta x$ that a particle moving at speed $|u|$ crosses in one timestep\n",
+    "$\\Delta t$. Keeping $C \\lesssim 1$ ensures that particles do not skip over cells. See\n",
+    "[a priori estimation of `dt`](./tutorial_dt_integrators.ipynb#a-priori-estimation-of-dt) for the same\n",
+    "reasoning applied to a structured grid.\n",
+    "\n",
+    "To use the local cell width in a kernel, we add a `cell_width` field to the `FieldSet`. While there are a number of valid measures of the cell width of a triangle, here we use the shortest distance from a corner to the opposite edge. This is equivalent to twice the triangle's area divided by its longest edge. The cell width is defined on the faces and does not change in time or depth, so we broadcast it over the `time` and `zc` dimensions, build a `FieldSet` from it, and add that to the velocity `FieldSet`."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "14",
+   "id": "16",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# corner coordinates of each triangle, shape (n_face, 3, 2)\n",
+    "corners = np.stack([node_x[face_nodes], node_y[face_nodes]], axis=-1)\n",
+    "edges = np.roll(corners, -1, axis=1) - corners  # edge vectors\n",
+    "edge_lengths = np.hypot(edges[..., 0], edges[..., 1])\n",
+    "cell_area = 0.5 * np.abs(\n",
+    "    edges[:, 0, 0] * edges[:, 1, 1] - edges[:, 0, 1] * edges[:, 1, 0]\n",
+    ")\n",
+    "cell_width = (\n",
+    "    2.0 * cell_area / edge_lengths.max(axis=1)\n",
+    ")  # smallest height of each triangle\n",
+    "\n",
+    "cell_width_da = xr.DataArray(\n",
+    "    np.broadcast_to(cell_width, (U.sizes[\"time\"], zc.size, cell_width.size)),\n",
+    "    dims=(\"time\", \"zc\", \"n_face\"),\n",
+    "    coords={\"time\": U[\"time\"]},\n",
+    ")\n",
+    "uxds_cell_width = ux.UxDataset(\n",
+    "    xr.Dataset(\n",
+    "        {\"cell_width\": cell_width_da}, coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc)}\n",
+    "    ),\n",
+    "    uxgrid=uxgrid,\n",
+    ")\n",
+    "fieldset = fieldset + parcels.FieldSet.from_ugrid_conventions(uxds_cell_width)\n",
+    "fieldset.describe()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17",
+   "metadata": {},
+   "source": [
+    "The `CFLLimitedDt` kernel below shrinks each particle's `dt` so that it moves at most `CFL_target`\n",
+    "local cell widths per step: $\\Delta t \\le C\\,\\Delta x/|u|$, with $C$ = `CFL_target` and\n",
+    "$\\Delta x$ = `cell_width`. To expose the `CFL_target` as a value in the kernel, we use the `fieldset.add_context(name,value)` method to store a scalar value under the name `CFL_target`. The `dt` passed to `pset.execute()` then acts as the *maximum* timestep, and\n",
+    "Parcels sub-steps each particle until it reaches the next output time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def CFLLimitedDt(particles, fieldset):\n",
+    "    \"\"\"Shrink each particle's dt so it moves at most CFL_target cell widths per step.\"\"\"\n",
+    "    cell_width = fieldset.cell_width[particles]\n",
+    "    u, v = fieldset.UV[particles]\n",
+    "    dt_cfl = (\n",
+    "        fieldset.CFL_target * cell_width / np.maximum(np.sqrt(u * u + v * v), 1e-12)\n",
+    "    )\n",
+    "    particles.dt = np.sign(particles.dt) * np.minimum(np.abs(particles.dt), dt_cfl)\n",
+    "\n",
+    "\n",
+    "fieldset.add_context(\"CFL_target\", 1.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "19",
+   "metadata": {},
+   "source": [
+    "## Release particles and advect\n",
+    "Here, the initial particle positions are selected at the face centroids in the deep basin and in areas where the cells are expected to be small."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -290,19 +396,47 @@
     "idx = np.where(deep)[0]\n",
     "idx = idx[:: max(1, idx.size // n_max)][:n_max]\n",
     "\n",
-    "x, y = cx[idx], cy[idx]\n",
-    "z = np.full(x.size, 2.0)  # release at 2 m depth\n",
-    "print(f\"releasing {x.size} particles at z = 2 m in the deep basin\")\n",
     "\n",
+    "# release 500 more on faces where the CFL limit at dt = 5 min is below 1 (small cells, fast flow), at 1 m depth\n",
+    "surf_speed = np.hypot(U.isel(time=0, zf=0).values, V.isel(time=0, zf=0).values)\n",
+    "dt_cfl = cell_width / np.maximum(\n",
+    "    np.nan_to_num(surf_speed[face_nodes].mean(axis=1)), 1e-12\n",
+    ")\n",
+    "fine = np.where(dt_cfl < 300)[0]\n",
+    "idx_fine = np.random.default_rng(0).choice(fine, 500, replace=False)\n",
+    "\n",
+    "x = np.concatenate([cx[idx], cx[idx_fine]])\n",
+    "y = np.concatenate([cy[idx], cy[idx_fine]])\n",
+    "z = np.concatenate([np.full(idx.size, 2.0), np.full(idx_fine.size, 1.0)])\n",
+    "print(\n",
+    "    f\"releasing {idx.size} particles in the deep basin and {idx_fine.size} in fine nearshore cells\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "21",
+   "metadata": {},
+   "source": [
+    "When calling `pset.execute` we call the `CFLLimitedDt` kernel before the `parcels.kernels.AdvectionRK4`. With each execution loop, each particle's time step size is adjusted so that its local CFL estimate is at most `CFL_target` before being advected."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "pset = parcels.ParticleSet(fieldset=fieldset, pclass=SchismParticle, x=x, y=y, z=z)\n",
     "output_file = parcels.ParticleFile(\n",
     "    \"output-schism.parquet\", outputdt=np.timedelta64(30, \"m\")\n",
     ")\n",
     "\n",
     "pset.execute(\n",
-    "    [parcels.kernels.AdvectionRK4, StopBelowBed],\n",
+    "    [CFLLimitedDt, parcels.kernels.AdvectionRK4, StopBelowBed],\n",
     "    runtime=np.timedelta64(5, \"h\"),  # the dataset spans 5 hours\n",
-    "    dt=np.timedelta64(5, \"m\"),\n",
+    "    dt=np.timedelta64(5, \"m\"),  # maximum timestep\n",
     "    output_file=output_file,\n",
     "    verbose_progress=False,\n",
     ")\n",
@@ -311,7 +445,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "15",
+   "id": "23",
    "metadata": {},
    "source": [
     "## Plot the velocity field and trajectories\n",
@@ -325,7 +459,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "16",
+   "id": "24",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -388,7 +522,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "17",
+   "id": "25",
    "metadata": {},
    "source": [
     "The particles move with the SCHISM surface currents, and any that wander over shallow water where\n",
@@ -406,7 +540,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "18",
+   "id": "26",
    "metadata": {},
    "source": [
     "## Adapting this tutorial to your own SCHISM setup\n",
@@ -420,7 +554,8 @@
     "  $\\Delta x$ approximated by each triangle's shortest edge) is well above 1 in the narrow channel, so\n",
     "  particles can jump across several cells, or out of the mesh entirely, in a single step. Choose `dt`\n",
     "  from the *smallest* cells and *fastest* flows that your particles will visit, not from the domain\n",
-    "  average. See [choosing a timestep and integration method](./tutorial_dt_integrators.ipynb).\n",
+    "  average, or use the `CFLLimitedDt` kernel from this tutorial to adapt `dt` per particle. See\n",
+    "  [choosing a timestep and integration method](./tutorial_dt_integrators.ipynb).\n",
     "\n",
     "  ![CFL number for u and v near a SCHISM river inflow with dt = 5 minutes](../../_static/schism_inflow_cfl.png)\n",
     "\n",
@@ -437,15 +572,15 @@
     "\n",
     "* **Particles leaving the mesh.** Because of the two points above, some particles near the coast may\n",
     "  still exit the mesh. Decide how to handle them: release particles away from the boundaries (as done\n",
-    "  here with the `valid_levels` filter, whose threshold depends on your mesh), or catch the relevant\n",
-    "  status codes (e.g. add `parcels.StatusCode.ErrorGridSearching` to `OUT_OF_BOUNDS_STATES`) and stop\n",
-    "  those particles. See [working with status codes](./tutorial_statuscodes.md)."
+    "  here for the deep-basin set with the `valid_levels` filter, whose threshold depends on your mesh),\n",
+    "  or catch the relevant status codes (as done here with `OUT_OF_BOUNDS_STATES`) and stop those\n",
+    "  particles. See [working with status codes](./tutorial_statuscodes.md)."
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Parcels:default (3.14.8)",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -459,7 +594,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.14.8"
+   "version": "3.14.3"
   },
   "mystnb": {
    "execution_in_temp": true