Parcels also supports curvilinear grids such as those used in the NEMO models.
However, there is one important additional step that needs to be taken before running Parcels on these curvilear grids: create a file of the 'angles' between the grid direction and the zonal (and meridional) direction. This tutorial explains how to do just that.
We will be using the example data in the NemoCurvilinear_data/
directory. These fiels are a purely zonal flow on an aqua-planet (so zonal-velocity is 1 m/s and meridional-velocity is 0 m/s everywhere, and no land). However, because of the curvilinear grid, the U
and V
fields vary north of 26N.
First, import the usual modules, as well as the compute_curvilinearGrid_rotationAngles
module.
from parcels import FieldSet, ParticleSet, JITParticle, ParticleFile, plotTrajectoriesFile
from parcels import compute_curvilinearGrid_rotationAngles
from parcels import AdvectionRK4
import numpy as np
from datetime import timedelta as delta
Now the first step is to create a new NetCDF file with these angle fields, based on the mesh_mask.nc4
file:
mesh_filename = 'NemoCurvilinear_data/mesh_mask.nc4'
rotation_angles_filename = 'NemoCurvilinear_data/rotation_angles.nc'
compute_curvilinearGrid_rotationAngles(mesh_filename, rotation_angles_filename)
INFO: Generating rotation angles fields in file: NemoCurvilinear_data/rotation_angles.nc
Note that the compute_curvilinearGrid_rotationangles()
function assumes some default values for the dimensions of U
, V
and F
. If your mesh_file has different dimensions, you can set these with a dictionary such as the default one below, and add that as dimensions=dimensions
to the call of compute_curvilinearGrid_rotationangles()
.
dimensions = {'U': {'lon': 'glamu', 'lat': 'gphiu'},
'V': {'lon': 'glamv', 'lat': 'gphiv'},
'F': {'lon': 'glamf', 'lat': 'gphif'}}
Now, we can create a FieldSet
just like we do for normal grids. Note that we need to provide extra variables for the rotation angles fields that we just computed
filenames = {'U': 'NemoCurvilinear_data/U_purely_zonal-ORCA025_grid_U.nc4',
'V': 'NemoCurvilinear_data/V_purely_zonal-ORCA025_grid_V.nc4',
'cosU': rotation_angles_filename,
'sinU': rotation_angles_filename,
'cosV': rotation_angles_filename,
'sinV': rotation_angles_filename}
variables = {'U': 'U',
'V': 'V',
'cosU': 'cosU',
'sinU': 'sinU',
'cosV': 'cosV',
'sinV': 'sinV'}
dimensions = {'U': {'lon': 'nav_lon_u', 'lat': 'nav_lat_u'},
'V': {'lon': 'nav_lon_v', 'lat': 'nav_lat_v'},
'cosU': {'lon': 'glamu', 'lat': 'gphiu'},
'sinU': {'lon': 'glamu', 'lat': 'gphiu'},
'cosV': {'lon': 'glamv', 'lat': 'gphiv'},
'sinV': {'lon': 'glamv', 'lat': 'gphiv'}}
field_set = FieldSet.from_netcdf(filenames, variables, dimensions, mesh='spherical', allow_time_extrapolation=True)
WARNING: Casting lon data to np.float32 WARNING: Casting lat data to np.float32 WARNING: Casting time data to np.float64 WARNING: Casting depth data to np.float32
And we can plot the U
field. Note that we need to rescale the grid.lon
so that it is monotonic (needed for plotting purposes only; Parcels itself can handle non-monotonic grids)
field_set.U.grid.lon[field_set.U.grid.lon<73] += 360
field_set.U.show()
As you see above, the U
field indeed is 1 m/s south of 50N, but varies with longitude and latitude north of that
Now we can run particles as normal. Parcels will take care to rotate the U
and V
fields
# Start 20 particles on a meridional line at 30E
npart = 20
lonp = 30 * np.ones(npart)
latp = [i for i in np.linspace(-70, 88, npart)]
# Create a periodic boundary condition kernel
def periodicBC(particle, pieldSet, time, dt):
if particle.lon > 180:
particle.lon -= 360
pset = ParticleSet.from_list(field_set, JITParticle, lon=lonp, lat=latp)
pfile = ParticleFile("nemo_particles", pset)
kernels = pset.Kernel(AdvectionRK4) + periodicBC
pset.execute(kernels, runtime=delta(days=50), dt=delta(hours=6),
interval=delta(days=1), output_file=pfile)
INFO: Compiled JITParticleAdvectionRK4periodicBC ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/c77947732f5b8acf06a97ad5cd7c1995.so
And then we can plot these trajectories. As expected, all trajectories go exactly zonal and due to the curvature of the earth, ones at higher latitude move more degrees eastward (even though the distance in km is equal for all particles)
plotTrajectoriesFile("nemo_particles.nc")