diff --git a/tests/test_advection.py b/tests/test_advection.py
index b926c25eb..3395d2a16 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels import (
     AdvectionAnalytical,
@@ -98,7 +99,7 @@ def test_advection_meridional(lon, lat):
     pset = ParticleSet(fieldset, pclass=Particle, lon=np.linspace(-60, 60, npart), lat=np.linspace(0, 30, npart))
     delta_lat = np.diff(pset.lat)
     pset.execute(AdvectionRK4, runtime=timedelta(hours=2), dt=timedelta(seconds=30))
-    assert np.allclose(np.diff(pset.lat), delta_lat, rtol=1.0e-4)
+    assert_allclose(np.diff(pset.lat), delta_lat, rtol=1.0e-4)
 
 
 def test_advection_3D():
@@ -119,7 +120,7 @@ def test_advection_3D():
     )
     time = timedelta(hours=2).total_seconds()
     pset.execute(AdvectionRK4, runtime=time, dt=timedelta(seconds=30))
-    assert np.allclose(pset.depth * pset.time, pset.lon, atol=1.0e-1)
+    assert_allclose(pset.depth * pset.time, pset.lon, atol=1.0e-1)
 
 
 @pytest.mark.parametrize("direction", ["up", "down"])
@@ -161,8 +162,8 @@ def SubmergeParticle(particle, fieldset, time):  # pragma: no cover
     pset.execute(kernels, runtime=11.0, dt=1)
 
     if direction == "up" and wErrorThroughSurface:
-        assert np.allclose(pset.lon[0], 0.6)
-        assert np.allclose(pset.depth[0], 0)
+        assert_allclose(pset.lon[0], 0.6)
+        assert_allclose(pset.depth[0], 0)
     else:
         assert len(pset) == 0
 
@@ -227,7 +228,7 @@ def test_conversion_3DCROCO():
     for zi, z in enumerate(z_xroms):
         sigma[zi] = _croco_from_z_to_sigma_scipy(fieldset, 0, z, lat, lon, None)
 
-    assert np.allclose(sigma, s_xroms, atol=1e-3)
+    assert_allclose(sigma, s_xroms, atol=1e-3)
 
 
 @pytest.mark.v4alpha
@@ -246,8 +247,8 @@ 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)
-    assert np.allclose(pset.depth, Z.flatten(), atol=5)  # TODO lower this atol
-    assert np.allclose(pset.lon_nextloop, [x + runtime for x in X.flatten()], atol=1e-3)
+    assert_allclose(pset.depth, Z.flatten(), atol=5)  # TODO lower this atol
+    assert_allclose(pset.lon_nextloop, [x + runtime for x in X.flatten()], atol=1e-3)
 
 
 @pytest.mark.v4alpha
@@ -262,8 +263,8 @@ def test_advection_2DCROCO():
     pset = ParticleSet(fieldset=fieldset, pclass=Particle, lon=X, lat=Y, depth=Z)
 
     pset.execute([AdvectionRK4], runtime=runtime, dt=100)
-    assert np.allclose(pset.depth, Z.flatten(), atol=1e-3)
-    assert np.allclose(pset.lon_nextloop, [x + runtime for x in X], atol=1e-3)
+    assert_allclose(pset.depth, Z.flatten(), atol=1e-3)
+    assert_allclose(pset.lon_nextloop, [x + runtime for x in X], atol=1e-3)
 
 
 def create_periodic_fieldset(xdim, ydim, uvel, vvel):
@@ -315,8 +316,8 @@ def test_advection_periodic_zonal_meridional():
     fieldset.add_periodic_halo(zonal=True, meridional=True)
     assert len(fieldset.U.lat) == ydim + 10  # default halo size is 5 grid points
     assert len(fieldset.U.lon) == xdim + 10  # default halo size is 5 grid points
-    assert np.allclose(np.diff(fieldset.U.lat), fieldset.U.lat[1] - fieldset.U.lat[0], rtol=0.001)
-    assert np.allclose(np.diff(fieldset.U.lon), fieldset.U.lon[1] - fieldset.U.lon[0], rtol=0.001)
+    assert_allclose(np.diff(fieldset.U.lat), fieldset.U.lat[1] - fieldset.U.lat[0], rtol=0.001)
+    assert_allclose(np.diff(fieldset.U.lon), fieldset.U.lon[1] - fieldset.U.lon[0], rtol=0.001)
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=[0.4], lat=[0.5])
     pset.execute(AdvectionRK4 + pset.Kernel(periodicBC), runtime=timedelta(hours=20), dt=timedelta(seconds=30))
@@ -431,8 +432,8 @@ def test_stationary_eddy(fieldset_stationary, method, rtol, diffField):
 
     exp_lon = [truth_stationary(x, y, pset[0].time)[0] for x, y in zip(lon, lat, strict=True)]
     exp_lat = [truth_stationary(x, y, pset[0].time)[1] for x, y in zip(lon, lat, strict=True)]
-    assert np.allclose(pset.lon, exp_lon, rtol=rtol)
-    assert np.allclose(pset.lat, exp_lat, rtol=rtol)
+    assert_allclose(pset.lon, exp_lon, rtol=rtol)
+    assert_allclose(pset.lat, exp_lat, rtol=rtol)
 
 
 def test_stationary_eddy_vertical():
@@ -460,9 +461,9 @@ def test_stationary_eddy_vertical():
     exp_lon = [truth_stationary(x, z, pset[0].time)[0] for x, z in zip(lon, depth, strict=True)]
     exp_depth = [truth_stationary(x, z, pset[0].time)[1] for x, z in zip(lon, depth, strict=True)]
     print(pset, exp_lon)
-    assert np.allclose(pset.lon, exp_lon, rtol=1e-5)
-    assert np.allclose(pset.lat, lat, rtol=1e-5)
-    assert np.allclose(pset.depth, exp_depth, rtol=1e-5)
+    assert_allclose(pset.lon, exp_lon, rtol=1e-5)
+    assert_allclose(pset.lat, lat, rtol=1e-5)
+    assert_allclose(pset.depth, exp_depth, rtol=1e-5)
 
     data = {"U": fldzero, "V": fld2, "W": fld1}
     fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
@@ -471,9 +472,9 @@ def test_stationary_eddy_vertical():
     pset.execute(AdvectionRK4_3D, dt=dt, endtime=endtime)
     exp_depth = [truth_stationary(z, y, pset[0].time)[0] for z, y in zip(depth, lat, strict=True)]
     exp_lat = [truth_stationary(z, y, pset[0].time)[1] for z, y in zip(depth, lat, strict=True)]
-    assert np.allclose(pset.lon, lon, rtol=1e-5)
-    assert np.allclose(pset.lat, exp_lat, rtol=1e-5)
-    assert np.allclose(pset.depth, exp_depth, rtol=1e-5)
+    assert_allclose(pset.lon, lon, rtol=1e-5)
+    assert_allclose(pset.lat, exp_lat, rtol=1e-5)
+    assert_allclose(pset.depth, exp_depth, rtol=1e-5)
 
 
 def truth_moving(x_0, y_0, t):
@@ -536,8 +537,8 @@ def test_moving_eddy(fieldset_moving, method, rtol, diffField):
 
     exp_lon = [truth_moving(x, y, t)[0] for x, y, t in zip(lon, lat, pset.time, strict=True)]
     exp_lat = [truth_moving(x, y, t)[1] for x, y, t in zip(lon, lat, pset.time, strict=True)]
-    assert np.allclose(pset.lon, exp_lon, rtol=rtol)
-    assert np.allclose(pset.lat, exp_lat, rtol=rtol)
+    assert_allclose(pset.lon, exp_lon, rtol=rtol)
+    assert_allclose(pset.lat, exp_lat, rtol=rtol)
 
 
 def truth_decaying(x_0, y_0, t):
@@ -619,8 +620,8 @@ def test_decaying_eddy(fieldset_decaying, method, rtol, diffField):
 
     exp_lon = [truth_decaying(x, y, t)[0] for x, y, t in zip(lon, lat, pset.time, strict=True)]
     exp_lat = [truth_decaying(x, y, t)[1] for x, y, t in zip(lon, lat, pset.time, strict=True)]
-    assert np.allclose(pset.lon, exp_lon, rtol=rtol)
-    assert np.allclose(pset.lat, exp_lat, rtol=rtol)
+    assert_allclose(pset.lon, exp_lon, rtol=rtol)
+    assert_allclose(pset.lat, exp_lat, rtol=rtol)
 
 
 def test_analyticalAgrid():
@@ -669,6 +670,6 @@ def test_uniform_analytical(u, v, w, direction, tmp_zarrfile):
     ds = xr.open_zarr(tmp_zarrfile)
     times = (direction * ds["time"][:]).values.astype("timedelta64[s]")[0]
     timeref = np.arange(1, 5).astype("timedelta64[s]")
-    assert np.allclose(times, timeref, atol=np.timedelta64(1, "ms"))
+    assert_allclose(times, timeref, atol=np.timedelta64(1, "ms"))
     lons = ds["lon"][:].values
-    assert np.allclose(lons, x0 + direction * u * np.arange(1, 5))
+    assert_allclose(lons, x0 + direction * u * np.arange(1, 5))
diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py
index 389174c75..99df498bf 100644
--- a/tests/test_diffusion.py
+++ b/tests/test_diffusion.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 from scipy import stats
 
 from parcels import (
@@ -44,10 +45,10 @@ def test_fieldKh_Brownian(mesh):
     lons = pset.lon
 
     tol = 500 * mesh_conversion  # effectively 500 m errors
-    assert np.allclose(np.std(lats), expected_std_lat, atol=tol)
-    assert np.allclose(np.std(lons), expected_std_lon, atol=tol)
-    assert np.allclose(np.mean(lons), 0, atol=tol)
-    assert np.allclose(np.mean(lats), 0, atol=tol)
+    assert_allclose(np.std(lats), expected_std_lat, atol=tol)
+    assert_allclose(np.std(lons), expected_std_lon, atol=tol)
+    assert_allclose(np.mean(lons), 0, atol=tol)
+    assert_allclose(np.mean(lats), 0, atol=tol)
 
 
 @pytest.mark.parametrize("mesh", ["spherical", "flat"])
@@ -78,8 +79,8 @@ def test_fieldKh_SpatiallyVaryingDiffusion(mesh, kernel):
     lats = pset.lat
     lons = pset.lon
     tol = 2000 * mesh_conversion  # effectively 2000 m errors (because of low numbers of particles)
-    assert np.allclose(np.mean(lons), 0, atol=tol)
-    assert np.allclose(np.mean(lats), 0, atol=tol)
+    assert_allclose(np.mean(lons), 0, atol=tol)
+    assert_allclose(np.mean(lats), 0, atol=tol)
     assert stats.skew(lons) > stats.skew(lats)
 
 
@@ -106,7 +107,7 @@ def vertical_randomexponential(particle, fieldset, time):  # pragma: no cover
 
     depth = pset.depth
     expected_mean = 1.0 / fieldset.lambd
-    assert np.allclose(np.mean(depth), expected_mean, rtol=0.1)
+    assert_allclose(np.mean(depth), expected_mean, rtol=0.1)
 
 
 @pytest.mark.parametrize("mu", [0.8 * np.pi, np.pi])
@@ -134,8 +135,8 @@ def vonmises(particle, fieldset, time):  # pragma: no cover
 
     angles = np.array([p.angle for p in pset])
 
-    assert np.allclose(np.mean(angles), mu, atol=0.1)
+    assert_allclose(np.mean(angles), mu, atol=0.1)
     vonmises_mean = stats.vonmises.mean(kappa=kappa, loc=mu)
-    assert np.allclose(np.mean(angles), vonmises_mean, atol=0.1)
+    assert_allclose(np.mean(angles), vonmises_mean, atol=0.1)
     vonmises_var = stats.vonmises.var(kappa=kappa, loc=mu)
-    assert np.allclose(np.var(angles), vonmises_var, atol=0.1)
+    assert_allclose(np.var(angles), vonmises_var, atol=0.1)
diff --git a/tests/test_field.py b/tests/test_field.py
index 7d7ec1ce3..ad5ef4266 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -2,6 +2,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels import Field
 from parcels.tools.converters import (
@@ -22,8 +23,8 @@ def test_field_from_netcdf_variables():
     variable = {"U": "vozocrtx"}
     f3 = Field.from_netcdf(filename, variable, dims)
 
-    assert np.allclose(f1.data, f2.data, atol=1e-12)
-    assert np.allclose(f1.data, f3.data, atol=1e-12)
+    assert_allclose(f1.data, f2.data, atol=1e-12)
+    assert_allclose(f1.data, f3.data, atol=1e-12)
 
     with pytest.raises(AssertionError):
         variable = {"U": "vozocrtx", "nav_lat": "nav_lat"}  # multiple variables will fail
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index ceb99d162..c69bcf72f 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels import (
     AdvectionRK4,
@@ -89,8 +90,8 @@ def test_fieldset_from_data(xdim, ydim):
     assert fieldset.U._creation_log == "from_data"
     assert len(fieldset.U.data.shape) == 3
     assert len(fieldset.V.data.shape) == 3
-    assert np.allclose(fieldset.U.data[0, :], data["U"], rtol=1e-12)
-    assert np.allclose(fieldset.V.data[0, :], data["V"], rtol=1e-12)
+    assert_allclose(fieldset.U.data[0, :], data["U"], rtol=1e-12)
+    assert_allclose(fieldset.V.data[0, :], data["V"], rtol=1e-12)
 
 
 def test_fieldset_extra_syntax():
@@ -147,9 +148,9 @@ def test_fieldset_from_data_different_dimensions(xdim, ydim):
     assert len(fieldset.V.data.shape) == 3
     assert len(fieldset.P.data.shape) == 4
     assert fieldset.P.data.shape == (tdim, zdim, ydim / 2, xdim / 2)
-    assert np.allclose(fieldset.U.data, 0.0, rtol=1e-12)
-    assert np.allclose(fieldset.V.data, 1.0, rtol=1e-12)
-    assert np.allclose(fieldset.P.data, 2.0, rtol=1e-12)
+    assert_allclose(fieldset.U.data, 0.0, rtol=1e-12)
+    assert_allclose(fieldset.V.data, 1.0, rtol=1e-12)
+    assert_allclose(fieldset.P.data, 2.0, rtol=1e-12)
 
 
 @pytest.mark.parametrize("xdim", [100, 200])
@@ -163,8 +164,8 @@ def test_fieldset_from_parcels(xdim, ydim, tmpdir):
     fieldset = FieldSet.from_parcels(filepath)
     assert len(fieldset.U.data.shape) == 3  # Will be 4 once we use depth
     assert len(fieldset.V.data.shape) == 3
-    assert np.allclose(fieldset.U.data[0, :], data["U"], rtol=1e-12)
-    assert np.allclose(fieldset.V.data[0, :], data["V"], rtol=1e-12)
+    assert_allclose(fieldset.U.data[0, :], data["U"], rtol=1e-12)
+    assert_allclose(fieldset.V.data[0, :], data["V"], rtol=1e-12)
 
 
 def test_fieldset_from_modulefile():
@@ -392,7 +393,7 @@ def test_fieldset_write_curvilinear(tmpdir):
     assert fieldset2.dx._creation_log == "from_netcdf"
 
     for var in ["lon", "lat", "data"]:
-        assert np.allclose(getattr(fieldset2.dx, var), getattr(fieldset.dx, var))
+        assert_allclose(getattr(fieldset2.dx, var), getattr(fieldset.dx, var))
 
 
 def addConst(particle, fieldset, time):  # pragma: no cover
@@ -490,7 +491,7 @@ def UpdateU(particle, fieldset, time):  # pragma: no cover
     assert fieldset.U.data[0, 1, 0] == 11
 
     da = xr.open_dataset(str(tmp_zarrfile).replace(".zarr", "_0005U.nc"))
-    assert np.allclose(fieldset.U.data, da["U"].values, atol=1.0)
+    assert_allclose(fieldset.U.data, da["U"].values, atol=1.0)
 
 
 @pytest.mark.v4remove
@@ -517,12 +518,12 @@ def test_timestamps(datetype, tmpdir):
     fieldset3 = FieldSet.from_parcels(tmpdir.join("file*"))
     timestamps = [dims1["time"], dims2["time"]]
     fieldset4 = FieldSet.from_parcels(tmpdir.join("file*"), timestamps=timestamps)
-    assert np.allclose(fieldset3.U.grid.time_full, fieldset4.U.grid.time_full)
+    assert_allclose(fieldset3.U.grid.time_full, fieldset4.U.grid.time_full)
 
     for d in [0, 8, 10, 12]:
         fieldset3.computeTimeChunk(d * 86400.0, 1.0)
         fieldset4.computeTimeChunk(d * 86400.0, 1.0)
-        assert np.allclose(fieldset3.U.data, fieldset4.U.data)
+        assert_allclose(fieldset3.U.data, fieldset4.U.data)
 
 
 @pytest.mark.v4remove
@@ -610,10 +611,10 @@ def sampleTemp(particle, fieldset, time):  # pragma: no cover
         temp_theo = temp_vec[-1]
     elif dt_sign == -1:
         temp_theo = temp_vec[0]
-    assert np.allclose(temp_theo, pset.temp[0], atol=1e-5)
-    assert np.allclose(pset.u1[0], pset.u2[0])
-    assert np.allclose(pset.v1[0], pset.v2[0])
-    assert np.allclose(pset.d[0], 1.0)
+    assert_allclose(temp_theo, pset.temp[0], atol=1e-5)
+    assert_allclose(pset.u1[0], pset.u2[0])
+    assert_allclose(pset.v1[0], pset.v2[0])
+    assert_allclose(pset.d[0], 1.0)
 
 
 @pytest.mark.parametrize("tdim", [10, None])
@@ -654,9 +655,9 @@ def generate_dataset(xdim, ydim, zdim=1, tdim=1):
 
     pset.execute(AdvectionRK4, dt=1, runtime=10)
     if tdim == 10:
-        assert np.allclose(pset.lon_nextloop[0], 4.5) and np.allclose(pset.lat_nextloop[0], 10)
+        assert_allclose(pset.lon_nextloop[0], 4.5) and np.allclose(pset.lat_nextloop[0], 10)
     else:
-        assert np.allclose(pset.lon_nextloop[0], 5.0) and np.allclose(pset.lat_nextloop[0], 10)
+        assert_allclose(pset.lon_nextloop[0], 5.0) and np.allclose(pset.lat_nextloop[0], 10)
 
 
 @pytest.mark.v4alpha
@@ -695,16 +696,16 @@ def test_fieldset_from_data_gridtypes():
     plon = pset.lon
     plat = pset.lat
     # sol of  dx/dt = (init_depth+1)*x+0.1; x(0)=0
-    assert np.allclose(plon, [0.17173462592827032, 0.2177736932123214])
-    assert np.allclose(plat, [1, 1])
+    assert_allclose(plon, [0.17173462592827032, 0.2177736932123214])
+    assert_allclose(plat, [1, 1])
 
     # Rectilinear S grid
     dimensions["depth"] = depth_s
     fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
     pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4])
     pset.execute(AdvectionRK4, runtime=1.5, dt=0.5)
-    assert np.allclose(plon, pset.lon)
-    assert np.allclose(plat, pset.lat)
+    assert_allclose(plon, pset.lon)
+    assert_allclose(plat, pset.lat)
 
     # Curvilinear Z grid
     dimensions["lon"] = lonm
@@ -713,8 +714,8 @@ def test_fieldset_from_data_gridtypes():
     fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
     pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4])
     pset.execute(AdvectionRK4, runtime=1.5, dt=0.5)
-    assert np.allclose(plon, pset.lon)
-    assert np.allclose(plat, pset.lat)
+    assert_allclose(plon, pset.lon)
+    assert_allclose(plat, pset.lat)
 
     # Curvilinear S grid
     dimensions["depth"] = depth_s
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index c8033a658..8108f03cf 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels import (
     AdvectionRK4,
@@ -92,10 +93,10 @@ def test_fieldset_sample(fieldset):
     lat = np.linspace(-80, 80, ydim, dtype=np.float32)
     v_s = np.array([fieldset.UV[0, 0.0, 70.0, x][1] for x in lon])
     u_s = np.array([fieldset.UV[0, 0.0, y, -45.0][0] for y in lat])
-    assert np.allclose(
+    assert_allclose(
         v_s, lon, rtol=1e-5
     )  # Tolerances were rtol=1e-7, increased due to numpy v2 float32 changes (see #1603)
-    assert np.allclose(u_s, lat, rtol=1e-5)
+    assert_allclose(u_s, lat, rtol=1e-5)
 
 
 def test_fieldset_sample_eval(fieldset):
@@ -105,10 +106,10 @@ def test_fieldset_sample_eval(fieldset):
     lat = np.linspace(-80, 80, ydim, dtype=np.float32)
     v_s = np.array([fieldset.UV.eval(0, 0.0, 70.0, x)[1] for x in lon])
     u_s = np.array([fieldset.UV.eval(0, 0.0, y, 0.0)[0] for y in lat])
-    assert np.allclose(
+    assert_allclose(
         v_s, lon, rtol=1e-5
     )  # Tolerances were rtol=1e-7, increased due to numpy v2 float32 changes (see #1603)
-    assert np.allclose(u_s, lat, rtol=1e-5)
+    assert_allclose(u_s, lat, rtol=1e-5)
 
 
 @pytest.mark.v4remove
@@ -168,7 +169,7 @@ def test_pset_from_field():
     pdens = np.histogram2d(pset.lat, pset.lon, bins=[np.linspace(0.0, 1.0, ydim + 1), np.linspace(0.0, 1.0, xdim + 1)])[
         0
     ]
-    assert np.allclose(pdens / sum(pdens.flatten()), startfield / sum(startfield.flatten()), atol=1e-2)
+    assert_allclose(pdens / sum(pdens.flatten()), startfield / sum(startfield.flatten()), atol=1e-2)
 
 
 def test_nearest_neighbor_interpolation2D():
@@ -189,8 +190,8 @@ def test_nearest_neighbor_interpolation2D():
     xv, yv = np.meshgrid(np.linspace(0.0, 1.0, int(np.sqrt(npart))), np.linspace(0.0, 1.0, int(np.sqrt(npart))))
     pset = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten())
     pset.execute(SampleP, endtime=1, dt=1)
-    assert np.allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5)], 1.0, rtol=1e-5)
-    assert np.allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5)], 0.0, rtol=1e-5)
+    assert_allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5)], 1.0, rtol=1e-5)
+    assert_allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5)], 0.0, rtol=1e-5)
 
 
 def test_nearest_neighbor_interpolation3D():
@@ -215,8 +216,8 @@ def test_nearest_neighbor_interpolation3D():
     pset2 = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten(), depth=np.ones(npart))
     pset.add(pset2)
     pset.execute(SampleP, endtime=1, dt=1)
-    assert np.allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5) & (pset.depth > 0.5)], 1.0, rtol=1e-5)
-    assert np.allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5) & (pset.depth < 0.5)], 0.0, rtol=1e-5)
+    assert_allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5) & (pset.depth > 0.5)], 1.0, rtol=1e-5)
+    assert_allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5) & (pset.depth < 0.5)], 0.0, rtol=1e-5)
 
 
 @pytest.mark.parametrize("withDepth", [True, False])
@@ -249,7 +250,7 @@ def test_inversedistance_nearland(withDepth, arrtype):
     if arrtype == "rand":
         assert np.all((pset.p > 2) & (pset.p < 3))
     else:
-        assert np.allclose(pset.p, 1.0, rtol=1e-5)
+        assert_allclose(pset.p, 1.0, rtol=1e-5)
 
     success = False
     try:
@@ -298,17 +299,17 @@ def test_partialslip_nearland_zonal(boundaryslip, withW, withT):
     kernel = AdvectionRK4_3D if withW else AdvectionRK4
     pset.execute(kernel, endtime=2, dt=1)
     if boundaryslip == "partialslip":
-        assert np.allclose([p.lon for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1)
-        assert np.allclose([pset[0].lon, pset[-1].lon], 0.06)
-        assert np.allclose([pset[1].lon, pset[-2].lon], 0.08)
+        assert_allclose([p.lon for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1)
+        assert_allclose([pset[0].lon, pset[-1].lon], 0.06)
+        assert_allclose([pset[1].lon, pset[-2].lon], 0.08)
         if withW:
-            assert np.allclose([p.depth for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1)
-            assert np.allclose([pset[0].depth, pset[-1].depth], 0.06)
-            assert np.allclose([pset[1].depth, pset[-2].depth], 0.08)
+            assert_allclose([p.depth for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1)
+            assert_allclose([pset[0].depth, pset[-1].depth], 0.06)
+            assert_allclose([pset[1].depth, pset[-2].depth], 0.08)
     else:
-        assert np.allclose([p.lon for p in pset], 0.1)
+        assert_allclose([p.lon for p in pset], 0.1)
         if withW:
-            assert np.allclose([p.depth for p in pset], 0.1)
+            assert_allclose([p.depth for p in pset], 0.1)
 
 
 @pytest.mark.parametrize("boundaryslip", ["freeslip", "partialslip"])
@@ -340,17 +341,17 @@ def test_partialslip_nearland_meridional(boundaryslip, withW):
     kernel = AdvectionRK4_3D if withW else AdvectionRK4
     pset.execute(kernel, endtime=2, dt=1)
     if boundaryslip == "partialslip":
-        assert np.allclose([p.lat for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1)
-        assert np.allclose([pset[0].lat, pset[-1].lat], 0.06)
-        assert np.allclose([pset[1].lat, pset[-2].lat], 0.08)
+        assert_allclose([p.lat for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1)
+        assert_allclose([pset[0].lat, pset[-1].lat], 0.06)
+        assert_allclose([pset[1].lat, pset[-2].lat], 0.08)
         if withW:
-            assert np.allclose([p.depth for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1)
-            assert np.allclose([pset[0].depth, pset[-1].depth], 0.06)
-            assert np.allclose([pset[1].depth, pset[-2].depth], 0.08)
+            assert_allclose([p.depth for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1)
+            assert_allclose([pset[0].depth, pset[-1].depth], 0.06)
+            assert_allclose([pset[1].depth, pset[-2].depth], 0.08)
     else:
-        assert np.allclose([p.lat for p in pset], 0.1)
+        assert_allclose([p.lat for p in pset], 0.1)
         if withW:
-            assert np.allclose([p.depth for p in pset], 0.1)
+            assert_allclose([p.depth for p in pset], 0.1)
 
 
 @pytest.mark.parametrize("boundaryslip", ["freeslip", "partialslip"])
@@ -372,13 +373,13 @@ def test_partialslip_nearland_vertical(boundaryslip):
     )
     pset.execute(AdvectionRK4, endtime=2, dt=1)
     if boundaryslip == "partialslip":
-        assert np.allclose([p.lon for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1)
-        assert np.allclose([p.lat for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1)
-        assert np.allclose([pset[0].lon, pset[-1].lon, pset[0].lat, pset[-1].lat], 0.06)
-        assert np.allclose([pset[1].lon, pset[-2].lon, pset[1].lat, pset[-2].lat], 0.08)
+        assert_allclose([p.lon for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1)
+        assert_allclose([p.lat for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1)
+        assert_allclose([pset[0].lon, pset[-1].lon, pset[0].lat, pset[-1].lat], 0.06)
+        assert_allclose([pset[1].lon, pset[-2].lon, pset[1].lat, pset[-2].lat], 0.08)
     else:
-        assert np.allclose([p.lon for p in pset], 0.1)
-        assert np.allclose([p.lat for p in pset], 0.1)
+        assert_allclose([p.lon for p in pset], 0.1)
+        assert_allclose([p.lat for p in pset], 0.1)
 
 
 def test_fieldset_sample_particle():
@@ -396,11 +397,11 @@ def test_fieldset_sample_particle():
     lat = np.linspace(-80, 80, npart)
     pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.v, lon, rtol=1e-6)
+    assert_allclose(pset.v, lon, rtol=1e-6)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.u, lat, rtol=1e-6)
+    assert_allclose(pset.u, lat, rtol=1e-6)
 
 
 def test_fieldset_sample_geographic(fieldset_geometric):
@@ -412,11 +413,11 @@ def test_fieldset_sample_geographic(fieldset_geometric):
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.v, lon, rtol=1e-6)
+    assert_allclose(pset.v, lon, rtol=1e-6)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.u, lat, rtol=1e-6)
+    assert_allclose(pset.u, lat, rtol=1e-6)
 
 
 def test_fieldset_sample_geographic_noconvert(fieldset_geometric):
@@ -428,11 +429,11 @@ def test_fieldset_sample_geographic_noconvert(fieldset_geometric):
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0)
     pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.v, lon * 1000 * 1.852 * 60, rtol=1e-6)
+    assert_allclose(pset.v, lon * 1000 * 1.852 * 60, rtol=1e-6)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0)
     pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.u, lat * 1000 * 1.852 * 60, rtol=1e-6)
+    assert_allclose(pset.u, lat * 1000 * 1.852 * 60, rtol=1e-6)
 
 
 def test_fieldset_sample_geographic_polar(fieldset_geometric_polar):
@@ -444,11 +445,11 @@ def test_fieldset_sample_geographic_polar(fieldset_geometric_polar):
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.v, lon, rtol=1e-6)
+    assert_allclose(pset.v, lon, rtol=1e-6)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0)
     pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.u, lat, rtol=1e-2)
+    assert_allclose(pset.u, lat, rtol=1e-2)
 
 
 def test_meridionalflow_spherical():
@@ -559,27 +560,27 @@ def test_sampling_out_of_bounds_time(allow_time_extrapolation):
     pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=-1.0)
     if allow_time_extrapolation:
         pset.execute(SampleP, endtime=-0.9, dt=0.1)
-        assert np.allclose(pset.p, 0.0, rtol=1e-5)
+        assert_allclose(pset.p, 0.0, rtol=1e-5)
     else:
         with pytest.raises(RuntimeError):
             pset.execute(SampleP, endtime=-0.9, dt=0.1)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=0)
     pset.execute(SampleP, runtime=0.1, dt=0.1)
-    assert np.allclose(pset.p, 0.0, rtol=1e-5)
+    assert_allclose(pset.p, 0.0, rtol=1e-5)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=0.5)
     pset.execute(SampleP, runtime=0.1, dt=0.1)
-    assert np.allclose(pset.p, 0.5, rtol=1e-5)
+    assert_allclose(pset.p, 0.5, rtol=1e-5)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=1.0)
     pset.execute(SampleP, runtime=0.1, dt=0.1)
-    assert np.allclose(pset.p, 1.0, rtol=1e-5)
+    assert_allclose(pset.p, 1.0, rtol=1e-5)
 
     pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=2.0)
     if allow_time_extrapolation:
         pset.execute(SampleP, runtime=0.1, dt=0.1)
-        assert np.allclose(pset.p, 1.0, rtol=1e-5)
+        assert_allclose(pset.p, 1.0, rtol=1e-5)
     else:
         with pytest.raises(RuntimeError):
             pset.execute(SampleP, runtime=0.1, dt=0.1)
@@ -655,7 +656,7 @@ def test_sample(particle, fieldset, time):  # pragma: no cover
 
     kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample)
     pset.execute(kernels, runtime=10, dt=1)
-    assert np.allclose(pset.sample_var, 10.0)
+    assert_allclose(pset.sample_var, 10.0)
 
 
 @pytest.mark.parametrize("npart", [1, 10])
@@ -694,7 +695,7 @@ def test_sample(particle, fieldset, time):  # pragma: no cover
 
     kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample)
     pset.execute(kernels, runtime=10, dt=1)
-    assert np.allclose(pset.sample_var, 10.0)
+    assert_allclose(pset.sample_var, 10.0)
 
 
 @pytest.mark.parametrize("ugridfactor", [1, 10])
@@ -830,7 +831,7 @@ def Recover(particle, fieldset, time):  # pragma: no cover
     pset.execute(pset.Kernel(AdvectionRK4) + SampleP + Recover, runtime=1, dt=1)
     assert np.isclose(pset.lat[0], 0)
     assert np.isclose(pset.p[0], 999)
-    assert np.allclose(fieldset.UV[0][0, 0, 0, 0], [0.1, 0.2])
+    assert_allclose(fieldset.UV[0][0, 0, 0, 0], [0.1, 0.2])
 
 
 def test_fieldset_sampling_updating_order(tmp_zarrfile):
diff --git a/tests/test_grids.py b/tests/test_grids.py
index 21a3d2ea7..ca1f15096 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels import (
     AdvectionRK4,
@@ -243,7 +244,7 @@ def sampleTemp(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet.from_list(fieldset, MyParticle, lon=[lon], lat=[lat], depth=[bath_func(lon) * ratio])
 
     pset.execute(pset.Kernel(sampleTemp), runtime=1)
-    assert np.allclose(pset.temp[0], ratio, atol=1e-4)
+    assert_allclose(pset.temp[0], ratio, atol=1e-4)
 
 
 def test_rectilinear_s_grids_advect1():
@@ -286,7 +287,7 @@ def bath_func(lon):
     pset = ParticleSet.from_list(fieldset, Particle, lon=lon, lat=lat, depth=depth)
 
     pset.execute(AdvectionRK4_3D, runtime=10000, dt=500)
-    assert np.allclose(pset.depth / bath_func(pset.lon), ratio)
+    assert_allclose(pset.depth / bath_func(pset.lon), ratio)
 
 
 def test_rectilinear_s_grids_advect2():
@@ -330,7 +331,7 @@ def moveEast(particle, fieldset, time):  # pragma: no cover
     kernel = pset.Kernel(moveEast)
     for _ in range(10):
         pset.execute(kernel, runtime=100, dt=50)
-        assert np.allclose(pset.relDepth[0], depth / bath_func(pset.lon[0]))
+        assert_allclose(pset.relDepth[0], depth / bath_func(pset.lon[0]))
 
 
 def test_curvilinear_grids():
@@ -362,7 +363,7 @@ def sampleSpeed(particle, fieldset, time):  # pragma: no cover
 
     pset = ParticleSet.from_list(fieldset, MyParticle, lon=[400, -200], lat=[600, 600])
     pset.execute(pset.Kernel(sampleSpeed), runtime=1)
-    assert np.allclose(pset.speed[0], 1000)
+    assert_allclose(pset.speed[0], 1000)
 
 
 def test_nemo_grid():
@@ -627,10 +628,10 @@ def OutBoundsError(particle, fieldset, time):  # pragma: no cover
         assert pset.out_of_bounds[1] == 0
         assert pset.out_of_bounds[2] == 1
     else:
-        assert np.allclose(pset.zonal, 0.015)
-        assert np.allclose(pset.meridional, 0.01)
-        assert np.allclose(pset.vert, -0.01)
-        assert np.allclose(pset.tracer, 1)
+        assert_allclose(pset.zonal, 0.015)
+        assert_allclose(pset.meridional, 0.01)
+        assert_allclose(pset.vert, -0.01)
+        assert_allclose(pset.tracer, 1)
 
 
 @pytest.mark.parametrize("gridindexingtype", ["mitgcm", "nemo"])
@@ -716,7 +717,7 @@ def UpdateR(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet(fieldset, pclass=MyParticle, lon=0, lat=4e3, time=0)
 
     pset.execute(pset.Kernel(UpdateR) + AdvectionRK4, runtime=timedelta(hours=14), dt=timedelta(minutes=5))
-    assert np.allclose(pset.radius, pset.radius_start, atol=10)
+    assert_allclose(pset.radius, pset.radius_start, atol=10)
 
 
 @pytest.mark.parametrize("gridindexingtype", ["mitgcm", "nemo"])
@@ -792,7 +793,7 @@ def UpdateR(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet(fieldset, pclass=MyParticle, depth=4e3, lon=0, lat=0, time=0)
 
     pset.execute(pset.Kernel(UpdateR) + AdvectionRK4_3D, runtime=timedelta(hours=14), dt=timedelta(minutes=5))
-    assert np.allclose(pset.radius, pset.radius_start, atol=10)
+    assert_allclose(pset.radius, pset.radius_start, atol=10)
 
 
 @pytest.mark.parametrize("gridindexingtype", ["pop", "mom5"])
@@ -869,7 +870,7 @@ def UpdateR(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet(fieldset, pclass=MyParticle, depth=-9.995e3, lon=0, lat=0, time=0)
 
     pset.execute(pset.Kernel(UpdateR) + AdvectionRK4_3D, runtime=timedelta(hours=14), dt=timedelta(minutes=5))
-    assert np.allclose(pset.radius, pset.radius_start, atol=10)
+    assert_allclose(pset.radius, pset.radius_start, atol=10)
 
 
 @pytest.mark.parametrize("gridindexingtype", ["mom5"])  # TODO v4: add pop in params?
@@ -971,10 +972,10 @@ def VelocityInterpolator(particle, fieldset, time):  # pragma: no cover
 
         convfactor = 0.01 if gridindexingtype == "pop" else 1.0
         if pointtype in ["U", "V"]:
-            assert np.allclose(pset.Uvel[0], u * convfactor)
-            assert np.allclose(pset.Vvel[0], v * convfactor)
+            assert_allclose(pset.Uvel[0], u * convfactor)
+            assert_allclose(pset.Vvel[0], v * convfactor)
         elif pointtype == "W":
             if extrapolation:
-                assert np.allclose(pset.Wvel[0], 0, atol=1e-9)
+                assert_allclose(pset.Wvel[0], 0, atol=1e-9)
             else:
-                assert np.allclose(pset.Wvel[0], w * convfactor)
+                assert_allclose(pset.Wvel[0], w * convfactor)
diff --git a/tests/test_interaction.py b/tests/test_interaction.py
index 7c0e4fc08..32450cbed 100644
--- a/tests/test_interaction.py
+++ b/tests/test_interaction.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
 from parcels import Field, FieldSet, ParticleSet
 from parcels.application_kernels.advection import AdvectionRK4
@@ -56,7 +57,7 @@ def test_simple_interaction_kernel(fieldset_unit_mesh):
         interaction_distance=interaction_distance,
     )
     pset.execute(DoNothing, pyfunc_inter=DummyMoveNeighbor, endtime=2.0, dt=1.0)
-    assert np.allclose(pset.lat, [0.1, 0.2, 0.1, 0.0], rtol=1e-5)
+    assert_allclose(pset.lat, [0.1, 0.2, 0.1, 0.0], rtol=1e-5)
 
 
 @pytest.mark.parametrize("mesh", ["spherical", "flat"])
@@ -75,10 +76,10 @@ def test_zonal_periodic_distance(mesh, periodic_domain_zonal):
     )
     pset.execute(DoNothing, pyfunc_inter=DummyMoveNeighbor, endtime=2.0, dt=1.0)
     if periodic_domain_zonal:
-        assert np.allclose([pset[0].lat, pset[2].lat], 0.6)
-        assert np.allclose(pset[1].lat, 0.5)
+        assert_allclose([pset[0].lat, pset[2].lat], 0.6)
+        assert_allclose(pset[1].lat, 0.5)
     else:
-        assert np.allclose([p.lat for p in pset], 0.5)
+        assert_allclose([p.lat for p in pset], 0.5)
 
 
 def test_concatenate_interaction_kernels(fieldset_unit_mesh):
@@ -103,7 +104,7 @@ def test_concatenate_interaction_kernels(fieldset_unit_mesh):
     # The kernel results are only applied after all interactionkernels
     # have been executed, so we expect the result to be double the
     # movement from executing the kernel once.
-    assert np.allclose(pset.lat, [0.2, 0.4, 0.2, 0.0], rtol=1e-5)
+    assert_allclose(pset.lat, [0.2, 0.4, 0.2, 0.0], rtol=1e-5)
 
 
 def test_concatenate_interaction_kernels_as_pyfunc(fieldset_unit_mesh):
@@ -125,7 +126,7 @@ def test_concatenate_interaction_kernels_as_pyfunc(fieldset_unit_mesh):
     # The kernel results are only applied after all interactionkernels
     # have been executed, so we expect the result to be double the
     # movement from executing the kernel once.
-    assert np.allclose(pset.lat, [0.2, 0.4, 0.2, 0.0], rtol=1e-5)
+    assert_allclose(pset.lat, [0.2, 0.4, 0.2, 0.0], rtol=1e-5)
 
 
 def test_neighbor_merge(fieldset_unit_mesh):
@@ -216,7 +217,7 @@ def compare_results_by_idx(instance, particle_idx, ref_result, active_idx=None):
         assert neigh in active_idx
     assert set(cur_neigh) <= set(active_idx)
     neigh_by_coor, _ = instance.find_neighbors_by_coor(instance._values[:, particle_idx])
-    assert np.allclose(cur_neigh, neigh_by_coor)
+    assert_allclose(cur_neigh, neigh_by_coor)
 
     assert isinstance(cur_neigh, np.ndarray)
     assert set(ref_result) == set(cur_neigh)
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index 8de32ed96..b73d552b2 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
 import parcels
 from parcels import (
@@ -64,10 +65,10 @@ def SampleP(particle, fieldset, time):  # pragma: no cover
     if kernel_type == "update_dlon":
         assert np.isclose(lons[0], lons[1])
         assert np.isclose(ps[0], ps[1])
-        assert np.allclose(lons[0], 0)
+        assert_allclose(lons[0], 0)
     else:
         assert np.isclose(ps[0] - ps[1], 0.1)
-        assert np.allclose(lons[0], 0.2)
+        assert_allclose(lons[0], 0.2)
 
 
 @pytest.mark.parametrize(
@@ -87,7 +88,7 @@ def test_execution_endtime(fieldset_unit_mesh, start, end, substeps, dt):
         fieldset_unit_mesh, pclass=Particle, time=start, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart)
     )
     pset.execute(DoNothing, endtime=end, dt=dt)
-    assert np.allclose(pset.time_nextloop, end)
+    assert_allclose(pset.time_nextloop, end)
 
 
 @pytest.mark.parametrize(
@@ -109,7 +110,7 @@ def test_execution_runtime(fieldset_unit_mesh, start, end, substeps, dt):
     t_step = abs(end - start) / substeps
     for _ in range(substeps):
         pset.execute(DoNothing, runtime=t_step, dt=dt)
-    assert np.allclose(pset.time_nextloop, end)
+    assert_allclose(pset.time_nextloop, end)
 
 
 def test_execution_fail_python_exception(fieldset_unit_mesh):
@@ -126,7 +127,7 @@ def PythonFail(particle, fieldset, time):  # pragma: no cover
         pset.execute(PythonFail, endtime=20.0, dt=2.0)
     assert len(pset) == npart
     assert np.isclose(pset.time[0], 10)
-    assert np.allclose(pset.time[1:], 0.0)
+    assert_allclose(pset.time[1:], 0.0)
 
 
 def test_execution_fail_out_of_bounds(fieldset_unit_mesh):
@@ -160,8 +161,8 @@ def MoveLeft(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=lon, lat=lat)
     pset.execute([MoveRight, MoveLeft], endtime=11.0, dt=1.0)
     assert len(pset) == npart
-    assert np.allclose(pset.lon, lon, rtol=1e-5)
-    assert np.allclose(pset.lat, lat, rtol=1e-5)
+    assert_allclose(pset.lon, lon, rtol=1e-5)
+    assert_allclose(pset.lat, lat, rtol=1e-5)
 
 
 def test_execution_check_all_errors(fieldset_unit_mesh):
@@ -209,8 +210,8 @@ def MoveRight(particle, fieldset, time):  # pragma: no cover
 def test_kernel_add_no_new_variables(fieldset_unit_mesh):
     pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5])
     pset.execute(pset.Kernel(MoveEast) + pset.Kernel(MoveNorth), endtime=2.0, dt=1.0)
-    assert np.allclose(pset.lon, 0.6, rtol=1e-5)
-    assert np.allclose(pset.lat, 0.6, rtol=1e-5)
+    assert_allclose(pset.lon, 0.6, rtol=1e-5)
+    assert_allclose(pset.lat, 0.6, rtol=1e-5)
 
 
 def test_multi_kernel_duplicate_varnames(fieldset_unit_mesh):
@@ -226,7 +227,7 @@ def Kernel2(particle, fieldset, time):  # pragma: no cover
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5])
     pset.execute([Kernel1, Kernel2], endtime=2.0, dt=1.0)
-    assert np.allclose(pset.lon, 0.3, rtol=1e-5)
+    assert_allclose(pset.lon, 0.3, rtol=1e-5)
 
 
 def test_multi_kernel_reuse_varnames(fieldset_unit_mesh):
@@ -241,7 +242,7 @@ def MoveEast2(particle, fieldset, time):  # pragma: no cover
 
     pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5])
     pset.execute(pset.Kernel(MoveEast1) + pset.Kernel(MoveEast2), endtime=2.0, dt=1.0)
-    assert np.allclose(pset.lon, [0.9], rtol=1e-5)  # should be 0.5 + 0.2 + 0.2 = 0.9
+    assert_allclose(pset.lon, [0.9], rtol=1e-5)  # should be 0.5 + 0.2 + 0.2 = 0.9
 
 
 def test_combined_kernel_from_list(fieldset_unit_mesh):
@@ -303,4 +304,4 @@ def MoveWest(particle, fieldset, time):  # pragma: no cover
     pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5])
     pset.execute(pset.Kernel(MoveEast), endtime=1.0, dt=1.0)
     pset.execute(pset.Kernel(MoveWest), endtime=3.0, dt=1.0)
-    assert np.allclose(pset.lon, 0.3, rtol=1e-5)  # should be 0.5 + 0.1 - 0.3 = 0.3
+    assert_allclose(pset.lon, 0.3, rtol=1e-5)  # should be 0.5 + 0.1 - 0.3 = 0.3
diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py
index f6722ddc2..b45d692c2 100644
--- a/tests/test_kernel_language.py
+++ b/tests/test_kernel_language.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
 from parcels import (
     Field,
@@ -127,7 +128,7 @@ def kernel(particle, fieldset, time):  # pragma: no cover
             particle.p *= 2.0
 
     pset.execute(kernel, endtime=1.0, dt=1.0)
-    assert np.allclose(pset.p, 20.0, rtol=1e-12)
+    assert_allclose(pset.p, 20.0, rtol=1e-12)
 
 
 def test_nested_if():
@@ -144,7 +145,7 @@ def kernel(particle, fieldset, time):  # pragma: no cover
                 particle.p1 = -1
 
     pset.execute(kernel, endtime=10, dt=1.0)
-    assert np.allclose([pset.p0[0], pset.p1[0]], [0, 1])
+    assert_allclose([pset.p0[0], pset.p1[0]], [0, 1])
 
 
 def test_pass():
@@ -157,7 +158,7 @@ def kernel(particle, fieldset, time):  # pragma: no cover
         pass
 
     pset.execute(kernel, endtime=10, dt=1.0)
-    assert np.allclose(pset[0].p, -1)
+    assert_allclose(pset[0].p, -1)
 
 
 def test_dt_as_variable_in_kernel():
@@ -220,7 +221,7 @@ def kernel(particle, fieldset, time):  # pragma: no cover
             particle.p += 1
 
     pset.execute(kernel, endtime=1.0, dt=1.0)
-    assert np.allclose(pset.p, 7.0, rtol=1e-12)
+    assert_allclose(pset.p, 7.0, rtol=1e-12)
     return
 
 
@@ -306,7 +307,7 @@ def test_small_dt(dt, expectation):
 
     with expectation:
         pset.execute(DoNothing, dt=dt, runtime=dt * 101)
-        assert np.allclose([p.time for p in pset], dt * 100)
+        assert_allclose([p.time for p in pset], dt * 100)
 
 
 def test_TEOSdensity_kernels():
@@ -336,7 +337,7 @@ def generate_fieldset(xdim=2, ydim=2, zdim=2, tdim=1):
     pset = ParticleSet(fieldset, pclass=DensParticle, lon=5, lat=5, depth=1000)
 
     pset.execute(PolyTEOS10_bsq, runtime=1)
-    assert np.allclose(pset[0].density, 1022.85377)
+    assert_allclose(pset[0].density, 1022.85377)
 
 
 def test_EOSseawaterproperties_kernels():
@@ -351,18 +352,18 @@ def test_EOSseawaterproperties_kernels():
     )
     pset = ParticleSet(fieldset, pclass=PoTempParticle, lon=5, lat=5, depth=1000)
     pset.execute(PtempFromTemp, runtime=1)
-    assert np.allclose(pset[0].potemp, 36.89073)
+    assert_allclose(pset[0].potemp, 36.89073)
 
     TempParticle = Particle.add_variables(
         [Variable("temp", dtype=np.float32), Variable("pressure", dtype=np.float32, initial=10000)]
     )
     pset = ParticleSet(fieldset, pclass=TempParticle, lon=5, lat=5, depth=1000)
     pset.execute(TempFromPtemp, runtime=1)
-    assert np.allclose(pset[0].temp, 40)
+    assert_allclose(pset[0].temp, 40)
 
     pset = ParticleSet(fieldset, pclass=TempParticle, lon=5, lat=30, depth=7321.45)
     pset.execute(PressureFromLatDepth, runtime=1)
-    assert np.allclose(pset[0].pressure, 7500, atol=1e-2)
+    assert_allclose(pset[0].pressure, 7500, atol=1e-2)
 
 
 @pytest.mark.parametrize("pressure", [0, 10])
@@ -397,6 +398,6 @@ def generate_fieldset(p, xdim=2, ydim=2, zdim=2, tdim=1):
     pset.execute(UNESCODensity, runtime=1)
 
     if pressure == 0:
-        assert np.allclose(pset[0].density, 1005.9465)
+        assert_allclose(pset[0].density, 1005.9465)
     elif pressure == 10:
-        assert np.allclose(pset[0].density, 1006.4179)
+        assert_allclose(pset[0].density, 1006.4179)
diff --git a/tests/test_mpirun.py b/tests/test_mpirun.py
index cb242fb4a..e050ed76b 100644
--- a/tests/test_mpirun.py
+++ b/tests/test_mpirun.py
@@ -1,9 +1,9 @@
 import os
 from glob import glob
 
-import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 
 from parcels._compat import MPI
 from tests.utils import PROJECT_ROOT
@@ -35,7 +35,7 @@ def test_mpi_run(tmpdir, repeatdt, maxage, nump):
         for v in ds2.variables.keys():
             if v == "time":
                 continue  # skip because np.allclose does not work well on np.datetime64
-            assert np.allclose(ds1.variables[v][:], ds2.variables[v][:], equal_nan=True)
+            assert_allclose(ds1.variables[v][:], ds2.variables[v][:], equal_nan=True)
 
         for a in ds2.attrs:
             if a != "parcels_version":
diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py
index 994166e0b..c77c2b377 100755
--- a/tests/test_particlefile.py
+++ b/tests/test_particlefile.py
@@ -6,6 +6,7 @@
 import numpy as np
 import pytest
 import xarray as xr
+from numpy.testing import assert_allclose
 from zarr.storage import MemoryStore
 
 import parcels
@@ -102,7 +103,7 @@ def test_pfile_array_remove_all_particles(fieldset, chunks_obs, tmp_zarrfile):
     pfile.write(pset, 2)
 
     ds = xr.open_zarr(tmp_zarrfile).load()
-    assert np.allclose(ds["time"][:, 0], np.timedelta64(0, "s"), atol=np.timedelta64(1, "ms"))
+    assert_allclose(ds["time"][:, 0], np.timedelta64(0, "s"), atol=np.timedelta64(1, "ms"))
     if chunks_obs is not None:
         assert ds["time"][:].shape == chunks
     else:
@@ -173,7 +174,7 @@ def Update_v(particle, fieldset, time):  # pragma: no cover
     ofile = pset.ParticleFile(name=tmp_zarrfile, outputdt=0.1)
     pset.execute(pset.Kernel(Update_v), endtime=1, dt=0.1, output_file=ofile)
 
-    assert np.allclose(pset.v_once - time - pset.age * 10, 1, atol=1e-5)
+    assert_allclose(pset.v_once - time - pset.age * 10, 1, atol=1e-5)
     ds = xr.open_zarr(tmp_zarrfile)
     vfile = np.ma.filled(ds["v_once"][:], np.nan)
     assert vfile.shape == (npart,)
@@ -213,12 +214,12 @@ def IncrLon(particle, fieldset, time):  # pragma: no cover
     samplevar = ds["sample_var"][:]
     if type == "repeatdt":
         assert samplevar.shape == (runtime // repeatdt, min(maxvar + 1, runtime))
-        assert np.allclose(pset.sample_var, np.arange(maxvar, -1, -repeatdt))
+        assert_allclose(pset.sample_var, np.arange(maxvar, -1, -repeatdt))
     elif type == "timearr":
         assert samplevar.shape == (runtime, min(maxvar + 1, runtime))
     # test whether samplevar[:, k] = k
     for k in range(samplevar.shape[1]):
-        assert np.allclose([p for p in samplevar[:, k] if np.isfinite(p)], k + 1)
+        assert_allclose([p for p in samplevar[:, k] if np.isfinite(p)], k + 1)
     filesize = os.path.getsize(str(tmp_zarrfile))
     assert filesize < 1024 * 65  # test that chunking leads to filesize less than 65KB
     ds.close()
@@ -323,7 +324,7 @@ def Update_lon(particle, fieldset, time):  # pragma: no cover
     ofile = pset.ParticleFile(name=tmp_zarrfile, outputdt=0.05)
     pset.execute(pset.Kernel(Update_lon), endtime=0.12, dt=0.02, output_file=ofile)
 
-    assert np.allclose(pset.lon, 0.6)
+    assert_allclose(pset.lon, 0.6)
 
 
 def test_correct_misaligned_outputdt_dt(fieldset, tmp_zarrfile):
@@ -337,10 +338,8 @@ def Update_lon(particle, fieldset, time):  # pragma: no cover
     pset.execute(pset.Kernel(Update_lon), endtime=11, dt=2, output_file=ofile)
 
     ds = xr.open_zarr(tmp_zarrfile)
-    assert np.allclose(ds.lon.values, [0, 3, 6, 9])
-    assert np.allclose(
-        ds.time.values[0, :], [np.timedelta64(t, "s") for t in [0, 3, 6, 9]], atol=np.timedelta64(1, "ns")
-    )
+    assert_allclose(ds.lon.values, [0, 3, 6, 9])
+    assert_allclose(ds.time.values[0, :], [np.timedelta64(t, "s") for t in [0, 3, 6, 9]], atol=np.timedelta64(1, "ns"))
 
 
 def setup_pset_execute(*, fieldset: FieldSet, outputdt: timedelta, execute_kwargs, particle_class=Particle):
diff --git a/tests/test_particles.py b/tests/test_particles.py
index 172d2a429..9ff3fa8df 100644
--- a/tests/test_particles.py
+++ b/tests/test_particles.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
 from parcels import (
     AdvectionRK4,
@@ -40,9 +41,9 @@ def addOne(particle, fieldset, time):  # pragma: no cover
         particle.p_int += 1
 
     pset.execute(pset.Kernel(AdvectionRK4) + addOne, runtime=1.0, dt=1.0)
-    assert np.allclose([p.p_float for p in pset], 11.0, rtol=1e-12)
-    assert np.allclose([p.p_double for p in pset], 12.0, rtol=1e-12)
-    assert np.allclose([p.p_int for p in pset], 13, rtol=1e-12)
+    assert_allclose([p.p_float for p in pset], 11.0, rtol=1e-12)
+    assert_allclose([p.p_double for p in pset], 12.0, rtol=1e-12)
+    assert_allclose([p.p_int for p in pset], 13, rtol=1e-12)
 
 
 @pytest.mark.parametrize("type", ["np.int8", "mp.float", "np.int16"])
@@ -82,7 +83,7 @@ def test_variable_init_relative(fieldset, coord_type):
     # Adjust base variable to test for aliasing effects
     for p in pset:
         p.p_base += 3.0
-    assert np.allclose([p.p_base for p in pset], 13.0, rtol=1e-12)
-    assert np.allclose([p.p_relative for p in pset], 10.0, rtol=1e-12)
-    assert np.allclose([p.p_lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.p_lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.p_base for p in pset], 13.0, rtol=1e-12)
+    assert_allclose([p.p_relative for p in pset], 10.0, rtol=1e-12)
+    assert_allclose([p.p_lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.p_lat for p in pset], lat, rtol=1e-12)
diff --git a/tests/test_particlesets.py b/tests/test_particlesets.py
index e7c6f58df..a2d46c632 100644
--- a/tests/test_particlesets.py
+++ b/tests/test_particlesets.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
 from parcels import (
     CurvilinearZGrid,
@@ -32,8 +33,8 @@ def test_pset_create_lon_lat(fieldset):
     lon = np.linspace(0, 1, npart, dtype=np.float32)
     lat = np.linspace(1, 0, npart, dtype=np.float32)
     pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle)
-    assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.lat for p in pset], lat, rtol=1e-12)
 
 
 @pytest.mark.parametrize("lonlatdepth_dtype", [np.float64, np.float32])
@@ -44,8 +45,8 @@ def test_pset_create_line(fieldset, lonlatdepth_dtype):
     pset = ParticleSet.from_line(
         fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=Particle, lonlatdepth_dtype=lonlatdepth_dtype
     )
-    assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.lat for p in pset], lat, rtol=1e-12)
     assert isinstance(pset[0].lat, lonlatdepth_dtype)
 
 
@@ -66,9 +67,9 @@ def test_pset_create_list_with_customvariable(fieldset):
 
     v_vals = np.arange(npart)
     pset = ParticleSet.from_list(fieldset, lon=lon, lat=lat, v=v_vals, pclass=MyParticle)
-    assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
-    assert np.allclose([p.v for p in pset], v_vals, rtol=1e-12)
+    assert_allclose([p.lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.v for p in pset], v_vals, rtol=1e-12)
 
 
 @pytest.mark.parametrize("restart", [True, False])
@@ -95,10 +96,10 @@ def Kernel(particle, fieldset, time):  # pragma: no cover
     )
 
     for var in ["lon", "lat", "depth", "time", "p", "p2", "p3"]:
-        assert np.allclose([getattr(p, var) for p in pset], [getattr(p, var) for p in pset_new])
+        assert_allclose([getattr(p, var) for p in pset], [getattr(p, var) for p in pset_new])
 
     if restart:
-        assert np.allclose([p.id for p in pset], [p.id for p in pset_new])
+        assert_allclose([p.id for p in pset], [p.id for p in pset_new])
     pset_new.execute(Kernel, runtime=2, dt=1)
     assert len(pset_new) == 3 * len(pset)
     assert pset[0].p3.dtype == np.float64
@@ -159,11 +160,11 @@ def test_pset_create_with_time(fieldset):
     lat = np.linspace(1, 0, npart)
     time = 5.0
     pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle, time=time)
-    assert np.allclose([p.time for p in pset], time, rtol=1e-12)
+    assert_allclose([p.time for p in pset], time, rtol=1e-12)
     pset = ParticleSet.from_list(fieldset, lon=lon, lat=lat, pclass=Particle, time=[time] * npart)
-    assert np.allclose([p.time for p in pset], time, rtol=1e-12)
+    assert_allclose([p.time for p in pset], time, rtol=1e-12)
     pset = ParticleSet.from_line(fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=Particle, time=time)
-    assert np.allclose([p.time for p in pset], time, rtol=1e-12)
+    assert_allclose([p.time for p in pset], time, rtol=1e-12)
 
 
 def test_pset_create_outside_time():
@@ -181,20 +182,20 @@ def Addlon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += particle.dt  # noqa
 
     pset.execute(Addlon, dt=1, runtime=2)
-    assert np.allclose([p.lon_nextloop for p in pset], [2 - t for t in times])
+    assert_allclose([p.lon_nextloop for p in pset], [2 - t for t in times])
 
 
 def test_pset_repeated_release(fieldset):
     npart = 10
     time = np.arange(0, npart, 1)  # release 1 particle every second
     pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), pclass=Particle, time=time)
-    assert np.allclose([p.time for p in pset], time)
+    assert_allclose([p.time for p in pset], time)
 
     def IncrLon(particle, fieldset, time):  # pragma: no cover
         particle_dlon += 1.0  # noqa
 
     pset.execute(IncrLon, dt=1.0, runtime=npart + 1)
-    assert np.allclose([p.lon for p in pset], np.arange(npart, 0, -1))
+    assert_allclose([p.lon for p in pset], np.arange(npart, 0, -1))
 
 
 def test_pset_repeatdt_check_dt(fieldset):
@@ -204,7 +205,7 @@ def IncrLon(particle, fieldset, time):  # pragma: no cover
         particle.lon = 1.0
 
     pset.execute(IncrLon, dt=2, runtime=21)
-    assert np.allclose([p.lon for p in pset], 1)  # if p.dt is nan, it won't be executed so p.lon will be 0
+    assert_allclose([p.lon for p in pset], 1)  # if p.dt is nan, it won't be executed so p.lon will be 0
 
 
 def test_pset_repeatdt_custominit(fieldset):
@@ -213,7 +214,7 @@ def test_pset_repeatdt_custominit(fieldset):
     pset = ParticleSet(fieldset, lon=0, lat=0, pclass=MyParticle, repeatdt=1, sample_var=5)
 
     pset.execute(DoNothing, dt=1, runtime=21)
-    assert np.allclose([p.sample_var for p in pset], 5.0)
+    assert_allclose([p.sample_var for p in pset], 5.0)
 
 
 def test_pset_stop_simulation(fieldset):
@@ -233,8 +234,8 @@ def test_pset_access(fieldset):
     lat = np.linspace(1, 0, npart, dtype=np.float32)
     pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle)
     assert pset.size == 100
-    assert np.allclose([pset[i].lon for i in range(pset.size)], lon, rtol=1e-12)
-    assert np.allclose([pset[i].lat for i in range(pset.size)], lat, rtol=1e-12)
+    assert_allclose([pset[i].lon for i in range(pset.size)], lon, rtol=1e-12)
+    assert_allclose([pset[i].lat for i in range(pset.size)], lat, rtol=1e-12)
 
 
 def test_pset_custom_ptype(fieldset):
@@ -243,8 +244,8 @@ def test_pset_custom_ptype(fieldset):
 
     pset = ParticleSet(fieldset, pclass=TestParticle, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart))
     assert pset.size == npart
-    assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), atol=1e-5)
-    assert np.allclose([p.n - 2 for p in pset], np.zeros(npart), rtol=1e-12)
+    assert_allclose([p.p - 0.33 for p in pset], np.zeros(npart), atol=1e-5)
+    assert_allclose([p.n - 2 for p in pset], np.zeros(npart), rtol=1e-12)
 
 
 def test_pset_add_explicit(fieldset):
@@ -256,8 +257,8 @@ def test_pset_add_explicit(fieldset):
         particle = ParticleSet(pclass=Particle, lon=lon[i], lat=lat[i], fieldset=fieldset, lonlatdepth_dtype=np.float64)
         pset.add(particle)
     assert pset.size == npart
-    assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.lat for p in pset], lat, rtol=1e-12)
 
 
 def test_pset_add_shorthand(fieldset):
@@ -268,8 +269,8 @@ def test_pset_add_shorthand(fieldset):
     for i in range(npart):
         pset += ParticleSet(pclass=Particle, lon=lon[i], lat=lat[i], fieldset=fieldset)
     assert pset.size == npart
-    assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
-    assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
+    assert_allclose([p.lon for p in pset], lon, rtol=1e-12)
+    assert_allclose([p.lat for p in pset], lat, rtol=1e-12)
 
 
 def test_pset_add_execute(fieldset):
@@ -283,7 +284,7 @@ def AddLat(particle, fieldset, time):  # pragma: no cover
         pset += ParticleSet(pclass=Particle, lon=0.1, lat=0.1, fieldset=fieldset)
     for _ in range(4):
         pset.execute(pset.Kernel(AddLat), runtime=1.0, dt=1.0)
-    assert np.allclose(np.array([p.lat for p in pset]), 0.4, rtol=1e-12)
+    assert_allclose(np.array([p.lat for p in pset]), 0.4, rtol=1e-12)
 
 
 def test_pset_merge_inplace(fieldset):
@@ -355,7 +356,7 @@ def AddLat(particle, fieldset, time):  # pragma: no cover
     k_add = pset.Kernel(AddLat)
     for _ in range(n + 1):
         pset.execute(k_add, runtime=1.0, dt=1.0)
-    assert np.allclose([p.lat - n * 0.1 for p in pset], np.zeros(npart), rtol=1e-12)
+    assert_allclose([p.lat - n * 0.1 for p in pset], np.zeros(npart), rtol=1e-12)
 
 
 def test_pset_multi_execute_delete(fieldset):
@@ -370,7 +371,7 @@ def AddLat(particle, fieldset, time):  # pragma: no cover
     for _ in range(n + 1):
         pset.execute(k_add, runtime=1.0, dt=1.0)
         pset.remove_indices(-1)
-    assert np.allclose(pset.lat, n * 0.1, atol=1e-12)
+    assert_allclose(pset.lat, n * 0.1, atol=1e-12)
 
 
 @pytest.mark.parametrize("staggered_grid", ["Agrid", "Cgrid"])
@@ -416,7 +417,7 @@ def SampleMask(particle, fieldset, time):  # pragma: no cover
 
     pset = ParticleSet.from_field(fieldset, size=400, pclass=SampleParticle, start_field=FMask, time=0)
     pset.execute(SampleMask, dt=1, runtime=1)
-    assert np.allclose([p.mask for p in pset], 1)
+    assert_allclose([p.mask for p in pset], 1)
     assert (np.array([p.lon for p in pset]) <= 1).all()
     test = np.logical_or(np.array([p.lon for p in pset]) <= 0, np.array([p.lat for p in pset]) >= 51)
     assert test.all()