diff --git a/.gitignore b/.gitignore
index c34fabf59c..8f4c7fd1ca 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,7 +1,6 @@
 build/*
 docs/_build/*
 docs/_downloads
-gh-pages/*
 lib/
 bin/
 parcels_examples
@@ -19,15 +18,13 @@ out-*
 .DS_store
 
 .vscode
-.idea/*
 .env
 Profile.prof
 **/.ipynb_checkpoints/*
 .cache/*
-.pytest_cache/*
 .eggs/*
 parcels.egg-info/*
 dist/parcels*.egg
-parcels/examples/particle*.png
 parcels/_version_setup.py
-/.pytest_cache/
+.pytest_cache
+.coverage
diff --git a/docs/examples/example_dask_chunk_OCMs.py b/docs/examples/example_dask_chunk_OCMs.py
index c03d9174d7..27dd8c2332 100644
--- a/docs/examples/example_dask_chunk_OCMs.py
+++ b/docs/examples/example_dask_chunk_OCMs.py
@@ -17,7 +17,7 @@ def compute_nemo_particle_advection(fieldset, mode):
     lonp = 2.5 * np.ones(npart)
     latp = [i for i in 52.0 + (-1e-3 + np.random.rand(npart) * 2.0 * 1e-3)]
 
-    def periodicBC(particle, fieldSet, time):
+    def periodicBC(particle, fieldSet, time):  # pragma: no cover
         if particle.lon > 15.0:
             particle_dlon -= 15.0  # noqa
         if particle.lon < 0:
@@ -652,7 +652,7 @@ class MyParticle(ptype[mode]):
 
     pset = parcels.ParticleSet(fieldset, pclass=MyParticle, lon=2.5, lat=52)
 
-    def Sample2D(particle, fieldset, time):
+    def Sample2D(particle, fieldset, time):  # pragma: no cover
         particle.sample_var_curvilinear += fieldset.nav_lon[
             time, particle.depth, particle.lat, particle.lon
         ]
diff --git a/docs/examples/example_globcurrent.py b/docs/examples/example_globcurrent.py
index d310a9883a..292ed1bf7e 100755
--- a/docs/examples/example_globcurrent.py
+++ b/docs/examples/example_globcurrent.py
@@ -129,7 +129,7 @@ def test_globcurrent_time_periodic(mode, rundays):
             fieldset, pclass=MyParticle, lon=25, lat=-35, time=fieldset.U.grid.time[0]
         )
 
-        def SampleU(particle, fieldset, time):
+        def SampleU(particle, fieldset, time):  # pragma: no cover
             u, v = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
             particle.sample_var += u
 
@@ -255,7 +255,7 @@ def test_globcurrent_startparticles_between_time_arrays(mode, dt, with_starttime
 
     MyParticle = ptype[mode].add_variable("sample_var", initial=0.0)
 
-    def SampleP(particle, fieldset, time):
+    def SampleP(particle, fieldset, time):  # pragma: no cover
         particle.sample_var += fieldset.P[
             time, particle.depth, particle.lat, particle.lon
         ]
@@ -288,7 +288,7 @@ def test_globcurrent_particle_independence(mode, rundays=5):
     fieldset = set_globcurrent_fieldset()
     time0 = fieldset.U.grid.time[0]
 
-    def DeleteP0(particle, fieldset, time):
+    def DeleteP0(particle, fieldset, time):  # pragma: no cover
         if particle.id == 0:
             particle.delete()
 
@@ -353,7 +353,7 @@ def test_error_outputdt_not_multiple_dt(mode, tmpdir):
     pset = parcels.ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0])
     ofile = pset.ParticleFile(name=filepath, outputdt=timedelta(days=1))
 
-    def DoNothing(particle, fieldset, time):
+    def DoNothing(particle, fieldset, time):  # pragma: no cover
         pass
 
     with pytest.raises(ValueError):
diff --git a/docs/examples/example_mitgcm.py b/docs/examples/example_mitgcm.py
index 95fec26ed7..22e66faf55 100644
--- a/docs/examples/example_mitgcm.py
+++ b/docs/examples/example_mitgcm.py
@@ -29,7 +29,7 @@ def run_mitgcm_zonally_reentrant(mode: Literal["scipy", "jit"], path: Path):
     fieldset.add_periodic_halo(zonal=True)
     fieldset.add_constant("domain_width", 1000000)
 
-    def periodicBC(particle, fieldset, time):
+    def periodicBC(particle, fieldset, time):  # pragma: no cover
         if particle.lon < 0:
             particle_dlon += fieldset.domain_width  # noqa
         elif particle.lon > fieldset.domain_width:
diff --git a/docs/examples/example_moving_eddies.py b/docs/examples/example_moving_eddies.py
index 40f679539a..cff035ec5c 100644
--- a/docs/examples/example_moving_eddies.py
+++ b/docs/examples/example_moving_eddies.py
@@ -269,7 +269,7 @@ def test_periodic_and_computeTimeChunk_eddies(mode):
         fieldset=fieldset, pclass=ptype[mode], lon=[3.3, 3.3], lat=[46.0, 47.8]
     )
 
-    def periodicBC(particle, fieldset, time):
+    def periodicBC(particle, fieldset, time):  # pragma: no cover
         if particle.lon < fieldset.halo_west:
             particle_dlon += fieldset.halo_east - fieldset.halo_west  # noqa
         elif particle.lon > fieldset.halo_east:
@@ -279,7 +279,7 @@ def periodicBC(particle, fieldset, time):
         elif particle.lat > fieldset.halo_north:
             particle_dlat -= fieldset.halo_north - fieldset.halo_south
 
-    def slowlySouthWestward(particle, fieldset, time):
+    def slowlySouthWestward(particle, fieldset, time):  # pragma: no cover
         particle_dlon -= 5 * particle.dt / 1e5  # noqa
         particle_dlat -= 3 * particle.dt / 1e5  # noqa
 
diff --git a/docs/examples/example_nemo_curvilinear.py b/docs/examples/example_nemo_curvilinear.py
index e394573313..dd751b4336 100644
--- a/docs/examples/example_nemo_curvilinear.py
+++ b/docs/examples/example_nemo_curvilinear.py
@@ -47,7 +47,7 @@ def run_nemo_curvilinear(mode, outfile, advtype="RK4"):
         latp = np.linspace(-70, 70, npart)
         runtime = timedelta(days=15)
 
-    def periodicBC(particle, fieldSet, time):
+    def periodicBC(particle, fieldSet, time):  # pragma: no cover
         if particle.lon > 180:
             particle_dlon -= 360  # noqa
 
diff --git a/docs/examples/example_peninsula.py b/docs/examples/example_peninsula.py
index 2ecf888596..a4e61b4f3b 100644
--- a/docs/examples/example_peninsula.py
+++ b/docs/examples/example_peninsula.py
@@ -91,7 +91,7 @@ def peninsula_fieldset(xdim, ydim, mesh="flat", grid_type="A"):
     return fieldset
 
 
-def UpdateP(particle, fieldset, time):
+def UpdateP(particle, fieldset, time):  # pragma: no cover
     if time == 0:
         particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
diff --git a/docs/examples/example_stommel.py b/docs/examples/example_stommel.py
index 8fefaded17..d5c9d7583c 100755
--- a/docs/examples/example_stommel.py
+++ b/docs/examples/example_stommel.py
@@ -77,13 +77,13 @@ def stommel_fieldset(xdim=200, ydim=200, grid_type="A"):
     return fieldset
 
 
-def UpdateP(particle, fieldset, time):
+def UpdateP(particle, fieldset, time):  # pragma: no cover
     if time == 0:
         particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
-def AgeP(particle, fieldset, time):
+def AgeP(particle, fieldset, time):  # pragma: no cover
     particle.age += particle.dt
     if particle.age > fieldset.maxage:
         particle.delete()
diff --git a/parcels/application_kernels/EOSseawaterproperties.py b/parcels/application_kernels/EOSseawaterproperties.py
index d85b6633a2..962b945ab4 100644
--- a/parcels/application_kernels/EOSseawaterproperties.py
+++ b/parcels/application_kernels/EOSseawaterproperties.py
@@ -5,7 +5,7 @@
 __all__ = ["AdiabticTemperatureGradient", "PressureFromLatDepth", "PtempFromTemp", "TempFromPtemp", "UNESCODensity"]
 
 
-def PressureFromLatDepth(particle, fieldset, time):
+def PressureFromLatDepth(particle, fieldset, time):  # pragma: no cover
     """
     Calculates pressure in dbars from depth in meters and latitude.
 
@@ -28,7 +28,7 @@ def PressureFromLatDepth(particle, fieldset, time):
     particle.pressure = ((1 - C1) - math.pow(((math.pow((1 - C1), 2)) - (8.84e-6 * particle.depth)), 0.5)) / 4.42e-6
 
 
-def AdiabticTemperatureGradient(particle, fieldset, time):
+def AdiabticTemperatureGradient(particle, fieldset, time):  # pragma: no cover
     """Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
 
 
@@ -81,7 +81,7 @@ def AdiabticTemperatureGradient(particle, fieldset, time):
     )
 
 
-def PtempFromTemp(particle, fieldset, time):
+def PtempFromTemp(particle, fieldset, time):  # pragma: no cover
     """
     Calculates potential temperature as per UNESCO 1983 report.
 
@@ -188,7 +188,7 @@ def PtempFromTemp(particle, fieldset, time):
     particle.potemp = (th + (del_th - 2 * q) / 6) / 1.00024
 
 
-def TempFromPtemp(particle, fieldset, time):
+def TempFromPtemp(particle, fieldset, time):  # pragma: no cover
     """
     Calculates temperature from potential temperature at the reference
     pressure PR and in situ pressure P.
@@ -295,7 +295,7 @@ def TempFromPtemp(particle, fieldset, time):
     particle.temp = (th + (del_th - 2 * q) / 6) / 1.00024
 
 
-def UNESCODensity(particle, fieldset, time):
+def UNESCODensity(particle, fieldset, time):  # pragma: no cover
     # This is a kernel which calculates the UNESCO density
     # (https://link.springer.com/content/pdf/bbm%3A978-3-319-18908-6%2F1.pdf),
     # from pressure, temperature and salinity.
diff --git a/parcels/application_kernels/TEOSseawaterdensity.py b/parcels/application_kernels/TEOSseawaterdensity.py
index 0195b0d13e..8234dd52f6 100644
--- a/parcels/application_kernels/TEOSseawaterdensity.py
+++ b/parcels/application_kernels/TEOSseawaterdensity.py
@@ -5,7 +5,7 @@
 __all__ = ["PolyTEOS10_bsq"]
 
 
-def PolyTEOS10_bsq(particle, fieldset, time):
+def PolyTEOS10_bsq(particle, fieldset, time):  # pragma: no cover
     """Calculates density based on the polyTEOS10-bsq algorithm from Appendix A.2 of
     https://www.sciencedirect.com/science/article/pii/S1463500315000566
     requires fieldset.abs_salinity and fieldset.cons_temperature Fields in the fieldset
diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py
index b30513053c..b11d9912fa 100644
--- a/parcels/application_kernels/advection.py
+++ b/parcels/application_kernels/advection.py
@@ -14,7 +14,7 @@
 ]
 
 
-def AdvectionRK4(particle, fieldset, time):
+def AdvectionRK4(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
     (u1, v1) = fieldset.UV[particle]
     lon1, lat1 = (particle.lon + u1 * 0.5 * particle.dt, particle.lat + v1 * 0.5 * particle.dt)
@@ -27,7 +27,7 @@ def AdvectionRK4(particle, fieldset, time):
     particle_dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particle.dt  # noqa
 
 
-def AdvectionRK4_3D(particle, fieldset, time):
+def AdvectionRK4_3D(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
     (u1, v1, w1) = fieldset.UVW[particle]
     lon1 = particle.lon + u1 * 0.5 * particle.dt
@@ -47,7 +47,7 @@ def AdvectionRK4_3D(particle, fieldset, time):
     particle_ddepth += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * particle.dt  # noqa
 
 
-def AdvectionRK4_3D_CROCO(particle, fieldset, time):
+def AdvectionRK4_3D_CROCO(particle, fieldset, time):  # 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.
     """
@@ -92,14 +92,14 @@ def AdvectionRK4_3D_CROCO(particle, fieldset, time):
     ) / 6
 
 
-def AdvectionEE(particle, fieldset, time):
+def AdvectionEE(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using Explicit Euler (aka Euler Forward) integration."""
     (u1, v1) = fieldset.UV[particle]
     particle_dlon += u1 * particle.dt  # noqa
     particle_dlat += v1 * particle.dt  # noqa
 
 
-def AdvectionRK45(particle, fieldset, time):
+def AdvectionRK45(particle, fieldset, time):  # 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'
@@ -161,7 +161,7 @@ def AdvectionRK45(particle, fieldset, time):
         return StatusCode.Repeat
 
 
-def AdvectionAnalytical(particle, fieldset, time):
+def AdvectionAnalytical(particle, fieldset, time):  # 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).
@@ -232,14 +232,14 @@ def AdvectionAnalytical(particle, fieldset, time):
     else:
         dz = 1.0
 
-    c1 = fieldset.UV.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
-    c2 = fieldset.UV.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
-    c3 = fieldset.UV.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
-    c4 = fieldset.UV.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
+    c1 = i_u._geodetic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
+    c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
+    c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
+    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
-    dxdy = fieldset.UV.jacobian(py, px, eta, xsi) * meshJac
+    dxdy = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
 
     if withW:
         U0 = direction * fieldset.U.data[ti, zi + 1, yi + 1, xi] * c4 * dz
diff --git a/parcels/application_kernels/advectiondiffusion.py b/parcels/application_kernels/advectiondiffusion.py
index fc3c8594ad..461063ac78 100644
--- a/parcels/application_kernels/advectiondiffusion.py
+++ b/parcels/application_kernels/advectiondiffusion.py
@@ -10,7 +10,7 @@
 __all__ = ["AdvectionDiffusionEM", "AdvectionDiffusionM1", "DiffusionUniformKh"]
 
 
-def AdvectionDiffusionM1(particle, fieldset, time):
+def AdvectionDiffusionM1(particle, fieldset, time):  # 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`
@@ -47,7 +47,7 @@ def AdvectionDiffusionM1(particle, fieldset, time):
     particle_dlat += v * particle.dt + 0.5 * dKdy * (dWy**2 + particle.dt) + by * dWy  # noqa
 
 
-def AdvectionDiffusionEM(particle, fieldset, time):
+def AdvectionDiffusionEM(particle, fieldset, time):  # 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`
@@ -83,7 +83,7 @@ def AdvectionDiffusionEM(particle, fieldset, time):
     particle_dlat += ay * particle.dt + by * dWy  # noqa
 
 
-def DiffusionUniformKh(particle, fieldset, time):
+def DiffusionUniformKh(particle, fieldset, time):  # 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`.
diff --git a/parcels/field.py b/parcels/field.py
index c60447ad3c..a02b53a49e 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -1977,28 +1977,17 @@ def _check_grid_dimensions(grid1, grid2):
             and np.allclose(grid1.time_full, grid2.time_full)
         )
 
-    def dist(self, lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float):
-        if mesh == "spherical":
-            rad = np.pi / 180.0
-            deg2m = 1852 * 60.0
-            return np.sqrt(((lon2 - lon1) * deg2m * math.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2)
-        else:
-            return np.sqrt((lon2 - lon1) ** 2 + (lat2 - lat1) ** 2)
-
-    def jacobian(self, py: np.ndarray, px: np.ndarray, eta: float, xsi: float):
-        dphidxsi = [eta - 1, 1 - eta, eta, -eta]
-        dphideta = [xsi - 1, -xsi, xsi, 1 - xsi]
+    @deprecated_made_private  # TODO: Remove 6 months after v3.2.0
+    def dist(self, *args, **kwargs):
+        raise NotImplementedError
 
-        dxdxsi = np.dot(px, dphidxsi)
-        dxdeta = np.dot(px, dphideta)
-        dydxsi = np.dot(py, dphidxsi)
-        dydeta = np.dot(py, dphideta)
-        jac = dxdxsi * dydeta - dxdeta * dydxsi
-        return jac
+    @deprecated_made_private  # TODO: Remove 6 months after v3.2.0
+    def jacobian(self, *args, **kwargs):
+        raise NotImplementedError
 
     def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, applyConversion=True):
         grid = self.U.grid
-        (zeta, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
+        (_, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
 
         if grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]])
@@ -2014,10 +2003,10 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply
             px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:])
         xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3]
         assert abs(xx - x) < 1e-4
-        c1 = self.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
-        c2 = self.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
-        c3 = self.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
-        c4 = self.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
+        c1 = i_u._geodetic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
+        c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
+        c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
+        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))
         if grid.zdim == 1:
             if self.gridindexingtype == "nemo":
                 U0 = self.U.data[ti, yi + 1, xi] * c4
@@ -2049,7 +2038,7 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply
         else:
             meshJac = deg2m if grid.mesh == "spherical" else 1
 
-        jac = self.jacobian(py, px, eta, xsi) * meshJac
+        jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
 
         u = (
             (-(1 - eta) * U - (1 - xsi) * V) * px[0]
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 3f78eb4f41..a8f6d7ff90 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -300,7 +300,7 @@ def _cache_key(self):
 
     def add_scipy_positionupdate_kernels(self):
         # Adding kernels that set and update the coordinate changes
-        def Setcoords(particle, fieldset, time):
+        def Setcoords(particle, fieldset, time):  # pragma: no cover
             particle_dlon = 0  # noqa
             particle_dlat = 0  # noqa
             particle_ddepth = 0  # noqa
@@ -309,7 +309,7 @@ def Setcoords(particle, fieldset, time):
             particle.depth = particle.depth_nextloop
             particle.time = particle.time_nextloop
 
-        def Updatecoords(particle, fieldset, time):
+        def Updatecoords(particle, fieldset, time):  # pragma: no cover
             particle.lon_nextloop = particle.lon + particle_dlon  # type: ignore[name-defined] # noqa
             particle.lat_nextloop = particle.lat + particle_dlat  # type: ignore[name-defined] # noqa
             particle.depth_nextloop = particle.depth + particle_ddepth  # type: ignore[name-defined] # noqa
diff --git a/parcels/tools/interpolation_utils.py b/parcels/tools/interpolation_utils.py
index 03ebcda240..71d316e7cc 100644
--- a/parcels/tools/interpolation_utils.py
+++ b/parcels/tools/interpolation_utils.py
@@ -1,3 +1,4 @@
+import math
 from collections.abc import Callable
 from typing import Literal
 
@@ -144,6 +145,7 @@ def dphidxsi2D_lin(eta: float, xsi: float) -> tuple[list[float], list[float]]:
                   1-xsi]
 
     return dphideta, dphidxsi
+# fmt: on
 
 
 def dxdxsi2D_lin(
@@ -173,4 +175,22 @@ def interpolate(phi: Callable[[float], list[float]], f: list[float], xsi: float)
     return np.dot(phi(xsi), f)
 
 
-# fmt: on
+def _geodetic_distance(lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float) -> float:
+    if mesh == "spherical":
+        rad = np.pi / 180.0
+        deg2m = 1852 * 60.0
+        return np.sqrt(((lon2 - lon1) * deg2m * math.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2)
+    else:
+        return np.sqrt((lon2 - lon1) ** 2 + (lat2 - lat1) ** 2)
+
+
+def _compute_jacobian_determinant(py: np.ndarray, px: np.ndarray, eta: float, xsi: float) -> float:
+    dphidxsi = [eta - 1, 1 - eta, eta, -eta]
+    dphideta = [xsi - 1, -xsi, xsi, 1 - xsi]
+
+    dxdxsi = np.dot(px, dphidxsi)
+    dxdeta = np.dot(px, dphideta)
+    dydxsi = np.dot(py, dphidxsi)
+    dydeta = np.dot(py, dphideta)
+    jac = dxdxsi * dydeta - dxdeta * dydxsi
+    return jac
diff --git a/pyproject.toml b/pyproject.toml
index af6d4fb8c5..70d202689a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -149,7 +149,7 @@ module = [
     "mpi4py",
     "scipy.spatial",
     "sklearn.cluster",
-    "zarr>=2.11.0,!=2.18.0,<3",
+    "zarr",
     "cftime",
     "pykdtree.kdtree",
     "netCDF4",
diff --git a/tests/common_kernels.py b/tests/common_kernels.py
index 1eefd91115..449ea9039b 100644
--- a/tests/common_kernels.py
+++ b/tests/common_kernels.py
@@ -1,18 +1,18 @@
 """Shared kernels between tests."""
 
 
-def DoNothing(particle, fieldset, time):
+def DoNothing(particle, fieldset, time):  # pragma: no cover
     pass
 
 
-def DeleteParticle(particle, fieldset, time):
+def DeleteParticle(particle, fieldset, time):  # pragma: no cover
     if particle.state >= 50:  # This captures all Errors
         particle.delete()
 
 
-def MoveEast(particle, fieldset, time):
+def MoveEast(particle, fieldset, time):  # pragma: no cover
     particle_dlon += 0.1  # noqa
 
 
-def MoveNorth(particle, fieldset, time):
+def MoveNorth(particle, fieldset, time):  # pragma: no cover
     particle_dlat += 0.1  # noqa
diff --git a/tests/test_advection.py b/tests/test_advection.py
index 77705f0576..2ded597993 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -145,11 +145,11 @@ def test_advection_3D_outofbounds(mode, direction, wErrorThroughSurface):
     }
     fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
 
-    def DeleteParticle(particle, fieldset, time):
+    def DeleteParticle(particle, fieldset, time):  # pragma: no cover
         if particle.state == StatusCode.ErrorOutOfBounds or particle.state == StatusCode.ErrorThroughSurface:
             particle.delete()
 
-    def SubmergeParticle(particle, fieldset, time):
+    def SubmergeParticle(particle, fieldset, time):  # pragma: no cover
         if particle.state == StatusCode.ErrorThroughSurface:
             (u, v) = fieldset.UV[particle]
             particle_dlon = u * particle.dt  # noqa
@@ -246,7 +246,7 @@ def test_advection_3DCROCO(mode):
     pclass = ptype[mode].add_variable(Variable("w"))
     pset = ParticleSet(fieldset=fieldset, pclass=pclass, lon=X, lat=Y, depth=Z)
 
-    def SampleW(particle, fieldset, time):
+    def SampleW(particle, fieldset, time):  # pragma: no cover
         particle.w = fieldset.W[time, particle.depth, particle.lat, particle.lon]
 
     pset.execute([AdvectionRK4_3D, SampleW], runtime=runtime, dt=100)
@@ -279,7 +279,7 @@ def create_periodic_fieldset(xdim, ydim, uvel, vvel):
     return FieldSet.from_data(data, dimensions, mesh="spherical", transpose=True)
 
 
-def periodicBC(particle, fieldset, time):
+def periodicBC(particle, fieldset, time):  # pragma: no cover
     particle.lon = math.fmod(particle.lon, 1)
     particle.lat = math.fmod(particle.lat, 1)
 
diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py
index 9c6a93b37c..64237ec743 100644
--- a/tests/test_diffusion.py
+++ b/tests/test_diffusion.py
@@ -104,7 +104,7 @@ def test_randomexponential(mode, lambd):
         fieldset=fieldset, pclass=ptype[mode], lon=np.zeros(npart), lat=np.zeros(npart), depth=np.zeros(npart)
     )
 
-    def vertical_randomexponential(particle, fieldset, time):
+    def vertical_randomexponential(particle, fieldset, time):  # pragma: no cover
         # Kernel for random exponential variable in depth direction
         particle.depth = ParcelsRandom.expovariate(fieldset.lambd)
 
@@ -134,7 +134,7 @@ def test_randomvonmises(mode, mu, kappa):
         fieldset=fieldset, pclass=AngleParticle, lon=np.zeros(npart), lat=np.zeros(npart), depth=np.zeros(npart)
     )
 
-    def vonmises(particle, fieldset, time):
+    def vonmises(particle, fieldset, time):  # pragma: no cover
         particle.angle = ParcelsRandom.vonmisesvariate(fieldset.mu, fieldset.kappa)
 
     pset.execute(vonmises, runtime=1, dt=1)
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index cf8dbcb898..c48d822e96 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -523,7 +523,7 @@ def test_fieldset_cellareas(mesh):
             assert np.allclose(cell_areas[y, :], cell_areas[y, 0], rtol=1e-3)
 
 
-def addConst(particle, fieldset, time):
+def addConst(particle, fieldset, time):  # pragma: no cover
     particle.lon = particle.lon + fieldset.movewest + fieldset.moveeast
 
 
@@ -583,7 +583,7 @@ def test_add_second_vector_field(mode):
     UV2 = VectorField("UV2", fieldset2.U2, fieldset2.V2)
     fieldset.add_vector_field(UV2)
 
-    def SampleUV2(particle, fieldset, time):
+    def SampleUV2(particle, fieldset, time):  # pragma: no cover
         u, v = fieldset.UV2[time, particle.depth, particle.lat, particle.lon]
         particle_dlon += u * particle.dt  # noqa
         particle_dlat += v * particle.dt  # noqa
@@ -607,7 +607,7 @@ def test_fieldset_write(tmp_zarrfile):
 
     fieldset.U.to_write = True
 
-    def UpdateU(particle, fieldset, time):
+    def UpdateU(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[particle]
         fieldset.U.data[particle.ti, particle.yi, particle.xi] += 1
         fieldset.U.grid.time[0] = time
@@ -653,7 +653,7 @@ def advance(self):
     def perIterGC():
         gc.collect()
 
-    def periodicBoundaryConditions(particle, fieldset, time):
+    def periodicBoundaryConditions(particle, fieldset, time):  # pragma: no cover
         while particle.lon > 180.0:
             particle_dlon -= 360.0  # noqa
         while particle.lon < -180.0:
@@ -830,7 +830,7 @@ def temp_func(time):
             data, dimensions, mesh="flat", time_periodic=time_periodic, transpose=True, allow_time_extrapolation=True
         )
 
-    def sampleTemp(particle, fieldset, time):
+    def sampleTemp(particle, fieldset, time):  # pragma: no cover
         particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon]
         # test if we can interpolate UV and UVW together
         (particle.u1, particle.v1) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
@@ -955,7 +955,7 @@ def test_fieldset_initialisation_kernel_dask(time2, tmpdir):
         filepath, chunksize={"time": ("time_counter", 1), "depth": ("depthu", 1), "lat": ("y", 2), "lon": ("x", 2)}
     )
 
-    def SampleField(particle, fieldset, time):
+    def SampleField(particle, fieldset, time):  # pragma: no cover
         particle.u_kernel, particle.v_kernel = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
 
     SampleParticle = JITParticle.add_variables(
@@ -1114,7 +1114,7 @@ def test_deferredload_simplefield(mode, direction, time_extrapolation, tmpdir):
     SamplingParticle = ptype[mode].add_variable("p")
     pset = ParticleSet(fieldset, SamplingParticle, lon=0.5, lat=0.5)
 
-    def SampleU(particle, fieldset, time):
+    def SampleU(particle, fieldset, time):  # pragma: no cover
         particle.p, tmp = fieldset.UV[particle]
 
     runtime = tdim * 2 if time_extrapolation else None
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index ce5ed7e6f1..36b67005c2 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -30,15 +30,15 @@ def pclass(mode):
     )
 
 
-def SampleUV(particle, fieldset, time):
+def SampleUV(particle, fieldset, time):  # pragma: no cover
     (particle.u, particle.v) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
 
 
-def SampleUVNoConvert(particle, fieldset, time):
+def SampleUVNoConvert(particle, fieldset, time):  # pragma: no cover
     (particle.u, particle.v) = fieldset.UV.eval(time, particle.depth, particle.lat, particle.lon, applyConversion=False)
 
 
-def SampleP(particle, fieldset, time):
+def SampleP(particle, fieldset, time):  # pragma: no cover
     particle.p = fieldset.P[particle]
 
 
@@ -624,7 +624,7 @@ def test_sampling_3DCROCO(mode):
 
     SampleP = ptype[mode].add_variable("p", initial=0.0)
 
-    def SampleU(particle, fieldset, time):
+    def SampleU(particle, fieldset, time):  # pragma: no cover
         particle.p = fieldset.U[time, particle.depth, particle.lat, particle.lon, particle]
 
     pset = ParticleSet(fieldset, pclass=SampleP, lon=120e3, lat=50e3, depth=-0.4)
@@ -684,7 +684,7 @@ def test_sampling_multigrids_non_vectorfield_from_file(mode, npart, tmpdir, chs)
 
     pset = ParticleSet.from_line(fieldset, pclass=TestParticle, start=[0.3, 0.3], finish=[0.7, 0.7], size=npart)
 
-    def test_sample(particle, fieldset, time):
+    def test_sample(particle, fieldset, time):  # pragma: no cover
         particle.sample_var += fieldset.B[time, fieldset.sample_depth, particle.lat, particle.lon]
 
     kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample)
@@ -736,7 +736,7 @@ def test_sampling_multigrids_non_vectorfield(mode, npart):
 
     pset = ParticleSet.from_line(fieldset, pclass=TestParticle, start=[0.3, 0.3], finish=[0.7, 0.7], size=npart)
 
-    def test_sample(particle, fieldset, time):
+    def test_sample(particle, fieldset, time):  # pragma: no cover
         particle.sample_var += fieldset.B[time, fieldset.sample_depth, particle.lat, particle.lon]
 
     kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample)
@@ -865,7 +865,7 @@ def test_nestedfields(mode):
     P = NestedField("P", [P1, P2])
     fieldset.add_field(P)
 
-    def Recover(particle, fieldset, time):
+    def Recover(particle, fieldset, time):  # pragma: no cover
         if particle.state == StatusCode.ErrorOutOfBounds:
             particle_dlon = 0  # noqa
             particle_dlat = 0  # noqa
@@ -918,7 +918,7 @@ def calc_p(t, y, x):
     xv, yv = np.meshgrid(np.arange(0, 1, 0.5), np.arange(0, 1, 0.5))
     pset = ParticleSet(fieldset, pclass=pclass(mode), lon=xv.flatten(), lat=yv.flatten())
 
-    def SampleP(particle, fieldset, time):
+    def SampleP(particle, fieldset, time):  # pragma: no cover
         particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
     kernels = [AdvectionRK4, SampleP]
diff --git a/tests/test_grids.py b/tests/test_grids.py
index f26c413659..dfb31aa83f 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -76,7 +76,7 @@ def temp_func(lon, lat):
 
     fieldset = FieldSet(u_field, v_field, fields=other_fields)
 
-    def sampleTemp(particle, fieldset, time):
+    def sampleTemp(particle, fieldset, time):  # pragma: no cover
         # Note that fieldset.temp is interpolated at time=time+dt.
         # Indeed, sampleTemp is called at time=time, but the result is written
         # at time=time+dt, after the Kernel update
@@ -231,7 +231,7 @@ def bath_func(lon):
     other_fields["temp"] = temp_field
     fieldset = FieldSet(u_field, v_field, fields=other_fields)
 
-    def sampleTemp(particle, fieldset, time):
+    def sampleTemp(particle, fieldset, time):  # pragma: no cover
         particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon]
 
     MyParticle = ptype[mode].add_variable("temp", dtype=np.float32, initial=20.0)
@@ -321,7 +321,7 @@ def bath_func(lon):
 
     MyParticle = ptype[mode].add_variable("relDepth", dtype=np.float32, initial=20.0)
 
-    def moveEast(particle, fieldset, time):
+    def moveEast(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 5 * particle.dt  # noqa
         particle.relDepth = fieldset.relDepth[time, particle.depth, particle.lat, particle.lon]
 
@@ -356,7 +356,7 @@ def test_curvilinear_grids(mode):
     v_field = Field("V", v_data, grid=grid, transpose=False)
     fieldset = FieldSet(u_field, v_field)
 
-    def sampleSpeed(particle, fieldset, time):
+    def sampleSpeed(particle, fieldset, time):  # pragma: no cover
         u, v = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
         particle.speed = math.sqrt(u * u + v * v)
 
@@ -388,7 +388,7 @@ def test_nemo_grid(mode):
     # test ParticleSet.from_field on curvilinear grids
     ParticleSet.from_field(fieldset, ptype[mode], start_field=fieldset.U, size=5)
 
-    def sampleVel(particle, fieldset, time):
+    def sampleVel(particle, fieldset, time):  # pragma: no cover
         (particle.zonal, particle.meridional) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
 
     MyParticle = ptype[mode].add_variables(
@@ -449,7 +449,7 @@ def test_cgrid_uniform_2dvel(mode, time):
     fieldset.U.interp_method = "cgrid_velocity"
     fieldset.V.interp_method = "cgrid_velocity"
 
-    def sampleVel(particle, fieldset, time):
+    def sampleVel(particle, fieldset, time):  # pragma: no cover
         (particle.zonal, particle.meridional) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
 
     MyParticle = ptype[mode].add_variables(
@@ -502,7 +502,7 @@ def test_cgrid_uniform_3dvel(mode, vert_mode, time):
     fieldset.V.interp_method = "cgrid_velocity"
     fieldset.W.interp_method = "cgrid_velocity"
 
-    def sampleVel(particle, fieldset, time):
+    def sampleVel(particle, fieldset, time):  # pragma: no cover
         (particle.zonal, particle.meridional, particle.vertical) = fieldset.UVW[
             time, particle.depth, particle.lat, particle.lon
         ]
@@ -560,7 +560,7 @@ def test_cgrid_uniform_3dvel_spherical(mode, vert_mode, time):
     fieldset.V.interp_method = "cgrid_velocity"
     fieldset.W.interp_method = "cgrid_velocity"
 
-    def sampleVel(particle, fieldset, time):
+    def sampleVel(particle, fieldset, time):  # pragma: no cover
         (particle.zonal, particle.meridional, particle.vertical) = fieldset.UVW[
             time, particle.depth, particle.lat, particle.lon
         ]
@@ -601,11 +601,11 @@ def test_popgrid(mode, vert_discretisation, deferred_load):
 
     fieldset = FieldSet.from_pop(filenames, variables, dimensions, mesh="flat", deferred_load=deferred_load)
 
-    def sampleVel(particle, fieldset, time):
+    def sampleVel(particle, fieldset, time):  # pragma: no cover
         (particle.zonal, particle.meridional, particle.vert) = fieldset.UVW[particle]
         particle.tracer = fieldset.T[particle]
 
-    def OutBoundsError(particle, fieldset, time):
+    def OutBoundsError(particle, fieldset, time):  # pragma: no cover
         if particle.state == StatusCode.ErrorOutOfBounds:
             particle.out_of_bounds = 1
             particle_ddepth -= 3  # noqa
@@ -712,7 +712,7 @@ def calculate_UVR(lat, lon, dx, dy, omega, alpha):
     fieldset.U.interp_method = "cgrid_velocity"
     fieldset.V.interp_method = "cgrid_velocity"
 
-    def UpdateR(particle, fieldset, time):
+    def UpdateR(particle, fieldset, time):  # pragma: no cover
         if time == 0:
             particle.radius_start = fieldset.R[time, particle.depth, particle.lat, particle.lon]
         particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]
@@ -789,7 +789,7 @@ def populate_UVWR(lat, lon, depth, dx, dz, omega):
     fieldset.V.interp_method = "cgrid_velocity"
     fieldset.W.interp_method = "cgrid_velocity"
 
-    def UpdateR(particle, fieldset, time):
+    def UpdateR(particle, fieldset, time):  # pragma: no cover
         if time == 0:
             particle.radius_start = fieldset.R[time, particle.depth, particle.lat, particle.lon]
         particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]
@@ -867,7 +867,7 @@ def populate_UVWR(lat, lon, depth, dx, dz, omega):
     fieldset.V.interp_method = "bgrid_velocity"
     fieldset.W.interp_method = "bgrid_w_velocity"
 
-    def UpdateR(particle, fieldset, time):
+    def UpdateR(particle, fieldset, time):  # pragma: no cover
         if time == 0:
             particle.radius_start = fieldset.R[time, particle.depth, particle.lat, particle.lon]
         particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]
@@ -940,7 +940,7 @@ def test_bgrid_interpolation(gridindexingtype, mode, extrapolation):
     fieldset.U.units = UnitConverter()
     fieldset.V.units = UnitConverter()
 
-    def VelocityInterpolator(particle, fieldset, time):
+    def VelocityInterpolator(particle, fieldset, time):  # pragma: no cover
         particle.Uvel = fieldset.U[time, particle.depth, particle.lat, particle.lon]
         particle.Vvel = fieldset.V[time, particle.depth, particle.lat, particle.lon]
         particle.Wvel = fieldset.W[time, particle.depth, particle.lat, particle.lon]
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index 0aff20f6ab..bcf0e1b413 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -46,13 +46,13 @@ def test_execution_order(mode, kernel_type):
         {"U": [[0, 1], [2, 3]], "V": np.ones((2, 2))}, {"lon": [0, 2], "lat": [0, 2]}, mesh="flat"
     )
 
-    def MoveLon_Update_Lon(particle, fieldset, time):
+    def MoveLon_Update_Lon(particle, fieldset, time):  # pragma: no cover
         particle.lon += 0.2
 
-    def MoveLon_Update_dlon(particle, fieldset, time):
+    def MoveLon_Update_dlon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 0.2  # noqa
 
-    def SampleP(particle, fieldset, time):
+    def SampleP(particle, fieldset, time):  # pragma: no cover
         particle.p = fieldset.U[time, particle.depth, particle.lat, particle.lon]
 
     SampleParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0.0)
@@ -125,7 +125,7 @@ def test_execution_runtime(fieldset_unit_mesh, mode, start, end, substeps, dt):
 def test_execution_fail_python_exception(fieldset_unit_mesh, mode):
     npart = 10
 
-    def PythonFail(particle, fieldset, time):
+    def PythonFail(particle, fieldset, time):  # pragma: no cover
         if particle.time >= 10.0:
             raise RuntimeError("Enough is enough!")
         else:
@@ -145,7 +145,7 @@ def PythonFail(particle, fieldset, time):
 def test_execution_fail_out_of_bounds(fieldset_unit_mesh, mode):
     npart = 10
 
-    def MoveRight(particle, fieldset, time):
+    def MoveRight(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[time, particle.depth, particle.lat, particle.lon + 0.1, particle]
         particle_dlon += 0.1  # noqa
 
@@ -162,11 +162,11 @@ def MoveRight(particle, fieldset, time):
 def test_execution_recover_out_of_bounds(fieldset_unit_mesh, mode):
     npart = 2
 
-    def MoveRight(particle, fieldset, time):
+    def MoveRight(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[time, particle.depth, particle.lat, particle.lon + 0.1, particle]
         particle_dlon += 0.1  # noqa
 
-    def MoveLeft(particle, fieldset, time):
+    def MoveLeft(particle, fieldset, time):  # pragma: no cover
         if particle.state == StatusCode.ErrorOutOfBounds:
             particle_dlon -= 1.0  # noqa
             particle.state = StatusCode.Success
@@ -182,10 +182,10 @@ def MoveLeft(particle, fieldset, time):
 
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_execution_check_all_errors(fieldset_unit_mesh, mode):
-    def MoveRight(particle, fieldset, time):
+    def MoveRight(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[time, particle.depth, particle.lat, particle.lon, particle]
 
-    def RecoverAllErrors(particle, fieldset, time):
+    def RecoverAllErrors(particle, fieldset, time):  # pragma: no cover
         if particle.state > 4:
             particle.state = StatusCode.Delete
 
@@ -196,7 +196,7 @@ def RecoverAllErrors(particle, fieldset, time):
 
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_execution_check_stopallexecution(fieldset_unit_mesh, mode):
-    def addoneLon(particle, fieldset, time):
+    def addoneLon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 1  # noqa
 
         if particle.lon + particle_dlon >= 10:
@@ -214,7 +214,7 @@ def addoneLon(particle, fieldset, time):
 def test_execution_delete_out_of_bounds(fieldset_unit_mesh, mode):
     npart = 10
 
-    def MoveRight(particle, fieldset, time):
+    def MoveRight(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[time, particle.depth, particle.lat, particle.lon + 0.1, particle]
         particle_dlon += 0.1  # noqa
 
@@ -237,11 +237,11 @@ def test_kernel_add_no_new_variables(fieldset_unit_mesh, mode):
 def test_multi_kernel_duplicate_varnames(fieldset_unit_mesh, mode):
     # Testing for merging of two Kernels with the same variable declared
     # Should throw a warning, but go ahead regardless
-    def Kernel1(particle, fieldset, time):
+    def Kernel1(particle, fieldset, time):  # pragma: no cover
         add_lon = 0.1
         particle_dlon += add_lon  # noqa
 
-    def Kernel2(particle, fieldset, time):
+    def Kernel2(particle, fieldset, time):  # pragma: no cover
         add_lon = -0.3
         particle_dlon += add_lon  # noqa
 
@@ -254,11 +254,11 @@ def Kernel2(particle, fieldset, time):
 def test_multi_kernel_reuse_varnames(fieldset_unit_mesh, mode):
     # Testing for merging of two Kernels with the same variable declared
     # Should throw a warning, but go ahead regardless
-    def MoveEast1(particle, fieldset, time):
+    def MoveEast1(particle, fieldset, time):  # pragma: no cover
         add_lon = 0.2
         particle_dlon += add_lon  # noqa
 
-    def MoveEast2(particle, fieldset, time):
+    def MoveEast2(particle, fieldset, time):  # pragma: no cover
         particle_dlon += add_lon  # noqa
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=ptype[mode], lon=[0.5], lat=[0.5])
@@ -274,10 +274,10 @@ def test_combined_kernel_from_list(fieldset_unit_mesh):
     mixed functions and kernel objects.
     """
 
-    def MoveEast(particle, fieldset, time):
+    def MoveEast(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 0.1  # noqa
 
-    def MoveNorth(particle, fieldset, time):
+    def MoveNorth(particle, fieldset, time):  # pragma: no cover
         particle_dlat += 0.1  # noqa
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=JITParticle, lon=[0.5], lat=[0.5])
@@ -315,11 +315,11 @@ def test_combined_kernel_from_list_error_checking(fieldset_unit_mesh):
 def test_update_kernel_in_script(fieldset_unit_mesh, mode):
     # Testing what happens when kernels are updated during runtime of a script
     # Should throw a warning, but go ahead regardless
-    def MoveEast(particle, fieldset, time):
+    def MoveEast(particle, fieldset, time):  # pragma: no cover
         add_lon = 0.1
         particle_dlon += add_lon  # noqa
 
-    def MoveWest(particle, fieldset, time):
+    def MoveWest(particle, fieldset, time):  # pragma: no cover
         add_lon = -0.3
         particle_dlon += add_lon  # noqa
 
@@ -366,7 +366,7 @@ def test_explicit_ParcelsRandom(fieldset_unit_mesh, mode):
     """Testing `from parcels import ParcelsRandom` in kernel code"""
     from parcels import ParcelsRandom
 
-    def nudge_kernel(particle, fieldset, time):
+    def nudge_kernel(particle, fieldset, time):  # pragma: no cover
         dlat = ParcelsRandom.uniform(2, 3)
         particle_dlat += dlat  # noqa
 
@@ -379,7 +379,7 @@ def nudge_kernel(particle, fieldset, time):
 def test_parcels_dot_ParcelsRandom(fieldset_unit_mesh, mode):
     """Testing `parcels.ParcelsRandom` in kernel code"""
 
-    def nudge_kernel(particle, fieldset, time):
+    def nudge_kernel(particle, fieldset, time):  # pragma: no cover
         particle_dlat += parcels.ParcelsRandom.uniform(2, 3)  # noqa
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=ptype[mode], lon=[0.5], lat=[0.5])
@@ -391,7 +391,7 @@ def nudge_kernel(particle, fieldset, time):
 def test_parcels_dot_rng(fieldset_unit_mesh, mode):
     """Testing `parcels.rng` in kernel code."""
 
-    def nudge_kernel(particle, fieldset, time):
+    def nudge_kernel(particle, fieldset, time):  # pragma: no cover
         dlat = parcels.rng.uniform(2, 3)
         particle_dlat += dlat  # noqa
 
@@ -405,7 +405,7 @@ def test_custom_ParcelsRandom_alias(fieldset_unit_mesh, mode):
     """Testing aliasing ParcelsRandom to another name."""
     from parcels import ParcelsRandom as my_custom_name
 
-    def nudge_kernel(particle, fieldset, time):
+    def nudge_kernel(particle, fieldset, time):  # pragma: no cover
         particle_dlat += my_custom_name.uniform(2, 3)  # noqa
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=ptype[mode], lon=[0.5], lat=[0.5])
@@ -429,7 +429,7 @@ def test_outdated_kernel(fieldset_unit_mesh):
     """
     pset = ParticleSet(fieldset_unit_mesh, pclass=JITParticle, lon=[0.5], lat=[0.5])
 
-    def outdated_kernel(particle, fieldset, time, dt):
+    def outdated_kernel(particle, fieldset, time, dt):  # pragma: no cover
         particle.lon += 0.1
 
     with pytest.raises(ValueError) as e:
diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py
index f2d2861dab..863a97728c 100644
--- a/tests/test_kernel_language.py
+++ b/tests/test_kernel_language.py
@@ -130,7 +130,7 @@ def test_while_if_break(mode):
     TestParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0)
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=[0], lat=[0])
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         while particle.p < 30:
             if particle.p > 9:
                 break
@@ -150,7 +150,7 @@ def test_nested_if(mode):
     )
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=0, lat=0)
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         if particle.p1 >= particle.p0:
             var = particle.p0
             if var + 1 < particle.p1:
@@ -166,7 +166,7 @@ def test_pass(mode):
     TestParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0)
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=0, lat=0)
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         particle.p = -1
         pass
 
@@ -178,7 +178,7 @@ def kernel(particle, fieldset, time):
 def test_dt_as_variable_in_kernel(mode):
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=ptype[mode], lon=0, lat=0)
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         dt = 1.0  # noqa
 
     pset.execute(kernel, endtime=10, dt=1.0)
@@ -188,10 +188,10 @@ def test_parcels_tmpvar_in_kernel():
     """Tests for error thrown if variable with 'tmp' defined in custom kernel."""
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=JITParticle, lon=0, lat=0)
 
-    def kernel_tmpvar(particle, fieldset, time):
+    def kernel_tmpvar(particle, fieldset, time):  # pragma: no cover
         parcels_tmpvar0 = 0  # noqa
 
-    def kernel_pnum(particle, fieldset, time):
+    def kernel_pnum(particle, fieldset, time):  # pragma: no cover
         pnum = 0  # noqa
 
     for kernel in [kernel_tmpvar, kernel_pnum]:
@@ -207,13 +207,13 @@ def test_varname_as_fieldname():
     Particle = JITParticle.add_variable("speed")
     pset = ParticleSet(fset, pclass=Particle, lon=0, lat=0)
 
-    def kernel_particlename(particle, fieldset, time):
+    def kernel_particlename(particle, fieldset, time):  # pragma: no cover
         particle.speed = fieldset.speed[particle]
 
     pset.execute(kernel_particlename, endtime=1, dt=1.0)
     assert pset[0].speed == 10
 
-    def kernel_varname(particle, fieldset, time):
+    def kernel_varname(particle, fieldset, time):  # pragma: no cover
         vertical_speed = fieldset.vertical_speed  # noqa
 
     with pytest.raises(NotImplementedError):
@@ -224,7 +224,7 @@ def test_abs():
     """Tests for error thrown if using abs in kernel."""
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=JITParticle, lon=0, lat=0)
 
-    def kernel_abs(particle, fieldset, time):
+    def kernel_abs(particle, fieldset, time):  # pragma: no cover
         particle.lon = abs(3.1)
 
     with pytest.raises(NotImplementedError):
@@ -237,7 +237,7 @@ def test_if_withfield(fieldset_unit_mesh, mode):
     TestParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0)
     pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=[0], lat=[0])
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         u, v = fieldset.UV[time, 0, 0, 1.0]
         particle.p = 0
         if fieldset.U[time, 0, 0, 1.0] == u:
@@ -273,7 +273,7 @@ def test_print(fieldset_unit_mesh, mode, capfd):
     TestParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0)
     pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=[0.5], lat=[0.5])
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         particle.p = 1e-3
         tmp = 5
         print(f"{particle.id} {particle.p:f} {tmp:f}")
@@ -286,7 +286,7 @@ def kernel(particle, fieldset, time):
         abs(float(lst[0]) - pset.id[0]) < tol and abs(float(lst[1]) - pset.p[0]) < tol and abs(float(lst[2]) - 5) < tol
     )
 
-    def kernel2(particle, fieldset, time):
+    def kernel2(particle, fieldset, time):  # pragma: no cover
         tmp = 3
         print(f"{tmp:f}")
 
@@ -302,7 +302,7 @@ def kernel2(particle, fieldset, time):
 def test_fieldset_access(fieldset_unit_mesh, expectation, mode):
     pset = ParticleSet(fieldset_unit_mesh, pclass=ptype[mode], lon=0, lat=0)
 
-    def kernel(particle, fieldset, time):
+    def kernel(particle, fieldset, time):  # pragma: no cover
         particle.lon = fieldset.U.grid.lon[2]
 
     with expectation:
@@ -352,10 +352,10 @@ def test_random_kernel_concat(fieldset_unit_mesh, mode, concat):
     TestParticle = ptype[mode].add_variable("p", dtype=np.float32, initial=0)
     pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=0, lat=0)
 
-    def RandomKernel(particle, fieldset, time):
+    def RandomKernel(particle, fieldset, time):  # pragma: no cover
         particle.p += ParcelsRandom.uniform(0, 1)
 
-    def AddOne(particle, fieldset, time):
+    def AddOne(particle, fieldset, time):  # pragma: no cover
         particle.p += 1.0
 
     kernels = [RandomKernel, AddOne] if concat else RandomKernel
@@ -389,7 +389,7 @@ def func(U, lon, dt):
     else:
         c_include = os.path.join(os.path.dirname(__file__), "customed_header.h")
 
-    def ckernel(particle, fieldset, time):
+    def ckernel(particle, fieldset, time):  # pragma: no cover
         func("parcels_customed_Cfunc_pointer_args", fieldset.U, particle_dlon, particle.dt)  # noqa
 
     kernel = pset.Kernel(ckernel, c_include=c_include)
@@ -402,7 +402,7 @@ def test_dt_modif_by_kernel(mode):
     TestParticle = ptype[mode].add_variable("age", dtype=np.float32, initial=0)
     pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=[0.5], lat=[0])
 
-    def modif_dt(particle, fieldset, time):
+    def modif_dt(particle, fieldset, time):  # pragma: no cover
         particle.age += particle.dt
         particle.dt = 2
 
diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py
index cb3501eaa2..de04bde743 100755
--- a/tests/test_particlefile.py
+++ b/tests/test_particlefile.py
@@ -80,7 +80,7 @@ def test_pfile_set_towrite_False(fieldset, mode, tmp_zarrfile):
     pset.set_variable_write_status("lat", False)
     pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1)
 
-    def Update_lon(particle, fieldset, time):
+    def Update_lon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 0.1  # noqa
 
     pset.execute(Update_lon, runtime=10, output_file=pfile)
@@ -109,7 +109,7 @@ def test_pfile_array_remove_all_particles(fieldset, mode, chunks_obs, tmp_zarrfi
     pfile.write(pset, 1)
     pfile.write(pset, 2)
 
-    ds = xr.open_zarr(tmp_zarrfile)
+    ds = xr.open_zarr(tmp_zarrfile).load()
     assert np.allclose(ds["time"][:, 0], np.timedelta64(0, "s"), atol=np.timedelta64(1, "ms"))
     if chunks_obs is not None:
         assert ds["time"][:].shape == chunks
@@ -121,7 +121,7 @@ def test_pfile_array_remove_all_particles(fieldset, mode, chunks_obs, tmp_zarrfi
 
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_variable_write_double(fieldset, mode, tmp_zarrfile):
-    def Update_lon(particle, fieldset, time):
+    def Update_lon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 0.1  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
@@ -157,7 +157,7 @@ def test_write_dtypes_pfile(fieldset, mode, tmp_zarrfile):
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 @pytest.mark.parametrize("npart", [1, 2, 5])
 def test_variable_written_once(fieldset, mode, tmp_zarrfile, npart):
-    def Update_v(particle, fieldset, time):
+    def Update_v(particle, fieldset, time):  # pragma: no cover
         particle.v_once += 1.0
         particle.age += particle.dt
 
@@ -203,7 +203,7 @@ def test_pset_repeated_release_delayed_adding_deleting(type, fieldset, mode, rep
         )
     pfile = pset.ParticleFile(tmp_zarrfile, outputdt=abs(dt), chunks=(1, 1))
 
-    def IncrLon(particle, fieldset, time):
+    def IncrLon(particle, fieldset, time):  # pragma: no cover
         particle.sample_var += 1.0
         if particle.sample_var > fieldset.maxvar:
             particle.delete()
@@ -237,7 +237,7 @@ def test_pfile_chunks_repeatedrelease(fieldset, mode, repeatdt, nump, tmp_zarrfi
     chunks = (20, 10)
     pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1, chunks=chunks)
 
-    def DoNothing(particle, fieldset, time):
+    def DoNothing(particle, fieldset, time):  # pragma: no cover
         pass
 
     pset.execute(DoNothing, dt=1, runtime=runtime, output_file=pfile)
@@ -247,7 +247,7 @@ def DoNothing(particle, fieldset, time):
 
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_write_timebackward(fieldset, mode, tmp_zarrfile):
-    def Update_lon(particle, fieldset, time):
+    def Update_lon(particle, fieldset, time):  # pragma: no cover
         particle_dlon -= 0.1 * particle.dt  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lat=np.linspace(0, 1, 3), lon=[0, 0, 0], time=[1, 2, 3])
@@ -274,7 +274,7 @@ def test_write_xiyi(fieldset, mode, tmp_zarrfile):
         ]
     )
 
-    def Get_XiYi(particle, fieldset, time):
+    def Get_XiYi(particle, fieldset, time):  # 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.
@@ -284,7 +284,7 @@ def Get_XiYi(particle, fieldset, time):
         particle.pxi1 = particle.xi[1]
         particle.pyi = particle.yi[0]
 
-    def SampleP(particle, fieldset, time):
+    def SampleP(particle, fieldset, time):  # pragma: no cover
         if time > 5 * 3600:
             _ = fieldset.P[particle]  # To trigger sampling of the P field
 
@@ -324,7 +324,7 @@ def test_reset_dt(fieldset, mode, 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):
+    def Update_lon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 0.1  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
@@ -338,7 +338,7 @@ def Update_lon(particle, fieldset, time):
 def test_correct_misaligned_outputdt_dt(fieldset, mode, tmp_zarrfile):
     """Testing that outputdt does not need to be a multiple of dt."""
 
-    def Update_lon(particle, fieldset, time):
+    def Update_lon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += particle.dt  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
diff --git a/tests/test_particles.py b/tests/test_particles.py
index 9ba8f8a68d..6a8474ba41 100644
--- a/tests/test_particles.py
+++ b/tests/test_particles.py
@@ -39,7 +39,7 @@ def test_variable_init(fieldset, mode):
     TestParticle = TestParticle.add_variable("p_int", np.int32, initial=12.0)
     pset = ParticleSet(fieldset, pclass=TestParticle, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart))
 
-    def addOne(particle, fieldset, time):
+    def addOne(particle, fieldset, time):  # pragma: no cover
         particle.p_float += 1.0
         particle.p_double += 1.0
         particle.p_int += 1
diff --git a/tests/test_particlesets.py b/tests/test_particlesets.py
index 0678df36a9..2237c262f2 100644
--- a/tests/test_particlesets.py
+++ b/tests/test_particlesets.py
@@ -91,7 +91,7 @@ def test_pset_create_fromparticlefile(fieldset, mode, restart, tmp_zarrfile):
     pset = ParticleSet(fieldset, lon=lon, lat=lat, depth=[4] * len(lon), pclass=TestParticle, p3=np.arange(len(lon)))
     pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1)
 
-    def Kernel(particle, fieldset, time):
+    def Kernel(particle, fieldset, time):  # pragma: no cover
         particle.p = 2.0
         if particle.lon == 1.0:
             particle.delete()
@@ -189,7 +189,7 @@ def test_pset_not_multipldt_time(fieldset, mode):
     times = [0, 1.1]
     pset = ParticleSet(fieldset, lon=[0] * 2, lat=[0] * 2, pclass=ptype[mode], time=times)
 
-    def Addlon(particle, fieldset, time):
+    def Addlon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += particle.dt  # noqa
 
     pset.execute(Addlon, dt=1, runtime=2)
@@ -203,7 +203,7 @@ def test_pset_repeated_release(fieldset, mode):
     pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), pclass=ptype[mode], time=time)
     assert np.allclose([p.time for p in pset], time)
 
-    def IncrLon(particle, fieldset, time):
+    def IncrLon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 1.0  # noqa
 
     pset.execute(IncrLon, dt=1.0, runtime=npart + 1)
@@ -213,7 +213,7 @@ def IncrLon(particle, fieldset, time):
 def test_pset_repeatdt_check_dt(fieldset):
     pset = ParticleSet(fieldset, lon=[0], lat=[0], pclass=ScipyParticle, repeatdt=5)
 
-    def IncrLon(particle, fieldset, time):
+    def IncrLon(particle, fieldset, time):  # pragma: no cover
         particle.lon = 1.0
 
     pset.execute(IncrLon, dt=2, runtime=21)
@@ -234,7 +234,7 @@ def test_pset_repeatdt_custominit(fieldset, mode):
 def test_pset_stop_simulation(fieldset, mode):
     pset = ParticleSet(fieldset, lon=0, lat=0, pclass=ptype[mode])
 
-    def Delete(particle, fieldset, time):
+    def Delete(particle, fieldset, time):  # pragma: no cover
         if time == 4:
             return StatusCode.StopExecution
 
@@ -299,7 +299,7 @@ def test_pset_add_shorthand(fieldset, mode):
 def test_pset_add_execute(fieldset, mode):
     npart = 10
 
-    def AddLat(particle, fieldset, time):
+    def AddLat(particle, fieldset, time):  # pragma: no cover
         particle_dlat += 0.1  # noqa
 
     pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode])
@@ -364,7 +364,7 @@ def test_pset_remove_particle(fieldset, mode):
 def test_pset_remove_kernel(fieldset, mode):
     npart = 100
 
-    def DeleteKernel(particle, fieldset, time):
+    def DeleteKernel(particle, fieldset, time):  # pragma: no cover
         if particle.lon >= 0.4:
             particle.delete()
 
@@ -378,7 +378,7 @@ def test_pset_multi_execute(fieldset, mode):
     npart = 10
     n = 5
 
-    def AddLat(particle, fieldset, time):
+    def AddLat(particle, fieldset, time):  # pragma: no cover
         particle_dlat += 0.1  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=np.linspace(0, 1, npart), lat=np.zeros(npart))
@@ -393,7 +393,7 @@ def test_pset_multi_execute_delete(fieldset, mode):
     npart = 10
     n = 5
 
-    def AddLat(particle, fieldset, time):
+    def AddLat(particle, fieldset, time):  # pragma: no cover
         particle_dlat += 0.1  # noqa
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=np.linspace(0, 1, npart), lat=np.zeros(npart))
@@ -442,7 +442,7 @@ def test_from_field_exact_val(staggered_grid):
 
     SampleParticle = ptype["scipy"].add_variable("mask", initial=0)
 
-    def SampleMask(particle, fieldset, time):
+    def SampleMask(particle, fieldset, time):  # pragma: no cover
         particle.mask = fieldset.mask[particle]
 
     pset = ParticleSet.from_field(fieldset, size=400, pclass=SampleParticle, start_field=FMask, time=0)