The Use of Semigeostrophic Theory to Diagnose the Behaviour of an Atmospheric GCM

Abstract

A diagnostic method is presented for analysing the large-scale behaviour of the Met Office Unified Model, which is a comprehensive atmospheric model used for weather and climate prediction. Outside the boundary layer, on scales larger than the radius of deformation, semi-geostrophic theory will give an accurate approximation to the model evolution. In particular, the ageostrophic circulation required to maintain geostrophic and hydrostatic balance against prescribed forcing and a rate of change of the geostrophic pressure can be calculated. In the tropics, the balance condition degenerates to the weak temperature gradient approximation. Within the boundary layer, the semi-geotriptic approximation has to be used because friction and rotation are equally important. Assuming the calculated pressure tendency and ageotriptic circulation match the observed model behaviour, the influence of the large-scale state and the nature of the forcing on the model response can be deduced in a straightforward way. The capabilities of the diagnostic are illustrated by comparing predictions of the ageotriptic circulation from the theory and the model. It is then used to show that the effects of latent heat release can be included by modifying the static stability, and to show the effect of an idealised tropical heat source on the subtropical jet. Finally, the response of the ageotriptic flow to boundary layer heating in the tropics is demonstrated. These illustrations show that the model behaviour on large scales conforms with theoretical expectations, so that the results of the diagnostic can be used to aid the development of further improvements to the model, in particular investigating systematic errors and understanding the large-scale atmospheric response to forcing.

1. Introduction

In common with most other operational weather services and climate research centres, the Met Office uses a comprehensive atmospheric modelling system which is coupled to models of other parts of the Earth system, called the Unified Model (UM) [1]. In common with other centres, the Office undertakes a research programme aimed at improving the performance of the system. A particular focus is the study of systematic errors, where the Office uses a “seamless” approach which exploits the fact that many such errors are present at different model resolutions and time periods [2]. Typically, the errors have large spatial scales. However, correcting them is not straightforward because of the strong coupling between atmospheric variables on large scales resulting from the requirements of geostrophic and hydrostatic balance. Near the equator, the geostrophic requirement degenerates to the “weak temperature gradient” approximation (WTG) [3]. In the boundary layer, friction is also important, leading to “geotriptic” balance [4], where the horizontal pressure gradient is balanced by friction as well as the Coriolis term. This is particularly important in the tropics because it can support horizontal pressure and temperature gradients. The ageotriptic flow is strongly constrained by the need to maintain geotriptic and hydrostatic balance when the model is integrated in time. his constraint would be strictly enforced if the semi-geotriptic (SGT) equations were solved instead of the full governing equations. These generalise the semi-geostrophic (SG) equations by including friction in the lowest order balance. The result is that model changes which appear to be well-suited to reducing systematic errors may not have the intended effect.

In this paper, a diagnostic technique is proposed and illustrated which allows the geotriptic evolution and the ageotriptic circulation of the UM to be calculated for given data. The use of a decomposition of the UM fields into geotriptic and ageotriptic parts is only justified on scales larger than the deformation radius, where the potential vorticity can be well approximated by a function of pressure. In this case, the pressure can be used as the controlling scalar variable and geotriptic winds and hydrostatic temperatures calculated from it. In three-dimensional flow, the requirement that the horizontal scale is greater than the deformation radius becomes a requirement that the aspect ratio is less than the ratio of the Coriolis and Brunt–Väisälä frequencies. The Lagrangian Rossby number also has to be small, restricting the curvature of trajectories. This latter condition is satisfied by fronts and jet streams which are characterised by nearly straight trajectories. However, such flows are only stable to three-dimensional perturbations if the aspect ratio condition is satisfied (see [5], p. 214 and [6]). While these conditions appear very restrictive, they are also the conditions under which coherent anomalies are maintained. Near the equator, the decomposition is only justified outside the boundary layer for flows that are almost zonally symmetric and the WTG holds. These issues are discussed fully in [5]. Within the boundary layer, the usefulness of this method is greatly enhanced because horizontal pressure and temperature gradients can be supported. For instance, it was shown in [7] that using a SGT rather than a SG model allowed large-scale idealised cross-equatorial flows to be simulated.

The diagnostic works by solving the SGT system for a single timestep using UM pressure data and forcing fields from the physical parameterisations. Assuming that the UM solutions can indeed be well approximated by computed SGT solutions on large scales, systematic errors in the UM should be visible in the resulting large scale pressure tendencies and ageotriptic winds. Since the SGT evolution over a timestep is linear in the forcing fields, and is representative of the longer time evolution, it should be possible to relate the systematic errors to particular features of the forcing. More generally, the diagnostic can be used to study the large-scale response of the atmosphere to forcing by calculating the effect of changes to the forcing fields used or to changes in the basic state pressure field. This will aid the study of issues such as tropical–extratropical interaction, and its dependence on the large-scale flow.

The semi-geostrophic (SG) model was originally introduced by Eliassen [8] as an alternative description of large-scale flows to the quasi-geostrophic model of Ref. [9]. It is based on what was subsequently called the Type II geostrophic scaling by Phillips [10] and in particular allows full variation of the static stability and Coriolis parameter. It is reviewed with a broader perspective by Eliassen [11]. The model was then reintroduced by Hoskins [12] to study frontogenesis. In that study, a constant Coriolis parameter is used, which allows the equations to be transformed to a simpler form by using the geostrophic coordinate transformation. This version of the SG model can be solved for large times in isentropic coordinates by using optimal transport methods, as reviewed by Cullen [5]. Once the variable Coriolis parameter is included, the equations have to be solved in physical space. Formal arguments demonstrating how the equations can then be solved were given by Cullen et al. [13,14]. However, the only rigorous result so far is limited to short times, as given by Cheng et al. [15]. This result is sufficient to justify the diagnostic procedure used in this paper, where the equations are solved for a pressure tendency and ageostrophic wind at a single time in the manner of Ref. [16]. The extension to the semi-geotriptic (SGT) model was described by Beare et al. [4]. This is very important because it allows a realistic lower boundary condition. The need to use a nonlinear formulation of vertical momentum diffusion has meant that no rigorous treatment of this set of equations is yet available.

In the UM, the large-scale fields will be close to geotriptic and hydrostatic, so the ageotriptic circulation required to maintain this balance against forcing will be predicted by SGT theory. Thus, if the UM forcing is diagnosed and input into the calculation of a single SGT timestep, the resulting ageotriptic flow should form a part of the ageotriptic flow actually produced by the UM. The total UM circulation will also contain transients, not subject to this constraint. In studying systematic UM errors, it is useful to determine what aspects of the circulation are constrained by the balance requirement, and which parts are evolving independently.

This calculation is a three-dimensional generalisation of the widely used Sawyer–Eliassen equation (SEE) (e.g., [17]). An important feature of SEE is that the calculation of the ageostrophic circulation depends strongly on the static and inertial stability of the large-scale state. The effect of moisture on static stability can be included. Most of the applications of the SEE are in a vertical cross-section, where it is easy to justify its use. Examples specifically exploiting SG theory are found in [18,19]. Three dimensional solutions using SG were given by Wernli et al. [20], where the results are compared with solutions of the full equations using idealised data, and Ragone et al. [21], where the results are compared with other idealised models. Examples including the boundary layer using SGT theory were given by Beare et al. [4]. In three-dimensions, the validity of the SEE requires SG theory to be accurate, which implies the restrictions discussed above. This has limited its use.

In the next section, the analytical formulation of the diagnostic is presented. The method works from the pressure field, which implies geotriptic winds and hydrostatic potential temperatures. As noted above, this is appropriate on scales where the SGT approximation is valid. However, UM data will contain non-geotriptic motions which affect the pressure, so steps have to be taken to minimise their effect on the diagnostic calculation. The equation for the pressure tendency is an elliptic equation, and the solution procedure is described. The discretisation of the equation requires care over orography, and it is essential to follow the principles of the UM scheme.

Since systematic errors in the UM will typically have a large scale signal, it should be possible to reproduce them by applying the diagnostic to a large number of cases and averaging the results. Before doing so, however, it is necessary to demonstrate the usefulness of the diagnostic in individual cases. The results section therefore includes four examples illustrating the performance of the diagnostic. The first is a comparison of the actual ageotriptic circulation produced by the UM in the extratropics with that calculated using SGT theory. This shows that the procedure is viable with full UM data subject to the extra processing noted in the previous paragraph. The second example calculates the vertical motion required to maintain large scale balance in cloudy air assuming moist static stability, and then inferring the latent heat release from the difference between the potential temperature increment calculated using this vertical motion with moist and dry static stability. This should represent the “forced” part of the latent heating actually calculated by the UM physics. The third example shows the effect of an idealised tropical heat source on the extratropical flow following [22]. There is geostrophic adjustment in the north–south direction, as illustrated by Schubert et al. [23], which should be realistic. The diagnostic also adjusts in the east-west direction, but except for very long forcing timescales this would in reality generate a tropical wave, as in the simulations of Jin, et al. [22] using an idealised general circulation model. The final example shows that thermally forced ageotroptic circulations in the boundary layer can be simulated.

2. Materials and Methods

2.1. The SG Approximation to the UM Equations

The diagnostic procedure uses the SGT and hydrostatic approximation in the deep atmosphere compressible equations solved by the UM [24]. No other approximations are made, so that the diagnostic should extract the geotriptic part of the UM evolution.

We first write the equations in Lagrangian form. The equations are written in spherical polar and terrain-following coordinates $(\lambda ,\varphi ,\eta )$, where $\lambda ,\varphi$ are longitude and latitude, and $\eta$ is the model level number. The true radial coordinate is r. The boundary layer momentum mixing is included, with stability dependent diffusion coefficient $K_m$ as calculated by the UM. A derivative with respect to r is interpreted as a derivative with respect to $\eta$ multipled by $\partial \eta /\partial r$ with $(\lambda ,\varphi )$ constant.

$$\begin{array}{c}\hfill \frac{D\mathbf{u}}{Dt}+c_{pd}{\theta }_v\nabla \pi +2\mathrm{\Omega }\times \mathbf{u}=\mathbf{g}+\frac{\partial }{\partial r}\left(K_m\frac{\partial {\mathbf{u}}_h}{\partial r}\right)+S_{\mathbf{u}},\\ \hfill \frac{D{\theta }_{vd}}{Dt}=S_{{\theta }_{vd}},\\ \hfill \frac{D\mu }{Dt}=S_{\mu },\\ \hfill \frac{\partial {\rho }_d}{\partial t}+\nabla ·\left({\rho }_d\mathbf{u}\right)=0.\end{array}$$

(1)

where

$$\begin{array}{c}\hfill {\rho }_d=\left(\frac{p_0}{R_d}\right)\frac{{\pi }^{\frac{1-{\kappa }_d}{{\kappa }_d}}}{{\theta }_{vd}},\\ \hfill R_d={\kappa }_d{c}_{pd},\\ \hfill \mu =\mathrm{\Sigma }{m}_X,X=v,cl,cf\cdots ,\\ \hfill {\theta }_v=\frac{{\theta }_{vd}}{1+\mu }.\end{array}$$

(2)

In these equations, $\mathbf{u}$ is the vector velocity field and ${\theta }_v$ the virtual potential temperature with ${\theta }_{vd}$ the virtual potential temperature for dry air. ${\rho }_d$ is the density of dry air and $\pi$ the Exner pressure ${\left(\frac{p}{p_0}\right)}^{{\kappa }_d}$. $p_0$ is a reference pressure, $R_d$ is the gas constant for dry air and $c_{pd}$ is the specific heat of dry air at constant pressure. ${\kappa }_d$ is defined by the second part of Equation (2). $m_X$ indicates the mixing ratio of various moisture constituents denoted by X and $\mu$ is the sum of all these mixing ratios. $\mathrm{\Omega }$ is the Earth’s rotation vector and $\mathbf{g}$ is the gravitational and centrifugal vector. $S_x$ is a source term for variable x. The suffix h denotes horizontal components.

The lower boundary condition is no slip. The upper boundary condition is a rigid lid.

The semi-geotriptic approximation in spherical geometry assumes first that the shallow atmosphere hydrostatic approximation holds [25], so that the rotation term $2\mathrm{\Omega }\times \mathbf{u}$ is replaced by its horizontal components $(-fv,fu)$ where $f=2\mathrm{\Omega }\sin \varphi$, $\varphi$ is the latitude, and that the gravitational vector $\mathbf{g}$ only has a component in the local vertical. Then, the vertical momentum in Equation (1) becomes

$$c_{pd}{\theta }_v\frac{\partial \pi }{\partial r}=-g.$$

(3)

To understand the solutions, we subtract a reference state ${\pi }_0\left(r\right)$ satisfying

$$c_{pd}{\theta }_0\frac{\partial {\pi }_0}{\partial r}=-g,$$

(4)

where ${\theta }_0$ is a constant, from $\pi$, giving

$$c_{pd}{\theta }_v\frac{\partial {\pi }^{\prime }}{\partial r}=g\frac{{\theta }^{\prime }}{{\theta }_0},$$

(5)

where ${\theta }^{\prime }={\theta }_v-{\theta }_0$. The Brunt–Väisälä frequency N is then given by

$$N^2=\frac{g}{{\theta }_0}\frac{\partial {\theta }^{\prime }}{\partial r}.$$

(6)

Assume typical scalings P and $\mathrm{\Theta }$ for ${\pi }^{\prime }$ and ${\theta }^{\prime }$, ignore $\mu$ and assume typical length scales L and H in the horizontal and vertical. Then, $\partial {\pi }^{\prime }/\partial r\simeq P/H$ and the horizontal pressure gradient term ${\nabla }_h\varphi \simeq P/L$. In some regimes, the ratio of horizontal and vertical gradients of ${\pi }^{\prime }$ may be much less than $H/L$ if ${\pi }^{\prime }$ is close to hydrostatic balance with a spatially uniform stably stratified reference profile $\theta \left(r\right)$. This is not generally the case on large scales with real data, particularly because of large horizontal variations of the tropopause height.

The geostrophic regime is characterised by the terms $c_{pd}{\theta }_v{\nabla }_h{\pi }^{\prime }$ in the horizontal momentum equations having the same magnitude as $(-fv,fu)$. Equations (5) and (6) show that $c_{pd}{\theta }_v\partial {\pi }^{\prime }/\partial r\simeq N^{2}H$. Defining U as a horizontal velocity scale, the geostrophic regime requires that $fU\simeq N^2{H}^2/L$. Now, define the dimensionless Rossby and Froude numbers, $Ro$ and $Fr$, respectively, by

$$Ro=\frac{U}{fL},Fr=\frac{U}{NH}.$$

(7)

Then, geostrophic balance requires $Ro=F{r}^2$. The semi-geostrophic regime is thus characterised by a small parameter $ϵ=Ro=F{r}^2$. The radius of deformation $L_D$ is defined as $NH/f$, so that $L/L_D=\mathrm{O}\left({ϵ}^{-\frac{1}{2}}\right).$

Geotriptic balance requires that the horizontal pressure gradient is comparable to the friction term. If K is a typical magnitude of the momentum diffusion coefficient $K_m$, this requires that $N^2{H}^2/L\simeq KU/h^2$ where $h\ll H$ is a typical boundary layer depth. This scaling is discussed in [26].

Since the aim is to apply the diagnostic procedure directly to model data, the subtraction of a basic state as in Equation (4) is not used. Therefore, we define geotriptic and hydrostatic balance and the geotriptic wind ${\mathbf{u}}_e=(u_e,v_e)$ consistently with Equation (1) by

$$\begin{array}{c}{c}_{pd}{\theta }_v{\nabla }_h\pi -(f{v}_e,-f{u}_e)=\frac{\partial }{\partial r}\left(K_m\frac{\partial {\mathbf{u}}_e}{\partial r}\right),\hfill \\ c_{pd}{\theta }_v\frac{\partial \pi }{\partial r}=-g.\hfill \end{array}$$

(8)

Then, the first part of Equation (1) can be approximated to O(${ϵ}^2$) outside the boundary layer by

$$\begin{array}{c}\hfill \frac{D{\mathbf{u}}_e}{Dt}-f(v-v_e,u_e-u)=\frac{\partial }{\partial r}\left(K_m\frac{\partial ({\mathbf{u}}_e-{\mathbf{u}}_h)}{\partial r}\right)+S_{\mathbf{u}},\\ \hfill c_{pd}{\theta }_v\frac{\partial \pi }{\partial r}=-g.\end{array}$$

(9)

Note that the full velocity $\mathbf{u}$ is used in the Lagrangian time derivative $D/Dt$ as in the geostrophic momentum approximation. Equation (9) is only accurate to O($ϵ$) in the boundary layer, assuming $h/H$ is independent of $ϵ$. As discussed in [26], it is more typical for $h/H$ to decrease with $ϵ$; that paper suggests $h/H\simeq {ϵ}^{0.7}$. The SGT approximation to Equation (1) consists of Equations (8), (9) and Equation (1) after the first. The boundary conditions are the same as for Equation (1).

2.2. The Diagnostic Equations

The diagnostic equations are based on the original Equation (1). The equations are written in Eulerian form as in [16] to derive an Eulerian pressure tendency and an equation for the ageotriptic winds. The principle of the derivation is the same as that for the basic two-dimensional SEE (e.g., [27], p. 287). The first part of Equation (8) is written

$$\begin{array}{c}f{v}_e+\frac{\partial }{\partial r}\left(K_m\frac{\partial u_e}{\partial r}\right)=\frac{c_{pd}{\theta }_v}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right),\hfill \\ -f{u}_e+\frac{\partial }{\partial r}\left(K_m\frac{\partial v_e}{\partial r}\right)=\frac{c_{pd}{\theta }_v}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right).\hfill \end{array}$$

(10)

where the horizontal derivatives are computed on $\eta$ surfaces.

Equation (10) can be used to generate a single equation for $u_e$ (eliminating $v_e$ terms):

$$\begin{array}{c}\hfill f^2{u}_e+\frac{\partial }{\partial r}\left(K_m\frac{{\partial }^2}{\partial r^2}\left(K_m\frac{\partial u_e}{\partial r}\right)\right)=\\ \hfill \frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\left(\frac{c_{pd}{\theta }_v}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)\right)\right)-f\frac{c_{pd}{\theta }_v}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right).\end{array}$$

(11)

Given values of $\pi$, Equation (11) can be used to calculate $u_e$. The lower boundary condition is that $u_e=v_e=0$ at $\eta =0$. It is assumed that $K_m$ is zero above some $\eta$ level below the upper boundary, so that setting $K_m=0$ there will generate an upper boundary condition for Equation (11). For larger values of $\eta$, Equation (11) can be solved independently at each level.

The first part of Equation (9) and the second and third parts of Equation (1) can then be rewritten in Eulerian form as

$$\begin{array}{c}\hfill \frac{\partial u_e}{\partial t}+\mathbf{u}·\nabla u_e-f(v-v_e)=\frac{\partial }{\partial r}\left(K_m\frac{\partial (u_e-u)}{\partial r}\right)+S_u,\\ \hfill \frac{\partial v_e}{\partial t}+\mathbf{u}·\nabla v_e+f(u-u_e)=\frac{\partial }{\partial r}\left(K_m\frac{\partial (v_e-v)}{\partial r}\right)+S_v,\\ \hfill \frac{\partial {\theta }_v}{\partial t}+\mathbf{u}·\nabla {\theta }_v=\frac{1}{1+\mu }{S}_{{\theta }_{vd}}-\frac{{\theta }_{vd}}{{(1+\mu )}^2}{S}_{\mu }\equiv S_{{\theta }_v}.\end{array}$$

(12)

Equation (12) can be rewritten as

$$\mathbf{Q}\left(\begin{array}{c}{\displaystyle u-u_e}\\ {\displaystyle v-v_e}\\ {\displaystyle w}\end{array}\right)+\frac{\partial }{\partial t}\left(\begin{array}{c}{\displaystyle v_e}\\ {\displaystyle -u_e}\\ {\displaystyle {\theta }_v}\end{array}\right)=\mathbf{H},$$

(13)

where

$$\mathbf{Q}=\left(\begin{array}{ccc}{\displaystyle f+\frac{1}{r\cos \varphi }\frac{\partial v_e}{\partial \lambda }+\frac{u_e\tan \varphi }{r}}& {\displaystyle \frac{1}{r}\frac{\partial v_e}{\partial \varphi }+\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)}& {\displaystyle \frac{\partial v_e}{\partial r}}\\ {\displaystyle -\frac{1}{r\cos \varphi }\frac{\partial u_e}{\partial \lambda }+\frac{v_e\tan \varphi }{r}-\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)}& {\displaystyle f-\frac{1}{r}\frac{\partial u_e}{\partial \varphi }}& {\displaystyle -\frac{\partial u_e}{\partial r}}\\ {\displaystyle \frac{1}{r\cos \varphi }\frac{\partial {\theta }_v}{\partial \lambda }}& {\displaystyle \frac{1}{r}\frac{\partial {\theta }_v}{\partial \varphi }}& {\displaystyle \frac{\partial {\theta }_v}{\partial r}}\end{array}\right)$$

(14)

and

$$\mathbf{H}=\left(\begin{array}{c}{\displaystyle -{\mathbf{u}}_e·\nabla v_e+S_v}\\ {\displaystyle {\mathbf{u}}_e·\nabla u_e-S_u}\\ {\displaystyle -{\mathbf{u}}_e·\nabla {\theta }_v+S_{{\theta }_v}}\end{array}\right).$$

(15)

We now derive a single equation for the pressure tendency. Differentiating geotriptic balance (the second part of Equation (10)) with respect to time, and assuming $K_m$ is time independent, gives

$$-f\frac{\partial u_e}{\partial t}+\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\left(\frac{\partial v_e}{\partial t}\right)\right)=\frac{c_{pd}{\theta }_v}{r}\frac{\partial }{\partial t}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)+\frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)c_{pd}\frac{\partial {\theta }_v}{\partial t}.$$

(16)

Next, define

$$\mathbf{B}=\left(\begin{array}{ccc}f& -\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)& 0\\ \frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)& f& 0\\ 0& 0& g/{\theta }_v\end{array}\right).$$

(17)

Then, using Equation (10), the similar equation for $\partial \pi /\partial \lambda$ and the relation

$$c_{pd}{\theta }_v\frac{\partial }{\partial t}\left(\frac{\partial \pi }{\partial r}\right)=-\frac{g}{{\theta }_v}\frac{\partial {\theta }_v}{\partial t},$$

(18)

obtained by differentiating hydrostatic balance (the second part of Equation (9)) with respect to time, we can write

$$\mathbf{B}\left(\begin{array}{c}\frac{\partial v_e}{\partial t}\\ -\frac{\partial u_e}{\partial t}\\ -\frac{\partial {\theta }_v}{\partial t}\end{array}\right)=c_{pd}{\theta }_v\frac{\partial }{\partial t}\left(\begin{array}{c}{\displaystyle \frac{1}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)}\\ {\displaystyle \frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)}\\ {\displaystyle \frac{\partial \pi }{\partial r}}\end{array}\right)+\left(\begin{array}{c}{\displaystyle \frac{1}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)}\\ {\displaystyle \frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)}\\ {\displaystyle 0}\end{array}\right)c_{pd}\frac{\partial {\theta }_v}{\partial t}.$$

(19)

Using Equation (10), the second term on the right hand side can be replaced by

$$\begin{array}{c}\hfill \frac{1}{{\theta }_v}\mathbf{B}\left(\begin{array}{c}{\displaystyle v_e}\\ {\displaystyle -u_e}\\ {\displaystyle 0}\end{array}\right)\frac{\partial {\theta }_v}{\partial t}.\end{array}$$

Outside the boundary layer, this can be rewritten as

$$\mathbf{B}\left(\begin{array}{c}{\displaystyle \frac{v_e}{{\theta }_v}\frac{\partial {\theta }_v}{\partial t}}\\ {\displaystyle -\frac{u_e}{{\theta }_v}\frac{\partial {\theta }_v}{\partial t}}\\ {\displaystyle 0}\end{array}\right).$$

and so this term can be amalgamated with the left hand side giving

$$\mathbf{B}\left(\begin{array}{c}{\displaystyle \frac{\partial v_e}{\partial t}-\frac{v_e}{{\theta }_v}\frac{\partial {\theta }_v}{\partial t}}\\ {\displaystyle -\frac{\partial u_e}{\partial t}+\frac{u_e}{{\theta }_v}\frac{\partial {\theta }_v}{\partial t}}\\ {\displaystyle -\frac{\partial {\theta }_v}{\partial t}}\end{array}\right)\mathrm{or}\mathbf{B}\left(\begin{array}{c}{\displaystyle {\theta }_v\frac{\partial }{\partial t}\left(\frac{v_e}{{\theta }_v}\right)}\\ {\displaystyle -{\theta }_v\frac{\partial }{\partial t}\left(\frac{u_e}{{\theta }_v}\right)}\\ {\displaystyle -\frac{\partial {\theta }_v}{\partial t}}\end{array}\right).$$

Within the boundary layer, we note that the vertical variation of ${\mathbf{u}}_e/{\theta }_v$ will be dominated by the vertical variation of ${\mathbf{u}}_e$, so we use the amalgamated form of the left hand side of Equation (19) in the boundary layer also. This additional approximation is needed to allow the reduction to a single equation to proceed.

Using the third part of Equation (12), we can then rewrite Equations (13)–(15) as

$$\mathbf{B}{\mathbf{Q}}^{\prime }\left(\begin{array}{c}{\displaystyle u-u_e}\\ {\displaystyle v-v_e}\\ {\displaystyle w}\end{array}\right)+c_{pd}{\theta }_v\frac{\partial }{\partial t}\left(\begin{array}{c}{\displaystyle \frac{1}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)}\\ {\displaystyle \frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)}\\ {\displaystyle \frac{\partial \pi }{\partial r}}\end{array}\right)=\mathbf{B}{\mathbf{H}}^{\prime }$$

(20)

where

$${\mathbf{Q}}^{\prime }=\left(\begin{array}{ccc}{\displaystyle f+\frac{{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\left(\frac{v_e}{{\theta }_v}\right)+\frac{u_e\tan \varphi }{r}}& {\displaystyle \frac{{\theta }_v}{r}\frac{\partial }{\partial \varphi }\left(\frac{v_e}{{\theta }_v}\right)+\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)}& {\displaystyle {\theta }_v\frac{\partial }{\partial r}\left(\frac{v_e}{{\theta }_v}\right)}\\ {\displaystyle -\frac{{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\left(\frac{u_e}{{\theta }_v}\right)+\frac{v_e\tan \varphi }{r}-\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\right)}& {\displaystyle f-\frac{{\theta }_v}{r}\frac{\partial }{\partial \varphi }\left(\frac{u_e}{{\theta }_v}\right)}& {\displaystyle -{\theta }_v\frac{\partial }{\partial r}\left(\frac{u_e}{{\theta }_v}\right)}\\ {\displaystyle \frac{1}{r\cos \varphi }\frac{\partial {\theta }_v}{\partial \lambda }}& {\displaystyle \frac{1}{r}\frac{\partial {\theta }_v}{\partial \varphi }}& {\displaystyle \frac{\partial {\theta }_v}{\partial r}}\end{array}\right)$$

(21)

and

$${\mathbf{H}}^{\prime }=\left(\begin{array}{c}{\displaystyle -\left({\mathbf{u}}_e{\theta }_v\right)·\nabla \left(\frac{v_e}{{\theta }_v}\right)+S_v-\frac{v_e{S}_{{\theta }_v}}{{\theta }_v}}\\ {\displaystyle \left({\mathbf{u}}_e{\theta }_v\right)·\nabla \left(\frac{u_e}{{\theta }_v}\right)-S_u+\frac{u_e{S}_{{\theta }_v}}{{\theta }_v}}\\ {\displaystyle -{\mathbf{u}}_e·\nabla {\theta }_v+S_{{\theta }_v}}\end{array}\right).$$

(22)

We can then write ${\mathbf{BH}}^{\prime }$ as $\mathbf{G}$ where

$$\mathbf{G}=\left(\begin{array}{c}{\displaystyle -f\frac{u_e{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\frac{v_e}{{\theta }_v}-\frac{f{u}_e^2\tan \varphi }{r}-f\frac{v_e{\theta }_v}{r}\frac{\partial }{\partial \varphi }\frac{v_e}{{\theta }_v}+f{S}_v-\frac{f{v}_e{S}_{{\theta }_v}}{{\theta }_v}-\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\left(\frac{u_e{S}_{{\theta }_v}}{{\theta }_v}\right)\right)}\\ {\displaystyle f\frac{u_e{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\frac{u_e}{{\theta }_v}-\frac{f{u}_e{v}_e\tan \varphi }{r}+f\frac{v_e{\theta }_v}{r}\frac{\partial }{\partial \varphi }\frac{u_e}{{\theta }_v}-f{S}_u+\frac{f{u}_e{S}_{{\theta }_v}}{{\theta }_v}-\frac{\partial }{\partial r}\left(K_m\frac{\partial }{\partial r}\left(\frac{v_e{S}_{{\theta }_v}}{{\theta }_v}\right)\right)}\\ {\displaystyle -\frac{g{u}_e}{r{\theta }_v\cos \varphi }\frac{\partial {\theta }_v}{\partial \lambda }-\frac{g{v}_e}{r{\theta }_v}\frac{\partial {\theta }_v}{\partial \varphi }+\frac{g{S}_{{\theta }_v}}{{\theta }_v}}\end{array}\right).$$

(23)

Since $\mathbf{B}$ and $\mathbf{Q}$ involve values at adjacent levels, it is impractical to eliminate $\mathbf{u}$ directly from Equation (20) to solve for $\frac{\partial \pi }{\partial t}$. To solve Equation (20) reasonably accurately, a preconditioning is required. To invert ${\mathbf{BQ}}^{\prime }$, we therefore approximate ${\mathbf{BQ}}^{\prime }$ by an invertible matrix $\mathbf{P}$, so that Equation (20) is replaced by

$$\mathbf{P}\left(\begin{array}{c}{\displaystyle u-u_e}\\ {\displaystyle v-v_e}\\ {\displaystyle w}\end{array}\right)+c_{pd}{\theta }_v\frac{\partial }{\partial t}\left(\begin{array}{c}{\displaystyle \frac{1}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)}\\ {\displaystyle \frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)}\\ {\displaystyle \frac{\partial \pi }{\partial r}}\end{array}\right)=\mathbf{G},$$

(24)

where $\mathbf{P}={\mathbf{P}}_1{\mathbf{P}}_2{\mathbf{P}}_3$. ${\mathbf{P}}_1$ is given by

$${\mathbf{P}}_1=\left(\begin{array}{ccc}{f}^2+\frac{\partial }{\partial r}\left(K_m\frac{{\partial }^2}{\partial r^2}\left(K_m\frac{\partial }{\partial r}\right)\right)& 0& 0\\ 0& f^2+\frac{\partial }{\partial r}\left(K_m\frac{{\partial }^2}{\partial r^2}\left(K_m\frac{\partial }{\partial r}\right)\right)& 0\\ 0& 0& 1\end{array}\right).$$

(25)

${\mathbf{P}}_2$ is given by

$${\mathbf{P}}_2=\left(\begin{array}{ccc}{\displaystyle f^2+\frac{f{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\frac{v_e}{{\theta }_v}+\frac{f{u}_e\tan \varphi }{r}}& {\displaystyle \frac{f{\theta }_v}{r}\frac{\partial }{\partial \varphi }\frac{v_e}{{\theta }_v}}& {\displaystyle f{\theta }_v\frac{\partial }{\partial r}\frac{v_e}{{\theta }_v}}\\ {\displaystyle -\frac{f{\theta }_v}{r\cos \varphi }\frac{\partial }{\partial \lambda }\frac{u_e}{{\theta }_v}+\frac{f{v}_e\tan \varphi }{r}}& {\displaystyle f^2-\frac{f{\theta }_v}{r}\frac{\partial }{\partial \varphi }\frac{u_e}{{\theta }_v}}& {\displaystyle -f{\theta }_v\frac{\partial }{\partial r}\frac{u_e}{{\theta }_v}}\\ {\displaystyle \frac{g}{r\cos \varphi {\theta }_v}\frac{\partial {\theta }_v}{\partial \lambda }}& {\displaystyle \frac{g}{r{\theta }_v}\frac{\partial {\theta }_v}{\partial \varphi }}& {\displaystyle \frac{g}{{\theta }_v}\frac{\partial {\theta }_v}{\partial r}}\end{array}\right),$$

(26)

and

$${\mathbf{P}}_3=\left(\begin{array}{ccc}{f}^{-2}& 0& 0\\ 0& f^{-2}& 0\\ 0& 0& 1\end{array}\right).$$

(27)

We next have to eliminate $\mathbf{u}$ using the continuity equation and equation of state, which are the fourth and fifth parts of Equation (1). In $(\lambda ,\varphi ,\eta )$ coordinates, we have $\mathbf{u}=(u,v,\dot{\eta })$ and

$$\dot{\eta }\frac{\partial r}{\partial \eta }=w-\frac{u}{r\cos \varphi }\frac{\partial r}{\partial \lambda }-\frac{v}{r}\frac{\partial r}{\partial \varphi }.$$

(28)

This can be rewritten as

$$\mathbf{R}\left(\begin{array}{c}{\displaystyle u}\\ {\displaystyle v}\\ {\displaystyle w}\end{array}\right)=\left(\begin{array}{c}{\displaystyle u}\\ {\displaystyle v}\\ {\displaystyle \dot{\eta }}\end{array}\right).$$

The continuity equation then takes the form

$$\frac{\partial }{\partial t}\left(r^2{\rho }_d\frac{\partial r}{\partial \eta }\right)+\frac{1}{\cos \varphi }\left(\frac{\partial }{\partial \lambda }\left(r{\rho }_{d}u\frac{\partial r}{\partial \eta }\right)+\frac{\partial }{\partial \varphi }\left(r{\rho }_{d}v\cos \varphi \frac{\partial r}{\partial \eta }\right)+\frac{\partial }{\partial \eta }\left(r^2{\rho }_d\dot{\eta }\frac{\partial r}{\partial \eta }\right)\right)=0$$

(29)

Now, applying the operator

$$\nabla ·r^2{\rho }_d\frac{\partial r}{\partial \eta }\mathbf{R}$$

to Equation (24) and using Equation (29) gives

$$\begin{array}{c}-\frac{\partial }{\partial t}\left(r^2{\rho }_d\frac{\partial r}{\partial \eta }\right)+\nabla ·\left(r^2{\rho }_d{c}_{pd}{\theta }_v\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\frac{\partial }{\partial t}\left({\nabla }_r\pi \right)\right)=\hfill \\ \nabla ·\left(r^2{\rho }_d\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\mathbf{G}\right)+\frac{1}{\cos \varphi }\left(\frac{\partial }{\partial \lambda }\left(r{\rho }_d{u}_e\frac{\partial r}{\partial \eta }\right)+\frac{\partial }{\partial \varphi }\left(r{\rho }_d{v}_e\cos \varphi \frac{\partial r}{\partial \eta }\right)\right).\hfill \end{array}$$

(30)

Next, differentiating the equation of state with respect to time gives

$$\frac{{\kappa }_d-1}{{\kappa }_d}{\pi }^{-\frac{1}{{\kappa }_d}}\frac{\partial \pi }{\partial t}{\rho }_d{\theta }_{vd}+{\pi }^{\frac{{\kappa }_d-1}{{\kappa }_d}}\frac{\partial {\rho }_d}{\partial t}{\theta }_{vd}+{\pi }^{\frac{{\kappa }_d-1}{{\kappa }_d}}{\rho }_d\frac{\partial {\theta }_{vd}}{\partial t}=0.$$

(31)

Removing common factors, and using Equation (18) and the final part of Equation (2) to evaluate the final term of Equation (31) gives

$$\frac{{\kappa }_d-1}{{\kappa }_d}\frac{\partial \pi }{\partial t}{\rho }_d{\theta }_{vd}+\pi \frac{\partial {\rho }_d}{\partial t}{\theta }_{vd}-\pi {\rho }_d\frac{c_{pd}{\theta }_v^2}{g}\frac{\partial }{\partial t}\left(\frac{1}{(1+\mu )}\frac{\partial \pi }{\partial r}\right)=0.$$

(32)

We now substitute Equation (32) into Equation (30) to give

$$\begin{array}{c}\frac{{\kappa }_d-1}{{\kappa }_d}\frac{{\rho }_d}{\pi }{r}^2\frac{\partial r}{\partial \eta }\frac{\partial \pi }{\partial t}+r^2\frac{\partial r}{\partial \eta }\frac{{\rho }_d{c}_{pd}{\theta }_v^2}{g{\theta }_{vd}}\frac{\partial }{\partial t}\left(\frac{1}{(1+\mu )}\frac{\partial \pi }{\partial r}\right)+\nabla ·\left(r^2{\rho }_d{c}_{pd}{\theta }_v\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\frac{\partial }{\partial t}\left({\nabla }_r\pi \right)\right)=\hfill \\ \nabla ·\left(r^2{\rho }_d\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\mathbf{G}\right)+\frac{1}{\cos \varphi }\left(\frac{\partial }{\partial \lambda }\left(r{\rho }_d{u}_e\frac{\partial r}{\partial \eta }\right)+\frac{\partial }{\partial \varphi }\left(r{\rho }_d{v}_e\cos \varphi \frac{\partial r}{\partial \eta }\right)\right).\hfill \end{array}$$

(33)

Apart from the small term $\frac{\partial }{\partial t}\frac{1}{(1+\mu )}$, which needs to be diagnosed separately, this is an elliptic equation for $\partial \pi /\partial t$, which can be solved. Boundary conditions are required at the top and bottom. The boundary conditions inherited from Equation (1) imply that $\dot{\eta }=0$ at the top and bottom, which implies a zero gradient of

$$r^2{\rho }_d{c}_{pd}{\theta }_v\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\frac{\partial }{\partial t}\left({\nabla }_r\pi \right).$$

(34)

Equation (28) and the no slip lower boundary condition then implies that $w=0$ at the lower boundary also. Then, Equations (20) and (22) reduce to

$$c_{pd}{\theta }_v\frac{\partial \pi }{\partial r}=S_{{\theta }_v}.$$

(35)

We can then calculate the ageotriptic velocity $u-u_e,v-v_e$ and w by applying ${\mathbf{P}}^{-1}$ to Equation (24).

2.3. Application

The right hand side of Equation (33) shows that the evolution is driven by the divergence of the forcing vector $\mathbf{G}$ and the divergence of the geotriptic wind $(u_e,v_e)$. Changes to the forcing with zero divergence, in the sense of the first term on the right hand side of Equation (33), will not affect the evolution. Such changes will simply alter the ageostrophic wind vector $(u-u_e,v-v_e,w)$ so that Equation (20) is satisfied. This demonstrates that changes to the forcing in the UM may have a much smaller effect on the large-scale evolution than might be expected.

Using the definition of $\mathbf{G}$ in Equation (23), we see that the first term on the right hand side of Equation (33) includes geotriptic advection of the momentum and potential temperature, the dynamical forcing, together with the external forcing terms S which are generated by the UM physics. The second term on the right hand side would be just the divergence of the geostrophic wind if the boundary later were not included. This would be zero in pressure coordinates with a constant Coriolis parameter. There is a contribution from variations in the Coriolis parameter, which gives the effect of Rossby wave propagation on the pressure tendency. Since the divergence calculation used to derive Equation (33) is in terrain-following coordinates, the orography will also contribute to this term through the horizontal variations of $\frac{\partial r}{\partial \eta }$. This shows how the orography can have a strong effect on the large-scale flow.

Moisture will interact strongly with the diagnostic procedure through the reduction of the effective static stability by latent heat release. This effect has been studied, for instance, by Booth et al. [28]. In the UM, the thermodynamic equation, which is the second part of Equation (1), includes the terms

$$w\frac{\partial {\theta }_{vd}}{\partial r}=S_{LH},$$

where $S_{LH}$ is the part of the forcing term resulting from latent heat release. This term is dominated by the term $\frac{D{q}_{sat}}{Dt}$, where $q_{sat}$ is the saturated mixing ratio which is largely a function of temperature. In a simple form, this can be represented by calculating a cloud fraction weighted static stability, where $\partial {\theta }_v/\partial r$ is the stability of the unsaturated air allowing for the virtual temperature effect and $\partial {\theta }_e^{*}/\partial z$ is the moist stability, where ${\theta }_e^{*}$ is the equivalent potential temperature calculated assuming saturation. w is not known in advance, so the idea that latent heating only occurs with upward motion cannot be used as the problem becomes nonlinear in w. The use of cloud fraction weighting avoids this problem, as evaporation will occur if w is negative in cloudy regions. If $\alpha$ is the cloud fraction, we set

$${\left(\frac{\partial {\theta }_v}{\partial r}\right)}_{eff}=\frac{\partial (\alpha {\theta }_e^{*}+(1-\alpha ){\theta }_v))}{\partial r}.$$

(36)

A more sophisticated approach would be to use a formula set out in [29] for wave propagation in inhomogeneous media. Applying this to a wet/dry mixture gives a modified Brunt–Väisälä frequency as

$$N_{eff}^2={\left(\frac{\alpha N_d^2+(1-\alpha )N_w^2}{N_w^2{N}_d^2}\right)}^{-1},$$

(37)

where

$$N_w^2=\frac{g}{\theta }\frac{\partial {\theta }_e^{*}}{\partial r},N_d^2=\frac{g}{\theta }\frac{\partial {\theta }_v}{\partial r}.$$

(38)

We thus replace $\frac{\partial {\theta }_v}{\partial r}$ in ${\mathbf{Q}}^{\prime }$, Equation (21), by

$${\left(\frac{\partial {\theta }_v}{\partial r}\right)}_{eff}={\left(\frac{\alpha \frac{\partial {\theta }_v}{\partial r}+(1-\alpha )\frac{\partial {\theta }_e^{*}}{\partial r}}{\frac{\partial {\theta }_e^{*}}{\partial r}\frac{\partial {\theta }_v}{\partial r}}\right)}^{-1}.$$

(39)

The model value of cloud fraction $\alpha$ can be used, provided again that small space and time scales are excluded.

If the temperature increment $S_L$ from the large scale precipitation scheme is added to S, the solution of Equation (20) will become

$$\mathbf{B}{\mathbf{Q}}^{\prime }\left(\begin{array}{c}{\displaystyle u-u_e}\\ {\displaystyle v-v_e}\\ {\displaystyle w}\end{array}\right)+c_{pd}{\theta }_v\frac{\partial }{\partial t}\left(\begin{array}{c}{\displaystyle \frac{1}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)}\\ {\displaystyle \frac{1}{r}\left(\frac{\partial \pi }{\partial \varphi }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \varphi }\right)}\\ {\displaystyle \frac{\partial \pi }{\partial r}}\end{array}\right)=\mathbf{B}\left({\mathbf{H}}^{\prime }+\left(\begin{array}{c}{\displaystyle \mathbf{0}}\\ {\displaystyle \mathbf{0}}\\ {\displaystyle {\mathbf{S}}_{\mathbf{L}}}\end{array}\right)\right).$$

(40)

The solution of Equation (40) for $\frac{\partial \pi }{\partial t}$ is identical to that of Equation (20) if the element $\frac{\partial {\theta }_v}{\partial r}$ in ${\mathbf{Q}}^{\prime }$ in that equation is replaced by ${\left(\frac{\partial {\theta }_v}{\partial r}\right)}_{eff}$ as in Equation (39) and

$$S_L=w\left(\frac{\partial {\theta }_v}{\partial r}-{\left(\frac{\partial {\theta }_v}{\partial r}\right)}_{eff}\right).$$

(41)

This will be approximately true if the latent heat release is dominated by the effects of vertical motion. This is illustrated in Section 3.3.

2.4. Computational Aspects

The elliptic equation (33) is solved by a preconditioned generalised conjugate residual method which requires explicit evaluation of the elliptic operator on the left hand side for a given estimate of $\partial \pi /\partial t$, and an evaluation of the right hand side. The preconditioning is done by a vertical matrix inversion, which uses a tridiagonal solver, and a horizontal smoother.

Given $\pi ,K_m$ from model data, solve Equations (10) and (11) for $(u_e,v_e)$. Care has to be taken in calculating the right hand side of Equation (10) because the two terms cancel strongly over orography. The scheme, based on that used in the UM and exploiting hydrostatic balance, is

$$\frac{c_{pd}{\theta }_v}{r\cos \varphi }\left(\frac{\partial \pi }{\partial \lambda }-\frac{\partial \pi }{\partial r}\frac{\partial r}{\partial \lambda }\right)\simeq c_{pd}{{\overline{\theta }}_v}^{\lambda }{\delta }_{\lambda }\pi +g{\delta }_{\lambda }r,$$

(42)

where ${\delta }_{\lambda }$ indicates a central difference in the $\lambda$ direction divided by $2r\cos \varphi \delta \lambda$. A similar scheme is used in the $\varphi$ direction. Solving Equation (11) then involves inverting a $N\times N$ penta-diagonal matrix at each horizontal gridpoint, where N is the number of model levels and $K_m$ is nonzero.

Next, use $(u_e,v_e)$ to evaluate the components of ${\mathbf{P}}_2$, Equation (26). Central differencing has to be used because the upwind direction is not known in advance. The vertical stability is evaluated more accurately as a vertical second difference of $\pi$. The values are filtered towards a zonal mean close to the equator in order to remain within the validity of the theory. If Equation (36) or (39) is used, the relative humidity, moisture content and cloud fraction are used to calculate the effective static stability.

Next, calculate the eigenvalues of ${\mathbf{P}}_2$ at each point. Modify ${\mathbf{P}}_2$ as necessary to remove negative eigenvalues. This can be done by first replacing any negative diagonal terms by suitable reference values, such as ${10}^{-10}$ in the first two rows and ${10}^{-6}$ in the third, and then rescaling the off-diagonal terms as required.

Next, construct the vertical pre-conditioner. If ${\mathbf{P}}_{2ij}$ is the value of the $(i,j)$ component of ${\mathbf{P}}_2$, calculate the rms value ${\mathbf{PI}}_{rms}$ of ${\mathbf{P}}_{233}^{-1}$ as a function of model level. Then, at each horizontal gridpoint, we precondition Equation (33) by replacing the left hand side by

$$\mathbf{A}\frac{\partial \pi }{\partial t}\equiv C_1\ast \frac{{\kappa }_d-1}{{\kappa }_d}\frac{{\rho }_d}{\pi }{r}^2\frac{\partial r}{\partial \eta }\frac{\partial \pi }{\partial t}+\nabla ·\left(r^2{\rho }_d{c}_{pd}{\theta }_v\frac{\partial r}{\partial \eta }{\mathbf{RPI}}_{\mathbf{v}}\frac{\partial }{\partial t}\left({\nabla }_r\pi \right)\right).$$

(43)

$C_1$ is set to 100 in the examples illustrated. ${\mathbf{PI}}_v$ is a tridiagonal matrix at each gridpoint defined as

$${\mathbf{PI}}_v={\mathbf{C}}_2{\mathbf{PI}}_{rms}+({\mathbf{P}}_{211}^{-1}+{\mathbf{P}}_{222}^{-1}).$$

(44)

${\mathbf{C}}_2$ is a matrix $\left(\begin{array}{ccc}1& 0& 0\\ 0& C_2& 0\\ 0& 0& 1\end{array}\right)$, with $C_2$ set to 1.5 in the examples illustrated.

Now, solve Equation (33) for $\frac{\partial \pi }{\partial t}$ preconditioned with ${\mathbf{A}}^{-1}$ as defined in Equation (43) using a generalised conjugate residual method. Back substitute in Equation (33) with the left hand side modified to

$$\begin{array}{c}\hfill \frac{{\kappa }_d-1}{{\kappa }_d}\frac{{\rho }_d}{\pi }{r}^2\frac{\partial r}{\partial \eta }\frac{\partial \pi }{\partial t}+r^2\frac{\partial r}{\partial \eta }\frac{{\rho }_d{c}_{pd}{\theta }_v^2}{g{\theta }_{vd}}\frac{\partial }{\partial t}\left(\frac{1}{(1+\mu )}\frac{\partial \pi }{\partial r}\right)+\\ \hfill C_3\nabla ·\left(r^2{\rho }_d{c}_{pd}{\theta }_v\frac{\partial r}{\partial \eta }{\mathbf{RP}}^{-1}\frac{\partial }{\partial t}\left({\nabla }_r\pi \right)\right).\end{array}$$

(45)

$C_3$ is a constant set to 1.2 in the examples illustrated. Carry out iterations until the convergence slows (25 in the examples illustrated).

Now, back substitute in Equation (33) and calculate a residual. Repeat the inner iterations. This outer iteration was run 10 times in the examples illustrated. Finally, use Equation (24) to obtain $\mathbf{u}$. It is not currently practical to iterate this procedure further to solve Equation (20) exactly.

3. Results

3.1. Experimental Setup

The diagnostic can be calculated from UM data with whatever resolution is to be studied. In the examples quoted here, this was on a latitude longitude grid with $640\times 480$ points, with a north–south gridlength of about 40 km. The data have 70 levels extending to 80 km. The diagnostic is run on a grid of $160\times 120$ points, giving a north–south gridlength of about 125 km, but the same levels as the UM. The illustrations use a single case from February 2014.

This procedure allows reasonable computational accuracy in computing scales of 500 km or greater. However, it will still capture variability on scales smaller than where the SGT approximation is valid. Therefore, the forcing fields are further smoothed. The coefficients of the ${\mathbf{Q}}^{\prime }$ matrix calculated from the data will also emphasise smaller scales because of the derivatives used to define the coefficients. These are also smoothed. Near the equator an additional filter is applied to the geostrophic winds which relaxes them towards the zonal mean. This is needed because of the severe restrictions on the validity of SGT in the tropics outside the boundary layer.

3.2. Comparison of Diagnostic and Model-Derived Ageotriptic Winds

Root mean square values of the geotriptic wind are about 21 ms${}^{-1}$ for this data, and root mean square values of the ageotriptic wind are about 6.3 ms${}^{-1}$, giving a global averaged Rossby number of about 0.3. This would suggest a 10% expected error in the difference between the model wind and the total wind deduced from Equation (20). The actual global rms difference between the diagnosed wind and the model wind is 9.7 ms${}^{-1}$ in a overall rms wind of 22 ms${}^{-1}$. The global rms difference between the calculated geostrophic wind and the model wind is 7.3 ms${}^{-1}$.

To illustrate this, Figure 1 shows a mid-tropospheric zonal geotriptic wind calculated using Equations (10) and (11) from UM pressure and potential temperature fields. The area illustrated is the extratropical North Atlantic, so away from significant topography. It would be reasonable to expect the ageotriptic wind to be accurate over such an area. Figure 1 shows the ageotriptic wind calculated from Equation (20) and the difference between the UM wind and the geotriptic wind deduced from the pressure and temperature. There is a reasonable match between the two estimates of the ageotroptic wind, indicating that SGT theory is working reasonably accurately for the UM fields over the area chosen. The rms ageotriptic wind calculated from Equation (20) is 5.1 ms${}^{-1}$ over the area and level shown in Figure 1. The rms difference between the UM wind and the geotriptic wind is is 4.7 ms${}^{-1}$ in an overall rms wind of 22 ms${}^{-1}$, so the proportional error in the geotriptic approximation is 20% rather than 45% globally. The rms difference between the two estimates of ageotriptic wind illustrated in Figure 1 is 3.8 ms${}^{-1}$. If these fields are smoothed horizontally on the scale of the diagnostic grid, about 150 km, this difference reduces to 3.0 ms${}^{-1}$. The correlation coefficient between the two estimates of the ageotriptic wind is 0.69 over the area and level illustrated. This both illustrates the applicability of the theory and the ability of the UM to reproduce it.

Figure 1 shows that the diagnosed ageotriptic wind is larger than the model value in some regions. This is typical of situations where the space or time scale is becoming too small for SGT theory to be accurate, and leads to the rms difference between the diagnosed wind and the model wind being no smaller than that between the calculated geotriptic wind and the model wind.

3.3. Use of a Modified Static Stability to Represent Latent Heat Release

In Section 2.3, it is shown that the latent heating resulting from maintenance of geotriptic balance in cloudy air can be deduced by solving Equation (20) with a modified static stability as in Equation (39). The resulting w can then be used in Equation (41) to deduce the latent heating. This can be compared with the temperature increments from the schemes representing cloud and precipitation, both convective and dynamic, in the UM.

The results are shown in Figure 2 for the same case as used in Figure 1. The area chosen is again the extratropical North Atlantic, so it can be expected that the vertical motion required to maintain geotriptic balance will be similar to that generated by the UM, and that precipitation processes, which generate most of the latent heating, will be strongly linked to the vertical velocity. In the first two panels, the vertical velocity w diagnosed from Equation (20) with Equation (39) is compared with that directly output from the UM. As in the comparison of ageotriptic winds in Figure 1, there is reasonable agreement. Comparing with Figure 1a, which shows the position of the strongest zonal winds, we see that there is systematic upward motion at the downstream end of the region of strongest zonal winds, and downward motion in the centre and rearward parts of this region. This is consistent with standard synoptic experience. The diagnosed vertical velocity has more coherent large scale structure, while the direct model output includes more small scale features, as would be expected. In the second pair of panels, the latent heating deduced from Equation (41) is compared with the temperature increments from the large-scale cloud scheme and the convection scheme combined. There is a strong correspondence between the diagnosed latent heating and the diagnosed vertical motion, as would be expected from Equation (39). The correspondence between the model vertical motion and the model latent heating does not look as strong. However, the correlation coefficient between them for the plotted area and level is again 0.69. The visual mismatch is primarily because there is no latent heating in areas of downward motion. The diagnostic procedure will generate latent cooling in cloudy air, but this would result in rapid dissipation of the cloud and not correspond to significant latent cooling in the UM.

3.4. Effect of Tropical Heating on the Subtropical Jet

The diagnostic procedure can be applied using artificial physical forcing as well as that derived from the model. There has been considerable interest in the effect of tropical heating on the extratropical circulation. An example is the paper by Matthews et al. [30], where time dependent heating deduced from reanalyses is used to determine the response to forcing of a simple general circulation model (GCM) using a climatological mean state. This is compared with the observed evolution deduced from satellite data.

The diagnostic presented here calculates the “instant” response of the geostrophic flow to forcing at a particular time. It is thus easiest to compare it with the earlier study of [22]. In that paper, the forcing was maintained over a long timescale, and an actual atmospheric state as well as a climatological state could be used to determine the response to the forcing. The atmospheric state was maintained in time by using an artificial forcing term. We apply a forcing similar to [22] and calculate the instant geostrophic response to this forcing using the atmospheric state illustrated in the rest of the paper. In [22], Figure 17c shows the divergent response to the forcing, which is set up very quickly. This should be reproduced by the diagnostic. Figure 18c also shows a significant downstream response at about 30${}^{\circ }$ N in the subtropical jet after five days. Since Figure 17c shows that the divergent response does not propagate much downstream, the formula for the source term plotted in Figure 18c shows that this must represent a perturbation vorticity. Thus, in our diagnostic, we look for a localised perturbation to the geotriptic flow in the subtropical jet, which would then propagate downstream over the five-day period illustrated by Jin et al. [22].

The heating in the mid troposphere is illustrated in Figure 3a. It is centred on the dateline as in [22], Figure 17c. The heating extends through the troposphere. The zonal geostrophic wind is shown in Figure 3b in the upper troposphere, noting that strong filtering has been applied near the equator to eliminate values for which SGT theory is inappropriate. The diagnosed zonal geostrophic wind tendency is shown in Figure 3c. This peaks at the same longitude as the forcing, so would propagate the jet forwards. The values peak at about 0.5 ${}^{\circ }$K day${}^{-1}$. If this is applied over 15 days, this gives a similar impact to that shown in Figure 12a of [22]. There are also impacts at lower latitudes, but these are not likely to be physically correct as they are outside the validity of SGT.

The reason for the strong interaction with the subtropical jet can be seen in Figure 3d. This plots the second diagonal element of the matrix ${\mathbf{BQ}}^{\prime }$ defined in Equations (17) and (21). In the absence of the non-trivial model state this would take the value $f^2$, about $0.5\times {10}^{-8}$ at 30${}^{\circ }$N. It can be seen that on the southern flank of the jet, where there is strong anticyclonic shear, the values drop to well below this. The effect can be seen from Equation (20). If the matrix ${\mathbf{Q}}^{\prime }$ has a small eigenvalue, there will be a very strong response to forcing in this direction, so that in the present case there is a very efficient transfer of information from the heating to the subtropical jet. Comparing Figure 3a,c shows that the largest impact is to the west of the forcing longitude, where the matrix coefficient is smallest. The diagnostic can thus pick up sensitivity of the response to forcing to the structure of the basic state, which is a topic of wide interest.

3.5. Effect of Boundary Layer Heating

A feature of the diagnostic procedure is the inclusion of a proper treatment of the boundary layer. As shown in Equation (35), the boundary layer heating projects strongly onto the temperature, as the vertical motion is suppressed near the surface. This heating will also allow a horizontal temperature gradient to develop in the boundary layer. Since there is no horizontal pressure gradient above the boundary layer, hydrostatic balance will generate a heat low near the surface, resulting in convergence of the geotriptic wind which flows down the pressure gradient. Hence, upward motion will be generated at the top of the boundary layer. This effect forms part of the explanation of the heat lows that often form over tropical land masses [31].

This is particularly important in the tropics, where it allows representation of the return circulation consistent with the upward motion in the deep tropics. The example shown in Figure 4 illustrates the response to boundary layer heating over Borneo. The local time is around 1000 (0200UTC). There is strong convergence over both Borneo and Celebes to the east, leading to upward motion (not shown), and also anticyclonic rotation.

4. Discussion

The paper has developed a diagnostic procedure based on semi-geotriptic theory that can be applied successfully to data from a comprehensive global model. It is necessary to filter the data to large horizontal scales, particularly in the tropics, as expected. This means that it can be used to extract the “balanced” response to forcing, either by the model dynamics and physics, or artificially imposed forcing. The first example shows that the SGT approximation gives results which match the large-scale behaviour of the UM sufficiently closely to be useful. In particular, the predicted ageostrophic flow correlates quite well with the UM’s ageostrophic flow. Greater accuracy cannot realistically be expected given the wide range of scales and forcing mechanisms present in real and model-simulated flow. An alternative view is to say that the UM appears to reproduce the balanced response to forcing as calculated by the diagnostic.

The diagnostic calculates the large-scale response of the atmosphere to forcing, and thus allows studying the effect of both changes to the forcing and changes to the background flow. This will contribute to the understanding of large-scale circulation patterns, such as monsoons. In particular, it should aid the study of systematic errors, which usually have a large-scale signal. The predicted circulation is strongly dependent on the model state, so that errors in the resulting circulation can come from errors in the model state as well as errors in the forcing. In this application, the results of the diagnostic procedure would need to be averaged over a large number of cases, so that the noise created by using data on the borderline for applicability of SGT would be filtered out. The second example supports this by showing that much of the extratropical precipitation signal can be reproduced using the diagnosed vertical motion and the cloud fraction. This will help in identifying the causes of systematic errors in precipitation. The third example shows that the effect of tropical–extratropical interaction is strongly dependent on the background state assumed; the causes of errors in this could be distinguished by replacing the model state with a reanalysis, while still using the same forcing. The final example shows the leading order effect of boundary layer forcing on the circulation, again allowing the causes of systematic errors to be identified.