Skip to content

Commit

Permalink
Held-Suarez example, NoTypes for surface fluxes (#514)
Browse files Browse the repository at this point in the history
* Split into 2D/3D examples

* NoTypes for surface fluxes

* Held-Suarez example added

* SpeedyTransforms docs with Makie too
  • Loading branch information
milankl authored Apr 5, 2024
1 parent c90d5b8 commit 35314c0
Show file tree
Hide file tree
Showing 16 changed files with 405 additions and 190 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ about dos and don'ts. Just express your interest to contribute and we'll be happ
## Example use

For a more comprehensive tutorial with several examples, see
[Examples](https://speedyweather.github.io/SpeedyWeather.jl/dev/examples/) in the documentation.
[Examples](https://speedyweather.github.io/SpeedyWeather.jl/dev/examples_2D/) in the documentation.
The interface to SpeedyWeather.jl consist of 5 steps: define the grid, create model components,
construct the model, initialize, run

Expand Down
3 changes: 2 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ makedocs(
],
"Running SpeedyWeather" => [
"How to run SpeedyWeather"=>"how_to_run_speedy.md",
"Examples"=>"examples.md",
"Examples 2D"=>"examples_2D.md",
"Examples 3D"=>"examples_3D.md",
"Analysis"=>"analysis.md",
"Tree structure"=>"structure.md",
"Particle advection"=>"particles.md",
Expand Down
29 changes: 25 additions & 4 deletions docs/src/convection.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,17 @@ precipitating as the condensed humidity forms cloud droplets that
eventually fall down as convective precipitation.
See also [Large-scale condensation](@ref) in comparison.

## Simplified Betts-Miller convective adjustment
## Convection implementations

Currently implemented are
```@example convection
using InteractiveUtils # hide
using SpeedyWeather
subtypes(SpeedyWeather.AbstractConvection)
```
which are described in the following.

## [Simplified Betts-Miller convection](@id BettsMiller)

We follow the simplification of the Betts-Miller convection scheme
[^Betts1986][^BettsMiller1986] as studied by Frierson, 2007 [^Frierson2007].
Expand All @@ -24,9 +34,9 @@ convection as an adjustment towards a (pseudo-) moist adiabat
reference profile and its associated humidity profile. Meaning that
conceptually for every vertical column in the atmosphere we

1. Diagnose the vertical temperature and humidity profile of the environment relative to the surface and calculate the adiabat to the level of zero buoyancy.
1. Diagnose the vertical temperature and humidity profile of the environment relative to the adiabat up to the level of zero buoyancy.
2. Decide whether convection should take place and whether it is deep (precipitating) or shallow (non-precipitating).
3. Relax temperature and humidity towards (adjusted) profiles from 1.
3. Relax temperature and humidity towards (corrected) profiles from 1.

## Reference profiles

Expand Down Expand Up @@ -152,7 +162,8 @@ Q_{ref} &= \int_{p_0}^{p_{LZB}} -q_{ref} dp \\
f_q &= 1 - \frac{\Delta q}{Q_ref} \\
q_{ref, 2} &= f_q q_{ref} \\
\Delta T &= \frac{1}{\Delta p} \int_{p_0}^{p_{LZB}} -(T - T_{ref}) dp \\
T_{ref,2} = T_{ref} - \Delta T
T_{ref,2} &= T_{ref} - \Delta T
\end{aligned}
```

## Corrected relaxation
Expand Down Expand Up @@ -182,6 +193,16 @@ P = -\int \frac{\Delta t}{g \rho} \delta q dp
In the shallow convection case ``P=0`` due to the correction even though in
the first guess relaxation ``P<0`` was possible, but for deep convection ``P>0`` by definition.

## Dry convection

In the primitive equation model with humidity the [Betts-Miller convection scheme](@ref BettsMiller)
as described above is defined. Without humidity, a dry version reduces to the
[Shallow convection](@ref) case. The two different shallow convection schemes in
Frierson 2007[^Frierson2007], the "shallower" shallow convection scheme and the "qref"
(as implemented here in [Shallow convection](@ref)) in that case also reduce to
the same formulation. The dry Betts-Miller convection scheme is the default
in the primitive equation model without humidity.

## References

[^Betts1986]: Betts, A. K., 1986: A new convective adjustment scheme. Part I: Observational and theoretical basis. Quart. J. Roy. Meteor. Soc.,112, 677-691. DOI: [10.1002/qj.49711247307](https://doi.org/10.1002/qj.49711247307)
Expand Down
124 changes: 18 additions & 106 deletions docs/src/examples.md → docs/src/examples_2D.md
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
# Examples
# [Examples 2D](@id Examples)

The following is a collection of model setups, starting with an easy setup
The following is a collection of example model setups, starting with an easy setup
of the [Barotropic vorticity equation](@ref barotropic_vorticity_model) and
continuing with more complicated setups.
continuing with the [shallow water model](@ref shallow_water_model).

See also [Examples 3D](@ref) for examples with the primitive equation models.

## 2D turbulence on a non-rotating sphere

Expand Down Expand Up @@ -137,9 +139,12 @@ nothing # hide
You see that in comparison the unicode plot heavily coarse-grains the simulation, well it's unicode after all!
Here, we have unpacked the netCDF file using [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl)
and then plotted via `heatmap(lon, lat, vor)`. While you can do that to give you more control
on the plotting, SpeedyWeather.jl also defines an extension for Makie.jl, see [Extensions](@ref)
because if our matrix `vor` here was an `AbstractGrid` (see [RingGrids](@ref)) then all
its geographic information (which grid point is where) is encoded in the type. So we can also do
on the plotting, SpeedyWeather.jl also defines an extension for Makie.jl, see [Extensions](@ref).
Because if our matrix `vor` here was an `AbstractGrid` (see [RingGrids](@ref)) then all
its geographic information (which grid point is where) would directly be encoded in the type.
From the netCDF file you need to use the longitude and latitude dimensions.

So we can also just do
```@example galewsky_setup
vor_grid = FullGaussianGrid(vor)
Expand All @@ -152,15 +157,15 @@ nothing # hide

Note that here you need to know which grid the data comes on (an error is thrown if `FullGaussianGrid(vor)`
is not size compatible). By default the output will be on the FullGaussianGrid, but if you
play around with other grids, you'd need to change this here.
play around with other grids, you'd need to change this here,
see [NetCDF output on other grids](@ref output_grid).

We did want to showcase the usage of [NetCDF output](@ref) here, but from now on
we will use `heatmap` to plot data on our grids directly, without storing output first.
So for our current simulation, that means at time = 12 days, vorticity on the grid
is stored in the diagnostic variables and can be visualised with

```@example galewsky_setup
t = ds.dim["time"]
vor = simulation.diagnostic_variables.layers[1].grid_variables.vor_grid
heatmap(vor, title="Relative vorticity [1/s]")
save("galewsky2.png", ans) # hide
Expand All @@ -176,9 +181,9 @@ Let's try it out! We create an `EarthOrography` struct like so
orography = EarthOrography(spectral_grid)
```

It will read the orography from file as shown, and there are some smoothing options too, but let's not change them.
Same as before, create a model, initialize into a simulation, run. This time directly for 12 days so that we can
compare with the last plot
It will read the orography from file as shown (only at `initialize!(model)`), and there are some smoothing
options too, but let's not change them. Same as before, create a model, initialize into a simulation, run.
This time directly for 12 days so that we can compare with the last plot

```@example galewsky_setup
model = ShallowWaterModel(; spectral_grid, orography, initial_conditions)
Expand Down Expand Up @@ -304,107 +309,14 @@ Hb = model.orography.orography
η = simulation.diagnostic_variables.surface.pres_grid
h = @. η + H - Hb # @. to broadcast grid + scalar - grid
heatmap(h, title="Dynamic layer thickness h")
heatmap(h, title="Dynamic layer thickness h", colormap=:oslo)
save("gravity_waves.png", ans) # hide
nothing # hide
```
![Gravity waves pyplot](gravity_waves.png)

Can you spot the Himalayas or the Andes?

## Jablonowski-Williamson baroclinic wave

```@example jablonowski
using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, nlev=8, Grid=FullGaussianGrid, dealiasing=3)
orography = ZonalRidge(spectral_grid)
initial_conditions = InitialConditions(
vordiv = ZonalWind(),
temp = JablonowskiTemperature(),
pres = ZeroInitially())
model = PrimitiveDryModel(; spectral_grid, orography, initial_conditions, physics=false)
simulation = initialize!(model)
model.feedback.verbose = false # hide
run!(simulation, period=Day(9))
nothing # hide
```

The Jablonowski-Williamson baroclinic wave test case[^JW06] using the
[Primitive equation model](@ref primitive_equation_model) particularly the dry model,
as we switch off all physics with `physics=false`.
We want to use 8 vertical levels, and a lower resolution of T31 on a
[full Gaussian grid](@ref FullGaussianGrid).
The Jablonowski-Williamson initial conditions are `ZonalWind` for vorticity and divergence
(curl and divergence of ``u, v``), `JablonowskiTemperature` for temperature and
`ZeroInitially` for pressure. The orography is just a `ZonalRidge`.
There is no forcing and the initial conditions are baroclinically unstable which kicks
off a wave propagating eastward. This wave becomes obvious when visualised with

```@example jablonowski
using CairoMakie
vor = simulation.diagnostic_variables.layers[end].grid_variables.vor_grid
heatmap(vor, title="Surface relative vorticity")
save("jablonowski.png", ans) # hide
nothing # hide
```
![Jablonowski pyplot](jablonowski.png)

## Aquaplanet

```@example aquaplanet
using SpeedyWeather
# components
spectral_grid = SpectralGrid(trunc=31, nlev=5)
ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273)
land_sea_mask = AquaPlanetMask(spectral_grid)
orography = NoOrography(spectral_grid)
# create model, initialize, run
model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography)
simulation = initialize!(model)
model.feedback.verbose = false # hide
run!(simulation, period=Day(50))
nothing # hide
```

Here we have defined an aquaplanet simulation by
- creating an `ocean::AquaPlanet`. This will use constant sea surface temperatures that only vary with latitude.
- creating a `land_sea_mask::AquaPlanetMask` this will use a land-sea mask with `false`=ocean everywhere.
- creating an `orography::NoOrography` which will have no orography and zero surface geopotential.

All passed on to the model constructor for a `PrimitiveWetModel`, we have now a model with humidity
and physics parameterization as they are defined by default (typing `model` will give you an overview
of its components). We could have change the `model.land` and `model.vegetation` components too,
but given the land-sea masks masks those contributions to the surface fluxes anyway, this is not
necessary. Note that neither sea surface temperature, land-sea mask
or orography have to agree. It is possible to have an ocean on top of a mountain.
For an ocean grid-cell that is (partially) masked by the land-sea mask, its value will
be (fractionally) ignored in the calculation of surface fluxes (potentially leading
to a zero flux depending on land surface temperatures).

Now with the following we visualize the surface humidity after the 50 days of
simulation. We use 50 days as without mountains it takes longer for the initial conditions to
become unstable. The surface humidity shows small-scale patches in the tropics, which is a result
of the convection scheme, causing updrafts and downdrafts in both humidity and temperature.

```@example aquaplanet
using CairoMakie
humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid
heatmap(humid, title="Surface specific humidity [kg/kg]")
save("aquaplanet.png", ans) # hide
nothing # hide
```
![Aquaplanet pyplot](aquaplanet.png)



## References

[^G04]: Galewsky, Scott, Polvani, 2004. *An initial-value problem for testing numerical models of the global shallow-water equations*, Tellus A. DOI: [10.3402/tellusa.v56i5.14436](https://doi.org/10.3402/tellusa.v56i5.14436)
[^JW06]: Jablonowski, C. and Williamson, D.L. (2006), A baroclinic instability test case for atmospheric model dynamical cores. Q.J.R. Meteorol. Soc., 132: 2943-2975. [10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12)
[^G04]: Galewsky, Scott, Polvani, 2004. *An initial-value problem for testing numerical models of the global shallow-water equations*, Tellus A. DOI: [10.3402/tellusa.v56i5.14436](https://doi.org/10.3402/tellusa.v56i5.14436)
Loading

0 comments on commit 35314c0

Please sign in to comment.