🖥️ AIFS tutorial#
Atmospheric dispersion studies are frequently driven by reanalysis or forecast products such as ERA5 or the (A)IFS. These products are often interpolated towards an A-grid, on which the U, V and W velocity components are defined on the same nodes. In this tutorial we advect particles with data from the AIFS forecast, a data-driven forecast product delivered by the European Centre for Medium-Range Weather Forecasts (ECMWF) and openly available from the ECMWF open data portal. Here we use a version that has been regridded onto an A-grid and provided by the MeteoCloud of the Jülich Supercomputing Centre. The workflow described below also applies to reanalysis products such as ERA5.
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import parcels
import parcels.tutorial
Building a meteorological dataset#
Below, we use the pre-computed AIFS data to build a meteorological dataset for Parcels. AIFS output is provided every six hours in separate NetCDF files. Note that in this dataset, we have already stacked these files along a new time dimension and assign each one the corresponding timestamp. You may need to do this yourself if you are using a different dataset. The xarray library can help you with this task.
ds = parcels.tutorial.open_dataset("AIFSmeteo/data")
Adapting the dataset dimensions#
The raw NetCDF dimensions (latitude, longitude and pressure) need a few adjustments before they can be used with Parcels. First, the latitude axis is stored in descending order, so we reverse both the coordinate and the data along it. Second, the pressure levels are stored in descending order as well; we reverse them so that pressure increases monotonically from the top of the atmosphere down to the surface. Pressure is also converted from Pa to hPa.
In 1644 Torricelli wrote “We live submerged at the bottom of an ocean of air”. Because pressure decreases with height and increases as we descend towards the surface, it plays a similiar role as depth in an ocean model. Therefore we use pressure as the vertical (“depth”) coordinate of the atmosphere. Concretely, this simply means renaming the plev coordinate to depth.
With pressure as the vertical coordinate, the vertical velocity is the pressure tendency \(\omega\).
Note
The AIFS data used here is provided on pressure levels, but reanalysis and forecast products are also available on other vertical coordinates with other vertical velocities. This has to be taken into account when preparing other datasets for atmospheric Parcels simulations.
# Reverse the latitude axis (and its data) to ascending order
ds = ds.isel(lat=slice(None, None, -1))
# Reverse the pressure levels so that pressure increases downwards
ds = ds.isel(plev=slice(None, None, -1))
# Convert pressure from Pa to hPa and the pressure tendency from Pa/s to hPa/s
ds["w"] = ds["w"] / 100.0
ds["plev"] = ds["plev"] / 100.0
# Rename plev to depth, since pressure acts as the vertical coordinate
ds = ds.rename({"plev": "depth"})
# Rename u, v, w to match Parcels' expected variable names
ds = ds.rename({"u": "U", "v": "V", "w": "W"})
Adding grid metadata#
Next we add the grid faces and the remaining properties that Parcels needs to interpret the AIFS data. The pressure, longitude and latitude faces are computed as the midpoints between neighbouring full levels.
We then declare the dataset conventions as CF-1.6, SGRID-0.3 and tag each coordinate with the axis it represents (X, Y, Z or T), so that Parcels knows how the coordinates map onto the model dimensions. From this information we build the grid topology, attach it to the u, v and w variables, and mark those velocities as being located on the nodes.
# Compute grid faces as the midpoints between neighbouring full levels
lon_face = 0.5 * (ds["lon"].values[:-1] + ds["lon"].values[1:])
ds["lon_face"] = xr.DataArray(lon_face, dims=["lon_face"])
lat_face = 0.5 * (ds["lat"].values[:-1] + ds["lat"].values[1:])
ds["lat_face"] = xr.DataArray(lat_face, dims=["lat_face"])
depth_face = 0.5 * (ds["depth"].values[:-1] + ds["depth"].values[1:])
ds["depth_face"] = xr.DataArray(depth_face, dims=["depth_face"])
# Declare the dataset conventions and the axis of each coordinate
ds.attrs["Conventions"] = "CF-1.6, SGRID-0.3"
ds["lon"].attrs["axis"] = "X"
ds["lat"].attrs["axis"] = "Y"
ds["depth"].attrs["axis"] = "Z"
ds["time"].attrs["axis"] = "T"
# Describe the grid topology
ds["grid"] = xr.DataArray(
data=np.int32(0),
attrs={
"cf_role": "grid_topology",
"topology_dimension": 2,
"node_dimensions": "lon lat",
"node_coordinates": "lon lat",
"face_dimensions": "lon_face: lon (padding: none) lat_face: lat (padding: none)",
"face_coordinates": "lon_face lat_face",
"vertical_dimensions": "depth_face: depth (padding: none)",
},
)
The resulting dataset looks as follows:
display(ds)
<xarray.Dataset> Size: 42MB
Dimensions: (time: 5, depth: 13, lat: 201, lon: 201, lon_face: 200,
lat_face: 200, depth_face: 12)
Coordinates:
* time (time) datetime64[ns] 40B 2025-07-01 ... 2025-07-02
* depth (depth) float64 104B 50.0 100.0 150.0 ... 850.0 925.0 1e+03
* lat (lat) float64 2kB -70.0 -69.75 -69.5 ... -20.5 -20.25 -20.0
* lon (lon) float64 2kB 25.0 25.25 25.5 25.75 ... 74.5 74.75 75.0
* lon_face (lon_face) float64 2kB 25.12 25.38 25.62 ... 74.38 74.62 74.88
* lat_face (lat_face) float64 2kB -69.88 -69.62 -69.38 ... -20.38 -20.12
* depth_face (depth_face) float64 96B 75.0 125.0 175.0 ... 775.0 887.5 962.5
Data variables:
U (time, depth, lat, lon) float32 11MB dask.array<chunksize=(1, 13, 201, 201), meta=np.ndarray>
t (time, depth, lat, lon) float32 11MB dask.array<chunksize=(1, 13, 201, 201), meta=np.ndarray>
W (time, depth, lat, lon) float32 11MB dask.array<chunksize=(1, 13, 201, 201), meta=np.ndarray>
V (time, depth, lat, lon) float32 11MB dask.array<chunksize=(1, 13, 201, 201), meta=np.ndarray>
grid int32 4B 0
Attributes:
CDI: Climate Data Interface version 2.4.4 (https://mpimet.mpg.de...
Conventions: CF-1.6, SGRID-0.3
institution: European Centre for Medium-Range Weather Forecasts
history: Mon Sep 14 16:55:01 2026: cdo -s -O -f nc4 selname,u,v,w,t ...
CDO: Climate Data Operators version 2.4.4 (https://mpimet.mpg.de...Creating the FieldSet#
With the dataset prepared, we can build the parcels.FieldSet using parcels.FieldSet.from_sgrid_conventions(). We pass mesh="spherical" because the velocities are given in local (spherical) coordinates. For the interpolation of the 2D and 3D fields we choose linear interpolation, which is the method most commonly used in atmospheric Lagrangian transport models.
We also speed up the I/O by using the ChunkCachedArray backend - see the squeezing performance out of Parcels tutorial for more details.
# Build the FieldSet and use linear interpolation
fieldset = parcels.FieldSet.from_sgrid_conventions(ds)
fieldset.to_chunk_cached_arrays()
fieldset.describe()
| Name | Type | Grid number | Interp method / value | Parcels backend |
|:-------|:------------|--------------:|:------------------------|:------------------|
| U | Field | 0 | XLinear(...) | ChunkCachedArray |
| V | Field | 0 | XLinear(...) | ChunkCachedArray |
| W | Field | 0 | XLinear(...) | ChunkCachedArray |
| t | Field | 0 | XLinear(...) | ChunkCachedArray |
| UV | VectorField | 0 | XLinear_Velocity(...) | - |
| UVW | VectorField | 0 | XLinear_Velocity(...) | - |
mesh: SphericalMesh(radius=6366707.019493707)
time interval: (np.datetime64('2025-07-01T00:00:00.000000000'), np.datetime64('2025-07-02T00:00:00.000000000'))
Creating the initial particles#
We initialize the particles on a three-dimensional grid with NumPy’s meshgrid, releasing them at three pressure levels (900 hPa, 500 hPa and 100 hPa) over a small region of interest. The parcels.ParticleSet is then created with Parcels’ own constructor.
pressure, lat, lon = np.meshgrid(
np.array([900, 500, 100]),
np.linspace(-65, -25, 10),
np.linspace(30, 70, 10),
indexing="ij",
)
pset = parcels.ParticleSet(fieldset=fieldset, x=lon, y=lat, z=pressure)
Running the simulation#
The simulation is run with ParticleSet.execute, which applies a list of kernels over a given runtime (here one day) with a fixed integration time step (here 900 s).
First we need a three-dimensional advection kernel; we use Runge-Kutta 2 through parcels.kernels.AdvectionRK2_3D. Second, because we simulate trajectories within a regional lat-lon box, we define a kernel that removes particles leaving the domain. The DeleteOutOfBounds kernel below deletes any parcel whose status is ErrorOutOfBounds. For global simulations one would instead add a periodic boundary kernel that wraps longitude at the dateline and handle the lower and upper pressure boundaries, for example by reflecting particles or letting them slide along the maximum and minimum pressure levels.
We store the output in a ParticleFile; here we write a Parquet file with an output frequency of 15 minutes.
# Define a kernel that deletes particles that leave the domain
def DeleteOutOfBounds(particles, fieldset):
"""Delete particles that have left the domain boundaries."""
out_of_bounds = particles.state == parcels.StatusCode.ErrorOutOfBounds
particles[out_of_bounds].state = parcels.StatusCode.Delete
outfile = parcels.ParticleFile(
"atmosphere_parcels.parquet", outputdt=np.timedelta64(15, "m"), mode="w"
)
pset.execute(
[parcels.kernels.AdvectionRK2_3D, DeleteOutOfBounds],
runtime=np.timedelta64(1, "D"),
dt=np.timedelta64(900, "s"),
output_file=outfile,
)
Plotting the trajectories#
Finally, we study the trajectories with a Matplotlib plotting routine. We use parcels.read_particlefile() to load the output and plot the particle positions in three dimensions for the full day. The trajectories show an anticyclone at the surface and fast zonal winds in the higher layers.
from mpl_toolkits.mplot3d import Axes3D
df = parcels.read_particlefile("atmosphere_parcels.parquet")
fig = plt.figure(figsize=(13, 8))
ax = plt.axes(projection="3d")
ax.view_init(azim=-145, elev=25)
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot3D(traj["x"], traj["y"], traj["z"], ".-")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_zlabel("Pressure (mbar)")
ax.set_zlim(df["z"].max(), 0)
plt.show()
