Global climate simulation

This example configures a global ocean–sea ice simulation at 1ᵒ horizontal resolution with realistic bathymetry and a few closures including the "Gent-McWilliams" IsopycnalSkewSymmetricDiffusivity. The atmosphere is represented by a 4-layer SpeedyWeather simulation on the T63 spectral grid (this grid has approximately 1.875ᵒ resolution).

The atmosphere is initialized with the Jablonowski and Williamson (2006) initial conditions, which consist of a zonal wind centered at mid-latitudes and higher altitudes and a temperature profile that is baroclinically unstable. The surface pressure is adjusted by topography for approximately globally constant mean-sea level pressure. The initial specific humidity is calculated from temperature for a constant relative humidity everywhere. The ocean and sea ice are initialized by ocean temperature, salinity, sea ice concentration, and sea ice thickness from the ECCO state estimate.

For this example, we need Oceananigans.HydrostaticFreeSurfaceModel (the ocean), ClimaSeaIce.SeaIceModel (the sea ice) and SpeedyWeather.PrimitiveWetModel (the atmosphere). All these are coupled and orchestrated by the NumericalEarth.EarthSystemModel (the coupled system).

The ConservativeRegridding.jl package is used to regrid fields between the atmosphere and ocean–sea ice components. Both the ocean–sea ice and atmosphere components run on a CUDA-enabled GPU.

using Oceananigans, SpeedyWeather, NumericalEarth, ConservativeRegridding
using CUDA
using NCDatasets, CairoMakie
using Oceananigans.Units
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Grids: φnode
using Printf, Statistics, Dates

[ Info: Oceananigans will use 8 threads

Ocean and sea-ice model configuration

The ocean and sea-ice are a simplified versions of the one-degree ocean-sea ice example.

The first step is to create the grid with realistic bathymetry.

Nx = 360
Ny = 180
Nz = 60
z = ExponentialDiscretization(Nz, -6000, 0; mutable=true)
grid = TripolarGrid(Oceananigans.GPU(); size=(Nx, Ny, Nz), z, halo=(5, 5, 4))

We regrid the bathymetry.

bottom_height = regrid_bathymetry(grid; minimum_depth=10, major_basins=2, interpolation_passes=10)
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)

[ Info: Loading cached field from /storage5/github-action-runners/julia-depot/scratchspaces/904d977b-046a-4731-8b86-9235c0d1ef02/field_cache/field_bottom_height_360x180x60_0f93908a.jld2

Now we can specify the numerical details and the closures for the ocean simulation. We use a biharmonic horizontal viscosity with a timescale of 15 days and set a background vertical diffusivity following Henyey et al. (1986).

momentum_advection = WENOVectorInvariant(order=5)
tracer_advection   = WENO(order=5)
free_surface       = SplitExplicitFreeSurface(grid; substeps=70)
catke_closure      = NumericalEarth.Oceans.default_ocean_closure()
eddy_closure       = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=500, κ_symmetric=200)

@inline νhb(i, j, k, grid, timescale) = Oceananigans.Operators.Azᶜᶜᶜ(i, j, k, grid)^2 / timescale
ν = Oceananigans.Field{Center, Center, Center}(grid)
Oceananigans.set!(ν, KernelFunctionOperation{Center, Center, Center}(νhb, grid, 15days))
horizontal_viscosity = HorizontalScalarBiharmonicDiffusivity(; ν)

νz = 1e-5
# κz(φ) = max(2e-6, 3e-5 * |sin(φ)|)
@inline henyey_diffusivity(i, j, k, grid) = max(2e-6, 3e-5 * abs(sind(φnode(i, j, k, grid, Center(), Center(), Center()))))
κz = Oceananigans.Field{Center, Center, Center}(grid)
Oceananigans.set!(κz, KernelFunctionOperation{Center, Center, Center}(henyey_diffusivity, grid))
vertical_diffusivity = VerticalScalarDiffusivity(ν=νz, κ=κz)

closures = (catke_closure, eddy_closure, horizontal_viscosity, vertical_diffusivity)

The ocean simulation, complete with initial conditions for temperature and salinity from ECCO on Jan 1st, 1992.

ocean = ocean_simulation(grid;
                         momentum_advection,
                         tracer_advection,
                         free_surface,
                         closure = closures)

ecco_set = MetadataSet(:temperature, :salinity,
                       :sea_ice_thickness, :sea_ice_concentration;
                       dataset = ECCO4Monthly(),
                       date = DateTime(1992, 1, 1))
Oceananigans.set!(ocean.model, ecco_set)   # T, S

The sea-ice simulation, complete with initial conditions for sea-ice thickness and sea-ice concentration from ECCO.

sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=5))

Oceananigans.set!(sea_ice.model, ecco_set)   # h, ℵ

Atmosphere model configuration

The atmosphere is provided by SpeedyWeather.jl. Here, we configure a T63L4 model with a 3-hour output interval. The atmosphere_simulation function takes care of building an atmosphere model with appropriate hooks so that NumericalEarth can compute inter-component fluxes.

nlayers = 4
output_interval = 3hours
stop_time = 30days
spectral_grid = SpeedyWeather.SpectralGrid(; NF=Float64, truncation=64, nlayers, Grid=FullClenshawGrid, architecture=SpeedyWeather.GPU())
time_stepping = SpeedyWeather.Leapfrog(spectral_grid; Δt_at_T32=Minute(40)) # gives Δt = 20 min at truncation=64
atmosphere = atmosphere_simulation(spectral_grid; output_interval, time_stepping, stop_time)
atmosphere.model.feedback.verbose = false  # disable SpeedyWeather's progress bar in favor of the callback defined below

Precalculate forward matrix:   3%|▋                      |  ETA: 0:01:09Precalculate forward matrix:  19%|████▎                  |  ETA: 0:00:21Precalculate forward matrix:  35%|████████               |  ETA: 0:00:13Precalculate forward matrix:  51%|███████████▊           |  ETA: 0:00:09Precalculate forward matrix:  67%|███████████████▍       |  ETA: 0:00:05Precalculate forward matrix:  83%|███████████████████    |  ETA: 0:00:03Precalculate forward matrix:  99%|██████████████████████▊|  ETA: 0:00:00Precalculate forward matrix: 100%|███████████████████████| Time: 0:00:15
Precalculate backward matrix:  20%|████▍                 |  ETA: 0:00:08

The atmosphere model already includes some initial conditions as described above:

atmosphere.model.initial_conditions

(vordiv = SpeedyWeather.ZonalWind{Float64} <: SpeedyWeather.AbstractInitialConditions
├ η₀::Float64 = 0.252
├ u₀::Float64 = 35.0
├ perturb_lat::Float64 = 40.0
├ perturb_lon::Float64 = 20.0
├ perturb_uₚ::Float64 = 1.0
└ perturb_radius::Float64 = 0.1, pres = SpeedyWeather.PressureOnOrography <: SpeedyWeather.AbstractInitialConditions
, temp = SpeedyWeather.JablonowskiTemperature{Float64} <: SpeedyWeather.AbstractInitialConditions
├ η₀::Float64 = 0.252
├ σ_tropopause::Float64 = 0.2
├ u₀::Float64 = 35.0
├ ΔT::Float64 = 0.0
└ Tmin::Float64 = 200.0, humid = SpeedyWeather.ConstantRelativeHumidity{Float64} <: SpeedyWeather.AbstractInitialConditions
└ relhumid_ref::Float64 = 0.7)

The coupled model

We are now ready to blend everything together. Here we set the timestep to be the same across all models.

Δt = convert(eltype(grid), atmosphere.model.time_stepping.Δt)

We build the complete coupled earth_model and the coupled simulation. NumericalEarth still computes turbulent (sensible heat, water vapor) fluxes here using its own bulk-formula machinery and writes them back to Speedy. Top-level radiation, however, is not yet wired up against SpeedyWeather's upwelling-LW path, so we leave radiation = nothing for now.

earth_model = EarthSystemModel(; atmosphere, sea_ice, ocean)

EarthSystemModel{GPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: SpeedyWeather.Simulation
├── land: Nothing
├── sea_ice: SeaIceModel{CUDAGPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── ocean: HydrostaticFreeSurfaceModel{CUDAGPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfaces

Building and running the simulation

We are ready to build and run the coupled simulation. But before we do, we add callbacks to write outputs to disk every 3 hours.

earth = Oceananigans.Simulation(earth_model; Δt, stop_time)
outputs = merge(ocean.model.velocities, ocean.model.tracers)
sea_ice_fields = merge(sea_ice.model.velocities, sea_ice.model.dynamics.auxiliaries.fields,
                       (; h=sea_ice.model.ice_thickness, ℵ=sea_ice.model.ice_concentration))

ocean.output_writers[:free_surf] = JLD2Writer(ocean.model, (; η=ocean.model.free_surface.displacement);
                                              overwrite_files=true,
                                              schedule=TimeInterval(output_interval),
                                              filename="ocean_free_surface.jld2")

ocean.output_writers[:surface] = JLD2Writer(ocean.model, outputs;
                                            overwrite_files=true,
                                            schedule=TimeInterval(output_interval),
                                            filename="ocean_surface_fields.jld2",
                                            indices=(:, :, grid.Nz))

sea_ice.output_writers[:fields] = JLD2Writer(sea_ice.model, sea_ice_fields;
                                             overwrite_files=true,
                                             schedule=TimeInterval(output_interval),
                                             filename="sea_ice_fields.jld2")

𝒬ᵀᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat
𝒬ᵛᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat
τˣᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.x_momentum
τʸᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.y_momentum
𝒬ᵀᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat
𝒬ᵛᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat
τˣᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.x_momentum
τʸᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.y_momentum
𝒬ⁱᵒ  = earth.model.interfaces.net_fluxes.sea_ice.bottom.heat
Jˢⁱᵒ  = earth.model.interfaces.sea_ice_ocean_interface.fluxes.salt
fluxes = (; 𝒬ᵀᵃᵒ, 𝒬ᵛᵃᵒ, τˣᵃᵒ, τʸᵃᵒ, 𝒬ᵀᵃⁱ, 𝒬ᵛᵃⁱ, τˣᵃⁱ, τʸᵃⁱ, 𝒬ⁱᵒ, Jˢⁱᵒ)

ocean.output_writers[:fluxes] = JLD2Writer(earth.model.ocean.model, fluxes;
                                           overwrite_files=true,
                                           schedule=TimeInterval(output_interval),
                                           filename="intercomponent_fluxes.jld2")

JLD2Writer scheduled on TimeInterval(3 hours):
├── filepath: intercomponent_fluxes.jld2
├── 10 outputs: (𝒬ᵀᵃᵒ, 𝒬ᵛᵃᵒ, τˣᵃᵒ, τʸᵃᵒ, 𝒬ᵀᵃⁱ, 𝒬ᵛᵃⁱ, τˣᵃⁱ, τʸᵃⁱ, 𝒬ⁱᵒ, Jˢⁱᵒ)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

We also add a callback function that prints out a helpful progress message while the simulation runs.

wall_time = Ref(time_ns())
model_time = Ref(0.0)

function progress(sim)
    atmos = sim.model.atmosphere
    ocean = sim.model.ocean

    ua, va     = atmos.variables.dynamics.u_mean_grid, atmos.variables.dynamics.v_mean_grid
    uo, vo, wo = ocean.model.velocities

    # RingGrids only defines the one-argument `maximum`, so `abs` must be broadcast first to stay on the GPU
    uamax = (maximum(abs.(ua)), maximum(abs.(va)))
    uomax = (maximum(abs, uo), maximum(abs, vo), maximum(abs, wo))

    step_time = 1e-9 * (time_ns() - wall_time[])
    sypd = (time(sim) - model_time[]) / step_time / 365

    msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim))
    msg2 = @sprintf(", max|ua|: (%.1e, %.1e) m s⁻¹", uamax...)
    msg3 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", uomax...)
    msg4 = @sprintf(", wall time: %s", prettytime(step_time))
    msg5 = @sprintf(", SYPD: %.3f \n", sypd)

    @info msg1 * msg2 * msg3 * msg4 * msg5

    wall_time[] = time_ns()
    model_time[] = time(sim)

     return nothing
end

add_callback!(earth, progress, TimeInterval(2days))

Let's run the coupled model!

Oceananigans.run!(earth)

[ Info: Initializing simulation...
[ Info: time: 0 seconds, iter: 0, max|ua|: (0.0e+00, 0.0e+00) m s⁻¹, max|uo|: (0.0e+00, 0.0e+00, 0.0e+00) m s⁻¹, wall time: 1.069 minutes, SYPD: 0.000 
[ Info:     ... simulation initialization complete (17.622 seconds)
[ Info: Executing initial time step...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters /storage5/github-action-runners/julia-depot/packages/Oceananigans/lCP9j/src/OutputWriters/jld2_writer.jl:259
[ Info:     ... initial time step complete (13.764 minutes).
[ Info: time: 2 days, iter: 144, max|ua|: (3.0e+01, 9.1e+00) m s⁻¹, max|uo|: (1.4e+00, 1.4e+00, 1.3e-02) m s⁻¹, wall time: 15.525 minutes, SYPD: 0.508 
[ Info: time: 4 days, iter: 288, max|ua|: (3.5e+01, 1.9e+01) m s⁻¹, max|uo|: (1.6e+00, 1.0e+00, 1.1e-02) m s⁻¹, wall time: 1.781 minutes, SYPD: 4.430 
[ Info: time: 6 days, iter: 432, max|ua|: (3.3e+01, 2.1e+01) m s⁻¹, max|uo|: (1.6e+00, 1.4e+00, 1.1e-02) m s⁻¹, wall time: 1.807 minutes, SYPD: 4.366 
[ Info: time: 8 days, iter: 576, max|ua|: (3.6e+01, 1.9e+01) m s⁻¹, max|uo|: (1.6e+00, 1.1e+00, 1.1e-02) m s⁻¹, wall time: 1.815 minutes, SYPD: 4.346 
[ Info: time: 10 days, iter: 720, max|ua|: (3.8e+01, 1.8e+01) m s⁻¹, max|uo|: (2.0e+00, 1.2e+00, 1.1e-02) m s⁻¹, wall time: 1.826 minutes, SYPD: 4.322 
[ Info: time: 12 days, iter: 864, max|ua|: (3.9e+01, 1.8e+01) m s⁻¹, max|uo|: (1.6e+00, 1.0e+00, 1.1e-02) m s⁻¹, wall time: 1.798 minutes, SYPD: 4.389 
[ Info: time: 14 days, iter: 1008, max|ua|: (4.2e+01, 1.8e+01) m s⁻¹, max|uo|: (1.6e+00, 1.3e+00, 1.1e-02) m s⁻¹, wall time: 1.800 minutes, SYPD: 4.383 
[ Info: time: 16 days, iter: 1152, max|ua|: (4.1e+01, 1.9e+01) m s⁻¹, max|uo|: (1.6e+00, 1.3e+00, 1.1e-02) m s⁻¹, wall time: 1.812 minutes, SYPD: 4.355 
[ Info: time: 18 days, iter: 1296, max|ua|: (4.7e+01, 2.2e+01) m s⁻¹, max|uo|: (1.6e+00, 1.5e+00, 1.1e-02) m s⁻¹, wall time: 1.813 minutes, SYPD: 4.353 
[ Info: time: 20 days, iter: 1440, max|ua|: (3.9e+01, 2.0e+01) m s⁻¹, max|uo|: (1.6e+00, 1.2e+00, 1.1e-02) m s⁻¹, wall time: 1.805 minutes, SYPD: 4.371 
[ Info: time: 22 days, iter: 1584, max|ua|: (4.0e+01, 1.8e+01) m s⁻¹, max|uo|: (1.5e+00, 1.1e+00, 1.0e-02) m s⁻¹, wall time: 1.808 minutes, SYPD: 4.363 
[ Info: time: 24 days, iter: 1728, max|ua|: (4.1e+01, 1.6e+01) m s⁻¹, max|uo|: (1.6e+00, 8.6e-01, 1.0e-02) m s⁻¹, wall time: 1.815 minutes, SYPD: 4.347 
[ Info: time: 26 days, iter: 1872, max|ua|: (4.3e+01, 1.8e+01) m s⁻¹, max|uo|: (1.5e+00, 1.2e+00, 9.6e-03) m s⁻¹, wall time: 1.823 minutes, SYPD: 4.328 
[ Info: time: 28 days, iter: 2016, max|ua|: (4.5e+01, 1.7e+01) m s⁻¹, max|uo|: (1.5e+00, 9.6e-01, 9.2e-03) m s⁻¹, wall time: 1.843 minutes, SYPD: 4.280 
[ Info: Simulation is stopping after running for 41.233 minutes.
[ Info: Simulation time 30 days equals or exceeds stop time 30 days.
[ Info: time: 30 days, iter: 2160, max|ua|: (4.2e+01, 1.8e+01) m s⁻¹, max|uo|: (1.6e+00, 9.6e-01, 9.0e-03) m s⁻¹, wall time: 1.887 minutes, SYPD: 4.181 

Visualizing the results

We can visualize some of the results. Here, we plot the surface speeds in the atmosphere, ocean, and sea-ice as well as the 2-meter temperature in the atmosphere, the sea surface temperature, and the sensible and latent heat fluxes at the atmosphere-ocean interface. SpeedyWeather outputs are stored in a NetCDF file located in the run_0001 folder, while ocean and sea-ice outputs are stored in JLD2 files that can be read by Oceananigans.jl using the FieldTimeSeries type.

SWO = Dataset("run_0001/output.nc")

Ta = reverse(SWO["temp"][:, :, nlayers, :], dims=2)
ua = reverse(SWO["u"][:, :, nlayers, :],    dims=2)
va = reverse(SWO["v"][:, :, nlayers, :],    dims=2)
sp = sqrt.(ua.^2 + va.^2)

SST = FieldTimeSeries("ocean_surface_fields.jld2", "T")
SSU = FieldTimeSeries("ocean_surface_fields.jld2", "u")
SSV = FieldTimeSeries("ocean_surface_fields.jld2", "v")

SIU = FieldTimeSeries("sea_ice_fields.jld2", "u")
SIV = FieldTimeSeries("sea_ice_fields.jld2", "v")
SIA = FieldTimeSeries("sea_ice_fields.jld2", "ℵ")

𝒬ᵀᵃᵒ = FieldTimeSeries("intercomponent_fluxes.jld2", "𝒬ᵀᵃᵒ")
𝒬ᵛᵃᵒ = FieldTimeSeries("intercomponent_fluxes.jld2", "𝒬ᵛᵃᵒ")

Nt = min(length(sp[1, 1, :]), length(𝒬ᵀᵃᵒ))
times = 𝒬ᵀᵃᵒ.times

uotmp = Oceananigans.Field{Face, Center, Nothing}(SST.grid)
votmp = Oceananigans.Field{Center, Face, Nothing}(SST.grid)

uitmp = Oceananigans.Field{Face,   Center, Nothing}(SST.grid)
vitmp = Oceananigans.Field{Center, Face,   Nothing}(SST.grid)
atmp  = Oceananigans.Field{Center, Center, Nothing}(SST.grid)

sotmp = Oceananigans.Field(sqrt(uotmp^2 + votmp^2))
sitmp = Oceananigans.Field(sqrt(uitmp^2 + vitmp^2) * atmp)

iter = Observable(1)

san = @lift sp[:, :, $iter]
son  = @lift begin
    parent(uotmp) .= parent(SSU[$iter])
    parent(votmp) .= parent(SSV[$iter])
    Oceananigans.compute!(sotmp)
    Oceananigans.interior(sotmp, :, :, 1)
end
ssn  = @lift begin
    Oceananigans.set!(uitmp, SIU[$iter])
    Oceananigans.set!(vitmp, SIV[$iter])
    Oceananigans.set!(atmp,  SIA[$iter])
    Oceananigans.compute!(sitmp)
    Oceananigans.interior(sitmp, :, :, 1)
end

fig = Figure(size = (666, 1000))

ax1 = Axis(fig[1, 1], title = "Surface speed, atmosphere")
ax2 = Axis(fig[2, 1], title = "Surface speed, ocean")
ax3 = Axis(fig[3, 1], title = "Surface speed, sea-ice")

hm1 = heatmap!(ax1, san; colormap = :deep,  nan_color=:lightgray, colorrange = (0, 35))
hm2 = heatmap!(ax2, son; colormap = :magma, nan_color=:lightgray, colorrange = (0, 0.6))
hm3 = heatmap!(ax3, ssn; colormap = :ice,   nan_color=:lightgray, colorrange = (0, 0.6))

Colorbar(fig[1, 2], hm1, label="(m s⁻¹)")
Colorbar(fig[2, 2], hm2, label="(m s⁻¹)")
Colorbar(fig[3, 2], hm3, label="(m s⁻¹)")

for ax in (ax1, ax2, ax3)
    hidedecorations!(ax)
end

title = @lift string(prettytime(times[$iter] - times[1]))
Label(fig[0, :], title, fontsize=18)

CairoMakie.record(fig, "surface_speeds.mp4", 1:Nt, framerate=8) do i
    iter[] = i
end

Tan = @lift Ta[:, :, $iter]
Ton = @lift interior(SST[$iter], :, :, 1)
𝒬ᵀn = @lift interior(𝒬ᵀᵃᵒ[$iter], :, :, 1)
𝒬ᵛn = @lift interior(𝒬ᵛᵃᵒ[$iter], :, :, 1)

fig = Figure(size = (700, 1400))

ax1 = Axis(fig[1, 1], title = "2m Temperature, atmosphere")
ax2 = Axis(fig[2, 1], title = "Sea Surface Temperature")
ax3 = Axis(fig[3, 1], title = "Sensible heat flux")
ax4 = Axis(fig[4, 1], title = "Latent heat flux")

hm1 = heatmap!(ax1, Tan; colormap = :plasma, nan_color=:lightgray, colorrange = (-45, 30))
hm2 = heatmap!(ax2, Ton; colormap = :plasma, nan_color=:lightgray, colorrange = (-2, 32))
hm3 = heatmap!(ax3, 𝒬ᵀn; colormap = :balance, colorrange = (-200, 200),  nan_color=:lightgray)
hm4 = heatmap!(ax4, 𝒬ᵛn; colormap = :balance, colorrange = (-200, 200),  nan_color=:lightgray)

Colorbar(fig[1, 2], hm1, label="(ᵒC)")
Colorbar(fig[2, 2], hm2, label="(ᵒC)")
Colorbar(fig[3, 2], hm3, label="(W/m²)")
Colorbar(fig[4, 2], hm4, label="(W/m²)")

for ax in (ax1, ax2, ax3, ax4)
    hidedecorations!(ax)
end

title = @lift string(prettytime(times[$iter] - times[1]))
Label(fig[0, :], title, fontsize=18)

CairoMakie.record(fig, "surface_temperature_and_heat_flux.mp4", 1:Nt, framerate=8) do i
    iter[] = i
end

This page was generated using Literate.jl.