the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Threedimensional normal mode functions: openaccess tools for their computation in isobaric coordinates (p3DNMF.v1)
Carlos A. F. Marques
Martinho MartaAlmeida
A free software package for the computation of the threedimensional normal modes of an hydrostatic atmosphere is presented. This software performs the computations in isobaric coordinates and was developed for two userfriendly languages: MATLAB and Python. The software can be used to expand the global atmospheric circulation onto the 3D normal modes. This expansion allows the computation of a 3D energetic scheme, which partitions the energy reservoirs and energy interactions between 3D 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 3D normal mode expansion can be used for spatialscale filtering of atmospheric motion, filtering of balanced motion and mass unbalanced motions, and barotropic and baroclinic components. Fixing the meridional scale, the 3D 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.
The use of the threedimensional normal mode functions (3D 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 3D 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 (inertiogravity 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 3D normal mode energy scheme, including the exchanges of kinetic and available potential energy (APE) among the different threedimensional scales and types of motions. Recently, Marques and Castanheira (2012) extended the methodology to analyse the conversion of APE into kinetic energy.
The 3D 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 lowfrequency 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 opensource software package for the projection of global horizontal wind and mass fields onto 3D 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 3D NMFs using isobaric coordinates. This software allows for the analysis of the full atmospheric energy cycle and was developed for two userfriendly languages: MATLAB and Python. The methodology described in this study to obtain the 3D 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.1 The basic equations
In the traditional shallow atmosphere approximation, a set of primitive equations in isobaric coordinates may be written as follows:
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_{0}(p) and ϕ_{0}(p), respectively, and the static stability parameter, S_{0}, is given by
Multiplying Eq. (5) by R∕(pS_{0}), 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
The righthand side of Eqs. (1), (2), and (7) contain the nonlinear terms, frictional forces, and the diabatic heat sources. By setting the righthand sides of those equations to zero, the linearised system of equations for the three dependent variables u, v and ϕ is written as
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),
and substituting into Eq. (10), one obtains
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)
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
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
The vertical structure functions (VSFs), Ψ_{k}(p), can be normalised to satisfy the following orthonormality condition
Using this orthonormality condition, one can define a vertical transform
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
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 finitedifference method that the derivatives of the vertical structure function can be calculated by analytical differentiation. The input data are the temperature profile T_{0}, 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_{0} 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_{0} and the base 10 logarithm of the stability parameter S_{0}. The profiles were computed based on the 3D temperature field of the ERAInterim reanalysis (Dee et al., 2011), covering the period from 1 January 1979 to 31 December 2010. The 3D temperature field was obtained with a $\mathrm{1.5}{}^{\circ}\phantom{\rule{0.25em}{0ex}}\text{lat.}\times \mathrm{1.5}{}^{\circ}\phantom{\rule{0.25em}{0ex}}\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_{0} and stability parameter S_{0} represented in Fig. 1.
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, Ψ_{0}, 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 massweighted vertical mean of that variable. The vertical structure function Ψ_{0} in Fig. 2 is not strictly constant. Nevertheless, since Ψ_{0} 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):
where
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
The dimensionless parameter, α_{k}, in the linear differential matrix operator L is defined as
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)
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)\phantom{\rule{0.125em}{0ex}}{e}^{i\phantom{\rule{0.125em}{0ex}}n\phantom{\rule{0.125em}{0ex}}\mathit{\lambda}}$ are known as the horizontal structure functions, and the meridional structures Θ_{nlk} are the Hough vector functions (LonguetHiggins, 1968; Swarztrauber and Kasahara, 1985). Substituting Eq. (27) into Eq. (21) gives
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 highfrequency eastward and westwardpropagating gravity inertia waves, and the second kind describes lowfrequency westwardpropagating rotational waves of the Rossby–Haurwitz type (LonguetHiggins, 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
where the superscript ∗ denotes a conjugate transpose and the righthand 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
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)\phantom{\rule{0.125em}{0ex}}{e}^{i\phantom{\rule{0.125em}{0ex}}n\phantom{\rule{0.125em}{0ex}}\mathit{\lambda}}$ in Eq. (29), one can see that the Hough vector functions satisfy the orthonormal condition given by
Using orthonormal condition (29), a complete set of Fourier–Hough transforms may be constructed as
In order to distinguish each wave type, i.e. the westwardpropagating Rossby wave and the westward and eastwardpropagating 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 westwardpropagating gravity inertia and rotational waves (Fig. 4a) and eastwardpropagating 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{\u03f5}}$ of Swarztrauber and Kasahara (1985). The frequencies are plotted as a function of α^{−1} 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 westwardpropagating 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.
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 α^{−1} 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 eastwardpropagating 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.
As mentioned above, the horizontal modes include Kelvin, westward and eastward inertiogravity (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_{0}≃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_{1}≃6.1 km and k=10, with h_{10}≃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 (LonguetHiggins, 1968; Cohn and Dee, 1989; Žagar et al., 2015b).
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 nondivergent, 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 nondivergent, 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.
2.5 3D normal mode functions
The threedimensional (3D) 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 3D 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):
where
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:
with $\mathbf{W}={\left[u,\phantom{\rule{0.125em}{0ex}}v,\phantom{\rule{0.125em}{0ex}}\mathit{\varphi}\right]}^{T}$.
2.6 3D normal mode energetics scheme
Using the 3D 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 3D normal mode energetics (NME) scheme, which combines three onedimensional spectral energetics in domains of the zonal wavenumber, n; meridional mode number, l; and vertical mode number, k. The 3D NME scheme therefore complements the standard energetics in the zonal wavenumber domain of Saltzman (1957), since it can diagnose the 3D 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 3D 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:
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 nonlinear 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
These vectors are made dimensionless using the inverse of scaling matrix Y_{k}, which is given by
Applying the Fourier–Hough transform (33) to (37), yields
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
and the total energy, $E=K+A$, is a conserved quantity of the full nonlinear Eqs. (1), (2), and (7) under frictionless and adiabatic conditions (see Tanaka and Kung, 1988).
In the transformed 3D normal mode space, the total energy E_{nlk} associated with each mode (nlk) is
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 3D normal modes are finally obtained as
where
The terms I_{nlk},J_{nlk}, and S_{nlk} that contribute to the time change of E_{nlk} are, respectively, associated with nonlinear 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 ERAInterim 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}\phantom{\rule{0.25em}{0ex}}\text{lat.}\times \mathrm{1.5}{}^{\circ}\phantom{\rule{0.25em}{0ex}}\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 (6hourly) and then averaged over the 30year 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 planetaryscale 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.
A disadvantage of the 3D 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 nonlinear 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)
where
and
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 ERAInterim reanalysis mentioned above for the 3D 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 (6hourly) and then averaged over the 34year 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_{0} 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_{0}=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 nonlinear 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 massweighted 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 nonslipping condition $(u,v{)}_{\mathrm{s}}=\mathrm{0}$ at p=p_{s}, though this may seem contradictory with the absence of friction.
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 nonlinear term vector due to wind field, and of the nonlinear 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 spatialscale 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).

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_{0}, 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 userdefined 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 ϕ_{0}(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 precomputed 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 3D spectrum of total energy (E_{nlk}) and of the nonlinear 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 spatialscale 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).
The exact version of the software code used to produce the results presented in this paper is archived on Zenodo (https://doi.org/10.5281/zenodo.3631310, Marques et al., 2020). A stepbystep tutorial is included in the repository for both Python and MATLAB. Also included are the mean vertical profiles of temperature and geopotential, computed from 32 years (1979–2010) of the ERAInterim reanalysis data as described in the text (see Sect. 2.2). These two vertical profiles, along with the horizontal wind (u, v), pressure velocity (ω), temperature, and geopotential fields that may be obtained from the ERAInterim data server with $\mathrm{1.5}{}^{\circ}\phantom{\rule{0.25em}{0ex}}\text{lat.}\times \mathrm{1.5}{}^{\circ}\phantom{\rule{0.25em}{0ex}}\text{long.}$ grid resolution for all of the provided 37 isobaric levels from 1000 to 1 hPa, at time intervals of 6 h, are the only datasets needed to obtain all of the figures in the paper. The two exceptions are Figs. 4 and 5, for which we have used the same values as Swarztrauber and Kasahara (1985); see, for example, their Table 5. Note that our α corresponds to their $\mathrm{1}/\sqrt{\mathit{\u03f5}}$.
CAFM conceived the idea, wrote the initial version of the software code for MATLAB, and tested the final version of the software for both MATLAB and Python. MMA wrote the software code for Python and the final version of the code for MATLAB and collaborated in the tests of the software. JMC elaborated the structure of the software and collaborated in writing the initial version of the software code for MATLAB and in the software tests. All authors contributed to writing the manuscript.
The authors declare that they have no conflict of interest.
The code is made publicly available without any warranty.
Carlos A. F. Marques was supported by the Portuguese Foundation for Science and Technology (FCT) under grant no. SFRH/BPD/76232/2011. The ERAInterim data have been obtained from the ECMWF data server. Thanks are due to FCT/MCTES for the financial support to CESAM (UIDP/50017/2020+UIDB/50017/2020) through national funds.
This research has been supported by the Portuguese Foundation for Science and Technology (FCT) (grant no. SFRH/BPD/76232/2011) and by FCT/MCTES through national funds for the financial support to CESAM (UIDP/50017/2020+UIDB/50017/2020).
This paper was edited by David Ham and reviewed by two anonymous referees.
Blaauw, M. and Žagar, N.: Multivariate analysis of Kelvin wave seasonal variability in ECMWF L91 analyses, Atmos. Chem. Phys., 18, 8313–8330, https://doi.org/10.5194/acp1883132018, 2018. a
Castanheira, J. M. and Marques, C. A. F.: Convectively Coupled Equatorial Waves Diagnosis using 3–D Normal Modes, Q. J. Roy. Meteor. Soc., 141, 2776–2792, https://doi.org/10.1002/qj.2563, 2015. a, b, c
Castanheira, J. M. and Marques, C. A. F.: The energy cascade associated with daily variability of the North Atlantic Oscillation, Q. J. Roy. Meteor. Soc., 145, 197–210, https://doi.org/10.1002/qj.3422, 2019. a, b, c, d, e, f
Castanheira, J. M., DaCamara, C. C., and Rocha, A.: Numerical solutions of the vertical structure equation and associated energetics, Tellus, 51A, 337–348, 1999. a, b
Castanheira, J. M., Graf, H.F., DaCamara, C. C., and Rocha, A.: Using a Physical Reference Frame to Study Global Circulation Variability, J. Atmos. Sci., 59, 1490–1501, https://doi.org/10.1175/15200469(2002)059<1490:UAPRFT>2.0.CO;2, 2002. a
Charney, J. G.: Geostrophic turbulence, J. Atmos. Sci., 28, 1087–1095, 1971. a
Cohn, S. E. and Dee, D. P.: An analysis of the vertical structure equation for arbitrary thermal profiles, Q. J. Roy. Meteor. Soc., 115, 143–171, 1989. a, b, c
Daley, R.: Atmospheric data analysis, Cambridge University Press, Cambridge, UK, 1991. a, b, c
Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., MongeSanz, B. M., Morcrette, J.J., Park, B.K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.N., and Vitart, F.: The ERA–Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a, b
Franzke, C. L. E., Jelic, D., Lee, S., and Feldstein, S. B.: Systematic decomposition of the MJO and its Northern Hemispheric extratropical response into Rossby and inertiogravity components, Q. J. Roy. Meteor. Soc., 145, 1147–1164, https://doi.org/10.1002/qj.3484, 2019. a
Gehne, M. and Kleeman, R.: Spectral Analysis of Tropical Atmospheric Dynamical Variables Using a Linear ShallowWater Modal Decomposition, J. Atmos. Sci., 69, 2300–2316, 2012. a
Kasahara, A.: Normal modes of ultralong waves in the atmosphere, Mon. Weather Rev., 104, 669–690, 1976. a
Kasahara, A.: Numerical Integration of the Global Barotropic Primitive Equations with Hough Harmonic Expansions, J. Atmos. Sci., 34, 687–701, 1977. a
Kasahara, A.: The Linear Response of a Stratified Global Atmosphere to Tropical Thermal Forcing, J. Atmos. Sci., 41, 2217–2237, 1984. a, b
Kasahara, A. and Puri, K.: Spectral representation of threedimensional global data by expansion in normal mode functions, Mon. Weather Rev., 109, 37–51, 1981. a, b, c, d, e, f, g, h
Kitsios, V., O'Kane, T. J., and Žagar, N.: A ReducedOrder Representation of the Madde–Julian Oscillation Based on Reanalyzed Normal Mode Coherences, J. Atmos. Sci., 76, 2463–2480, https://doi.org/10.1175/JASD180197.1, 2019. a
LonguetHiggins, M. S.: The eigenfunctions of Laplace's tidal equations over a sphere, Philos. T. Roy. Soc. London, A262, 511–607, 1968. a, b, c
Marques, C. A. F. and Castanheira, J. M.: A detailed normal mode energetics of the general circulation of the atmosphere, J. Atmos. Sci., 69, 2718–2732, 2012. a, b, c, d, e, f, g, h, i, j, k, l, m
Marques, C. A. F. and Castanheira, J. M.: Barotropic decelerations of the southern stratospheric polar vortex, Q. J. Roy. Meteor. Soc., 143, 744–755, https://doi.org/10.1002/qj.2961, 2017. a, b, c, d, e, f
Marques, C. A. F. and Castanheira, J. M.: Diagnosis of Free and Convectively Coupled Equatorial Waves, Math. Geosci., 50, 585–606, https://doi.org/10.1007/s110040189729y, 2018. a, b
Marques, C. A. F., Castanheira, J. M., and Rocha, A.: Changes in the normal mode energetics of the general atmospheric circulation in a warmer climate, Clim. Dynam., 42, 1887–1903, https://doi.org/10.1007/s0038201317508, 2014. a, b
Marques, C. A. F., MartaAlmeida, M., and Castanheira, J. M.: NMF3D – tridimensional normal mode functions, Zenodo, https://doi.org/10.5281/zenodo.3631310, 2020. a
Nastrom, G. D. and Gage, K. S.: A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft, J. Atmos. Sci., 42, 950–960, 1985. a
Saltzman, B.: Equations governing the energetics of the larger scales of atmospheric turbulence in the domain of wave number, J. Meteorol., 14, 513–523, 1957. a
Shigehisa, Y.: Normal modes of the shallow water equations for zonal wavenumber zero, J. Meteorol. Soc. Jpn., 61, 479–494, 1983. a, b
Swarztrauber, P. N. and Kasahara, A.: The vector harmonic analysis of Laplace's tidal equations, SIAM J. Sci. Stat. Comput., 6, 464–491, 1985. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u
Tanaka, H. L.: Global energetics analysis by expansion into threedimensional normalmode functions during the FGGE winter, J. Meteorol. Soc. Jpn., 63, 180–200, 1985. a, b, c, d, e, f, g
Tanaka, H. L. and Kung, E. C.: Normalmode energetics of the general circulation during the FGGE year, J. Atmos. Sci., 45, 3723–3736, 1988. a, b, c, d, e, f, g, h
Tanaka, H. L. and Seki, S.: Development of a ThreeDimensional Spectral Linear Baroclinic Model and its Application to the Baroclinic Instability Associated with Positive and Negative Arctic Oscillation Indices, J. Meteorol. Soc. Jpn. Ser. II, 91, 193–213, https://doi.org/10.2151/jmsj.2013207, 2013. a
Tanaka, H. L. and Tokinaga, H.: Baroclinic Instability in High Latitudes Induced by Polar Vortex: A Connection to the Arctic Oscillation, J. Atmos. Sci., 59, 69–82, https://doi.org/10.1175/15200469(2002)059<0069:BIIHLI>2.0.CO;2, 2002. a
Tanaka, H. L., Kung, E. C., and Baker, W. E.: Energetics analysis of the observed and simulated general circulation using threedimensional normal mode expansions, Tellus, 38A, 412–428, 1986. a
Terasaki, K. and Tanaka, H. L.: An Analysis of the 3D Atmospheric Energy Spectra and Interactions Using Analytical Vertical Structure Functions and Two Reanalyses, J. Meteorol. Soc. Jpn., 85, 785–796, 2007. a
Yamagami, A. and Tanaka, H. L.: Analysis of Unstable Solution with an MJO Structure Using a ThreeDimensional Spectral Linear Baroclinic Model, SOLA, 10, 103–107, https://doi.org/10.2151/sola.2014021, 2014. a
Yamagami, A. and Tanaka, H. L.: Characteristics of the JRA55 and ERAInterim Datasets by Using the ThreeDimensional Normal Mode Energetics, SOLA, 12, 27–31, https://doi.org/10.2151/sola.2016006, 2016. a
Yang, G.Y., Hoskins, B., and Slingo, J.: Convectively Coupled Equatorial Waves: A New Methodology for Identifying Wave Structures in Observational Data, J. Atmos. Sci., 60, 1637–1654, 2003. a
Žagar, N., Buizza, R., and Tribbia, J.: A ThreeDimensional Multivariate Modal Analysis of Atmospheric Predictability with Application to the ECMWF Ensemble, J. Atmos. Sci., 72, 4423–4444, https://doi.org/10.1175/JASD150061.1, 2015a. a
Žagar, N., Kasahara, A., Terasaki, K., Tribbia, J., and Tanaka, H.: Normalmode function representation of global 3D data sets: openaccess software for the atmospheric research community, Geosci. Model Dev., 8, 1169–1195, https://doi.org/10.5194/gmd811692015, 2015b. a, b, c, d
Žagar, N., Kosovelj, K., Manzini, E., Martin, H., and Castanheira, J.: An assessment of scaledependent variability and bias in global prediction models, Clim. Dynam., 54, 287–306, https://doi.org/10.1007/s0038201905001x, 2020. a, b