diff --git a/docs/make.jl b/docs/make.jl
index 68e47a329..ed33d7c27 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -9,6 +9,7 @@ makedocs(
     pages=["Home"=>"index.md",
             "Installation"=>"installation.md",
             "How to run SpeedyWeather.jl"=>"how_to_run_speedy.md",
+            "Model setups"=>"setups.md",
             "Spherical Harmonic Transform"=>"spectral_transform.md",
             "Grids"=>"grids.md",
             "Barotropic model"=>"barotropic.md",
diff --git a/docs/src/barotropic.md b/docs/src/barotropic.md
index 5294d6d96..e5b72459a 100644
--- a/docs/src/barotropic.md
+++ b/docs/src/barotropic.md
@@ -89,7 +89,7 @@ Now loop over
 
 ## [Horizontal diffusion](@id diffusion)
 
-In SpeedyWeather.jl we use hyerdiffusion through an ``n``-th power Laplacian ``(-1)^{n+1}\nabla^{2n}`` (hyper when ``n>1``) which
+In SpeedyWeather.jl we use hyperdiffusion through an ``n``-th power Laplacian ``(-1)^{n+1}\nabla^{2n}`` (hyper when ``n>1``) which
 can be implemented as a multiplication of the spectral coefficients ``\Psi_{lm}`` with ``(-l(l+1))^nR^{-2n}`` (see
 spectral [Laplacian](@ref)) It is therefore computationally not more expensive to apply hyperdiffusion over diffusion
 as the ``(-l(l+1))^nR^{-2n}`` can be precomputed. Note the sign change ``(-1)^{n+1}`` here is such that the dissipative nature of the diffusion operator is retained for ``n`` odd and even.
@@ -162,7 +162,7 @@ A scaling for vorticity ``\zeta`` and stream function ``\Psi`` is used that is
 \tilde{\zeta} = \zeta R, \tilde{\Psi} = \Psi R^{-1}.
 ```
 This is also convenient as vorticity is often ``10^{-5}\text{ s}^{-1}`` in the atmosphere,
-but the streamfunction more like ``10^5\text{ m}^2\text{ s}^{-1}`` and so this scaling
+but the stream function more like ``10^5\text{ m}^2\text{ s}^{-1}`` and so this scaling
 brings both closer to 1 with a typical radius of the Earth of 6371km.
 The inversion of the Laplacians in order to obtain ``\Psi`` from ``\zeta`` therefore becomes
 ```math
@@ -184,7 +184,7 @@ with
 So scaling with the radius squared means we can use dimensionless operators, however, this comes at the
 cost of needing to deal with both a time step in seconds as well as a scaled time step in seconds per
 meter, which can be confusing. Furthermore, some constants like Coriolis or the diffusion coefficient
-need to be scaled too during initialisation, which may be confusing too because values are not
+need to be scaled too during initialization, which may be confusing too because values are not
 what users expect them to be. SpeedyWeather.jl follows the logic that the scaling to the prognostic
 variables is only applied just before the time integration and variables are unscaled for output
 and after the time integration finished. That way, the scaling is hidden as much as possible from
@@ -211,7 +211,7 @@ For the Leapfrog time integration two time steps of the prognostic variables hav
 to ``i-1`` in a step that also swaps the indices for the next time step ``i \to i-1`` and ``i+1 \to i``,
 so that no additional memory than two time steps have to be stored at the same time.
 
-The Leapfrog time integration has to be initialised with an Euler forward step in order
+The Leapfrog time integration has to be initialized with an Euler forward step in order
 to have a second time step ``i+1`` available when starting from ``i`` to actually leapfrog over.
 SpeedyWeather.jl therefore does two initial time steps that are different from
 the leapfrog time steps that follow and that have been described above.
diff --git a/docs/src/extending.md b/docs/src/extending.md
index b6242c83c..eb2bb207b 100644
--- a/docs/src/extending.md
+++ b/docs/src/extending.md
@@ -1,3 +1,3 @@
-# New model setups
+# Extending SpeedyWeather
 
 more to come...
\ No newline at end of file
diff --git a/docs/src/grids.md b/docs/src/grids.md
index fe4b9b623..7e3e197be 100644
--- a/docs/src/grids.md
+++ b/docs/src/grids.md
@@ -3,19 +3,15 @@
 The spectral transform (the [Spherical Harmonic Transform](@ref)) in SpeedyWeather.jl supports any ring-based
 equi-longitude grid. Several grids are already implemented but other can be added. The following pages will
 describe an overview of these grids and but let's start but how they can be used
-```julia
-julia> spectral_grid = SpectralGrid(Grid = FullGaussianGrid)
-SpectralGrid:
- Spectral:   T31 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
- Grid:       4608-element, 48-ring FullGaussianGrid{Float32} (quadratic)
- Resolution: 333km (average)
- Vertical:   8-level SigmaCoordinates
+```@example grids
+using SpeedyWeather
+spectral_grid = SpectralGrid(Grid = FullGaussianGrid)
 ```
 The life of every SpeedyWeather.jl simulation starts with a `SpectralGrid` object which defines the
 resolution in spectral and in grid-point space. The generator `SpectralGrid()` can take as a keyword
 argument `Grid` which can be any of the grids described below. The resolution of the grid, however,
 is not directly chosen, but determined from the spectral resolution `trunc` and the `dealiasing`
-factor. More in [Matching spectral and grid resolution](@ref).
+factor. More in [SpectralGrid](@ref) and [Matching spectral and grid resolution](@ref).
 
 !!! info "RingGrids is a module too!"
     While RingGrids is the underlying module that SpeedyWeather.jl uses for data structs
@@ -179,7 +175,7 @@ Exactness here means that only rounding errors occur in the transform, meaning t
 transform error is very small compared to other errors in a simulation.
 This property arises from that property of the
 [Gauss-Legendre quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature),
-which is used in the [Spherical harmonic transform](@ref) with a full Gaussian grid.
+which is used in the [Spherical Harmonic Transform](@ref) with a full Gaussian grid.
 
 On the full Gaussian grid there are in total ``N_\phi N_\theta`` grid points, which are squeezed
 towards the poles, making the grid area smaller and smaller following a cosine.
diff --git a/docs/src/how_to_run_speedy.md b/docs/src/how_to_run_speedy.md
index c9729ee06..67ac98564 100644
--- a/docs/src/how_to_run_speedy.md
+++ b/docs/src/how_to_run_speedy.md
@@ -1,245 +1,212 @@
 # How to run SpeedyWeather.jl
 
-Creating a SpeedyWeather.jl simulation and running it consists conceptually of 4 steps
+Creating a SpeedyWeather.jl simulation and running it consists conceptually of 4 steps.
+In contrast to many other models, these steps are bottom-up rather then top-down.
+There is no monolithic interface to SpeedyWeather.jl, instead all options that a
+user may want to adjust are chosen and live in their respective model components.
 
-1. Create a `SpectralGrid` which defines the grid and spectral resolution
-2. Create a model
-3. Initialize a model to obtain a `Simulation`.
-4. Run the simulation.
+1. Create a [SpectralGrid](@ref) which defines the grid and spectral resolution.
+2. [Create model components](@ref create_model_components) and combine to a [model](@ref create_model).
+3. [Initialize the model to create a simulation](@ref initialize).
+4. [Run that simulation](@ref run).
 
-In the following we will describe these steps in more detail,
-but let's start with some examples first.
-
-## Example 1: 2D turbulence on a non-rotating sphere
-
-We want to use the barotropic model to simulate some free-decaying 2D turbulence
-on the sphere without rotation. We start by defining the `SpectralGrid` object.
-To have a resolution of about 200km, we choose a spectral resolution of
-T63 (see [Available horizontal resolutions](@ref)) and `nlev=1` vertical levels.
-The `SpectralGrid` object will provide us with some more information
-```@example howtorun
-using SpeedyWeather
-spectral_grid = SpectralGrid(trunc=63,nlev=1)
-```
-
-We could have specified further options, but let's ignore that for now.
-Next step we create a planet that's like Earth but not rotating
-```@example howtorun
-still_earth = Earth(rotation=0)
-```
-There are other options to create a planet but they are irrelevant for the
-barotropic vorticity equations. We also want to specify the initial conditions,
-randomly distributed vorticity is already defined
-```@example howtorun
-using Random # hide
-Random.seed!(1234) # hide
-initial_conditions = StartWithRandomVorticity()
-```
-By default, the power of vorticity is spectrally distributed with ``k^{-3}``, ``k`` being the
-horizontal wavenumber, and the amplitude is ``10^{-5}\text{s}^{-1}``.
-
-Now we want to construct a `BarotropicModel`
-with these
-```@example howtorun
-model = BarotropicModel(;spectral_grid, initial_conditions, planet=still_earth)
-nothing # hide
-```
-The `model` contains all the parameters, but isn't initialized yet, which we can do
-with and then run it. The `run!` command will always return the prognostic variables, which, by default, are 
-plotted for surface relative vorticity with a unicode plot. The resolution of the plot
-is not necessarily representative but it lets us have a quick look at the result
-```@example howtorun
-simulation = initialize!(model)
-run!(simulation,n_days=20)
-```
-
-Woohoo! I can see turbulence! You could pick up where this simulation stopped by simply
-doing `run!(simulation,n_days=50)` again. We didn't store any output, which
-you can do by `run!(simulation,output=true)`, which will switch on NetCDF output
-with default settings. More options on output in [NetCDF output](@ref).
-
-## Example 2: Shallow water with mountains
-
-As a second example, let's investigate the Galewsky et al.[^1] test case for the shallow
-water equations with and without mountains. As the shallow water system has also only
-one level, we can reuse the `SpectralGrid` from Example 1.
-```@example howtorun
-spectral_grid = SpectralGrid(trunc=63,nlev=1)
-```
-Now as a first simulation, we want to disable any orography, so we create a `NoOrography`
-```@example howtorun
-orography = NoOrography(spectral_grid)
-```
-Although the orography is zero, you have to pass on `spectral_grid` so that it can
-still initialize zero-arrays of the right size and element type. Awesome.
-This time the initial conditions should be set the the Galewsky et al.[^1] zonal
-jet, which is already defined as
-```@example howtorun
-initial_conditions = ZonalJet()
-```
-The jet sits at 45˚N with a maximum velocity of 80m/s and a perturbation as described in their paper.
-Now we construct a model, but this time a `ShallowWaterModel`
-```@example howtorun
-model = ShallowWaterModel(;spectral_grid, orography, initial_conditions)
-simulation = initialize!(model)
-run!(simulation,n_days=6)
-```
-Oh yeah. That looks like the wobbly jet in their paper. Let's run it again for another 6 days
-but this time also store [NetCDF output](@ref).
-```@example howtorun
-run!(simulation,n_days=6,output=true)
-```
-The progress bar tells us that the simulation run got the identification "0001"
-(which just counts up, so yours might be higher), meaning that
-data is stored in the folder `/run_0001`, so let's plot that data properly (and not just using UnicodePlots).
-```@example howtorun
-using PythonPlot, NCDatasets
-ioff() # hide
-ds = NCDataset("run_0001/output.nc")
-ds["vor"]
-```
-Vorticity `vor` is stored as a lon x lat x vert x time array, we may want to look at the first time step,
-which is the end of the previous simulation (time=6days) which we didn't store output for.
-```@example howtorun
-time = 1
-vor = Matrix{Float32}(ds["vor"][:,:,1,time]) # convert from Matrix{Union{Missing,Float32}} to Matrix{Float32}
-lat = ds["lat"][:]
-lon = ds["lon"][:]
-pcolormesh(lon,lat,vor')
-xlabel("longitude")
-ylabel("latitude")
-tight_layout() # hide
-savefig("galewsky1.png") # hide
-nothing # hide
-```
-![Galewsky jet pyplot](galewsky1.png)
-
-You see that in comparison the unicode plot heavily coarse-grains the simulation, well it's unicode after all!
-And now the last time step, that means time = 12days is
-
-```@example howtorun
-time = 25
-vor = Matrix{Float32}(ds["vor"][:,:,1,time])
-pcolormesh(lon,lat,vor')
-xlabel("longitude")
-ylabel("latitude")
-tight_layout() # hide
-savefig("galewsky2.png") # hide
-nothing # hide
-```
-![Galewsky jet pyplot](galewsky2.png)
-
-The jet broke up into many small eddies, but the turbulence is still confined to the northern hemisphere, cool!
-How this may change when we add mountains (we had `NoOrography` above!), say Earth's orography, you may ask?
-Let's try it out! We create an `EarthOrography` struct like so
-
-```@example howtorun
-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
-
-```@example howtorun
-model = ShallowWaterModel(;spectral_grid, orography, initial_conditions)
-simulation = initialize!(model)
-run!(simulation,n_days=12,output=true)
-```
-
-This time the run got a new run id, `0002` in our case, which you can always check
-after the `run!` call (the automatic run id is only determined just before the main time loop starts)
-with `model.output.id`, but otherwise we do as before.
-
-```@example howtorun
-ds = NCDataset("run_0002/output.nc")
-time = 49
-vor = Matrix{Float32}(ds["vor"][:,:,1,time])
-pcolormesh(lon,lat,vor')
-xlabel("longitude")
-ylabel("latitude")
-tight_layout() # hide
-savefig("galewsky3.png") # hide
-nothing # hide
-```
-![Galewsky jet pyplot](galewsky3.png)
-
-Interesting! The initial conditions have zero velocity in the southern hemisphere, but still, one can see
-some imprint of the orography on vorticity. You can spot the coastline of Antarctica; the Andes and
-Greenland are somewhat visible too. Mountains also completely changed the flow after 12 days,
-probably not surprising!
+In the following we will describe these steps in more detail.
+If you want to start with some examples first, see [Model setups](@ref)
+which has easy and more advanced examples of how to set up
+SpeedyWeather.jl to run the simulation you want.
 
 ## SpectralGrid
 
 The life of every SpeedyWeather.jl simulation starts with a `SpectralGrid` object.
-We have seen some examples above, now let's look into the details
-
-```julia
-help?> SpectralGrid
-search: SpectralGrid
-
-  Defines the horizontal spectral resolution and corresponding grid and the
-  vertical coordinate for SpeedyWeather.jl. Options are
-
-    •  NF::Type{<:AbstractFloat}: number format used throughout the model
-
-    •  trunc::Int64: horizontal resolution as the maximum degree of
-       spherical harmonics
-
-    •  Grid::Type{<:SpeedyWeather.RingGrids.AbstractGrid}: horizontal
-       grid used for calculations in grid-point space
-
-    •  dealiasing::Float64: how to match spectral with grid resolution:
-       dealiasing factor, 1=linear, 2=quadratic, 3=cubic grid
-
-    •  radius::Float64: radius of the sphere [m]
-
-    •  nlat_half::Int64: number of latitude rings on one hemisphere
-       (Equator incl)
-
-    •  npoints::Int64: total number of grid points in the horizontal
-
-    •  nlev::Int64: number of vertical levels
-
-    •  vertical_coordinates::SpeedyWeather.VerticalCoordinates:
-       coordinates used to discretize the vertical
-
-  nlat_half and npoints should not be chosen but are derived from trunc, Grid
-  and dealiasing.
-```
-
-Say we wanted double precision (`Float64`), a spectral resolution of T42 on
-a regular longitude-latitude grid (`FullClenshawGrid`) with cubic truncation
-(`dealiasing=3`) and 4 vertical levels, we would do this by
-
-```@example howtorun
-spectral_grid = SpectralGrid(NF=Float64, trunc=42, Grid=FullClenshawGrid, dealiasing=3, nlev=4)
-```
-
-We don't specify the resolution of the grid (its `nlat_half` parameter) directly,
-instead we chose a spectral truncation `trunc`
-and through the `dealiasing` factor a grid resolution will be automatically chosen. 
-Here T42 will be a matched with a 192x95 regular longitude-latitude grid
-that has 18240 grid points in total. For details see [Matching spectral and grid resolution](@ref).
-
-We could have also defined `SpectralGrid` on a smaller sphere than Earth,
-or with a different vertical spacing
-```@example howtorun
-vertical_coordinates = SigmaCoordinates(0:0.2:1)
-```
-These are regularly spaced [Sigma coordinates](@ref), defined through their half levels.
-```@example howtorun
-spectral_grid = SpectralGrid(;vertical_coordinates,radius=1e6)
-```
-
-In the end, a `SpectralGrid` defines the physical domain of the simulation and its discretization.
+A `SpectralGrid` defines the physical domain of the simulation and its discretization.
 This domain has to be a sphere because of the spherical harmonics, but it can have a different radius.
 The discretization is for spectral, grid-point space and the vertical as this determines the size of many
 arrays for preallocation, for which als the number format is essential to know.
 That's why `SpectralGrid` is the beginning of every SpeedyWeather.jl simulation and that is why
 it has to be passed on to (most) model components.
 
-## References
+The default `SpectralGrid` is
+
+```@example howto
+using SpeedyWeather
+spectral_grid = SpectralGrid()
+```
+You can also get the help prompt by typing `?SpectralGrid`.
+Let's explain the details: The spectral resolution is T31, so the largest
+wavenumber in spectral space is 31, and all the complex spherical harmonic
+coefficients of a given 2D field (see [Spherical Harmonic Transform](@ref))
+are stored in a [`LowerTriangularMatrix`](@ref lowertriangularmatrices)
+in the number format Float32. The radius of the sphere is
+6371km by default, which is the radius of the Earth.
+
+This spectral resolution is combined with an
+[octahedral Gaussian grid](@ref OctahedralGaussianGrid) of 3168 grid points.
+This grid has 48 latitude rings, 20 longitude points around the poles
+and up to 96 longitude points around the Equator. Data on that
+grid is also stored in Float32. The resolution is therefore on average about 400km.
+In the vertical 8 levels are used, using [Sigma coordinates](@ref).
+
+The resolution of a SpeedyWeather.jl simulation is adjusted using the
+`trunc` argument, this defines the spectral resolution and the grid
+resolution is automatically adjusted to keep the aliasing between
+spectral and grid-point space constant (see [Matching spectral and grid resolution](@ref)).
+```@example howto
+spectral_grid = SpectralGrid(trunc=85)
+```
+Typical values are 31,42,63,85,127,170,... although you can technically
+use any integer, see [Available horizontal resolutions](@ref) for details.
+Now with T85 (which is a common notation for `trunc=85`) the grid
+is of higher resolution too. You may play with the `dealiasing` factor,
+a larger factor increases the grid resolution that is matched with a given
+spectral resolution. You don't choose the resolution of the grid directly,
+but using the `Grid` argument you can change its type (see [Grids](@ref))
+```@example howto
+spectral_grid = SpectralGrid(trunc=85,dealiasing=3,Grid=HEALPixGrid)
+```
+
+## Vertical coordinates and resolution
+
+The number of vertical layers or levels (we use both terms often interchangeably)
+is determined through the `nlev` argument. Especially for the
+`BarotropicModel` and the `ShallowWaterModel` you want to set this to
+```@example howto
+spectral_grid = SpectralGrid(nlev=1)
+```
+For a single vertical level the type of the vertical coordinates does not matter,
+but in general you can change the spacing of the sigma coordinates
+which have to be discretized in ``[0,1]``
+```@example howto
+vertical_coordinates = SigmaCoordinates(0:0.2:1)
+```
+These are regularly spaced [Sigma coordinates](@ref), defined through their half levels.
+The cell centers or called full levels are marked with an ×.
+
+## [Creating model components](@id create_model_components)
+
+SpeedyWeather.jl deliberately avoids the namelists that are the main
+user interface to many old school models. Instead, we employ a modular approach
+whereby every non-default model component is created with its respective
+parameters to finally build the whole model from these components.
+
+If you know which components you want to adjust you would create them right away,
+however, a new user might first want to know which components there are,
+so let's flip the logic around and assume you know you want to run a `ShallowWaterModel`.
+You can create a default `ShallowWaterModel` like so and inspect its components
+```@example howto
+model = ShallowWaterModel()
+```
+
+So by default the `ShallowWaterModel` uses a [Leapfrog `time_stepping`](@ref leapfrog),
+which you can inspect like so
+```@example howto
+model.time_stepping
+```
+
+Model components often contain parameters from the `SpectralGrid` as they are needed
+to determine the size of arrays and other internal reasons. You should, in most cases,
+just ignore those. But the `Leapfrog` time stepper comes with `Δt_at_T31` which
+is the parameter used to scale the time step automatically. This means at a spectral
+resolution of T31 it would use 30min steps, at T63 it would be ~half that, 15min, etc.
+Meaning that if you want to have a shorter or longer time step you can create a new
+`Leapfrog` time stepper. But remember that every model component depends on a
+`SpectralGrid` as first argument.
+```@example howto
+spectral_grid = SpectralGrid(trunc=63,nlev=1)
+time_stepping = Leapfrog(spectral_grid,Δt_at_T31=15)
+```
+The actual time step at the given resolution (here T63) is then `Δt_sec`, there's
+also `Δt` which is a scaled time step used internally, because SpeedyWeather.jl
+[scales the equations](@ref scaled_swm) with the radius of the Earth,
+but this is largely hidden (except here) from the user. With this new 
+`Leapfrog` time stepper constructed we can create a model by passing
+on the components (they are keyword arguments so either use `;time_stepping`
+for which the naming must match, or `time_stepping=my_time_stepping` with
+any name)
+```@example howto
+model = ShallowWaterModel(;spectral_grid,time_stepping)
+```
+This logic continues for all model components. See the [Model setups](@ref)
+for examples. All model components are also subtype (i.e. `<:`) of
+some abstract supertype, this feature can be made use of to redefine
+not just constant parameters of existing model components but also
+to define new ones. This is more elaborated in [Extending SpeedyWeather](@ref).
+
+A model in SpeedyWeather.jl is understood as a collection of model components that
+roughly belong into one of three groups. 
+
+1. Components (numerics, dynamics, or physics) that have parameters, possibly contain some data (immutable, precomputed) and some functions associated with them. For example, a forcing term may be defined through some parameters, but also require precomputed arrays, or data to be loaded from file and a function that instructs how to calculate this forcing on every time step.
+2. Components that are merely a collection of parameters that conceptually belong together. For example, `Earth` is an `AbstractPlanet` that contains all the orbital parameters important for weather and climate (rotation, gravity, etc).
+3. Components for output purposes. For example, `OutputWriter` controls the NetCDF output, and `Feedback` controls the information printed to the REPL and to file.
+
+## [Creating a model](@id create_model)
+
+SpeedyWeather.jl implements (currently) 4 different models
+
+1. `BarotropicModel`, see [Barotropic vorticity model](@ref)
+2. `ShallowWaterModel`, see [Shallow water model](@ref)
+3. `PrimitiveDryModel`, see [Primitive equation model](@ref) but with zero humidity.
+4. `PrimitiveWetModel`, see [Primitive equation model](@ref).
+
+Overall they are kept quite similar, but there are still fundamental
+differences arising from the different equations. For example,
+the barotropic and shallow water models do not have any physical
+parameterizations. Conceptually you construct these different models with
 
-[^1] 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)
\ No newline at end of file
+```julia
+spectral_grid = SpectralGrid(trunc=...,...)
+component1 = SomeComponent(spectral_grid,parameter1=...,...)
+component2 = SomeOtherComponent(spectral_grid,parameter2=...,...)
+model = BarotropicModel(;spectral_grid,all_other_components...,...)
+```
+or `model = ShallowWaterModel(;spectral_grid,...)`, etc.
+
+## [Model initialization](@id initialize)
+
+In the previous section the model was created, but this is conceptually
+just gathering all its components together. However, many components
+need to be initialized. This step is used to precompute arrays,
+load necessary data from file or to communicate those between components.
+Furthermore, prognostic and diagnostic variables are allocated.
+It is (almost) all that needs to be done before the model can be run
+(exception being the output initialization). Many model components
+have a `initialize!` function associated with them that it executed here.
+However, from a user perspective all that needs to be done here is
+```@example howto
+simulation = initialize!(model)
+```
+and we have initialized the `ShallowWaterModel` we have defined earlier.
+After this step you can continue to tweak your model setup but note that
+some model components are immutable, or that your changes may not be
+propagated to other model components that rely on it. But you can, for
+example, change the output time step (in hours) like so
+```@example howto
+simulation.model.output.output_dt = 1
+```
+Now, if there's output, it will be every hour. Furthermore the initial
+conditions can be set with the `initial_conditions` model component
+which are then set during `initialize!(::ModelSetup)`, but you can also
+change them now, before the model runs
+```@example howto
+simulation.prognostic_variables.layers[1].timesteps[1].vor[1] = 0
+```
+So with this we have set the zero mode of vorticity of the first (and only)
+layer in the shallow water to zero. Because the leapfrogging is a 2-step
+time stepping scheme we set here the first. As it is often tricky
+to set the initial conditions in spectral space, it is generally advised
+to do so through the `initial_conditions` model component.
+
+## [Run a simulation](@id run)
+
+After creating a model, initializing it to a simulation, all that is left
+is to run the simulation.
+```@example howto
+run!(simulation)
+```
+By default this runs for 10 days without output. These are the options left
+to change, so with
+```@example howto
+model.output.id = "test" # hide
+run!(simulation,n_days=5,output=true)
+```
+You would continue this simulation (the previous `run!` call already integrated
+10 days!) for another 5 days and storing default [NetCDF output](@ref).
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 567642dce..ea73a8a6a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -26,6 +26,7 @@ See the following pages of the documentation for more details
 
 - [Installation](installation.md)
 - [How to run SpeedyWeather.jl](how_to_run_speedy.md)
+- [Model setups](setups.md)
 - [Spherical harmonic transform](spectral_transform.md)
 - [Grids](grids.md)
 - [Barotropic model](barotropic.md)
diff --git a/docs/src/primitiveequation.md b/docs/src/primitiveequation.md
index aa6ff991c..797148af2 100644
--- a/docs/src/primitiveequation.md
+++ b/docs/src/primitiveequation.md
@@ -2,7 +2,7 @@
 
 The [primitive equations](https://en.wikipedia.org/wiki/Primitive_equations) are a hydrostatic approximation
 of the compressible Navier-Stokes equations for an ideal gas on a rotating sphere. We largely follow
-the idealised spectral dynamical core developed by GFDL[^1] and documented therein[^2].
+the idealised spectral dynamical core developed by GFDL[^GFDL1] and documented therein[^GFDL2].
 
 The primitive equations solved by SpeedyWeather.jl for relative vorticity ``\zeta``, divergence ``\mathcal{D}``,
 logarithm of surface pressure ``\ln p_s``, temperature ``T`` and specific humidity ``q`` are
@@ -57,7 +57,7 @@ p_w &= \rho_w R_w T \\
 \end{aligned}
 ```
 with the respective specific gas constants ``R_d = R/m_d`` and ``R_w = R/m_w``
-obtained from the univeral gas constant ``R`` divided by the molecular masses
+obtained from the universal gas constant ``R`` divided by the molecular masses
 of the gas. The total pressure ``p`` in the air parcel is
 ```math
 p = p_d + p_w = (\rho_d R_d + \rho_w R_w)T
@@ -67,7 +67,7 @@ using the ideal gas law, with the temperature ``T``, so that we never have
 to calculate the density explicitly. However, in order to not deal with
 two densities (dry air and water vapour) we would like to replace
 temperature with a virtual temperature that includes the effect of
-humidity on the density. So, whereever we use the ideal gas law
+humidity on the density. So, wherever we use the ideal gas law
 to replace density with temperature, we would use the virtual temperature,
 which is a function of the absolute temperature and specific humidity,
 instead. A higher specific humidity in an air parcel lowers 
@@ -78,7 +78,7 @@ think of the virtual temperature as the temperature dry air would need to have
 to be as light as moist air.
 
 Starting with the last equation, with some manipulation we can write
-the ideal gas law as total density ``rho`` times a gas constant
+the ideal gas law as total density ``\rho`` times a gas constant
 times the virtual temperature that is supposed to be a function
 of absolute temperature, humidity and some constants
 ```math
@@ -133,7 +133,7 @@ The chain rule lets us differentiate ``\Psi`` with respect to ``z`` or ``\eta``
 \frac{\partial \Psi}{\partial z} = \frac{\partial \Psi}{\partial \eta}\frac{\partial \eta}{\partial z},
 \qquad \frac{\partial \Psi}{\partial \eta} = \frac{\partial \Psi}{\partial z}\frac{\partial z}{\partial \eta}
 ```
-But for derivatives with respect to ``x,y,t`` we have to apply the multivariable
+But for derivatives with respect to ``x,y,t`` we have to apply the multi-variable
 chain-rule as both ``\Psi`` and ``\eta`` depend on it. So a derivative with respect to
 ``x`` on ``\eta`` levels (where ``\eta`` constant) becomes
 ```math
@@ -311,7 +311,7 @@ for layer ``k``
 - (M_{k+\tfrac{1}{2}}T_{k+\tfrac{1}{2}} - M_{k-\tfrac{1}{2}}T_{k-\tfrac{1}{2}})
 ```
 which can be through the gradient product rule, and using the conservation of mass
-(see [Vertical velocity](@ref)) transformed into an advective form. In sigma coordinates this simpifies to
+(see [Vertical velocity](@ref)) transformed into an advective form. In sigma coordinates this simplifies to
 ```math
 \frac{\partial T_k}{\partial t} = - \mathbf{u}_k \cdot \nabla T_k
 - \frac{1}{\Delta \sigma_k}\left(\dot{\sigma}_{k+\tfrac{1}{2}}(T_{k+\tfrac{1}{2}} - T_k) - \dot{\sigma}_{k-\tfrac{1}{2}}(T_k - T_{k-\tfrac{1}{2}})\right)
@@ -420,7 +420,7 @@ From left to right: The pressure gradient force in ``z``-coordinates; in pressur
 and in sigma coordinates. Each denoted with the respective subscript on gradients. 
 SpeedyWeather.jl uses the latter.
 In sigma coordinates we may drop the ``\sigma`` subscript on gradients, but still meaning
-that the gradient is evaluted on a surface of our vertical coordinate.
+that the gradient is evaluated on a surface of our vertical coordinate.
 In vorticity-divergence formulation of the momentum equations the ``\nabla_\sigma \Phi``
 drops out in the vorticity equation (``\nabla \times \nabla \Phi = 0``),
 but becomes a ``-\nabla^2 \Phi`` in the divergence equation,
@@ -488,7 +488,7 @@ with ``\mathbf{u}_\perp = (v,-u)`` the rotated velocity vector, see [Barotropic
 ## Humidity equation
 
 The dynamical core treats humidity as an (active) tracer, meaning that after the physical
-parameterizations for humidity ``\matchal{P}`` are calculated in grid-point space,
+parameterizations for humidity ``\mathcal{P}`` are calculated in grid-point space,
 humidity is only advected with the flow. The only exception is the [Virtual temperature](@ref)
 as high levels of humidity will lower the effective density, which is why we use the
 virtual instead of the absolute temperature. The equation to be solved for humidity
@@ -636,7 +636,7 @@ So for every linear term in ``N_I`` we have two options corresponding to two sid
 2. Or, evaluate it at the current time step ``i`` as ``N(V_i)``, but then move it back to the previous time step ``i-1`` by adding (in spectral space) the linear operator ``N_I`` evaluated with the difference between the two time steps.
 
 If there is a tendency that is easily evaluated in spectral space it is easier to follow 1. However,
-a term that is costly to evalute in grid-point space should usually follow the latter. The reason is that
+a term that is costly to evaluate in grid-point space should usually follow the latter. The reason is that
 the previous time step is generally not available in grid-point space (unless recalculated through a costly transform
 or stored with additional memory requirements) so it is easier to follow 2 where the ``N_I`` is
 available in spectral space. For the adiabatic conversion term in the [Temperature equation](@ref)
@@ -738,7 +738,7 @@ be recomputed. Otherwise the algorithm for the semi-implicit scheme is as follow
 Then for every time step
 
 1. Compute the uncorrected tendencies evaluated at the current time step for the explicit terms and the previous time step for the implicit terms.
-2. Excpetion in SpeedyWeather.jl is the adiabatic conversion term, which is, using ``\mathbf{L}`` moved afterwards from the current ``i`` to the previous time step ``i-1``.
+2. Exception in SpeedyWeather.jl is the adiabatic conversion term, which is, using ``\mathbf{L}`` moved afterwards from the current ``i`` to the previous time step ``i-1``.
 3. Compute the combined tendency ``G`` from the uncorrected tendencies ``G_\mathcal{D}``, ``G_T``, ``G_{\ln p_s}``.
 4. With the inverted operator get the corrected tendency for divergence, ``\delta \mathcal{D} = \mathbf{S}^{-1}G``.
 5. Obtain the corrected tendencies for temperature ``\delta T`` and surface pressure ``\delta \ln p_s`` from ``\delta \mathcal{D}``.
@@ -804,8 +804,8 @@ Now loop over
 
 ## References
 
-[^1]: Geophysical Fluid Dynamics Laboratory, [Idealized models with spectral dynamics](https://www.gfdl.noaa.gov/idealized-models-with-spectral-dynamics/)
-[^2]: Geophysical Fluid Dynamics Laboratory, [The Spectral Dynamical Core](https://www.gfdl.noaa.gov/wp-content/uploads/files/user_files/pjp/spectral_core.pdf)
+[^GFDL1]: Geophysical Fluid Dynamics Laboratory, [Idealized models with spectral dynamics](https://www.gfdl.noaa.gov/idealized-models-with-spectral-dynamics/)
+[^GFDL2]: Geophysical Fluid Dynamics Laboratory, [The Spectral Dynamical Core](https://www.gfdl.noaa.gov/wp-content/uploads/files/user_files/pjp/spectral_core.pdf)
 [^Vallis]: GK Vallis, 2006. [Atmopsheric and Ocean Fluid Dynamics](http://vallisbook.org/), Cambridge University Press.
 [^SB81]: Simmons and Burridge, 1981. *An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates*, Monthly Weather Review. DOI: [10.1175/1520-0493(1981)109<0758:AEAAMC>2.0.CO;2](https://doi.org/10.1175/1520-0493(1981)109<0758:AEAAMC>2.0.CO;2).
 [^HS75]: Hoskins and Simmons, 1975. *A multi-layer spectral model and the semi-implicit method*, Quart. J. R. Met. Soc. DOI: [10.1002/qj.49710142918](https://doi.org/10.1002/qj.49710142918)
\ No newline at end of file
diff --git a/docs/src/setups.md b/docs/src/setups.md
new file mode 100644
index 000000000..e0e76bd5f
--- /dev/null
+++ b/docs/src/setups.md
@@ -0,0 +1,361 @@
+# Model setups
+
+The following is a collection of model setups, starting with an easy setup
+of the [Barotropic vorticity equation](@ref) and continuing with more
+complicated setups.
+
+## 2D turbulence on a non-rotating sphere
+
+!!! info "Setup script"
+    ```julia
+    using SpeedyWeather
+    spectral_grid = SpectralGrid(trunc=63,nlev=1)
+    still_earth = Earth(rotation=0)
+    initial_conditions = StartWithRandomVorticity()
+    model = BarotropicModel(;spectral_grid, initial_conditions, planet=still_earth)
+    simulation = initialize!(model)
+    run!(simulation,n_days=20)
+    ```
+
+We want to use the barotropic model to simulate some free-decaying 2D turbulence
+on the sphere without rotation. We start by defining the `SpectralGrid` object.
+To have a resolution of about 200km, we choose a spectral resolution of
+T63 (see [Available horizontal resolutions](@ref)) and `nlev=1` vertical levels.
+The `SpectralGrid` object will provide us with some more information
+```@example barotropic_setup
+using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=63,nlev=1)
+```
+Next step we create a planet that's like Earth but not rotating
+```@example barotropic_setup
+still_earth = Earth(rotation=0)
+```
+There are other options to create a planet but they are irrelevant for the
+barotropic vorticity equations. We also want to specify the initial conditions,
+randomly distributed vorticity is already defined
+```@example barotropic_setup
+using Random # hide
+Random.seed!(1234) # hide
+initial_conditions = StartWithRandomVorticity()
+```
+By default, the power of vorticity is spectrally distributed with ``k^{-3}``, ``k`` being the
+horizontal wavenumber, and the amplitude is ``10^{-5}\text{s}^{-1}``.
+
+Now we want to construct a `BarotropicModel`
+with these
+```@example barotropic_setup
+model = BarotropicModel(;spectral_grid, initial_conditions, planet=still_earth)
+nothing # hide
+```
+The `model` contains all the parameters, but isn't initialized yet, which we can do
+with and then run it. The `run!` command will always return the prognostic variables, which, by default, are 
+plotted for surface relative vorticity with a unicode plot. The resolution of the plot
+is not necessarily representative but it lets us have a quick look at the result
+```@example barotropic_setup
+simulation = initialize!(model)
+run!(simulation,n_days=20)
+```
+
+Woohoo! Something is moving! You could pick up where this simulation stopped by simply
+doing `run!(simulation,n_days=50)` again. We didn't store any output, which
+you can do by `run!(simulation,output=true)`, which will switch on NetCDF output
+with default settings. More options on output in [NetCDF output](@ref).
+
+## Shallow water with mountains
+
+!!! info "Setup script"
+    ```julia
+    using SpeedyWeather
+    spectral_grid = SpectralGrid(trunc=63,nlev=1)
+    orography = NoOrography(spectral_grid)
+    initial_conditions = ZonalJet()
+    model = ShallowWaterModel(;spectral_grid, orography, initial_conditions)
+    simulation = initialize!(model)
+    run!(simulation,n_days=6)
+    ```
+
+As a second example, let's investigate the Galewsky et al.[^G04] test case for the shallow
+water equations with and without mountains. As the shallow water system has also only
+one level, we can reuse the `SpectralGrid` from Example 1.
+```@example galewsky_setup
+using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=63,nlev=1)
+```
+Now as a first simulation, we want to disable any orography, so we create a `NoOrography`
+```@example galewsky_setup
+orography = NoOrography(spectral_grid)
+```
+Although the orography is zero, you have to pass on `spectral_grid` so that it can
+still initialize zero-arrays of the correct size and element type. Awesome.
+This time the initial conditions should be set the the Galewsky et al.[^G04] zonal
+jet, which is already defined as
+```@example galewsky_setup
+initial_conditions = ZonalJet()
+```
+The jet sits at 45˚N with a maximum velocity of 80m/s and a perturbation as described in their paper.
+Now we construct a model, but this time a `ShallowWaterModel`
+```@example galewsky_setup
+model = ShallowWaterModel(;spectral_grid, orography, initial_conditions)
+simulation = initialize!(model)
+run!(simulation,n_days=6)
+```
+Oh yeah. That looks like the wobbly jet in their paper. Let's run it again for another 6 days
+but this time also store [NetCDF output](@ref).
+```@example galewsky_setup
+run!(simulation,n_days=6,output=true)
+```
+The progress bar tells us that the simulation run got the identification "0001"
+(which just counts up, so yours might be higher), meaning that
+data is stored in the folder `/run_0001`, so let's plot that data properly (and not just using UnicodePlots).
+```@example galewsky_setup
+using PythonPlot, NCDatasets
+ioff() # hide
+ds = NCDataset("run_0001/output.nc")
+ds["vor"]
+```
+Vorticity `vor` is stored as a lon x lat x vert x time array, we may want to look at the first time step,
+which is the end of the previous simulation (time=6days) which we didn't store output for.
+```@example galewsky_setup
+t = 1
+vor = Matrix{Float32}(ds["vor"][:,:,1,t]) # convert from Matrix{Union{Missing,Float32}} to Matrix{Float32}
+lat = ds["lat"][:]
+lon = ds["lon"][:]
+
+fig,ax = subplots(1,1,figsize=(10,6))
+ax.pcolormesh(lon,lat,vor')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Relative vorticity")
+tight_layout() # hide
+savefig("galewsky1.png") # hide
+nothing # hide
+```
+![Galewsky jet pyplot1](galewsky1.png)
+
+You see that in comparison the unicode plot heavily coarse-grains the simulation, well it's unicode after all!
+And now the last time step, that means time = 12days is
+
+```@example galewsky_setup
+t = ds.dim["time"]
+vor = Matrix{Float32}(ds["vor"][:,:,1,t])
+ax.pcolormesh(lon,lat,vor')
+savefig("galewsky2.png") # hide
+nothing # hide
+```
+![Galewsky jet pyplot2](galewsky2.png)
+
+The jet broke up into many small eddies, but the turbulence is still confined to the northern hemisphere, cool!
+How this may change when we add mountains (we had `NoOrography` above!), say Earth's orography, you may ask?
+Let's try it out! We create an `EarthOrography` struct like so
+
+```@example galewsky_setup
+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
+
+```@example galewsky_setup
+model = ShallowWaterModel(;spectral_grid, orography, initial_conditions)
+simulation = initialize!(model)
+run!(simulation,n_days=12,output=true)
+```
+
+This time the run got a new run id, which you see in the progress bar, but can also always check
+after the `run!` call (the automatic run id is only determined just before the main time loop starts)
+with `model.output.id`, but otherwise we do as before.
+```@example galewsky_setup
+id = model.output.id
+```
+
+```@example galewsky_setup
+ds = NCDataset("run_$id/output.nc")
+time = 49
+vor = Matrix{Float32}(ds["vor"][:,:,1,time])
+
+fig,ax = subplots(1,1,figsize=(10,6))
+ax.pcolormesh(lon,lat,vor')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Relative vorticity")
+tight_layout() # hide
+savefig("galewsky3.png") # hide
+nothing # hide
+```
+![Galewsky jet pyplot3](galewsky3.png)
+
+Interesting! The initial conditions have zero velocity in the southern hemisphere, but still, one can see
+some imprint of the orography on vorticity. You can spot the coastline of Antarctica; the Andes and
+Greenland are somewhat visible too. Mountains also completely changed the flow after 12 days,
+probably not surprising!
+
+## Polar jet streams in shallow water
+
+Setup script:
+```@example jet_stream_setup
+using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=63,nlev=1)
+forcing = JetStreamForcing(spectral_grid,latitude=60)
+drag = QuadraticDrag(spectral_grid)
+output = OutputWriter(spectral_grid,ShallowWater,output_dt=6,output_vars=[:u,:v,:pres,:orography])
+model = ShallowWaterModel(;spectral_grid,output,drag,forcing)
+simulation = initialize!(model)
+model.feedback.verbose = false # hide
+run!(simulation,n_days=20)   # discard first 20 days   
+run!(simulation,n_days=20,output=true)
+nothing # hide
+```
+
+We want to simulate polar jet streams in the shallow water model. We add a `JetStreamForcing`
+that adds momentum at 60˚N to inject kinetic energy into the model. This energy needs to be removed
+(the [diffusion](@ref diffusion) is likely not sufficient) through a drag, we have implemented
+a `QuadraticDrag` and use the default drag coefficient. Outputting ``u,v,\eta`` (called `:pres`,
+as it is the pressure equivalent in the shallow water system) we run 20 days without output
+to give the system some time to adapt to the forcing. And visualize zonal wind after another
+20 days with
+
+```@example jet_stream_setup
+using PythonPlot, NCDatasets
+ioff() # hide
+
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+timestep = ds.dim["time"]
+u = Matrix{Float32}(ds["u"][:,:,1,timestep])
+lat = ds["lat"][:]
+lon = ds["lon"][:]
+
+fig,ax = subplots(1,1,figsize=(10,6))
+q = ax.pcolormesh(lon,lat,u')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Zonal wind [m/s]")
+colorbar(q,ax=ax)
+tight_layout() # hide
+savefig("polar_jets.png") # hide
+nothing # hide
+```
+![Polar jets pyplot](polar_jets.png)
+
+
+## Gravity waves on the sphere
+
+Setup script:
+```@example gravity_wave_setup
+using Random # hide
+Random.seed!(1234) # hide
+using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=127,nlev=1)
+time_stepping = SpeedyWeather.Leapfrog(spectral_grid,Δt_at_T31=30)
+implicit = SpeedyWeather.ImplicitShallowWater(spectral_grid,α=0.5)
+orography = EarthOrography(spectral_grid,smoothing=false)
+initial_conditions = SpeedyWeather.RandomWaves()
+output = OutputWriter(spectral_grid,ShallowWater,output_dt=12,output_vars=[:u,:pres,:div,:orography])
+model = ShallowWaterModel(;spectral_grid,orography,output,initial_conditions,implicit,time_stepping)
+simulation = initialize!(model)
+model.feedback.verbose = false # hide
+run!(simulation,n_days=2,output=true)
+nothing # hide
+```
+
+How are gravity waves propagating around the globe? We want to use the shallow water model
+to start with some random perturbations of the interface displacement (the "sea surface height")
+but zero velocity and let them propagate around the globe. We set the ``\alpha`` parameter
+of the [semi-implicit time integration](@ref implicit_swm) to ``0.5`` to have a centred
+implicit scheme which dampens the gravity waves less than a backward implicit scheme would do.
+But we also want to keep orography, and particularly no smoothing on it, to have the orography
+as rough as possible. The initial conditions are set to `RandomWaves` which set the spherical
+harmonic coefficients of ``\eta`` to between given wavenumbers to some random values
+```@example gravity_wave_setup
+SpeedyWeather.RandomWaves()
+```
+so that the amplitude `A` is as desired, here 2000m. Our layer thickness is by default
+```@example gravity_wave_setup
+model.atmosphere.layer_thickness
+```
+8.5km so those waves are with an amplitude of 2000m quite strong.
+But the semi-implicit time integration can handle that even with fairly large time steps of
+```@example gravity_wave_setup
+model.time_stepping.Δt_sec
+```
+seconds. Note that the gravity wave speed here is ``\sqrt{gH}`` so almost 300m/s.
+Let us also output divergence, as gravity waves are quite pronounced in that variable.
+But given the speed of gravity waves we don't have to integrate for long.
+Visualise with
+
+
+```@example gravity_wave_setup
+using PythonPlot, NCDatasets
+ioff() # hide
+
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+timestep = ds.dim["time"]
+div = Matrix{Float32}(ds["div"][:,:,1,timestep])
+lat = ds["lat"][:]
+lon = ds["lon"][:]
+
+fig,ax = subplots(1,1,figsize=(10,6))
+ax.pcolormesh(lon,lat,div')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Divergence")
+tight_layout() # hide
+savefig("gravity_waves.png") # 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 = ZonalWind()
+model = PrimitiveDryModel(;spectral_grid,orography,initial_conditions,physics=false)
+simulation = initialize!(model)
+model.feedback.verbose = false # hide
+run!(simulation,n_days=9,output=true)
+nothing # hide
+```
+
+The Jablonowski-Williamson baroclinic wave test case[^JW06] using the [Primitive equation model](@ref)
+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 in `ZonalWind`, 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 PythonPlot, NCDatasets
+ioff() # hide
+
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+timestep = ds.dim["time"]
+surface = ds.dim["lev"]
+vor = Matrix{Float32}(ds["vor"][:,:,surface,timestep])
+lat = ds["lat"][:]
+lon = ds["lon"][:]
+
+fig,ax = subplots(1,1,figsize=(10,6))
+ax.pcolormesh(lon,lat,vor')
+ax.set_xlabel("longitude")
+ax.set_ylabel("latitude")
+ax.set_title("Surface relative vorticity")
+tight_layout() # hide
+savefig("jablonowski.png") # hide
+nothing # hide
+```
+![Jablonowski pyplot](jablonowski.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)
\ No newline at end of file
diff --git a/docs/src/shallowwater.md b/docs/src/shallowwater.md
index c9d4499d3..a9f43614b 100644
--- a/docs/src/shallowwater.md
+++ b/docs/src/shallowwater.md
@@ -106,7 +106,7 @@ In spectral space ``\nabla^2`` is a diagonal operator, meaning that there is no
 harmonics and its inversion is therefore easily done on a mode-by-mode basis of the harmonics.
 
 This can be made use of when facing time stepping constraints with explicit schemes, where
-ridiculuously small time steps to resolve fast waves would otherwise result in a horribly slow simulation. 
+ridiculously small time steps to resolve fast waves would otherwise result in a horribly slow simulation. 
 In the shallow water system there are gravity waves that propagate at a wave speed of ``\sqrt{gH}``
 (typically 300m/s), which, in order to not violate the
 [CFL criterion](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition)
@@ -219,7 +219,7 @@ The idea of the semi-implicit time stepping is now as follows:
 Some notes on the semi-implicit time stepping
 
 - The inversion of the semi-implicit time stepping depends on ``\delta t``, that means every time the time step changes, the inversion has to be recalculated.
-- You may choose ``\alpha = 1/2`` to dampen gravity waves but initialisation shocks still usually kick off many gravity waves that propagate around the sphere for many days.
+- You may choose ``\alpha = 1/2`` to dampen gravity waves but initialization shocks still usually kick off many gravity waves that propagate around the sphere for many days.
 - With increasing ``\alpha > 1/2`` these waves are also slowed down, such that for ``\alpha = 1`` they quickly disappear in several hours.
 - Using the [scaled shallow water equations](@ref scaled_swm) the time step ``\delta t`` has to be the scaled time step ``\tilde{\Delta t} = \delta t/R`` which is divided by the radius ``R``. Then we use the normalized eigenvalues ``-l(l+1)`` which also omit the ``1/R^2`` scaling, see [scaled shallow water equations](@ref scaled_swm) for more details.
 
diff --git a/docs/src/spectral_transform.md b/docs/src/spectral_transform.md
index a24afd01d..9825a6737 100644
--- a/docs/src/spectral_transform.md
+++ b/docs/src/spectral_transform.md
@@ -227,11 +227,11 @@ Some remarks on this table
 - This assumes the default quadratic truncation, you can always adapt the grid resolution via the `dealiasing` option, see [Matching spectral and grid resolution](@ref)
 - `nlat` refers to the total number of latitude rings, see [Grids](@ref). With non-Gaussian grids, `nlat` will be one one less, e.g. 47 instead of 48 rings.
 - `nlon` is the number of longitude points on the [Full Gaussian Grid](@ref FullGaussianGrid), for other grids there will be at most these number of points around the Equator.
-- ``\Delta x`` is the horizontal resolution. For a spectral model there are many ways of estimating this[^9]. We use here the square root of the average area a grid cell covers, see [Effective grid resolution](@ref)
+- ``\Delta x`` is the horizontal resolution. For a spectral model there are many ways of estimating this[^Randall2021]. We use here the square root of the average area a grid cell covers, see [Effective grid resolution](@ref)
 
 ## Effective grid resolution
 
-There are many ways to estimate the effective grid resolution of spectral models[^9].
+There are many ways to estimate the effective grid resolution of spectral models[^Randall2021].
 Some of them are based on the wavelength a given spectral resolution allows to
 represent, others on the total number of real variables per area.
 However, as many atmospheric models do represent a considerable amount of physics
@@ -354,7 +354,7 @@ number ``m`` times imaginary ``i``.
 ### Meridional derivative
 
 The meridional derivative of the spherical harmonics is a derivative of the Legendre polynomials for which the following
-recursion relation applies[^10],[^11]
+recursion relation applies[^Randall2021],[^Durran2010],[^GFDL],[^Orszag70]
 
 ```math
 \cos\theta \frac{dP_{l,m}}{d\theta} = -l\epsilon_{l+1,m}P_{l+1,m} + (l+1)\epsilon_{l,m}P_{l-1,m}.
@@ -394,7 +394,7 @@ degree ``l``, but scalar quantities should not make use of it. Equivalently, the
 set to zero before the time integration, which only advances scalar quantities.
 
 
-In SpeedyWeather.jl vector quantities
+In SpeedyWeather.jl, vector quantities
 like ``u,v`` use therefore one more meridional mode than scalar quantities such as vorticity ``\zeta`` or stream
 function ``\Psi``. The meridional derivative in SpeedyWeather.jl also omits the ``1/R``-scaling as explained for
 the [Zonal derivative](@ref) and in [Radius scaling](@ref scaling).
@@ -404,7 +404,7 @@ the [Zonal derivative](@ref) and in [Radius scaling](@ref scaling).
 The meridional gradient as described above can be applied to scalars, such as ``\Psi`` and ``\Phi`` in the conversion
 to velocities ``(u,v) = \nabla^\bot\Psi + \nabla\Phi``, however, the operators curl ``\nabla \times`` and divergence
 ``\nabla \cdot`` in spherical coordinates involve a ``\cos\theta`` scaling _before_ the meridional gradient is applied.
-How to translate this to spectral coefficients has to be derived separately[^10],[^11].
+How to translate this to spectral coefficients has to be derived separately[^Randall2021],[^Durran2010].
 
 The spectral transform of vorticity ``\zeta`` is
 ```math
@@ -414,7 +414,7 @@ P_{l,m}(\sin\theta) e^{im\lambda} d\lambda \cos\theta d\theta
 Given that ``R\zeta = \cos^{-1}\partial_\lambda v - \cos^{-1}\partial_\theta (u \cos\theta)``,
 we therefore have to evaluate a meridional integral of the form
 ```math
-\int P_{l,m} \frac{1}{\cos \theta} \partial_\theta(u \cos\theta)) \cos \theta d\theta
+\int P_{l,m} \frac{1}{\cos \theta} \partial_\theta(u \cos\theta) \cos \theta d\theta
 ```
 which can be solved through integration by parts. As ``u\cos\theta = 0`` at ``\theta = \pm \tfrac{\pi}{2}`` only the integral
 ```math
@@ -443,7 +443,7 @@ The spectral Laplacian is easily applied to the coefficients ``\Psi_{lm}`` of a
 as the spherical harmonics are eigenfunctions of the Laplace operator ``\nabla^2`` in spherical
 coordinates with eigenvalues ``-l(l+1)`` divided by the radius squared ``R^2``, i.e.
 ``\nabla^2 \Psi`` becomes ``\tfrac{-l(l+1)}{R^2}\Psi_{lm}`` in spectral space. For example,
-vorticity ``\zeta`` and streamfunction ``\Psi`` are related by ``\zeta = \nabla^2\Psi``
+vorticity ``\zeta`` and stream function ``\Psi`` are related by ``\zeta = \nabla^2\Psi``
 in the barotropic vorticity model. Hence, in spectral space this is equivalent for every
 spectral mode of degree ``l`` and order ``m`` to
 
@@ -502,3 +502,4 @@ as further described in [Radius scaling](@ref scaling).
 [^GFDL]: Geophysical Fluid Dynamics Laboratory, [The barotropic vorticity equation](https://www.gfdl.noaa.gov/wp-content/uploads/files/user_files/pjp/barotropic.pdf).
 [^FFT]: Depending on the implementation of the Fast Fourier Transform ([Cooley-Tukey algorithm](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm), or or the [Bluestein algorithm](https://en.wikipedia.org/wiki/Chirp_Z-transform#Bluestein.27s_algorithm)) *easily Fourier-transformable* can mean different things: Vectors of the length ``n`` that is a power of two, i.e. ``n = 2^i`` is certainly easily Fourier-transformable, but for most FFT implementations so are ``n = 2^i3^j5^k`` with ``i,j,k`` some positive integers. In fact, [FFTW](http://fftw.org/) uses ``O(n \log n)`` algorithms even for prime sizes.
 [^Bourke72]: Bourke, W. An Efficient, One-Level, Primitive-Equation Spectral Model. Mon. Wea. Rev. 100, 683–689 (1972). doi:[10.1175/1520-0493(1972)100<0683:AEOPSM>2.3.CO;2](https://doi.org/10.1175/1520-0493(1972)100<0683:AEOPSM>2.3.CO;2)
+[^Orszag70]: Orszag, S. A., 1970: Transform Method for the Calculation of Vector-Coupled Sums: Application to the Spectral Form of the Vorticity Equation. J. Atmos. Sci., 27, 890–895, [10.1175/1520-0469(1970)027<0890:TMFTCO>2.0.CO;2](https://doi.org/10.1175/1520-0469(1970)027<0890:TMFTCO>2.0.CO;2). 
\ No newline at end of file
diff --git a/docs/src/speedytransforms.md b/docs/src/speedytransforms.md
index 9da816e8e..47cea3e60 100644
--- a/docs/src/speedytransforms.md
+++ b/docs/src/speedytransforms.md
@@ -38,7 +38,7 @@ map = gridded(alms)
 ```
 By default, the `gridded` transforms onto a [`FullGaussianGrid`](@ref FullGaussianGrid) unravelled here
 into a vector west to east, starting at the prime meridian, then north to south, see [RingGrids](@ref).
-We can visualize `map` quickly with a unicodeplot via `plot` (see [Visualising RingGrid data](@ref))
+We can visualize `map` quickly with a UnicodePlot via `plot` (see [Visualising RingGrid data](@ref))
 ```@example speedytransforms
 plot(map)
 ```
@@ -108,7 +108,7 @@ Note that because we chose `dealiasing=3` (cubic truncation) we now match a T5 s
 with a 12-ring octahedral Gaussian grid, instead of the 8 rings as above. So going from
 `dealiasing=2` (default) to `dealiasing=3` increased our resolution on the grid while the
 spectral resolution remains the same. The `SpectralTransform` also has options for the
-recomputation or precomputation of the Legendre polynomials, fore more information see
+recomputation or pre-computation of the Legendre polynomials, fore more information see
 [(P)recompute Legendre polynomials](@ref).
 
 Passing on `S` the `SpectralTransform` now allows us to transform directly on the grid
@@ -164,7 +164,7 @@ m = Matrix(map) # hide
 m
 ```
 You hopefully know which grid this data comes on, let us assume it is a regular
-latitde-longitude grid, which we call the `FullClenshawGrid`. Note that for the spectral
+latitude-longitude grid, which we call the `FullClenshawGrid`. Note that for the spectral
 transform this should not include the poles, so the 96x47 matrix size here corresponds
 to 
 
diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl
index 491cc713b..5a91e7a53 100644
--- a/src/dynamics/models.jl
+++ b/src/dynamics/models.jl
@@ -21,7 +21,7 @@ whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
 Base.@kwdef struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Barotropic
     "dictates resolution for many other components"
-    spectral_grid::SpectralGrid = SpectralGrid()
+    spectral_grid::SpectralGrid = SpectralGrid(nlev=1)
 
     # DYNAMICS
     "contains physical and orbital characteristics"
@@ -77,7 +77,7 @@ whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
 Base.@kwdef struct ShallowWaterModel{NF<:AbstractFloat, D<:AbstractDevice} <: ShallowWater
     "dictates resolution for many other components"
-    spectral_grid::SpectralGrid = SpectralGrid()
+    spectral_grid::SpectralGrid = SpectralGrid(nlev=1)
 
     # DYNAMICS
     "contains physical and orbital characteristics"
diff --git a/src/dynamics/prognostic_variables.jl b/src/dynamics/prognostic_variables.jl
index 5d3113e8c..b79c840b4 100644
--- a/src/dynamics/prognostic_variables.jl
+++ b/src/dynamics/prognostic_variables.jl
@@ -13,6 +13,12 @@ Base.@kwdef mutable struct Clock
     n_timesteps::Int = 0  
 end
 
+function Base.show(io::IO,C::Clock)
+    println(io,"$(typeof(C))")
+    keys = propertynames(C)
+    print_fields(io,C,keys)
+end
+
 """
 $(TYPEDSIGNATURES)
 Initialize the clock with the time step `Δt` in the `time_stepping`."""
diff --git a/src/dynamics/vertical_coordinates.jl b/src/dynamics/vertical_coordinates.jl
index 14558f41c..d8d8d701f 100644
--- a/src/dynamics/vertical_coordinates.jl
+++ b/src/dynamics/vertical_coordinates.jl
@@ -15,15 +15,15 @@ SigmaCoordinates(σ_half::AbstractVector) = SigmaCoordinates(nlev=length(σ_half
 SigmaCoordinates(σ_half::AbstractRange) = SigmaCoordinates(collect(σ_half))
 
 function Base.show(io::IO,σ::SigmaCoordinates)
-    println("$(σ.nlev)-level SigmaCoordinates")
+    println(io,"$(σ.nlev)-level SigmaCoordinates")
     nchars = length(string(σ.nlev))
     format = Printf.Format("%$(nchars)d")
     for k=1:σ.nlev
-        println("├─ ", @sprintf("%1.4f",σ.σ_half[k]),"  k = ",Printf.format(format,k-1),".5")
+        println(io,"├─ ", @sprintf("%1.4f",σ.σ_half[k]),"  k = ",Printf.format(format,k-1),".5")
         σk = (σ.σ_half[k] + σ.σ_half[k+1])/2
-        println("│× ",@sprintf("%1.4f",σk),"  k = ",Printf.format(format,k))
+        println(io,"│× ",@sprintf("%1.4f",σk),"  k = ",Printf.format(format,k))
     end
-    print("└─ ",@sprintf("%1.4f",σ.σ_half[end]),"  k = ",Printf.format(format,σ.nlev),".5")
+    print(io,"└─ ",@sprintf("%1.4f",σ.σ_half[end]),"  k = ",Printf.format(format,σ.nlev),".5")
 end
 
 """Coefficients of the generalised logistic function to describe the vertical coordinate.
diff --git a/src/output/feedback.jl b/src/output/feedback.jl
index c4ee82214..ff6e15f35 100644
--- a/src/output/feedback.jl
+++ b/src/output/feedback.jl
@@ -29,6 +29,18 @@ mutable struct Feedback <: AbstractFeedback
     nars_detected::Bool                     
 end
 
+function Base.show(io::IO,F::AbstractFeedback)
+    println(io,"$(typeof(F)) <: AbstractFeedback")
+    keys = propertynames(F)
+    print_fields(io,F,keys)
+end
+
+function Base.show(io::IO,P::ProgressMeter.Progress)
+    println(io,"$(typeof(P)) <: ProgressMeter.AbstractProgress")
+    keys = propertynames(P)
+    print_fields(io,P,keys)
+end
+
 """
 $(TYPEDSIGNATURES)
 Generator function for a Feedback struct."""
@@ -59,8 +71,9 @@ Initializes the a `Feedback` struct."""
 function initialize!(feedback::Feedback,clock::Clock,model::ModelSetup)
 
     # reinitalize progress meter, minus one to exclude first_timesteps! which contain compilation
-    (;enabled, showspeed, desc) = feedback.progress_meter
-    feedback.progress_meter = ProgressMeter.Progress(clock.n_timesteps-1;enabled,showspeed, desc)
+    (;showspeed, desc) = feedback.progress_meter
+    (;verbose) = feedback
+    feedback.progress_meter = ProgressMeter.Progress(clock.n_timesteps-1;enabled=verbose, showspeed, desc)
     
     # set to false to recheck for NaRs
     feedback.nars_detected = false