24 Basic concepts of numerics
A weather model is a copy of the atmosphere. In it, all fields and differential operators are replaced by discretized versions, since you can only store a finite amount of information. One writes a substitution for a differential operator $D$
\[ \begin{align} D\to D', \end{align} \]
analogous for atmospheric conditions
\[ \begin{align} Z\to Z'. \end{align} \]
Then you define a dynamic, i.e. a procedure
\[ \begin{align} Z'\left(t\right)\stackrel{D'}{\to}Z'\left(t + \Delta t\right) \end{align} \]
to determine the time evolution of the model atmosphere with the aim of bringing the state trajectory of the model as close as possible to that of the real atmosphere. To do this, you define a distance
\[ \begin{align} \left|Z - Z'\right|\geq 0. \end{align} \]
Weather forecast has the following workflow:
- Observation of the state of the atmosphere.
- Determination of the initial state of the model.
- Integration of the model.
Errors arise for the following reasons:
- The initial state was not fully observed.
- The initial state was incorrectly observed.
- The governing equations are equations of statistical physics and therefore do not apply exactly.
- The model can store less information than the real atmosphere.
- The dynamics of the model are different from those of the real atmosphere.
- There are rounding errors due to the finite computer precision.
- The computer is making errors.
24.1 Basic concepts of discretizations
Let the solution of a differential equation be denoted by $F$ and the solution of a discretized version of the equation by $F_d$. A numerical procedure is called consistent if and only if the solution of the scheme converges to the solution of the differential equation as $\Delta x$ and $\Delta t$ become smaller, i.e. $\lim\limits_{\Delta\to 0}F_d = F$ in the above terms. A discretization is called self-consistent if all conserved quantities of the continuous system are also conserved quantities of the discretized system. This can also be transferred to approximations. A scheme is called stable if and only if the calculated solution is bounded for all times, i.e. if a $C>0$ exists with $|F_d|
24.2 Finite differences
Let an infinitely often differentiable function $f:\mathbb{R}\to\mathbb{R}$ be given. You can write using a Taylor series
\[ \begin{align} f\left(x + \Delta x\right) = f\left(x\right) + \Delta xf'\left(x\right) + \frac{f''\left(x\right)}{2!}\Delta x^2 + \frac{f'''\left(x\right)}{3!}\Delta x^3 + \dotsc\tag{24.5}\label{eq:taylor_links} \end{align} \]
and
\[ \begin{align} f\left(x - \Delta x\right) = f\left(x\right) - \Delta xf'\left(x\right) + \frac{f''\left(x\right)}{2!}\Delta x^2 - \frac{f'''\left(x\right)}{3!}\Delta x^3 + \dotsc\tag{24.6}\label{eq:taylor_rechts} \end{align} \]
If you subtract both equations, you get
\[ \begin{align} f\left(x + \Delta x\right) - f\left(x - \Delta x\right) = 2\Delta xf'\left(x\right) + R', \end{align} \]
where the remainder $R'$ represents the third and higher order terms. This results in a change
\[ \begin{align} f'\left(x\right) = \frac{f\left(x + \Delta x\right) - f\left(x - \Delta x\right)}{2\Delta x} + R, \tag{24.8}\label{eq:approx_derivative_zentral} \end{align} \]
where $R$ is a polynomial in $\Delta x$ with the lowest power 2, one says that (24.8) gives a second-order approximation of the derivative and instead of $R$ writes the symbol $\mathcal{O}\left(\Delta x^2\right)$. This way of approximating derivatives is called central difference quotients. A one-sided difference quotient, for example
\[ \begin{align} f'\left(x\right) \approx\frac{f\left(x + \Delta x\right) - f\left(x\right)}{\Delta x} \end{align} \]
is also a suitable approximation, but this only converges to the first order, as can be easily done by rearranging Eq. (24.5) can be seen. The higher the order, the better, since the accuracy of a Taylor expansion increases with the order. For this reason, central spatial difference quotients are used wherever possible. You can also derive an approximation for the second derivative of a function from the equations (24.5) and (24.6):
\[ \begin{align} f\left(x + \Delta x\right) + f\left(x - \Delta x\right)&= 2f\left(x\right) + f''\left(x\right)\Delta x^2 + \mathcal{O}\left(\Delta x^4\right)\nonumber\\ \Leftrightarrow f''\left(\Delta x\right) &= \frac{f\left(x + \Delta x\right) - 2f\left(x\right) + f\left(x - \Delta x\right)}{\Delta x^2} + \mathcal{O}\left(\Delta x^2\right) \end{align} \]
For example, if you were given a function $f\left(x\right) = A\sin\left(kx\right)$, then $f' = Ak\cos\left(kx\right)$, Eq. However, (24.8) gives
\[ \begin{align} f' &\approx& \frac{A}{2\Delta x}\left(\sin\left(kx + k\Delta x\right) - \sin\left(kx - k\Delta x\right)\right) = \frac{A}{\Delta x}\cos\left(kx\right)\sin\left(k\Delta x\right) = \frac{\sin\left(k\Delta x\right)}{k\Delta x}f'. \end{align} \]
This shows that the approximation, as expected, converges to the correct value of the derivative if $\Delta x$ approaches zero (we have $\lim\limits_{\alpha\to 0}\frac{\sin\left(\alpha\right)}{\alpha} = 1$ according to L'Hospital's rule). If you write $k = \frac{2\pi}{\lambda}$ with $\lambda$ as the wavelength and insert $\Delta x = \lambda/2$, the central difference quotient always results in $f' = 0$. However, sufficiently long waves ($\lambda\gg 2\Delta x$) are resolved well.
24.3 Plane spectral method
24.4 Spherical spectral method
Let $\psi = \psi\left(\varphi, \lambda, z\right)$ be a scalar field. Since the spherical surface functions presented in Sect. C.5 form a complete orthonormal basis on spherical shells, it is possible to express the horizontal dependence of $\psi$ by this set of functions:
\[ \begin{align} \psi\left(\varphi, \lambda, z\right) = \sum_{l = 0}^\infty\sum_{m = -l}^l\newtilde{\psi}_{l, m}\left(z\right)Y_{l, m}\left(\varphi, \lambda\right) \end{align} \]
The definition of the spherical area functions Eq. (C.155) is formulated here in geographical coordinates:
\[ \begin{align} Y_{l, m}\left(\varphi, \lambda\right) = \sqrt{\frac{2l + 1}{4\pi}\frac{\left(l - m\right)!}{\left(l + m\right)!}}P_{l, m}\left(\sin\left(\varphi\right)\right)\exp\left(im\lambda\right) \end{align} \]