SpeedyWeather.jl

SpeedyWeather.jl is a global atmospheric climate model designed for research with a focus on speed, interactivity, and extensibility. It is built in Julia to allow users to interactively run, extend, and visualize models in the REPL, notebooks (Jupyter/Pluto), or via scripts.

Key Features

Dynamics and Physics

• Supports various physical equations: barotropic vorticity, shallow water, and primitive equations (with or without humidity).
• Includes 2D particle advection and customizable 2D/3D tracer advection.
• Physics parameterizations for convection, precipitation, boundary layer, radiation, etc.
• Includes a slab ocean, thermodynamic sea ice, and a 2-layer land bucket model.

Numerics and Computing

• Grids: Full and octahedral, Gaussian and Clenshaw-Curtis, HEALPix, and OctaHEALPix.
• Resolutions: T31 to T1023 and higher (approx. 400km to 10km).
• Arithmetics: Supports Float32 (default), Float64, and experimental BFloat16 with stochastic rounding.
• Transforms: Uses the SpeedyTransforms library for fast spherical harmonics transforms.

User Interface and Extensibility

• Visualization: 2D and 3D interactive visualization powered by Makie.
• Modularity: Models are constructed from components; non-default components are passed as arguments.
• Callbacks: Users can inject code after every time step for custom output or event handling.

In SpeedyWeather, atmospheric prognostic equations are scaled with the planet's radius R. Most tendencies (parameterizations, forcing, and drag) are automatically scaled by the dynamical core. If you are adding a custom term to force temperature, you can simply write it in units of Kelvin per second without manual scaling.

When manual scaling is required: You must manually apply scaling only when forcing the vorticity or divergence equations directly (in spectral space) and the forcing is not proportional to the variable being forced.

If the forcing term is proportional to the variable (e.g., $\propto \zeta^n$ where $n=1$), the dynamical core handles the scaling automatically. If the forcing is non-linear or constant, you must apply manual scaling with $R^{1-n}$:

• Constant forcing ($n=0$): Scale the forcing manually with $R$.
• Quadratic forcing ($n=2$): Scale the forcing manually with $1/R$.
• General case: Scale with $R^{1-n}$.

You can write multiple ensemble members into a single Zarr store by adding an ensemble dimension. This is designed for parallel ensemble members running as separate processes.

Each process should construct a ZarrOutput with the same path, id, run_number, and filename, but a unique ensemble_index.

Options:

• ensemble_index::Int: The index of this member (starting at 1). Set to 0 (default) to disable ensemble output.
• ensemble_size::Int: Total number of ensemble members. Must satisfy ensemble_size ≥ ensemble_index.
• ensemble_timeout::Int: Seconds a member waits for member 1 to create the shared store before erroring (default 600).

Note: This feature is currently specific to ZarrOutput and is not supported by NetCDFOutput.

using SpeedyWeather, Zarr

# In a parallel process environment:
# member is the current process's ensemble index (1..ensemble_size)
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
output = ZarrOutput(spectral_grid, PrimitiveWet;
    ensemble_index = member,
    ensemble_size = 10,
    interval = Hour(6),
)
model = PrimitiveWetModel(spectral_grid; output)
simulation = initialize!(model)
run!(simulation, period=Day(10), output=true)

For high-performance scenarios where you interpolate between the same two grids repeatedly, or interpolate onto the same number of points with changing locations, reuse an interpolator object to avoid recomputing indices and weights.

Reusing an Interpolator between two grids

1. Create the interpolator: interp = RingGrids.interpolator(grid_out, grid_in).
2. Use interpolate! with the precomputed interpolator: interpolate!(field_out, field_in, interp).

Reusing an Interpolator for moving coordinates

If you have a fixed number of points but their locations change, use an AnvilInterpolator:

1. Initialize: interp = AnvilInterpolator(input_grid, npoints, NF=Float32).
2. Update locations: RingGrids.update_locator!(interp, lons, lats).
3. Interpolate: interpolate!(output_vec, field_in, interp).

Precision Control

You can specify the interpolation precision using the NF keyword argument (e.g., NF=Float32) when creating the interpolator. This allows you to perform interpolation in higher or lower precision than the input/output data types.

# Reusing interpolator for two grids
grid_in = HEALPixGrid(4)
grid_out = FullClenshawGrid(6)
interp = RingGrids.interpolator(grid_out, grid_in)

field_in = rand(grid_in)
field_out = zeros(grid_out)
interpolate!(field_out, field_in, interp)

# High precision interpolation (Float32) with Float16 data
field_in = randn(OctahedralGaussianField{Float16}, 24)
field_out = zeros(FullClenshawField{Float16}, 24)
interp = RingGrids.interpolator(field_out, field_in, NF=Float32)
interpolate!(field_out, field_in, interp)

# Moving coordinates with AnvilInterpolator
npoints = 10
grid = HEALPixGrid(24)
interp = AnvilInterpolator(grid, npoints, NF=Float32)

londs = collect(-10.0:2.0:8.0)
latds = collect(0.0:5.0:45.0)
RingGrids.update_locator!(interp, londs, latds)

output_vec = zeros(10)
field_in = rand(grid)
interpolate!(output_vec, field_in, interp)

SpeedyWeather represents surface fluxes as the exchange of momentum, heat, and humidity/moisture between the ocean/land and the lowermost atmospheric layer.

Key implementation details:

• Momentum Fluxes: Represent drag experienced by the boundary layer wind.
• Heat Fluxes: Represent sensible heat flux between the surface and the atmosphere.
• Humidity Fluxes: Represent evaporation (ocean) or evapotranspiration/condensation (land).

Important: Interdependence of surface flux computations By default, model.boundary_layer.surface_condition extrapolates atmospheric variables (wind, temperature, humidity) to the surface. The model.boundary_layer also centrally computes a drag coefficient used for momentum, heat, and humidity fluxes.

Warning: Setting model.boundary_layer.surface_condition = nothing or model.boundary_layer = nothing will disable all other surface fluxes unless they have use_boundary_layer_drag = false (which allows them to use their own independent drag coefficient).

using SpeedyWeather
# To see available momentum flux implementations:
subtypes(SpeedyWeather.AbstractSurfaceMomentumFlux)

# To see available heat flux implementations:
subtypes(SpeedyWeather.AbstractSurfaceHeatFlux)

# To see available humidity flux implementations:
subtypes(SpeedyWeather.AbstractSurfaceHumidityFlux)

# To see available boundary layer implementations (handles drag/extrapolation):
subtypes(SpeedyWeather.AbstractBoundaryLayer)