Public Documentation

<!– Title –> <h1 align="center"> ClimaOcean.jl </h1>

<!– description –> <p align="center"> <strong>🌎 A framework for realistic ocean-only and coupled ocean + sea-ice simulations driven by prescribed atmospheres and based on <a href=https://github.com/CliMA/Oceananigans.jl>Oceananigans</a> and <a href=https://github.com/CliMA/ClimaSeaIce.jl>ClimaSeaIce</a></strong>. </p>

Installation

ClimaOcean is a registered package. To install from a Julia REPL:

julia> using Pkg

julia> Pkg.add("ClimaOcean")

julia> Pkg.instantiate()

Use Pkg.add("url=https://github.com/CliMA/ClimaOcean.jl.git", rev="main") to install the latest version of ClimaOcean. For more information, see the documentation for Pkg.jl.

Why? What's the difference between ClimaOcean and Oceananigans?

Oceananigans is a general-purpose library for ocean-flavored fluid dynamics. ClimaOcean implements a framework for driving realistic Oceananigans simulations with prescribed atmospheres, and coupling them to prognostic sea ice simulations.

A core abstraction: ClimaOcean.OceanSeaIceModel

Our system for realistic modeling is anchored by ClimaOcean.OceanSeaIceModel, which encapsulates the ocean simulation, sea ice simulation, prescribed atmospheric state, and specifies how the three communicate. To illustrate how OceanSeaIceModel works we set up a simulation on a grid with 10 vertical levels and 1/4-degree horizontal resolution:

using Oceananigans
using Oceananigans.Units
using Dates, CFTime
import ClimaOcean

arch = GPU()
grid = LatitudeLongitudeGrid(arch,
                             size = (1440, 560, 10),
                             halo = (7, 7, 7),
                             longitude = (0, 360),
                             latitude = (-70, 70),
                             z = (-3000, 0))

bathymetry = ClimaOcean.regrid_bathymetry(grid) # builds gridded bathymetry based on ETOPO1
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bathymetry))

# Build an ocean simulation initialized to the ECCO state estimate on Jan 1, 1993
ocean = ClimaOcean.ocean_simulation(grid)
date  = DateTimeProlepticGregorian(1993, 1, 1)
set!(ocean.model, T = ClimaOcean.ECCOMetadata(:temperature; date),
                  S = ClimaOcean.ECCOMetadata(:salinity; date))

# Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis
atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch)
coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere)
simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days)
run!(simulation)

The simulation above achieves approximately 8 simulated years per day of wall time on an Nvidia H100 GPU.

Since ocean.model is an Oceananigans.HydrostaticFreeSurfaceModel, we can leverage Oceananigans features in our scripts. For example, to plot the surface speed at the end of the simulation we write

u, v, w = ocean.model.velocities
speed = Field(sqrt(u^2 + v^2))
compute!(speed)

using GLMakie
heatmap(view(speed, :, :, ocean.model.grid.Nz), colorrange=(0, 0.5), colormap=:magma, nan_color=:lightgray)

which produces

Additional features: a utility for ocean_simulations and data wrangling

A second core abstraction in ClimaOcean is ocean_simulation. ocean_simulation configures an Oceananigans model for realistic simulations including temperature and salinity, the TEOS-10 equation of state, boundary conditions to store computed air-sea fluxes, the automatically-calibrated turbulence closure CATKEVerticalDiffusivity, and the WENOVectorInvariant advection scheme for mesoscale-turbulence-resolving simulations.

ClimaOcean also provides convenience features for wrangling datasets of bathymetry, ocean temperature, salinity, ocean velocity fields, and prescribed atmospheric states.

ClimaOcean is built on top of Oceananigans and ClimaSeaIce, so it's important that ClimaOcean users become proficient with Oceananigans. Note that though ClimaOcean is currently focused on hydrostatic modeling with Oceananigans.HydrostaticFreeSurfaceModel, realistic nonhydrostatic modeling is also within the scope of this package.

@distribute communicator for i in iterable
    ...
end

Distribute a for loop among different ranks in communicator. If communicator is not provided, MPI.COMM_WORLD is used.

@handshake communicator exs...

perform exs on all ranks in communicator, but only one rank at a time, where ranks r2 > r1 wait for rank r1 to finish before executing exs. If communicator is not provided, MPI.COMM_WORLD is used.

@onrank communicator rank exs...

Perform exp only on rank rank (0-based index) in communicator. Other ranks will wait for rank rank to finish before continuing. The expression is run anyways if MPI in not initialized. If communicator is not provided, MPI.COMM_WORLD is used.

@root communicator exs...

Perform exs only on rank 0 in communicator, otherwise known as the "root" rank. Other ranks will wait for the root rank to finish before continuing. If communicator is not provided, MPI.COMM_WORLD is used.

MixedLayerDepthField(grid, tracers, buoyancy_model; Δb = 3e-4, field_kw...)

Return a reduced Field{Center, Center, Nothing} that represents mixed layer depth for model, based on a buoyancy differential criterion. The mixed layer depth is defined as the depth $h$ for which

\[b(z=0) - b(z=-h) = Δb\]

This criterion is solved by integrating downwards and linearly interpolating to find h, assuming that $b$ decreases with depth.

Keyword arguments

• Δb: Buoyancy differential used to calculate mixed layer depth
• field_kw: Keyword arguments passed to Field.

Example

h = MixedLayerDepth(model)
compute!(h) # compute mixed layer depth

struct ECCOMetadata{D, V}

Metadata information for an ECCO dataset:

• name: The name of the dataset.
• dates: The dates of the dataset, in an AbstractCFDateTime format.
• version: The version of the dataset, could be ECCO2Monthly, ECCO2Daily, or ECCO4Monthly.
• dir: The directory where the dataset is stored.

ECCOMetadata(name::Symbol; 
             dates = DateTimeProlepticGregorian(1993, 1, 1),
             version = ECCO4Monthly(),
             dir = download_ECCO_cache)

Construct an ECCOMetadata object with the specified parameters.

Arguments

• name::Symbol: The name of the metadata.

Keyword Arguments

• dates: The date(s) of the metadata. Note this can either be a single date, representing a snapshot, or a range of dates, representing a time-series. Default: DateTimeProlepticGregorian(1993, 1, 1).

• version: The data version. Supported versions are ECCO2Monthly(), ECCO2Daily(), or ECCO4Monthly().

• dir: The directory of the data file. Default: download_ECCO_cache.

ECCORestoring(variable_name::Symbol, [ arch_or_grid = CPU(), ];
              version = ECCO4Monthly(),
              dates = all_ECCO_dates(version),
              time_indices_in_memory = 2, 
              time_indexing = Cyclical(),
              mask = 1,
              rate = 1,
              inpainting = NearestNeighborInpainting(Inf),
              cache_inpainted_data = true)

Build a forcing term that restores to values stored in an ECCO field time series. The restoring is applied as a forcing on the right hand side of the evolution equations calculated as

\[Fψ = r μ (ψ_{ECCO} - ψ)\]

where $μ$ is the mask, $r$ is the restoring rate, $ψ$ is the simulation variable, and the ECCO variable $ψ_ECCO$ is linearly interpolated in space and time from the ECCO dataset of choice to the simulation grid and time.

Arguments

• variable_name: The name of the variable to restore. Choices include: * :temperature, * :salinity, * :u_velocity, * :v_velocity, * :sea_ice_thickness, * :sea_ice_area_fraction.

                Note that `ECCOMetadata` may be provided as the first argument instead
                of `variable_name`. In this case the `version` and `dates` kwargs (described below)
                cannot be provided.

• arch_or_grid: Either the architecture of the simulation, or a grid on which the ECCO data is pre-interpolated when loaded. If an architecture is provided, such as arch_or_grid = CPU() or arch_or_grid = GPU(), ECCO data will be interpolated on-the-fly when the forcing tendency is computed. Default: CPU().

Keyword Arguments

• version: The version of the ECCO dataset. Default: ECCO4Monthly().

• dates: The dates to use for the ECCO dataset. Default: all_ECCO_dates(version).

• time_indices_in_memory: The number of time indices to keep in memory. The number is chosen based on a trade-off between increased performance (more indices in memory) and reduced memory footprint (fewer indices in memory). Default: 2.

• time_indexing: The time indexing scheme for the field time series.

• mask: The mask value. Can be a function of (x, y, z, time), an array, or a number.

• rate: The restoring rate, i.e., the inverse of the restoring timescale (in s⁻¹).

• inpainting: inpainting algorithm, see inpaint_mask!. Default: NearestNeighborInpainting(Inf).

• cache_inpainted_data: If true, the data is cached to disk after inpainting for later retrieving. Default: true.

LinearlyTaperedPolarMask(; northern = (70,   75),
                           southern = (-75, -70),
                           z = (-20, 0))

Build a mask that is linearly tapered in latitude between the northern and southern edges. The mask is constant in depth between the z and equals zero everywhere else. The mask is limited to lie between (0, 1). The mask has the following functional form:

n = 1 / (northern[2] - northern[1]) * (φ - northern[1])
s = 1 / (southern[1] - southern[2]) * (φ - southern[2])

valid_depth = (z[1] < z < z[2])

mask = valid_depth ? clamp(max(n, s), 0, 1) : 0

ECCO_field(metadata::ECCOMetadata;
           architecture = CPU(),
           inpainting = nothing,
           mask = nothing,
           horizontal_halo = (7, 7),
           cache_inpainted_data = false)

Return a Field on architecture described by ECCOMetadata with horizontal_halo size. If not nothing, the inpainting method is used to fill the cells within the specified mask. mask is set to ECCO_mask for non-nothing inpainting.

ECCO_immersed_grid(metadata, architecture = CPU())

Compute the ImmersedBoundaryGrid for metadata with a bottom height field that is defined by the first non-missing value from the bottom up.

ECCO_mask(architecture = CPU(); minimum_value = Float32(-1e5))

A boolean field where true represents a missing value in the ECCO dataset.

regrid_bathymetry(target_grid;
                  url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", 
                  height_above_water = <none>,
                  minimum_depth = 0,
                  dir = download_cache,
                  filename = "ETOPO_2022_v1_60s_N90W180_surface.nc")

Regrid bathymetry associated with the NetCDF file at path = joinpath(dir, filename) to target_grid. If path does not exist, then a download is attempted from joinpath(url, filename).

Arguments

• target_grid: grid to interpolate onto

Keyword Arguments

• height_above_water: limits the maximum height of above-water topography (where h > 0). If nothing the original topography is retained

• minimum_depth: minimum depth for the shallow regions, defined as a positive value. h > - minimum_depth will be considered land

• dir: directory of the bathymetry-containing file

• filename: file containing bathymetric data. Must be netcdf with fields: (1) lat vector of latitude nodes (2) lon vector of longitude nodes (3) z matrix of depth values

• interpolation_passes: regridding/interpolation passes. The bathymetry is interpolated in interpolation_passes - 1 intermediate steps. With more steps the final bathymetry will be smoother.

  Example: interpolating from a 400x200 grid to a 100x100 grid in 4 passes involves:

  • 400x200 → 325x175
  • 325x175 → 250x150
  • 250x150 → 175x125
  • 175x125 → 100x100

  If coarsening the original grid, linear interpolation in passes is equivalent to applying a smoothing filter, with more passes increasing the strength of the filter. If refining the original grid, additional passes will not help and no intermediate steps will be performed.

• major_basins: Number of "independent major basins", or fluid regions fully encompassed by land, that are retained by remove_minor_basins!. Basins are removed by order of size: the smallest basins are removed first. major_basins=1 will retain only the largest basin. Default: Inf, which does not remove any basins.

retrieve_bathymetry(grid, filename; kw...)

Retrieve the bathymetry data from a file or generate it using a grid and save it to a file.

Arguments

• grid: The grid used to generate the bathymetry data.
• filename: The name of the file to read or save the bathymetry data.
• kw...: Additional keyword arguments.

Returns

• bottom_height: The retrieved or generated bathymetry data.

If the specified file exists, the function reads the bathymetry data from the file. Otherwise, it generates the bathymetry data using the provided grid and saves it to the file before returning it.