The barotropic vorticity model describes the evolution of a 2D non-divergent flow with
velocity components \mathbf{u} = (u, v)
through self-advection, forces and dissipation.
Due to the non-divergent nature of the flow, it can be described by (the vertical component)
of the relative vorticity \zeta = \nabla \times \mathbf{u}
.
The dynamical core presented here to solve the barotropic vorticity equations largely follows the idealized models with spectral dynamics developed at the Geophysical Fluid Dynamics Laboratory1: A barotropic vorticity model2.
Many concepts of the [Shallow water model](@ref shallow_water_model) and the [Primitive equation model](@ref primitive_equation_model) are similar, such that for example [horizontal diffusion](@ref diffusion) and the [Time integration](@ref leapfrog) are only explained here.
The barotropic vorticity equation is the prognostic equation that describes the time evolution of
relative vorticity \zeta
with advection, Coriolis force, forcing and diffusion in a
single global layer on the sphere.
We denote timet
, velocity vector \mathbf{u} = (u, v)
, Coriolis parameter f
,
and hyperdiffusion (-1)^{n+1} \nu \nabla^{2n} \zeta
(n
is the hyperdiffusion order, see [Horizontal diffusion](@ref diffusion)).
We also include possible forcing terms
F_\zeta, \mathbf{F}_\mathbf{u} = (F_u, F_v)
which act on the vorticity and/or on the
zonal velocity u
and the meridional velocity v
and hence the curl
\nabla \times \mathbf{F}_\mathbf{u}
is a tendency for relative vorticity \zeta
.
See Extending SpeedyWeather how to define these.
Starting with some relative vorticity \zeta
, the Laplacian is
inverted to obtain the stream function \Psi
The zonal velocity u
and meridional velocity v
are then the (negative) meridional gradient
and zonal gradient of \Psi
which is described in Derivatives in spherical coordinates. Using u
and v
we can then
advect the absolute vorticity \zeta + f
. In order to avoid to calculate both the curl and the
divergence of a flux we rewrite the barotropic vorticity equation as
with \mathbf{u}_\perp = (v, -u)
the rotated velocity vector, because
-\nabla\cdot\mathbf{u} = \nabla \times \mathbf{u}_\perp
. This is the form that is solved
in the BarotropicModel
, as outlined in the following section.
We briefly outline the algorithm that SpeedyWeather.jl uses in order to integrate the barotropic vorticity equation. As an initial step
0. Start with initial conditions of \zeta_{lm}
in spectral space and
transform this model state to grid-point space:
- Invert the Laplacian of vorticity
\zeta_{lm}
to obtain the stream function\Psi_{lm}
in spectral space - obtain zonal velocity
(\cos(\theta)u)_{lm}
through a Meridional derivative - obtain meridional velocity
(\cos(\theta)v)_{lm}
through a Zonal derivative - Transform zonal and meridional velocity
(\cos(\theta)u)_{lm}
,(\cos(\theta)v)_{lm}
to grid-point space - Unscale the
\cos(\theta)
factor to obtainu, v
- Transform
\zeta_{lm}
to\zeta
in grid-point space
Now loop over
- Compute the forcing (or drag) terms
F_\zeta, \mathbf{F}_\mathbf{u}
- Multiply
u, v
with\zeta+f
in grid-point space - Add
A = F_u + v(\zeta + f)
andB = F_v - u(\zeta + f)
- Transform these vector components to spectral space
A_{lm}
,B_{lm}
- Compute the curl of
(A, B)_{lm}
in spectral space, add toF_\zeta
to accumulate the tendency of\zeta_{lm}
- Compute the [horizontal diffusion](@ref diffusion) based on that tendency
- Compute a leapfrog time step as described in [Time integration](@ref leapfrog) with a Robert-Asselin and Williams filter
- Transform the new spectral state of
\zeta_{lm}
to grid-pointu, v, \zeta
as described in 0. - Possibly do some output
- Repeat from 1.
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). 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.
In SpeedyWeather.jl the diffusion is applied implicitly. For that, consider a
leapfrog scheme with time step \Delta t
of
variable \zeta
to obtain from time steps i-1
and i
, the next time step i+1
with d\zeta
being some tendency evaluated from \zeta_i
. Now we want to add a diffusion term
(-1)^{n+1}\nu \nabla^{2n}\zeta
with coefficient \nu
, which however, is implicitly calculated from \zeta_{i+1}
, then
As the application of (-1)^{n+1}\nu\nabla^{2n}
is, for every spectral mode, equivalent to a multiplication of
a constant, we can rewrite this to
and expand the numerator to
Hence the diffusion can be applied implicitly by updating the tendency d\zeta
as
which only depends on \zeta_{i-1}
. Now let D_\text{explicit} = (-1)^{n+1}\nu\nabla^{2n}
be the explicit part and
D_\text{implicit} = 1 - (-1)^{n+1} 2\Delta t \nu\nabla^{2n}
the implicit part. Both parts can be precomputed and are
D_\text{implicit} = 1 - 2\Delta t \nu\nabla^{2n}
the implicit part. Both parts can be precomputed and are
only an element-wise multiplication in spectral space. For every spectral harmonic l, m
we do
Hence 2 multiplications and 1 subtraction with precomputed constants.
However, we will normalize the (hyper-)Laplacians as described in the following. This also will take care of
the alternating sign such that the diffusion operation is dissipative regardless the power n
.
In physics, the Laplace operator \nabla^2
is often used to represent diffusion due to viscosity in a fluid or
diffusion that needs to be added to retain numerical stability. In both cases,
the coefficient is \nu
of units \text{m}^2\text{s}^{-1}
and the full operator reads as \nu \nabla^2
with units (\text{m}^2\text{s}^{-1})(\text{m}^{-2}) = \text{s}^{-1}
.
This motivates us to normalize the Laplace operator by a constant of units \text{m}^{-2}
and the coefficient
by its inverse such that it becomes a damping timescale of unit \text{s}^{-1}
. Given the application in
spectral space we decide to normalize by the largest eigenvalue -l_\text{max}(l_\text{max}+1)
such that
all entries in the discrete spectral Laplace operator are in [0, 1]
. This also has the effect that the
alternating sign drops out, such that higher wavenumbers are always dampened and not amplified.
The normalized coefficient \nu^* = l_\text{max}(l_\text{max}+1)\nu
(always positive) is
therefore reinterpreted as the (inverse) time scale at which the highest wavenumber is dampened
to zero due to diffusion. Together we have
and the hyper-Laplacian of power n
follows as
and the implicit part is accordingly D^\text{implicit, n}_{l, m} = 1 - 2\Delta t D^\text{explicit, n}_{l, m}
.
Note that the diffusion time scale \nu^*
is then also scaled by the radius, see next section.
Similar to a non-dimensionalization of the equations, SpeedyWeather.jl scales the barotropic vorticity
equation with R^2
to obtain normalized gradient operators as follows.
A scaling for vorticity \zeta
and stream function \Psi
is used that is
This is also convenient as vorticity is often 10^{-5}\text{ s}^{-1}
in the atmosphere,
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
where the dimensionless gradients simply omit the scaling with 1/R
, \tilde{\nabla} = R\nabla
.
The Barotropic vorticity equation scaled with R^2
is
with
\tilde{t} = tR^{-1}
, the scaled timet
\mathbf{u} = (u, v)
, the velocity vector (no scaling applied)\tilde{f} = fR
, the scaled Coriolis parameterf
\tilde{\mathbf{F}} = R\mathbf{F}
, the scaled forcing vector\mathbf{F}
\tilde{\nu} = \nu^* R
, the scaled diffusion coefficient\nu^*
, which itself is normalized to a damping time scale, see Normalization of diffusion.
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 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 the user. In hopefully many other cases it is clearly denoted that a variable or constant is scaled.
SpeedyWeather.jl is based on the Leapfrog time integration,
which, for relative vorticity \zeta
, is
in its simplest form
meaning we step from the previous time step i-1
, leapfrogging over the current time stepi
to the next time step i+1
by evaluating the tendencies on the right-hand side RHS
at the current time step i
. The time stepping is done in spectral space.
Once the right-hand side RHS
is evaluated, leapfrogging is a linear operation, meaning
that its simply applied to every spectral coefficient \zeta_{lm}
as one would evaluate
it on every grid point in grid-point models.
For the Leapfrog time integration two time steps of the prognostic variables have to be stored,
i-1
and i
. Time step i
is used to evaluate the tendencies which are then added
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 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.
- an Euler forward step with
\Delta t/2
, then - one leapfrog time step with
\Delta t
, then - leapfrog with
2 \Delta t
till the end
This is particularly done in a way that after 2. we have t=0
at i-1
and t=\Delta t
at i
available so that 3. can start the leapfrogging without any offset from the intuitive spacing
0, \Delta t, 2\Delta t, 3\Delta t, ...
. The following schematic can be useful
time at step i-1 |
time at step i |
time step at i+1 |
|
---|---|---|---|
Initial conditions | t = 0 |
||
1: Euler | (T) \quad t = 0 |
t=\Delta t/2 |
|
2: Leapfrog with \Delta t |
t = 0 |
(T) \quad t = \Delta t/2 |
t = \Delta t |
3 to n : Leapfrog with 2\Delta t |
t-\Delta t |
(T) \qquad \quad \quad t |
t+\Delta t |
The time step that is used to evaluate the tendencies is denoted with (T). It is always the time step furthest in time that is available.
The standard leapfrog time integration is often combined with a Robert-Asselin filter34
to dampen a computational mode. The idea is to start with a standard leapfrog step to obtain
the next time step i+1
but then to correct the current time step i
by applying a filter
which dampens the computational mode. The filter looks like a discrete Laplacian in time
with a (1, -2, 1)
stencil, and so, maybe unsurprisingly, is efficient to filter out
a "grid-scale oscillation" in time, aka the computational mode. Let v
be the unfiltered
variable and u
be the filtered variable, F
the right-hand side tendency,
then the standard leapfrog step is
Meaning we start with a filtered variable u
at the previous time step i-1
, evaluate
the tendency F(v_i)
based on the current time step i
to obtain an
unfiltered next time step v_{i+1}
. We then filter the current time step i
(which will become i-1
on the next iteration)
by adding a discrete Laplacian with coefficient \tfrac{\nu}{2}
to it, evaluated
from the available filtered and unfiltered time steps centred around i
:
v_{i-1}
is not available anymore because it was overwritten by the filtering
at the previous iteration, u_i, u_{i+1}
are not filtered yet when applying
the Laplacian. The filter parameter \nu
is typically chosen between 0.01-0.2,
with stronger filtering for higher values.
Williams5 then proposed an additional filter step to regain accuracy
that is otherwise lost with a strong Robert-Asselin filter67.
Now let w
be unfiltered, v
be once filtered, and u
twice filtered, then
with the Williams filter parameter \alpha \in [0.5, 1]
. For \alpha=1
we're back with the Robert-Asselin filter (the first two lines).
The Laplacian in the parentheses is often called a displacement,
meaning that the filtered value is displaced (or corrected) in the direction
of the two surrounding time steps. The Williams filter now also applies
the same displacement, but in the opposite direction to the next time
step i+1
as a correction step (line 3 above) for a once-filtered
value v_{i+1}
which will then be twice-filtered by the Robert-Asselin
filter on the next iteration. For more details see the referenced publications.
The initial Euler step (see [Time integration](@ref leapfrog), Table) is not filtered. Both the the Robert-Asselin and Williams filter are then switched on for all following leapfrog time steps.
Footnotes
-
Geophysical Fluid Dynamics Laboratory, Idealized models with spectral dynamics ↩
-
Geophysical Fluid Dynamics Laboratory, The barotropic vorticity equation. ↩
-
Robert, André. "The Integration of a Low Order Spectral Form of the Primitive Meteorological Equations." Journal of the Meteorological Society of Japan 44 (1966): 237-245. ↩
-
ASSELIN, R., 1972: Frequency Filter for Time Integrations. Mon. Wea. Rev., 100, 487-490, doi:10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2 ↩
-
Williams, P. D., 2009: A Proposed Modification to the Robert-Asselin Time Filter. Mon. Wea. Rev., 137, 2538-2546, 10.1175/2009MWR2724.1. ↩
-
Amezcua, J., E. Kalnay, and P. D. Williams, 2011: The Effects of the RAW Filter on the Climatology and Forecast Skill of the SPEEDY Model. Mon. Wea. Rev., 139, 608-619, doi:10.1175/2010MWR3530.1. ↩
-
Williams, P. D., 2011: The RAW Filter: An Improvement to the Robert-Asselin Filter in Semi-Implicit Integrations. Mon. Wea. Rev., 139, 1996-2007, doi:10.1175/2010MWR3601.1. ↩