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

Abstract. A free software package for the computation of the 3-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 partition the energy reservoirs and energy interactions 5 between 3-D spatial scales, and between barotropic and baroclinic components as well as between 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, and 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 these features are useful both in data analysis and in 10 assessment of general circulation atmospheric models.


Introduction
The use of the three-dimensional normal mode functions (3-D NMFs) of the linearized 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 energy of the global motion in two kinds of 15 motions: mass balanced rotational modes (Rossby/Hawitz waves) and divergent modes (inertio-gravity waves).
The original work of Kasahara and Puri (1981) was developed for terrain following sigma-coordinates. Four years later, Tanaka (1985) extended the methodology to isobaric coordinates. 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 3-dimensional scales and types of motions. Recently, Marques and Castanheira (2012) 20 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 extratropi-the specific data set 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 data sets (such as reanalysis and outputs from climate models) and requiring the minimum user interaction as we found possible.

The basic equations
In the traditional shallow-atmosphere approximation, a set primitive equations in isobaric coordinates may be written as ∂u ∂t − 2Ωv sin θ + 1 a cos θ ∂φ ∂λ = −V · ∇u − ω ∂u ∂p + tan θ a uv + F u , ∂v ∂t + 2Ωu sinθ + 1 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. VariableQ 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 (5) by R/(pS 0 ), then calculating the derivatives with respect to p, and finally using on the resulting equation the hydrostatic and continuity equations ( (3) and (4)), the thermodynamic energy equation takes the form The right-hand side of equations (1), (2) and (7) contain the nonlinear terms, frictional forces and the diabatic heat sources.

25
The VSE (13) with appropriate boundary conditions form an eigenvalue problem, with discrete eigensolutions Ψ k (p) and eigenvalues −1/gh k (Cohn and Dee, 1989). The software here presented 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 = dz/dt 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 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 normalized 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 circulations variables 15 (u, v, φ) T , may be expanded as where In this software package, the VSE is solved by using the spectral method introduced by Kasahara (1984) (see also Castanheira et al. 20 (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 is 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 + 20. The vertical integration is performed by Gaussian quadrature interpolating T 0 into G L = 2J − 1 Gaussian levels. The maximum number of VSFs saved in the output is equal to the number of pressure levels (N k ). Fig. 1 shows the vertical profiles of temperature T 0 and base-10 logarithm of the stability parameter S 0 . 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 1.5 • lon. × 1.5 • lat. grid resolution for 5 the 37 isobaric levels from 1 to 1000 hPa, at time intervals of six-hours. The first 12 vertical structures associated with the 12 smallest eigenvalues ξ 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 figures 2 and 3, i.e. the vertical structures with one or several nodes are close to each other. 10 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 correspond exactly to the mass weighted 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.

15
Applying the vertical transform (18) to the system of equations (8)-(10), a dimensionless equation is obtained in the following vector form (Tanaka, 1985;Marques and Castanheira, 2012): wherê t = 2 Ω t,  Figure 1. Left: Vertical profile of temperature at the ERA-interim pressure levels (black dots) and interpolated into GL gaussian levels (red line). Right: stability parameter S0 calculated from T0 at original pressure levels (black dots) and from interpolated T0 at gaussian levels (red line). The stability parameter S0 was made dimensionless multiplying eq. (6) by p 2 s /gH * with H * = 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.  (14) and (15).
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 Since equation (21) is a linear system with respect to λ andt, the solution W k can be expressed as zonal waves (Swarztrauber and Kasaha 1985) where σ nlk are the dimensionless frequencies for the free waves with zonal wave number n and meridional structures Θ nlk (θ).
The functions H nlk (λ, θ) = Θ nlk (θ) e i n λ are known as the horizontal structure functions, and the meridional structures Θ nlk 5 are the Hough vectors functions (Longuet-Higgins, 1968;Swarztrauber and Kasahara, 1985). Substituting (27) into (21) gives and thus the horizontal structure functions are obtained as free eigenvalue problem. Equation (28) is referred 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 Kasah 1985;Žagar et al., 2015b).
It can be shown that for real α k all the frequencies σ of L are real and that any modes H corresponding to distinct frequencies 5 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 right-hand side is unity if n = n ′ and l = l ′ , 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, 15 1985). Substituting eq. (30) into H nlk (λ, θ) = Θ nlk (θ) e i n λ 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 westward-propagating 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).

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 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 5 n and meridional indices l. The computational method was developed by Swarztrauber and Kasahara (1985) and consists of expanding the eigensolutions of (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).  Swarztrauber and Kasahara (1985), because we have used the same values of α has they did to compute the frequencies σ.
The α used here corresponds to 1/ √ ǫ of Swarztrauber and Kasahara (1985). The frequencies are plotted as a function of α −1 for 10 meridional indices (l = 0, . . ., 9), five symmetric (continous lines) and five antisymmetric (dashed lines). The logarithm scale was used in both axis 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 15 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 gravityinertia 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) 20 are zero for n = 0 (Swarztrauber and Kasahara, 1985). Therefore, instead of saving zeros for all 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 ).
Left panel of Figure 5 shows the curves of asymptotic dimensionless frequencies σ a in the case of n → 0. The frequencies 25 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 r = −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 r = −1, which is similar to the eastward propagating Kelvin wave in Figure 4 (right panel). For this reason, the l r = −1 mode is saved as the eastward gravity mode l e = 0 as in Swarztrauber and Kasahara (1985). The curves for the modes l r = 1, 2, . . ., 9 behave analogously to those of rotational modes  12 https://doi.org/10.5194/gmd-2019-340 Preprint. Discussion started: 5 March 2020 c Author(s) 2020. CC BY 4.0 License.
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 subsection 2.2, using boundary conditions (14) and (15) (the first 12 VSFs are represented in Figure 2), and for 43 wavenumbers (n = 0, . . ., 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 5 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). Increasing the vertical index k (i.e. reducing the equivalent height h k ) or increasing the zonal wavenumber n for the same vertical index k, the horizontal structure functions become more confined around the equator (Longuet-Higgins, 1968; 10 Cohn and Dee, 1989;Žagar et al., 2015b).

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 15 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 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). Fig. 8 shows the same meridional profiles as in Fig. 6, but for the Haurwitz modes.

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 Π nlk (λ, θ, p), form a complete orthogonal basis, allowing therefore to expand the horizontal wind and the geopotential fields of the global atmosphere (Kasahara and Puri, 1981;Tanaka, 25 1985;Daley, 1991): where Π nlk (λ, θ, p) = Ψ k (p) Θ nlk (θ) e i n λ . 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:

3-D normal mode energetics scheme 5
Using the 3-D NMFs of the linearized primitive equations ( (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 onedimensional spectral energetics in domains of zonal wavenumber, n, meridional mode number, l, and vertical mode number, k.
The 3-D NME scheme complements therefore 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 10 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 summarized as follows: Applying the vertical transform (18) to equations (1), (2) and (7), a dimensionless equation is obtained in the following vector form: where subscript k denotes the k th component of the vertical transform andt, W k and L are given by (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 ,  nlk and s nlk are the Fourier-Hough transforms of the vector variables W k , I k , J k 5 and S k , respectively.
The kinetic energy, K, and available potential energy, A, are defined respectively by 10 And the total energy, E = K + A, is conserved a quantity of the full nonlinear equations (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 where ǫ n = 2 for the zonal (n = 0) modes, and ǫ n = 1 for the eddy (n > 0) modes.

15
Substituting (42) into the time derivatives of (45), the energy balance equations for the 3-D normal modes are finally obtained as d dt where The terms I nlk , J nlk , and S nlk that contribute to the time change of E nlk are respectively associated with nonlinear interac-25 tions of kinetic and available potential energies and with an energy source or sink due to diabatic heating and dissipation.The spetra of these terms can therefore be analysed separately for the zonal mean and eddy components, for the barotropic and baroclinic modes, and also for the Rossby and gravity waves (Tanaka and Kung, 1988).
As another illustration of the capabilities of the current software, in Figure 9 it is shown the spectra for total energy E (top panels), interactions of available potential energy J (middle panels) and interactions of kinetic energy I (bottom panels). The left column of panels contain the wavenumber spectra whereas the middle and right columns show the vertical spectra for the zonal mean and eddy components, respectively. In each panel it is represented the total spectra along with the contributions of the Rossby and Gravity components. These spectra were computed using 34 years of ERA-interim reanalysis (Dee et al.,5 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 1.5 • lon. × 1.5 • lat. grid resolution for the 37 isobaric levels from 1000-to 1 hPa, at time intervals of six-hours. 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 years period.

10
A detailed analysis of the spectra in figure 9 is beyond the scope of this study, but some of its characteristics are mentioned.
The wavenumber energy spectrum of Rossby modes follows approximately the -3 power law over the synoptic and mesoscale region, while that of the gravity modes has approximately a -5/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 contibution by 15 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 20 energy due to the gravity waves is in the opposite direction.
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, also the conversion rate of available potential energy into kinetic energy cannot be acessed. 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) 25 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 also 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 30 energy (A) is given by (see Marques and Castanheira, 2012, for details) where 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 5 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 equations (50) and (51). The calculated energetics terms may then be used to construct an extended energy cycle diagram, as in Fig (13)) were obtained using boundary conditions (14) and (15). The energetics terms were computed at each time step (6-hourly) and then averaged over the 34 years 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 figure 10.

15
The energy cycle diagram of figure 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 20 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)

25
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 later, 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 30 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 can not be derived from de eigenvalues of the VSE (13). In this case the smallest eigenvalue is null, and corresponds to the separation constant in eq. (12) to be null. For such case, it is not possible to define an equivalent  height as in eq. (12). In order to maintain 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. Performing in this way the equations (50) and (51), with its individual terms given by equations (52)-(58), may be also used to construct an extended energy cycle diagram as illustrated in figure 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 figure 10. The estimates in 5 this diagram were based on the same data as those of figure 10. The solutions of the VSE were also approximated using 57 Legendre polynomials, being retained 37 vertical modes. It was also considered the same zonal wavenumbers as those for the estimates in figure 10.
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 be shown easily that the total energy is exactly conserved in by full non-linear equations (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 5 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 guaranty the conservation of energy for frictionless adiabatic motions it is necessary to impose a non-slipping condition (u, v) s = 0 at p = p s . And this may seem contradictory with the absence of friction.

10
The software here presented 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) 15 and the computation of complex expansion coefficients w nlk , ı nlk and  nlk , which are the vertical-Fourier-Hough transforms of the dependent variable vector [u, v, φ] 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 equations (50) and (51) and its ex-20 pressions are given in equations (52)-(64). With these vertical transforms the user may compute the global energy cycle in the wavenumber and vertical domains (as in Castanheira, 2012, 2017;Castanheira and Marques, 2019), with the dissipation and generation terms computed as residuals from (50) and (51).
Other applications of this software include the spatial scale filtering 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,25 2015; Žagar et al., 2019). Furthermore, the solutions of HSE, i.e. the Hough vectors 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 these features are useful both in data analysis and in assessment of general circulation atmospheric models (e.g. Žagar et al., 2019).
(h k ). The VSE may be solved using either one of the two pairs of boundary conditions as described in section 2.2. These pairs of boundary conditions are controled 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 later case the first vertical structure is strictly constant and the associated equivalent height is infinite.
-Horizontal structure 10 The solution of HSE (21) is obtained with function called hough_functions. It uses the computational method developed by Swarztrauber and Kasahara (1985) which consists of expanding the eigensolutions of (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 vectors functions (Θ nlk (θ)) 15 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 controled internaly by 20 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  nlk may be obtained with function called expansion_coeffs. The inputs are the vertical structure functions (Ψ k (p)), the equivalent heights (h k ), the Hough vectors functions (Θ nlk (θ)) and a data structure with appropriate fields. For the expansion coefficients w nlk it is required a data structure with fields 25 [u, v, φ] (recall that field φ corresponds to the deviations from a hydrostatic reference state φ 0 (p)). The computation of expansion coefficients ı nlk and  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 (59), (60) and (61), respectively (Note that in (61) the vertical derivative is outside the square brackets, and this is accounted in the expansion_coeffs function). The returned coefficients, w nlk , ı nlk and  nlk , may then be used to compute the 3-D 30 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 equations (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 equation (34), can be obtained with function inv_expansion_coeffs. By choosing an appropriate subset of modes, this function can be used for spatial scale filtering of atmospheric motion, and filtering of balanced and mass unbalanced motions, as well as of 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 as inputs the 5 vertical structure functions (Ψ k (p)), the equivalent heights (h k ), the Hough vectors functions (Θ nlk (θ)), the complex expansion coefficients (w nlk ), the pressure levels (in hPa), the longitudes (in degrees) and the wavenumber, meridional and vertical indices. In addition, the user can choose to compute all three fields, u, v and φ (the default), or only a subset of those, which is controled by the input argument uvz.
-Vertical and Hough transforms  (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 the 20 terms of the normal mode energetics formulation represented by balance equations (50) and (51), may then be easily obtained using equations (52)-(56), with the remainder dissipation and generation terms computed as residuals from (50) and (51).
Code and data availability. The exact version of the software code used to produce the results presented in this paper is archived on Zenodo (Marques et al., 2020). A step by step tutorial is included in the repository for both Python and MATLAB. Also included are the data sets 25 with the mean vertical profiles of temperature and geopotential, computed from 32 years  of the ERA-Interim reanalysis data as described in the text (see section 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 ERA-Interim data server with 1.5 • lon. × 1.5 • lat. grid resolution for all the provided 37 isobaric levels from 1000-to 1 hPa, at time intervals of six-hours, are the only data sets needed to obtain all the figures in the paper. The two exceptions are figures 4 and 5 for which we have used the same values as Swarztrauber and Kasahara (1985) (see for example their Table