Downslope windstorms in the Isthmus of Tehuantepec during Tehuantepecer events: a numerical study with WRF high-resolution simulations

Tehuantepecers or Tehuanos are extreme winds produced in the Isthmus of Tehuantepec, blowing south through the Chivela pass, the mountain gap across the isthmus, from the Gulf of Mexico into the Pacific Ocean. They are the result of the complex interaction between large-scale meteorological conditions and local orographic forcings around Chivela pass and occur mainly in winter months, due to cold air damming in the wake of cold fronts that reach as far south. Even though the name refers 5 mostly to the intense mountain gap outflow, Tehuantepecer episodes can also generate other localized extreme wind conditions across the region, such as downslope windstorms and hydraulic jumps, strong turbulent flows that have a direct effect on the Pacific side of the isthmus and the Gulf of Tehuantepec. This study focuses on investigating these phenomena using high horizontal and vertical resolution WRF (Weather Research and Forecasting) model simulations. In particular, we employ a four nested grid configuration, with up to 444 m horizontal 10 spacing in the innermost domain and 70 hybrid-sigma vertical levels, 8 of which lie within the first 200 m above ground. We select one 36 hour period in December 2013, when favorable conditions for a strong gap wind situation were observed. The high-resolution WRF experiment reveals a significant fine-scale structure in the strong Tehuano wind flow, beyond the well known gap jet. Depending on the Froude number upstream the topographic barrier, different downslope windstorm conditions and hydraulic jumps with rotor circulations develop simultaneously at different locations east of Chivela pass with varied 15 crest height. A comparison with observations suggests that the model accurately represents the spatially heterogeneous intense downslope windstorm and the formation of mountain wave clouds for several hours, with low errors in wind speed, wind direction, and temperature. Copyright statement. TEXT


Introduction
The Isthmus of Tehuantepec, in the Mexican state of Oaxaca, is the narrowest stretch of land separating the Gulf of Mexico from the Pacific Ocean. The Sierra Madre mountains cross the isthmus from east to west, leaving, however, a pronounced gap in the middle (Chivela Pass), coinciding with the point of shortest distance between the two sea masses, of only 200 km. The elevation of the Chivela pass is 224 m, whereas mountain peaks in the side sierras reach 2000 m, creating ideal conditions to 5 generate a powerful wind corridor (Romero-Centeno et al., 2003). In winter, cold high-pressure systems originating in North America move over the Gulf of Mexico in the wake of south-reaching cold fronts, and large pressure differences develop across the isthmus between the bay of Campeche and the Gulf of Tehuantepec, on the Pacific side. This pressure gradient results in a northerly wind situation, in which the flow is accelerated southward by cold air damming, traveling through Chivela Pass to finally blow violently outward into the Pacific Ocean, as it fans out and curves anticyclonically. This gap outflow extends for thus, the gap wind does not encounter any ambient large-scale flow to merge with, describing a very-close-to-inertial trajectory with a large radius, also as a consequence of the diminished Coriolis force (Steenburgh et al., 1998). In the Gulf of Tehuantepec, the strong sea surface wind stress generates intense upwelling and vertical mixing in the upper ocean (Hong et al., 2018). The 15 largest number of these events tend to occur in December, with a mean duration of 48h (Brennan et al., 2010). These powerful mountain gap winds are called Tehuantepecers or Tehuanos, and have been the focus of several previous studies (Brennan et al., 2010;McCreary et al., 1989;Jaramillo and Borja, 2004) detailing their general setting, drivers and dynamics (Steenburgh et al., 1998) of strong wind situations in the isthmus. Little knowledge exists, however, on the fine-scale structure of the Tehuantepecer flow elsewhere across the Isthmus, away from Chivela pass. There is evidence from observations and earlier numerical studies 20 (Steenburgh et al., 1998) that the low elevation topography of Chivela pass can excite mountain waves, and also that these are not only restricted to the pass itself but extend into the much higher mountain crests to the west and especially to the east, as the cold air pool is often thick enough to surpass them. This gravity wave activity can potentially result in downslope windstorms (DSWS hereafter) producing severe turbulent phenomena such as rotors and hydraulic jumps (HJ hereafter) (Durran, 1986a;Sheridan and Vosper, 2006) on the Pacific side of the isthmus, to the lee of local orography, which would explain several 25 accidents related to strong winds reported by the Oaxaca´s Civil Protection Commission (Santiago, 2018;Hernández, 2018;Televisa, 2018;Rodríguez, 2018) every year during some Tehuano occurrences. The ability to understand and forecast these events is very relevant, since the Isthmus has been an important development site for wind farms since the 2000's (Coldwell et al., 2017). Currently this region allocates 76.8% of the wind power capacity installed in Mexico, with approximately 2360 MW (Baxter et al., 2017), which is expected to double to 5076 MW by 2020 (AMDEE ,2018).

30
The main goal of the present work is to study the variability of flow behavior, beyond the well-known mountain gap jet, in Tehuano wind episodes across the isthmus of Tehuantepec, depending on topographic barrier height and thermodynamic conditions of the air mass, using high-resolution simulations with the Weather Research and Forecast (WRF) model. Many studies have successfully employed WRF to analyze this kind of mountain-flow events in other parts of the world. In the US, for example, B. Pokharel et al. Pokharel et al. (2017b) study a DSWS and HJ to the lee of the Medicine Bow Mountains in southeast Wyoming. Another DSWS in this same area is also investigated by Grubišić et al. (2015) with WRF. Other studies such as Cao and Fovell (2016), Pokharel et al. (2017a), Prtenjak and Belusic (2009), Jung-Hoon and Chung (2006), Ágústsson and Ólafsson (2014), Priestley et al. (2017 also use WRF at high resolution to analyze mountain wave flows in different locations. To the best of our knowledge, there is not any previous work that studies in detail lee wave phenomena in the 5 Mexican state of Oaxaca.
The high-resolution WRF simulations employed in the present study allow us to obtain a more complete knowledge of these events, from the synoptic scale to the small scale, focusing on the downslope winds and the HJ that develop along the mountain ranges neighboring Chivela pass. The article is organized as follows: in section II the methodology is explained in detail, from the climatology of the region to the model configuration. In section III, the primary results obtained are shown, divided into 10 synoptic-mesoscale situation and upstream-downstream structure of the phenomena. Finally, in section IV the conclusions reached are discussed.

Mountain wave phenomena, hydraulic analog and relation with the Froude number
DSWSs result from the intense flow acceleration occurring on the lee slope of a mountain under certain cross-barrier wind conditions. They resemble the behavior of water flow in an open channel when encountering an obstacle, such as when a 15 relatively slow river increases its speed when flowing in a thin layer over a mill's dam or over a rock. As in the case of the river, DSWSs often end abruptly with a return to the state upstream of the obstacle through a turbulent HJ somewhere downstream.
The similarity in both fluid's behavior suggests that the physical processes behind are also alike, and that shallow-water theory could be applied somehow to the atmosphere (Long, 1953). However, the complexity of the unbounded atmospheric flow, without a free surface, makes it difficult to make the analog so simple, because gravity waves in the atmosphere propagate 20 vertically in addition to horizontally as in shallow water, and non-linear effects are important.
Considerable observational and numerical experiments have been developed to elucidate the dynamics of atmospheric mountain lee flows (see Durran (1990) for a review). The current view is that in the atmosphere, DSWSs are observed when stable air at low levels flows, similarly to shallow water, as like having a free surface somewhere above the obstacle preventing vertical energy dissipation. This occurs when gravity waves do not propagate vertically and break because of the presence of a critical 25 level  or due to overturning for having too large amplitude (Clark and Peltier, 1977), or when without wave breaking, there is an interface separating highly stable lower layers from less stable air above (Durran, 1986b;Klemp and Durran, 1987;Bacmeister and Pierrehumbert, 1988). Wave breaking creates a well-mixed layer to the lee of the obstacle, which generates a dividing streamline separating undisturbed flow aloft and trapped energy and flow analogous to hydraulics in the lower surficial branch (Smith, 1985b). A highly stable lower layer topped by less stable air results in reflecting and 30 decaying waves aloft, enhanced non-linear effects, and the atmosphere also flowing like a two-layer fluid and behaving like shallow water when encountering the barrier.
In shallow water theory, the Froude number (Fr), which is the ratio of the mean speed to the intrinsic gravity wave phase speed in the fluid, determines its behavior when encountering the obstacle, depending on whether gravity is balanced mostly by acceleration (Fr>1) or pressure gradient forces (Fr<1). HJs and significant flow acceleration to the lee of the obstacle occur when the fluid transitions from a subcritical (Fr <1) to supercritical (Fr>1) regime at the top of the barrier. It is not straightforward to define a Froude number to determine the regime of atmospheric flows (Sheridan and Vosper, 2006;Smith, 1989Smith, , 1985a as it is for shallow water, because there is no clear analog to flow depth controlling gravity wave speed and pressure perturbations. Key factors determining wave properties and flow characteristics in this case are stratification, wind 5 speed and barrier height, or rather, the relative values among them. A quantity that combines these three variables and is referred to as Froude number for shallow atmospheric flows over mountains by many authors, is the following (Smith, 1989): where U is the flow speed, N , the Brunt-Väisäla frequency (Equation (2)), and H, the mountain height.
with g the acceleration of gravity, dθ/dt, the potential temperature gradient in the stable layer, and θ 0 , the potential temperature at the base of this layer. To avoid confusion with the classical Froude number, the inverse of Fr defined as above is often used instead and called the non-dimensional mountain height. Fr is a measure of flow deceleration and stagnation upwind from the mountain (Baines, 1987). It can also be regarded as an estimate of nonlinearity. When Fr«1 there are significant nonlinear effects and blocking, whereas for Fr»1 the opposite occurs (Smolarkiewicz and Rotunno, 1989). Fr around 1 indicates a 15 transitional regime between the two states and favorable conditions for the formation of DSWS and HJs.

WRF Configuration
We use the Advanced Research WRF (ARW) model (Skamarock et al., 2008) version 3.9 (WRFV3.9) to perform the simulations. Based on a fully compressible and non-hydrostatic dynamic core, WRFV3.9 is a limited-area mesoscale and microscale 20 model, with a terrain-following hydrostatic-pressure vertical coordinate, designed for operational forecasting, as well as research. For the experiments, we employ a nested domain configuration, in order to achieve sufficiently high resolution in the innermost grids to capture the small scale structure of the flow, while reproducing the synoptic phenomenology conducive to local DSWS in the parent one ( Figure 1a).
D03 covers the whole isthmus area, with Chivela pass approximately in its center, while the highest resolution domains d04 We simulate a 36-hour period, from 2013-12-23 12:00 to 2013-12-25 00:00 UTC, which registers DSWS conditions in the observational data. Regarding the main physics options, the simulations use the tropical suite configuration (Table 1), introduced in WRF version 3.9., except for the planetary boundary layer, which is parametrized by the Shin-Hong scale-aware scheme (S-H) (Shin and Hong, 2015). The next table summarizes this physics configuration used.
The S-H planetary boundary layer option is more suitable for the high resolution of the innermost domain (444 m) because it helps to mitigate a double counting effect of the small-scale processes in gray-zone resolutions. Apart from this, this scheme provides a turbulent kinetic energy (TKE) diagnostic variable useful for our analyses.

Global model and real data
Global Forecast System (GFS) analysis data from the National Centers for Environmental Prediction (NCEP) is used as initial and boundary conditions for the WRF model, with a 3-h update interval. The horizontal resolution of this dataset for all variables is 0.5 • x 0.5 • , with 32 vertical levels ranging from 1000 to 10 hPa. The observational data used in this work is provided by the Mexican National Laboratory of remote sensors (https://clima.inifap.gob.mx/lnmysr/Estaciones/MapaEstaciones), col-10 lected every 15 minutes at two meteorological stations whose location is presented in Table 2  Model data are extrapolated from their native sigma levels to the height of the meteorological station using equation (3), which relates wind speed with friction wind speed, which is also a diagnostic variable in the model.
Where U * in the friction velocity, K is the von Kármán constant, z is the height and z 0 is the surface roughness (m).
3 Results and discussion   results suggest that they also occur in the mountains west and especially east of Chivela pass. The potential temperature field on the same sigma=3 level at the enhanced resolution (4km) of the next nested grid D03 (Figure 3c) illustrates how the stable cold air mass to the north surmounts the lower hills neighboring the pass, particularly those to the east where elevation increases more gradually. As shown in the following section 3.2, acceleration on the top of these mountains and to their lee is related to gravity wave activity, which results in the development of strong downslope winds, rotors and HJs. The model sounding at location NP on the Gulf of Mexico side of the Isthmus (Figure 3d), shows in the temperature profile a stable lower boundary 5 layer capped by a very stable isothermal layer from 850hPa up to about 2000 to 2500m, defining the depth of the cold air pool, indeed above the aforementioned mountain tops. A weak inversion associated with the subsidence within the high pressure system aloft is also clearly apparent just below 500hPa. Above this level, winds are weak and veer from being southeasterly to southwesterly in the upper troposphere. Below 500 hPa, winds back from an easterly to a northeasterly direction at about 800hPa, and more strongly in lower levels, becoming westerly at the surface, indicating intense cold air advection. There is a 10 pronounced reverse wind shear in the lower troposphere.
As mentioned earlier, the only available observations in the area are from stations MET1 and MET2, whose positions are marked in Figure 3c. occurrence. At about 0 UTC on the 29 th , the episode decays and winds go back to local breeze regimes. The extent of the simulated period covers the first 36h of the event (shaded in Fig.4), corresponding with its highest intensity in both locations.
These observations away from Chivela pass show evidence, as the simulations suggest, that the neighboring mountains may 20 also induce strong winds, even though they certainly do not reach out as far as the mountain gap outflow does.

Upstream-Downstream structure
In this section, we focus on the fine scale structure of the flow acceleration across the Isthmus, and more precisely, on that occurring on the westernmost section of the Sierra Madre de Chiapas mountains bordering Chivela Pass, apart from the wellknown strong gap wind jet. This is the area covered by the d04 domain of 1.3 km resolution (highlighted in green in Fig. 3c), encompassing with 137 km from north to south the wind flow path before and after crossing the mountains. Figure 5 depicts 5 from d04, latitudinal cross-sections (red lines in Figure 1b) facing east (south to the left, north to the right) of different variables at the longitude of the two validation points MET1 (Fig. 5a-f) and MET2 (Fig. 5g-l).
The tight isentropes on the windward side of the mountains (to the right in Fig. 5a, d, g and j), indicate a quite similar stable stratification in the whole lower tropospheric column in both cases; stronger in the layers below about 2.5 km and weaker above. The temperature profile in Fig. 3d evidences that the 2500m height level corresponds to the depth of the cold air. Winds 10 are northerly in these lower layers to the north of the mountains, with somewhat higher speeds above 20m/s upstream from MET1 than further east, north of MET2. In both cases, winds above the cool pool are much weaker, and back to an easterly component at about 4000m and above. The latter height is therefore a critical level, since the cross-mountain flow component becomes null; thus, triggered gravity waves will break and dissipate when they approach it. Markowski and Richardson (2010) outline seven conditions conducive to DSWS, albeit not all of them absolutely necessary. These conditions are: (1) asymmetric 15 mountain with steeper lee than windward side (Fig 5), (2) crossed by strong winds (> 15m/s) (Fig 5), (3) mostly normal to the barrier (Fig.3, 5). (4) A stable layer above the top and less stable above that (Fig 3d, 6c), (5) with cold air advection and large-scale subsidence to maintain the stability (Fig 3d, 6c). Apart from this, (6) reverse wind shear above (Fig 3d, 6c) and (7) no cool pool in the lee (Fig 3a, 6b), is also desirable. These conditions are all perfectly met for both locations analyzed, as discussed previously, and indeed intense DSWSs occur in both cases.

20
The stably stratified cross barrier flow displays wave activity from early on, and wave breaking, as the vertical isentropes suggest ( Fig. 5a and d), enhances turbulent mixing ( Fig. 5c and f) and yields a region of weak stability and reverse flow immediately downwind from the mountain crests (Fig 5g and j). In both cases, isentropes on the windward side sink sharply under these layers of low stability on the lee side, much more pronouncedly for the tallest mountain (Fig. 5j). Encompassing the well mixed region to the lee, a split streamline develops (Smith, 1985b), and below its lower branch there is flow thinning 25 and a significant increase in wind speed. The particular features existent on the lee side differ, however, depending of the height of the topographic obstacle. The strong accelerated flow confined by layers of strong wind shear and turbulence extends for many kilometers downwind from the lowest mountain (Fig 5i), while ending with a HJ and rotors at the foot when the barrier is higher (Fig 5j). The formation of either of these lee wave events is related to the Froude number upstream (Smith, 1989).
We note that only the turbulence produced by the existent wind shear in the column is well captured in the simulation (Fig.   30 5i and l). The subgrid scale turbulence associated with gravity waves due to rotors and non-local turbulent advection or with instability resulting from wave overturning is not represented and accounted for, since a much higher horizontal and vertical resolution of the order of tens of meters would be needed in order to explicitly resolve these features (Vosper et al. 2018).  during nighttime and even higher at some other times during the day, indicating supercritical conditions. The Brunt-Väisälä frequency around mountain top is between 0.020 and 0.025 s −1 . The generated mountain waves have relatively short wavelength and small amplitude and a modest HJ that propagates downstream can be seen at the initial stages of the episode at 15 5 UTC December 23 (Fig. 5a). 12 hours later, during nighttime, the aforementioned strong jet extending for tens of kilometers downwind is fully formed, with values above 35 m/s at about 750m above ground and strong turbulence at the surface and in the layers above the jet, where stability is much reduced. The temperature profile upwind (at 3 UTC December 24 and location UP1 in Fig 5) shows in the temperature profile a stable lower boundary layer capped by a very stable isothermal layer from 850hPa up to about 800hPa, or 2500m, defining the depth of the cold air pool. At the observation location MET1, downwind 5 from the mountain, the lowest layers up to about 900 m are well mixed (slightly lower, up to 750m further downstream), the result of the intense surface turbulence. Above the latter height and up to about 1500 m elevation, a strongly stable layer exists corresponding with the aforementioned packing of the isentropes (and streamlines). This is where the highest wind speeds are found. Stability is sharply reduced aloft, in the region encompassed by the dividing streamline, and winds are rather weak and have an easterly component for the most part, parallel to the mountains. from the jump to the location of observation station MET2 is also evident in Fig. 5l. The temperature profile upstream is very similar to that of the MET1 case, but downwind from the mountain at location MET2, lacks the strongly stratified layer present in the case of the lower mountain. We note that the soundings in Figure 6 show complex temperature profiles with 20 discontinuous stratification in the planetary boundary layer, which underscores the importance of using high vertical resolution in simulations for these kind of studies.
Results from the innermost nested grid d05 with the finest resolution (Fig 7) suggest that trapped lee waves develop in the MET1 case within the high stability layer where the strongest winds are found, just below the low stability region and the critical level aloft that prevent their vertical propagation. This wavelike pattern is a common feature in DSWS periods (Pokharel 25 et al., 2017b;Hertenstein and Kuettner, 2005) and fully formed 12 h later (Fig. 7c), extends for more than 100 km downstream aligned with the general orientation in the northwest-southeast direction of the Sierra in the region. Waves are absent further east in the MET2 cross section, where the topographic barrier is higher and a stationary HJ forms instead, as discussed above.
Our results suggest that during Tehuano wind events, the Pacific side of the Isthmus of Tehuantepec east of Chivela pass is very prone to host extreme winds' phenomena. The formation of DSWSs in the area increases the impacts of the already strong 30 mountain gap winds.

Validation
Finally, we contrast our simulation results with the very few data available for validation at meteorological stations MET1 and MET2. The two plots in Figure 8 compare the simulated wind speed (in d04 and d05) with observations from both stations.  Table 2). The similarity in the low mean error between both simulations and their high correlation throughout the period are due to the nature of the event in that area, an intense and mostly steady jet that the d04 domain resolution (1.3 km) is already capable of resolving accurately. However, results in MET2, downwind from the strong HJ, present more differences between d04 and d05 ( Figure 8b) and there is a 5 significant improvement in d05 with respect to its parent domain d04. Wind speeds in d04 are overestimated (Mean error ME = 2.76 m/s), and present a daily cycle that is absent or very subtle in the observations. The complexity and fast variability of HJs formation in this area are better resolved in the higher resolution grid, which perhaps reproduces more accurately the stagnant flow and rotor formations downstream from the HJ. With regard to wind direction, errors are small for MET1 and significantly higher for MET2, due to the same reasons. As in the case of wind speeds, wind direction results from the finer grid d05 are also better than in d04. Temperature errors are equal or below 1 K in both locations and domains, hence the surface thermal evolution is well captured.   Table 3. Wind speed (WS), wind direction (WD), and temperature (T) mean errors (ME) and mean absolute errors (MAE) at MET1 and MET2 locations during the simulated period and for the two higher resolution domains (d04 and d05).
Lee waves can promote orographic cloud formation at different scales, depending on the amplitude of the wave and elevation (Armi and Mayr, 2011;Szmyd, 2016). Model results in d05 suggest that lenticular clouds form at the crests of the trapped lee waves depicted in Figure 7.

Conclusions
In the present work, we studied lee wave phenomena occurring during Tehuano events on the Pacific side of the Isthmus of Tehuantepec using WRF high-resolution simulations. Orographic forcings at different scales result in the well-known gap wind jet off Chivela pass, but also in downslope windstorms and hydraulic jumps in the neighboring mountains. We analyzed these phenomena in an episode in December 2013 having the typical genesis of Tehuantepecer wind events. An Arctic air mass in 5 North America pushed as far south as the bay of Campeche due to cold air damming east of the Rockies continuing to the east of the Sierra Madre Oriental range in Mexico. The displacement of the associated high-pressure system on the wake of the cold front created large pressure differences across the Isthmus of Tehuantepec, ultimately producing the strong mountain gap winds through the low elevation of Chivela pass. The model simulates these intense winds, blowing with speeds at the surface of more than 25 m/s that extend for many kilometers from the mountain gap, fanning out well into the Gulf of Tehuantepec, as 10 it is commonly observed in Tehuano events (Steenburgh et al., 1998). The depth of the cold surge on the coastal plains of the Gulf of Mexico side of the Isthmus is about 2500m, therefore thick enough to surmount the lower elevations to the west, and especially to the east of Chivela pass. The flow over these mountains results in intense DSWSs to their lee, with the generation of intense turbulence, HJs and rotors, depending on the particular height of the topography. We focus on two locations of different barrier elevation from where there are surface observations 15 downwind: one of 963m closer to Chivela pass and another further east, with height increasing to 1736m The thermodynamic characteristics of the air mass are rather uniform upwind of both mountains, with strong stability within the cool pool and weaker above, and intense northerly winds that back and weaken aloft to a more easterly component parallel to the barrier. The critical level where the cross-mountain wind component becomes zero, inhibiting wave propagation, is about 4000m. Mountain waves are generated in both cases, with smaller amplitudes for the lower mountain, where the Brunt-Väisälä about half. Wave breaking produces mixing and generates a region of low stability to the lee of the mountains, which is deeper where the waves have higher amplitude. The Froude number is around 2.5 at crest height in the lower barrier and the flow presents a supercritical behavior. The region of low stability to the lee of the mountain lies above about 1500m and leads to a packing of the isentropes and streamlines underneath, resulting in strong stability and flow acceleration. An intense jet develops with wind speeds of 35 m/s at about 750m above ground extending for tens of kilometers downwind from the mountains. Wind 5 speeds are reduced closer to the surface due to intense turbulence. Trapped lee waves form at about 1500m, just below the well mixed layer aloft that prevents their vertical propagation. The Froude number decreases to about 1 further east as elevation rises and the flow presents a critical regime. Isentropes on the windward side of the mountain sink much more pronouncedly under the wider mixed layer generated by wave breaking to the lee, and are tightly packed in a shallow layer above the surface.
This generates an intense wind storm on the lee slope of the mountain, with surface wind speeds up to 35 m/s. The accelerated 10 flow down the mountain ends abruptly at its foot, where the flow turns to subcritical state and a marked stationary HJ forms, with vertical velocities of 6 m/s. A rotor circulation develops further downstream from the jump.
Only limited observations are available to validate our model results. Errors in surface wind speeds, directions and temperature are small at the only two stations available on the Pacific coastal plain downwind from the mountains. In addition, lenticular clouds similar in location and pattern to those produced by the model, are apparent in satellite imagery of the day of 15 the event, providing valuable indication that the mountain lee wave phenomena simulated indeed correspond to a real scenario.
Our model results suggest that when the cold air mass intruding from the north on the Gulf of Mexico side of the Isthmus is thick enough to surmount the mountains, extreme wind events develop in the area during Tehuano events beyond the gap wind jet. These include DSWSs and HJs, which are intense and highly turbulent flows that can have a substantial impact on the existent wind farm industry in the region.
was performed at the Centro de Supercomputacion de Galicia (CESGA) (http://www.cesga.es/). Their computer facilities and support have been indispensable to carry out this project. Finally, we would like to acknowledge the Non-Lineal Physics Group of Universidade de Santiago de Compostela (http://www.usc.es/en/investigacion/grupos/gfnl), which is where this whole project has been developed.