# CROCO 3D tutorial#

This tutorial will show how to run a 3D simulation with output from the CROCO model.

## Example setup#

We start with loading the relevant modules and the data.

```
[1]:
```

```
import os
import parcels
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
example_dataset_folder = parcels.download_example_dataset("CROCOidealized_data")
file = os.path.join(example_dataset_folder, "CROCO_idealized.nc")
```

The simulation will be a very simple, idealised flow: a purely zonal flow over a sloping bottom. This flow (which is somewhat unrealistic of course) nicely showcases that particles stay on their initial depth levels, even though the sigma-layers slope down.

This flow has been created by first running the example from the Shelf front example on the CROCO website. Then, we took the restart file are manually replaced all the `u`

-velocities with `1`

m/s and all the `v`

-velocities with `0`

m/s. This way we get a purely zonal flow. We then started a new simulation from the restart file, and CROCO then automatically calculated the `w`

velocities to
match the new zonal field. We saved the `time=0`

snapshot from this new run and use it below.

Now we create a FieldSet object using the `FieldSet.from_croco()`

method. Note that CROCO is a C-grid (with similar indexing at MITgcm), so we need to provide the longitudes and latitudes of the \(\rho\)-points of the grid (`lon_rho`

and `lat_rho`

). We also need to provide the sigma levels at the depth points (`s_w`

). Finally, it is important to also provide the bathymetry field (`h`

), which is needed to convert the depth levels of the particles to sigma-coordinates.

**Note** that in the code below we use the `w`

velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an `omega`

field, which may be more appropriate to use. The idealised simulation below only works when using `w`

, though. In other simulations, it is recommended to test whether `omega`

provides more realistic results. See OceanParcels/Parcels#1728 for more information.

```
[2]:
```

```
variables = {"U": "u", "V": "v", "W": "w", "H": "h"}
lon_rho = "x_rho" # Note, this would be "lon_rho" for a dataset on a spherical grid
lat_rho = "y_rho" # Note ,this would be "lat_rho" for a dataset on a spherical grid
dimensions = {
"U": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
"V": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
"W": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
"H": {"lon": lon_rho, "lat": lat_rho},
}
fieldset = parcels.FieldSet.from_croco(
file,
variables,
dimensions,
allow_time_extrapolation=True, # Note, this is only needed for this specific example dataset, that has only one snapshot
mesh="flat", # Note, this is only needed for this specific example dataset, that has been created on a 'flat' mesh (i.e. in km instead of in degrees)
)
```

```
/var/folders/x0/1qxj679n55zcybylvdsb4bxh0000gq/T/ipykernel_99569/3687013997.py:12: FieldSetWarning: Note that it is unclear which vertical velocity ('w' or 'omega') to use in 3D croco fields.
See https://docs.oceanparcels.org/en/latest/examples/tutorial_croco_3D.html for more information
fieldset = parcels.FieldSet.from_croco(
```

Now we can use this Fieldset to advect particles as we would normally do. Note that the particle depths should be provided in (negative) meters, not in sigma-coordinates.

```
[3]:
```

```
X, Z = np.meshgrid(
[40e3, 80e3, 120e3],
[100, -10, -130, -250, -400, -850, -1400, -1550],
)
Y = np.ones(X.size) * fieldset.U.grid.lat[25]
def DeleteParticle(particle, fieldset, time):
if particle.state >= 50:
particle.delete()
pset = parcels.ParticleSet(
fieldset=fieldset, pclass=parcels.JITParticle, lon=X, lat=Y, depth=Z
)
outputfile = pset.ParticleFile(name="croco_particles3D.zarr", outputdt=5000)
pset.execute(
[parcels.AdvectionRK4_3D, DeleteParticle],
runtime=5e4,
dt=100,
output_file=outputfile,
)
```

```
INFO: Output files are stored in croco_particles3D.zarr.
100%|██████████| 50000.0/50000.0 [00:00<00:00, 79485.48it/s]
```

Now we plot the particle trajectories below. Note that the particles stay on their initial depth levels, even though the sigma-layers slope down. Also note that particles released above the surface (where depth >0) or below the bathymetry are not advected (due to the `DeleteParticle`

kernel).

```
[4]:
```

```
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ds = xr.open_zarr("croco_particles3D.zarr")
ax.plot(ds.lon.T / 1e3, ds.z.T, ".-")
dsCROCO = xr.open_dataset(file)
for z in dsCROCO.s_w.values:
ax.plot(fieldset.H.lon / 1e3, fieldset.H.data[0, 25, :] * z, "k", linewidth=0.5)
ax.set_xlabel("X [km]")
ax.set_xlim(30, 170)
ax.set_ylabel("Depth [m]")
ax.set_title("Particles in idealized CROCO velocity field using 3D advection")
plt.tight_layout()
plt.show()
```

### A CROCO simulation with no vertical velocities#

It may be insightful to compare this 3D run with the `AdvectionRK4_3D`

kernel with a run where the vertical velocity (`W`

) is set to zero. In that case, the particles will not stay on their initial depth levels but instead follow sigma-layers.

```
[5]:
```

```
import copy
fieldset_noW = copy.copy(fieldset)
fieldset_noW.W.data[:] = 0.0
pset_noW = parcels.ParticleSet(
fieldset=fieldset_noW, pclass=parcels.JITParticle, lon=X, lat=Y, depth=Z
)
outputfile = pset.ParticleFile(name="croco_particles_noW.zarr", outputdt=5000)
pset_noW.execute(
[parcels.AdvectionRK4_3D, DeleteParticle],
runtime=5e4,
dt=100,
output_file=outputfile,
)
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ds = xr.open_zarr("croco_particles_noW.zarr")
ax.plot(ds.lon.T / 1e3, ds.z.T, ".-")
dsCROCO = xr.open_dataset(file)
for z in dsCROCO.s_w.values:
ax.plot(fieldset.H.lon / 1e3, fieldset.H.data[0, 25, :] * z, "k", linewidth=0.5)
ax.set_xlabel("X [km]")
ax.set_xlim(30, 170)
ax.set_ylabel("Depth [m]")
ax.set_title("Particles in idealized CROCO velocity field with W=0")
plt.tight_layout()
plt.show()
```

```
INFO: Output files are stored in croco_particles_noW.zarr.
100%|██████████| 50000.0/50000.0 [00:00<00:00, 82825.22it/s]
```

## The algorithms used#

When using `FieldSet.from_croco()`

, Parcels knows that depth needs to be converted to sigma-coordinates, before doing any interpolation. This is done under the hood, using code for interpolation (in this case a `T`

Field) like

```
sigma = particle.depth / fieldset.H[time, particle.depth, particle.lat, particle.lon]
temp = fieldset.T[time, sigma, particle.lat, particle.lon]
```

For the `AdvectionRK4_3D`

kernel, Parcels will replace the kernel with `AdvectionRK4_3D_CROCO`

, which works slightly different from the normal 3D advection kernel because it converts the vertical velocity in sigma-units.

In particular, the following algorithm is used (note that the RK4 version is slightly more complex than this Euler-Forward version, but the idea is identical)

```
# calculate local sigma level of particle, by scaling depth by local ocean depth H
sigma = particle.depth / fieldset.H[time, particle.depth, particle.lat, particle.lon]
(u, v, w) = fieldset.UVW[time, particle.depth, particle.lat, particle.lon, particle]
# scaling the w with the sigma level of the particle
w_sigma = w * sigma / fieldset.H[time, particle.depth, particle.lat, particle.lon]
lon_new = particle.lon + u*particle.dt
lat_new = particle.lat + v*particle.dt
# calculating new sigma level
sigma_new = sigma + w_sigma*particle.dt
# Converting back from sigma to depth, at _new_ location
depth_new = sigma_new * fieldset.H[time, particle.depth, lat_new, lon_new]
```