Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -67,3 +67,4 @@ CondaPkg.toml

# claude
.claude
test/Manifest.toml
2 changes: 1 addition & 1 deletion docs/src/developers/slab_ocean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ set!(sea_ice.model, h=Metadatum(:sea_ice_thickness, dataset=ECCO4Monthly()),
ℵ=Metadatum(:sea_ice_concentration, dataset=ECCO4Monthly()))

interfaces = ComponentInterfaces(atmosphere, slab_ocean, sea_ice; exchange_grid=grid)
coupled_model = NumericalEarth.EarthSystemModel(atmosphere, slab_ocean, sea_ice; interfaces)
coupled_model = NumericalEarth.EarthSystemModel(; atmosphere, sea_ice, ocean=slab_ocean, interfaces)

simulation = Simulation(coupled_model, Δt=60minutes, stop_time=120days)
run!(simulation)
Expand Down
16 changes: 12 additions & 4 deletions docs/src/earth_system_model.md
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ model = OceanOnlyModel(ocean)

# output
EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
Expand Down Expand Up @@ -51,6 +52,7 @@ model

# output
EarthSystemModel{CPU}(time = 1 hour, iteration = 3)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
Expand Down Expand Up @@ -87,10 +89,11 @@ a sea ice component as positional arguments:
```jldoctest esm
ocean = ocean_simulation(grid, timestepper = :QuasiAdamsBashforth2)
sea_ice = FreezingLimitedOceanTemperature()
model = OceanSeaIceModel(ocean, sea_ice)
model = OceanSeaIceModel(sea_ice, ocean)

# output
EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
Expand All @@ -108,9 +111,12 @@ Oceananigans `Simulation`, would also work.
EarthSystemModel
```

The full constructor takes positional arguments `(atmosphere, ocean, sea_ice)` and
gives access to every knob: radiation parameters, reference densities, heat capacities,
and -- most importantly -- the `interfaces` keyword, which controls how fluxes are computed.
The full constructor takes positional arguments
`(radiation, atmosphere, land, sea_ice, ocean)` -- the components in struct
order, top to bottom -- and gives access to every knob: reference densities,
heat capacities, and -- most importantly -- the `interfaces` keyword, which
controls how fluxes are computed. Pass `nothing` for components that are
absent.

## Customizing flux formulations

Expand All @@ -130,6 +136,7 @@ model = OceanOnlyModel(ocean; interfaces)

# output
EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
Expand All @@ -148,6 +155,7 @@ model = OceanOnlyModel(ocean; interfaces)

# output
EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
Expand Down
9 changes: 5 additions & 4 deletions examples/generate_surface_fluxes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,12 @@ S_metadata = ECCOMetadatum(:salinity; date=DateTime(1993, 1, 1))
set!(ocean.model; T=T_metadata, S=S_metadata)

# Finally, we construct a coupled model, which will compute fluxes during construction.
# We omit `sea_ice` so the model is ocean-only, and use the default `Radiation()` that
# uses the two-band shortwave (visible and UV) + longwave (mid and far infrared)
# decomposition of the radiation spectrum.
# We omit `sea_ice` so the model is ocean-only, and pair the JRA55 atmosphere with a
# matching `JRA55PrescribedRadiation` that supplies downwelling shortwave and
# longwave radiation as well as ocean / sea-ice surface properties.

coupled_model = OceanOnlyModel(ocean; atmosphere, radiation=Radiation(eltype))
radiation = JRA55PrescribedRadiation(; backend = JRA55NetCDFBackend(2))
coupled_model = OceanOnlyModel(ocean; atmosphere, radiation)

# Now that the surface fluxes are computed, we can extract and visualize them.
# The turbulent fluxes are stored in `coupled_model.interfaces.atmosphere_ocean_interface.fluxes`.
Expand Down
8 changes: 5 additions & 3 deletions examples/global_climate_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -98,10 +98,12 @@ nothing #hide
nothing #hide

# We build the complete coupled `earth_model` and the coupled simulation.
# Since radiation is idealized in this example, we set the emissivities to zero.
# 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.

radiation = Radiation(ocean_emissivity=0.0, sea_ice_emissivity=0.0)
earth_model = EarthSystemModel(atmosphere, ocean, sea_ice; radiation)
earth_model = EarthSystemModel(; atmosphere, sea_ice, ocean)

# ## Building and running the simulation
#
Expand Down
13 changes: 9 additions & 4 deletions examples/idealized_single_column_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,20 @@ qᵃᵗ = 0.01 # specific humidity
ℐꜜˢʷ = 400 # shortwave radiation (W m⁻², positive means heating right now)

# Build the atmosphere
radiation = Radiation(ocean_albedo=0.1)
atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
atmosphere_times = range(0, 1days, length=3)
atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times)

# Build the radiation component (lives on the same grid + times as the atmosphere
# in this single-column setup) and prescribe a constant shortwave forcing.
# Override the default ocean albedo (0.05) but keep the default emissivity (0.97).
radiation = PrescribedRadiation(atmosphere_grid, atmosphere_times;
ocean_surface = SurfaceRadiationProperties(albedo=0.1))

parent(atmosphere.tracers.T) .= Tᵃᵗ # K
parent(atmosphere.velocities.u) .= u₁₀ # m/s
parent(atmosphere.tracers.q) .= qᵃᵗ # mass ratio
parent(atmosphere.downwelling_radiation.shortwave) .= ℐꜜˢʷ # W
parent(radiation.downwelling_shortwave) .= ℐꜜˢʷ # W

# Build ocean model at rest with initial temperature stratification
grid = RectilinearGrid(size=20, z=(-100, 0), topology=(Flat, Flat, Bounded))
Expand All @@ -40,8 +45,8 @@ Tᵢ(z) = T₀ + dTdz * z
set!(ocean.model, T=Tᵢ, S=S₀)

atmosphere_ocean_fluxes = SimilarityTheoryFluxes(stability_functions=nothing)
interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes)
model = OceanOnlyModel(ocean; atmosphere, interfaces)
interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes, radiation)
model = OceanOnlyModel(ocean; atmosphere, radiation, interfaces)

𝒬ᵛ = model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat
𝒬ᵀ = model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat
Expand Down
8 changes: 4 additions & 4 deletions examples/meridional_heat_transport_ecco.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,10 @@ ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset)
set!(ocean.model, T=ecco_temperature, S=ecco_salinity)
set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration)

radiation = Radiation(arch)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80),
include_rivers_and_icebergs = false)
esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
jra55_backend = JRA55NetCDFBackend(80)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend)
radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend)
esm = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation)

simulation = Simulation(esm; Δt=20minutes, stop_time=5*365days)

Expand Down
20 changes: 5 additions & 15 deletions examples/near_global_ocean_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -88,24 +88,14 @@ set!(ocean.model, T=Metadatum(:temperature, dataset=ECCO4Monthly()),

# ### 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)

# The 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 [`InterfaceComputations`](@ref NumericalEarth.EarthSystemModels.InterfaceComputations).
# The number of snapshots that are loaded into memory is determined by
# the `backend`. Here, we load 41 snapshots at a time into memory.
# Next we build a prescribed atmosphere state and radiation component,
# which together drive the ocean simulation. The atmospheric data and
# downwelling shortwave / longwave radiation are both prescribed using JRA55.
# The number of snapshots loaded into memory is set by the backend.

jra55_backend = JRA55NetCDFBackend(41)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend)
radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend)
land = JRA55PrescribedLand(arch; backend=jra55_backend)

# ## The coupled simulation
Expand Down
8 changes: 6 additions & 2 deletions examples/one_degree_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -95,9 +95,13 @@ set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration)
# ### Atmospheric forcing

# We force the simulation with a JRA55-do atmospheric reanalysis.
radiation = Radiation(arch)
jra55_backend = JRA55NetCDFBackend(80)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend)
# Use a latitude-dependent ocean albedo (Large & Yeager 2009); keep the
# default ocean emissivity (0.97) and sea-ice surface (albedo 0.7,
# emissivity 1.0).
radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend,
ocean_surface = SurfaceRadiationProperties(albedo = LatitudeDependentAlbedo()))
land = JRA55PrescribedLand(arch; backend=jra55_backend)

# ### Coupled simulation
Expand All @@ -107,7 +111,7 @@ land = JRA55PrescribedLand(arch; backend=jra55_backend)

# With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 20 minutes.

coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, land, radiation)
coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, land, radiation)
simulation = Simulation(coupled_model; Δt=20minutes, stop_time=365days)

# ### A progress messenger
Expand Down
10 changes: 7 additions & 3 deletions examples/single_column_os_papa_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,11 @@ atmosphere = JRA55PrescribedAtmosphere(longitude = λ★,
end_date = DateTime(1990, 1, 31), # Last day of the simulation
backend = InMemory())

radiation = JRA55PrescribedRadiation(longitude = λ★,
latitude = φ★,
end_date = DateTime(1990, 1, 31),
backend = InMemory())

# This builds a representation of the atmosphere on the small grid

atmosphere.grid
Expand Down Expand Up @@ -101,7 +106,6 @@ lines!(axq, t_days, qa)
current_figure()

# We continue constructing a simulation.
radiation = Radiation()
coupled_model = OceanOnlyModel(ocean; atmosphere, radiation)
simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days)

Expand Down Expand Up @@ -205,8 +209,8 @@ ua = atmosphere.velocities.u
va = atmosphere.velocities.v
Ta = atmosphere.tracers.T
qa = atmosphere.tracers.q
ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave
ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave
ℐꜜˡʷ = radiation.downwelling_longwave
ℐꜜˢʷ = radiation.downwelling_shortwave
Pr = atmosphere.freshwater_flux.rain
Ps = atmosphere.freshwater_flux.snow

Expand Down
7 changes: 3 additions & 4 deletions examples/veros_ocean_forced_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,11 @@ ocean.set_forcing = set_forcing_tke_only
# radiation, as well as freshwater fluxes.

atmos = JRA55PrescribedAtmosphere(; backend = JRA55NetCDFBackend(10))
radiation = JRA55PrescribedRadiation(; backend = JRA55NetCDFBackend(10))

# The coupled ocean--atmosphere model.
# We use the default radiation model and we do not couple an ice model for simplicity.
# The coupled ocean--atmosphere model. We do not couple an ice model for simplicity.

radiation = Radiation()
coupled_model = OceanSeaIceModel(ocean, nothing; atmosphere=atmos, radiation)
coupled_model = OceanSeaIceModel(nothing, ocean; atmosphere=atmos, radiation)
simulation = Simulation(coupled_model; Δt = 1800, stop_time = 60days)

# We set up a progress callback that will print the current time, iteration, and maximum velocities
Expand Down
4 changes: 2 additions & 2 deletions experiments/arctic_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -90,13 +90,13 @@ set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset),
#####

atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40))
radiation = Radiation()
radiation = JRA55PrescribedRadiation(arch; backend=JRA55NetCDFBackend(40))

#####
##### Arctic coupled model
#####

arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
arctic = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation)
arctic = Simulation(arctic, Δt=5minutes, stop_time=365days)

# Sea-ice variables
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -79,14 +79,17 @@ salinity = ECCOMetadata(:salinity; dir="./")
ice_thickness = ECCOMetadata(:sea_ice_thickness; dir="./")
ice_concentration = ECCOMetadata(:sea_ice_concentration; dir="./")

atmosphere = JRA55PrescribedAtmosphere(arch, backend=JRA55NetCDFBackend(20))
radiation = Radiation(ocean_albedo = LatitudeDependentAlbedo(), sea_ice_albedo=0.6)
jra55_backend = JRA55NetCDFBackend(20)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend)
radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend,
ocean_surface = SurfaceRadiationProperties(LatitudeDependentAlbedo(), 0.97),
sea_ice_surface = SurfaceRadiationProperties(0.6, 1.0))
Comment thread
glwagner marked this conversation as resolved.

set!(ocean.model, T=temperature, S=salinity)
set!(sea_ice.model.ice_thickness, ice_thickness, inpainting=NearestNeighborInpainting(1))
set!(sea_ice.model.ice_concentration, ice_concentration, inpainting=NearestNeighborInpainting(1))

earth_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
earth_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation)
earth = Simulation(earth_model; Δt=30minutes, stop_iteration=10, stop_time=30days)

u, v, _ = ocean.model.velocities
Expand Down
2 changes: 1 addition & 1 deletion experiments/flux_climatology/flux_climatology.jl
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(1000))
##### A prescribed earth...
#####

earth_model = OceanOnlyModel(ocean; atmosphere, radiation = Radiation(ocean_emissivity=0, ocean_albedo=1))
earth_model = OceanOnlyModel(ocean; atmosphere) # radiation off (radiation=nothing default)
earth = Simulation(earth_model, Δt=3hours, stop_time=365days)

wall_time = Ref(time_ns())
Expand Down
5 changes: 3 additions & 2 deletions experiments/one_degree_simulation/one_degree_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surfac
set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
S=ECCOMetadata(:salinity; dates=first(dates)))

radiation = Radiation(arch)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
jra55_backend = JRA55NetCDFBackend(41)
atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend)
radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend)
coupled_model = OceanOnlyModel(ocean; atmosphere, radiation)

simulation = Simulation(coupled_model; Δt=10minutes, stop_iteration=100)
Expand Down
6 changes: 3 additions & 3 deletions ext/NumericalEarthReactantExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,9 @@ const OceananigansReactantExt = Base.get_extension(
Oceananigans, :OceananigansReactantExt
)

const ReactantOSIM{I, A, L, O, F, C} = Union{
EarthSystemModel{I, A, L, O, F, C, <:ReactantState},
EarthSystemModel{I, A, L, O, F, C, <:Distributed{ReactantState}},
const ReactantOSIM{R, A, L, I, O, F, C} = Union{
EarthSystemModel{R, A, L, I, O, F, C, <:ReactantState},
EarthSystemModel{R, A, L, I, O, F, C, <:Distributed{ReactantState}},
}

reconcile_state!(model::ReactantOSIM) = nothing
Expand Down
24 changes: 6 additions & 18 deletions src/Atmospheres/interpolate_atmospheric_state.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,6 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c
atmosphere_tracers = (T = atmosphere.tracers.T.data,
q = atmosphere.tracers.q.data)

ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave
ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave
downwelling_radiation = (shortwave=ℐꜜˢʷ.data, longwave=ℐꜜˡʷ.data)
freshwater_flux = map(ϕ -> ϕ.data, atmosphere.freshwater_flux)
atmosphere_pressure = atmosphere.pressure.data

Expand All @@ -44,14 +41,12 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c

# Simplify NamedTuple to reduce parameter space consumption.
# See https://github.com/CliMA/NumericalEarth.jl/issues/116.
atmosphere_data = (u = atmosphere_fields.u.data,
v = atmosphere_fields.v.data,
T = atmosphere_fields.T.data,
p = atmosphere_fields.p.data,
q = atmosphere_fields.q.data,
ℐꜜˢʷ = atmosphere_fields.ℐꜜˢʷ.data,
ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data,
Jᶜ = atmosphere_fields.Jᶜ.data)
atmosphere_data = (u = atmosphere_fields.u.data,
v = atmosphere_fields.v.data,
T = atmosphere_fields.T.data,
p = atmosphere_fields.p.data,
q = atmosphere_fields.q.data,
Jᶜ = atmosphere_fields.Jᶜ.data)

kernel_parameters = interface_kernel_parameters(grid)

Expand All @@ -72,7 +67,6 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c
atmosphere_velocities,
atmosphere_tracers,
atmosphere_pressure,
downwelling_radiation,
freshwater_flux,
atmosphere_backend,
atmosphere_time_indexing)
Expand All @@ -99,7 +93,6 @@ end
atmos_velocities,
atmos_tracers,
atmos_pressure,
downwelling_radiation,
prescribed_freshwater_flux,
atmos_backend,
atmos_time_indexing)
Expand All @@ -121,9 +114,6 @@ end
qᵃᵗ = interp_atmos_time_series(atmos_tracers.q, atmos_args...)
pᵃᵗ = interp_atmos_time_series(atmos_pressure, atmos_args...)

ℐꜜˢʷ = interp_atmos_time_series(downwelling_radiation.shortwave, atmos_args...)
ℐꜜˡʷ = interp_atmos_time_series(downwelling_radiation.longwave, atmos_args...)

# Usually precipitation
Mh = interp_atmos_time_series(prescribed_freshwater_flux, atmos_args...)

Expand All @@ -138,8 +128,6 @@ end
surface_atmos_state.T[i, j, 1] = Tᵃᵗ
surface_atmos_state.p[i, j, 1] = pᵃᵗ
surface_atmos_state.q[i, j, 1] = qᵃᵗ
surface_atmos_state.ℐꜜˢʷ[i, j, 1] = ℐꜜˢʷ
surface_atmos_state.ℐꜜˡʷ[i, j, 1] = ℐꜜˡʷ
surface_atmos_state.Jᶜ[i, j, 1] = Mh
end
end
Expand Down
Loading
Loading