Near-global ocean simulation
This example sets up and runs a near-global ocean simulation using the Oceananigans.jl and ClimaOcean.jl. The simulation covers latitudes from 75°S to 75°N with a horizontal resolution of 1/4 degree and 40 vertical levels.
The simulation's results are visualized with the CairoMakie.jl package.
Initial setup with package imports
We begin by importing the necessary Julia packages for visualization (CairoMakie), ocean modeling (Oceananigans, ClimaOcean), and handling dates and times (CFTime, Dates). These packages provide the foundational tools for setting up the simulation environment, including grid setup, physical processes modeling, and data visualization.
using ClimaOcean
using Oceananigans
using Oceananigans.Units
using CairoMakie
using CFTime
using Dates
using PrintfGrid configuration
We define a global grid with a horizontal resolution of 1/4 degree and 40 vertical levels. The grid is a LatitudeLongitudeGrid spanning latitudes from 75°S to 75°N. We use an exponential vertical spacing to better resolve the upper-ocean layers. The total depth of the domain is set to 6000 meters. Finally, we specify the architecture for the simulation, which in this case is a GPU.
arch = GPU()
Nx = 1440
Ny = 600
Nz = 40
depth = 6000meters
z_faces = exponential_z_faces(; Nz, depth)
grid = LatitudeLongitudeGrid(arch;
size = (Nx, Ny, Nz),
halo = (7, 7, 7),
z = z_faces,
latitude = (-75, 75),
longitude = (0, 360))1440×600×40 LatitudeLongitudeGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.GPU with 7×7×7 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [0.0, 360.0) regularly spaced with Δλ=0.25
├── latitude: Bounded φ ∈ [-75.0, 75.0] regularly spaced with Δφ=0.25
└── z: Bounded z ∈ [-6000.0, 0.0] variably spaced with min(Δz)=8.0258, max(Δz)=645.588Bathymetry and immersed boundary
We use regrid_bathymetry to derive the bottom height from ETOPO1 data. To smooth the interpolated data we use 5 interpolation passes. We also fill in (i) all the minor enclosed basins except the 3 largest major_basins, as well as (ii) regions that are shallower than minimum_depth.
bottom_height = regrid_bathymetry(grid;
minimum_depth = 10meters,
interpolation_passes = 5,
major_basins = 3)
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)1440×600×40 ImmersedBoundaryGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.GPU with 7×7×7 halo:
├── immersed_boundary: GridFittedBottom(mean(z)=-2513.5, min(z)=-6000.0, max(z)=0.0)
├── underlying_grid: 1440×600×40 LatitudeLongitudeGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.GPU with 7×7×7 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [0.0, 360.0) regularly spaced with Δλ=0.25
├── latitude: Bounded φ ∈ [-75.0, 75.0] regularly spaced with Δφ=0.25
└── z: Bounded z ∈ [-6000.0, 0.0] variably spaced with min(Δz)=8.0258, max(Δz)=645.588Let's see what the bathymetry looks like:
h = grid.immersed_boundary.bottom_height
fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-depth, 0))
Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false)
save("bathymetry.png", fig)
Ocean model configuration
We build our ocean model using ocean_simulation,
ocean = ocean_simulation(grid)Simulation of HydrostaticFreeSurfaceModel{GPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: Inf days
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entrieswhich uses the default ocean.model,
ocean.modelHydrostaticFreeSurfaceModel{GPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 1440×600×40 ImmersedBoundaryGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.GPU with 7×7×7 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S, e)
├── closure: CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
├── free surface: Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurfaces.SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── substepping: FixedTimeStepSize(20.254 seconds)
├── advection scheme:
│ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction
│ ├── T: FluxFormAdvection(x=WENO(order=7), y=WENO(order=7), z=Centered(order=2))
│ ├── S: FluxFormAdvection(x=WENO(order=7), y=WENO(order=7), z=Centered(order=2))
│ └── e: Nothing
└── coriolis: Oceananigans.Coriolis.HydrostaticSphericalCoriolis{Oceananigans.Coriolis.ActiveCellEnstrophyConserving, Float64}We initialize the ocean model with ECCO2 temperature and salinity for January 1, 1993.
date = DateTimeProlepticGregorian(1993, 1, 1)
set!(ocean.model, T=ECCOMetadata(:temperature; dates=date),
S=ECCOMetadata(:salinity; dates=date))Prescribed atmosphere and radiation
Next we build a prescribed atmosphere state and radiation model, which will drive the ocean simulation. We use the default Radiation model,
# The radiation model specifies an ocean albedo emissivity to compute the net radiative
# fluxes. The default ocean albedo is based on Payne (1982) and depends on cloud cover
# (calculated from the ratio of maximum possible incident solar radiation to actual
# incident solar radiation) and latitude. The ocean emissivity is set to 0.97.
radiation = Radiation(arch)Radiation{Float64}:
├── stefan_boltzmann_constant: 5.67e-8
├── emission: SurfaceProperties
│ ├── ocean: 0.97
│ └── sea_ice: 0.97
└── reflection: SurfaceProperties
├── ocean: LatitudeDepedendentAlbedo{Float64}: 0.069 - 0.011 ⋅ cos(2φ)
└── sea_ice: 0.7The atmospheric data is prescribed using the JRA55 dataset. The JRA55 dataset provides atmospheric data such as temperature, humidity, and winds to calculate turbulent fluxes using bulk formulae, see CrossRealmFluxes. The number of snapshots that are loaded into memory is determined by the backend. Here, we load 41 snapshots at a time into memory.
atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))640×320×1×2920 PrescribedAtmosphere{Float32} on Oceananigans.Grids.LatitudeLongitudeGrid:
├── times: 2920-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
├── reference_height: 10.0
└── boundary_layer_height: 600.0The coupled simulation
Next we assemble the ocean, atmosphere, and radiation into a coupled model,
coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)OceanSeaIceModel{GPU}(time = 0 seconds, iteration = 0)
├── ocean: HydrostaticFreeSurfaceModel{GPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── atmosphere: 640×320×1×2920 PrescribedAtmosphere{Float32}
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
└── fluxes: OceanSeaIceSurfaceFluxesWe then create a coupled simulation. We start with a small-ish time step of 90 seconds. We run the simulation for 10 days with this small-ish time step.
simulation = Simulation(coupled_model; Δt=90, stop_time=10days)Simulation of OceanSeaIceModel{GPU}(time = 0 seconds, iteration = 0)
├── Next time step: 1.500 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 10 days
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u_ocean on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesWe define a callback function to monitor the simulation's progress,
wall_time = Ref(time_ns())
function progress(sim)
ocean = sim.model.ocean
u, v, w = ocean.model.velocities
T = ocean.model.tracers.T
Tmax = maximum(interior(T))
Tmin = minimum(interior(T))
umax = (maximum(abs, interior(u)),
maximum(abs, interior(v)),
maximum(abs, interior(w)))
step_time = 1e-9 * (time_ns() - wall_time[])
msg = @sprintf("Iter: %d, time: %s, Δt: %s", iteration(sim), prettytime(sim), prettytime(sim.Δt))
msg *= @sprintf(", max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s",
umax..., Tmax, Tmin, prettytime(step_time))
@info msg
wall_time[] = time_ns()
end
simulation.callbacks[:progress] = Callback(progress, TimeInterval(5days))Callback of progress on TimeInterval(5 days)Set up output writers
We define output writers to save the simulation data at regular intervals. In this case, we save the surface fluxes and surface fields at a relatively high frequency (every day). The indices keyword argument allows us to save only a slice of the three dimensional variable. Below, we use indices to save only the values of the variables at the surface, which corresponds to k = grid.Nz
outputs = merge(ocean.model.tracers, ocean.model.velocities)
ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs;
schedule = TimeInterval(1days),
filename = "near_global_surface_fields",
indices = (:, :, grid.Nz),
with_halos = true,
overwrite_existing = true,
array_type = Array{Float32})JLD2OutputWriter scheduled on TimeInterval(1 day):
├── filepath: near_global_surface_fields.jld2
├── 6 outputs: (T, S, e, u, v, w)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 13.8 MiBSpinning up the simulation
We spin up the simulation with a small-ish time-step to resolve the "initialization shock" associated with starting from ECCO2 initial conditions that are both interpolated and also satisfy a different dynamical balance than our simulation.
run!(simulation)[ Info: Initializing simulation...
[ Info: Iter: 0, time: 0 seconds, Δt: 1.500 minutes, max|u|: (0.00e+00, 0.00e+00, 0.00e+00) m s⁻¹, extrema(T): (31.25, -1.87) ᵒC, wall time: 5.411 minutes
[ Info: ... simulation initialization complete (6.931 seconds)
[ Info: Executing initial time step...
┌ Warning: Simulation stopped during initialization.
└ @ Oceananigans.Simulations /central/scratch/esm/slurm-buildkite/climaocean-examples/1664/depot/default/packages/Oceananigans/Y6Lmd/src/Simulations/run.jl:135
[ Info: ... initial time step complete (2.945 minutes).
[ Info: Iter: 4800, time: 5 days, Δt: 1.500 minutes, max|u|: (1.75e+00, 2.27e+00, 2.17e-02) m s⁻¹, extrema(T): (32.24, -1.88) ᵒC, wall time: 47.479 minutes
[ Info: Simulation is stopping after running for 1.513 hours.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.
[ Info: Iter: 9600, time: 10 days, Δt: 1.500 minutes, max|u|: (2.02e+00, 2.05e+00, 1.19e-02) m s⁻¹, extrema(T): (33.27, -1.88) ᵒC, wall time: 43.331 minutes
Running the simulation for real
After the initial spin up of 10 days, we can increase the time-step and run for longer.
simulation.stop_time = 60days
simulation.Δt = 10minutes
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (1.349 ms)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (8.009 seconds).
[ Info: Iter: 10320, time: 15 days, Δt: 10 minutes, max|u|: (1.93e+00, 1.95e+00, 1.62e-02) m s⁻¹, extrema(T): (34.62, -1.88) ᵒC, wall time: 6.812 minutes
[ Info: Iter: 11040, time: 20 days, Δt: 10 minutes, max|u|: (1.74e+00, 1.83e+00, 1.83e-02) m s⁻¹, extrema(T): (35.96, -1.88) ᵒC, wall time: 6.793 minutes
[ Info: Iter: 11760, time: 25 days, Δt: 10 minutes, max|u|: (1.64e+00, 2.35e+00, 1.77e-02) m s⁻¹, extrema(T): (37.37, -1.88) ᵒC, wall time: 6.944 minutes
[ Info: Iter: 12480, time: 30 days, Δt: 10 minutes, max|u|: (2.11e+00, 2.17e+00, 2.19e-02) m s⁻¹, extrema(T): (38.07, -1.88) ᵒC, wall time: 6.844 minutes
[ Info: Iter: 13200, time: 35 days, Δt: 10 minutes, max|u|: (2.26e+00, 2.29e+00, 3.03e-02) m s⁻¹, extrema(T): (38.22, -1.88) ᵒC, wall time: 6.777 minutes
[ Info: Iter: 13920, time: 40 days, Δt: 10 minutes, max|u|: (1.93e+00, 2.26e+00, 2.80e-02) m s⁻¹, extrema(T): (38.65, -1.88) ᵒC, wall time: 6.812 minutes
[ Info: Iter: 14640, time: 45 days, Δt: 10 minutes, max|u|: (2.03e+00, 2.05e+00, 2.57e-02) m s⁻¹, extrema(T): (38.88, -1.88) ᵒC, wall time: 6.806 minutes
[ Info: Iter: 15360, time: 50 days, Δt: 10 minutes, max|u|: (2.28e+00, 2.38e+00, 2.72e-02) m s⁻¹, extrema(T): (36.86, -1.88) ᵒC, wall time: 6.771 minutes
[ Info: Iter: 16080, time: 55 days, Δt: 10 minutes, max|u|: (2.90e+00, 3.14e+00, 2.64e-02) m s⁻¹, extrema(T): (37.51, -1.88) ᵒC, wall time: 6.757 minutes
[ Info: Simulation is stopping after running for 1.133 hours.
[ Info: Simulation time 60 days equals or exceeds stop time 60 days.
[ Info: Iter: 16800, time: 60 days, Δt: 10 minutes, max|u|: (1.94e+00, 1.64e+00, 2.46e-02) m s⁻¹, extrema(T): (37.62, -1.88) ᵒC, wall time: 6.739 minutes
A pretty movie
It's time to make a pretty movie of the simulation. First we load the output we've been saving on disk and plot the final snapshot:
u = FieldTimeSeries("near_global_surface_fields.jld2", "u"; backend = OnDisk())
v = FieldTimeSeries("near_global_surface_fields.jld2", "v"; backend = OnDisk())
T = FieldTimeSeries("near_global_surface_fields.jld2", "T"; backend = OnDisk())
e = FieldTimeSeries("near_global_surface_fields.jld2", "e"; backend = OnDisk())
times = u.times
Nt = length(times)
n = Observable(Nt)
land = interior(T.grid.immersed_boundary.bottom_height) .>= 0
Tn = @lift begin
Tn = interior(T[$n])
Tn[land] .= NaN
view(Tn, :, :, 1)
end
en = @lift begin
en = interior(e[$n])
en[land] .= NaN
view(en, :, :, 1)
end
un = Field{Face, Center, Nothing}(u.grid)
vn = Field{Center, Face, Nothing}(v.grid)
s = @at (Center, Center, Nothing) sqrt(un^2 + vn^2) # compute √(u²+v²) and interpolate back to Center, Center
s = Field(s)
sn = @lift begin
parent(un) .= parent(u[$n])
parent(vn) .= parent(v[$n])
compute!(s)
sn = interior(s)
sn[land] .= NaN
view(sn, :, :, 1)
end
title = @lift string("Near-global 1/4 degree ocean simulation after ",
prettytime(times[$n] - times[1]))
λ, φ, _ = nodes(T) # T, e, and s all live on the same grid locations
fig = Figure(size = (1000, 1500))
axs = Axis(fig[1, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
axT = Axis(fig[2, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
axe = Axis(fig[3, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
hm = heatmap!(axs, λ, φ, sn, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray)
Colorbar(fig[1, 2], hm, label = "Surface Speed (m s⁻¹)")
hm = heatmap!(axT, λ, φ, Tn, colorrange = (-1, 30), colormap = :magma, nan_color=:lightgray)
Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)")
hm = heatmap!(axe, λ, φ, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray)
Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)")
Label(fig[0, :], title)
save("snapshot.png", fig)
And now we make a movie:
record(fig, "near_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn
n[] = nn
endThis page was generated using Literate.jl.