🖥️ Using the `gsw` toolbox to compute density#

This tutorial shows how to use the gsw toolbox (the Gibbs SeaWater Oceanographic Toolbox of TEOS-10) within Parcels to compute density from temperature and salinity fields. The gsw toolbox can be installed via conda install gsw or pixi add gsw.

First, load the necessary libraries and the data:

import gsw  # The package to calculate density from temperature, salinity and depth
import numpy as np

import parcels
import parcels.tutorial

# Load the CopernicusMarine data in the Agulhas region from the example_datasets
ds_fields = parcels.tutorial.open_dataset(
    "CopernicusMarine_data_for_Argo_tutorial/data"
)

# Convert to SGRID-compliant dataset and create FieldSet
fields = {
    "U": ds_fields["uo"],
    "V": ds_fields["vo"],
    "thetao": ds_fields["thetao"],
    "so": ds_fields["so"],
}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

# Convert the FieldSet to windowed arrays for better performance
fieldset = fieldset.to_windowed_arrays()

Now, define a custom Particle class that includes temperature, salinity, and density as variables, and create a ParticleSet with one particle at a known location:

GSWParticle = parcels.Particle.add_variable(
    [
        parcels.Variable("temp", dtype=np.float32, initial=np.nan),
        parcels.Variable("salt", dtype=np.float32, initial=np.nan),
        parcels.Variable("density", dtype=np.float32, initial=np.nan),
    ]
)

# Initiate one Argo float in the Agulhas Current
pset = parcels.ParticleSet(
    fieldset=fieldset,
    pclass=GSWParticle,
    x=[32],
    y=[-31],
    z=[200],
)

Now (as the core part of this tutorial) define a custom kernel that uses the gsw toolbox to compute density from temperature and salinity:

@parcels.validate_kernel
def ParcelsGSW(particles, fieldset):
    particles.temp = fieldset.thetao[particles]
    particles.salt = fieldset.so[particles]
    pressure = gsw.p_from_z(-particles.z, particles.y)
    particles.density = gsw.density.rho(particles.salt, particles.temp, pressure)

Finally, run the ParcelsGSW Kernel for one timestep and check (for Continuous Integration purposes) that the computed density is as expected:

pset.execute(
    ParcelsGSW,
    runtime=np.timedelta64(1, "s"),
    dt=np.timedelta64(1, "s"),
    verbose_progress=False,
)

np.testing.assert_allclose(pset.density, [1026.8281], rtol=1e-5)

print(
    f"Temperature: {pset.temp[0]:.2f}, Salinity: {pset.salt[0]:.2f}, Density: {pset.density[0]:.2f}"
)

Temperature: 16.32, Salinity: 35.52, Density: 1026.83

🖥️ Using the gsw toolbox to compute density#

🖥️ Using the `gsw` toolbox to compute density#