Determining which ROMS grid particles are in when using nested fields

[jwongala]

Hello!

Is there a way to know what ROMS a particle is located in when using a nested field? I've used this code fieldset.U[1].grid.lon before when adding in displacement to define the lon and lat. And when using the below code I get this output but am unsure what it means. Are the ROMS references in the fieldset based on the order they are nested in? For example, UNested = NestedField('U', [Ufineres, Umhi] would mean that fieldset.U[0] is the fineres ROMS?

Thanks for any assistance!

[erikvansebille]

Hi @jwongala, what output do you get? Please provide more detailed information.

And does the last section of the Nested Fields tutorial, on keeping track on the Fields interpolated not provide some ideas?

[jwongala]

Hi @erikvansebille thank you for providing the link for keeping track of Fields interpolated, that does give me ideas!

[jwongala]

Hi @erikvansebille I used the keeping track on the Fields interpolated to create a grid of 1's for the coarse ROMS and 2's for the fine resolution ROMS I am using for my simulations. I nested the grids of 1's and 2's together and added them to the fieldset with no problems. But once I tried to execute the model I got the following error (screenshot) and I am unsure what it means. Any/all help is much appreciated! I have also attached my script as a .py for you if that is helpful for you. Thank you!

import numpy as np
import numpy.ma as ma
from netCDF4 import Dataset
import xarray as xr
import pandas as pd
import math
import scipy.integrate as spi
from scipy import interpolate

import matplotlib.pyplot as plt

from parcels import FieldSet, ParticleSet, JITParticle, ScipyParticle, DiffusionUniformKh, Variable, Field, GeographicPolar, Geographic, ErrorCode, plotTrajectoriesFile, NestedField, AdvectionRK4
from datetime import timedelta as delta

file_name="SINKTEST_toots_OSS_2023.zarr"
path2 = '/Users/wongalaj/TOOTS/TOOTS_parcels/' # give path to where .nc file is
feather_path = 'SINKTEST_toots_OSS_2023.feather' # assign name of feather file

startDate = '2023-01-01' # YYYY-MM-DD
endDate = '2023-01-03' # YYYY-MM-DD

run_days = 1   # day 2018: 240, 2019: 365, 2020: 366, 2021: 365
kh = 10   # This is the eddy diffusivity in m2/s
pld = 3   # in days

file1 = 'https://pae-paha.pacioos.hawaii.edu/erddap/griddap/roms_hiomsg' # OAHU SS

variables = {'U': 'u', 'V': 'v'}
dimensions = {'time': 'time', 'depth': 'depth','lat': 'latitude', 'lon': 'longitude'}

ds1 = xr.open_dataset(file1)  # this puts the opendap data into a xarray dataset

myDat1 = ds1.sel(**{'time': slice(startDate,endDate)})

dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time'}
dimensions = {'time': 'time', 'depth': 'depth','lat': 'latitude', 'lon': 'longitude'} #

U1=Field.from_xarray(myDat1.u, 'U', dimensions, allow_time_extrapolation=True)
V1=Field.from_xarray(myDat1.v, 'V', dimensions, allow_time_extrapolation=True)
S1=Field.from_xarray(myDat1.salt, 'sF', dimensions, allow_time_extrapolation=True)
T1=Field.from_xarray(myDat1.temp, 'tF', dimensions, allow_time_extrapolation=True)

file2 = 'https://pae-paha.pacioos.hawaii.edu/erddap/griddap/roms_hiig' # MHI ROMS (BASE)

ds2 = xr.open_dataset(file2)  # this puts the opendap data into a xarray dataset
myDat2 = ds2.sel(**{'time': slice(startDate,endDate)}) # subset based on time and depth layer

variables = {'U': 'u', 'V': 'v'}
dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time', 'depth': 'depth'}

U=Field.from_xarray(myDat2.u, 'U', dimensions, allow_time_extrapolation=True)
V=Field.from_xarray(myDat2.v, 'V', dimensions, allow_time_extrapolation=True)
S=Field.from_xarray(myDat2.salt, 'sC', dimensions, allow_time_extrapolation=True)
T=Field.from_xarray(myDat2.temp, 'tC', dimensions, allow_time_extrapolation=True)

# MHI (coarse model)
Umhi = U
Vmhi = V
Smhi = S
Tmhi = T

# OSS (fineres model)
Ufineres = U1
Vfineres = V1
Sfineres = S1
Tfineres = T1

# depth array for MHI ROMS
file_mhi = "/Users/wongalaj/TOOTS/TOOTS_parcels/depth_ncdf4MHI.nc"
ds_mhi = xr.open_dataset(file_mhi)  # this puts the opendap data into a xarray dataset

dimensions = {'lon': 'lon', 'lat': 'lat', 'depth': 'depth_layers'}
Dmhi = Field.from_xarray(ds_mhi.depth, 'depthC', dimensions)

# depth array for OSS ROMS
file_fineres = "/Users/wongalaj/TOOTS/TOOTS_parcels/depth_ncdf4OSS.nc"
ds_fineres = xr.open_dataset(file_fineres)  # this puts the opendap data into a xarray dataset

dimensions = {'lon': 'lon', 'lat': 'lat', 'depth': 'depth_layers'}
Dfineres = Field.from_xarray(ds_fineres.depth, 'depthF', dimensions)

# grid ID array for MHI
file_mhi = "/Users/wongalaj/TOOTS/TOOTS_parcels/nested_grid_identifier_ncdf4MHI.nc"
ds_mhi = xr.open_dataset(file_mhi)  # this puts the opendap data into a xarray dataset

dimensions = {'lon': 'lon', 'lat': 'lat'}
Fmhi = Field.from_xarray(ds_mhi.grid_def, 'gridC', dimensions)

# grid ID array for OSS
file_fineres = "/Users/wongalaj/TOOTS/TOOTS_parcels/nested_grid_identifier_ncdf4OSS.nc"
ds_fineres = xr.open_dataset(file_fineres)  # this puts the opendap data into a xarray dataset

dimensions = {'lon': 'lon', 'lat': 'lat'}
Ffineres = Field.from_xarray(ds_fineres.grid_def, 'gridF', dimensions)

UNested = NestedField('U', [Ufineres, Umhi])
VNested = NestedField('V', [Vfineres, Vmhi])

SNested = NestedField('S', [Sfineres, Smhi])
TNested = NestedField('T', [Tfineres, Tmhi])

DNested = NestedField('depth', [Dfineres, Dmhi])

FNested = NestedField('gridID', [Ffineres, Fmhi])

U = UNested
V = VNested
S = SNested
T = TNested
D = DNested
F = FNested

fieldset = FieldSet(U, V)

fieldset.add_field(D)
fieldset.add_field(S)
fieldset.add_field(T)
fieldset.add_field(F)

# kernel to update salt and temp at each time step

def update_ST_kernel(particle, fieldset, time):

    print("particle.time is", particle.time)

    print("particle.age is", particle.age)

    salt1 = fieldset.S[particle]
    salt_new = fieldset.S[time + particle.dt, particle.depth, particle.lat, particle.lon, particle]
    particle.salt = salt_new

    temp1 = fieldset.T[particle]
    temp_new = fieldset.T[time + particle.dt, particle.depth, particle.lat, particle.lon, particle]
    particle.temp = temp_new

    field1 = fieldset.gridID[particle]
    field_new = fieldset.gridID[particle.lat, particle.lon, particle]
    particle.fieldID = field_new

# kernel to sinkin particle based on what salt and what ROMS particle is in
def sinking_kernel(particle, fieldset, time):

    if particle.salt >=26 & particle.fieldID == 1.:
        u1 = fieldset.U[particle]
        v1 = fieldset.V[particle]
        grid_m = 4000 # grid res of MHI ROMS in meters

        w_oc = ((u1/grid_m)+(v1/grid_m)*particle.depth) # nmeters
        w_vel = 0.00011173 # value was calculated based on info from shapiro et al paper (all particle sunk in 20 mins after rolling for 24 hours in 450 ml jar. estimated height of the jar) 0.000004976 # original 0.00049769
        w_tot = w_oc + w_vel

        dep1 = particle.depth + (w_tot*particle.dt)

        particle.depth = dep1

    elif particle.salt >=26 & particle.fieldID==2.:
        u1 = fieldset.U[particle]
        v1 = fieldset.V[particle]
        grid_m = 200 # grid res of MHI ROMS in meters

        w_oc = ((u1/grid_m)+(v1/grid_m)*particle.depth) # nmeters
        w_vel = 0.00011173 # value was calculated based on info from shapiro et al paper (all particle sunk in 20 mins after rolling for 24 hours in 450 ml jar. estimated height of the jar) 0.000004976 # original 0.00049769
        w_tot = w_oc + w_vel

        dep1 = particle.depth + (w_tot*particle.dt)

        particle.depth = dep1

    else:

        particle.depth = particle.depth

# kernle to resuspend particle based on time and depth
def resuspension_kernel(particle, fieldset, time):

    if particle.time <= 86400 and particle.depth <= 25:
        w_vel = -0.00010 # random value in m/s
        dep1 = particle.depth + (w_vel*particle.dt)
        particle.depth = dep1
    elif particle.time >= 216000 and particle.time <= 259200:
        w_vel = -0.00010 # random value in m/s
        dep1 = particle.depth + (w_vel*particle.dt)
        particle.depth = dep1
    else:
        particle.depth = particle.depth


def make_landmask1(fielddata):
    """Returns landmask where land = 1 and ocean = 0
    fielddata is a netcdf file.
    """
    datafile = Dataset(fielddata)

    landmask = datafile.variables['u'][0, 0]
    landmask = np.ma.masked_invalid(landmask)
    landmask = landmask.mask.astype('int')

    return landmask


file_path1 = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_OAHU_SS/roms_hiomsg_2018_05_05_ffdb_1c31_10ba.nc'
file_path2 = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_MHI/roms_hiig_2018_05_05_685a_c38d_b3fe.nc'

landmask_fineres = make_landmask1(file_path1) # OAHU SS
landmask_mhi = make_landmask1(file_path2)

def get_coastal_nodes(landmask):
    """Function that detects the coastal nodes, i.e. the ocean nodes directly
    next to land. Computes the Laplacian of landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the coastal nodes, the coastal nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap -= 4*landmask
    coastal = np.ma.masked_array(landmask, mask_lap > 0)
    coastal = coastal.mask.astype('int')

    return coastal


def get_shore_nodes(landmask):
    """Function that detects the shore nodes, i.e. the land nodes directly
    next to the ocean. Computes the Laplacian of landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the shore nodes, the shore nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap -= 4*landmask
    shore = np.ma.masked_array(landmask, mask_lap < 0)
    shore = shore.mask.astype('int')

    return shore


def get_coastal_nodes_diagonal(landmask):
    """Function that detects the coastal nodes, i.e. the ocean nodes where
    one of the 8 nearest nodes is land. Computes the Laplacian of landmask
    and the Laplacian of the 45 degree rotated landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the coastal nodes, the coastal nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap += np.roll(landmask, (-1,1), axis=(0,1)) + np.roll(landmask, (1, 1), axis=(0,1))
    mask_lap += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (1, -1), axis=(0,1))
    mask_lap -= 8*landmask
    coastal = np.ma.masked_array(landmask, mask_lap > 0)
    coastal = coastal.mask.astype('int')

    return coastal


def get_shore_nodes_diagonal(landmask):
    """Function that detects the shore nodes, i.e. the land nodes where
    one of the 8 nearest nodes is ocean. Computes the Laplacian of landmask
    and the Laplacian of the 45 degree rotated landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the shore nodes, the shore nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap += np.roll(landmask, (-1,1), axis=(0,1)) + np.roll(landmask, (1, 1), axis=(0,1))
    mask_lap += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (1, -1), axis=(0,1))
    mask_lap -= 8*landmask
    shore = np.ma.masked_array(landmask, mask_lap < 0)
    shore = shore.mask.astype('int')

    return shore


coastal_fineres = get_coastal_nodes_diagonal(landmask_fineres)
shore_fineres = get_shore_nodes_diagonal(landmask_fineres)

coastal_mhi = get_coastal_nodes_diagonal(landmask_mhi)
shore_mhi = get_shore_nodes_diagonal(landmask_mhi)


def create_displacement_field(landmask, double_cell=False):
    """Function that creates a displacement field 1 m/s away from the shore.

    - landmask: the land mask dUilt using `make_landmask`.
    - double_cell: Boolean for determining if you want a double cell.
      Default set to False.

    Output: two 2D arrays, one for each camponent of the velocity.
    """
    shore = get_shore_nodes(landmask)
    shore_d = get_shore_nodes_diagonal(landmask) # bordering ocean directly and diagonally
    shore_c = shore_d - shore                    # corner nodes that only border ocean diagonally

    Ly = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0) # Simple derivative
    Lx = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)

    Ly_c = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0)
    Ly_c += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (-1,1), axis=(0,1)) # Include y-component of diagonal neighbours
    Ly_c += - np.roll(landmask, (1,-1), axis=(0,1)) - np.roll(landmask, (1,1), axis=(0,1))

    Lx_c = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)
    Lx_c += np.roll(landmask, (-1,-1), axis=(1,0)) + np.roll(landmask, (-1,1), axis=(1,0)) # Include x-component of diagonal neighbours
    Lx_c += - np.roll(landmask, (1,-1), axis=(1,0)) - np.roll(landmask, (1,1), axis=(1,0))

    v_x = -Lx*(shore)
    v_y = -Ly*(shore)

    v_x_c = -Lx_c*(shore_c)
    v_y_c = -Ly_c*(shore_c)

    v_x = v_x + v_x_c
    v_y = v_y + v_y_c

    magnitude = np.sqrt(v_y**2 + v_x**2)
    # the coastal nodes between land create a problem. Magnitude there is zero
    # I force it to be 1 to avoid problems when normalizing.
    ny, nx = np.where(magnitude == 0)
    magnitude[ny, nx] = 1

    v_x = v_x/magnitude
    v_y = v_y/magnitude

    return v_x, v_y


v_x_c, v_y_c = create_displacement_field(landmask_mhi)
v_x_f, v_y_f = create_displacement_field(landmask_fineres)


def distance_to_shore(landmask, dx=1):
    """Function that computes the distance to the shore. It is based in the
    the `get_coastal_nodes` algorithm.

    - landmask: the land mask dUilt using `make_landmask` function.
    - dx: the grid cell dimension. This is a crude approxsimation of the real
    distance (be careful).

    Output: 2D array containing the distances from shore.
    """
    ci = get_coastal_nodes(landmask) # direct neighbours
    dist = ci*dx                     # 1 dx away

    ci_d = get_coastal_nodes_diagonal(landmask) # diagonal neighbours
    dist_d = (ci_d - ci)*np.sqrt(2*dx**2)       # sqrt(2) dx away

    return dist+dist_d

# make for both coarse and fine
d_2_s_f = distance_to_shore(landmask_fineres)
d_2_s_c = distance_to_shore(landmask_mhi)



u_displacement_c = v_x_c
u_displacement_f = v_x_f
v_displacement_c = v_y_c
v_displacement_f = v_y_f

# fieldset for U fine
fieldset.add_field(Field('dispUF', data=u_displacement_f,
                         lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat,
                         mesh='spherical')) # have to index to choose which field we want to base it off of; 1 is choosing coarser

# fieldset for U coarse
fieldset.add_field(Field('dispUC', data=u_displacement_c,
                         lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat,
                         mesh='spherical'))

# fieldset for V fine
fieldset.add_field(Field('dispVF', data=v_displacement_f,
                         lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat,
                         mesh='spherical'))

# fieldset for V coarse
fieldset.add_field(Field('dispVC', data=v_displacement_c,
                          lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat,
                          mesh='spherical'))

fieldset.dispUF.units = GeographicPolar()
fieldset.dispUC.units = GeographicPolar()
fieldset.dispVF.units = Geographic()
fieldset.dispVC.units = Geographic()


fieldset.add_field(Field('landmask_fine', landmask_fineres,
                         lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat,
                         mesh='spherical'))

fieldset.add_field(Field('landmask_coarse', landmask_mhi,
                         lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat,
                         mesh='spherical'))

fieldset.add_field(Field('distance2shore_fine', d_2_s_f,
                         lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat,
                         mesh='spherical'))

fieldset.add_field(Field('distance2shore_coarse', d_2_s_c,
                         lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat,
                         mesh='spherical'))


# Add even diffusivity to the fieldset
fieldset.add_constant_field('Kh_zonal', kh, mesh='spherical')
fieldset.add_constant_field('Kh_meridional', kh, mesh='spherical')

# Add PLD to fieldset
fieldset.add_constant('pld', (pld*86400))


class DisplacementParticle(ScipyParticle):
    dU = Variable('dU', dtype=np.int32)
    dV = Variable('dV', dtype=np.int32)
    d2s = Variable('d2s', initial=1e3)
    temp = Variable('temp', dtype=np.float32)
    salt = Variable('salt', dtype=np.float32)
    age = Variable('age', dtype=np.float32, initial=0.)
    releaseSite = Variable('releaseSite', dtype=np.int32)
    IslandReleaseSite = Variable('IslandReleaseSite', dtype=np.int32)
    fieldID = Variable('fieldID', dtype=np.float32)


def set_displacement(particle, fieldset, time):
    # try saying in between max and min

    if  particle.lon <= fieldset.max_lon and particle.lon >= fieldset.min_lon and particle.lat <= fieldset.max_lat and particle.lat >= fieldset.min_lat:
        particle.d2s = fieldset.distance2shore_fine[time, particle.depth,
                               particle.lat, particle.lon]
        if  particle.d2s < 0.5:
            dispUab = fieldset.dispUF[time, particle.depth, particle.lat,
                               particle.lon]
            dispVab = fieldset.dispVF[time, particle.depth, particle.lat,
                               particle.lon]

            particle.dU = dispUab
            particle.dV = dispVab
        else:
            particle.dU = 0.
            particle.dV = 0.
    else:
        particle.d2s = fieldset.distance2shore_coarse[time, particle.depth,
                               particle.lat, particle.lon]

        if  particle.d2s < 0.5:
            dispUab = fieldset.dispUC[time, particle.depth, particle.lat,
                               particle.lon]
            dispVab = fieldset.dispVC[time, particle.depth, particle.lat,
                               particle.lon]


            particle.dU = dispUab
            particle.dV = dispVab
        else:
            particle.dU = 0.
            particle.dV = 0.


def displace(particle, fieldset, time):
    if  particle.d2s < 0.5:
        particle.lon += particle.dU*particle.dt
        particle.lat += particle.dV*particle.dt


def DeleteParticle(particle, fieldset, time):
    print('deleted particle')
    particle.delete()


def Ageing(particle, fieldset, time):
    particle.age += particle.dt
    if particle.age >= fieldset.pld:
        particle.delete()

# add min lon of base
fieldset.add_constant('min_lon', -163.83070)

# add max lon of base
fieldset.add_constant('max_lon', -152.51930)

# add min lat of base
fieldset.add_constant('min_lat', 17.01843)

# add max lat of base
fieldset.add_constant('max_lat', 23.98238)

# add min lon of fine
fieldset.add_constant('min_lon', -158.1218)

#add max lon of fine
fieldset.add_constant('max_lon', -157.7903)

#add min lat of fine
fieldset.add_constant('min_lat', 21.23194)

#add max lat of fine
fieldset.add_constant('max_lat', 21.32355)


habilon = [-157.8181, -157.9064, -158.1666] # lon_fin
habilat = [21.61089, 21.26621, 21.65717] # lat_fin
habisite = [42, 37, 34] # site_fin
islandsite = [3, 3, 3] # islandsite_fin
# density = [1023.6, 1023.6, 1023.6] # density_fin
# diameter = [0.00024, 0.00024, 0.00024] # diameter_fin
# par_rho = [20, 20, 20] # parrho_fin
depth_layer = [0.25, 0.25, 0.25] # depth_fin

# initialize salt, temp , and gridID for pset
salt_fin = np.array([])
temp_fin = np.array([])
f_fin = np.array([])

for i in range(0,3):
    salt_init = fieldset.S[0, depth_layer[i], habilat[i], habilon[i]]
    temp_init = fieldset.T[0, depth_layer[i], habilat[i], habilon[i]]
    f_init = fieldset.gridID[0, depth_layer[i],habilat[i], habilon[i]]


    salt_fin = np.append(salt_fin, salt_init)
    temp_fin = np.append(temp_fin, temp_init)
    f_fin = np.append(f_fin, f_init)


f_init = f_fin
temp_init = temp_fin # [25.1495826 , 25.20291293, 25.42359876] # temp_fin
salt_init = salt_fin # [34.48447797, 34.24041612, 34.5609318]# salt_fin


release_days = 4 # how often to release particles (original == 15 days)

release_int = release_days*86400 # 15 days converted to seconds

# Define the pset
pset = ParticleSet.from_list(fieldset=fieldset, #parameter found in ParticleSet
                             pclass=DisplacementParticle, #parameter found in ParticleSet
                             lon=habilon, #parameter found in ParticleSet
                             lat=habilat, #parameter found in ParticleSet
                             releaseSite=habisite,
                             IslandReleaseSite=islandsite,
                             depth = depth_layer,
                             fieldID = f_init,
                             temp = temp_init,
                             salt  = salt_init,
                             repeatdt=release_int) #parameter found in ParticleSet


kernels = pset.Kernel(displace) + pset.Kernel(AdvectionRK4) + pset.Kernel(sinking_kernel) + pset.Kernel(resuspension_kernel) + pset.Kernel(update_ST_kernel) + pset.Kernel(DiffusionUniformKh) + pset.Kernel(set_displacement) + pset.Kernel(Ageing) #


output_file = pset.ParticleFile(name=file_name, outputdt=delta(hours=0.25))


pset.execute(kernels,
             runtime=delta(days=run_days),
             dt=delta(minutes=15),
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
             output_file=output_file)

# now stop the repeated release
pset.repeatdt = None

# now continue running for the remaining length of the PLD
pset.execute(kernels,
            runtime=delta(days=pld+1),
            dt=delta(minutes=15),
            recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
            output_file=output_file)

[erikvansebille]

Hi @jwongala,

Looking at the error message (Error: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''), I think the issue is with your use of & in your custom kernels.

Try replacing statements such as

if particle.salt >=26 & particle.fieldID == 1.:

with

if particle.salt >=26 and particle.fieldID == 1.:

or possibly

if particle.salt >=26 && particle.fieldID == 1.:

[jwongala]

Hi @erikvansebille

This fixed my issue! Thank you very much!