Orography
Orography (in height above the surface) forms the surface boundary of the lowermost layer in SpeedyWeather.
In the shallow-water equations the orography $H_b$ enters the equations when computing the layer thickness $h = \eta + H_0 - H_b$ for the volume fluxes $\mathbf{u}h$ in the continuity equation. Here, the orography is used in meters above the surface which shortens $h$ over mountains. The orography here is needed in grid-point space.
In the primitive equations the orography enters the equations when computing the Geopotential. So actually required here is the surface geopotential $\Phi_s = gz_s$ where $z_s$ is the orography height in meters as used in the shallow-water equations too $z_s = H_b$. However, the primitive equations require the orography in spectral space as the geopotential calculation is a linear operation in the horizontal and can therefore be applied in either grid-point or spectral space. The latter is more convenient as SpeedyWeather solves the equations to avoid additional transforms.
In the current formulation of the barotropic vorticity equations there is no orography. In fact, the field model.orography is not defined for model::BarotropicModel.
Orographies implemented
Currently implemented are
using SpeedyWeather
subtypes(SpeedyWeather.AbstractOrography)3-element Vector{Any}:
EarthOrography
NoOrography
ZonalRidgewhich are
- $\Phi_s = z_s = H_b = 0$ for
NoOrography - For
ZonalRidgethe zonal ridge from the Jablonowski and Williamson initial conditions, see Jablonowski-Williamson baroclinic wave - For
EarthOrographya high-resolution orography is loaded and interpolated to the resolution as defined byspectral_grid.
all orographies need to be created with spectral_grid::SpectralGrid as the first argument, so that the respective fields for geopot_surf, i.e. $\Phi_s$ and orography, i.e. $H_b$ can be allocated in the right size and number format.
Earth's orography
Earth's orography can be created with (here we use a resolution of T85, about 165km globally)
using SpeedyWeather
spectral_grid = SpectralGrid(trunc=85)
orography = EarthOrography(spectral_grid)EarthOrography{Float32, OctahedralGaussianGrid{Float32}} <: SpeedyWeather.AbstractOrography
├ path::String = SpeedyWeather.jl/input_data
├ file::String = orography.nc
├ file_Grid::UnionAll = FullGaussianGrid
├ scale::Float64 = 1.0
├ smoothing::Bool = true
├ smoothing_power::Float64 = 1.0
├ smoothing_strength::Float64 = 0.1
├ smoothing_fraction::Float64 = 0.05
└── arrays: orography, geopot_surfbut note that only allocates the orography, it does not actually load and interpolate the orography which happens at the initialize! step. Visualised with
model = PrimitiveDryModel(;spectral_grid, orography)
initialize!(orography, model) # happens also in simulation = initialize!(model)
using CairoMakie
heatmap(orography.orography, title="Earth's orography at T85 resolution, no smoothing")
typing ?EarthOrography shows the various options that are provided. An orogaphy at T85 resolution that is as smooth as it would be at T42 (controlled by the smoothing_fraction, the fraction of highest wavenumbers which are the top half here, about T43 to T85) for example can be created with
orography = EarthOrography(spectral_grid, smoothing=true, smoothing_fraction=0.5)
initialize!(orography, model)
heatmap(orography.orography, title="Earth's orography at T85 resolution, smoothed to T42")
Load orography from file
The easiest to load another orography from a netCDF file is to reuse the EarthOrography, e.g.
mars_orography = EarthOrography(spectal_grid,
path="path/to/my/orography",
file="mars_orography.nc",
file_Grid=FullClenshawGrid)the orography itself need to come on one of the full grids SpeedyWeather defines, i.e. FullGaussianGrid or FullClenshawGrid (a regular lat-lon grid, see FullClenshawGrid), which you can specify. Best to inspect the correct orientation with plot(mars_orography.orography) (where the first is whatever name you chose here). You can use smoothing as above.
Changing orography manually
You can also change orography manually, that means by mutating the elements in either orography.orography (to set it for the shallow-water model) or orography.geopot_surf (for the primitive equations, but this is in spectral space, advanced!). After the orography has been initialised (for testing to initialize!(orography, model) but in most cases when you do simulation = initialize!(model)), for example you can do
sort!(orography.orography)to move all mountains to the south pole, or
orography.orography[1] = 100100to set the grid point 1 (at 0˚E, on the first ring around the north pole) to a height of 100m. Whichever way you tweak the orography manually you want to reflect this in the surface geopotential geopot_surf which is used in the primitive equations by
spectral!(orography.geopot_surf, orography.orography, model.spectral_transform)
orography.geopot_surf .*= model.planet.gravity
spectral_truncation!(orography.geopot_surf)87×86 LowerTriangularMatrix{ComplexF32}:
8118.62+0.0im 0.0+0.0im … 0.0+0.0im
-11522.2+0.0im 0.113159-63.0705im 0.0+0.0im
9690.33+0.0im 0.0933886+91.7965im 0.0+0.0im
-6662.48+0.0im -0.342513-88.9813im 0.0+0.0im
4177.93+0.0im 1.11322+71.2606im 0.0+0.0im
-3201.76+0.0im -1.34881-65.8504im … 0.0+0.0im
2599.81+0.0im 1.71177+62.2399im 0.0+0.0im
-2118.58+0.0im -1.35034-58.1854im 0.0+0.0im
1481.34+0.0im 1.27353+45.976im 0.0+0.0im
-1144.17+0.0im -0.404195-40.3456im 0.0+0.0im
⋮ ⋱ ⋮
11.74+0.0im 0.186785+3.88613im 0.0+0.0im
-12.93+0.0im -0.0695883-5.16978im 0.0+0.0im
9.48157+0.0im 0.229045+2.73734im … 0.0+0.0im
-11.3232+0.0im 0.491396-4.96867im 0.0+0.0im
9.37732+0.0im -0.718931+3.21253im 0.0+0.0im
-12.9287+0.0im 1.86125-5.65274im 0.0+0.0im
11.8514+0.0im -2.70572+4.21834im 0.0+0.0im
-13.7203+0.0im 4.03883-6.85891im … 0.0599032-0.0341593im
0.0+0.0im 0.0+0.0im 0.0+0.0imIn the first line, the surface geopotential is still missing the gravity, which is multiplied in the second line. The spectral_truncation! only removes the $l_{max}+1$ degree of the spherical harmonics that scalar fields do not use, see One more degree for spectral fields.
Defining a new orography type
You can also define a new orography like we defined ZonalRidge or EarthOrography. The following explains what's necessary for this. The new MyOrography has to be defined as (mutable or not, but always with @kwdef)
@kwdef struct MyOrography{NF, Grid<:RingGrids.AbstractGrid{NF}} <: SpeedyWeather.AbstractOrography
# optional, any parameters as fields here, e.g.
constant_height::Float64 = 100
# add some other parameters with default values
# mandatory, every <:AbstractOrography needs those (same name, same type)
orography::Grid # in grid-point space [m]
geopot_surf::LowerTriangularMatrix{Complex{NF}} # in spectral space *gravity [m^2/s^2]
endfor convenience with a generator function is automatically defined for all AbstractOrography
my_orography = MyOrography(spectral_grid, constant_height=200)Main.MyOrography{Float32, OctahedralGaussianGrid{Float32}} <: SpeedyWeather.AbstractOrography
├ constant_height::Float64 = 200.0
└── arrays: orography, geopot_surfNow we have to extend the initialize! function. The first argument has to be ::MyOrography i.e. the new type we just defined, the second argument has to be ::ModelSetup although you could constrain it to ::ShallowWater for example but then it cannot be used for primitive equations.
function SpeedyWeather.initialize!(
orog::MyOrography, # first argument as to be ::MyOrography, i.e. your new type
model::ModelSetup, # second argument, use anything from model read-only
)
(; orography, geopot_surf) = orog # unpack
# maybe use lat, lon coordinates (in degree or radians)
(; latds, londs, lats, lons) = model.geometry
# change here the orography grid [m], e.g.
orography .= orography.constant_height
# then also calculate the surface geopotential for primitive equations
# given orography we just set
spectral!(geopot_surf, orography, model.spectral_transform)
geopot_surf .*= model.planet.gravity
spectral_truncation!(geopot_surf)
return nothing
end