🖥️ MOM6 tutorial#
This tutorial shows how to run a surface (2D) simulation with output from the MOM6 ocean model. MOM6 uses a curvilinear Arakawa C-grid: U is defined on the eastern cell faces (yh, xq), V on the northern cell faces (yq, xh), and the grid geometry is given by the cell corners geolon_c/geolat_c on (yq, xq). Parcels has built-in functionality to deal with C-grid data, so the velocities can be used on the native grid without regridding or rotating them. See also the NEMO tutorial and MITgcm tutorial for other examples of C-grid data.
Note
The MOM6-to-SGRID converter currently only works with 2-dimensional (surface) velocity fields. If you have a dataset with 3-dimensional velocity fields, please raise an Issue on the Parcels GitHub repository, specifying the details of your dataset and then we can provide guidance on how to handle it.
The example data is a small subset of the CEFI Northeast Pacific (NEP) MOM6 hindcast.
import matplotlib.pyplot as plt
import numpy as np
import parcels
import parcels.tutorial
U = parcels.tutorial.open_dataset("MOM6_example_data/U")["ssu"]
V = parcels.tutorial.open_dataset("MOM6_example_data/V")["ssv"]
coords = parcels.tutorial.open_dataset("MOM6_example_data/static")
We can use a combination of parcels.convert.mom6_to_sgrid() and parcels.FieldSet.from_sgrid_conventions() to read in the data. The converter needs the velocity fields and the MOM6 static grid file (ocean_static.nc), which contains the corner coordinates.
MOM6 can write output in a symmetric layout (one more corner point than cell centres along each axis, as in this regional dataset) or a non-symmetric layout (equal sizes). mom6_to_sgrid() detects the layout from the grid sizes and sets the SGRID padding accordingly. The time coordinate must be in ascending order, so we sort it first.
ds_fset = parcels.convert.mom6_to_sgrid(
fields={"U": U.sortby("time"), "V": V.sortby("time")}, coords=coords
)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
Our example MOM6 grid is curvilinear and rotated relative to longitude/latitude, we select release points from the model grid itself (every few corner points, away from the domain edges) so that all particles start inside the domain.
lon, lat = fieldset.UV.grid.lon, fieldset.UV.grid.lat # 2D, shape (yq, xq)
iy = np.linspace(2, lon.shape[0] - 3, 10).astype(int)
ix = np.linspace(2, lon.shape[1] - 3, 10).astype(int)
X = lon[np.ix_(iy, ix)]
Y = lat[np.ix_(iy, ix)]
Z = np.zeros_like(X)
def DeleteParticle(particles, fieldset):
any_error = particles.state >= 50
particles[any_error].state = parcels.StatusCode.Delete
pset = parcels.ParticleSet(fieldset=fieldset, x=X, y=Y, z=Z)
outputfile = parcels.ParticleFile(
path="mom6_particles.parquet",
outputdt=np.timedelta64(3, "h"),
)
pset.execute(
[parcels.kernels.AdvectionRK2, DeleteParticle],
runtime=np.timedelta64(2, "D"),
dt=np.timedelta64(30, "m"),
output_file=outputfile,
)
And make a simple plot of the particle trajectories, together with the grid corner points to show the curvilinear model grid:
df = parcels.read_particlefile("mom6_particles.parquet")
fig, ax = plt.subplots()
ax.plot(coords["geolon_c"], coords["geolat_c"], "k.", markersize=1, alpha=0.3)
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"], traj["y"], ".-")
ax.set_xlabel("Longitude [°]")
ax.set_ylabel("Latitude [°]")
plt.show()
