TRAM: A new non‐hydrostatic fully compressible numerical model suited for all kinds of regional atmospheric predictions

A new limited‐area numerical model (TRAM, for Triangle‐based Regional Atmospheric Model) has been built using a non‐hydrostatic and fully compressible version of the Navier–Stokes equations. Advection terms are solved using a Reconstruct–Evolve–Average (REA) strategy over the computational cells. These cells consist of equilateral triangles in the horizontal. The classical z‐coordinate is used in the vertical, allowing arbitrary stretching (e.g., higher resolution in the Planetary Boundary Layer, PBL). Proper treatment of terrain slopes in the bottom boundary conditions allows for accurately representing the orographic forcing. To gain computational efficiency, time splitting is used to integrate fast and slow terms separately and acoustic modes in the vertical are solved implicitly. For real cases on the globe, the Lambert map projection is applied, and all Coriolis and curvature terms are retained. No explicit filters are needed. The first part of the manuscript describes the dynamical core of the model and provides its thorough validation using a variety of benchmark tests (mostly in two dimensions) in the context of a dry‐adiabatic atmosphere. In the second part, TRAM is reformulated for a moist atmosphere and is completed with a proper set of physical parametrizations of cloud microphysics, cumulus convection, short and long‐wave radiation, PBL processes and surface fluxes. Various examples of the great versatility offered by this full version will be presented, with special emphasis on Mediterranean case studies. In summary, TRAM performs as well as state‐of‐the‐art numerical models and is suitable for simulating circulations ranging from small‐scale thermal bubbles (≈100 m scale) to synoptic‐scale baroclinic cyclones (>1000 km size), including orographic circulations, thermally driven flows, squall lines, supercells, all kinds of precipitation systems and medicanes. Besides opening a myriad of academic and research applications, TRAM regional forecasts at different resolutions are being disseminated in the web (see https://meteo.uib.es/tram).


INTRODUCTION
When numerically modelling the atmosphere for both research and weather forecasting, it has often been desirable to work with approximate versions of the governing equations.By filtering out certain kinds of motions, these approximate forms will be more tractable and easier to solve numerically (Lauritzen et al., 2011).Two of the most common approximations (e.g., Durran, 1999) are the hydrostatic approximation (only appropriate for horizontal scales greater than about 10 km) that neglects the acceleration term in the vertical momentum equation, and the anelastic approximation (it responds well even on scales smaller than 10 km) that neglects the elasticity of the fluid by using a special form of the mass continuity equation.Among other effects, the anelastic equations do not support acoustic modes and the hydrostatic equations do not support internal acoustic modes.Despite the very weak energy of these modes, their high frequency would make it expensive or complicated to retain them in the model.The multiscale nature of atmospheric dynamics, ranging from planetary waves to boundary-layer turbulence, needs to be represented as best as possible in a numerical model when pursuing accurate weather predictions at all scales.Consider, for instance, a critical but highly elusive surface variable as precipitation.Processes that lead to rainfall are linked to atmospheric circulations as diverse as synoptic-scale cyclones, fronts, thunderstorms and small cumuli (Hayhoe et al., 2017).Unfortunately, neither the hydrostatic nor the anelastic approximation is valid on all horizontal scales.Consequently, state-of-the-art atmospheric models, designed to work from global scales down to sub-kilometre scales, use the non-hydrostatic fully compressible (NHFC) equations (Lauritzen et al., 2011).
On the other hand, and continuing with our example, generation of precipitation-size drops will ultimately rely on complex microphysical transformations in each individual cloud.These clouds will evolve in association with the explicit dynamics but also through mutual interactions with other subgrid-scale physical processes, such as radiative forcing, Planetary Boundary Layer (PBL) mixing and surface fluxes (Stensrud, 2009).Thus, weather predictions will not only benefit from a formulation of the model dynamical core using some version of the primitive NHFC equations, but also from any improvement in the realism of the physical parametrizations and their complex interactions.Progress in both facets with a view to convection-permitting resolutions does not come without a price, namely the need of special numerical methods and an increase in computation requirements, but it does mean a stimulating challenge, intellectually and scientifically, that would certainly pay off in terms of improved performance in numerical weather prediction (NWP) and its applications.This is the main driver of our research proposal.
Development and maintenance of ambitious numerical models like the one pursued here are prevalently framed within big institutional projects dealing with fundamental research or practical applications, including commercial use.But occasionally there have been modest (in terms of size) research units or university groups with outstanding contributions in this field.For instance, a meteorology group at CNR-ISAC (Italy) of similar proportions to ours, is internationally known for developing the very versatile suite of GLOBO-BOLAM-MOLOCH models (see http://www.isac.cnr.it/dinamica/projects/forecasts/; Davolio et al., 2020 and references therein).In these contexts, models are originally devised and coded mainly for research and academic purposes.Without the pressure imposed by real-time applications, simulation of mesoscale flows over a region of interest and design of ad hoc idealized experiments can be naturally addressed, even at very high grid resolutions.Over time and at the sacrifice of some resolution, these modelling systems easily migrate towards forecasting tasks as well, provided they are suited to the spatial scales of interest (at least synoptic and mesoscale).
Our TRAM model was built from scratch, and this required a sequential completion of steps, from the formulation of the advection scheme at the very beginning (first in one dimension and then in two and three dimensions) to the inclusion and mutual coordination of the physical parametrizations as the last step.The following sections basically reproduce the actual chronology of milestones reached during TRAM development, describing with a didactic approach the sequence of conceptual/methodological building blocks of the model along with the corresponding validation tests (many of them extracted from the literature).A new model must not only successfully pass benchmark tests but also involve original or uncommon aspects in its formulation, such as our horizontal discretization of the equations using triangles (this giving rise to the acronym TRAM, for 'Triangle-based Regional Atmospheric Model').I will stress these and other innovative aspects in the description while also justifying other important choices of the formulation based on the experience with companion NHFC numerical models (e.g., MM5, Dudhia, 1993, Grell et al., 1995;WRF, Skamarock et al., 2008;andCM1, Bryan, 2002, Bryan &Fritsch, 2002).
The paper is organized in two parts as follows: part 1 focuses on the dynamical core of the model for a dry-adiabatic atmosphere, describing the horizontal mesh and advection scheme (Section 2.1), the time integration scheme (Section 2.2), the basic equations of the model (Section 2.3), the handling of orography and lateral boundary conditions (Section 2.4), the vertical stretching of the height coordinate and consequential use of a semi-implicit scheme (Section 2.5), and validation tests not included in previous sections (Section 2.6); part 2 focuses on the adaptation of the model to the realistic and moist atmosphere, and therefore presents the full form of its equations with details of the included physical parametrizations (Section 3.1), the description of the time-marching algorithm (Section 3.2), and a wide range of validation tests, structured as mesoscale-idealized simulations (Section 3.3) and synoptic-real-case simulations (Section 3.4).As is customary, a last section will present the main conclusions of the work, as well as ideas for future improvement and extension of TRAM capabilities.

Horizontal mesh and advection scheme
Advection terms are the genuinely nonlinear elements of any set of geophysical flow equations and deserve special treatment to avoid spurious or poor numerical solutions.For the moment we start with the linear advection equation in its simplest form, through modelling the transport of a scalar quantity q along the one-dimensional domain under the action of a prescribed current of constant value (e.g., a westerly wind ū):  t q = −ū x q.This equation can be naturally expressed in flux form as  t q = − x ( ūq).
Our advection scheme follows a Reconstruct-Evolve-Average (REA) philosophy, inspired by the finite-volume methods point of view (see Leveque, 2002).In the classical formulation of a finite-volume method, the latter equation would be integrated forward in time over a time step Δt at each grid cell of size Δx by approximating in some way -as function of known information at surrounding grid points -the time-integrated fluxes of field q at right and left boundaries of the cell.But a time integration with formulation of fluxes is equivalent to the application of a REA algorithm with piecewise linear reconstruction (Leveque, 2002), as follows: (i) for each grid cell or 'volume' (in fact, a segment in one dimension), a linear profile of the field is Reconstructed using the central grid point value and an adequate choice for the slope or gradient in that direction (R-step); (ii) this profile Evolves conservatively over time interval Δt (i.e., the segment is translated along the domain, and if wind was not spatially constant, also stretches/shrinks) by the action of the velocity found at segment edges (E-step); and (iii) the new grid point values from the advection process are calculated by Averaging within each cell's limits all segments (i.e., profiles) that totally or partially remain or enter the cell after the previous evolutions (A-step).
What exactly defines the numerical behaviour of the REA scheme is the particular choice of the slope for the R-step.For instance, a zero slope yields a scheme equivalent to the first-order upwind method, known for correctly preserving the monotonicity of the solution but inducing an unacceptable degree of dissipation.Slopes consistent with popular second (or higher)-order methods can also be defined, but then other characteristic problems of these schemes, like phase error and oscillatory contamination of the solution, can easily emerge.The previous pathologies can be greatly reduced by using slope limiter methods (a synonym term for the flux-limiting strategy distinctive of the high-resolution finite-volume methods, Leveque, 2002).Among the large family of these specially defined -and local solution dependent -slopes, and guided by several advection tests (see two-dimensional tests at the end of this section), we found the monotonized central-difference limiter (MC limiter;Van Leer, 1977) as the most suitable choice.The MC slope at each cell i for a given time step n is defined as the flattest slope (or zero if not all three of the same sign) of these three values: the centred slope (Q n i+1 − Q n i−1 )/2Δx, double the upstream slope 2(Q n i − Q n i−1 )/Δx, and double the downstream slope 2(Q n i+1 − Q n i )/Δx.(Note we use capital letters to refer to the discrete numerical solution.) Exactly the same REA approach can be applied to the more general, non-flux form of the equation where the advective wind is spatially heterogeneous ( t q = −u(x) x q) or even part of the pursued numerical solution: u(x,t).A good -but still simple -prototype of this type of nonlinear complexities, characteristic of the partial differential equations for the atmosphere, is found in the well-known Burgers' equation (Burgers, 1948):  t u = −u x u.Specifically, what is needed for the scheme is some reasonable assignment of velocity values at the edges, U n i−1/2 and U n i+1/2 , at the intermediate E-step.The most obvious choice is the average of the two neighbouring values, such as )/2 at the eastern edge, leading to a displacement of length U n i+1/2 Δt.However, the displacement method would better account for the spatial variabilities of the wind field if, instead of using a constant speed for the edge point, we incorporate a linear profile to its speed, leading to solving the kinematic equation ̇x = U n i+1/2 + Ax, where A is the wind gradient found at the cell penetrated by the edge (in our example, that would be cell i + 1 for a westerly wind, cell i for an easterly wind).This differential equation has a simple analytical solution for the total path covered by the edge point during interval Δt, providing a better estimate than the previous solution with A = 0. Indeed, multiple advection tests (illustrative two-dimensional examples are shown below) proved the benefits of this kind of dynamic treatment of the E-step, especially for large time steps.
Before proceeding with two and three dimensions, note that the formulation of the advection scheme assumes no staggering of variables; that is, all intervening fields (winds and scalars) are discretized over a common mesh of points.This is the case, for the moment, of our TRAM model.Moreover, without further criteria or limitation on the MC slope, strictly speaking the scheme would not be positive definite for naturally discontinuous variables such as water species; but in practice, this problem is almost insignificant owing to the strong non-oscillatory character of the scheme near acute gradients of the advected field.The last comment concerns the numerical stability of the scheme.The Lagrangian character of the E-step, where the field values along the segments are conserved during their movement, guarantees the stability of the results as long as the final A-step for any given cell is effectively completed for all possible segments intruding into that cell.This is no more than the Courant-Friedrichs-Lewy (CFL)-type condition for the stability of forward schemes applied to hyperbolic problems (Leveque, 2002): the domain of dependence of the numerical scheme (range of scanned neighbouring cells in our case) must contain the true physical domain of dependence (maximum possible distance travelled over time Δt).Nevertheless, for an easier coding of the advection algorithm and to prevent the loss of accuracy for larger time steps, we will limit Δt to the maximum value that guarantees that all cells are only influenced by the immediately adjacent neighbours.A flexibilization of this requirement will be introduced later, just for the vertical, when considering segments of varying size (Section 2.5).
The advection scheme would naturally extend to the three-dimensional space by implementing the so-called dimensional splitting (Leveque, 2002), where the method is sequentially applied in the three dimensions; for example, first X, then Y and finally Z.However, in recognition of the observed kinematic properties at many spatial scales, characterized by strongly stratified flows dominated by magnitudes, divergences, vorticities and deformations of the wind predominantly along the horizontal plane (see Bluestein, 1992), rather than acting independently on the X and Y directions the model will solve horizontal advection in-bloc.By limiting the dimensional splitting to only two steps (first horizontal and then vertical) the splitting error of the standard three-steps method would also be attenuated.
Figure 1 illustrates two important and unique characteristics of the TRAM model.In the first place, the topology of grid cells in the horizontal plane, structured as equilateral triangles of two types (T, with a vertex at the top; B, with a vertex at the bottom) filling the space.All model fields are defined at the barycentres of T and B. One could argue that the horizonal resolution of the mesh is defined by the side length Δx of the triangles, but referring to the traditional square-based mesh, this resolution could equivalently be identified by 2/3 Δx (i.e., the side length of a square covering approximately the same area as T/B).Note that the two-dimensional gradient of a field at a point T is naturally defined by the value of the field at the three points B bordering the cell in its NW, NE and S directions.Analogously, any gradient at a B cell is provided by the three surrounding T points located to its N, SW and SE.
In the second place, Figure 1 summarizes the implementation of the REA approach on the triangle-based mesh for solving the horizontal advection in-bloc.(Vertical advection implementation is no longer discussed as it is fully equivalent to the one-dimensional problem discussed above, because our model uses geometric height as vertical coordinate; thus, volumetrically the model uses triangular prismatic cells of constant thickness Δz, the parameter that defines the vertical resolution.) As a generalization of the one-dimensional scheme, for the R-step the MC slope (two-dimensional gradient in this case) is determined for each T/B cell by comparing the magnitude of the gradient at the cell with the magnitudes -once doubled -at the upgradient and downgradient cells.These up-and downgradient opposing cells will be two of the three surrounding B/T cells respectively.Remember that the MC criterium will choose the two-dimensional gradient of lowest magnitude among the three compared vectors, but this gradient will be zero if these vectors do not project positively onto each other.Again, as an extension of the one-dimensional algorithm, in the E-step the planes reconstructed for each cell (bluish and reddish areas in the example of Figure 1, originally coincident with a pair of adjacent T/B cells) will evolve advectively over time Δt and affect their neighbours.By mimicking the one-dimensional scheme, this evolution is 'dynamic', as illustrated by the curved green arrows in Figure 1.That is, rather than using the six-cell wind vector averages at triangle corners and shift the vertices in a straight line, we solve the two-dimensional kinematic equations for the motion of the vertices subjected to a linearly dependent wind field (wind field gradients are likewise extracted from the penetrated cells).Finally, the A-step for the advected field is completed by contributing to the new cell averages at time t + t.In the schematic example, the T cell will contribute to eight neighbour cells of the 13 possible (including itself), while the B cell will affect five of the 13 possible cells (Figure 1).
A validation of the described advection scheme was conducted for diverse experimental designs, involving both uniform and spatially variable winds and checking the behaviour for both smooth and sharp gradients of the advected field.In the two-dimensional example of Figure 2 we display the results of a rotating field, a common test for new algorithms (see for instance the overview of Smolarkiewicz, 2005).A solid-body rotational wind, with the axis over the southern Iberian Peninsula and an angular velocity of one revolution/day, is imposed on a scalar field built by adding a main gaussian with two others of half width shifted to its east and southwest (see these profiles in Figure 2a).For this particular test the size of the triangular cells is Δx = 20 km, the time step is Δt = 60 s and we do not incorporate the 'dynamic' improvement of winds in the E-step.The excellent behaviour of the MC-based method is observed after four days of simulation in Figure 2d, clearly contrasting with the much poorer performance of a variant of the REA approach where a zero slope is chosen in the R-step (essentially the first-order upwind method, Figure 2b) or where this slope is systematically taken from the downstream cell (a variant corresponding to the second-order Lax-Wendroff method; Figure 2c).The benefits of activating the 'dynamic' calculation of the triangle's new arrangements for the E-step are evidenced by a significant reduction of simulation error; specifically, Figure 3 intercompares the error of the previous MC slope experiment with the error of a theoretically degraded simulation (according to the integration time step, tripled to 180 s) but which is actually better owing to the 'dynamic' treatment of winds.
Similar conclusions about the good performance of the TRAM advection scheme could be drawn for the rotating cylinder, highlighting the lack of significant oscillations or negative values around the structure after several revolutions, as well as an arrangement of the simulation error basically axisymmetric and confined to the cylinder walls (experiment not shown).These kinds of analyses were also extended to the three-dimensional framework, requiring the coordination of horizonal and vertical REA methods as described above.Finally, before coupling the advection scheme to more complex systems of equations (next sections) the stability and fluent functioning of the scheme in the nonlinear context were also confirmed based on experiments with the two-dimensional Burgers' equation (e.g., evolution of a vortex embedded in zonal flow; not shown).It should be noted that the REA scheme carries a small but perceptible amount of computational diffusion, which for complex applications should be considered a beneficial property towards maintaining small-scale noise under control.

Time integration scheme (shallow-water model)
A good prototype of more general atmospheric flow equations is offered by the shallow-water (SW) model (Vreugdenhil, 1986).Without loss of generality, we can consider that we are modelling phenomena such as tsunami propagation, accurately described by a 'thin' layer of fluid of constant density in hydrostatic balance.The corresponding nonlinear SW equations over a flat bottom topography are: where h is the height of the free surface and (u,v) the two-dimensional velocity.Note we are not expressing the equations in any of the admitted conservation forms and also the fact that differential operators are fully developed into separate terms; the reason is to conform to the norm of the TRAM model.Despite its simplicity, the SW model captures the nonlinear interplay between mass and momentum variables, very much like it is found in the full atmospheric equations (next section).Additionally, we incorporate the possible role in the experiments of the rotating earth (f is the Coriolis parameter) and frictional/diffusion forces (through coefficients b and ).Like the gravity (g), all these parameters, if activated, are constant over the domain.
The basic terms of the SW dynamics are the advection terms (marked in red) and the 'forcing' terms responsible for the propagation of surface gravity waves (marked in green/magenta).Typically, flow velocities would reach a few m⋅s −1 at most, while gravity wave speed is about √ gH (H is the representative depth of the water (that is, its difference with respect to the true solution of Figure 2a); (b) numerical error for the same type of simulation but increasing the time step to Δt = 180 s while incorporating the 'dynamic' treatment of winds.In both maps the contour interval is five units, starting at −2.5 for negative values (dashed line) and at +2.5 for positive values (continuous line).
column).This distinction between 'slow' and 'fast' terms respectively, would allow a time-splitting strategy (see next section) but, for the moment, a single time step Δt is used for all terms.Ruled by the CFL condition for stability, the maximum operational value for Δt is set by the maximum propagation speed among advective/wave modes.
Regarding the numerical implementation of the model, we split the contribution of the different terms.REA integration of the advective terms on the triangle-based mesh is implemented for h, u and v at the end of the time step as described in the last section (red terms).Before that final step, the green/magenta terms are processed to provisionally update the three predictive variables, but as described in Durran (2010), a joint forward integration of these terms would make the numerical scheme unstable.Therefore, a forward-backward scheme is applied (mass variable through the green term is updated first, then followed by velocity variables through the magenta terms).To improve numerical accuracy the combined forward-backward integration is actually repeated twice, by proceeding through the midpoint of the time interval according to a second-order Runge-Kutta (RK2) cycle; alternatively, the Strong Stability-Preserving RK2 method (SSPRK2; see Durran, 2010) was also implemented: .
Factual differences in the results among the schemes, also including formulations of a third-order RK3, were too modest to recommend the necessary extra calculations.Thus, the TRAM model in the next section focuses exclusively on the RK2 cycle for the forward-backward integration of the mass-wind equations.Regarding the horizontal derivatives appearing in the green/magenta terms, they pose no problem in the triangular-based mesh.They are simply the X and Y components of the gradient vectors at T/B cells calculated using their neighbours (recall last section).This formulation corresponds to a second-order discretization of the spatial derivatives.Finally, some of the remaining terms in the SW equations (marked in grey) could be activated in the simulation.In that case, they would be just processed together with the magenta terms as part of the RK2 or SSPRK2 cycle (for the diffusion term, an equivalent second-order Laplacian operator on the triangular-based mesh can be easily formulated).
As an example of the correct performance of the time integration scheme, Figure 4 shows the results of the partial dam break problem.The same test problem can be found in Delis and Katsaounis (2005), except we rather consider meteorologically relevant scales by simulating a 2000 × 1400 km domain.Water depths on each side of the dam, before 'breaking' the structure along the off-centred breach, are 10 and 1 m.The h field shown corresponds to the state after 17 hours of the breakdown.This simulation does not include Coriolis, drag and diffusion forces and was performed with the SSPRK2 scheme using Δx = 5 km and Δt = 180 s; rigid-wall conditions were imposed along the physical boundaries of the dam and at the lateral limits of the domain.This type of tests was successfully complemented with a broad sampling of simulation parameters (including applications in one dimension; Romero et al., 2019) and with experiments of other classical shallow-water problems, like the spread of several kinds of smooth or steep water bumps (results not shown).

Non-hydrostatic fully compressible equations and numerical implementation
This section presents a first version of the TRAM model, still very limited in the scope of applicability but with the same dynamical core and numerical processing as the full model of Part 2. Specifically, we consider one of the traditional equation sets that govern the dynamics of a non-hydrostatic and fully compressible atmosphere under dry and adiabatic conditions (see Giraldo & Restelli, 2008).By combining the fundamental principles and introducing the Exner pressure, π = (P/P 0 ) R/cp , and potential temperature,  = T/π, the Euler equations in cartesian coordinates can be written as: where π' and ' are the deviations from a predefined basic state (overbar variables, function of z only) that satisfies the hydrostatic balance: The remaining predictive variables consist of the three components of wind, (u,v,w).Note we are omitting subscript d ('dry' air) for the gas constant (R) and specific heats (c p and c v ).The other parameters in the equations have the same meaning as in the previous section, but the Coriolis force is completed with its other components through f ̂= 2Ω t cosLat.The reader will note the same colour convention that was used for the shallow-water model.This colour coding has the same meaning in connection with the applied integration scheme.Recall, green/magenta terms are integrated over time step Δt in a forward-backward sequence and in two legs (RK2 cycle), using centred differences for the horizontal and vertical derivatives.In contrast, we now realize and take advantage of the broad range of waves and motions that are possible in the above set of unfiltered atmospheric equations, ranging from slow advective flows to fast acoustic waves (>300 m⋅s −1 ).Specifically, a time-splitting strategy has been implemented (e.g., Wicker & Skamarock, 1998).Instead of updating advection terms (vertical advection first, followed by horizontal advection) with the REA method every short time step Δt, they are updated at longer intervals, typically every 6-10 time steps (we will refer to this multiplier as Nstep).This is an interesting solution for the red terms, since the REA method is computationally intensive.
CFL conditions for numerical stability, considering the fastest acoustic modes, allows a maximum stable time step somewhat above Δt≈2Δx(Δz), with time in seconds and grid lengths in kilometres.The shortest values of Δx and Δz found across the domain should be used in this expression, but since the numerical experiments of this section use homogeneous grids with Δx = Δz, this type of consideration does not enter into play for the moment.In addition, by mimicking other authors, most of the validation tests were done with a two-dimensional version of the model operating in the (x,z) plane.This version simply neglects  y derivatives in the above equations, the REA algorithm and zonal/vertical derivatives are fully equivalent in both coordinates and, in the case of no rotation (i.e., Coriolis parameters set to zero), the v equation is also omitted.The corresponding rule of thumb for the maximum time step relaxes to Δt ≈ 3Δx(Δz).
Finally, diffusion terms with constant  coefficient have been added for consistency in the equations.These grey-coloured terms do not aim to represent the viscous processes in flows governed by the full set of Navier-Stokes equations, since atmospherically relevant circulations have high Reynolds numbers and are nearly inviscid.These terms would express the possible need of some type of filters for scale-selective dissipation, acting on the shortest and least reliable wavelengths of the numerical solution in order to prevent nonlinear instability (see Durran, 2010 for relevant discussions on this issue).However, as noted at the end of Section 2.1, this kind of explicit filters becomes unnecessary for the TRAM model (even in the absence of any turbulence parametrization; this will be treated in Part 2) thanks to the small-scale damping implicitly associated with the REA application.
We subjected our model to several benchmark tests found in the literature, displaying the same fields and output times for an easier validation of results.Some of these tests involve the simulation of thermal bubbles embedded in calm and neutrally stable environments, with no earth rotation.In two dimensions these tests can be designed at high resolution with low computational cost.The first example (Bryan & Fritsch, 2002) considers the evolution of a circular (2 km radius) axisymmetric warm anomaly located near the ground, with ' defined by a cos 2 profile and reaching +2 K at its centre.Figure 5a depicts the shape of the buoyant anomaly after 17 min, as in the reference work, and confirms the correct development of two rotors on the sides of the thermal and a large temperature gradient in its upper part.This simulation was performed with Δx = Δz = 50 m, Δt = 0.125 s and Nstep = 10.The subsequent effects of Kelvin-Helmholtz instability along the periphery of the thermal during its ascent lead to an irreversible loss of softness and the production of complex ' patterns (e.g., Robert, 1993).Figure 5b displays one of such attractive patterns, corresponding to the same type of simulation as Figure 5a but performed at double resolution and shown after 33 min.A second validation test is shown in Figure 6, which contains the same four stages of the exercise proposed by Robert (1993), regarding the interaction of a rising large warm bubble and a descending small cold bubble introduced again in a calm and neutrally stable environment.This experiment was run using Δx = Δz = 2.5 m, Δt = 0.0625 s and Nstep = 10.Times shown are 0, 4, 7 and 10 min, as in Figure 9 of the reference paper, and confirm the ability of TRAM to perfectly replicate the details of the interacting bubbles.Another well-known test for NHFC models is the density current simulation (Straka et al., 1993), initialized with a strong cold anomaly (dimensions 8 × 4 km, ' = −15 K at its core) embedded in the same kind of environment as the previous experiments.Our Figure 7   Finally, we successfully reproduced the inertia-gravity wave experiment of Giraldo and Restelli (2008) in Figure 8.This experiment uses Δx = Δz = 125 m, Δt = 0.3125 s and Nstep = 10 over a domain 300 × 10 km in size.The waves are forced by introducing a very small thermal perturbation (' = +0.01K at the maximum point in Figure 8a) in a non-rotating environment characterized by uniform wind and stability.Specifically, initial background wind is from the west, u = 20 m⋅s −1 , and vertical stratification is given by N = 0.01 s −1 (N is the Brunt-Väisälä frequency, defined as N 2 = gd[ln]/dz).After 50 min of simulation (Figure 8b) the TRAM model replicates in detail the results of the reference study, except near the lateral boundaries owing to the fact we have still used rigid-wall conditions in all the above experiments.Proper treatment of boundary conditions in our model will be discussed in the next section.

Inclusion of orography and lateral boundary conditions
Mesoscale applications over complex terrain demand effective inclusion of orographic forcing in the numerical model.Terrain-following vertical coordinates such as sigma and hybrid coordinates facilitate, in principle, this purpose, but these schemes are not free of serious problems over steep slopes owing to the decomposition of the pressure gradient force into two terms in the formulation of the equations (see Schär et al., 2002).Since we chose the z-coordinate for the TRAM model, a proper representation of the stepwise terrain and particularly of the slope-generated forcing will be necessary to guarantee realistic results with this terrain-intersecting coordinate.
Consider the database of terrain elevation for the region of interest (black silhouette in the background of Figure 9).This digital information will normally be available at much finer detail than the model grid resolution Δx ⨯ Δz (the grid cell size in the two-dimensional scheme of Figure 9 is indicated as the light-orange square near the top-left corner).First, the average terrain height at the grid resolution is calculated (dark-orange lines) and, based on this, the terrain mask is built: grid cells intersecting these dark-orange elevations and those located directly underneath will conform the mask.That is, the light-orange cells along the bottom of Figure 9 are considered ground cells, while those cells located above define the atmospheric domain.Only over oceanic regions these atmospheric columns would extend down to the lowest cell.
What is crucial in our treatment of orographic forcing is the terrain slope that is assigned to the ground-atmosphere transition cell at the bottom of each atmospheric column.Instead of assigning the 'grid' slope (a slope calculated using the stepwise height profile drawn by the light-orange cells) we will assign the 'true' slope (a slope calculated using the more gradual height profile provided by actual terrain elevations, i.e., the dark-orange values).Figure 9 exemplifies this important distinction between both slope formulations and their effects on the flow equations.While the grid slope would lead, from left to right, to an unnatural shift from null mechanical forcing to excessive forcing, the true slope would induce a more gradual and realistic change of the orographic forcing along the domain.With our approach we will avoid anomalous flow simulations such as those observed by Gallus and Klemp (2000) when using step-terrain coordinates, problems that are especially evident with poor vertical resolution.Note also that the procedure illustrated in Figure 9 for the two-dimensional case is easily generalized for a three-dimensional domain.This would end with the construction of the terrain mask on the triangular mesh of Figure 1 and the calculation of the slope vector (slope) at each triangle using the true elevations at its three neighbours.The topographic bottom boundary condition is written for the first ground cell (gnd) using information of the lowest atmospheric cell (atm) as a terrain-following flow condition: which is complemented by a zero-order extrapolation of other prognostic variables: The performance of the TRAM model under this treatment of orography was tested for the density current experiment of the last section by introducing in the domain analytically defined mountains of significant amplitude (results not shown), but before moving to more traditional tests involving mountain waves, a gravity wave-absorbing layer had to be formulated in the model.In effect, the finite vertical size of the domain and the presence of an artificial 'rigid-wall' boundary at the top would induce a significant downward reflection of any topographically generated energy, severely contaminating the numerical results.Following the ideas proposed by Klemp et al. (2008) for non-hydrostatic model equations that are solved using split-explicit time integration techniques, as ours, we simply add a Rayleigh damping term to the vertical velocity equation: where the maximum value of the damping coefficient, τ max , is 0.1 s −1 and this damping is applied over a depth Z D = 10 km from the model top, located at Z T .These are standard values that might change for certain experiments.
In addition, if the vertical dimension of the domain is less than 20 km, Rayleigh damping is forced to remain active only above 10 km altitude.As we will see, this simple method is an effective way of absorbing upper propagating gravity wave energy in the simulations and avoiding spurious results in the tropospheric layers.
On the other hand, we arm the TRAM system with a more versatile scheme for the lateral boundary conditions that has been proven to display good behaviour in many other models (e.g., Grell et al., 1995;Skamarock et al., 2008), avoiding spurious reflections and a smooth transition of the interior solution towards external data.These external data (referred to below as ϕ LS ) would usually consist of a fixed vertical profile of the meteorological variables coincident with the initial sounding used in idealized experiments or, more generally (Part 2), of spatially and temporally variable large-scale fields provided by a reanalysis or a global model.The relaxation scheme (Davies & Turner, 1977;Marbaix et al., 2003) consists of 'nudging' the model predicted variables ϕ mod towards ϕ LS by the action of Newtonian and diffusion terms: where this adjustment is applied after every time step at the four outermost grid cells of the domain, with a weight decreasing linearly with distance, from a value of 1 at the lateral boundary, to 0 at the first strictly interior cell.
In addition, the coefficients F and G are taken as 1/10t and 1/50t for idealized experiments such as the following examples, and five times greater when forcing with large-scale analyses or forecasts.
Our model successfully passed several kinds of mountain wave tests.In this section we include simulations of linear mountain waves as in Klemp et al. (2008) and replication of the mountain wave test of Schär et al. (2002).The linear mountain waves are forced by combining an initial environment possessing uniform westerly wind and stability (u = 10 m s −1 and N = 0.01 s −1 ) with a 10 m high bell-shaped mountain in three different cases regarding the mountain half width: 10, 50 and 2 km (recall, u = 10 m s −1 and N = 0.01 s −1 ) to the presence of a complex mountain of 250 m height.This mountain is defined by an analytical profile composed of a bell-shaped envelope combined with shorter scale cos 2 -shaped components.Consequently, the solution is very interesting as it contains a rich mixture of a larger-scale hydrostatic wave, with deep propagation in the vertical, and smaller-scale non-hydrostatic waves, rapidly decaying with height (see Schär et al., 2002).The results with TRAM in Figure 11 were obtained with Δx = Δz = 250 m, Δt = 0.75 s and Nstep = 10 using a domain of 200 × 20 km size.The standard value of τ max was used.This output and additional versions with changed resolution (not shown) can be directly compared with other model simulations and with the analytical solution for the stationary wave (see fig. 13 in Schär et al., 2002).The clear conclusion of this comparison is that TRAM is perfectly suitability for capturing any kind of mountain-induced disturbance, totally free of numerical artefacts associated with the use of the height coordinate and stepwise terrain.The reader will note in Figure 11 the apparent lack of a mountainous terrain at the bottom of the simulation domain; this is explained by our definition criteria for the terrain mask in a case in which Δz coincides with the mountain height; nevertheless, the proper treatment of mountain slopes in our method produces a numerical solution virtually identical to the analytical wave.

Vertical stretching and semi-implicit scheme
At this point the TRAM model still uses a uniform resolution Δz in the vertical.This is a disadvantage considering that a proper parametrization of turbulent mixing in the PBL (Part 2) will demand enhanced resolution at low levels.Finer resolution close to the surface would also lead to a better representation of highly variable terrains.In response, while maintaining the number of computational levels within reasonable limits, we improve the model by allowing a heterogeneous resolution in the vertical: levels are brought closer together in the lower troposphere, where more resolution is needed, and gradually stretched towards the top of the domain, where poorer resolution is less critical.The degree of stretching is flexible and user-defined by the simple choice of two parameters: the 'mean' resolution Δz m and the stretching parameter stretch.Then, the minimum thickness or highest resolution at sea level is given by Δz min = Δz m /stretch, and the maximum thickness or poorest resolution at the top of the domain is given by Δz max = Δz m + (Δz m − Δz min ).Intermediate resolutions for the other vertical nodes are provided by a cos profile connecting Δz min and Δz max values, although any other transition function could be easily formulated.The reader will easily note that the special case stretch = 1 corresponds to the previous scheme with uniform resolution, Δz = Δz m .
Several actions had to be undertaken to adapt the numerical scheme to the vertically stretched height coordinate.First, all vertical derivatives had to be reformulated in view of the variable resolution and this was done keeping the centred, second-order character in the discretized expressions.But the major change in the model deals with the practical implications of now using very short Δz in some regions of the domain.This ultrafine grid resolution would severely restrict the maximum time step compatible with numerical stability according to the CFL condition.
To circumvent this problem, acoustic vertical modes linked to the  z w and  z ' terms in the first and last equations of Section 2.3 respectively, are stabilized through the semi-implicit formulation of these derivatives in the finite difference form of the prognostic equations.This is a common practice in numerical codes of elastic models (e.g., Skamarock et al., 2008) and basically means that the CFL criterion behind the maximum allowable time step Δt will now be determined exclusively in terms of the horizontal grid length Δx.That is, Δt ≈ 2x for the three-dimensional model and Δt ≈ 3x for the two-dimensional version, without taking into account the vertical resolution.Of course, the stabilization of fast vertical modes comes at the price of degrading their accuracy in the numerical representation, but since these modes are of little meteorological significance, the loss of accuracy is basically inconsequential for the spatial and temporal scales of interest (Durran, 2010).Specifically, the integration of the green/magenta terms in the above equations now proceeds as follows i : where the problematic vertical derivates have been split into explicit (i.e., at step n) and implicit (at step n + 1) components using an off-centred scheme (α = 0.3, β = 0.7, as in Skamarock et al., 2008).The resulting expression for ' n+1 is substituted in the w n+1 development and -following the red arrow -an implicit expression for w n+1 is obtained that depends on their own first and second derivatives, apart from explicit information (A, B, C, D or their derivatives).Since the discrete expressions for  z w and  zz w are formulated with a centred scheme, we end with a tridiagonal problem for the vertical velocity at the klevels (again, a, b, c and f coefficients are fully explicit, i.e., they only depend on fields at step n).The tridiagonal problem is numerically solved for w n+1 with any of the standard techniques (e.g., the code available in appendix A.2.1 of Durran, 2010) and applying in this process the topographic and zero-velocity boundary conditions at the domain limits.Finally, by resorting back to the ' n+1 expression, the perturbation pressure field is also updated.With reference to the previous scheme, recall that the time integration over the short time step Δt concludes with the update of u n+1 and v n+1 using the forward-backward scheme.
Additional optimizations were introduced at this stage of model re-development.In short: (i) Flexibilization of the REA scheme for the vertical advection, in the sense of allowing cells to affect computational layers beyond their immediate neighbours; this is necessary to accommodate the action of strong updrafts or downdrafts in layers of high vertical resolution, once Δt is no longer ruled by the vertical grid length.(ii) Updating the Coriolis horizontal components and the entire ' equation, considered slow processes, solely in the Nstep cycle of the model, as the advection terms.(iii) Finally, in the rare instance that the grey-coloured terms in the equations of Section 2.3 need to be activated (they also are used in the Nstep cycle), application of an implicit numerical scheme for the vertical diffusion, another consequence of using a comparatively long Δt in this optimized version of the TRAM model.
The correct implementation of all the above improvements was verified by repeating, first, the same tests of Sections 2.3 and 2.4 under a stretch parameter equal to one.Virtually indistinguishable results were obtained.In the second place, as an interesting sensitivity analysis we repeated the Schär mountain wave test of Figure 11 for several degrees of stretching (stretch = 1, 5, 10, 20, 30) but keeping all other simulation parameters unaltered (we highlight the use of Δz m = 250 m with 81 vertical levels, Δt = 0.75 s and Nstep = 10).The stationary mountain wave remained nearly insensitive to the magnitude of the stretching, but when the vertical resolution was degraded (Δz m = 500 m with 41 vertical levels, and Δz m = 1000 m with 21 vertical levels) the wave lost appreciable entity unless some degree of stretching was activated in the model (stretch > 2 and stretch > 4 for these coarse-resolution experiments respectively; figures not shown).
Crucial illustrative examples of the definitive progress achieved in the design of the dynamical core of the model deal with the simulation of intense mountain waves using realistic soundings.On the one hand, the well-known 11 January 1972 Boulder windstorm event (Klemp & Lilly, 1978) was successfully simulated with the same initial data and domain characteristics as in Doyle et al. (2000).In fact, our simulation is fully compatible with the results of any of the 11 models intercompared in the referred study (TRAM fields not shown).On the other hand, we replicated the experiment of Doyle et al. (2011) and simulated the mountain wave excited by Sierra Nevada range (CA, USA) under the impingement of a flow that was initialized using a real upstream sounding, taken during the T-REX observational campaign (March 2006) and kindly provided by the author.A realistic topographic profile was also used for the simulation (Figure 12), which was run under Δx = 500 m, Δz m = 100 m, stretch = 5, Δt = 1.5 s and Nstep = 6 over a domain of 500 × 25 km.Non-linearities and transient features are highly influential in this case and obviously no stationary state is reached, but persistent features in the simulation are the downslope windstorm conditions and hydraulic jump found at low levels, along the eastern slope of the ridge, and the profound wave breaking occurring aloft.In Figure 12 these structures have been represented after four hours of simulation exactly in the same way as in the model intercomparison composite contained in fig. 5 of Doyle et al. (2011), which again analyses 11 different models.We can easily conclude that the TRAM model performs at least as well as state-of-the-art non-hydrostatic and NHFC modelling systems.

Additional validation tests
We close this series of dry-adiabatic simulations with the dynamical core of TRAM by considering the simulations of the reference study, highlighting the two vortices formed to the lee of the topographic obstacle, especially for the highest mountain.As expected, if the obstacle is elongated along the N-S axis, normal to the background wind, the splitting of the flow on the upstream side of the mountain is accentuated and the lee vortices become more robust (a nice example after tripling the meridional half width is included in Figure 13a).In the second place, asymmetries in the problem were introduced, either by considering non-circular mountains or, exactly in the same way as in Schär and Durran (1997), by imposing a small asymmetry in the initial temperature field.As in their simulations, these perturbations produce a relatively rapid transition to the asymmetric vortex shedding regime (see the TRAM example after six hours in Figure 13b).Oscillating Von Kármán vortex streets are developed in the wake of the obstacle, with the vortices persisting for many hours as they drift downstream (our simulations covered up to 48 hours).These striking simulated features are in good agreement with the real effects of mountainous islands, as occasionally revealed by satellite cloud imagery (e.g., Etling, 1989).

Physical parametrizations and a new form of equations
Many formulations of the subgrid-scale processes, including radiative and moisture-related forcings, have been developed during the last decades for inclusion in short-range to global numerical models (see Stensrud, 2009 for a comprehensive review).They come with different degrees of realism or hierarchy levels, which is tantamount to saying computational cost.It could be argued that any physical parametrization, eventually expressed as a piece of code, represents in itself a scientific topic that grows out of a fruitful combination of theory, observations and experiments.The development of new parametrizations for the TRAM model would be very arduous and is clearly beyond the scope of this work.We followed a pragmatic solution, consisting of adapting a set of parametrizations of intermediate complexity directly from the MM5 model (described in Dudhia, 1993 andGrell et al., 1995).This task required the arrangement of appropriate interfaces to connect the physical modules with our dynamical core, as well as significant reformulations of the Fortran source programs, considering that the prognostic variables are unique for each model and that we use the height coordinate (e.g., the vertical index of the field arrays increases upwards) while MM5 uses the sigma coordinate (the vertical index increases downwards).The adopted physical parametrizations form the same combination of schemes that was prioritized by Romero et al. (2014) on the basis of MM5 mesoscale simulations of a mediterranean convective storm.Their key characteristics and supporting references are summarized below: • Microphysics (Reisner graupel or Reisner 2 scheme; Reisner et al., 1998) -Extends the more primitive mixed-phase (Reisner 1) scheme by adding equations predicting graupel and ice number concentration and all relevant microphysical processes to the original set of explicit equations for the liquid phase (cloud and rain water fields) and ice phase (cloud ice and snow).The scheme is suitable for cloud-resolving models.
• Cumulus (Kain-Fritsch 2 scheme; Kain & Fritsch, 1993;Kain, 2004) -A newer version of the Kain-Fritsch scheme that includes shallow convection.This type of scheme works similarly to the Fritsch-Chappell strategy that relaxes to a profile due to properties of updraft, downdraft and subsidence region (see Fritsch & Chappell, 1980, for details) but using a sophisticated cloud-mixing scheme to determine entrainment/detrainment, and removing all available buoyant energy in the relaxation time.This scheme predicts both updraft and downdraft properties and also detrains cloud and precipitation.Shear effects on precipitation efficiency are also considered.
• Radiation (cloud radiation scheme; Benjamin, 1983) -This parametrization scheme is sophisticated enough to account for long-wave and short-wave interactions with explicit cloud and clear air.In addition to atmospheric temperature tendencies, the scheme provides surface radiation fluxes.It may be moderately expensive (this is typical of the radiation calculation component for any model) but it requires little memory.
• Planetary Boundary Layer (MRF PBL or Hong-Pan PBL scheme; Hong & Pan, 1996;Troen & Mahrt, 1986) -This scheme is suitable for high-resolution PBL (e.g., five layers in the lowest km and the surface layer less than 100 m thick).It is an efficient scheme based on a Troen-Mahrt representation of the countergradient term and K profile in the well-mixed PBL (see Hong & Pan, 1996 for details).Vertical diffusion uses an implicit scheme to allow longer time steps.
• Surface (five-layer soil model scheme; Dudhia, 1996) -Temperature on land is predicted in 1-, 2-, 4-, 8-and 16-cm layers (approximately), with fixed substrate below, using a vertical-diffusion equation.Thermal inertia is formulated as in the force/restore scheme (Blackadar, 1979), but by vertically resolving diurnal temperature variation, the scheme allows for a more rapid response of surface temperature.Additionally, moisture availability in the soil varies with time, particularly in response to rainfall and evaporation rates.Seasonal roughness, moisture content, and radiative and thermal properties are prescribed for each location according to the USGS classification of up to 24 different land-use types.It should be noted that on water bodies the surface temperature is kept fixed during the simulation.
Endowing the model with the physics package implies the addition of new prognostic fields, namely the mixing ratios of six water species (water vapour Q v , cloud water Q c , rainwater Q r , cloud ice Q i , snow Q s and graupel Q g ) plus ice number concentration (NC); in the notation we group these seven variables as Q χ .In addition, the parametrizations interact with each other and with the dynamics through their output tendencies, or forcings, for temperature (F T ), horizontal wind components (F u and F v ), water vapour (F Qv ) and the remaining water variables (F Qχ ).Additional fields delivered by the physical schemes include, among others: non-convective and convective rainfall; near-surface wind, temperature and moisture; temperature and moisture at ground level; and the several components of the surface energy budget.
The complete form of the model equations, now accounting for the new prognostic variables and forcing terms, i.e., applicable to the real moist atmosphere fully interacting with the earth surface and external radiation, are given below.The main modifications with respect to the dry-adiabatic version of Part 1 are highlighted in blue.In this form of the equations all effects of moisture on pressure (first equation) and thermodynamics (second equation) have been retained, as in the CM1 numerical model (Bryan, 2021).In fact, since R m and the specific heats c pm and c vm will not differ too much from their dry counterparts (R d , c p and c v ), the moist coefficients coupled to the forcing terms in these two equations contain only small corrections with respect to the dry atmosphere or would even tend to zero (divergence term and last term in the ' equation).Some NHFC models legitimately avoid these second-order effects on their prognostic equations (such as ARPS, MM5, WRF and Klemp-Wilhelmson model; see Bryan, 2021 for details) but we found a slightly positive influence on the validation tests and thus prefer to retain all of them.A further necessary modification for correctly modelling the moist atmosphere is the use of density potential temperature,  ρ , instead of potential temperature in the momentum equations.The reader will also note the inclusion in the w equation of the drag force linked to the presence of liquid (Q liq ) or solid (Q ice ) condensate.
Finally, it should be noted that the above equations are presented with no approximations and for the most general applications, intended to simulate regional or large domains over the spheric earth (of radius a).Specifically, all Coriolis and curvature terms have been retained (note f = 2Ω t sinLat and f ̂= 2Ω t cosLat) and a map projection (currently the Lambert conformal conic, quite appropriate for mid-latitudes; Haltiner & Williams, 1980) has been applied (note m is the map scale factor and α is the angle between the local meridian and x = constant).We refer to these kinds of problems as 'Synoptic-real case' (Section 3.4).On the other hand, one might be interested in more academic, small-scale simulations performed over a generic computational box, without regard to the real earth.We will refer to these problems as 'Mesoscale-idealized' (Section 3.3).The corresponding equations are easily expressed by considering a → ∞ in the above set and replacing m and α by the constant values 1 and 0 respectively.

Time-marching algorithm
Before proceeding with the validation tests, the sequence of steps that make up the time-marching procedure of the full model should be presented.Indeed, it is not obvious in which order the physical tendencies should be calculated and where and how these forcings be incorporated in the above prognostic equations.First, we note that physical forcings belong to the 'slow' part of the model, just as advection, and therefore they can be called with low frequency to save computational time, that is, in the Nstep cycle.Apart from that, the results exhibit some sensitivity to the specific order of the physical parametrizations in the global algorithm and to their mode of coupling with the dynamical core.Some simulations might even destabilize when repeatedly exposed to some sort of physical shock if all parametrizations are acting together.Numerous sensitivity tests for the kind of problems analysed in next sections led to the conclusion that the optimal time-marching algorithm for implementing TRAM consists of the following steps: From top to bottom in the FLOWCHART, and clearly distinguishing fast processes (accounted for in the short cycle of time step dt; blue) from slow process (long cycle of dtNstep; orange), the numerical integration obviously starts with the ingestion of the initial conditions for all fields (I.C.).In the short cycle, mass and wind equations are integrated in combination exactly as in Part 1, that is, in a forward-backward sequence using the RK2 method and the semi-implicit scheme for the vertical coupling of ' and w.In this part of the integration, not only the genuinely slow processes (advection and physical forcings F) are skipped, but also some additional terms.While w incorporates all other terms, and ' both the basic state influence and the divergence term, the u/v equations omit all Coriolis and curvature terms, that is, they only incorporate the pressure gradient force in the RK2-cycle.The Coriolis and curvature terms are added later, at the end of the short cycle (last line of the blue diagram).Just before, the ' equation is integrated based on the basic and divergence terms exclusively.The very minor influences of F T and F Qv on ' are also incorporated at this point, using the last available values of these forcings.At the end of the short cycle the boundary conditions for all fields (B.C.) are conveniently updated and a new integration over dt starts (n + 1 → n).
Each time the user-defined Nstep multiplier of dt is reached, the model enters the long cycle, or orange diagram of the displayed scheme, and slow processes over the long time step dtNstep are integrated.The order of calculations is again indicated by the grey-shaded arrow, from top to bottom.First, the PBL/Surface schemes are called and their forcings (i.e., mixing) are accumulated over u, v, Q χ and '.Surface variables such as ground temperature, T g , etc, are also predicted at this step.Other artificially introduced sources of diffusion, either horizontally (explicitly calculated) or vertically (an implicit scheme is used owing to the small grid lengths near ground), or both, would be activated next, but we emphasize again these are needless components of the model.Before calling the other physical packages (Radiation, ii Cumulus and Microphysics, in this order; see orange diagram) the advective tendencies on all prognostic fields are applied.These tendencies emerge from the REA algorithm and are calculated/applied in two steps, first vertically (REA-V) and then horizontally (REA-H), as explained in Part 1. Focusing on radiation (next step after Advection), this scheme will not only update the surface radiation fluxes but will also provide a first contribution to the atmospheric temperature tendency (stored in F T ); this tendency will be completed (+F T ) by the cumulus and microphysics packages (see scheme).Likewise, cumulus parametrization provides a first contribution to water species tendencies (F Qχ ) which are later completed by the microphysics package (i.e., +F Qχ ).All these tendencies are finally applied on the corresponding fields, Q χ and ', to conclude the Nstep cycle.Since thermodynamic coefficients strictly depend on moisture variables (cf. the last section) their values are also updated at the end of the long cycle.The reader will also note that, once this cycle is finished, the model is screened for possible output of primitive and diagnostic meteorological fields, according to the frequency requested by the user.Among the long list of output variables, we highlight the convective (R c ) and non-convective (R nc ) components of the total rainfall field, yielded respectively by the cumulus and microphysics parametrizations.

Mesoscale-idealized simulations
We start by exploring the capability of TRAM for simulating thermally driven circulations influenced by local orography.An ideal laboratory for these kinds of mesoscale circulations is the island of Mallorca, located at the centre of the western Mediterranean.With an approximate size of 80 × 80 km, the island develops a well-defined sea/land breeze regime almost daily during the prevalently anticyclonic warm season.The diurnal sea breezes are particularly regular and intense along the three major bays and nearby lowlands and, at the mature state, converge from all directions at the centre of the island.This circulation was for the first time systematically described by Jansà and Jaume (1946) and has been simulated in detail by numerical models, both hydrostatic (Ramis & Romero, 1995) and non-hydrostatic (Cuxart et al., 2014).Ramis and Romero (1995) also showed the reinforcement role exerted by the typical dryness of the soil and, in some sectors, by the overlap of upslope winds diurnally forced by the Tramuntana northern mountain range.The diurnal cycle was simulated under ideal circumstances (i.e., calm synoptic wind and horizontally homogeneous fields initially) starting from the Palma de Mallorca sounding of 0000 UTC 30 August 2004.Climatological values for sea surface temperature and soil moisture availability were assumed: 25 • C and 10% respectively.We used Δx = 1.5 km (i.e., square-based resolution is 1 km), Δz m = 400 m, stretch = 20, Δt = 3 s and Nstep = 10 in a computational domain of 180 × 155 × 16 km; note the high stretch parameter (meaning a refinement up to Δz = 20 m at sea level) in this type of simulations, since PBL physics is better represented using rich vertical resolution.The model reproduces with superb detail the characteristics of the full diurnal cycle, comprising the onset of a light land breeze during nocturnal hours and early morning, the subsequent development of a neat sea breeze during the day which ultimately affects the whole island, and the transition in the last hours of the 30-hour-long simulation towards new land breezes along the coasts.We include as an example the simulated wind field at mature state (Figure 14a).As in previous studies, the pre-eminence of the flows entering through the main bays, the concomitant inland wind convergences (even shaping a cyclonic vortex at the centre of the island), and the role of the Tramuntana range (both hampering the breeze from the north and deflecting/enhancing the breezes in the west-central plains) are all very clear.We also highlight the role of the breeze-induced convergence zones as a triggering mechanism for afternoon convection under propitious environments (Jansà & Jaume, 1946).These favourable environments are most commonly found in late August and September.In fact, the considered sounding featured a convective available potential energy (CAPE) value exceeding 1200 J kg −1 and a convective inhibition of about −100 J kg −1 ; in those circumstances and following the progressive warming and mixing of the PBL over land, the mature breeze serves as an ideal mesoscale mechanism to release the convective instability by low-level convergence.On 30 August 2004 heavy rainfall of this kind occurred in the interior of Mallorca, even leading to some local floods.Our simulation (Figure 14b) succeeds in capturing the genesis of short-lived afternoon downpours along the convergence areas and mountain slopes, confirming a proper implementation and interconnection of the physical parametrization schemes in TRAM.Additionally, the role of a non-zero synoptic wind was investigated by prescribing geostrophically adjusted basic currents of different speeds and directions in a new set of idealized initial conditions, and also by ingesting the actual wind profile of the Palma de Mallorca sounding.As expected, for slight to moderate values of the basic synoptic wind, its combination with the thermally driven forcing results in a downstream shift of the breezes, convergence zones and precipitation structures (results not shown).
Further, we have simulated with TRAM the genesis and evolution of two remarkable types of long-lived, organized convective storms: squall lines and supercells, problems recurring very frequently in severe weather studies embracing NHFC models.We resort again to simplified contexts for the simulations, consisting of idealized vertical soundings of temperature, humidity and wind for defining horizontally homogeneous initial states and fixed boundary conditions.As in the reference studies, we focus on the leading mechanisms that control the dynamics of these small-scale storms and, therefore, a non-rotating domain is used (i.e., Coriolis force is zero) and we omit radiation, PBL physics and cumulus parametrizations.That is, the cloud microphysics scheme is the only physical package turned on in these fine-grid simulations and moist convection is expected to be fully explicitly resolved.In addition, flat-bottom domains are used in both cases, and given the lack of triggering mechanisms in the initial fields, the convective storms are promoted by imposing some kind of localized thermal perturbation at the initial time.
The squall lines are simulated in a computational domain of 600 × 260 × 16 km using Δx = 1.5 km (square-based resolution 1 km), Δz m = 200 m, stretch = 10, Δt = 3 s and Nstep = 5.The classical WK82 thermodynamic sounding (Weisman & Klemp, 1982) is employed to initialize the vertical profiles of temperature and moisture.This sounding yields a surface-based CAPE of 2200 J kg −1 .Apart from instability, squall line development and maintenance demand favourable vertical wind shear; the balance between the storm-induced cold pool, the rear-inflow current and the environmental shear allows achieving a steady state and a great longevity (Weisman, 1993).Specifically, three different TRAM experiments that ingest the same three idealized wind profiles defined in Weisman et al. (1997) are launched.These profiles vary from weak to moderate and to deep shear, affecting the intensity and speed of the squall line and its degree of upshear orientation over time.In all three cases the linear convection is kicked off by a meridionally extended cold pool introduced along the western fringe of the domain.This ' perturbation linearly decreases from −8 K at surface to zero at a height of 2.5 km.Convective updrafts are readily forced along the leading edge of the cold pool, then environmental shear favours their continuous upright regeneration, and, in less than two hours, the eastward-moving storms are fully organized as squall lines.
Figure 15 displays a W-E central cross-section of the simulated squall line for the moderate shear case, five hours after its initialization.This environment prescribes a westerly wind linearly increasing from zero at surface to 17.5 m s −1 at a height of 2.5 km, and then kept uniform above this level (Weisman et al., 1997).The mature convective line moves precisely at that speed, approximately, so we subtract the value of 17.5 m s −1 from the model simulated u-wind in order to obtain the storm-relative circulation of Figure 15a.It is clear that TRAM reproduces with a wealth of details the known features of vivid squall lines and fulfils the conceptual model of these convective systems (Houze et al., 1989).In particular: (i) the morphologies and relative positions of cloud top, cloud base (including signs of a shelf cloud) and radar echo boundary are well reproduced; (ii) the model reflects the coexistence, from east to west, of incipient updrafts, vigorous convective towers, and old cells; (iii) both the ascending front-to-rear flow and the descending rear inflow become very distinctive features in the simulated storm-relative circulation; and (iv) the thermal influence of the massive latent heat release in the core of the storm and the development of a strong low-level cold pool with a gust front, in phase with the system, are also neatly reproduced by TRAM.The supercells are simulated in a computational domain of 300 × 260 × 20 km using Δx = 0.75 km (squarebased resolution 500 m), Δz m = 400 m, stretch = 20, Δt = 1.5 s and Nstep = 5.Our initial configuration essentially replicates the conditions of the idealized WK82 supercell simulated by Potvin and Flora (2015).That is, temperature and moisture are initialized again using the severe convective sounding of Weisman and Klemp (1982) and, as wind profile, the quarter-circle hodograph is used.As in the reference study, convection is triggered by a thermal bubble with ' = 3 K at its centre (located 1.5 km above the ground) and with radii of 10 and 1.5 km in the horizontal and vertical respectively.This initial thermal bubble is placed near the western boundary and meridionally centred, since the simulated storms will move predominantly eastwards.
Hodographs turning clockwise with a height as the one used (the 0-3 km storm-relative helicity, SRH, exceeds 200 m 2 ⋅s −2 ) should induce the intensification (weakening) of the right-(left)-moving supercells following the splitting of the original cell (Klemp, 1987).Figure 16 displays the TRAM results for some key fields after 4.5 hours of evolution.In effect, the right-mover storm (in the lower portion of Figure 16a) evolves into a supercell storm with a clear hook echo structure, while the left-mover has led to more disorganized convection, composed of several cells forced along spreading outflow boundaries (upper half of Figure 16a).The surface wind associated with the supercell (Figure 16b) reflects the severe character of these storms; in addition, the simulated pattern is fully compatible with the conceptual scheme of the mature circulation (Lemon & Doswell III, 1979), which highlights the mesocyclone around the updraft and the gust fronts associated with the forward-flank (FFD) and rear-flank (RFD) downdrafts.Finally, the role of the storms for stabilizing the environment (leaving cold air at low levels in their wake) and as efficient rainfall producers, is visibly displayed in Figure 16c,d respectively.As expected, when the initial hodograph is forced to adopt a straight profile (by setting to zero the v-component at all heights), the splitting of the updraft leads to the formation of right-and left-moving supercells with mirror-image symmetry (Markowski & Richardson, 2010) (simulation not shown).

Synoptic-real-case simulations
Virtually all kinds of atmospheric evolutions are of interest for the validation of a numerical model.Countless real cases, each with their own specificities, could be listed for this purpose.For the sake of brevity, we focus the validation task of TRAM on examining its performance for a few prototypical high-impact weather situations that affected the Mediterranean region and nearby countries.Specifically, we simulate a historical extratropical cyclone (the 'Superstorm' event), an explosive cyclogenesis event of sub-synoptic size (storm Hugo), a Mediterranean hurricane (medicane Zorbas), a heavy precipitation episode linked to a mid-tropospheric cut-off low, or 'Dana' in Spanish, and the extreme winds and abundant rainfall produced by the extraordinary storm Gloria.The quoted names of these storms were at the time adopted by the National Weather Services of the affected countries.All TRAM simulations use a domain centred in the Balearic Islands (western Mediterranean) but with the size and grid spacing depending on the event (see below).The Superstorm of 10-12 November 2001 (Romero, 2008) was the result of a deep baroclinic development in the western Mediterranean basin.The devastation produced by wind and rain in coastal regions like the Balearic Islands, and the loss of human lives in Algeria, have no parallel in recent decades.The simulations are initialized at 0000 UTC 9 November and forced at the lateral boundaries using the National Centers for Environmental Prediction meteorological grid analyses, available at 12-hour intervals.The examples shown in Figure 17a,b correspond to two independent simulations performed at coarse (Δx = 50 km) and high (Δx = 9 km) resolution respectively.The short time steps are set in proportion (Δt = 75 s and 15 s) and Nstep = 6 in both cases.Both simulations use a vertical domain of 16 km, with Δz m = 200 m and stretch = 1 (thus, 81 computational levels are included).The synoptic-scale domain correctly captures the intense cyclonic development that followed the mere initial presence of a large-scale trough along the western flank of Europe.At 0000 UTC 11 November, the peak moment of the storm, TRAM effectively places the cyclone centre to the southeast of the Balearic Islands, producing extreme sea-level pressure gradients over broad areas of the Mediterranean (Figure 17a).These features, as well as the characteristic asymmetric cloud structure wrapping around the low, are in good agreement with the verification analyses and satellite images (not shown).This is interpreted as a good unfolding of the baroclinic growth mechanism in the model (linked to the optimal coupling of upper-tropospheric potential vorticity advection and low-level thermal advection for this event) as well as a realistic modulation of the simulated storm by the diabatic and sub-grid scale processes.But even the mesoscale simulation, which relies strongly on correct propagation of boundary conditions into the inner domain, generates a deep Mediterranean cyclone within the same span of 48 hours (Figure 17b).Obviously, the much higher horizontal resolution of this experiment endows the surface circulation with additional features, many of these associated with the prominent regional orography.
Storm Hugo (23-24 March 2018;Ruiz et al., 2018) represents a special challenge to the TRAM model for two reasons.First, the precursor disturbance consisted, at surface, of a weak sub-synoptic trough rapidly approaching from the open Atlantic Ocean, and thus a feature necessarily handled by the western boundary conditions for any simulation domain of reasonable size.Second, this low-pressure system experienced an extreme deepening on its way to Europe, ending up as a relatively small, but intense cyclone that inflicted severe winds on the Cantabrian Sea (Spain) and western coasts of France.On the other hand, the medicane Zorbas in the central Mediterranean (28-29 September 2018;Portmann et al., 2020) presents, as all medicanes, similar challenges regarding the sub-synoptic dimension of these disturbances and the subtle characteristics of the larger-scale precursor disturbances.While these tropical-like cyclones are geographically more autonomous than Atlantic storms and most baroclinic lows, since they fully evolve over the Mediterranean Sea and nearby coasts, its prediction is further complicated by the dominant role of surface fluxes and mixing in the maritime boundary layer, moist processes and diabatic heating (e.g., Romero & Emanuel, 2013;Tous et al., 2013).
Hugo and Zorbas were simulated under identical conditions, using a mid-resolution domain (Δx = 25 km) and considering Δz m = 200 m (with 81 vertical levels), stretch = 10, Δt = 45 s and Nstep = 5.These 90-hour-long simulations were nested in the ERA5 meteorological grid reanalyses, updated at six-hour intervals.Starting from initially weak disturbances TRAM succeeds in fully developing both the explosive Atlantic cyclone (Figure 18a, corresponding to a 72-hour forecast) and the medicane (Figure 18b, a 48-hour forecast).The size (≈ 500 km in diameter), intensity and trajectory of both cyclones during their whole life cycle are reasonably well captured, based on inspection of the ERA5 analyses and satellite images (not shown).At the times displayed, Hugo is correctly placed off the Brittany coast, and Zorbas, while heading towards the Ionian Sea, is about to impact southern Greece, where most of the damage was actually inflicted.Additional simulations switching off the latent heat release in the model prove the crucial role of the diabatic factor for the genesis of Zorbas, but also its substantial contribution to the intensification of Hugo (results not shown).
Heavy precipitation and flash flooding are among the most devastating natural hazards in the countries of the Mediterranean basin (see for instance Drobinski et al., 2014).The selected Dana episode of 11-14 September 2019 (Romero-Díaz & Pérez-Morales, 2021) is a perfect example of the typical synoptic setting conducive to catastrophic floods in Valencia and Murcia (eastern Spain) during the late summer and fall.A cut-off low in the mid and upper troposphere combines with a rather shallow depression at surface that continuously supplies warm and moist air from the east to the coastal environment.Cold air aloft and Mediterranean characteristics at low levels guarantee the development of conditional or potential instability in a synoptic setting governed by upward dynamic forcing and water vapour flux convergence (Doswell III et al., 1998).This convective instability will be readily released by the supplementary contribution of mesoscale circulations and by direct mechanical lifting over the coastal slopes under the impinging easterly flow.
Despite the enhanced predictability offered by the complex orography of the Mediterranean basin (e.g., Romero et al., 2005), quantitative precipitation forecasts and the potential for flash flooding will nevertheless remain problematic owing to the dependence of convective systems on multiple -and often uncertain -spatial and temporal scales (Roebber et al., 2004) and Nstep = 5.Once again, this simulation was initialized and fed at the boundaries with the six-hour apart ERA5 reanalyses.The results for the total accumulated precipitation are simply outstanding and in close agreement with the observations, both spatially and in quantitative terms (see Figure 19).The model correctly forecasts amounts exceeding 50 mm in the entire Southeast region of the Iberian Peninsula, and highlights to an acceptable level the flood risk in coastal areas of Valencia and Murcia with two embedded maxima in the range 200-500 mm.
Finally, storm Gloria represents a paradigmatic case in the recent meteorological history of Mediterranean Spain for the broad class of adverse effects it produced on 19-23 January 2020 (Amores et al., 2020).It brought heavy snow at low elevations, excessive precipitation with floods, violent winds and, at sea, storm surges, wind waves and currents of extreme amplitude.The widespread damages in coastal infrastructures and the erosion of beaches and deltas are unparalleled over the last decades.Again, we focus the analysis on the most critical but less predictable variable, precipitation, using the same domain, resolution, input data and numerical configuration as in the previous case.Figure 20 compares the TRAM accumulated precipitation (in 138 hours, starting at 0000 UTC 18 January) against the estimated distribution from data of the AEMET rain-gauge network.TRAM captures the substantial accumulations of this episode along the eastern flank of the Iberian Peninsula and Balearic Islands and reproduces with reasonable accuracy the main torrential centres, where rainfall exceeded 200 mm.Other key aspects of Gloria, like the cyclogenesis that took place to the southwest of the Balearic Islands and the outbreak of intense winds over the maritime areas, are likewise well captured (fields not shown).These types of accomplishments are clear symptoms of an adequate cascade (from synoptic to meso-and local scales) of all relevant kinematic and thermodynamical ingredients in the TRAM simulations.

CONCLUSIONS
At the beginning of this project, a few years ago, we had in mind the development from scratch of a novel atmospheric numerical model aimed at a wide range of time-space scales.Encouraged by the positive results successively achieved throughout all stages of model-building -those steps precisely presented and discussed in this paper -we ended up producing a state-of-the-art numerical model and, undoubtedly, upgrading the modelling capabilities of Meteo-UIB.The combination of the NHFC dynamical set under minimal simplifications, sophisticated numerical integration methods, and physical parametrizations of intermediate complexity, resulted in a versatile modelling system.Indeed, TRAM is suited to simulate processes ranging from thermal bubbles to extratropical baroclinic cyclones, that is, circulations differing in size and life cycle by several orders of magnitude.In between, the model accurately represents the full range of atmospheric waves, flow perturbations and instability types linked to both internal dynamics and external factors (e.g., orography), circulations associated with differential heating (e.g., sea/land breezes and slope/valley winds), cyclonic disturbances whatever their origin -thus including tropical-like cyclones -and the many facets resulting from atmospheric moisture, particularly its crucial role in the genesis of high-impact convective and mesoscale weather systems.
The main technical characteristics of the new TRAM model can be listed as follows: • The dynamical core of the model consists of a classical NHFC version of Euler's equations for the atmosphere which poses predictive equations for the three velocity components and for the perturbations of Exner pressure and potential temperature.Since these equations are not written in flux form, mass and energy cannot be strictly conserved by the numerical scheme.This apparent limitation is generally not problematic in the context of short-to medium-range weather predictions, although predictions of surface pressure might be influenced in some situations (see Thuburn, 2008).
• A mesh of equilateral triangles is used in the horizontal (three-dimensional version of the model), with no staggering of variables (i.e., all predictive fields defined at the barycentres of the triangular cells).Horizontal advection avoids dimensional splitting and is formulated using the REA strategy; the MC slope limiter is applied in the reconstruction phase.
• The classical height coordinate is used in the vertical, but allowing arbitrary stretching of computational levels (e.g., higher resolution in the PBL).The fields are also not staggered in the vertical direction, and an analogous one-dimensional REA method is used to solve advection.A proper treatment of terrain slopes and bottom boundary conditions allows correctly incorporating the effects of the complex orography.
• Time splitting is applied for the time integration of the model, that is, a short time step is used for the fast terms (e.g., terms responsible for the gravity waves and acoustic modes) and a 5-10 times longer time step is used for the slow terms (e.g., advection and physical parametrizations forcings).Additionally, the fast terms are integrated using a second-order Runge-Kutta (RK2) cycle and, in the vertical, instead of explicitly, they are solved semi-implicitly in order to relax the CFL stability condition (i.e., the time step is ruled by the horizontal grid length exclusively).
• Fully coupled with the dynamical core and considering up to six water species, the model includes a realistic set of physical parametrizations of the effects of cloud microphysics, cumulus convection, short-and long-wave radiation, PBL processes and surface fluxes.All moist effects on pressure and thermodynamics are retained in the equations.
• Applications on the real earth use the Lambert map projection.All Coriolis and curvature terms are retained in the corresponding equations.
• No explicit filters are needed in the model to control possible degrading effects of linear and non-linear numerical instabilities.
Despite having reached the most crucial milestones as regards the effective development of a NHFC numerical model for atmospheric applications, the new TRAM model cannot be considered complete.First, we would like to expand the options and complexity of available physical schemes as, currently, only one scheme is incorporated for each parametrized process.Besides initial/boundary conditions, physical parametrizations are recognized as a major source of uncertainty in meso-and convective-scale numerical simulations and, thus, more options available would allow promoting predictability studies and the design of multiphysics ensemble forecasting applications (e.g., Amengual et al., 2021).Another step towards widening the operational scope of the model, by being able to define autonomously well-adapted initial states from observations and analyses, including ensemble members, will consist of implementing some of the different flavours of data assimilation methods (see Kalnay, 2002).More straightforward additions in TRAM, such as the capability of running nested domains in a cascade of horizontal resolutions with either one-or two-way coupling (Madhulatha et al., 2021), should also be accomplished to broaden its applicability and operational perspective.
While we are aware that TRAM improvements, testing and debugging will remain ongoing for an extended period, and that at the end of the process an open-source version -and user guide -for external use should also be released, the model is already an excellent tool to foster locally at Meteo-UIB educational and research activities.It is being run by master's students in a course on 'Geophysical Fluid Simulation' and it has been successfully applied to investigate the synoptic and mesoscale mechanisms that drove recent Mediterranean flooding episodes (papers in preparation).Finally, major scientific-technical efforts have been combined to have TRAM running automatically, twice daily, over three independent domains (continental, regional and insular; square-based resolutions of 17, 6 and 2 km respectively) nested in the Global Forecast System (GFS-NCEP) forecast fields.These real-time applications become, in fact, the best possible testbed to assess in the future the performance of successive versions of the model for all kinds of weather situations and at different scales.The produced high-resolution TRAM forecasts are available to the public and key stakeholders, as they are continuously updated at http://meteo.uib.es/tram.

F
Schematic depiction of the triangle-based horizonal mesh used in the TRAM model.It is composed of two types of equilateral triangular cells, T and B, whose dimensions are defined by the side length dx.A hypothetical advective evolution of two contiguous cells is shown: green trajectories define the new position of the vertexes after the time interval Δt; bluish and reddish triangles are the new configurations of the original T-B pair.These new triangles will affect in varying proportions the neighbouring cells as part of the Reconstruct-Evolve-Average (REA) scheme (see text for details).
Results of the two-dimensional advection test, consisting of the initial scalar field shown in (a) in blue (contour interval 10 units, starting at 1) subjected to a wind field (light blue) described by solid-body rotation (angular velocity is 1 rev/day).Results are shown after four revolutions for (b) Upwind, (c) Lax-Wendroff and (d) monotonized central-difference (MC) slopes in the Reconstruct-Evolve-Average (REA) scheme.Simulations use Δt = 60 s and 'static' winds.The geography of the western Mediterranean is superimposed to appreciate the spatial scale of the test.
a) Numerical error of the simulation by the monotonized central-difference (MC) method shown in Figure2d

F
Numerical simulation of the partial dam break problem, considering an initial state with water depth upstream of the dam of 10 m (dark blue) and, downstream, of 1 m (green).See text for the details of the experiment.
Numerical simulation of the rising thermal bubble (see text for the configuration of both experiments).Results are shown after (a) 17 min in a 'normal' resolution experiment, and (b) 33 min in a 'double' resolution experiment.Note the displayed frames are zoomed in on the rising thermal perturbation, over a partial domain of 20 × 13 km.
shows again the ' evolution (half domain only) at exactly the same moments of the Straka et al. experiment (0, 5, 10 and 15 min, compare with their figure 1) according to a configuration with Δx = Δz = 100 m, F I G U R E 6 Numerical simulation of the interaction of a large warm bubble and a small cold bubble, whose initial characteristics are indicated in the upper-left panel.Following this initial configuration (see text for additional details), results of the interacting thermal perturbations after 4, 7 and 10 min are shown in the other panels.F I G U R E 7 Numerical simulation of the density current, with the sequence of plots showing the evolution of the cold perturbation at 0, 5, 10 and 15 min.See text for the specific characteristics of this experiment.Note that in this representation only the right half of the domain is shown (dimensions: 20 × 8 km).
Δt = 0.25 s and Nstep = 10.The simulation captures in detail the different phases and continuous reshaping of the density current as the negatively buoyant thermal perturbation hits the surface and spreads horizontally.
Numerical simulation of inertial-gravity waves (see text for details).Warm and cold perturbations are shown initially in (a) and after 50 min in (b).The domain size is 300 × 10 km.F I G U R E 9 Schematic depiction of the treatment of orography in the TRAM model.In this two-dimensional example the grid cell size is indicated as the light-orange square near the top-left corner.The resulting terrain mask (same colour) is shown along the bottom part of the figure.In its central part, the distinction between 'grid' slope and 'true' slope is indicated.Dark-orange lines account for the average terrain elevation (see text for details).

F
I G U R E 10 Numerical simulation of linear mountain waves for three different values of the mountain half-width: (a) 10 km; (b) 50 km; and (c) 2 km.Vertical velocity contours are plotted (ascending motion in yellow-red, descending motion in green-blue) with intervals of 0.002, 0.0003, and 0.006 m⋅s −1 respectively.Horizontal dimension of the subdomain shown is 80, 600 and 24 km respectively, while the vertical dimension is 20 km in all cases.(Figure10).These three simulations are performed with the two-dimensional model using Δz = 200 m, Δt = 0.6 s and Nstep = 10 over a domain of 20 km in vertical extent.Other simulation parameters are adapted to each experiment, like the domain horizontal dimension (80, 300 and 24 km respectively), horizontal resolution Δx (1 km, 5 km and 500 m) and τ max (0.1, 0.4 and 0.025 s −1 ).In addition, for the simulation with the widest mountain of 50 km the Coriolis parameter is activated in the model (f = 10 −4 s −1 ) as rotational influences cannot be neglected at these spatial scales.The stationary state is rapidly reached in all these experiments.By comparing our results (Figure10a-c) with identical types of plots contained inKlemp et al., 2008 (see their figs 4, 6 and 7 respectively) the correct performance of the TRAM model in both the hydrostatic and non-hydrostatic limits is confirmed, as revealed by the outstanding reproduction of the analytical solution for the three linear cases.Regarding the Schär mountain wave test (Figure11), this well-known two-dimensional experiment analyses the stationary response of the same previous westerly current F I G U R E 11 Numerical simulation of the Schär mountain wave.Vertical velocity contours are plotted (ascending motion in orange-red, descending motion in turquoise-magenta) with an interval of 0.05 m⋅s −1 (zero contour omitted; see scale).The subdomain shown corresponds to a zoomed area of 50 × 10 km.

F
I G U R E 12 Numerical simulation of the T-REX intense mountain wave.Fields shown after a simulation time of four hours are: the horizontal perturbation wind component (colour field according to scale, with an interval of 5 m⋅s −1 ) and the potential temperature field (black contours, with an interval of 10 K).The domain shown corresponds to an area of 400 × 25 km.
Numerical simulation of two cases of flows past isolated topography: (a) elliptic mountain of 3 km height, using homogeneous initial conditions; and (b) circular mountain of 3 km height, introducing an initial temperature asymmetry.Terrain height is indicated with black contours (interval of 1000 m, starting at 500 m).Fields shown after a simulation time of six hours are: surface wind (as a reference, unperturbed wind vectors correspond to 10 m⋅s −1 ) and surface relative vorticity (colours, ranging from −3 × 10 −3 s −1 to +3 × 10 −3 s −1 from magenta to brown, with yellow corresponding to zero vorticity).Both domains encompass an area of 300 × 139 km.three-dimensional simulations of vortex formation and vortex shedding in continuously stratified flows past isolated topography, similarly toSchär and Durran (1997).These tests are run under Δx = 2 km, Δz m = 500 m, stretch = 2, Δt = 4 s and Nstep = 10 over a domain of 300 × 139 × 20 km.The topographic response is forced through the initial impingement of a uniform flow from the west (u = 10 m⋅s −1 ) with constant vertical stratification (N = 0.01 s −1 ).Coriolis effects are neglected.First, a circular bell-shaped mountain of half width 10 km and height 1.5 km or 3 km was analysed (i.e., inverse Froude numbers are 1.5 and 3 respectively;Smolarkiewicz & Rotunno, 1989).The results are entirely consistent with U R E 14 Numerical simulation of the sea breeze on the island of Mallorca: (a) wind field (vectors every other point; speed in m⋅s −1 , according to scale) at the mature state, 1300 UTC; and (b) breeze-induced daily rainfall (in mm, according to scale).Tramuntana mountain range along the northern coast is indicated using 500-m and 1000-m height contours.Simulation domain is 180 × 155 km.
Numerical simulation of a squall line five hours after its initialization, corresponding to the case with moderate vertical wind shear and: (a) reflectivity (colour field, with an interval of 5 dbZ starting at 5 dbZ) and storm relative wind (vector field); (b) potential temperature perturbation (colours, ranging from −10 K to +10 K from magenta to brown, with yellow indicating zero perturbation) and absolute wind (vector field).In both panels the cloud envelope of the storm is also displayed by means of the grey contour.Regarding the wind fields, it should be noted that a low density of grid values has been used for the clarity of the representation; as reference of vector magnitudes, the u-component's size of the horizontal cell would correspond to 35 m⋅s −1 , the w-component's size of the vertical cell to 5 m⋅s −1 .The domain shown corresponds to a central W-E-oriented vertical cross-section of 600 × 15 km.
Numerical simulation of a supercell 4.5 hours after its initialization, corresponding to the case with the quarter-circle hodograph.(a) Reflectivity at a height of 2 km (in dbZ, according to scale); (b) surface wind field (speed in m⋅s −1 , according to scale); (c) surface potential temperature perturbation (in K, according to scale); (d) rainfall since the beginning of the simulation (in mm, according to scale).As reference, the last three panels include the reflectivity contour of 35 dBZ (grey line).The simulation domain covers an area of 300 × 260 km.

F
to two domains/resolutions, 48 hours after the initialization.(a) Mean sea level pressure (contours, in hPa) and vertically-integrated cloud content (kg⋅m −2 , according to scale); (b) surface wind field (speed in m⋅s −1 , according to scale).
. Here we test TRAM's ability to simulate the full Dana episode over an ultrahigh-resolution domain (Δx = 4.5 km), starting one day before, at 00 UTC 10 September, and for a forecast horizon of 90 hours.The rest of parameters are: Δz m = 200 m (81 vertical levels), stretch = 10, Δt = 9 s cyclonic storm Hugo, and (b) medicane Zorbas.It is shown the mean sea level pressure (contours, in hPa) and vertically integrated cloud content (kg⋅m −2 , according to scale), 72 and 48 hours after the initialization respectively.

F
I G U R E 19 (a) Observed total precipitation (mm) in Valencia and Murcia during the Dana episode of September 2019 ; and (b) numerical simulation with the TRAM model (in mm, according to scale).(source:Spanish agency AEMET) total precipitation (mm) during the Gloria episode of January 2020 ; and (b) numerical simulation with the TRAM model (in mm, according to scale).(source:AEMET)