Particles stuck on staggered C grid

[ckleeaj]

Question

Question

Hi everyone,

I am computing lagrangian velocity with my CROCO eulerian output. CROCO uses a staggered c grid and I have tried creating fieldset according to it. However, the particles still get stuck in my simulation.

I have 12 output files, which include data of every different month, and I need to advect particles for one year.

Below is a snapshot of my simulation, which is at day 47 (the output frequency is 1 day) and the timestep used to advect the particles is 5 mins.

Supporting code/error messages

I have tried generating a fieldset with from_c_grid_dataset, from_nemo and from_mitgcm. However, the same error (index 0 is out of bounds for axis 0 with size 0) popped out.

# Import the mesh mask (grid file)
mesh_mask = r'E:\CROCO\MD_2020_hmin\croco_grd.nc'

# Import the CROCO output files
dataset_folder = "E:/CROCO/MD_2020_hmin/updated_files2/"
data_files = sorted(glob.glob(f"{dataset_folder}/*.nc"))
fname = f"{dataset_folder}/*.nc"

# Create Fieldset
filenames = {"U": {"lon": mesh_mask, "lat": mesh_mask, "data": fname}, 
             "V": {"lon": mesh_mask, "lat": mesh_mask, "data": fname}}

# Get velocity with land mask applied
variables = {"U": "masked_u", "V": "masked_v",}

# Get the lon and lat from f_point (psi)
dimensions = {"lat": "lat_psi", "lon": "lon_psi", "time": "time"}

fieldset_cgrid = parcels.FieldSet.from_c_grid_dataset(filenames, variables, dimensions)

{
	"name": "IndexError",
	"message": "index 0 is out of bounds for axis 0 with size 0",
	"stack": "---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[4], line 15
     11 variables = {\"U\": \"masked_u\", \"V\": \"masked_v\",}
     13 dimensions = {\"lat\": \"lat_psi\", \"lon\": \"lon_psi\", \"time\": \"time\"}
---> 15 fieldset_cgrid = parcels.FieldSet.from_c_grid_dataset(filenames, variables, dimensions)

File c:\\Users\\leeka\\anaconda3\\envs\\parcels\\Lib\\site-packages\\parcels\\fieldset.py:673, in FieldSet.from_c_grid_dataset(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, gridindexingtype, chunksize, **kwargs)
    670 if 'creation_log' not in kwargs.keys():
    671     kwargs['creation_log'] = 'from_c_grid_dataset'
--> 673 return cls.from_netcdf(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    674                        allow_time_extrapolation=allow_time_extrapolation, interp_method=interp_method,
    675                        chunksize=chunksize, gridindexingtype=gridindexingtype, **kwargs)

File c:\\Users\\leeka\\anaconda3\\envs\\parcels\\Lib\\site-packages\\parcels\\fieldset.py:454, in FieldSet.from_netcdf(cls, filenames, variables, dimensions, indices, fieldtype, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, chunksize, **kwargs)
    452                     dFiles = fields[procvar].dataFiles
    453                     break
--> 454     fields[var] = Field.from_netcdf(paths, (var, name), dims, inds, grid=grid, mesh=mesh, timestamps=timestamps,
    455                                     allow_time_extrapolation=allow_time_extrapolation,
    456                                     time_periodic=time_periodic, deferred_load=deferred_load,
    457                                     fieldtype=fieldtype, chunksize=varchunksize, dataFiles=dFiles, **kwargs)
    459 u = fields.pop('U', None)
    460 v = fields.pop('V', None)

File c:\\Users\\leeka\\anaconda3\\envs\\parcels\\Lib\\site-packages\\parcels\\field.py:464, in Field.from_netcdf(cls, filenames, variable, dimensions, indices, grid, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, **kwargs)
    459     raise RuntimeError('Multiple files given but no time dimension specified')
    461 if grid is None:
    462     # Concatenate time variable to determine overall dimension
    463     # across multiple files
--> 464     time, time_origin, timeslices, dataFiles = cls.collect_timeslices(timestamps, data_filenames,
    465                                                                       _grid_fb_class, dimensions,
    466                                                                       indices, netcdf_engine, netcdf_decodewarning)
    467     grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
    468     grid.timeslices = timeslices

File c:\\Users\\leeka\\anaconda3\\envs\\parcels\\Lib\\site-packages\\parcels\\field.py:305, in Field.collect_timeslices(timestamps, data_filenames, _grid_fb_class, dimensions, indices, netcdf_engine, netcdf_decodewarning)
    302 time = time_origin.reltime(time)
    304 if not np.all((time[1:] - time[:-1]) > 0):
--> 305     id_not_ordered = np.where(time[1:] < time[:-1])[0][0]
    306     raise AssertionError(f'Please make sure your netCDF files are ordered in time. First pair of non-ordered files: {dataFiles[id_not_ordered]}, {dataFiles[id_not_ordered + 1]}')
    307 return time, time_origin, timeslices, dataFiles

IndexError: index 0 is out of bounds for axis 0 with size 0"
}

Then I used from_xarray_dataset to load my fieldset and it had no problem with generating the fieldset. I have also specified the interpolation method as ‘cgrid_velocity’. But the particles were still stuck on the shore so I also tried ‘freeslip’ which did not solve my problem either.

### Create fieldset

mesh_mask = r'E:\CROCO\MD_2020_hmin\croco_grd.nc'

dataset_folder = "E:/CROCO/MD_2020_hmin/updated_files3/"

data_files = sorted(glob.glob(f"{dataset_folder}/*.nc"))

# Open all files as a single dataset
ds = xr.open_mfdataset(data_files, combine='nested', concat_dim='time')

ds = ds.set_coords(['lon_psi', 'lat_psi'])

# Select a single time slice (using the first time step)
lon_2d = ds.lon_psi.isel(time=0).values
lat_2d = ds.lat_psi.isel(time=0).values

# Ensure these are 2D arrays
assert lon_2d.ndim == 2 and lat_2d.ndim == 2, "Longitude and Latitude should be 2D arrays"

# Update the dataset with these 2D arrays
ds['lon_psi_2d'] = (('eta_psi', 'xi_psi'), lon_2d)
ds['lat_psi_2d'] = (('eta_psi', 'xi_psi'), lat_2d)

# Set as coordinates
ds = ds.set_coords(['lon_psi_2d', 'lat_psi_2d'])

# Create FieldSet with C-grid specifications
fieldset = parcels.FieldSet.from_xarray_dataset(ds, 
                                                variables={'U': 'masked_u', 'V': 'masked_v'},
                                                dimensions={'U': {'lat': 'lat_psi_2d', 'lon': 'lon_psi_2d', 'time': 'time'},
                                                            'V': {'lat': 'lat_psi_2d', 'lon': 'lon_psi_2d', 'time': 'time'}},
                                                mesh='spherical',
                                                interp_method={'U': 'cgrid_velocity', 'V': 'cgrid_velocity'},
                                                # interp_method={'U': 'freeslip', 'V': 'freeslip'},
                                                allow_time_extrapolation=True)

And here is how I set my particles and the kernels. For the size of timestep, I have tried 1 minute and 5 minutes and opted for 5 minutes.

pset = parcels.ParticleSet.from_line(
    fieldset=fieldset,
    pclass=parcels.JITParticle,
    size=5000,  # releasing 5000 particles

    # Sanmen Island
    start= (114.6750, 22.4798),  # releasing on a line: the start longitude and latitude
    finish=(114.6750, 22.4108),  # releasing on a line: the end longitude and latitude
)

### Deal with out of boundary particles

from parcels import StatusCode

def DeleteErrorParticle(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()

### Execute trajectory calculation and save the results into a new zarr file

output_file = pset.ParticleFile(
    name="E:/CROCO/particles_file/test2_dt5.zarr", outputdt = timedelta(days=1)
)
pset.execute(
    [parcels.AdvectionRK4, DeleteErrorParticle],
    runtime = timedelta(days=370),
    dt = timedelta(minutes=5),
    output_file=output_file,
)

Thanks for the help and please let me know if I need to provide any further information.

Kay