OZO v . 1 . 0 : Software for solving a generalized omega equation and the Zwack-Okossi height tendency equation using WRF model output

A software package (OZO, Omega-Zwack-Okossi) was developed to diagnose the processes that affect vertical motions and geopotential height tendencies in weather systems simulated by the Weather Research and Forecasting (WRF) model. First, this software solves a generalized omega equation to calculate the vertical motions associated with different physical forcings: vorticity advection, thermal advection, friction, diabatic heating, and an imbalance term between vorticity and temperature tendencies. After this, the corresponding height tendencies are calculated with the Zwack-Okossi tendency 5 equation. The resulting height tendency components thus contain both the direct effect from the forcing itself and the indirect effects (related to the vertical motion induced by the same forcing) of each physical mechanism. This approach has an advantage compared with previous studies with the Zwack-Okossi equation, in which vertical motions were used as an independent forcing but were typically found to compensate the effects of other forcings. The software is currently tailored to use input from WRF simulations with Cartesian geometry. As an illustration, results for 10 an idealized 10-day baroclinic wave simulation are presented. An excellent agreement is found between OZO and the direct WRF output for both the vertical motion (correlation 0.95 in the midtroposphere) and the height tendency fields(correlation 0.95-0.97 in the whole troposphere). The individual vertical motion and height tendency components demonstrate the importance of both adiabatic and diabatic processes for the simulated cyclone. OZO is an open source tool for both research and education, and the distribution of the software will be supported by the authors. 15


Introduction
Today, high-resolution atmospheric reanalyses provide a three-dimensional (3-D) view on the evolution of synopticscale weather systems (Dee et al., 2011;Rienecker et al., 2011). On the other hand, simulations by atmospheric models allow for exploring the sensitivity of both real-world and idealised weather systems to factors such as the initial state (e.g. Leutbecher et al., 2002;Hoskins and Coutinho, 2005), lower boundary conditions (e.g. Elguindi et al., 2005;Hirata et al., 2016), and representation of sub-grid scale processes (e.g. Wernli et al., 2002;Liu et al., 2004;Beare, 2007). Nevertheless, the complexity of atmospheric dynamics often makes the physical interpretation of reanalysis data and model output far from simple. Therefore, there is also a need for diagnostic methods that help to separate the effects of individual dynamical and physical processes on the structure and evolution of weather systems.
Two variables that are of special interest in the study of synoptic-scale weather systems are the geopotential height tendency and vertical motion (Holton and Hakim, 2012). Height tendencies are directly related to the movement and intensification or decay of low-and high-pressure systems. Vertical motions affect atmospheric humidity, cloudiness, and precipitation. They also play a crucial role in atmospheric dynamics by inducing adiabatic temperature changes, by generating cyclonic or anticyclonic vorticity, and by converting available potential energy to kinetic energy (Lorenz, 1955;Holton and Hakim, 2012).
For the need of diagnostic tools, some software packages have been developed to separate the contributions of each forcing to the vertical motion and height tendency. DIONYSOS (Caron et al., 2006), a tool for analysing numerically simulated weather systems, provides currently online daily diagnostics for the output of numerical weather prediction models. RIP4 (Stoelinga, 2009) can calculate Q-vectors (Hoskins et al., 1978;Holton and Hakim, 2012) and a quasigeostrophic (QG) vertical motion (Holton and Hakim, 2012) from Weather Forecast and Research (WRF) model output, but the division of ω into contributions from various atmospheric processes is not possible in RIP4. Furthermore, many research groups have developed tools for their own needs but do not have resources to distribute the software.
Here, we introduce a software package Omega-Zwack-Okossi (OZO) that can be used for diagnosing the contributions of different dynamical and physical processes to atmospheric vertical motions and height tendencies. OZO calculates vertical motion from a quasi-geostrophic and a generalised omega equation (Räisänen, 1995) while height tendencies are calculated using the Zwack-Okossi tendency equation (Zwack and Okossi, 1986;Lupo and Smith, 1998). The current, first version of OZO has been tailored to use output from the WRF model (Wang et al., 2007;Shamarock et al., 2008) simulations run with idealised Cartesian geometry. Due to the wide use of the WRF model, we expect OZO to be a useful open-source tool for both research and education.
In the following, we first introduce the equations solved by OZO: the two forms of the omega equation in Sect. 2.1 and the Zwack-Okossi height tendency equation in Sect. 2.2. The numerical techniques used in solving these equations are described in Sect. 3. We have tested the software using output from an idealised 25 km resolution WRF simulation described in Sect. 4. Section 5 provides some computational aspects of the software. The next two sections give an overview of the vertical motion (Sect. 6) and height tendency (Sect. 7) calculations for this simulation. Software limitations and plans for future development are presented in Sect. 8, and the conclusions are given in Sect. 9. Finally, information about the data and code availability is given in Sect. 10.

Omega equation
The omega equation is a diagnostic tool for estimating atmospheric vertical motions and studying their physical and dynamical causes. Its well-known QG form, obtained by com-bining the QG vorticity and thermodynamic equations, infers vertical motion from geostrophic advection of absolute vorticity and temperature (Holton, 1992): where and the two right-hand side (RHS) terms are (4) (the notation is conventional, see Table 1 for an explanation of the symbols.) Qualitatively, the QG omega equation indicates that cyclonic (anticyclonic) vorticity advection increasing with height and a maximum of warm (cold) advection should induce rising (sinking) motion in the atmosphere. However, when deriving this equation, ageostrophic winds, diabatic heating, and friction are neglected. In addition, hydrostatic stability is treated as a constant and several terms in the vorticity equation are omitted. Although Eq. (1) often provides a reasonable estimate of synoptic-scale vertical motions at extratropical latitudes (Räisänen, 1995), these approximations inevitably deteriorate the accuracy of the QG omega equation solution.
The omega equation can be generalised by relaxing the QG approximations (e.g. Krishnamurti, 1968;Pauley and Nieman, 1992;Räisänen, 1995). Here we use the formulation where is a generalised form of Eq.
(2) and the RHS terms have the expressions Apart from the reorganisation of the terms in Eq. (6), this generalised omega equation is identical with the one used by Räisänen (1995). It follows directly from the isobaric primitive equations, which assume hydrostatic balance but omit the other approximations in the QG theory. The first four terms on the RHS represent the effects of vorticity advection (F V ), thermal advection (F T ), friction (F F ), and diabating heating (F Q ). The last term (F A ) describes imbalance between the temperature and vorticity tendencies. For constant f and constant R, the imbalance term is directly proportional to the pressure derivative of the ageostrophic vorticity tendency (Räisänen, 1995).
Because the operators L QG and L are linear, the contributions of the various RHS terms to ω can be calculated separately if homogeneous boundary conditions (ω = 0 at horizontal and vertical boundaries) are used. These contributions will be referred to as ω X , where X identifies the corresponding forcing term. The contribution of orographic vertical motions could be added by using a non-zero lower boundary condition (Krishnamurti, 1968) but is not included in the current version of OZO.
The geostrophic vorticity tendency at level p L is obtained from the equation where p s (here 1000 hPa) and p t (here 100 hPa) are the lower and the upper boundaries of the vertical domain. The vorticity tendency ∂ζ ∂t is calculated from the vorticity equation (Eq. 2 in Räisänen, 1997) and the temperature tendency ∂T ∂t from the thermodynamic equation (Eq. 3.6 in Holton and Hakim, 2012). The ageostrophic vorticity tendency ∂ζ ag ∂t is estimated from time series of vorticity and geostrophic vorticity (Eq. 12) using centred time differences ∂ζ ag ∂t ≈ For the calculations shown in this paper, t = 1800 s was used.
In Eq. (14), the first integral gives the mass-weighted vertical average of the geostrophic vorticity tendency between the levels p s and p t . The difference between the geostrophic vorticity tendency at level p L and this mass-weighted average is obtained from the second double integral. In this latter integral, hydrostatic balance is assumed to link the pressure derivative of the geostrophic vorticity tendency to the Laplacian of the temperature tendency.
Analogously with the vertical motion, the height tendency is divided in OZO to the contributions of different physical and dynamical processes as ii. Thermal advection (T ) and diabatic heating (Q) have a direct effect on the temperature tendency in Eq. (14).
iii. The ageostrophic vorticity tendency in Eq. (14) is attributed to the imbalance term (A).
iv. All five terms also affect the vorticity and temperature tendencies indirectly via vertical motions, which are calculated for each of them separately with the generalised omega equation.
This results in the following new expressions: The equation system used in this study has been adopted partly from Räisänen (1997). However, whereas Räisänen (1997) used the vorticity equation and the non-linear balance equation to obtain height tendencies, the Zwack-Okossi equation is used here. The main advantage of this choice is its smaller sensitivity to numerical errors. Our method produces quite smooth vertical profiles of height tendencies because the tendencies at neighbouring levels are bound to each other by the vertical integration in Eq. (14). On the other hand, our method differs from earlier applications of the Zwack-Okossi equation (e.g. Zwack and Okossi, 1986;Lupo and Smith, 1998) because the use of the generalised omega equation eliminates vertical motion as an independent height tendency forcing. This is an important advantage because these earlier studies have shown a tendency of compensation between vertical motions and the other forcing terms.
Earlier diagnostic tools have come close to our technique. The most similar approach is probably used in the DIONYSOS (Caron et al., 2006) tool. Regardless of many similarities, there are still three major differences between DIONYSOS and OZO. First, DIONYSOS eliminates the ageostrophic vorticity tendency as an independent forcing using an iterative procedure. Second, DIONYSOS uses simple parameterisations of diabatic heating and friction, whereas OZO relies directly on model output. Third, DIONYSOS uses the method of Räisänen (1997) to convert vorticity tendencies to height tendencies, whereas the Zwack-Okossi method is used in OZO.

Vorticity and temperature advection by non-divergent and divergent winds
Following the Helmholtz theorem, the horizontal wind can be divided to non-divergent (V ψ ) and divergent (V χ ) parts. Their contributions to vorticity advection and temperature advection can be separated as and the same applies to the corresponding ω and height tendency contributions. This separation between V ψ and V χ contributions is included in OZO because Räisänen (1997) found it to be important for the height tendencies associated with vorticity advection. OZO calculates the divergent part of the wind (V χ ) from the gradient of the velocity potential χ (Eq. 24), which is derived from the horizontal divergence by inverting the Laplacian in Eq. (25).
The non-divergent wind is obtained as the difference between V and V χ . The output of OZO explicitly includes the vorticity advection and temperature advection terms of the ω and height tendency equations due to both the full wind field V and the divergent wind V χ . The contributions associated with the non-divergent wind V ψ can be calculated as their residual.

Numerical methods
The first version of the OZO software package is tailored to use output from WRF simulations in idealised Cartesian geometry. The computational domain is periodic in the zonal direction, whereas symmetric boundary conditions are used at the northern and southern boundaries. Before the calculation, the WRF data need to be interpolated to pressure coordinates.

Calculation of right-hand side terms
All of the right-hand side terms of the omega equation (Eq. 5) and the Zwack-Okossi equation (Eq. 16) are calculated in grid point space. Horizontal and vertical derivatives are approximated with two-point central differences with the exception at the meridional and vertical boundaries, where onesided differences are used. In the calculation of the imbalance term of the omega and Zwack-Okossi equations (Eqs. 11 and 21), also tendencies of T , ζ and Z are needed. These tendencies are approximated by central differences of the corresponding variables.
Because the calculations are done in pressure coordinates, the lower boundary of the domain does not correspond to the actual surface. To mitigate the impact of this, vorticity and temperature advection, friction, diabatic heating, and the ageostrophic vorticity tendency are all attenuated below the actual surface by multiplying them by a factor varying from 0 to 1. The multiplication factor at each level depends on how far down the level is below the ground. For example, for a surface pressure of 950 hPa and vertical mass-centred grid spacing of 50 hPa, the multiplication factor is 0 at 1000 hPa, 0.5 at 950 hPa, and 1 at 900 hPa and all higher levels.

Inversion of left-hand side operators
The omega equation is solved using a multigrid algorithm (Fulton et al., 1986). Each multigrid cycle starts from the original (finest) grid, denoted below with the superscript (1). A rough solution in this grid ( ω (1) ) is found using ν 1 iterations of simultaneous under relaxation, starting either from the previous estimate of ω or (in the first cycle) ω = 0. The residuals, are then upscaled to a coarser grid (superscript (2)), in which the number of points is halved in all three dimensions. In this grid, the equation is then roughly solved replicating the method used in the first grid. The residual of this calculation is fed to the next coarser grid, and the procedure is continued until the grid only has five points on its longest (meridional) axis. Thus, for an idealised baroclinic wave simulation with horizontal resolution of 25 km, the meridional axis has 320 grid points. That makes six coarser resolutions (160, 80, 40, 20, 10, and 5 points on the meridional axis) in addition to the original one.
Having obtained the estimate ω (N) for the coarsest grid, a new estimate for the second coarsest grid is formed as where α is a relaxation coefficient. To reduce the effect of regridding errors, this estimate is refined using ν 2 iterations of simultaneous under relaxation. The result, ω NEW2 , is then substituted to the next finer grid and the sequence is repeated until the original grid is reached. After each multigrid cycle, the maximum difference between the new and the previous estimate of ω in the original (finest) grid is computed. If this difference exceeds a userdefined threshold (by default 5 × 10 −5 Pa s −1 ), the multigrid cycle is repeated. Typically, several hundreds of cycles are required to achieve the desired convergence.
OZO has four parameters for governing the multigrid algorithm, with the following default values: the under relaxation coefficient (α = 0.017), the number of sub-cycle iterations in the descending (ν 1 = 30) and ascending phases of the multigrid cycle (ν 2 = 4), and the threshold for testing convergence (toler = 5 × 10 −5 Pa s −1 ). All these parameters can be adjusted via a name list. The mentioned default values of α, ν 1 and ν 2 have been optimised for a 25 km grid resolution. At lower resolution, α can be increased and ν 1 and ν 2 reduced to speed up the algorithm.
In the Zwack-Okossi equation, geostrophic vorticity tendencies are converted to geopotential height tendencies using Eq. (13), which is a 2-D Poisson's equation. To solve this equation computationally effectively, we utilise Intel's MKL (Math Kernel Library) fast Poisson solver routines, which employ the DFT (discrete Fourier transform) method. Intel's MKL is widely available and free to download, although registration is required.

Boundary conditions
In the omega equation, homogeneous boundary conditions (ω = 0) at both the meridional and the lower and the upper boundaries are used. For the Zwack-Okossi equation, a slightly more complicated procedure is used to ensure that the area means of the individual height tendency components are consistent with the corresponding temperature tendencies. First, for all of V , T , F , Q and A, the height tendency is initially solved from Eq. (13) using homogenous boundary conditions ( ∂Z ∂t = 0) at the northern and southern boundaries. Then, the area mean temperature tendencies for these five terms are calculated, taking into account both the direct effect of temperature advection (for T ) and diabatic heating (for Q) and the adiabatic warming-cooling associated with the corresponding omega component. These temperature tendencies are substituted to the hypsometric equation to calculate the corresponding area mean height tendencies. In the vertical integration of the hypsometric equation, the area mean height tendency at the lower boundary (1000 hPa) is set to 0, following the expectation that the total atmospheric mass in the model domain is conserved. Horizontally homogeneous adjustments are then made to the initial height tendencies for V , T , F , Q and A, forcing their area means to agree with those derived from the hypsometric equation.

The WRF set-up
WRF is a non-hydrostatic model and can generate atmospheric simulations using real data or idealised configurations (Wang et al., 2007;Shamarock et al., 2008). The calculations presented in this paper used input data from an idealised moist baroclinic wave simulation, which simulates the evolution of a baroclinic wave within a baroclinically unstable jet in the Northern Hemisphere under the f-plane approximation (Blázquez et al., 2013). The value of the Coriolis parameter was set to 10 −4 s −1 in the whole model domain.
The simulation presented in Sects. 6-7 was run for 10 days with a 30 min output interval, in a domain of 4000 km×8000 km×16 km, with 25 km horizontal grid spacing. The horizontal grid was Cartesian and staggered, and 64σ levels were used in the vertical direction. The boundary conditions were symmetric in the meridional direction and periodic in the zonal direction. Cloud microphysics is parameterised using the WRF single-moment 3-class scheme (Hong et al., 2004). Cumulus convection is treated with the Kain-Fritsch scheme (Kain and Fritsch, 1993) and boundary layer turbulence with the YSU scheme (Hong et al., 2006). The radiation scheme was switched off in our model integration.
After running the simulation, data were interpolated from model levels to 19 evenly spaced pressure levels (1000, 950, . . . , 100 hPa). The interpolation was done with WRF utility wrf_interp, which is freely available from the WRF website. During the interpolation, the horizontal data grid was unstaggered to mass points, thus having 160 grid points in the zonal and 320 grid points in the meridional direction.
The model output data contained all the variables needed in solving the generalised omega equation and the Zwack-Okossi equation: temperature, wind components, geopotential height, surface pressure, and parameterised diabatic heating and friction components. Diabatic heating and friction in WRF included contributions from various physical processes, such as cumulus convection, boundary layer physics and microphysics. The physical tendencies are not in the default WRF output, and need to be added by modifying the WRF registry file.
To study the performance of OZO at a resolution that is more similar to that used in several earlier diagnostic studies of synoptic-scale dynamics (Räisänen, 1995(Räisänen, , 1997Caron et al., 2006;Stepanyuk et al., 2017), WRF was also run in the same domain at 100 km grid spacing. Results for this simulation are presented in the supplementary material. An intermediate simulation at 50 km resolution was also made to study the resolution dependence of the computational performance. The inversion of the left-hand side operator of the omega equation (Eq. 6) is computationally quite a heavy process. Hence, in our previously described test runs, the calculation of the five plus two ω components took almost all of the total computing time. For the height tendencies, the inversion of the horizontal Laplacian (Eq. 13) is much more straightforward to do and, moreover, MKL fast Poisson solver routines are employed. For that reason, only a small fraction of the computing time is spent for the height tendency calculation. Table 2 provides information on the absolute computing times and how they depend on the number of grid points in the model domain. It is our goal to improve the computational performance by utilising better scalable solver routines for the omega equation in a future version of the software.
6 Results -vertical motion 6.1 Comparison between calculated and WRF-simulated vertical motions Figure 1 compares the solution of the generalised omega equation (ω TOT ) with ω as being obtained directly from the WRF output (ω WRF ), at the 700 hPa level after 118 h of simulation when the cyclone is approaching its maximum intensity. The agreement is excellent. A strong maximum of rising motion (ω ≈ −3 Pa s −1 ) occurs near the occlusion point to the east of the surface low in both cases, with a somewhat weaker ascent along the cold frontal zone to the south-west and in the north-eastern sector of the low. Descent prevails further east of the low and behind the cold front. However, lots of mesoscale details are visible in both the simulated and the calculated ω, although with some not perfect agreement between these two. The difference between ω TOT and ω WRF (Fig. 1c) is noisy, although it suggests slightly more discrepancy behind the cold front where shallow convection takes place.
The QG omega equation (Eq. 1) also captures the largescale patterns of rising and sinking motion reasonably well (Fig. 1d). However, it substantially underestimates the ascent to the east of the low and along the cold front, and there are two bands of strong descent (behind the cold front and to the east of the low) that are much weaker in ω WRF and ω TOT . Furthermore, many of the mesoscale details shared by ω WRF and ω TOT are missing in the QG solution.
The majority of the differences between ω TOT and ω WRF most likely result from numerical errors. One source of these errors is the approximation of the time derivatives in F A in Eq. (11) with central differences. As shown in Fig. 2, the agreement between ω TOT and ω WRF gradually improves when the half-span of these time differences is reduced from 2 h to 1 min (i.e. the time step in the WRF simulation). The 30 min output interval used for the other calculations shown in this paper is a compromise between accuracy and the need to control the data volume of the WRF output for the 10-day simulation.
A more comprehensive statistical evaluation of the calculated vertical motions is given in Figs. 3 and 4, by using data from the whole model area and the 8 last days of the simulation. The first 2 days, when both the cyclone and the vertical motions are still weak, are neglected. Figure 3 shows the time-averaged spatial correlation between ω WRF and various omega equation solutions. The correlation between ω TOT (line VTFQA) and ω WRF is excellent, reaching 0.95 in the mid-troposphere and exceeding 0.88 at all levels from 250 to 850 hPa. However, leaving out ω A , which requires non-synoptic information from the time derivatives of temperature and vorticity, deteriorates the correlation substantially (line VTFQ). The solution of the QG omega equation only correlates with ω WRF at r ∼ 0.6 in the mid-troposphere (line VT(QG)), although the correlation approaches 0.65 at 300 hPa where diabatic heating is less important than at lower levels. Caron et al. (2006) reported excellent correlations (at most levels, from 0.9 to 0.96) between the vertical motions calculated by the DIONYSOS tool and their original model output (their Figs. 2a and 3a), for two cases and nine synoptic times at 3 h intervals for both. These correlations are similar to or even higher than those shown in Fig. 3. However, in these tests the horizontal resolution of DIONYSOS was 100 km, whereas we used 25 km resolution in our WRF simulation. In fact, the correlations for OZO are also improved and reach 0.97 in the mid-troposphere for a 100 km resolution WRF simulation (see the Supplement, Fig. S3). Naturally, the performance may also depend on the synoptic case studied. Figure 4 compares the root-mean-square (RMS) amplitudes of ω WRF and ω TOT . RMS(ω TOT ) is typically about 5 % smaller than RMS(ω WRF ), which is presumably due to the truncation errors that occur when the derivatives in the forcing terms are approximated with second-order central differences. An exception to this is the lower troposphere where vertical motions calculated by OZO(ω TOT ) are slightly stronger compared to the direct WRF output (ω WRF ). The RMS amplitudes of the individual ω components will be discussed in the next subsection. Figure 5 shows the contributions of the five individual forcing terms to ω(700 hPa) for the situation studied in Fig. 1. Vorticity advection and thermal advection both contribute substantially to the vertical motions ( Fig. 5a-b), but the maximum of ascent is further east for ω T (Fig. 5b) than ω V (Fig. 5a). Due to this phase shift, there is a cancellation between rising motion from ω V and sinking motion from ω T behind the cold front. This cancellation effect is well-known and typically occurs just behind the cold front, where cold advection and increasing cyclonic vorticity advection with height take place (Lackmann, 2011). On the other hand, these two terms both induce rising motion to the south-east of the centre of the low. Diabatic heating, which is dominated by latent heat release, strongly enforces the ascent along the main frontal zones of the cyclone (Fig. 5d). Compensating subsidence prevails in the surrounding areas, except for spots of localised ascent associated with convective precipitation well behind the cold front. The imbalance term ω A is remarkably large but noisy. Friction induces ascent close to the centre of the low and descent around and to the northeast of the surface high (Fig. 5c), but its contribution is quite weak at the 700 hPa level. Figures S1 and S2 in the Supplement present the actual forcing fields for the vertical motions fields in Fig. 5a-e.

Contributions of individual forcing terms to ω
In terms of the RMS amplitudes evaluated over the whole model area and the last 8 days of the simulation, temperature advection makes the largest contribution to the calculated ω (line T in Fig. 4). Vorticity advection (line V ), diabatic heating (line Q) and the imbalance term (line A) are all similar in magnitude in the mid-and upper troposphere. Near the surface, the imbalance term has a tendency to compensate the effects of temperature advection, as their RMS values are even larger than the RMS of the total vertical motion (line TOT). RMS(ω Q ) peaks at 850 hPa, which is mostly due to shallow convection behind the cold front (not shown). This peak is visible also in the RMS(ω TOT ). Our sensitivity experiments with different output intervals suggest that the half-hour tem-    poral resolution is too sparse for proper estimation of the imbalance term in the presence of fast-moving convection cells, which causes the over-estimation of total vertical motion in the lowest troposphere. RMS(ω F ) is at its maximum near the top of the boundary layer at 900 hPa but remains weak even at this level. These results are partly consistent with similar calculations made for observation-based analysis data (Räisänen, 1995) and for model data (Stepanyuk et al., 2017) at lower spatial resolution. However, the imbalance term was relatively small in the study by Räisänen (1995), because the 6 h time resolution of its input data was not sufficient for a proper estimation of this term. In addition, RMS(ω T ) is more dominant in our study than in Räisänen (1995) and in Stepanyuk et al. (2017). The main reason for this is the improved horizontal resolution. Räisänen (1995) found that the influence of thermal advection compared to vorticity advection increases when the horizontal resolution is improved. This is because vertical motions induced by thermal advection typically have a smaller horizontal scale than those induced by vorticity advection (Räisänen, 1995). In this paper, the horizontal resolution is almost 10 times higher than in the study by Räisänen (1995). Conversely, RMS(ω F ) is smaller in Fig. 4 than found for the midlatitudes in Räisänen (1995) and Stepanyuk et al. (2017). This may be at least in part because Räisänen (1995) and Stepanyuk et al. (2017) included both land and sea areas, whereas our WRF simulation was made for an idealised ocean surface.
A further division of ω V and ω T to contributions from advection by rotational and divergent wind reveals that they both are largely dominated by the rotational wind (not shown).

Comparison between calculated and WRF-simulated height tendencies
In this subsection, the calculated total height tendencies are compared with height tendencies from the WRF simulation. The latter were estimated as central differences from the 30 min time series of the simulated geopotential heights. Figure 6 shows the distributions of the calculated height tendency, the WRF height tendency and their difference slightly before the cyclone reaches its maximum intensity (t = 118 h). The values are shown at the 900 hPa level, which is sufficiently low to represent the processes affecting the low-level cyclogenesis. Negative (positive) height tendency  over the low centre indicates deepening (weakening) of the low. In general, a very close agreement between the fields is seen, although a somewhat more positive bias is visible behind the cold front. Figure 7 shows the time-averaged RMS error and correlation coefficient between the calculated and WRF height tendency as a function of height. The RMS error is quite constant from 950 hPa up to the 250 hPa level (Fig. 7a). Above this level, the error grows rapidly towards the stratosphere. This error growth is accompanied by a decrease of correlation coefficient at the same altitude. This deterioration is presumably at least partly due to the 50 hPa vertical resolution in the OZO, which is too coarse for an adequate representation of stratospheric dynamics.
The correlation between the calculated and WRF height tendency is highest in the upper troposphere, which is roughly 0.97. The correlation weakens slightly closer the surface, but still exceeds 0.95. Thus, the calculated height tendency is generally in very good agreement with the tendency diagnosed directly from the WRF output. These correlations are comparable but mostly slightly higher than those reported by Caron et al. (2006) for DIONYSOS (their Fig. 6a).

Contributions of individual terms during the mature stage
The contributions of the individual height tendency components at the 900 hPa level at 118 h are shown in Fig. 8. Vorticity advection (Fig. 8a) produces a wide and strongly positive height tendency behind the surface low (see Sect. 7.4 for further analysis of this term). Thermal advection (Fig. 8b) causes a positive height tendency in the area behind the surface low and a negative height tendency at the opposite side. This large negative height tendency ahead of the low is caused by warm air advection in the mid-and upper troposphere. In this baroclinic life cycle simulation, thermal advection is the main contributor to the movement of the cyclone, which is in agreement with the study of Räisänen (1997). Friction (Fig. 8c) always acts to damp synoptic-scale weather systems and is thereby inducing a positive (negative) height tendency over the surface low (high). Diabatic heating (Fig. 8d) is causing uniformly negative lower tropospheric height tendencies in the vicinity of the surface low. The largest negative height tendency due to diabatic heating is located south-east from the low centre, where strong latent heat release occurs in connection with frontal precipitation.
The imbalance term (Fig. 8e) shows more small-scale structure than the other terms. In general, however, it is in phase with the total height tendency near the centre of the low, with negative values to the east and positive values to the west. The reason for this feature is most probably the following. The conversion from geostrophic vorticity tendencies to height tendencies for the other terms was done by assuming a geostrophic balance according to Eq. (13). However, in cyclones the wind is typically sub-geostrophic (e.g. Holton and Hakim, 2012). Therefore, the tendency of geostrophic vorticity exceeds the actual vorticity tendency. This implies that height tendencies calculated from the actual vorticity tendency under the geostrophic assumption will be too small. The imbalance term takes care of this and makes the calculated total height tendency to correspond better to the actual change of the geopotential height field. Figure 9 shows the 900 hPa level height tendencies induced by the five individual terms in the cyclone centre during the deepening period. The low deepens vigorously roughly between 72 and 120 h of simulation, as shown by the negative total height tendency (black line) during this period. The total height tendency is also in a good agreement with the WRF height tendency (dotted line). The deepening is mostly due to vorticity advection (blue line) and the imbalance term (grey line). Later on, roughly from 96 h onward, diabatic heating (red line) and thermal advection (orange line) together with the imbalance term make the largest contributions to maintaining the intensity of the surface low. Friction (green line) systematically destroys cyclonic vorticity over the cyclone centre and thus produces a slightly positive height tendency during the whole life cycle. After circa 96 h, vorticity advection also acts a damping mechanism for the surface low (see also the next subsection). . Time series of individual height tendency components at the 900 hPa level from the cyclone centre during the deepening period. V is vorticity advection, T is thermal advection, F is friction, Q is diabatic heating, A is imbalance term, Tot is total and WRF is WRF height tendency. Height tendencies at the cyclone centre were averaged over all grid boxes in which the 900 hPa geopotential height was less than 5 m above its minimum. In addition, the values are 2.5 h moving averages. 7.4 The effect of vorticity and thermal advection by divergent and non-divergent winds Figure 10 presents the height tendencies associated with vorticity advection by V χ and V ψ separately. The divergent wind vorticity advection causes widespread and strong positive height tendency over and around the surface low (Fig. 10a). According to Räisänen (1997), the divergent wind transports anticyclonic vorticity from the surroundings of the surface low, and is thus acting to reduce the cyclonic vorticity at the centre of the low. In the case of the non-divergent wind component (Fig. 10b), positive (negative) height tendencies behind (ahead of) the low originate from the upper troposphere, where the non-divergent wind and thereby vorticity advection is the strongest. Cyclonic vorticity advection ahead of the trough produces negative height tendency in the same area, while anticyclonic vorticity advection ahead of the ridge does the opposite. Furthermore, the non-divergent vorticity advection is substantially contributing to the deepening of the cyclone, since the area of the negative height tendency reaches the centre of the low as well.

Height tendencies in the cyclone centre
In contrast to vorticity advection, thermal advection by divergent winds was found to cause a negligible height tendency (not shown).

Limitations and plans for future developments
The idealised baroclinic wave simulation provides an effective and widely used tool for studying cyclone dynamics in an easily controlled model environment. For this reason, we chose to begin the development work of OZO from a relatively simple Cartesian implementation. Nevertheless, the idealised Cartesian geometry obviously reduces the number of potential applications of OZO. We aim to extend OZO to more complex spherical grid applications in the future. However, in principle the use of OZO is possible with some limited-area real cases where spherical geometry has been used, if the data are regridded to Cartesian geometry afterwards. However, in this case, one must change the lateral boundary conditions of OZO since they are tailored for periodic model domain by default.
Another limiting factor in the current version of OZO is the weak scalability in the multigrid omega equation solver together with the lack of parallelisation of the source code. For high-resolution runs, this means significant slowing of the calculations (Table 2). To reduce the lengthy computing times, we are currently developing a parallel version of OZO, which uses a different solving method for the omega equation. We aim to release this parallel version by the end of the year 2017.
An issue associated with the physical interpretation of high-resolution simulations should also be noted. In OZO, the vertical motion and height tendency contributions of vorticity advection, thermal advection, friction, and diabatic heating are calculated assuming geostrophic balance between the vorticity and temperature tendencies, while the deviations from this balance are attributed to the imbalance term. As the balance assumption is increasingly violated at higher resolution, the imbalance term grows larger. This complicates the interpretation of the results particularly when the imbalance term opposes the other terms. Such compensation is indeed apparent in our results for vertical motions in the lower troposphere, where the contributions of thermal advection and diabatic heating are strongly opposed by the imbalance term (Fig. 4). At resolutions much higher than 25 km, issues like this are likely to become increasingly problematic.

Conclusions
In this paper, a software package called OZO is introduced. OZO is a tool for investigating the physical and dynamical factors that affect atmospheric vertical motions and geopotential height tendencies, tailored for WRF simulations with idealised Cartesian geometry. As input to OZO, the output of the WRF model interpolated to evenly spaced pressure levels is required.
The generalised omega equation diagnoses the contributions of different physical and dynamical processes to vertical motions: vorticity advection, thermal advection, friction, diabatic heating, and imbalance between temperature and vorticity tendencies. Then, analogously with the vertical motion, the height tendencies associated with these forcings are calculated. As an advance over traditional applications of the Zwack-Okossi equation (Zwack and Okossi, 1986;Lupo and Smith, 1998), the use of the generalised omega equation allows OZO to eliminate vertical motion as an independent forcing in the calculation of height tendencies.
The calculated total vertical motions and height tendencies in the test case are generally in excellent agreement with the vertical motions and height tendencies diagnosed directly from the WRF simulations. The time-averaged correlation between the calculated and the WRF height tendency was 0.95-0.97 in the troposphere. For the vertical motion as well, a correlation of 0.95 was found in the mid-troposphere. Our analysis further illustrates the importance of both adiabatic and diabatic processes to atmospheric vertical motions and the development of the simulated cyclone.
The OZO software is applicable to different types of WRF simulations, as far as Cartesian geometry with periodic boundary conditions in the zonal direction is used. One example of potential applications are simulations with increased sea surface temperatures as the lower boundary condition. Combined with OZO, such simulations provide a simple framework for studying the changes in cyclone dynamics in a warmer climate.

Data and code availability
The source code of OZO is freely available under MIT licence in GitHub (https://github.com/mikarant/ozo). OZO v.1.0 described in this manuscript is also archived at http: //doi.org/10.5281/zenodo.157188. In addition to the source code, the package also includes a makefile for compiling and running the program, a small sample input dataset for testing the functionality, and two README files containing the instructions for both generating input data with WRF and running the OZO program. The WRF model as well as the interpolation utility (wrf_interp) are downloadable from the WRF users page (http://www2.mmm.ucar.edu/wrf/users/ downloads.html). Intel's MKL can be downloaded after registration from their web page (https://software.intel.com/ en-us/articles/free-mkl). OZO v1.0 is guaranteed to work with WRF version 3.8.1.
The Supplement related to this article is available online at doi:10.5194/gmd-10-827-2017-supplement.