Three-dimensional normal mode functions: open-access tools for their computation in isobaric coordinates (p-3DNMF.v1)

A free software package for the computation of the three-dimensional normal modes of an hydrostatic atmosphere is presented. This software performs the computations in isobaric coordinates and was developed for two user-friendly languages: MATLAB and Python. The software can be used to expand the global atmospheric circulation onto the 3-D normal modes. This expansion allows the computation of a 3-D energetic scheme, which partitions the energy reservoirs and energy interactions between 3-D spatial scales, barotropic and baroclinic components, and balanced (rotational) and unbalanced (divergent) circulation fields. Moreover, by retaining only a subset of the expansion coefficients, the 3-D normal mode expansion can be used for spatial-scale filtering of atmospheric motion, filtering of balanced motion and mass unbalanced motions, and barotropic and baroclinic components. Fixing the meridional scale, the 3-D normal mode filtering can be used to isolate tropical components of the atmospheric circulation. All of these features are useful both in data analysis and in assessment of general circulation atmospheric models.

1 Introduction

The use of the three-dimensional normal mode functions (3-D NMFs) of the linearised primitive equations as a basis to expand the global horizontal wind and mass fields simultaneously was presented for the first time by Kasahara and Puri (1981). The expansion of both atmospheric fields onto the 3-D NMFs allows the partition of the energy of the global motion into two kinds of motion: mass balanced rotational modes (Rossby–Haurwitz waves) and divergent modes (inertio-gravity waves).

The original work of Kasahara and Puri (1981) was developed for terrain, following sigma coordinates. Tanaka (1985) extended the methodology to isobaric coordinates 4 years later. In these coordinates, Tanaka (1985) and Tanaka and Kung (1988) were able to develop a 3-D normal mode energy scheme, including the exchanges of kinetic and available potential energy (APE) among the different three-dimensional scales and types of motions. Recently, Marques and Castanheira (2012) extended the methodology to analyse the conversion of APE into kinetic energy.

The 3-D normal mode expansion has been applied in several type of studies, such as the comparisons of reanalysis datasets (Marques and Castanheira, 2012; Yamagami and Tanaka, 2016, and references therein), or the analysis of ensemble prediction system uncertainties (Žagar et al., 2015a). The methodology was also applied in studies of low-frequency extratropical variability, like the Arctic Oscillation–North Atlantic Oscillation (AO–NAO) (Castanheira et al., 2002; Tanaka and Tokinaga, 2002; Tanaka and Seki, 2013; Castanheira and Marques, 2019), as well as in tropical variability analysis (Yamagami and Tanaka, 2014; Castanheira and Marques, 2015; Marques and Castanheira, 2018; Blaauw and Žagar, 2018; Franzke et al., 2019; Kitsios et al., 2019).

Recently, Žagar et al. (2015b) developed an open-source software package for the projection of global horizontal wind and mass fields onto 3-D NMFs. Their software was developed for sigma coordinates but does not allow for the analysis of the full energy cycle. Here we present a software package for the projections of horizontal wind and mass fields onto 3-D NMFs using isobaric coordinates. This software allows for the analysis of the full atmospheric energy cycle and was developed for two user-friendly languages: MATLAB and Python. The methodology described in this study to obtain the 3-D NMFs has been used in previous works, Marques and Castanheira (e.g. 2012); Marques et al. (e.g. 2014); Castanheira and Marques (e.g. 2015); Marques and Castanheira (e.g. 2018); Castanheira and Marques (e.g. 2019), in the form of unstructured and disperse scripts written for the specific dataset used in those studies. The software presented here consists of a structured collection of functions, intended to be used for a wider and more general datasets (such as reanalysis and outputs from climate models) and requiring the minimum amount of user interaction as we found to be possible.

2 3-D normal mode functions

2.1 The basic equations

In the traditional shallow atmosphere approximation, a set of primitive equations in isobaric coordinates may be written as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial u}{\partial t}}-\mathrm{2}\mathrm{\Omega }v\sin \mathit{\theta }+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\lambda }}}=\\ & -\mathit{V}\cdot \mathrm{\nabla }u-\mathit{\omega }{\displaystyle \frac{\partial u}{\partial p}}+{\displaystyle \frac{\tan \mathit{\theta }}{a}}uv+F_u,\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial v}{\partial t}}+\mathrm{2}\mathrm{\Omega }u\sin \mathit{\theta }+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\theta }}}=\\ & -\mathit{V}\cdot \mathrm{\nabla }v-\mathit{\omega }{\displaystyle \frac{\partial v}{\partial p}}-{\displaystyle \frac{\tan \mathit{\theta }}{a}}{u}^{\mathrm{2}}+F_v,\end{array}\text{(3)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial p}}=-{\displaystyle \frac{RT}{p}},\text{(4)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{V}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}=\mathrm{0},\text{(5)}& {\displaystyle \frac{\partial T}{\partial t}}-{\displaystyle \frac{p{S}_{\mathrm{0}}}{R}}\mathit{\omega }=-\mathit{V}\cdot \mathrm{\nabla }T-\mathit{\omega }{\displaystyle \frac{\partial T}{\partial p}}+{\displaystyle \frac{\dot{Q}}{c_p}},\end{array}$$

where λ and θ are longitude and latitude, respectively. Variables u and v are the zonal and meridional components of horizontal wind vector V, and the vertical wind component in the isobaric system is represented by ω. The constants Ω and a denote the Earth's angular velocity and Earth's radius, respectively. Variable $\dot{Q}$ is the rate of diabatic heating per mass unity, c_p the specific heat at constant pressure (p), and R is the dry air gas constant. Temperature and geopotential, T and ϕ, are the deviations from a hydrostatic reference state T₀(p) and ϕ₀(p), respectively, and the static stability parameter, S₀, is given by

$$\begin{array}{}\text{(6)}& {\displaystyle }{S}_{\mathrm{0}}={\displaystyle \frac{R}{p}}\left({\displaystyle \frac{R{T}_{\mathrm{0}}}{p{c}_p}}-{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{0}}}{\mathrm{d}p}}\right).\end{array}$$

Multiplying Eq. (5) by R∕(pS₀), then calculating the derivatives with respect to p, and finally using the hydrostatic and continuity Eqs. (3) and (4) on the resulting equation, the thermodynamic energy equation takes the form

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& -{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial p}}\right)+\mathrm{\nabla }\cdot \mathit{V}=\\ & {\displaystyle \frac{\partial }{\partial p}}\left[{\displaystyle \frac{R}{p{S}_{\mathrm{0}}}}\left(-\mathit{V}\cdot \mathrm{\nabla }T-\mathit{\omega }{\displaystyle \frac{\partial T}{\partial p}}\right)\right]+{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{R\dot{Q}}{p{c}_p{S}_{\mathrm{0}}}}\right).\end{array}\end{array}$$

The right-hand side of Eqs. (1), (2), and (7) contain the non-linear terms, frictional forces, and the diabatic heat sources. By setting the right-hand sides of those equations to zero, the linearised system of equations for the three dependent variables u, v and ϕ is written as

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial u}{\partial t}}-\mathrm{2}\mathrm{\Omega }v\sin \mathit{\theta }+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\lambda }}}=\mathrm{0},\text{(9)}& {\displaystyle \frac{\partial v}{\partial t}}+\mathrm{2}\mathrm{\Omega }u\sin \mathit{\theta }+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\theta }}}=\mathrm{0},\text{(10)}& {\displaystyle }-{\displaystyle \frac{\partial }{\partial p}}\left[{\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}\right)\right]+\mathrm{\nabla }\cdot \mathit{V}=\mathrm{0}.\end{array}$$

The linearised system (8)–(10) describes small oscillations of an incompressible, homogeneous, hydrostatic and inviscid fluid over a rotating sphere, around a state at rest (Tanaka and Kung, 1988; Daley, 1991). Assuming the following separation of variables (Kasahara and Puri, 1981; Daley, 1991),

$$\begin{array}{}\text{(11)}& {\displaystyle }\left[u,v,\mathit{\varphi }\right]=\mathrm{\Psi }\left(p\right)\left[\stackrel{\mathrm{̃}}{u}(\mathit{\lambda },\mathit{\theta },t),\stackrel{\mathrm{̃}}{v}(\mathit{\lambda },\mathit{\theta },t),\stackrel{\mathrm{̃}}{\mathit{\varphi }}(\mathit{\lambda },\mathit{\theta },t)\right],\end{array}$$

and substituting into Eq. (10), one obtains

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Psi }}}{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathrm{\Psi }}{\partial p}}\right)={\displaystyle \frac{\mathrm{\nabla }\cdot \stackrel{\mathrm{̃}}{\mathit{V}}}{\partial \stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial t}}=-{\displaystyle \frac{\mathrm{1}}{gh}},\end{array}$$

where g is the constant gravity acceleration and 1∕(gh) is the separation constant.

Following this, the separable vertical structures Ψ of the solutions of the linearised system (8)–(10) are given by the eigensolution of the following vertical structure equation (VSE)

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathrm{\Psi }}{\partial p}}\right)+{\displaystyle \frac{\mathrm{1}}{gh}}\mathrm{\Psi }=\mathrm{0},\end{array}$$

satisfying an appropriate pair of upper and lower boundary conditions.

2.2 Vertical structure functions

The VSE (13), with appropriate boundary conditions, forms an eigenvalue problem with discrete eigensolutions Ψ_k(p) and eigenvalues $-\mathrm{1}/g{h}_k$ (Cohn and Dee, 1989). The software presented here was developed for two sets of boundary conditions. One set was derived from the linearised version of the thermodynamic equation with the assumption that the vertical p velocity, ω, vanishes at the top of the atmosphere and from the assumption that the linearised geometric vertical velocity, $w=\mathrm{d}z/\mathrm{d}t$, vanishes at a constant pressure, p_s, near the surface (Tanaka, 1985). Such boundary conditions may be written as

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_k}{\mathrm{d}p}}+{\displaystyle \frac{p{S}_{\mathrm{0}}}{R{T}_{\mathrm{0}}}}{\mathrm{\Psi }}_k=\mathrm{0}\text{at}p=p_s,\text{(15)}& {\displaystyle }\underset{p\to \mathrm{0}}{lim}{\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_k}{\mathrm{d}p}}=\mathrm{0}.\end{array}$$

Another set of boundary conditions assumes that the atmospheric mass is bounded by an isobaric surface p_s and that the pressure vertical velocity ω must also vanish at p_s. In this case, the lower boundary condition (14) becomes

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_k}{\mathrm{d}p}}=\mathrm{0}\text{at}p=p_s.\end{array}$$

The vertical structure functions (VSFs), Ψ_k(p), can be normalised to satisfy the following orthonormality condition

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{1}}{p_s}}\underset{\mathrm{0}}{\overset{p_s}{\int }}{\mathrm{\Psi }}_i\left(p\right){\mathrm{\Psi }}_j\left(p\right)\mathrm{d}p={\mathit{\delta }}_{ij}.\end{array}$$

Using this orthonormality condition, one can define a vertical transform

$$\begin{array}{}\text{(18)}& f_k={\displaystyle \frac{\mathrm{1}}{p_s}}\underset{\mathrm{0}}{\overset{p_s}{\int }}f\left(p\right){\mathrm{\Psi }}_k\left(p\right)\mathrm{d}p,\end{array}$$

where f(p) is an arbitrary function of pressure.

It can be shown (Cohn and Dee, 1989) that the VSFs form a complete orthonormal basis and the circulation variables $(u,v,\mathit{\varphi }{)}^T$ may be expanded as

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{\left(u,v,\mathit{\varphi }\right)}^T=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\left(u_k,v_k,{\mathit{\varphi }}_k\right)}^T{\mathrm{\Psi }}_k\left(p\right),\text{(20)}& {\displaystyle }& {\displaystyle }\text{where}{\left(u_k,v_k,{\mathit{\varphi }}_k\right)}^T={\displaystyle \frac{\mathrm{1}}{p_s}}\underset{\mathrm{0}}{\overset{p_s}{\int }}{\left(u,v,\mathit{\varphi }\right)}^T{\mathrm{\Psi }}_k\left(p\right)\mathrm{d}p.\end{array}$$

In this software package, the VSE is solved by using the spectral method introduced by Kasahara (1984) (see also Castanheira et al., 1999), which has the advantage over the finite-difference method that the derivatives of the vertical structure function can be calculated by analytical differentiation. The input data are the temperature profile T₀, which is assumed to be defined on N_k pressure levels. Following the computational results of Castanheira et al. (1999), the number J of Legendre polynomials used to approximate Ψ is defined as $J=N_k+\mathrm{20}$. The vertical integration is performed by Gaussian quadrature interpolating T₀ into $G_L=\mathrm{2}J-\mathrm{1}$ Gaussian levels. The maximum number of VSFs saved in the output is equal to the number of pressure levels (N_k).

Figure 1 shows the vertical profiles of temperature T₀ and the base 10 logarithm of the stability parameter S₀. The profiles were computed based on the 3-D temperature field of the ERA-Interim reanalysis (Dee et al., 2011), covering the period from 1 January 1979 to 31 December 2010. The 3-D temperature field was obtained with a $\mathrm{1.5}{}^{\circ }\text{lat.}\times \mathrm{1.5}{}^{\circ }\text{long.}$ grid resolution for the 37 isobaric levels from 1 to 1000 hPa at time intervals of 6 h. The first 12 vertical structures associated with the 12 smallest eigenvalues of ξ_k, which were obtained using boundary conditions (14) and (15), are shown in Fig. 2. The first 12 vertical structures obtained using boundary conditions (16) and (15) are illustrated in Fig. 3. Both sets of 37 VSFs were derived using the temperature profile T₀ and stability parameter S₀ represented in Fig. 1.

Figure 1 (a) Vertical profile of temperature at the ERA-Interim pressure levels (black dots) and interpolated to G_L Gaussian levels (red line). (b) Stability parameter S₀ calculated from T₀ at original pressure levels (black dots) and from interpolated T₀ at Gaussian levels (red line). The stability parameter S₀ was made dimensionless by multiplying Eq. (6) by $p_s^{\mathrm{2}}/g{H}_{*}$ with $H_{*}=\mathrm{8}$ km. The profiles were computed from the 3-D temperature field of the ERA-Interim reanalysis, covering the period from 1 January 1979 to 31 December 2010.

Figure 2 The first 12 vertical structure functions derived using T₀ and S₀ of Fig. 1 with boundary conditions (14) and (15).

Figure 3 The first 12 vertical structure functions derived using T₀ and S₀ of Fig. 1 with boundary conditions (16) and (15).

The baroclinic structures in both Figs. 2 and 3, i.e. the vertical structures with one or several nodes, are close to each other. In Fig. 3 the vertical structure function, Ψ₀, is a pure barotropic vertical structure, i.e. it is a constant, and the projection of an atmospheric variable onto the barotropic mode corresponds exactly to the mass-weighted vertical mean of that variable. The vertical structure function Ψ₀ in Fig. 2 is not strictly constant. Nevertheless, since Ψ₀ of Fig. 2 does not change sign as its counterpart in Fig. 3 and is approximately constant in the troposphere, it is also called as barotropic mode.

2.3 Computation of the horizontal structure functions

Applying the vertical transform (18) to the system of Eqs. (8)–(10), a dimensionless equation is obtained in the following vector form (Tanaka, 1985; Marques and Castanheira, 2012):

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial }{\partial \widehat{t}}}{\mathbf{W}}_k+\mathbf{L}{\mathbf{W}}_k=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }\widehat{t}=\mathrm{2}\mathrm{\Omega }t,\text{(23)}& {\displaystyle }& {\displaystyle }{\mathbf{W}}_k={\mathbf{X}}_k^{-\mathrm{1}}{\left[\begin{array}{c}u\\ v\\ \mathit{\varphi }\end{array}\right]}_k,\text{(24)}& {\displaystyle }& {\displaystyle }\mathbf{L}=\left[\begin{array}{ccc}\mathrm{0}& -\sin \mathit{\theta }& \frac{{\mathit{\alpha }}_k}{\cos \mathit{\theta }}\frac{\partial }{\partial \mathit{\lambda }}\\ \sin \mathit{\theta }& \mathrm{0}& {\mathit{\alpha }}_k\frac{\partial }{\partial \mathit{\theta }}\\ \frac{{\mathit{\alpha }}_k}{\cos \mathit{\theta }}\frac{\partial }{\partial \mathit{\lambda }}& \frac{{\mathit{\alpha }}_k}{\cos \mathit{\theta }}\frac{\partial }{\partial \mathit{\theta }}\left(\cos \mathit{\theta }\left(\cdot \right)\right)& \mathrm{0}\end{array}\right].\end{array}$$

The time was made dimensionless by multiplying with 2 Ω, and the vectors are made dimensionless by multiplying by the inverse of the scaling matrix X_k, which is given by

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathbf{X}}_k=\left[\begin{array}{ccc}\sqrt{g{h}_k}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sqrt{g{h}_k}& \\ \mathrm{0}& \mathrm{0}& g{h}_k\end{array}\right].\end{array}$$

The dimensionless parameter, α_k, in the linear differential matrix operator L is defined as

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathit{\alpha }}_k={\displaystyle \frac{\sqrt{g{h}_k}}{\mathrm{2}\mathrm{\Omega }a}}.\end{array}$$

Equation (21) is a dimensionless form of the linearised shallow water equations.

Since Eq. (21) is a linear system with respect to λ and $\widehat{t}$, the solution W_k can be expressed as zonal waves (Swarztrauber and Kasahara, 1985)

$$\begin{array}{}\text{(27)}& {\displaystyle }{\mathbf{W}}_{nlk}\left(\mathit{\lambda },\mathit{\theta },\widehat{t}\right)={\mathbf{\Theta }}_{nlk}\left(\mathit{\theta }\right)e^{i\left(n\mathit{\lambda }-{\mathit{\sigma }}_{nlk}\widehat{t}\right)},\end{array}$$

where σ_nlk is the dimensionless frequencies for the free waves with zonal wave number n and meridional structures Θ_nlk(θ). The functions ${\mathbf{H}}_{nlk}(\mathit{\lambda },\mathit{\theta })={\mathbf{\Theta }}_{nlk}\left(\mathit{\theta }\right)e^{in\mathit{\lambda }}$ are known as the horizontal structure functions, and the meridional structures Θ_nlk are the Hough vector functions (Longuet-Higgins, 1968; Swarztrauber and Kasahara, 1985). Substituting Eq. (27) into Eq. (21) gives

$$\begin{array}{}\text{(28)}& {\displaystyle }-i{\mathit{\sigma }}_{nlk}{\mathbf{H}}_{nlk}\left(\mathit{\lambda },\mathit{\theta }\right)+\mathbf{L}{\mathbf{H}}_{nlk}\left(\mathit{\lambda },\mathit{\theta }\right)=\mathrm{0},\end{array}$$

and thus the horizontal structure functions are obtained as a free eigenvalue problem. Equation (28) is referred to as the horizontal structure equation (HSE).

For n>0, there are two kinds of motion (i.e. solutions) with distinct frequencies for system (21). The first kind describes high-frequency eastward- and westward-propagating gravity inertia waves, and the second kind describes low-frequency westward-propagating rotational waves of the Rossby–Haurwitz type (Longuet-Higgins, 1968; Kasahara and Puri, 1981; Swarztrauber and Kasahara, 1985; Žagar et al., 2015b).

It can be shown that for real α_k all of the frequencies σ of L are real and that any modes H corresponding to distinct frequencies are orthogonal (e.g. Swarztrauber and Kasahara, 1985; Žagar et al., 2015b). For n>0, the frequencies are distinct and thus the modes are orthogonal. However, for zonal wavenumber n=0, the rotational modes are not necessarily orthogonal because their frequencies are all zero. Nevertheless, it is possible to derive an orthogonal set of rotational modes for n=0, which are also orthogonal to the modes for n>0 (Kasahara and Puri, 1981; Shigehisa, 1983; Swarztrauber and Kasahara, 1985). Therefore, all horizontal structure functions H satisfy the orthonormal condition given by

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{\mathbf{H}}_{nlk}^{\ast }\cdot {\mathbf{H}}_{n^{\prime }{l}^{\prime }k}\cos \mathit{\theta }\mathrm{d}\mathit{\lambda }\mathrm{d}\mathit{\theta }={\mathit{\delta }}_{n{n}^{\prime }}{\mathit{\delta }}_{l{l}^{\prime }},\end{array}$$

where the superscript ∗ denotes a conjugate transpose and the right-hand side is unity if $n=n^{\prime }$ and $l=l^{\prime }$ and zero otherwise.

The Hough vector function, Θ_nlk(θ), has three components, zonal velocity, meridional velocity, and height, and is given by

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathbf{\Theta }}_{nlk}\left(\mathit{\theta }\right)=\left[\begin{array}{c}{\mathrm{U}}_{nlk}\left(\mathit{\theta }\right)\\ i{\mathrm{V}}_{nlk}\left(\mathit{\theta }\right)\\ {\mathrm{\Phi }}_{nlk}\left(\mathit{\theta }\right)\end{array}\right],\end{array}$$

where the factor i in front of V_nlk is introduced to account for a phase shift of π∕2 (Kasahara, 1977; Swarztrauber and Kasahara, 1985). Substituting Eq. (30) into ${\mathbf{H}}_{nlk}(\mathit{\lambda },\mathit{\theta })={\mathbf{\Theta }}_{nlk}\left(\mathit{\theta }\right)e^{in\mathit{\lambda }}$ in Eq. (29), one can see that the Hough vector functions satisfy the orthonormal condition given by

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& \underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}{\mathbf{\Theta }}_{nlk}^{\ast }\cdot {\mathbf{\Theta }}_{n{l}^{\prime }k}\cos \mathit{\theta }\mathrm{d}\mathit{\theta }=\\ & \underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}({\mathrm{U}}_{nlk}{\mathrm{U}}_{n{l}^{\prime }k}+{\mathrm{V}}_{nlk}{\mathrm{V}}_{n{l}^{\prime }k}+{\mathrm{\Phi }}_{nlk}{\mathrm{\Phi }}_{n{l}^{\prime }k})\cos \mathit{\theta }\mathrm{d}\mathit{\theta }\\ & ={\mathit{\delta }}_{l{l}^{\prime }}.\end{array}\end{array}$$

Using orthonormal condition (29), a complete set of Fourier–Hough transforms may be constructed as

$$\begin{array}{}\text{(32)}& {\displaystyle }& {\displaystyle }{\mathbf{W}}_k\left(\mathit{\lambda },\mathit{\theta },t\right)=\sum _{l=\mathrm{0}}^{\mathrm{\infty }}\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{w}_{nlk}\left(t\right){\mathbf{H}}_{nlk}\left(\mathit{\lambda },\mathit{\theta }\right),\text{(33)}& {\displaystyle }& {\displaystyle }{w}_{nlk}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{\mathbf{H}}_{nlk}^{\ast }\left(\mathit{\lambda },\mathit{\theta }\right)\cdot {\mathbf{W}}_k\left(\mathit{\lambda },\mathit{\theta },t\right)\cos \mathit{\theta }\mathrm{d}\mathit{\lambda }\mathrm{d}\mathit{\theta }.\end{array}$$

In order to distinguish each wave type, i.e. the westward-propagating Rossby wave and the westward- and eastward-propagating gravity inertia waves, the meridional index l is defined as a sequence of three distinct modes, l_r, l_w and l_e, respectively (e.g. Kasahara and Puri, 1981; Marques and Castanheira, 2012).

2.4 Solutions of the horizontal structure equation

As shown in Swarztrauber and Kasahara (1985), the meridional modal functions Θ_nlk(θ) for the eastward and westward gravity inertia modes are symmetric (antisymmetric) with respect to the Equator for the even (odd) meridional index l. For the rotational modes, the functions Θ_nlk(θ) are symmetric (antisymmetric) for odd (even) l. For each vertical mode k, the frequencies σ_nlk and the corresponding meridional modal functions Θ_nlk(θ) are computed for a range of zonal wavenumbers n and meridional indices l. The computational method was developed by Swarztrauber and Kasahara (1985) and consists of expanding the eigensolutions of Eq. (28) in terms of spherical vector harmonics. The modes for n=0 are determined following the approach suggested by Shigehisa (1983), as described in Swarztrauber and Kasahara (1985).

Figure 4 shows the dimensionless frequencies σ of westward-propagating gravity inertia and rotational waves (Fig. 4a) and eastward-propagating gravity inertia waves (Fig. 4b), for the case n=1. This figure corresponds to Fig. 2 in Swarztrauber and Kasahara (1985) because we have used the same values of α as they did to compute the frequencies σ. The α used here corresponds to $\mathrm{1}/\sqrt{\mathit{ϵ}}$ of Swarztrauber and Kasahara (1985). The frequencies are plotted as a function of α⁻¹ for 10 meridional indices ($l=\mathrm{0},\mathrm{\dots },\mathrm{9}$), 5 symmetric (continuous lines), and 5 antisymmetric (dashed lines). The logarithm scale was used on both axes and the y axis for the westward-propagating waves has been reversed. The rotational mode l_r=0 is referred to as a mixed Rossby–gravity wave because it behaves like a rotational mode for large values of α but behaves like a gravity mode for small values of α. The symmetric mode l_e=0 is referred to as a Kelvin wave.

Figure 4 Dimensionless frequencies σ for westward-propagating gravity waves and rotational waves (a) and for eastward-propagating gravity waves (b), plotted as a function of α⁻¹ and for zonal wavenumber n=1.

The meaning of eastward and westward propagations is lost in the case of n=0. However, since the frequencies of gravity inertia motion appear as pairs of positive and negative values with the same magnitudes, we also use the term eastward (westward) to indicate positive (negative) frequency, as adopted by Swarztrauber and Kasahara (1985). The frequencies of the rotational motion, along with the frequencies of the gravity motion corresponding to the lowest meridional index (l=0), are zero for n=0 (Swarztrauber and Kasahara, 1985). Therefore, instead of saving zeros for all of the rotational frequencies and for the gravity frequencies corresponding to the lowest meridional index, the software computes the asymptotic rate (σ_a) at which those frequencies go to zero (see Eqs. 4.14 and 4.18 in Swarztrauber and Kasahara, 1985, for the definition and computational formula of σ_a).

Figure 5a shows the curves of asymptotic dimensionless frequencies σ_a in the case of n→0. The frequencies are plotted as a function of α⁻¹ for 10 meridional indices, using a log–log scale with the y axis reversed. In this case, the lowest symmetric mode is identified as $l_{\mathrm{r}}=-\mathrm{1}$ because the symmetric equations begin at index −1 (see Swarztrauber and Kasahara, 1985, for details;). The values of σ_a are all negative, except for the mode $l_{\mathrm{r}}=-\mathrm{1}$, which is similar to the eastward-propagating Kelvin wave in Fig. 4b. For this reason, the $l_{\mathrm{r}}=-\mathrm{1}$ mode is saved as the eastward gravity mode l_e=0 as in Swarztrauber and Kasahara (1985). The curves for the modes $l_{\mathrm{r}}=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{9}$ behave analogously to those of rotational modes in Fig. 4a. The absolute values of dimensionless frequencies σ in the case of n=0 of the gravity waves, are illustrated in Fig. 5b for $l=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{9}$. Figure 5a and b correspond to Figs. 1 and 3 of Swarztrauber and Kasahara (1985), respectively.

Figure 5 Asymptotic dimensionless frequencies σ_a in the case of n→0 (a) and dimensionless frequencies $\left|\mathit{\sigma }\right|$ for gravity waves and zonal wavenumber n=0 (b), plotted as a function of α⁻¹.

As mentioned above, the horizontal modes include Kelvin, westward and eastward inertio-gravity (WIG and EIG, respectively), Rossby, and mixed Rossby–gravity (MRG) waves. As an illustration of the software outputs, the horizontal structure functions were calculated for each equivalent height computed in Sect. 2.2 using boundary conditions (14) and (15) (the first 12 VSFs are represented in Fig. 2), and for 43 wavenumbers ($n=\mathrm{0},\mathrm{\dots },\mathrm{42}$) and 80 meridional modes (20 WIG, 20 EIG and 40 Rossby modes). The meridional profiles of the Hough functions corresponding to two Rossby modes l_r=1 and l_r=2 are shown in Fig. 6 for the barotropic mode (k=0, with h₀≃9.8 km) and for two zonal wavenumbers (n=0 and n=10). Figure 7 shows the meridional profiles corresponding to the Kelvin mode (l_e=0) and to the Rossby mode (l_r=1) for two baroclinic modes (k=1, with h₁≃6.1 km and k=10, with h₁₀≃65 m) and for two zonal wavenumbers (n=0 and n=10). By increasing the vertical index k (i.e. reducing the equivalent height h_k) or the zonal wavenumber n for the same vertical index k, the horizontal structure functions become more confined around the Equator (Longuet-Higgins, 1968; Cohn and Dee, 1989; Žagar et al., 2015b).

Figure 6 Meridional structures of the barotropic (k=0, h₀≃9.8 km) Hough harmonics for two Rossby modes, l_r=1 (a, b) and l_r=2 (c, d), and two zonal wavenumbers, n=0 (a, c) and n=10 (b, d). The equivalent height h₀ was obtained by solving the vertical structure equation with boundary conditions (14) and (15).

Figure 7 Meridional structures of two baroclinic Hough harmonics (k=1, h₁≃6.1 km and k=10, h₁₀≃65 m) for the Kelvin mode l_e=0 (left column), the Rossby mode, l_r=1 (right column), and for two zonal wavenumbers, n=0 and n=10. The equivalent heights h_k were obtained by solving the vertical structure equation with boundary conditions (14) and (15).

2.4.1 The Haurwitz waves

The Hough vector functions define modes in the transformed variables as given by Eq. (23). If the equivalent height h_k is infinite, then transform (23) is no longer applicable. This is the case for the barotropic mode when the lower boundary condition (16) is used. However, an infinite equivalent height means that the separation constant in Eq. (12) is null, and the horizontal motion is non-divergent, and there must be no vertical dependence in the linear system (8)–(10). Consistently, the respective vertical structure is a constant. In this case, the horizontal modes are computed in the untransformed variables u,v, and ϕ as in Swarztrauber and Kasahara (1985), and are called the Rossby–Haurwitz waves. Additionally, because the linear system (8)–(10) is non-divergent, there are no WIG and EIG modes (and also no MRG mode for the zonal mean component n=0). Figure 8 shows the same meridional profiles as in Fig. 6 but for the Haurwitz modes.

Figure 8 The same as in Fig. 6 but for the Haurwitz modes (h₀=∞).

2.5 3-D normal mode functions

The three-dimensional (3-D) normal mode functions (NMFs) are obtained by the product of the vertical normal modes Ψ_k(p) (the eigensolutions of the VSE) and the horizontal normal modes H_nlk(λ,θ) – the eigensolutions of the HSE (e.g. Kasahara, 1976; Kasahara and Puri, 1981). The 3-D NMFs, denoted by ${\mathbf{\Pi }}_{nlk}\left(\mathit{\lambda },\mathit{\theta },p\right)$, form a complete orthogonal basis, therefore allowing us to expand the horizontal wind and the geopotential fields of the global atmosphere (Kasahara and Puri, 1981; Tanaka, 1985; Daley, 1991):

$$\begin{array}{}\text{(34)}& {\left(u,v,\mathit{\varphi }\right)}^T=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\sum _{l=\mathrm{0}}^{\mathrm{\infty }}\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{w}_{nlk}\left(t\right){\mathbf{X}}_k{\mathbf{\Pi }}_{nlk}\left(\mathit{\lambda },\mathit{\theta },p\right),\end{array}$$

where

$$\begin{array}{}\text{(35)}& {\mathbf{\Pi }}_{nlk}\left(\mathit{\lambda },\mathit{\theta },p\right)={\mathrm{\Psi }}_k\left(p\right){\mathbf{\Theta }}_{nlk}\left(\mathit{\theta }\right)e^{in\mathit{\lambda }}.\end{array}$$

The expansion coefficients w_nlk are obtained by means of the vertical projection onto the vertical structure functions, followed by the horizontal projection onto the horizontal structure functions:

$$\begin{array}{}\text{(36)}& \begin{array}{rl}& w_{nlk}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{p}_s}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{p_s}{\int }}{\mathbf{\Pi }}_{nlk}^{\ast }\left(\mathit{\lambda },\mathit{\theta },p\right)\\ & \cdot {\mathbf{X}}_k^{-\mathrm{1}}\mathbf{W}\cos \mathit{\theta }\mathrm{d}p\mathrm{d}\mathit{\lambda }\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

with $\mathbf{W}={\left[u,v,\mathit{\varphi }\right]}^T$.

2.6 3-D normal mode energetics scheme

Using the 3-D NMFs of the linearised primitive Eqs. (8)–(10) as a basis to expand the global circulation field, Tanaka (1985) (see also Tanaka and Kung, 1988) developed a 3-D normal mode energetics (NME) scheme, which combines three one-dimensional spectral energetics in domains of the zonal wavenumber, n; meridional mode number, l; and vertical mode number, k. The 3-D NME scheme therefore complements the standard energetics in the zonal wavenumber domain of Saltzman (1957), since it can diagnose the 3-D spectral distribution of energy and energy interactions, the energetics characteristics of Rossby waves and gravity inertia waves, and the energy interaction between the barotropic and baroclinic modes (Tanaka and Kung, 1988).

Derivation of the 3-D NME may be summarised as follows.

Applying the vertical transform (18) to Eqs. (1), (2), and (7), a dimensionless equation is obtained in the following vector form:

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\partial }{\partial \widehat{t}}}{\mathbf{W}}_k+\mathbf{L}{\mathbf{W}}_k={\mathbf{I}}_k+{\mathbf{J}}_k+{\mathbf{S}}_k,\end{array}$$

where subscript k denotes the kth component of the vertical transform and $\widehat{t}$, W_k, and L are given by Eqs. (22), (23), and (24), respectively. The dimensionless non-linear term vectors for the wind and mass fields, I_k and J_k, and the energy source or sink term due to diabatic heating and dissipation, S_k, are given by

$$\begin{array}{}\text{(38)}& {\displaystyle }& {\displaystyle }{\mathbf{I}}_k={\mathbf{Y}}_k^{-\mathrm{1}}{\left[\begin{array}{c}-\mathit{V}\cdot \mathrm{\nabla }u-\mathit{\omega }\frac{\partial u}{\partial p}+\frac{\tan \mathit{\theta }}{a}uv\\ -\mathit{V}\cdot \mathrm{\nabla }v-\mathit{\omega }\frac{\partial v}{\partial p}-\frac{\tan \mathit{\theta }}{a}{u}^{\mathrm{2}}\\ \mathrm{0}\end{array}\right]}_k,\text{(39)}& {\displaystyle }& {\displaystyle }{\mathbf{J}}_k={\mathbf{Y}}_k^{-\mathrm{1}}{\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \frac{\partial }{\partial p}\left[\frac{R}{p{S}_{\mathrm{0}}}\left(-\mathit{V}\cdot \mathrm{\nabla }T-\mathit{\omega }\frac{\partial T}{\partial p}\right)\right]\end{array}\right]}_k,\text{(40)}& {\displaystyle }& {\displaystyle }{\mathbf{S}}_k={\mathbf{Y}}_k^{-\mathrm{1}}{\left[\begin{array}{c}{F}_u\\ F_v\\ \frac{\partial }{\partial p}\left(\frac{qR}{p{c}_p{S}_{\mathrm{0}}}\right)\end{array}\right]}_k.\end{array}$$

These vectors are made dimensionless using the inverse of scaling matrix Y_k, which is given by

$$\begin{array}{}\text{(41)}& {\displaystyle }{\mathbf{Y}}_k=\left[\begin{array}{ccc}\mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}& \\ \mathrm{0}& \mathrm{0}& \mathrm{2}\mathrm{\Omega }\end{array}\right].\end{array}$$

Applying the Fourier–Hough transform (33) to (37), yields

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}\widehat{t}}}{w}_{nlk}+i{\mathit{\sigma }}_{nlk}{w}_{nlk}={\mathit{ı}}_{nlk}+j_{nlk}+s_{nlk},\end{array}$$

where the complex variables w_nlk, ı_nlk, j_nlk, and s_nlk are the Fourier–Hough transforms of the vector variables W_k, I_k, J_k, and S_k, respectively.

The kinetic energy, K, and available potential energy, A, are defined, respectively, by

$$\begin{array}{}\text{(43)}& {\displaystyle }K={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(u^{\mathrm{2}}+v^{\mathrm{2}}\right),\text{(44)}& {\displaystyle }A={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{S}_{\mathrm{0}}}}{\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial p}}\right)}^{\mathrm{2}},\end{array}$$

and the total energy, $E=K+A$, is a conserved quantity of the full non-linear Eqs. (1), (2), and (7) under frictionless and adiabatic conditions (see Tanaka and Kung, 1988).

In the transformed 3-D normal mode space, the total energy E_nlk associated with each mode (nlk) is

$$\begin{array}{}\text{(45)}& {\displaystyle }{E}_{nlk}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{ϵ}}_n}}{p}_s{h}_k\mathrm{|}{w}_{nlk}{\mathrm{|}}^{\mathrm{2}},\end{array}$$

where ϵ_n=2 for the zonal (n=0) modes, and ϵ_n=1 for the eddy (n>0) modes.

Substituting Eq. (42) into the time derivatives of Eq. (45), the energy balance equations for the 3-D normal modes are finally obtained as

$$\begin{array}{}\text{(46)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{E}_{nlk}=I_{nlk}+J_{nlk}+S_{nlk},\end{array}$$

where

$$\begin{array}{}\text{(47)}& {\displaystyle }{I}_{nlk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\left[w_{nlk}^{\ast }{\mathit{ı}}_{nlk}+w_{nlk}{\mathit{ı}}_{nlk}^{\ast }\right],\text{(48)}& {\displaystyle }{J}_{nlk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\left[w_{nlk}^{\ast }{j}_{nlk}+w_{nlk}{j}_{nlk}^{\ast }\right],\text{(49)}& {\displaystyle }{S}_{nlk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\left[w_{nlk}^{\ast }{s}_{nlk}+w_{nlk}{s}_{nlk}^{\ast }\right].\end{array}$$

The terms I_nlk,J_nlk, and S_nlk that contribute to the time change of E_nlk are, respectively, associated with non-linear interactions of kinetic and available potential energies and with an energy source or sink due to diabatic heating and dissipation. The spectra of these terms can therefore be analysed separately for the zonal mean and eddy components, the barotropic and baroclinic modes, and the Rossby and gravity waves (Tanaka and Kung, 1988).

As another illustration of the capabilities of the current software, in Fig. 9 the spectra for total energy E (Fig. 9a, b, c), interactions of available potential energy J (Fig. 9d, e, f), and interactions of kinetic energy I are shown (Fig. 9g, h, i). Figure 9a, d, and g contain the wavenumber spectra, whereas Fig. 9b, e, and h and Fig. 9c, f, and i show the vertical spectra for the zonal mean and eddy components, respectively. In each panel the total spectra are represented, along with the contributions of the Rossby and gravity components. These spectra were computed using 34 years of ERA-Interim reanalysis (Dee et al., 2011) data covering the period from 1 January 1979 to 31 December 2012. The horizontal wind (u,v), pressure velocity (ω), temperature, and geopotential fields were obtained with $\mathrm{1.5}{}^{\circ }\text{lat.}\times \mathrm{1.5}{}^{\circ }\text{long.}$ grid resolution for the 37 isobaric levels from 1000 to 1 hPa at time intervals of 6 h. The solutions of the VSE (13), with boundary conditions (16) and (15), were approximated using 57 Legendre polynomials and 37 vertical modes were retained. The zonal wavenumber has been truncated at n=42. The spectra were computed at each time step (6-hourly) and then averaged over the 30-year period.

A detailed analysis of the spectra in Fig. 9 is beyond the scope of this study, but some of its characteristics are mentioned. The wavenumber energy spectrum of Rossby modes approximately follows the −3 power law over the synoptic and mesoscale region, while that of the gravity modes has approximately a $-\mathrm{5}/\mathrm{3}$ slope, which is in line with the literature (e.g. Charney, 1971; Nastrom and Gage, 1985; Tanaka, 1985; Tanaka et al., 1986; Terasaki and Tanaka, 2007). The zonal mean available potential energy is transferred into eddy available potential energy essentially due to the Rossby waves, with a small contribution by planetary-scale gravity waves, as indicated by the positive values in the wavenumber spectra of J. On the other hand, the wavenumber spectra of I indicate that the Rossby waves transfer energy out of the eddy kinetic energy reservoir, whereas the smaller quantity of energy transferred by the gravity waves is in the opposite direction. The vertical spectra for I shows that the eddy kinetic energy contained in the Rossby baroclinic modes is transferred into the Rossby barotropic modes of both zonal mean and eddy components, and also into baroclinic zonal mean kinetic energy. Again, it is seen that the transfer of kinetic energy due to the gravity waves is in the opposite direction.

Figure 9 Spectra of total energy E (a, b, c) and of non-linear interactions of available potential J (d, e, f) and kinetic I (g, h, i) energies, in the zonal wavenumber (a, d, g) and the vertical index, separately, for the zonal (b, e, h) (0) and eddy (c, f, i) (n) components. In each panel the contributions of the Rossby (R) and gravity (G) components to the total (T) spectra of energy or energy interactions are represented.

A disadvantage of the 3-D NME is that one cannot separate the available potential and kinetic energies, only the total energy (E_nlk) associated with each mode can be calculated. Consequently, the conversion rate of available potential energy into kinetic energy also cannot be accessed. However, the separation of the available potential and kinetic energies can be performed if one uses only the expansions in the vertical and wavenumber domains. This reasoning led Marques and Castanheira (2012) to present a normal mode energetics formulation that performs an explicit evaluation of the available potential and kinetic energies, as well as the conversion rates between them. In addition, the generation and dissipation rates and the non-linear interactions of each energy form can be performed in both the zonal wavenumber and vertical mode domains. Using the vertical normal mode basis that results from applying lower boundary condition (14) and (15), and then performing the vertical and the zonal wavenumber (Fourier) expansions, the set of balance equations for the kinetic energy (K) and available potential energy (A) is given by (see Marques and Castanheira, 2012, for details)

$$\begin{array}{}\text{(50)}& {\displaystyle \frac{\mathrm{d}{K}_{nk}}{\mathrm{d}t}}=C_{nk}+I_{nk}-D_{nk},\text{(51)}& {\displaystyle \frac{\mathrm{d}{A}_{nk}}{\mathrm{d}t}}=-C_{nk}+J_{nk}+G_{nk},\end{array}$$

where

$$\begin{array}{}\text{(52)}& {\displaystyle }& {\displaystyle }{K}_{nk}={\displaystyle \frac{p_s{h}_k}{\mathrm{2}{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left({\widehat{u}}_{nk}^{\ast }{\widehat{u}}_{nk}+{\widehat{v}}_{nk}^{\ast }{\widehat{v}}_{nk}\right)\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\text{(53)}& {\displaystyle }& {\displaystyle }{A}_{nk}={\displaystyle \frac{p_s{h}_k}{\mathrm{2}{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left({\widehat{\mathit{\varphi }}}_{nk}^{\ast }{\widehat{\mathit{\varphi }}}_{nk}\right)\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}$$

$$\begin{array}{}\text{(54)}& \begin{array}{rl}& C_{nk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left\{{\displaystyle \frac{in{\mathit{\alpha }}_k}{\cos \mathit{\theta }}}\left({\widehat{u}}_{nk}{\widehat{\mathit{\varphi }}}_{nk}^{\ast }-{\widehat{u}}_{nk}^{\ast }{\widehat{\mathit{\varphi }}}_{nk}\right)\right.\\ & -\left.\phantom{{\displaystyle \frac{in{\mathit{\alpha }}_k}{\cos \mathit{\theta }}}}{\mathit{\alpha }}_k\left({\widehat{v}}_{nk}^{\ast }{\displaystyle \frac{\partial {\widehat{\mathit{\varphi }}}_{nk}}{\partial \mathit{\theta }}}+{\widehat{v}}_{nk}{\displaystyle \frac{\partial {\widehat{\mathit{\varphi }}}_{nk}^{\ast }}{\partial \mathit{\theta }}}\right)\right\}\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

$$\begin{array}{}\text{(55)}& {\displaystyle }\begin{array}{rl}& I_{nk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left\{\left({\widehat{u}}_{nk}^{\ast }{\left({\widehat{I}}_{\mathrm{1}}\right)}_{nk}+{\widehat{u}}_{nk}{\left({\widehat{I}}_{\mathrm{1}}\right)}_{nk}^{\ast }\right)\right.\\ & +\left.\phantom{{\left({\widehat{I}}_{\mathrm{1}}\right)}_{nk}^{\ast }}\left({\widehat{v}}_{nk}^{\ast }{\left({\widehat{I}}_{\mathrm{2}}\right)}_{nk}+{\widehat{v}}_{nk}{\left({\widehat{I}}_{\mathrm{2}}\right)}_{nk}^{\ast }\right)\right\}\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}\text{(56)}& {\displaystyle }\begin{array}{rl}& J_{nk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left({\widehat{\mathit{\varphi }}}_{nk}^{\ast }{\left({\widehat{J}}_{\mathrm{3}}\right)}_{nk}\right.\\ & \left.+{\widehat{\mathit{\varphi }}}_{nk}{\left({\widehat{J}}_{\mathrm{3}}\right)}_{nk}^{\ast }\right)\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

$$\begin{array}{}\text{(57)}& {\displaystyle }\begin{array}{rl}& D_{nk}=-{\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left\{\left({\widehat{u}}_{nk}^{\ast }{\left({\widehat{F}}_u\right)}_{nk}+{\widehat{u}}_{nk}{\left({\widehat{F}}_u\right)}_{nk}^{\ast }\right)\right.\\ & +\left.\phantom{\left({\widehat{F}}_u\right)}\left({\widehat{v}}_{nk}^{\ast }{\left({\widehat{F}}_v\right)}_{nk}+{\widehat{v}}_{nk}{\left({\widehat{F}}_v\right)}_{nk}^{\ast }\right)\right\}\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}\text{(58)}& {\displaystyle }\begin{array}{rl}& G_{nk}={\displaystyle \frac{p_s\mathrm{\Omega }{h}_k}{{\mathit{ϵ}}_n}}\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\left({\widehat{\mathit{\varphi }}}_{nk}^{\ast }{\left({\widehat{N}}_{\mathrm{3}}\right)}_{nk}\right.\\ & \left.+{\widehat{\mathit{\varphi }}}_{nk}{\left({\widehat{N}}_{\mathrm{3}}\right)}_{nk}^{\ast }\right)\cos \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

and

$$\begin{array}{}\text{(59)}& {\displaystyle }& {\displaystyle }{\left({\widehat{I}}_{\mathrm{1}}\right)}_{nk}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}}}{\left[-\mathit{V}\cdot \mathrm{\nabla }u-\mathit{\omega }{\displaystyle \frac{\partial u}{\partial p}}+{\displaystyle \frac{\tan \mathit{\theta }}{a}}uv\right]}_{nk},\text{(60)}& {\displaystyle }& {\displaystyle }{\left({\widehat{I}}_{\mathrm{2}}\right)}_{nk}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}}}{\left[-\mathit{V}\cdot \mathrm{\nabla }v-\mathit{\omega }{\displaystyle \frac{\partial v}{\partial p}}-{\displaystyle \frac{\tan \mathit{\theta }}{a}}{u}^{\mathrm{2}}\right]}_{nk},\end{array}$$

$$\begin{array}{}\text{(61)}& {\displaystyle }& {\displaystyle }({\widehat{J}}_{\mathrm{3}}{)}_{nk}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }}}{\left({\displaystyle \frac{\partial }{\partial p}}\left[{\displaystyle \frac{R}{p{S}_{\mathrm{0}}}}\left(-\mathit{V}\cdot \mathrm{\nabla }T-\mathit{\omega }{\displaystyle \frac{\partial T}{\partial p}}\right)\right]\right)}_{nk},\text{(62)}& {\displaystyle }& {\displaystyle }({\widehat{N}}_{\mathrm{3}}{)}_{nk}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }}}{\left[{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{qR}{p{c}_p{S}_{\mathrm{0}}}}\right)\right]}_{nk},\end{array}$$

$$\begin{array}{}\text{(63)}& {\displaystyle }& {\displaystyle }{\left({\widehat{F}}_u\right)}_{nk}={\displaystyle \frac{{\left(F_u\right)}_{nk}}{\mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}}},\text{(64)}& {\displaystyle }& {\displaystyle }{\left({\widehat{F}}_v\right)}_{nk}={\displaystyle \frac{{\left(F_v\right)}_{nk}}{\mathrm{2}\mathrm{\Omega }\sqrt{g{h}_k}}}.\end{array}$$

Figure 10 Extended energy cycle diagram describing the flow of energy among the zonal mean and eddy components and among the barotropic and baroclinic components for ERA-Interim in DJF (top values) and JJA (bottom values) climates. Units are 10⁵ J m⁻² for energy levels and W m⁻² for conversion and transfer rates and generation and dissipation rates.

Figure 11 The same as in Fig. 10 but replacing lower boundary condition (14) with (16).

The term C_nk represents the conversion of available potential energy into kinetic energy, whereas I_nk (J_nk) represents interactions or exchanges of kinetic (available potential) energy between different scales or types of motion. Finally, the dissipation of kinetic energy and the generation of available potential energy, due to diabatic processes, are represented by terms D_nk and G_nk, respectively, and may be computed as residuals from balance Eqs. (50) and (51). The calculated energetic terms may then be used to construct an extended energy cycle diagram, as in Fig. 7 of Marques and Castanheira (2012) (see also Fig. 6 in Marques et al., 2014). Such an extended energy cycle diagram is illustrated in Fig. 10, in which the boxes represent the levels of energy and the arrows the energy generation and dissipation rates and the energy conversion and transfer rates. The estimates in the diagram were based on the same 34 years of ERA-Interim reanalysis mentioned above for the 3-D spectra of E, J, and I. The solutions of the VSE (13) were obtained using boundary conditions (14) and (15). The energetic terms were computed at each time step (6-hourly) and then averaged over the 34-year period to obtain the corresponding mean values for the northern winter (DJF) and summer (JJA) seasons, which are represented as the top (black) and bottom (red) values in Fig. 10.

The energy cycle diagram of Fig. 10 describes the flow of energy among the zonal mean and eddy components, and also among the barotropic and baroclinic components. All terms in the diagram are decomposed into the zonal mean (n=0) and eddy (n≥1) components, which are denoted, respectively, by subscripts Z and E. Each one of these components is also decomposed into the barotropic (k=0) and baroclinic (k≥1) components by using the extra subscripts B and b, respectively. The two exceptions are the barotropic–baroclinic interactions of available potential energy, denoted by J_Bb, and the baroclinic–barotropic interactions of kinetic energy, which is represented as I_bB. These values represent the balances of the flows of kinetic energy and available potential energy between the barotropic and baroclinic components and both zonal mean and eddies, implying that the energy flows represented by the dotted lines in the diagram cannot by quantified individually. These flows are associated with eddy generation by barotropic instability and the barotropic decay of baroclinic eddies (see Marques and Castanheira, 2012, and references therein, for details).

As an alternative, Marques and Castanheira (2017) and Castanheira and Marques (2019) presented a similar NME scheme, but using the vertical normal mode basis that results from applying lower boundary condition (16) in place of (14). The main difference between the NME schemes presented in Marques and Castanheira (2012) and in Marques and Castanheira (2017) is that in the latter, the first vertical structure is strictly constant, which implies that its vertical derivative is null and there is no available potential energy associated with the barotropic mode (k=0). Therefore, the terms A_n0, C_n0, J_n0, and G_n0 are null in this variant of the energetics scheme. All baroclinic terms (k≥1) have identical expressions in both NME schemes. The terms associated with the kinetic energy of the barotropic component are calculated formally with identical expressions in both NME schemes. However, a constant h₀ cannot be derived from the eigenvalues of the VSE (13). In this case the smallest eigenvalue is null and corresponds to the separation constant in Eq. (12) being null. For such a case, it is not possible to define an equivalent height as it was in Eq. (12). In order to keep the expressions for the barotropic and baroclinic terms formally identical in both NME schemes, the software fixes h₀=1m, when the boundary condition (16) is used. Calculating the terms given by Eqs. (52)–(58) in this way,  (50) and (51) may also be used to construct an extended energy cycle diagram as illustrated in Fig. 11 (see also Fig. A2 in Marques and Castanheira, 2017). In this case, there is no barotropic branch in the available potential energy side, as compared to extended energy cycle diagram of Fig. 10. The estimates in this diagram were based on the same data as those of Fig. 10. The solutions of the VSE were also approximated using 57 Legendre polynomials, retaining 37 vertical modes. The same zonal wavenumbers as those for the estimates in Fig. 10 were also considered.

Although boundary condition (16) may be seen as somewhat unrealistic, it presents some useful features. First, using the boundary conditions (15) and (16), it can easily be shown that the total energy is exactly conserved by full non-linear Eqs. (1)–(5), in the case of frictionless adiabatic motions (see Tanaka and Kung, 1988). Moreover, an extra surface term, which was interpreted by Tanaka and Kung (1988) as a geopotential flux across the lower boundary, does not appear in the available potential energy as defined by Eq. (44), when the boundary condition (16) is used. Another feature of (16) is that the barotropic mode represents the mass-weighted average of the circulation field.

The use of the boundary condition (14) may seem more realistic, but in order to guarantee the conservation of energy for frictionless adiabatic motions it is necessary to impose a non-slipping condition $(u,v{)}_{\mathrm{s}}=\mathrm{0}$ at p=p_s, though this may seem contradictory with the absence of friction.

3 Conclusions

The software presented here consists of tools or functions for specific tasks and includes a Python and a MATLAB version. The Python version requires NumPy and SciPy modules and works with both Python 2 and 3. The netCDF4 library is recommended. The input–output data format for each function may be either the native format of the chosen version, i.e. “.mat” for MATLAB and “.npz” for Python, or netCDF for both versions.

The tasks to be achieved with the functions in the software include the solution of the VSE (13), the solution of HSE (28), and the computation of complex expansion coefficients w_nlk, ı_nlk, and j_nlk, which are the vertical Fourier–Hough transforms of the dependent variable vector $[u,v,\mathit{\varphi }{]}^T$, of the non-linear term vector due to wind field, and of the non-linear term due to mass field, respectively. In addition, the software also permits the computation of vertical transforms of the terms in the normal mode energetics formulation presented by Marques and Castanheira (2012) or its alternative presented by Marques and Castanheira (2017) and Castanheira and Marques (2019). These terms constitute the set of balance Eqs. (50) and (51), and its expressions are given in Eqs. (52)–(64). With these vertical transforms the user may compute the global energy cycle in the wavenumber and vertical domains (as in Marques and Castanheira, 2012, 2017; Castanheira and Marques, 2019), with the dissipation and generation terms computed as residuals from Eqs. (50) and (51).

Other applications of this software include the spatial-scale filtering of atmospheric motion and filtering of balanced motion and mass unbalanced motions by retaining an appropriate subset of terms in the expansion (34) (e.g. Castanheira and Marques, 2015; Žagar et al., 2020). Furthermore, the solutions of HSE, i.e. the Hough vector functions (Θ_nlk(θ)), can be used as meridional basis functions onto which dynamical data may be projected. This can be an alternative tool for investigating studies on tropical convection as in Yang et al. (2003) and Gehne and Kleeman (2012), which have used parabolic cylinder functions as meridional basis functions. All of these features are useful both in data analysis and in assessment of general circulation atmospheric models (e.g. Žagar et al., 2020).

Appendix A: Software functions

• Vertical structure.

  The vertical structure equation (VSE), (13), is solved with a function named vertical_structure that uses the spectral method introduced by Kasahara (1984). The function requires the mean vertical profile of temperature T₀, defined on pressure levels, and the pressure levels as inputs, returning the vertical structure functions (Ψ_k(p)) and equivalent heights (h_k). The VSE may be solved using either one of the two pairs of boundary conditions, as described in Sect. 2.2. These pairs of boundary conditions are controlled by the input argument ws0. When ws0=False (the default) the VSE is solved using boundary conditions (14) and (15), whereas when ws0=True the VSE is solved using boundary conditions (16) and (15). In the latter case the first vertical structure is strictly constant and the associated equivalent height is infinite.

• Horizontal structure.

  The solution of HSE (21) is obtained with a function called hough_functions. It uses the computational method developed by Swarztrauber and Kasahara (1985), which consists of expanding the eigensolutions of Eq. (28) in terms of spherical vector harmonics. The function inputs are the equivalent heights (h_k) and the user-defined parameters for the maximum zonal wavenumber, the total number of Rossby modes and half the number of gravity modes used in the expansion, along with the latitude points, which may be linear (the default) or Gaussian. It returns the Hough vector functions (Θ_nlk(θ)) and the frequencies of the modes. In the output, the meridional modes are written in the following order: first the westward gravity modes, then the eastward gravity modes, and finally the (westward) Rossby modes. Both types of gravity modes are written following the symmetric–antisymmetric mode order, whereas the order for Rossby modes is reversed. If the VSE (13) is solved using boundary conditions (16) and (15), then the solutions for the barotropic mode returned by hough_functions are the Haurwitz waves (Swarztrauber and Kasahara, 1985). This is controlled internally by hough_functions which calls functions named hvf_baroclinic and hvf_barotropic for the Haurwitz modes.

• Expansion coefficients.

  The complex expansion coefficients w_nlk, ı_nlk, and j_nlk may be obtained with a function called expansion_coeffs. The inputs are the vertical structure functions (Ψ_k(p)), the equivalent heights (h_k), the Hough vector functions (Θ_nlk(θ)), and a data structure with appropriate fields. For the expansion coefficients w_nlk a data structure with fields $[u,v,\mathit{\varphi }]$ is required (recall that field ϕ corresponds to the deviations from a hydrostatic reference state ϕ₀(p)). The computation of expansion coefficients ı_nlk and j_nlk requires data structures with fields [I1,I2] and [J3], respectively. Each one of these fields need to be pre-computed by the user and correspond to the terms between square brackets of Eqs. (59), (60) and (61), respectively (note that in Eq. 61 the vertical derivative is outside the square brackets, and that this is accounted for in the expansion_coeffs function). The returned coefficients, w_nlk, ı_nlk and j_nlk, may then be used to compute the 3-D spectrum of total energy (E_nlk) and of the non-linear interactions of kinetic (I_nlk) and available potential (J_nlk) energies (see Fig. 9), using Eqs. (45), (47), and (48), respectively.

• The inverse transforms.

  The zonal and meridional wind, u and v, and geopotential perturbation ϕ (from the reference geopotential), as given by Eq. (34), can be obtained with the function inv_expansion_coeffs. By choosing an appropriate subset of modes, this function can be used for spatial-scale filtering of atmospheric motion, balanced and mass unbalanced motions, and barotropic and baroclinic components. Moreover, choosing the appropriate meridional indices the filtering can be used to isolate tropical components of the atmospheric circulation. The function requires the vertical structure functions (Ψ_k(p)); the equivalent heights (h_k); the Hough vector functions (Θ_nlk(θ)); the complex expansion coefficients (w_nlk); the pressure levels (in hPa); the longitudes (in degrees); and the wavenumber, meridional, and vertical indices as inputs. In addition, the user can choose to compute all three fields, u,v and ϕ (the default), or only a subset of those, which is controlled by the input argument uvz.

• Vertical and Hough transforms.

  All functions described above have the same name in both versions, MATLAB and Python, differing only in the file extension (“.m” and “.py”, respectively). The software tools presented here include two more functions used to compute the vertical and the Hough transforms. In the MATLAB version there are two separate functions, vertical_transform.m and hough_transform.m, whereas for Python both are included in module transforms.py, being invoked as transforms.vertical and transforms.hough. The vertical and Hough transforms are used by the expansion_coeffs function, but they may be also executed independently by the user. For example, the vertical_transform function (or transforms.vertical) may be used to compute the vertical transforms of u, v, ϕ, I1, I2, and J3. Again, I1, I2, and J3 correspond to the terms between square brackets of Eqs. (59), (60), and (61). The computation of these six vertical transforms constitute the essential task to obtain the global energy cycle in the wavenumber and vertical domains as in Marques and Castanheira (2012) or Marques and Castanheira (2017) (see also Castanheira and Marques, 2019). Having these vertical transforms, all of the terms of the normal mode energetics formulation, represented by balance Eqs. (50) and (51), may then be easily obtained using Eqs. (52)–(56), with the remainder dissipation and generation terms computed as residuals from Eqs. (50) and (51).