B Ribstein

C Millet

F Lott

A De

La Cámara
email: acamarai@ucm.es

Can we improve the realism of gravity wave parameterizations by imposing sources at all altitudes in the atmosphere?

Can we improve the realism of gravity wave parameterizations by imposing sources at all altitudes in

the atmosphere?

B Ribstein, C Millet, F Lott, A de La Cámara

Introduction

Atmospheric gravity waves (GWs) have long been observed with radio soundings [START_REF] Tsuda | Radiosonde observations of equatorial atmosphere dynamics over indonesia 2. characteristics of gravity waves[END_REF], using radars [START_REF] Love | Gravity wavemomentum flux in themesospheremeasured by vhf radar at davis, antarctica[END_REF][START_REF] Shibuya | A study of the dynamical characteristics of inertiagravity waves in the antarctic mesosphere combining the pansy radar and a non-hydrostatic general circulation model[END_REF], lidars [START_REF] Baumgarten | Inertia gravity wave in the stratosphere and mesosphere observed by doppler[END_REF][START_REF] Khaykin | Seasonal variation of gravity wave activity at midlatitudes from 7years of cosmic gps and rayleigh lidar temperature observations[END_REF], and satellites [START_REF] Alexander | Global and seasonal variations in three-dimensional gravity wave momentum flux from satellite limb-sounding temperatures[END_REF]. Their amplitude grows as they propagate vertically and they impact the middle atmosphere circulation when they break, even if they have relatively small amplitudes at the source level. The spatial scales of GWs are generally too small to be resolved by general circulation models (GCMs), and the generation, propagation and dissipation of these waves need to be parameterized for GCMs to produce a reasonable circulation. Such parameterizations were first introduced in the 1980's in models with a barely resolved stratosphere, and only highamplitude orographically-forced GWs were needed to be taken into account [START_REF] Palmer | Alleviation of a systematic westerly bias in general circulation and numerical weather prediction models through an orographic gravity wave drag parametrization[END_REF]. Nowadays, most models resolve the middle atmosphere requiring the parameterization of the effects of smaller-amplitude, non-orographic gravity waves [START_REF] Manzini | Impact of the doppler spread parameterization on the simulation of the middle atmosphere circulation using the ma/echam4 general circulation model[END_REF].

One way to parameterize non-orographic GWs consists in using the observational evidence that over a large number of realizations the GW fluctuations of vertical wind and temperature in the middle atmosphere follow a "universal" spectra, which shape is derived from radiosondes and satellite data (e.g., Cot, 2001;[START_REF] Zhao | Lidar observations of stratospheric gravity waves from 2011 to[END_REF]. These spectra are numerically robust [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], and various theories have been developed to explain them. Some involve wave breaking (e.g., [START_REF] Dewan | Saturation and the universal spectrum for vertical profiles of horizontal scalar winds in the atmosphere[END_REF], whereas other include nonlinear effects like triade interactions, Doppler spreading, and inverse cascades (e.g., [START_REF] Broutman | On the interaction of small-scale oceanic internal waves with near-inertial waves[END_REF][START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF][START_REF] Hines | Nonlinearity of gravity wave saturated spectra in the middle atmosphere[END_REF][START_REF] Métais | Inverse cascade in stably stratified rotating turbulence[END_REF]. Beyond the theoretical debate, a practical result is that the existence of a saturated spectra allows semi-theoretical integrations that ease the parameterization of non-orographic GWs (C. O. [START_REF] Hines | Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. part 2: Broad and quasi monochromatic spectra, and implementation[END_REF][START_REF] Warner | On the propagation and dissipation of gravity wave spectra through a realistic middle atmosphere[END_REF][START_REF] Manzini | Impact of the doppler spread parameterization on the simulation of the middle atmosphere circulation using the ma/echam4 general circulation model[END_REF]. To a certain extent, this approach is challenged by the recent balloon observations showing that the GW field is very intermittent, and is often dominated by rather well-defined GWs packets [START_REF] Hertzog | Estimation of gravity wave momentum flux and phase speeds from quasilagrangian stratospheric balloon flights. part ii: Results from the vorcore campaign in antarctica[END_REF][START_REF] Wright | Global observations of gravity wave intermittency and its impact on the observed momentum flux morphology[END_REF][START_REF] Alexander | Global and seasonal variations in three-dimensional gravity wave momentum flux from satellite limb-sounding temperatures[END_REF]. This intermittent nature makes that in each model gridbox and at each time step, the number of GWs packets is not large enough to fulfill the law of large numbers underlying the construction of spectra out of individual realisations. More generally, this concern is related to that of statistical equilibrium and is central in the recent development of stochastic parameterizations [START_REF] Berner | Stochastic parameterization toward a new view of weather and climate models[END_REF]. In the case of the GWs this implies that the "universal" spectral shape should be checked a posteriori and over a large number of days, rather than being realised every time. This leads to an alternative approach to parameterize non-orographic GWs, which is based on representing the GW field with a Fourier series in the horizontal and temporal directions (e.g., [START_REF] Alexander | A Spectral Parameterization of Mean-Flow Forcing due to Breaking Gravity Waves[END_REF]. In this approach the individual harmonics are a crude representation of the individual wave packets and the intermittency is taken into account by launching stochastically a few harmonics each timestep [START_REF] Eckermann | Explicitly Stochastic Parameterization of Nonorographic Gravity Wave Drag[END_REF][START_REF] Lott | A stochastic parameterization of nonorographic gravity waves: Formalism and impact on the equatorial stratosphere[END_REF]. The challenge is then to reconcile the two types of schemes ("spectral" versus "multiwave") and the two types of observations (stationary universal "vertical spectra" versus intermittent "wave packets"). As we shall see, this can be done by showing that the ensemble average of the periodograms associated with superposed harmonics can reproduce the observational "universal" spectra (de la [START_REF] De La Cámara | Intermittency in a stochastic parameterization of nonorographic gravity waves[END_REF]. In route to realize this objective, the study of [START_REF] Souprayen | Atmospheric gravity-wave spectra: A stochastic approach[END_REF]) is encouraging since it shows that the wave filtering by the large scale flow and the breaking of individual GWs packets can yield realistic spectra in the upper stratosphere and mesosphere.

In some modeling centers the amplitude of the parameterized GWs is related to their non-orographic sources, i.e. convection and fronts and/or flow imbalance, and the GWs are launched from a single source level in the troposphere. This last choice is not well justified for several reasons. One is that the presence of unbalanced flows that emit GWs is not restricted to the troposphere: [START_REF] Dörnbrack | Gravity waves excited during a minor sudden stratospheric warming[END_REF] and [START_REF] Sato | Gravity wave generation around the polar vortex in the stratosphere revealed by 3-hourly radiosonde observations at syowa station[END_REF] provide observational evidence for GW generation in the stratospheric polar vortex, and [START_REF] Polichtchouk | Spontaneous inertia-gravity wave emission from a nonlinear critical layer in the stratosphere[END_REF] discuss the generation from the critical layer of the stratospheric polar night jet using an idealized numerical model. A second reason is that observations often show that the GWs in the middle atmosphere have rather long periods (often around 6 hours and longer in [START_REF] Reichert | Retrieval of intrinsic mesospheric gravity wave parameters using lidar and airglow temperature and meteor radar wind data[END_REF], near the inertial period in [START_REF] Gelinas | Observations of an inertial peak in the intrinsic wind spectrum shifted by rotation in the antarctic vortex[END_REF], [START_REF] Bellenger | Tropospheric turbulence over the tropical open ocean: Role of gravity waves[END_REF], [START_REF] Shibuya | A study of the dynamical characteristics of inertiagravity waves in the antarctic mesosphere combining the pansy radar and a non-hydrostatic general circulation model[END_REF][START_REF] Shibuya | A study of the dynamical characteristics of inertiagravity waves in the antarctic mesosphere combining the pansy radar and a non-hydrostatic general circulation model[END_REF]Vincent and[START_REF] Vincent | Balloon-borne observations of short vertical wavelength gravity waves and interaction with qbo winds[END_REF]). The presence of these slow waves is difficult to justify if the GW sources are only in the troposphere, as waves with small intrinsic frequency have short vertical wavelengths and saturate more easily than faster waves. This process is often referred to as dynamical filtering and occurs systematically when the waves approach a critical level. For GWs, this dynamical filtering is very effective, explaining most of the relation between the wind speed and GW amplitude in balloon measurements [START_REF] Plougonven | On the relation between gravity waves and wind speed in the lower stratosphere over the southern ocean[END_REF].

The present paper analyses if the "multiwave" non-orographic GW parameterization due to fronts and jets of de la Cámara and Lott (2015) can be adapted to include sources from all levels, small intrinsic phase speed waves (including near inertial waves), and if these modifications can help to reproduce the universal spectral shape systematically enforced in the "spectral schemes". Although some of these concerns could be adressed with any other scheme, this one has few characteristics that makes it suitable to treat all of them. The first and most important is that this scheme is based on a spontaneous emission theory [START_REF] Lott | Gravity waves generated by sheared three-dimensional potential vorticity anomalies[END_REF], which is a theory that partly explains the GW emission seen in quite sophisticated high resolution simulations [START_REF] Polichtchouk | Spontaneous inertia-gravity wave emission from a nonlinear critical layer in the stratosphere[END_REF], and which is not limited to the troposphere. The second is that it points to PV anomalies as a source of GWs, which is coherent with the fact that processes such as "classical" geostrophic adjustement or re-emission are associated with the presence of PV anomalies. The third is that it is operational in the Institute Pierre Simon Laplace (IPSL) Earth system model [START_REF] Hourdin | Impact of the lmdz atmospheric grid configuration on the climate and sensitivity of the ipsl-cm5a coupled model[END_REF], so it is routinely tested. The fourth is about observational constraints, in the sense that it qualitatively produces the observed intermittency of the nonorographic GW field [START_REF] Hertzog | On the Intermittency of Gravity Wave Momentum Flux in the Stratosphere[END_REF][START_REF] Wright | Global observations of gravity wave intermittency and its impact on the observed momentum flux morphology[END_REF][START_REF] De La Cámara | Intermittency in a stochastic parameterization of nonorographic gravity waves[END_REF][START_REF] Alexander | Global and seasonal variations in three-dimensional gravity wave momentum flux from satellite limb-sounding temperatures[END_REF]. Indeed, de la Cámara, [START_REF] De La Cámara | On the gravity wave forcing during the southern stratospheric final warming in lmdz[END_REF] demonstrated that a good representation of the GW intermittency can be beneficial for models, helping the IPSL model to better simulate the timing of the Southern Hemisphere stratospheric final warming. A final (and much less positive) reason is that in its operational version, de la Cámara and Lott (2015) use a Gaussian distribution of intrinsic phase speeds with standard deviation near 40ms -1 . This large value is used because it helps the waves to propagate up to the upper mesosphere without attenuation, but it contradicts the fact that frontal waves resulting from spontaneous adjustment have small intrinsic phase speeds near their source [START_REF] Lott | Gravity waves generated by sheared three-dimensional potential vorticity anomalies[END_REF].

A central assumption made in this paper is that this bias toward larger than expected intrinsic phase speeds is common to other schemes, and that trying to correct it in one scheme could be indicative of what could be done in other schemes. In fact the corrections we test in this study have a general character: we include the Coriolis force because it can be significant at low intrinsic frequencies, and we place sources at all levels rather than in the troposphere only. Finally, it is worth noting that the shortcomings we try to deal with are today considered as priorities in the community (see discussion about low phase speed waves in [START_REF] Alexander | Observational validation of parameterized gravity 2 waves from tropical convection in the whole 3 atmosphere community climate model (waccm)[END_REF]), but are not the only ones. Some authors place more emphasis on including three-dimensional propagation of gravity waves in parameterizations (e.g., [START_REF] Muraschko | On the application of wentzel-kramer-brillouin theory for the simulation of the weakly nonlinear dynamics of gravity waves[END_REF]Ribstein & Achatz, 2016;[START_REF] Amemiya | A new gravity wave parameterization including three-dimensional propagation[END_REF]. Including lateral propagation in highly parallelized code is extremely challenging computationally, so it is worthwhile to test if some improvements can be obtained through other routes.

The goal of this paper is to present modifications so that the frontal GWs parameterization can include slow intrinsic phase speed waves, rotation and GWs sources at all levels in the atmosphere. Section 2 describes the modifications we propose to the frontal GWs parameterization used in LMDz. In section 3, offline tests are performed to test if the scheme predicts realistic vertical energy spectra. Section 4 presents online results obtained with the LMDz GCM, using the standard parameterization and the updates.

Section 5 gives the main conclusions.

2 Non-orographic gravity waves due to fronts and flow imbalance

General formalism

We next summarize the formalism of the stochastic parameterization used in LMDz (de la [START_REF] De La Cámara | A parameterization of gravity waves emitted by fronts and jets[END_REF], and emphasize the modifications introduced in this study.

The horizontal wind and temperature disturbances (u , T ) due to GWs are represented by a stochastic Fourier series of J monochromatic waves,

(u , T )(x, z, t) = J j=1 C j ûj (z), Tj (z) e i(kj •bf x-ωj t) (1)
whose horizontal wavevector k j and absolute frequency ω j are chosen randomly, and the complex amplitudes ûj (z), Tj (z) vary in the vertical direction measured by the logpressure coordinate z. The intermittency coefficient C j measures the probability of the presence of the corresponding wave at a given horizontal location within the gridbox. In previous versions of the parameterization we had always assumed equiprobability for simplicity and taken

J j=1 C 2 j = 1, i.e. C j = 1/ √ J. (2)
As we shall see, the value of this parameter can be changed as to represent sources located at different vertical levels.

To evaluate the wave amplitude, we adapt Lott, Plougonven, and Vanneste (2012)'s analytical estimate of the GW momentum flux emitted by a potential vorticity (PV) anomaly in a vertically sheared flow, and consider that a given model level z l of thickness dz l emits a vertical Eliassen-Palm (EP) flux with value,

F(z l ) = G 0 ∆zdz l 4f 0 ρ(z l )N (z l ) f tanh ζ(z l ) f 2 e -π N (z l ) Uz (z l ) , (3)
where U z = |U z | is the modulus of the vertical shear, N is the Brunt-Väisälä frequency, ρ = e -z/H , and H = 7km is the characteristic vertical scale of density decay. Still in equation ( 3), G 0 is a tuning parameter, f 0 = 10 -4 s -1 is a characteristic value of f , ∆z a characteristic depth of the subgrid scale PV anomalies. To derive equation ( 3) from theory we made the approximation that the subgrid scale vorticity equals the gridscale one ζ. The only novelty at this stage is the hyperbolic tangent term that is used to limit the relative vorticity to values below |f |. This is a reasonable assumption since flows with relative vorticity larger than |f | are likely to be strongly unstable, and it also moderates the emission intensity in the tropics.

de la Cámara and Lott (2015) assumed that the momentum fluxes essentially come from the troposphere, but integrated the contributions in equation ( 3) over all model levels L, to calculate a total emitted flux,

F = L l=1 F(z l ). (4)
It is then imposed that this momentum radiates from a specified launching level z La in the troposphere, distributed over the ensemble of J monochromatic waves in (1). Among the randomly chosen parameters, the direction of each harmonic F j of the total EP-flux F is specified through the random horizontal wavector k j , following the rule

F j (z La ) = - k j k j F. (5)
To derive ( 5) from ( 4) we have chosen by convention that the sign of the intrinsic phase speed at the launch level

ωj (z La ) = ω j -k j • U(z La ) > 0. (6)
To evaluate the vertical profile of the flux above the launching level we consider that from one vertical level to the next above, the flux is (i) reduced by a small diffusivity µ/ρ 0 (z), (ii) limited by that of a saturated wave (e.g., [START_REF] Lindzen | A note on the limits of rossby wave amplitudes[END_REF], and (iii) set to zero immediately above an inertial level [START_REF] Lott | Inertia-gravity waves in inertially stable and unstable shear flows[END_REF]:

F j (z l+1 ) = - k j k j Θ ωj (z l+1 ) 2 -f 2 × min F j (z l ) exp 2 µm 3 j (z l ) ρ 0 ωj (z l ) dz l , ρ r S 2 c N (z l+1 ) 2 |m j (z l+1 )| 3 k 2 min k j . (7)
In equation ( 7), Θ is the Heaviside function to handle inertial levels, k min is the horizontal wave number associated with the largest resolved GW in the model (it is related to the model's horizontal resolution), and S c is a tunable parameter that controls the saturation amplitude. Still in equation ( 7), the vertical wavenumber m j and intrinsic frequency ωj are given by

m j (z) = - N (z) k j ωj (z) 2 -f 2 , where ωj (z) = ω j -k j • U(z), (8)
where the minus sign in the definition of m j ensures upward propagation above the launching level. A novelty is that we have included the Coriolis force in equations ( 7) and ( 8).

At this stage, the emitted flux is equidistributed among all possible horizontal directions, which somehow contradicts the theory saying that the preferential emission is in the direction opposite to the wind shear. In practice, however, waves emitted with phase speeds in the direction of the shear get their intrinsic phase speeds decreased and their vertical wave numbers increased as they are evaluated at the next vertical level above.

According to (7) the fluxes for these waves are much more reduced than for the waves in the direction opposed to it. In the scheme we also do exclude the emission of highly saturated waves, which also avoid imposing a huge drag just above the launching level, by systematically replacing the launching value of the flux by its value at the level above (in z La+1 = z La +dz La ). This naturally tends to reduce emissions of waves with phase speeds in the direction of the shear. Finally, below the launching altitude the fluxes are kept constant, which technically allows to define F j (z) at all model levels. This last choice is consistent with the fact that in the model we do not extract momentum from the source's surroundings region to balance the emitted wave drag, a shortcoming that is justified by the rapid decrease of air density with altitude: the corrections to the fields and to the tendencies would presumably be of small amplitude around and below the launch level.

Once F j (z) is evaluated at all vertical levels, we use the Wentzel-Kramers-Brillouin (WKB) formalism to relate the EP flux to the disturbance fields in (1):

ûj (z) = -k j ωj -if e z × k j ω2 j -f 2 φ j (z)e i z 0 mj (z )dz +iξj , (9)
Tj (z) = im j H R φ j (z)e i z 0 mj (z )dz +iξj (10)
where the amplitude of the geopotential is

φ j (z) = F j N ρm j k j , (11)
and ξ j is a phase with no effect on the EP fluxes, which is chosen randomly when computing the physical fields offline. For completeness, note that the GW drag is computed as

∂u ∂t | GW = J j=1 C 2 j ∂F j ∂z . ( 12
)
Thermal effects are taken into account by evaluating the work performed against the wind:

∂T ∂t | GW = -u • ∂u ∂t | GW /C p , (13)
where C p is the heat capacity of dry air.

At this point, we have applied the described equations to J waves emitted at a single launching level, z La . To adapt our formalism in order to work with J waves emitted from different launching levels we choose the launching level randomly, z j , together with the horizontal wavevector k j and intrinsic frequency ωj (see Fig. 1). Then the launching flux in ( 5) is replaced by

F j (z j ) = - k j k j F(z j ), (14)
the emitted amplitude F at the corresponding level being directly given by (3). The various profiles are then evaluated above and below z j following equations ( 7) to ( 11), but with z La replaced by z j for each waves. But now that we have J different launching altitudes, with the possibility that J = L, we need to take for intermittency parameter:

C 2 j = L/J. (15)

Model and reanalysis

In eterizations [START_REF] Lott | A new subgrid-scale orographic drag parametrization: Its formulation and testing[END_REF][START_REF] Lott | A stochastic parameterization of the gravity waves due to convection and its impact on the equatorial stratosphere[END_REF], the different setting of the parameterization of the GWs due to fronts and jet imbalances are detailed in the next subsection. The reader is referred to [START_REF] Hourdin | Lmdz6a: The atmospheric component of the ipsl climate model with improved and better tuned physics[END_REF] for a comprehensive description of the LMDz model equations and specifications of its grid.

Gravity waves parameterization setup

In all experiments discussed in this paper, we consider J = 48 waves each physical time step in the model as well as in the offline reconstructions. We also choose the horizontal wavenumber amplitude k j = ||k j || with uniform probability in the interval k min ≤ k j ≤ k max , with 2π/k max = 6.3km and 2π/k min = 315km crudely bounding the smallest horizontal wavelengths that can be attributed to GWs, and the largest horizontal wavelength that cannot be represented in the model gridbox. The direction of horizontal propagation θ j (k j = k j (cos θ j i + sin θ j j)) is chosen uniformly within the interval 0 ≤ θ j ≤ 2π. The attribution of frequency is indirect since we first select the wave intrinsic phase speed at the launch level z l from a half-normal distribution with a standard deviation of c ϕ , and in the direction given by k j . The parameter c ϕ is key, and tuning it in different experimental setups requires adjustments in the launched momentum flux amplitude and saturation parameters G 0 and S c in (3) and ( 7) respectively.

In this study we present results from three different setups of the nonorographic GWs parameterizations. The strategy adopted consists of (1) decreasing the phase speed drastically and ( 2) introduce multiple level sources adapting the other parameters to give performances that are quite comparable to the existing operationnal scheme. We therefore target rather neutral effects on the model climate, which is in itself a task that demands a substantial amount of trial simulations.

In the first experimental setup, we proceed as in de la Cámara and Lott ( 2015 we see that we are essentially taking into account long waves (since c ϕ is large), and we will refer to this setup as LW.

In the second experimental setup we keep a single launching level but decreases the intrinsic phase speed down to c ϕ = 10ms -1 . This requires a slight increase in G up to G = 5, but since the characteristic vertical scale is now much shorter, i.e. 2πc ϕ /N ≈ 4km, the saturation parameters need to be increased substantially up to S c = 2.5 to keep the saturated flux of about the same amplitude within the middle atmosphere. Note that the characteristic vertical scale is closer to the wavenumber m * introduced by [START_REF] Warner | On the propagation and dissipation of gravity wave spectra through a realistic middle atmosphere[END_REF]. These changes in the wave parameters imply parameterized waves of shorter vertical scales, so this setup will be referred to as SW-1L.

In the third experimental setup, we consider waves emitted from multiple model levels. Nevertheless we realized that as the amplitude of emission in equation ( 3) is very sensitive to the Richardson number and to the relative vorticity ζ, there is often one level that dominates the sum of EP fluxes in equation ( 5). Finding this level out of J cases is equivalent to average over J levels, so there is no need to change substantially the tuning parameters G and S c compared to the previous case, and we take for the new values G = 4.2 and S C = 3. This configuration will be referred to as SW-ML.

Parameterized gravity wave spectra

The performance of the three configurations LW, SW-1L, and SW-ML is addressed in off-line runs using the meteorological fields from MERRA2 as input. 2015)).

Next we analyze the ability of LW, SW-1L and SW-ML to reproduce the empirical GW energy spectra. The vertical profiles of the wind and temperature are obtained from the parameterized profile of vertical EP flux as in equations ( 9), ( 10), and (11). More specifically, and to build spectra out of individual realisations we first construct 1000 monochromatic waves, each corresponding to randomly chosen values of phase speed and horizontal wavenumber. The wind and Temperature profiles of each wave are then sampled every 100m, which is a much higher vertical resolution than that used in MERRA2. We therefore linearly interpolate the EP flux to the target grid, together with the large scale fields of wind and Temperature needed in Eqs ( 9) and (10). To construct 1 realization out of these 1000 monochromatic waves, we pick J = 48 of them randomly, choose the phase ξ i of each randomly and then sum over the J waves. Let u (z) be one such realization of the horizontal wind disturbance, we use here the same notation as in ( 1), we then perform a Fourier analysis of this realisation which gives the periodogram ûû * , û being the Fourier coefficients and the stars indicating conjugation. To avoid numerical artifacts in the boundaries, a tapered cosine window is used together before the fast Fourier transform. The spectra are then estimated at each horizontal places by doing an ensemble average of individual periodograms ûû * constructing an ensemble of 20 independent realisations. Although we found our results about spectral shape to be little sensitive to the procedure, we could have averaged at a given place over different days, or average the same day at different places, we adopt this one because it permits to compare without ambiguity regions with presumably large GW emission to regions with small GW emission (for instance the locations A and B in Figures 2a-b and the same day).

In the following we restrict the discussion to the horizontal wind spectra, but we verified that the temperature spectra exhibit a similar shape (not shown). The spectra for altitudes lower than 25 km and 65 km are displayed in Figure 3 at the two specific locations A and B and still for the 22 january 2012. For both locations the energy spectra obtained with the configurations LW and SW-1L are proportional to m -4.5 , the shape obtained with vertically distributed GW sources (i.e. SW-ML) is characterized by a m -3 tail, which is suggestive that saturation occurs much more systematically under this configuration (e.g., [START_REF] Dewan | Saturation and the universal spectrum for vertical profiles of horizontal scalar winds in the atmosphere[END_REF]. Besides, the vertical spectra for SW-1L and SW-ML are shifted toward smaller wavelengths, which is consistent with the fact that GWs are launched with slower phase velocities.

To summarize, when smaller intrinsic phase speed are imposed (i.e. in SW-1L and SW-ML), there is more energy concentrated at shorter wavelengths, this energy corresponds more often to saturated waves, and the effect of the saturation on decreasing the wave amplitude is compensated in SW-ML by launching waves from all altitudes. These results demonstrate that multiwave schemes with small intrinsic phase speed can fairly reproduce the observed GW energy spectra, bridging a gap with the spectral schemes that include sources but prescribe a saturated spectrum [START_REF] Bushell | Parameterized gravity wave momentum fluxes from sources related to convection and large-scale precipitation processes in a global atmosphere model[END_REF].

4 Impacts on the simulation of the stratosphere

We next evaluate the performance of the three configurations of the frontal GW parameterization in 15-year runs with the climate model LMDz. Figures 4a-b show latitudeheight cross-sections of zonal-mean zonal wind in LW-LMDz (black contours), seasonally averaged for December-January-February (DJF) and June-July-August (JJA). Wellknown features of the zonal mean structure of the troposphere and middle atmosphere are captured in LW-LMDz. The subtropical jets in the upper troposphere are stronger in the winter than in the summer hemisphere, and displaced further poleward in the summer hemisphere. In the middle atmosphere, there are westerly winds in winter (i.e., the polar night jet) and easterlies in summer. The color shading in Figures 4a-b represents the difference of zonal wind between MERRA2 (period 1996-2010) and LW-LMDz. In the tropics, there are differences above 40-km height in all seasons that are related to the representation of the amplitude of the semiannual oscillation [START_REF] Lott | A stochastic parameterization of the gravity waves due to convection and its impact on the equatorial stratosphere[END_REF][START_REF] Smith | The Semiannual Oscillation of the Tropical Zonal Wind in the Middle Atmosphere Derived from Satellite Geopotential Height Retrievals[END_REF]. In the extratropics, the strongest bias takes place in the upper stratosphere and mesosphere in JJA in the SH (Figure 4b), where the weaker winds in LW-LMDz than in MERRA2 are related to a polar night jet in the model that does not tilt equatorward with height as compared to reanalysis. In the NH in DJF (Figure 4a), the westerly winds in LW-LMDz are also systematically weaker than in MERRA2, which implies a weaker polar vortex.

The middle and bottom panels of Figure 4 show the corresponding zonal wind profiles (black contours) for SW-1L-LMDz (Figures 4c-d) and SW-ML-LMDz (Figures 4ef), with the color shading displaying differences with respect to LW-LMDz. Both SW-1L-LMDz and SW-ML-LMDz produce a stronger westerly jet in DJF in the upper stratosphere and mesosphere than LW-LMDz, something consistent with MERRA2. In general, the differences between SW-1L-LMDz and LW-LMDz do not necessarily imply reduced biases (cf. Figures 4a-b), but SW-ML-LMDz performs qualitatively better. This is particularly evident in the case of the polar night jet in the SH in JJA (Figure 4f), which tilts towards the tropics with altitude and the wind differences with LW-LMDz are similar to those between MERRA2 and LW-LMDz (Figure 4b).

To test if these changes can be associated with changes in GW drag, Fig. 5 shows the cross-sections of zonal mean non-orographic GW drag in SW-1L-LMDz and SW-ML-LMDz (black contours), and the corresponding differences with respect to the control run LW-LMDz (shading). The non-orographic GW drag is more effective at higher altitudes in the upper stratosphere and mesosphere, where we find positive drag in the summer hemisphere and negative drag in the winter hemisphere, contributing to decelerate the easterlies and westerlies, respectively. For SW-1L and SW-ML the differences of GW drag and zonal mean zonal wind with respect to LW are generally consistent with this view since regions of slower (larger) mesospheric winds roughly correspond to regions of smaller (larger) GW drag in Figure 5 (see the summer mesosphere between 60km< z < 80km and 80 • S-60 • S in Figs. 4c and5c). As we have modified the phase speeds substantially, this correspondence between GW drag and zonal wind differences could be simply related to changes in dynamical filtering at high altitude. In the regions mentioned before, this does not seem to be the case: weaker winds in summer should allow the propagation of more waves with positive intrinsic phase speeds that would strengthen the net positive drag, which is not what is found in summer in the upper mesosphere. Therefore, it would seem that the new sets of parameters are simply producing less GW drag and acceleration.

Dynamical filtering is more evident in the differences between the SW-ML run and the LW run and at the places where the differences in zonal wind are the more pronounced -i.e. at altitudes between 60km and 80km and around 30 • in each hemisphere (Figures. 4( between the two sides of the shear zone is essential, it reveals that we do not put globally more or less drag, as was the case in summer for SW-1L, but that we distribute it differently according to the wind speed. Another important thing to notice is that the latitudes around 30 • also correspond to the latitudes of the winter tropospheric jet center at z ≈ 15km. We can therefore speculate that since the wind shear is negative on the upper flank of the jet, positive phase speed waves are produced, they create acceleration right above at mesospheric levels, but are not much efficient further above because they have small intrinsic phase speeds and are easily filtered out by the wind they produces. Our explanation is therefore that the stratospheric jet tilt in the SH, which is reproduced in SW-ML-LMDz, is in part supported by GWs with small intrinsic phase speeds coming from the upper troposphere and lower stratosphere.This can be realized by a GW scheme with sources that are not confined to the troposphere.

Overall, the mean middle atmospheric circulation with the two new settings of the frontal GW parameterization (i.e. SW-1L and SW-ML) is reasonable, and presents specific improvements in the shape of the austral winter jet that look promising to correct long-standing biases. To gain further confidence into the performances of the model in these new settings, we next evaluate the interannual variability of the different model runs. Figure 6 shows the annual cycle of the zonal mean temperature at 80 • N and the variability over the period 1996-2010 quantified through the 5th and 95th percentiles (shading). Starting with the control run LW-LMDz (red line and shading in Figure 6a), the model has a clear warm bias at 10 hPa as compared to MERRA2 (black line and shading in Figure 6a, and Fig. 6d), which is consistent with a weak wind bias (Figure 4a).

This temperature bias is reduced using the alternative settings SW-1L and SW-ML (Fig-

ures 6b-c), both in the mean cycle and the variability given by the percentiles. However, the three model simulations clearly overestimate temperature fluctuations in early and late winter.

Conclusions

Specifying the emission of GWs based on the model grid-scale dynamics, and using our observational knowledge of the GW field to constrain tunable parameters in the process, are two major challenges for improving parameterizations of GW drag [START_REF] Alexander | Recent developments in gravity-wave effects in climate models and the global distribution of gravity-wave momentum flux from observations and models[END_REF][START_REF] Plougonven | How does knowledge of atmospheric gravity waves guide their parameterizations? Q[END_REF], particularly for GWs generated by fronts and jet imbalance [START_REF] Plougonven | Internal gravity waves from atmospheric jets and fronts[END_REF]. The parameterization introduced in de la Cámara and Lott (2015) was based on a formalism that can be used to address these challenges.

For this purpose we introduce the following modifications to this parameterization: 1)

we reduce the horizontal intrinsic phase speeds of the launched waves following sugges-tions from observations and high-resolution model runs [START_REF] Plougonven | On the relation between gravity waves and wind speed in the lower stratosphere over the southern ocean[END_REF] (SW-1L configuration); and 2) we launch GWs from all model levels (SW-ML configuration) instead of launching from only one tropospheric level. This can result in more saturated waves everywhere in the middle atmosphere.

A technical result of our study is that launching GWs at different model levels can be achieved at a moderate numerical cost, in our case the cost of the stochastic parameterization of non-orographic GWs used in LMDz. From a more scientific point of view, some results are worth highlighting. The first is that decreasing the intrinsic phase speeds and launching GWs from all levels shifts the slope of the vertical energy spectra of parametrized GWs toward the observed value of -3. Second, the middle atmospheric circulation in climate model simulations responds reasonably well to the applied changes. There are even some indications of a weak reduction of model biases, such as an improved equatorward tilt of the austral polar night jet and stratospheric polar temperatures over the Arctic. These bias reductions may be model dependent, or could have been obtained in the same model with further refined tuning of the initial scheme. However, we obtain these results using configurations with about the same amount of wave stress launched (see Figure 2c) and only moderate alterations of the parameters of the scheme. This may indicate that the improvements obtained are not simply a signature of refined tuning.

Despite the modifications performed to the GW parameterization, placing inertial gravity waves in the middle atmosphere suggested by observations [START_REF] Gelinas | Observations of an inertial peak in the intrinsic wind spectrum shifted by rotation in the antarctic vortex[END_REF] has proven difficult (not shown). This may call for substantial changes in the formalism, perhaps in the source amplitude specification, the latter only including spontaneous adjustment. "Classical" geostrophic adjustment of the kind arising after short gravity waves breaking (secondary emission) could also be considered [START_REF] Vadas | Mechanism for the generation of secondary waves in wave breaking regions[END_REF][START_REF] Lott | Large-scale flow response to short gravity waves breaking in a rotat-ing shear flow[END_REF].

In any case, the choice of including sources at all model levels stays a valid option. Secondary emission, for instance, generates waves at altitudes where primary GWs break, and this usually takes place well in the middle atmosphere.

Current efforts of GW paramaterization development are generally focused on adding complexity in the way parameterizations treat wave propagation and dissipation (e.g., [START_REF] Kim | Towards transient subgrid-scale gravity wave representation in atmospheric models. Part II: Wave intermittency simulated with convective sources[END_REF]Plougonven et al., 2020, and references therein). Our results demonstrate that modest changes in source-specifications guided by GW observations can have desirable effects in both the realism of the energy spectra and the simulated middle at-mospheric circulation in a climate model. These changes are much less computationally demanding than the development of schemes allowing lateral propagation (e.g., Bölöni et al., 2016). Since an argument for including lateral propagation is that columnar propagation hinders GWs to propagate to specific regions where they are needed, our approach provides a stopgap to represent the desirable waves.

de la Cámara, A., Lott, F., & Abalos, M. (2016). Climatology of the middle atmosphere in lmdz: Impact of source-related parameterizations of gravity wave .

  offline mode we use daily fields of temperature and horizontal winds from the Modern-ERa Retrospective Analysis for Research and Applications version 2 (MERRA2) (Gelaro et al., 2017) at 1 • × 1 • resolution and from the ground to 0.01hPa. In online mode we use the stratospheric version of the Laboratoire de Meteorologie Dynamique Zoom (LMDz) model with 142×144 uniform latitude longitude grid and L = 80 vertical levels, the model top being at 0.01hPa. The simulations are forced with the observed seasonal cycle of sea surface temperatures and sea-ice from the CMIP5 database for the period 1980-1995, and the ozone climatology is built from the ACC/SPARC ozone database. All runs have the same settings of the orographic and convective GWs param-

  ) and choose the launching altitude at z La = 500 m, and take c ϕ = 50 m•s -1 , G = 4, and S c = 0.6. These values stay reasonably close to those used in previous studies (e.g., de la Cámara, Lott, & Abalos, 2016) considering that we now include the Coriolis force and bound the disturbance vorticity to values below f . If we consider that the characteristic vertical scale of the waves produced with this setup is 2π/m ≈ 2πc ϕ /N ≈ 20km,

  Figures 2a-b show the spatial distribution on a given day (22 January 2012) of the daily average of vertical EP flux amplitude in SW-ML at the 500 hPa level. Peak values of around 50 mPa are found in the vicinity of the subtropical jets in both hemispheres, reflecting the direct relation between the emitted EP flux and the grid-scale relative vorticity, as indicated by equation (3). Figure 2c compares the latitudinal distribution of the zonal-mean EP flux for LW, SW-1L, and SW-ML. The three curves present very similar features with very low values in the tropics, a maximum around the subtropical jets, and a gradual descent towards the poles. In general, we see that SW-ML represents slightly larger fluxes at all latitudes than the other two setups, and that SW-1L consistently works with smaller fluxes. It is important to emphasize that although the amplitude of the EP flux is calculated deterministically, the probability density function of the launched fluxes qual-itatively follows the observed log-normal distribution in the lower stratosphere (not shown but see de la Cámara and Lott (

  ). The difference in GW drag at those places (Figures. 5(e)-(f)) shows a weaker drag in SW-ML below the jet shear zone, and a stronger drag above. The difference in sign

Figure 1 .

 1 Figure 1. Schematic representation of the scheme used in (de la Cámara & Lott, 2015) (left) and of its modification to include waves emitted from all model levels (right).

Figure 2 .

 2 Figure 2. (a,b) Horizontal distribution of the total EP flux (mPa) at pressure level 500 hPa in the SW-ML setup. The geographical locations A and B are the locations where the energy spectra are constructed as shown in Figure 3. (c) Zonal mean EP flux at 500 hPa for the three GWs configurations.

Figure 3 .

 3 Figure 3. Zonal wind spectra from realizations of the GWs fields obtained using parameterizations LW, SW-1L and SW-ML, respectively, for locations A (a,c) and B (b,d) (see figure 2) and altitude ranges 0-65 km (a,b) and 0-25 km (c,d).

Figure 4 .

 4 Figure 4. Cross-sections of zonal mean zonal wind from the three experiments (black contours, interval: 10m•s -1 ): a) DJF LW-LMDz, b) JJA LW-LMDz, c) DJF SW-1L-LMDz, d) JJA SW-1L-LMDz, e) DJF SW-ML-LMDz, and f) JJA SW-ML-LMDz. The color shading corresponds to the differences between LW-LMDz minus MERRA in a) and b), SW-1L-LMDz minus LW-LMDz in c), and d), and SW-ML-LMDz minus LW-LMDz in e) and f). Units are in m•s -1

Acknowledgments

This work was supported by the Laboratoire de Recherche Conventionn Yves Rocard, a collaborative unit between CEA and Ecole Normale Supérieure. BR acknowledges funding from the European project ARISE 2 (Horizon 2020, GAN653980), FL from VESRI Schmidt Future project "DataWave", and AC from the Santander-UCM grant StratWars (PR87/19-22679).