OutOfBounds error despite periodic boundaries

[MiriamSterl]

I'm running particle tracking simulations on a flat square grid, with both x and y running from -L/2 to L/2 with L=7e5 m. I've added a halo to my fieldset via

fieldset.add_constant('halo_west', fieldset.U.grid.lon[0])
fieldset.add_constant('halo_east', fieldset.U.grid.lon[-1])
fieldset.add_constant('halo_south', fieldset.U.grid.lat[0])
fieldset.add_constant('halo_north', fieldset.U.grid.lat[-1])
fieldset.add_periodic_halo(zonal=True,meridional=True)

and added a kernel that implements periodic boundary conditions:

def periodicBC(particle, fieldset, time):
    if particle.lon < fieldset.halo_west:
        particle.lon += fieldset.halo_east - fieldset.halo_west
    elif particle.lon > fieldset.halo_east:
        particle.lon -= fieldset.halo_east - fieldset.halo_west
    if particle.lat < fieldset.halo_south:
        particle.lat += fieldset.halo_north - fieldset.halo_south
    elif particle.lat > fieldset.halo_north:
        particle.lat -= fieldset.halo_north - fieldset.halo_south

The simulation runs well at first, but at the end I get an OutOfBounds error:

OutOfBoundsError: 0
Particle P[26](lon=-197467.718750, lat=-349982.031250, depth=0.000000, u=-0.855570, v=-4.064669, dy=2070851.000000, time=103478400.000000)
Time: 103478400.0,	timestep dt: 3600.000000
Out-of-bounds sampling by particle at (-197467.718750, -349982.031250, 0.000000)

Why do I suddenly get this error despite the periodic boundaries? The only thing I can think of is that the particle that gets the error is very close to the southern boundary -L/2 = -350000, and the meridional velocity v is quite strongly negative, so the particle will move out of the southern boundary. But why is that not caught by the periodic BC Kernel?
NB: when I run the simulation with a smaller timestep (dt = 1800 s instead of 3600 s), I don't get the error.

(Small note: my flow fields are 2D, hence why the depth of the particles remains 0.)

[erikvansebille]

Thanks for the question, @MiriamSterl. This could be because you are using AdvectionRK4 (see code below), which requires interpolation of the field at intermediate time steps
https://github.com/OceanParcels/parcels/blob/54ede99e63f97e14faca12de23b02d713ac0762d/parcels/application_kernels/advection.py#L10-L22
One solution would be to create your own Advection Kernel, based on the kernel above, but where checks are done when e.g., lat1, lat2 or lat3 are outside of the domain. Makes sense?

[MiriamSterl]

Thanks Erik! That does make sense, I'll give that a try.