Let's start with importing the relevant modules. The key ones are all in the `parcels`

package.

```
In [1]:
```%matplotlib inline
from parcels import FieldSet, ParticleSet, Variable, JITParticle, AdvectionRK4, plotTrajectoriesFile
import numpy as np
import math
from datetime import timedelta
from operator import attrgetter

`FieldSet`

object, which is simply a collection of hydrodynamic fields. In this first case, we use a simple flow of two idealised moving eddies. That field is saved in NetCDF format in the directory `examples/MovingEddies_data`

. Since we know that the files are in what's called Parcels FieldSet format, we can call these files using the function `FieldSet.from_parcels()`

.

```
In [2]:
```fieldset = FieldSet.from_parcels("MovingEddies_data/moving_eddies")

`fieldset`

can then be visualised with the `show()`

function. To show the zonal velocity (`U`

), give the following command

```
In [3]:
```fieldset.U.show()

```
```

`ParticleSet`

. In this case, we start 2 particles at locations (330km, 100km) and (330km, 280km) using the `from_list`

constructor method, that are advected on the `fieldset`

we defined above. Note that we use `JITParticle`

as `pclass`

, because we will be executing the advection in JIT (Just-In-Time) mode. The alternative is to run in `scipy`

mode, in which case `pclass`

is `ScipyParticle`

```
In [4]:
```pset = ParticleSet.from_list(fieldset=fieldset, # the fields on which the particles are advected
pclass=JITParticle, # the type of particles (JITParticle or ScipyParticle)
lon=[3.3e5, 3.3e5], # a vector of release longitudes
lat=[1e5, 2.8e5]) # a vector of release latitudes

Print the `ParticleSet`

to see where they start

```
In [5]:
```print(pset)

```
```

`not_yet_set`

, that is because we didn't specify a `time`

when we defined the `pset`

.

To plot the positions of these particles on the zonal velocity, use the following command

```
In [6]:
```pset.show(field=fieldset.U)

```
```

`ParticelSet`

. We run the particles using the `AdvectionRK4`

kernel, which is a 4th order Runge-Kutte implementation that comes with Parcels. We run the particles for 6 days (using the `timedelta`

function from `datetime`

), at an RK4 timestep of 5 minutes. We store the trajectory information at an interval of 1 hour in a file called `EddyParticles.nc`

. Because `time`

was `not_yet_set`

, the particles will be advected from the first date available in the `fieldset`

, which is the default behaviour.

```
In [7]:
```pset.execute(AdvectionRK4, # the kernel (which defines how particles move)
runtime=timedelta(days=6), # the total length of the run
dt=timedelta(minutes=5), # the timestep of the kernel
output_file=pset.ParticleFile(name="EddyParticles.nc", outputdt=timedelta(hours=1))) # the file name and the time step of the outputs

```
```

The code should have run, which can be confirmed by printing and plotting the `ParticleSet`

again

```
In [8]:
```print(pset)
pset.show(field=fieldset.U)

```
```

`U`

field have moved in the plot above. Also, the `time`

of the particles is now 518400 seconds, which is 6 days.

`EddyParticles.nc`

file. It can be quickly plotted using the `plotTrajectoriesFile`

function.

```
In [9]:
```plotTrajectoriesFile('EddyParticles.nc');

```
```

`plotTrajectoriesFile`

function can also be used to show the trajectories as an animation, by specifying that it has to run in `movie2d_notebook`

mode. If we pass this to our function above, we can watch the particles go!

```
In [10]:
```plotTrajectoriesFile('EddyParticles.nc', mode='movie2d_notebook')

```
Out[10]:
```

`plotTrajectoriesFile`

can also be used to display 2-dimensional histograms (`mode=hist2d`

) of the number of particle observations per bin. Use the `bins`

argument to control the number of bins in the longitude and latitude direction. See also the matplotlib.hist2d page.

```
In [11]:
```plotTrajectoriesFile('EddyParticles.nc', mode='hist2d', bins=[30, 20]);

```
```

*during execution*, which is great for debugging. To rerun the particles while plotting them on top of the zonal velocity field (`fieldset.U`

), first reinitialise the `ParticleSet`

and then re-execute. However, now rather than saving the output to a file, display a movie using the `moviedt`

display frequency, in this case with the zonal velocity `fieldset.U`

as background

```
In [12]:
```# THIS DOES NOT WORK IN THIS IPYTHON NOTEBOOK, BECAUSE OF THE INLINE PLOTTING.
# THE 'SHOW_MOVIE' KEYWORD WILL WORK ON MOST MACHINES, THOUGH
# pset = ParticleSet(fieldset=fieldset, size=2, pclass=JITParticle, lon=[3.3e5, 3.3e5], lat=[1e5, 2.8e5])
# pset.execute(AdvectionRK4,
# runtime=timedelta(days=6),
# dt=timedelta(minutes=5),
# moviedt=timedelta(hours=1),
# movie_background_field=fieldset.U)

Running particles in backward time is extremely simple: just provide a `dt`

< 0.

```
In [13]:
```pset.execute(AdvectionRK4,
dt=-timedelta(minutes=5), # negative timestep for backward run
runtime=timedelta(days=6), # the run time
output_file=pset.ParticleFile(name="EddyParticles_Bwd.nc", outputdt=timedelta(hours=1))) # the file name and the time step of the outputs

```
In [14]:
```print(pset)
pset.show(field=fieldset.U)

```
```

```
In [15]:
```def WestVel(particle, fieldset, time):
if time > 86400:
uvel = -2.
particle.lon += uvel * particle.dt

`ParticleSet`

again, and re-execute. Note that we have now changed `kernel`

to be `AdvectionRK4 + k_WestVel`

, where `k_WestVel`

is the `WestVel`

function as defined above cast into a `Kernel`

object (via the `pset.Kernel`

call).

```
In [16]:
```pset = ParticleSet.from_list(fieldset=fieldset, pclass=JITParticle, lon=[3.3e5, 3.3e5], lat=[1e5, 2.8e5])
k_WestVel = pset.Kernel(WestVel) # casting the WestVel function to a kernel object
pset.execute(AdvectionRK4 + k_WestVel, # simply add kernels using the + operator
runtime=timedelta(days=2),
dt=timedelta(minutes=5),
output_file=pset.ParticleFile(name="EddyParticles_WestVel.nc", outputdt=timedelta(hours=1)))

```
```

And now plot this new trajectory file

```
In [17]:
```plotTrajectoriesFile('EddyParticles_WestVel.nc');

```
```

`FieldSet.from_netcdf()`

function.

The `examples`

directory contains a set of GlobCurrent files of the region around South Africa.

`*`

) and the filenames for U and V can be the same (as in this case)

```
In [18]:
```filenames = {'U': "GlobCurrent_example_data/20*.nc",
'V': "GlobCurrent_example_data/20*.nc"}

`U`

and `V`

) and dimensions (`lon`

, `lat`

and `time`

; note that in this case there is no `depth`

because the GlobCurrent data is only for the surface of the ocean)

```
In [19]:
```variables = {'U': 'eastward_eulerian_current_velocity',
'V': 'northward_eulerian_current_velocity'}
dimensions = {'lat': 'lat',
'lon': 'lon',
'time': 'time'}

`FieldSet.from_netcdf`

function with the above-defined `filenames`

, `variables`

and `dimensions`

```
In [20]:
```fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)

```
```

`ParticleSet`

, in this case with 5 particle starting on a line between (28E, 33S) and (30E, 33S) using the `ParticleSet.from_line`

constructor method

```
In [21]:
```pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle,
size=5, # releasing 5 particles
start=(28, -33), # releasing on a line: the start longitude and latitude
finish=(30, -33)) # releasing on a line: the end longitude and latitude

And finally execute the `ParticleSet`

for 10 days using 4th order Runge-Kutta

```
In [22]:
```pset.execute(AdvectionRK4,
runtime=timedelta(days=10),
dt=timedelta(minutes=5),
output_file=pset.ParticleFile(name="GlobCurrentParticles.nc", outputdt=timedelta(hours=6)))

```
```

`plotParticles`

script again. Note you can plot the particles on top of one of the velocity fields using the `tracerfile`

, `tracerfield`

, etc keywords.

```
In [23]:
```plotTrajectoriesFile('GlobCurrentParticles.nc',
tracerfile='GlobCurrent_example_data/20020101000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc',
tracerlon='lon',
tracerlat='lat',
tracerfield='eastward_eulerian_current_velocity');

```
```

`P`

) field, using `extra_fields={'P': 'P'}`

. Note that, because this flow does not depend on time, we need to set `allow_time_extrapolation=True`

when reading in the fieldset.

```
In [24]:
```fieldset = FieldSet.from_parcels("Peninsula_data/peninsula", extra_fields={'P': 'P'}, allow_time_extrapolation=True)

`Particle`

class that has an extra `Variable`

: the pressure. We initialise this by sampling the `fieldset.P`

field.

```
In [25]:
```class SampleParticle(JITParticle): # Define a new particle class
p = Variable('p', initial=fieldset.P) # Variable 'p' initialised by sampling the pressure

`ParticleSet`

using the `from_line`

method also used above in the GlobCurrent data. Plot the `pset`

and print their pressure values `p`

```
In [26]:
```pset = ParticleSet.from_line(fieldset=fieldset, pclass=SampleParticle,
start=(3000, 3000), finish=(3000, 46000), size=5, time=0)
pset.show(field='vector')
print('p values before execution:', [p.p for p in pset])

```
```

`fieldset.P`

field at the particle location. Cast this function to a `Kernel`

.

```
In [27]:
```def SampleP(particle, fieldset, time): # Custom function that samples fieldset.P at particle location
particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
k_sample = pset.Kernel(SampleP) # Casting the SampleP function to a kernel.

`pset`

with a combination of the `AdvectionRK4`

and `SampleP`

kernels, plot the `pset`

and print their new pressure values `p`

```
In [28]:
```pset.execute(AdvectionRK4 + k_sample, # Add kernels using the + operator.
runtime=timedelta(hours=20),
dt=timedelta(minutes=5))
pset.show(field=fieldset.P, show_time=0)
print('p values after execution:', [p.p for p in pset])

```
```

`p`

are (within roundoff errors) the same as the pressure values before the execution of the kernels. The particles thus stay on isobars!

`Particle`

class that includes three extra variables. The `distance`

variable will be written to output, but the auxiliary variables `prev_lon`

and `prev_lat`

won't be written to output (can be controlled using the `to_write`

keyword)

```
In [29]:
```class DistParticle(JITParticle): # Define a new particle class that contains three extra variables
distance = Variable('distance', initial=0., dtype=np.float32) # the distance travelled
prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
initial=attrgetter('lon')) # the previous longitude
prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
initial=attrgetter('lat')) # the previous latitude.

`TotalDistance`

that calculates the sum of Euclidean distances between the old and new locations in each RK4 step

```
In [30]:
```def TotalDistance(particle, fieldset, time):
# Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)
lat_dist = (particle.lat - particle.prev_lat) * 1.11e2
# Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth
lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi / 180)
# Calculate the total Euclidean distance travelled by the particle
particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2))
particle.prev_lon = particle.lon # Set the stored values for next iteration.
particle.prev_lat = particle.lat

*Note:* here it is assumed that the latitude and longitude are measured in degrees North and East, respectively. However, some datasets (e.g. the `MovingEddies`

used above) give them measured in (kilo)meters, in which case we must *not* include the factor `1.11e2`

.

`TotalDistance`

function on a `ParticleSet`

containing the five particles within the `GlobCurrent`

fieldset from above. Note that `pclass=DistParticle`

in this case.

```
In [31]:
```filenames = {'U': "GlobCurrent_example_data/20*.nc",
'V': "GlobCurrent_example_data/20*.nc"}
variables = {'U': 'eastward_eulerian_current_velocity',
'V': 'northward_eulerian_current_velocity'}
dimensions = {'lat': 'lat',
'lon': 'lon',
'time': 'time'}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
pset = ParticleSet.from_line(fieldset=fieldset,
pclass=DistParticle,
size=5, start=(28, -33), finish=(30, -33))

Again define a new kernel to include the function written above and execute the `ParticleSet`

.

```
In [32]:
```k_dist = pset.Kernel(TotalDistance) # Casting the TotalDistance function to a kernel.
pset.execute(AdvectionRK4 + k_dist, # Add kernels using the + operator.
runtime=timedelta(days=6),
dt=timedelta(minutes=5),
output_file=pset.ParticleFile(name="GlobCurrentParticles_Dist.nc", outputdt=timedelta(hours=1)))

```
```

`EddyParticles_Dist.nc`

file)

```
In [33]:
```print([p.distance for p in pset]) #the distances in km travelled by the particles

```
```