Particles don't delete at domain boundary

[JonathanC33]

I am currently running ocean parcels using a gridded particle field with new releases every hour. When I execute the particle set, the particles advect towards the boundary of the domain, stay stationary next to it, and then are advected back into the center of the domain once the tidal direction changes. I've included the error.outofbounds delete particle definition but for some reason, it is not working in this script. The script is included below.

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation, writers, PillowWriter

import matplotlib
matplotlib.interactive(False)
from parcels import FieldSet, ParticleSet, Variable, JITParticle, AdvectionRK4, plotTrajectoriesFile, ErrorCode
import numpy as np
import math
from scipy import stats,interpolate,integrate
from datetime import timedelta
from operator import attrgetter
import xarray as xr

import cartopy.crs as ccrs
import cartopy
from shapely import geometry


#In this script I am creating a Ocean Parcels output to match the TRAP core time step
#An OP output will be created for each ToE listed below

ToI = {}
#Define the times of interest
ToI['Fall_Freshet'] = {}
ToI['Fall_Freshet'] = {
    2019: ['20191104','20191105'],
    2020: ['20201022','20201023']}
ToI['Winter_Low'] = {}
ToI['Winter_Low'] = {
    2019: ['20190301','20190302'],
    2020: ['20200222','20200223']}
ToI['Spring_High'] = {}
ToI['Spring_High'] = {
    2019: ['20190425','20190426'],
    2020: ['20200508','20200509']}
ToI['Summer_Low'] = {}
ToI['Summer_Low'] = {
    2019: ['20190801','20190802'],
    2020: ['20200715','20200716']}


#Import the Saint Johns small domain U and V velocity fields
filepath = '/home/joc000/Documents/Data/ocean_drifters_100m/'
filepath_const = '/home/stt001/gpfs4/CONSTANTS/SJ100/'
coord = filepath_const + 'Bathymetry_SJ100.nc'
mesh_mask = filepath_const + 'mesh-mask_SJ100.nc'
filepath_depth = filepath + 'NEMO_RPN_1h_grid_U_20191225-20191225.nc'

filepath_ocean_u = filepath+'NEMO_RPN_1h_grid_U_20*.nc'
filepath_ocean_v = filepath+'NEMO_RPN_1h_grid_V_20*.nc'


#Import the land mask
pfile = xr.open_dataset(mesh_mask)
land_mask = pfile.tmask.values[0,0].astype(float)
land_mask[land_mask == 0] = np.nan
lon = pfile.nav_lon.values
lat = pfile.nav_lat.values
pfile.close()

#Names for the U and V files, opens multiple
filenames = { 
'U': {'lon': coord, 'lat': coord, 'depth':filepath_depth, 'data': filepath_ocean_u},
'V': {'lon': coord, 'lat': coord, 'depth':filepath_depth, 'data': filepath_ocean_v}
}

#Define the dictionary of the variables and the dimensions, no depth in this case
variables = { 
'U': 'uo',
'V': 'vo'
}

dimensions = { 
'lat': 'nav_lat',
'lon': 'nav_lon',
'depth': 'depthu',
'time': 'time_instant'
}

#Take the depth at 0m (approximately, it's at 14.68m)
indices = {
'U': {'depth': [0]},
'V': {'depth': [0]}
}

#Fieldset using FieldSet.from_netcdf using the filenames, variables, and dimensions
fieldset_ocean = FieldSet.from_netcdf(
	filenames, 
	variables, 
	dimensions,
    indices=indices,
    deferred_load=True,
	allow_time_extrapolation=True) 


#Run through each of the ToEs, setting up particle grids and repeating every 3 hours
#Run for 2 days total to cover the entire time of interest
#Mask out land for the lons and lats used later for particle locations 
#Every 5th space, approx. 600m spacing 
lons = (lon*land_mask)[::5,::5]
lats = (lat*land_mask)[::5,::5]
#Release new particles every x hour(s) 
repeatdt = timedelta(hours=1)


for i in ToI.keys():
    for ii in ToI[i].keys():
        
        #Take the first day of each ToI
        iii = ToI[i][ii][0]
        
        
        #Define the ToE
        OP_start = np.datetime64(iii[:4]+'-'+iii[4:6]+'-'+iii[6:])
        #Define Particle Locations, and time according to ToE
        pset = ParticleSet.from_list(fieldset=fieldset_ocean, pclass=JITParticle,
        lon=lons.flatten()[~np.isnan(lons.flatten())],
        lat=lats.flatten()[~np.isnan(lats.flatten())],
        time=OP_start,
        repeatdt=repeatdt)
        
        #Don't save the depth of the particles (always 0)
        pset.set_variable_write_status('depth', False)
        
        def DeleteParticle(particle, fieldset, time):
            particle.delete()
            
        #Run the OP output for 2 days
        output_file = pset.ParticleFile(name='/gpfs/fs4/dfo/dpnm/dfo_dpnm/joc000/OceanParcelsOutput/ToI/'+\
                                        i+'/ToI-GriddedRelease_'+\
                                        i+'_'+iii+\
                                        '_100m-res_1hr-repeat.nc',
                                        outputdt=timedelta(minutes=60))
            
        kernels = pset.Kernel(AdvectionRK4)
        pset.execute(kernels,
        runtime=timedelta(hours=24*2),
        dt=timedelta(minutes=1),
        recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
        output_file=output_file)
        output_file.export()

[erikvansebille]

Dear @JonathanC33. Sorry, but it's very difficult to help you with a statement "[...] but for some reason, it is not working in this script." What si not working? How is it not working? You will need to provide more information if you want us to help

[JonathanC33]

Hi @erikvansebille, sorry for the lack of clarity in the question. Upon further discussing the issue with the team behind the modelling output, I found out about a set of non-physical zero values around the edge of the domain and this was the cause for the particles not leaving the domain.