Dynamical downscaling of historical climate over CORDEX Central America domain with a regionally coupled atmosphere–ocean model

The climate in Mexico and Central America is influenced by the Pacific and the Atlantic oceanic basins and atmospheric conditions over continental North and South America. These factors and important ocean–atmosphere coupled processes make the region’s climate a great challenge for global and regional climate modeling. We explore the benefits that coupled regional climate models may introduce in the representation of the regional climate with a set of coupled and uncoupled simulations forced by reanalysis and global model data. Uncoupled simulations tend to stay close to the large-scale patterns of the driving fields, particularly over the ocean, while over land they are modified by the regional atmospheric model physics and the improved orography representation. The regional coupled model adds to the reanalysis forcing the air–sea interaction, which is also better resolved than in the global model. Simulated fields are modified over the ocean, improving the representation of the key regional structures such as the Intertropical Convergence Zone and the Caribbean Low Level Jet. Higher resolution leads to improvements over land and in regions of intense air–sea interaction, e.g., off the coast of California. The coupled downscaling improves the representation of the Mid Summer Drought and the meridional rainfall distribution in southernmost Central America. Over the regions of humid climate, the coupling corrects the wet bias of the uncoupled runs and alleviates the dry bias of the driving model, yielding a rainfall seasonal cycle similar to that in the reanalysis-driven experiments.

interaction from the Pacific Ocean, the Gulf of Mexico and the Caribbean Sea. The regional weather and climate, largely driven by the effect of the vast surrounding water bodies, results from the junction of various processes at different time and spatial scales. Rainfall is featured by a contrasting Pacific-Caribbean distribution characterized by the Mid Summer Drought (MSD, Magaña et al. 1999) on the Pacific slope. The Inter-Tropical Convergence Zone (ITCZ), deep convection, low level moisture transport, cyclone activity and mid-latitude air intrusions are pointed as key drivers for regional rainfall (Amador et al. 2016). Interannual variability is mainly modulated by ENSO, which is known to be linked with drought-prone conditions as demonstrated by the analysis of observational data (Sánchez-Murillo et al. 2017). Meanwhile, the regional impact of other low frequency modes such as PDO (Pacific Decadal Oscillation) and AMO (Atlantic Multidecadal Oscillation) is related with the modulation of dry and wet extremes depending on the phasing of the signals. This behavior has been detected from paleorecords for northern Mesoamerica (Park et al. 2017) in agreement with findings of the PDO and AMO modulating rainfall over Mexico and Central America (Méndez and Magaña 2010) and the PDO-ENSO relationship over Mexico Pavia et al. 2006). The region is highly vulnerable to climatic fluctuations (Bárcena et al. 2010;Ahmed et al. 2009) and a well known climate change "hot spot" (Giorgi 2006). In fact, the latter work identified this region as the most responsive tropical region to climate change. This area is well known for its extreme weather, including major droughts (Stahle et al. 2011), impact of hurricanes (Amador et al. 2010) and increasing temperature trends (Aguilar et al. 2005).
To date, different dynamic downscaling modeling frameworks have been applied to the analysis of the climate in the Mesoamerican domain for Central America. For example, Rivera and Amador (2009) showed a general overestimation of rainfall with the MM5 model using ECHAM4.5 boundary conditions and 90 and 30 km resolutions, and suggested further studies to evaluate the sensitivity to boundary layer and radiation schemes. Using PRECIS RCM, Karmalkar et al. (2011) analyzed two sets of experiments with 25 km grid resolution, and found a reasonable representation of the key regional climate features for recent climate conditions, while future projections suggested a drying pattern in most of Central America. Diro et al. (2012) used RegCM4 to conduct an analysis over the Central America CORDEX domain at a 50 km grid resolution to evaluate the sensitivity of rainfall to land and sea surface schemes. Their results conclude the simulations capture fairly well the rainfall distribution. However, their study identified an underestimation and misplacement in the location of the CLLJ. Moreover, systematic biases were found at sub-diurnal scales, with little differences across the land schemes. The sensitivity to resolution, convection schemes, and surface flux parameterizations in the CORDEX Central America domain was investigated by Fuentes-Franco et al. (2016). Using the Reg-CM4 model, they showed the relevance of resolving coastal topography for a better simulation of tropical cyclone activity in the region.
Considering the well known Pacific-Caribbean precipitation seesaw, a general drying trend for the region requires a more thorough analysis. Most modeling studies highlight that despite the general climate features been captured, the representation of the Mid Summer Drought (MSD, Magaña et al. 1999), the Caribbean Low Level Jet (CLLJ) structure and intensity as well as its related moisture transport , cold surges (Schultz et al. 1998) and tropical cyclones dynamics (Knutson et al. 2010) need to be improved. A common but important caveat is that most modeling studies performed in this region do not include ocean-atmosphere feedbacks at the regional scale. Typically, the regional atmospheric models are forced at their lower boundaries by prescribed sea surface temperature (SST), so that simulations cannot account for thermally ocean-driven intraseasonal interactions. Meanwhile, other studies based on GCMs projections, highlight the need for higher resolution and a better representation of ocean interaction processes (Nakaegawa et al. 2014). Recently, Sitz et al. (2017) have shown the usefulness of coupled regional models for a good representation of precipitation and SST in different CORDEX domains, among them the Central American domain.
In general terms, most of the added value of regional models over land areas in comparison to global coupled models (GCM) simulations can be traced back to a better representation of the orography (e. g. Di Luca et al. 2015;Xu et al. 2018;Teichmann et al. 2013). However, coupled regional models show a clear improvement with respect to both uncoupled atmospheric simulation and the driving global model in regions of intense air-sea interactions or of large-scale feedback (e.g., Cabos et al. 2017). In the last few years, a number of regional coupled atmosphere-ocean models have been used to explore a wide range of inter-annual and intra-seasonal processes in present climate simulations over different regions of the world (e.g. Zou et al. 2016;Ratman et al. 2009;Sein et al. 2014;Cabos et al. 2017), including the Mesoamerica region .
Another critical factor to improve the representation of the climate of the region in Regional Climate Models (RCM) simulations is the extension of the domain (Sein et al. 2014). For this reason, in our model domain the northern boundaries of the Central American domain defined in the framework of the Coordinated Regional Downscaling Experiment (CORDEX) project have been shifted away from the Mexican territory (similar to Fuentes-Franco et al. 2014). The domain targets to simulate adequately some of the key large-scale circulation features that affect regional climate. In particular, the extension of the CORDEX Central American domain allows for a better representation of the ITCZ movements (Schneider et al. 2014), the CLLJ dynamics (Amador 2008), the seasonal changes of the WHWP (Western Hemisphere Warm Pool, Wang and Enfield 2001), and the entrance of cold fronts as well as the large-scale circulation interactions that account for the regional rainfall climatology (see e.g., Englehart and Douglas 2001). Because of the difficulty in representing well all these processes, it is important to explore how higher horizontal resolution and ocean-atmosphere coupling may improve the representation of the climate in the region.
Here we perform a set of both ocean-atmosphere coupled and atmosphere-only simulations. The model systems considered are the atmosphere-only RCM REMO and the regionally coupled ocean-atmosphere model ROM (REMO-OASIS-MPIOM), which consists of a regional atmosphere coupled to a global ocean. The simulations are forced by 32 years of the ERA-Interim reanalysis and a 55-year integration of present-time climate with the Max-Planck Institute (MPI) Earth System Model (ESM), MPI-ESM. The atmospheric domain is common for ROM and REMO and extends the Central American CORDEX domain to include all the Mexican territory and is run in two different horizontal resolutions (~ 25 and ~ 50 km).
The objectives of this study can be summarized as follows: 1. Examine the added value that high-resolution REMO and ROM bring with respect to the driving global model in the area of study. 2. Examine the role of air-sea coupling processes in simulating the present climate by comparing coupled and uncoupled runs. 3. Asses the skills of both ROM and REMO in reproducing the observed regional climate over our extended COR-DEX Central American domain.
The remainder of the paper is organized as follows. In Sect. 2 we present the system design along with the data and methods used for the model validation. In Sect. 3 we show results from the climatology that the model system is able to capture, examining the impact that higher resolution and ocean-atmospheric coupling have on this ability. Finally, in Sect. 4 we discuss the improvements and limitations of coupled and non-coupled simulations and how this may impact the simulations of future climate change scenarios.

The modeling frameworks
In this study, we use the coupled system ROM (see Sein et al. 2015) and its atmospheric component, the regional climate model REMO (Jacob et al. 2001). REMO is a hydrostatic, three-dimensional regional climate atmospheric model that uses the physical package of the global circulation model ECHAM4 (Roeckner et al. 1996). In the vertical, variations of the prognostic variables are represented by a hybrid vertical coordinate system. In ROM, REMO is coupled to the global ocean model, MPIOM, that includes sea-ice and marine biogeochemistry modules (cf. Marsland et al. 2003). In addition, ROM includes a global hydrological discharge model (HD) that computes river runoff at 0.5° horizontal grid resolution and is coupled to the atmospheric and ocean components. The coupling between the atmosphere and the ocean is done through the OASIS coupler (Valcke et al. 2013), which we configured with a coupling period of 3-h model time. This coupling frequency is necessary to resolve well the diurnal cycle. In this work we switch off the marine biogeochemistry module and refine the computational grid of MPIOM to increase the horizontal resolution up to 10 km in the Gulf of Mexico and the Tropical North East Atlantic region (Fig. 1a). The refinement is achieved by placing the coordinate grid poles over North America and North West Africa. In the vertical the MPIOM model uses a z-type coordinate that in our configuration is discretized on 40 levels.
In REMO standalone, both surface and boundary conditions are taken from ERA-Interim or MPI-ESM data. In ROM, the surface fluxes (i.e., heat, momentum and freshwater fluxes) forcing the oceanic component outside the domain of the atmospheric model are taken from 6-h reanalysis or GCM data. The atmospheric component of ROM is forced laterally by 6-h fields taken from the same global dataset used for the oceanic component. The extent of the atmospheric domain is akin to the domains used in Central America CORDEX simulations, but extended northward so that the Mexican territory is located away from the domain boundaries (Diro et al. 2012;Fuentes-Franco et al. 2014). The use of a global ocean model prevents the complications associated with prescribing lateral boundary conditions in regional ocean models. Specifically, we avoid inconsistencies in spatial and temporal resolutions between the regional model solution and the lateral driving fields. Furthermore, a global ocean model allows coastal trapped waves originated outside the coupled domain to influence the barotropic sea level variability and the bottom pressure in the coupled domain (Sein et al. 2015).

Experiment design
In this work we compare observational datasets, a fully coupled earth system model simulation (MPI-ESM), and downscaled simulations with the regionally coupled model, ROM, and its standalone atmospheric component, REMO. We use ERA-interim and MPI-ESM as boundary condition to perform two sets of experiments with ROM and REMO.
In a first set of experiments, lateral atmospheric boundary conditions are taken from ERA-Interim (Dee et al. 2011) for both the coupled and stand-alone simulations, while surface boundary condition are generated by MPIOM in the coupled ROM simulations or taken from ERA-Interim in the standalone REMO simulations. These simulations cover the period 1980-2012. The second set of experiments, covering the period 1950-2005, is similar to the first set but uses CMIP5 MPI-ESM twentieth century simulations as boundary conditions. For the first (second) set of experiments, the spin-up of the coupled model, which is forced by ERA-Interim (MPI-ESM), is performed by running three consecutive times the period 1980-2012 (1950-2005), with each consecutive iteration starting from the final state of the previous run. Then, the final state of the third integration is used as initial condition for the final simulation.
In both sets of experiments, the atmospheric model uses a rotated grid with uniform horizontal resolution of about 50 km. However, to explore the impact of the horizontal resolution, we repeated the simulations of the first set using the same domain but with a horizontal grid of about 25 km. A summary of the experiments and naming convention is shown in Table 1.

Observational datasets
We use monthly mean values for the period 1980-2012 from two reanalysis datasets and various observational datasets, including time series of in-situ and satellite data. Sea level pressure (SLP), 2-m temperature (T2M), and both wind components at 925 hPa are taken from the commonly used ERA-Interim reanalysis (Dee et al. 2011), which provides fields at about 80 km resolution for the period 1979-2017. The sea surface height (SSH) and temperature (SST) are taken from the partially coupled Climate Forecast System Reanalysis (CFSR; Saha et al. 2010), which provides fields on a regular grid of about 100 km for the period 1979-2017. We choose the CFSR reanalysis because it is the closest dataset available to a coupled global reanalysis of the where coupling is computed and ocean/ice model grid (grid lines in black). b Orography of the Mesoamerican region (shading) and regions (rectangular boxes) used for the definition of the different indices and geographical points mentioned in the text. c Koppen-Geiger climate regions analyzed in this study (see Table 2 for details) atmosphere and the ocean. This coupled reanalysis provides fields that are physically more consistent than other atmosphere or ocean only datasets. In addition to the SST and SSH of CFSR, we use two more observational products: (1) the optimum interpolation SST (OISSTv2; Reynolds et al. 2002) analysis which is available for the period 1982-2014, and (2), the AVISO Altimetry SSH that is available for the period 1992-2010 and retrieved from https ://podaa c.jpl. nasa.gov/datas et/AVISO _L4_DYN_TOPO_1DEG_1MO). The simulated precipitation rate (P) is compared against the 12-month climatology of the TRMM 3B43 precipitation dataset (Huffman et al. 2007) and the Global Precipitation Climatology Project (GPCP) precipitation dataset (Kummerow et al. 2000). For both TRMM 3B43 and GPCP, the climatology is obtained as long-term mean of monthly values from January 1998 to December 2014. The TRMM 3B43 product combines ground rain-gauge measurements with microwave remote sensing to adjust infrared observations of geostationary satellites to build a 0.25° × 0.25° gridded dataset. The GPCP combines in-situ observations over land with satellite precipitation data to build a 2.5° × 2.5° global gridded datasets. The 925 hPa monthly wind climatology from IGRA radiosonde data located at Hato airport, Curaçao, is used to evaluate the CLLJ. To analyze the model skill to simulate T2M and P in region with a coherent seasonal cycle, we consider the Köppen-Geiger climate classification scheme (Kottek et al. 2006;Rubel et al. 2016), which aggregates the climate information in zones of similar T and P regimes. The climate type for each grid point within the Mesoamerican domain has been derived from the high-resolution 1986-2010 classification available at http:// koepp en-geige r.vu-wien.ac.at/prese nt.htm. It is based on the Climatic Research Unit (CRU) and the GPCC datasets. The features defining each Köppen-Geiger climate type can be found in Table 1 of Kottek et al. (2006). The climate types appearing in our study are listed in Table 2.

The atmospheric circulation climatology
The large-scale circulation is dominated by the North Atlantic Subtropical High (NASH). The NASH pattern has a strong seasonality that is characterized by a high-pressure system in the eastern Atlantic, which weakens during winter ( Fig. 2a) and strengthens during summer (Fig. 2b). The wintertime westward extension of the NASH (Fig. 2a) into North America has an important tropical connection, which is consistent with the sharp decrease of the Great Plains lowlevel jet (Weaver and Nigam 2008). This NASH behavior has been related with a thermal land forcing (Muñoz et al. 2008;Cook and Vizy 2010), and often occurs in conjunction with the intensification of the CLLJ in summertime. On the Pacific, the North Pacific High is stronger during summer. At higher latitudes, the Aleutian low deepens during winter ( Fig. 2a), fostering the intrusion of cold air and extratropical systems. The latter is of extreme importance as the associated frontal systems which provide a suitable environment for rainfall during the dry season over Central America (winter).
The regional atmospheric circulation is also strongly affected by the temperature seasonal cycle. Winter SSTs over the subtropics depict a pronounced meridional gradient ( Fig. 2g), which weakens in summer (Fig. 2h). The transition of the region of maximum SST from the Pacific (Fig. 2g) to the Caribbean and the Gulf of Mexico (Fig. 2h) clearly depicts the annual cycle of the WHWP. This transition occurs in association with the northernmost position of the ITCZ and brings the (summer) rainy season over the area (Fig. 2f). This migration of the ITCZ is connected with land processes that drive moisture from the Gulf of Mexico inland, contributing to the development of the socalled North America Monsoon. Furthermore, the rainfall intensification shows a marked contrast between the Pacific and Caribbean basins, remarking the horizontal gradient featured by the MSD. In contrast, rainfall is drastically reduced during winter, with the exception of the southern Central America and its Atlantic windward coast (Fig. 2e).

Model biases
We evaluate the climate of the extended CORDEX Central American domain as simulated with the global MPI-ESM and the regional climate model experiments. As ROM and REMO share the physics and domain, differences in the simulations driven by the same global model can be primarily attributed to the modified underlying SSTs due to the inclusion of the local air-sea coupling. With this aim, ERA-Interim and other observational datasets are used as a reference.  Figure 3 shows the SLP biases with respect to ERA-Interim for winter ( Fig. 3a-d, i-k) and summer ( Fig. 3e-h, l-n), for the global MPI-ESM, the regionally coupled ROM and the atmosphere-only REMO, respectively. Since the regional setups have a higher horizontal resolution (25 and 50 km), their data have been interpolated to the ERA-Interim resolution (~ 80 km), in order to prevent small-scale features from being spuriously interpreted as biases. The MPI-ESM  (Fig. 3a, e) shows a positive (negative) bias in winter (summer) over Mexico, the Caribbean and the Central America region, which is particularly pronounced over the mountainous regions of western US. Over the ocean, the positive bias in the Aleutian low in winter is noteworthy (cf. Fig. 3a), as well as the negative bias off Baja California (Fig. 3e). The latter is related to the atmospheric component of MPI-ESM (ECHAM), since it appears also in atmosphere-only simulations (see discussion of Figs. 4, 6). The sign of the large-scale bias pattern remains the very similar for both seasons and all the regional simulations. The coupled ROM simulations forced by ERA-Interim show the smallest SLP bias (Fig. 3d, h). While the biases over the ocean depend strongly on the coupling, their spatial pattern over land is rather sensitive to resolution and other details of the atmospheric component, and their magnitude is influenced by the global driving data. For the MPI-ESM forced simulations, the effect of the coupling is an overall increase in SLP over the ocean, modifying the circulation pattern of the driving model (Fig. 3b, f). The positive SLP bias over the eastern Pacific and the Caribbean is consistent with the strengthening of the NASH and colder SSTs (Fig. 6c, d), e.g. in line with Wang et al. (2008). Over land, an improved orography representation leads to a weaker SLP over the Sierra Madre, Rocky Mountains and Sierra Nevada, and stronger SLP over the lower regions of the Great Basin and east of the Rocky Mountains. This effect is clearer in the ERA-Interim-driven runs, especially in ROM25I (Fig. 3d, h). In this case, the surface pressure of the forcing can be considered as "perfect" and the biases over land seem to be associated to differences in orography (cf. also Sein et al. 2014). Therefore, the differences between ROM and ERA-Interim in this are can be primarily be attributed to the better representation of orography.

Sea level pressure
For uncoupled REMO simulations, the SLP biases over land have a spatial pattern similar to those of the coupled ROM, albeit with overall more negative values (Fig. 3i-n). The stronger negative biases are found in summer along the Mexican highlands (markedly in REMO50M), thus amplifying the MPI-ESM bias (Fig. 3l). The imposed SST introduces its own biases over the ocean, especially in summer (cp. Fig. 3i, m). For example, the MPI-ESM SST bias off California becomes stronger in REMO and seems to reinforce the bias over the Sierra Madre range. On the other hand, the ERA-Interim forced simulations (run with observed Hadley SST) do not show this bias. However, the differences between the driving model and the regional model are quite similar in spatial structure for both forcings. The negative SLP bias is mitigated in ROM25I with respect to ROM50I, and the coupling also seems helpful for the summer negative SLP bias over Mesoamerica, (cp. Fig. 3m, n with Fig. 3g, h), where ROM exhibits a better performance. Despite the land surface features (orography) being responsible for most of the detected bias, it is clear that including the air-sea coupling improves the representation of regional large-scale patterns.

Surface air temperature
A warm bias can be found over land for MPI-ESM (Fig. 4a, e), in both seasons, with values is below 2 °C in most of the region. This is clear over the central and coastal plains, the Yucatan peninsula, and most of the regions considered for South America. In winter, the biases are larger over southern Baja California and the Great Plains, while in summer they are stronger over South America and in northern Baja California. Conversely, the negative bias for both seasons are mainly found over northwestern Mexico and southwestern USA mountain ranges, particularly along the Sierra Nevada and the southern Rocky Mountains, and can probably be attributed to the different representation of the orography. Here, the dominance of low height stations in the assimilated data might play a role (cf. Karmalkar et al. 2011). Over the ocean, the strongest biases in winter can be found north of Cape Hatteras, and mainly arise from a late separation of the Gulf Stream. In summer, the largest biases off the coasts of California and Baja California are related to a poor representation of the California eastern boundary current.
When comparing the simulated MPI-ESM T2M over land areas (Fig. 4a, e) with the RCMs, both REMO ( Fig. 4i, l) and ROM (Fig. 4b, f) are generally colder over western North America and Mesoamerica. Over the ocean, the coupling has a cooling effect on T2M, reducing the positive bias that MPI-ESM exhibits along the coasts of Baja California and the Gulf Stream. Still, cold T2M biases develop over the Caribbean Sea and western North Atlantic over a colder SST, and in the Tropical Pacific due to a displaced ITCZ. It is noteworthy that REMO reproduces closely the driving model biases over the ocean. Generally, the ERA-Interim forced simulations show the smallest T2M bias (Fig. 4c,d,g,h,j,k,m,n), and higher resolution adds little to the largescale bias pattern. Unlike MPI-ESM, the downscaled T2M over land is generally warmer than the driving reanalysis, with the exception of the mountainous regions in western North America and Mexico, which is more noticeable in ROM50I and ROM25I (Fig. 4c, d, g, h). A warm bias is found over Central America in summer in both REMO50I and the coupled runs, with stronger magnitude in the latter. Over the ocean, the spatial bias pattern resembles the MPI-ESM forced runs, with more positive values. This bias has been linked to the misrepresentation of key processes in the gulf currents (Jungclaus et al. 2006;Kwon et al. 2010).

Precipitation
The simulated rainfall over the Mesoamerica region is compared with GPCP and TRMM observational precipitation data. Figure 2e, f show the spatial distribution of the average DJF and JJA daily mean rainfall. In DJF, the ITCZ is in its southernmost position and dry conditions prevail over most of Mesoamerica. These conditions are not spatially uniform, as southern Central America and its Caribbean coast present wetter conditions. In summer, the ITCZ reaches its northernmost position and rainfall along the eastern Pacific coast is enhanced by the ITCZ (Karnauskas and Busalacchi 2009). The North America Monsoon leads to strong precipitation in southwestern United States and northwestern Mexico (Gochis et al. 2006). MPI-ESM shows biases that reflect a poleward shift of a weakened ITCZ, such as intense rain in the eastern Pacific and northeastern México (Fig. 5a, e). Moreover, MPI-ESM overestimates the precipitation in the Caribbean, and simulates an unrealistic rainfall distribution between land and ocean (Stevens et al. 2013) resulting in dry biases over Central America and eastern Mexico. Thus, the precipitation bias pattern is stronger in the wet season (JJA), than in the dry season (DJF).
The REMO and ROM (Fig. 5b, f, i, l) improves substantially the precipitation bias of MPI-ESM except in summer over Central America, where the dry bias is similar and even slightly reinforced by ROM (Fig. 5f). Moreover, ROM simulations provide a better description of the meridional rainfall distribution in southernmost Central America despite the a drier ITCZ (Fig. 5b, f), as they can represent the Central American Dry Corridor (CADC, Hidalgo et al. 2017). Over the ocean, both regional models improve the ITCZ representation. Still, REMO (ROM) shows a stronger ITCZ and a wet (dry) bias in the Pacific coast of Central America. In general, both MPI-ESM and REMO simulations show the largest wet bias in summer (Fig. 5e, l-n), while for ROM a general improvement of the precipitation bias is found. This is particularly evident in summer (Fig. 5f, h). Nevertheless, this reduction of the precipitation leads to a dry bias in some regions. Concluding, the ocean-atmospheric coupling improves the representation of the ITCZ as it permits the air-sea interactions important to properly represent the energy and moisture transfers. These results are consistent with Kang et al. (2008).  Figure 6 shows the SST biases of MPI-ESM and RCM simulations with respect to the Reynolds OISSTV2 climatology (Fig. 2g, h). In both summer and winter, MPI-ESM exhibits a warm bias off the California Gulf, but with larger values in summer. Moderate warm biases are also present in the coastal regions of the Gulf of Mexico, while cold biases are found in the eastern part of the Gulf of Mexico and in the East Pacific Warm Pool region during winter. All the RCM simulations show colder biases over the open oceans and the Gulf Stream region, while the warm biases off the Californian Coast are alleviated. This effect is stronger in ROM50M. The cold bias in the Caribbean, Gulf of Mexico, and Gulf Stream is stronger in ROM50M than in ROM50I. This bias might be associated with stronger trade winds in ROM, leading to enhanced evaporation, mixing and Ekman pumping beneath cyclonic wind shear (de Szoeke and Xie 2008). Over the Atlantic Warm pool, the bias could be partially related to radiative flux errors (Lin 2007;Liu et al. 2012). ROM25I has the lowest bias, showing that the coupling and the higher resolution lead to more reliable SSTs, which may have a relevant role as the coupling enhances the representation of the SST distribution, hence improving the simulation of the WHWP.

Sea surface temperature
Following , we select four subdomains for the temporal analysis of oceanic variability: Gulf of Mexico (GM), Caribbean Sea (CS), Western Atlantic (WA) and Eastern tropical Pacific (EP; cf. Fig. 1b). In the simulations forced by ERA-Interim, we expect the model SSTs to be at least partially in phase with the observations, given that the large-scale climate signal propagates in the modeled region in two ways: (1) through the surface fluxes that force the ocean outside the region of coupling, and (2), through the lateral boundary conditions that force the regional atmospheric model. However, since part of the SSTs variability is locally-generated and thus independent from the large-scale signal, we do not expect a perfect match between modeled and observed SSTs. Moreover, we also consider the Niño3 region, which is outside of the coupled domain and thus do not include any ocean-atmosphere coupling. Figure 7 shows the time series of the area-averaged SST anomaly over Niño3, GM, CS, WA, and EP from observations and ROM runs. The simulated variability over the Niño3 region is in good agreement with the observations. In fact, the correlation between the time series of OISSTV2 and both ROM simulations is as high as 0.9, though the simulated standard deviation is higher than observed. This suggests that ROM is able to better represent the oceanic influence from outside of the coupling region. The ocean alone forced experiments show slightly different results, suggesting that the coupling region is able to influence the ocean also outside of this region. Generally, the magnitude of this remote influence depends on the domain location and size (e.g., Sein et al. 2014). In the EP box, the SST standard deviation (0.98 °C for observations and 0.88 °C for both ROM50I and ROM25I) is higher than over the Niño3 and the other selected regions (Fig. 8b). However, the temporal correlation with OISSTV2 is still very high (0.78 for ROM50I and 0.79 for ROM25I). For the other boxes, the agreement is lower. In particular, the SST standard deviation WA box (Fig. 7e) is 0.36 °C for observations and 0.40 and 0.42 °C for ROM50I and ROM25I simulations, and the correlation is 0.45 and 0.36, respectively.
It is noteworthy that the correlation for the EP box (Fig. 7b) is higher than estimated by . In this case, the box is close to the coupling domain boundary, and the small difference in the correlation can be due to the effects on ocean circulations introduced in ROARS by the oceanic boundary conditions scheme. The other boxes are located farther from the boundaries, and the smaller correlations could reflect the stronger internal variability. The magnitude of SST standard deviation is similar to the observed values, reflecting a correct simulation of the climate statistics. It is also interesting to note that in this case higher resolution does not bring a significant SST improvement in the selected regions.

Sea surface height
The standard deviation of the SSH from AVISO, MPI-ESM, and ROM is shown in Fig. 8. Since SSH variability strongly depends on resolution (c.f. Sein et al. 2017), we would expect that MPI-ESM and CFSR would be unable to capture the SSH variability as featured in AVISO (compare Fig. 8c, e to Fig. 8a), as both datasets have a relatively coarse resolution. Note that SSH data was not assimilated in the CFSR (Xue et al. 2014). In contrast, ROM is able to capture the main features of the spatial distribution of the SSH variance over these regions (Fig. 8b, d, f). Still, the ocean resolution is not eddy resolving in ROM, and the model fails to show the high SSH variance related to the characteristic gyres of the Gulf of Mexico circulation. Furthermore, ROM overestimates the variability in the Gulf Stream region, probably due to the instabilities generated by the erroneous east-west  Fig. 1c for region definitions). OISSTV2 in black, ROM50I simulation in green, and ROM25I in blue SST gradient across the Gulf Stream (see Fig. 6). This shortcoming must be taken into consideration when analyzing e.g. the easterly flow associated with eddies propagation and its influence on marine ecosystems, since they are highly dependent on the gyres, which are not very well resolved by the model.

Temperature and precipitation in the Köppen-Geiger areas
For a more detailed study of the seasonal cycle in the Mesoamerica region, we use the Köppen-Geiger climate types regionalization proposed by Rubel et al. (2017). From the 18 climate types present in Mexico and Central America, we choose only 9 (see Fig. 1c; Table 2) since the areas of the remaining are small and not well represented in the models or/and the observations. It must be noted that the Köppen-Geiger procedure considers exclusively the seasonal temperature and precipitation regimes, regardless of the physical agents shaping them at each location. This can result in areas of a common type that are actually governed by different processes. We believe this drawback to be partly reduced, on the one hand, due to the limited extent of our domain. On the other hand, we also consider that the widespread use of the Köppen-Geiger classification in many climate-related fields makes it interesting to examine the ability of our simulations in reproducing its features, implicitly including its spatial distribution. Figure 9 shows the seasonal cycle of T2M averaged over the selected regions for all simulations, ERA-Interim and CFSR reanalyses, and the CRU observational dataset. This last dataset is used in Fig. 10 as reference to construct the Taylor Diagrams for the time series presented in Fig. 9. In the equatorial climate regions (i.e., Af, Am, As and Aw), the weak annual cycle is well reproduced by most models but with generally higher variance. MPI-ESM deviates from the observed annual cycle and produces a second T2M maximum in summer. This bi-modal distribution is transferred to the ROM50M and REMO50M simulations, resulting in the artificial summer maximum stronger in the latter. These three runs tend to display lower (however still large) correlations in the corresponding Taylor diagrams (Fig. 10a-d), with ROM50M improving the MPI-ESM correlation for all regions with equatorial climate except Af. The pronounced annual cycle of temperature in the arid (BSh, BSk, and BWh) climates appears to be well captured in amplitude and timing in all runs, which is reflected in the good performance in the Taylor diagrams ( Fig. 10e-g). The peak temperatures seem to be slightly delayed (by about 1 month) in MPI-ESM but well simulated by ROM50M and REMO50M. In the warm temperate regions (Csa and Cwb), the models tend to display a somewhat stronger seasonal cycle than the one represented by CRU, resulting in higher variances Fig. 8 SSH variability, measured by its standard deviation for observations, AVISO, a global coupled reanalysis, CFSR, the driving model, MPI-ESM, and the coupled regional models, ROM50M,ROM50I,ROM25l in the Taylor diagrams (Fig. 10h, i). For these regions, all regional models show a better variance and correlation than MPI-ESM. For the Csa climate, T2M falls comparatively quickly, in all runs but MPI-ESM, from the early summer maximum to early autumn, as opposed to the higher persistence in CRU. It must be noted that CFSR and ERAI share the same behavior of the models. In Cwb, MPI-ESM again shows too high temperatures in late summer, departing from  Fig. 1c for all the REMO and ROM simulations, MPI-ESM and the ERA-Interim and CFSR reanalysis the observed annual cycle. With much lower mean temperature and a bit less noticeably, ROM50M and REMO50M reproduce this feature.
In Fig. 11 we compare the seasonal cycle of precipitation in the selected Köppen-Geiger areas, for the MPI-ESM, REMO, ROM simulations and for the ERA-Interim and CFSR reanalyses, against three observational datasets: CRU, CHIRPS, GPCP, and TRMM. We use these four different datasets for validation because it is known that the scarcity of stations and some uncertainties in the satellite algorithms used for precipitation estimates generate some differences in the resulting gridded precipitation. Figure 12 displays Taylor diagrams for the seasonal cycle in each region relative to the seasonal cycle featured by the CRU precipitation data. As could be expected, the seasonal cycle of total precipitation shows a high spread in most of the selected climate type regions. In Af, Am, and Aw equatorial climate regions, the precipitation presents a pattern typical of the MSD, with two maxima in June and September. In these regions MPI-ESM and ROM show less precipitation than both the observations and all the reanalysis while the simulated REMO precipitation tends to be stronger, leading to pronounced deviations Fig. 10 Taylor diagram of the seasonal cycle of the 2 m temperature in nine representative climate types for the Central America and Mexico area. The variance of the seasonal cycle in each region is represented against the correlation with CRU data from the CRU reference in the Taylor diagrams of Fig. 12. The three REMO experiments exhibit a higher variance in the Aw-type Taylor diagram (Fig. 12d), reflecting the wet bias in summer along the Pacific slope of Central America (Fig. 5l-n). This wet bias is also seen in the Af-type diagram (Fig. 12a), and only for REMO50M in the As-type (Fig. 12c), as its extension is smaller in the ERA-Interim driven simulations and does not affect the As-type area. It is noteworthy that CFSR also tends to show a stronger than observed precipitation, not only for the A types, but for the others as well, while the other reanalysis generally stays closer to observations. The coupled simulations, in agreement with the results in Fig. 5b-h, suffer from a dry bias in the A-type areas, except for the As-type, noticeable in the seasonal cycles of Fig. 11a, b, d and in the smaller variances in the corresponding Taylor diagrams (Fig. 12a, b, d). As the precipitation in regions depends heavily on the humidity transport from the Caribbean, the lower precipitation in Fig. 11 Seasonal cycle of the precipitation in nine representative climate types for the Central America and Mexico area. The precipitation is averaged over the corresponding regions represented in Fig. 1c for all the REMO and ROM simulations, MPI-ESM, the CFSR reanalysis and the CRU, TRMM, GPCP and CHIRPS observed datasets the coupled simulations can be explained by the suppressed evaporation caused by the colder SST in the Caribbean Sea. The smaller precipitation rates in the arid B-type areas lead to less spread in the seasonal cycles of Fig. 11e-g, though the Taylor diagrams reveal some deficiencies in the simulations. For the BSh type, the REMO50M experiment shows a much higher variance, due to the summer wet bias over eastern Mexico (Fig. 5l). The poorest performance in terms of correlation and variance is exhibited by the coupled simulations for the BWh type. This is attributable to the wet bias along the mainland coast of the Gulf of California which affects these runs and, more severely, also the REMO50M and the MPI-ESM experiments (Fig. 5a-d, i). For the warm temperate C types, the dry bias in summer over southern Mexico results in lower variance in the Taylor diagram of the Cwb type (Fig. 12i), correcting the excess in the uncoupled runs. The correlation is high, though, particularly for ROM25I, hinting at a correct timing in the seasonal march of precipitations in these mountainous regions. The reverse occurs for the Cwb type seasonal cycle, whose amplitude is well captured by the ROM experiments, despite a low correlation in the Taylor diagram of Fig. 12h. Overall, ROM Fig. 12 Taylor diagram of the seasonal cycle of the precipitation in nine representative climate types for the Central America and Mexico area. The variance of the seasonal cycle in each region is represented against the correlation with CRU data offers a good representation of the seasonal cycle in precipitation over the most extended A types, generally correcting the wet bias of the uncoupled runs and improving the MPI-ESM precipitation. In the wetter Af and Am areas, ROM50M partially mitigates the dry bias of MPI-ESM in the second half of the year, with a performance similar to that of the ERA-driven experiments (Fig. 11a, b). It also amends the early onset of the wet season in the As-type displayed by MPI-ESM in Fig. 11c, once more leading to a behavior comparable to ROM50I.

Mid summer drought
The MSD over Mexico and Central America (Magaña et al. 1999) is characterized by a bimodal distribution in the annual cycle of precipitation over southern Mexico and Central America with maximum values during June and September-October and a relative minimum during July and August.
The MSD has been related with fluctuations in the intensity and location of the eastern Pacific ITCZ (Magaña et al. 1999), and with the July maximum of the CLLJ (Amador 2008). The increase in the magnitude of the convective forcing causes the intensification of the trade winds over the Caribbean during July and August, which reinforces the transport of moisture from the Caribbean. This intensification in the moisture transport combines with the orographic forcing of the mountains over most of Central America to produce the maximum precipitation along the Caribbean coast, which contrasts with minimum precipitation along the Pacific coast of Central America ). On the Pacific side, changes in the divergent (convergent) low-level winds over the "warm pool" off the west coast of southern Mexico and Central America determines the evolution of the MSD (Small et al. 2007). Figure 13 shows the representation of the MSD pattern as the difference in mean precipitation between July and August (JA) minus June and September (JS) for GPCP, TRMM and the different coupled and uncoupled simulations. Thus, positive values occur in regions where July and August dominate the wet season and negative values where precipitation is reduced during these months. Although the patterns of the MSD in both observations agree in the large-scale, TRMM (Fig. 13b) shows a wetter signal over the Eastern Central America, the ITCZ region, and the North-western Mexican coast, the latter being linked to an enhanced monsoon over the western Sierra Madre mountain range. In contrast, the TRMM is dryer than GPCP over the Gulf of Mexico. In the higher resolution of TRMM, it is clearly seen a south-north wet-dry dipole that is related to fluctuations in the position of the eastern Pacific ITCZ, which is much weaker in GPCP.
The MPI-ESM reproduces well the spatial distribution of the large-scale MSD pattern and is quite similar to GPCP (compare Fig. 13a and Fig. 13c), although the magnitude is stronger, especially in the region located southwest of Baja California Peninsula where there are large latitudinal Fig. 13 Representation of the Mid-Summer Drought (July and August mean minus June and September mean) in a GPCP, b TRMM, c MPI-ESM, the uncoupled simulations, d REMO50M, e REMO50I, f REMO25I, and the coupled simulations, g ROM50M, h ROM50I, i ROM25I displacements of the ITCZ. The MPI-ESM fails to simulate the wet-dry precipitation pattern between eastern and western Central America and the northward extension of the MSD pattern in Northeastern Mexico. Also, because of its low resolution, the MPI-ESM is unable to simulate the finer scale MSD structures that can be seen in TRMM (Fig. 13b). In contrast, the higher resolution of all coupled and uncoupled regional simulations is able to resolve most of these features, especially over land in the coastal regions ( Fig. 13d-i). Moreover, the regionally coupled simulations represent better the east-west contrast in the precipitation regimes of Mexico, which is characterized by a bimodal distribution along the west coast and a single peak distribution along the east coast. In Fig. 13 this difference results in wet and dry contrast between the eastern and western part of Mexico.
It is noteworthy that this east-west contrast has an opposite sign in the southern Central America region, and that this feature is missing in the global simulation (i.e., MPI-ESM) and strongly underestimated by the REMO model forced by simulated SST (Fig. 13d). These results show the importance of the higher resolution but also of the appropriate representation of the SST, which is obtained by either prescribing the field (i.e., REMO50I and REMO25I) or coupling the atmosphere to an active ocean (i.e., ROM50M). Moreover, the importance of the coupling is shown by the smaller ROM biases over the ITCZ latitudinal displacement region and the coastal regions of Northwestern Mexico. Concluding, the results of Fig. 13 point to the importance of the coupling (especially over the ocean) and higher resolution (especially over mountainous regions) for accurately simulate the MSD. Figure 14 shows the annual cycle of precipitation over two domains. One is located in the East Pacific Warm Pool (WP), which covers the region 105W-95W and 10N-15N and represents the region with the warmest water in the eastern Pacific Ocean and is located directly adjacent to southern Mexico and Central America (see SST climatology in Fig. 2g, h). The other domain extends from 94W-87W and 14N-22N and covers southeastern Mexico (SEM), including land and waters from the Gulf of Mexico, Caribbean Sea and eastern Pacific Ocean (see Fig. 1b). In the WP region, ROM25I reproduces the amplitude and phase of the observed annual cycle of precipitation. This agreement decreases when lower spatial resolution is used (ROM50I), and it is even weaker when ROM is forced by . In this last case, ROM50M clearly underestimates the amplitude of the seasonal cycle of precipitation.
The uncoupled simulations and the global model, however, provide a large overestimation of the amplitude of the seasonal cycle of precipitation, illustrating the key role of both spatial resolution and coupling for an accurate simulation of the evolution of precipitation in the WP region. In the SEM region, however, the geographical location is more complex and the global model and the coupled simulations reproduce very well the phase of the seasonal cycle of precipitation but underestimate its amplitude, while the stand-alone simulations reproduce very well the phase and the amplitude. This seems to be related to the cooling bias of the coupled simulations in the Caribbean Sea and deserves a more detailed analysis to isolate the air-ocean-land processes influencing the amount of precipitation in the SEM region.

Caribbean Low Level Jet
The CLLJ (Amador 2008) is characterized by strong easterly zonal winds in the Caribbean region that varies semi-annually, with two maxima in the summer and winter and two minima in the fall and spring. During summer, the strengthening of the CLLJ is associated to a maximum of sea level pressure (SLP), a relative minimum of rainfall (the midsummer drought), and a minimum of tropical cyclogenesis in July in the Caribbean Sea (Wang 2007). The CLLJ acts as moisture conveyor from the tropical Atlantic into the Caribbean Sea, Gulf of Mexico, and the continental US, hence, influencing rainfall both locally in the Caribbean and Central America, and remotely in the US (see e.g., Durán-Quesada et al. 2017).
Figure15 shows the winter and summer climatology of 925 hPa winds from ERA-Interim, MPI-ESM and our simulations. During winter, the region with MPI-ESM winds in excess of 10 m/s along the northern coast of Colombia and Venezuela features a CLLJ that extends too far east when compared with the reanalysis. In the west, the MPI-ESM fails to simulate the zonal wind extension to the Costa Rican Dome (CRD; see e.g., Fiedler 2002), which is located over the Pacific Ocean in the region off of Costa Rica, which results from the wind funneling through the Papagayo Gap during the CLLJ winter maximum. In this case, this bias in the spatial representation of the zonal wind arises from a resolution issue, as the simulation even at its highest resolution is not able to resolve the Papagayo Gap. Hence, the winter intrusion of the CLLJ is not resolved. Both the east and west bias in the location of the CLLJ are improved by REMO and ROM that show a more realistic representation of the 925 hPa wind. This improvement seems to be largely independent from the driving field and more linked with the incorporation of the air-sea coupling.
During the boreal summer ( Fig. 15i-p), the MPI-ESM simulates a weaker CLLJ when compared with the reanalysis. This bias is reproduced by REMO50M, while ROM50M shows a more realistic CLLJ pointing to the importance of the SST gradients in this season. Moreover, REMO50I, which is forced by the ERA-Interim SST, simulates a stronger CLLJ, although weaker than ROM50I, pointing to the relevance of the SST gradients and ocean-atmosphere feedbacks rather than only the large-scale pressure forcing. A comparison of Fig. 15e, f and Fig. 15g, h for winter and Fig. 15m, n and Fig. 15o, p for summer shows that for both seasons, the impact of atmospheric resolution is higher in boreal winter, which presents a stronger CLLJ core. Hence, surface forcing is more relevant for the depiction of the CLLJ compared to the summer peak, supporting previous considerations on forcing of the CLLJ maximum during winter and summer .
Outside of the CLLJ region, the winter and summer climatologies of easterly winds over the Caribbean Sea and the northern tropical Atlantic simulated by MPIESM, ROM, and REMO are comparable to ERA-Interim, while they tend to be stronger in the eastern Pacific in winter. Considering the ITCZ position and intensity to be linked with the CLLJ (as suggested by Hidalgo et al. 2015), it is expected that deficiencies in the simulation of the ITCZ appear along with the misrepresentation of the CLLJ.
An index of the CLLJ intensity is computed for the models, the ERA-Interim and CFSR reanalysis. Following the definition by Wang (2007), the absolute value of the zonal wind in the region is averaged between 80W-70W and 12.5N-17.5N (see Fig. 16a). Figure 16b shows the performance of models and Reanalysis against the corresponding 925 hPa monthly wind climatology from IGRA radiosonde data over Plesman Field (Hato airport, Curaçao, located at 68.97W and 12.2N). The CLLJ is consistently stronger in CFSR than in ERA-Interim, probably due to both the resolution (~ 38 km for CFSR and ~ 80 km for ERA-Interim) and the interactive coupling between the ocean and the atmosphere. MPI-ESM simulates reasonably well the strength of the CLLJ in winter and spring, remaining close to CFSR while in summer and autumn it becomes weaker than in both reanalyses. The behavior of the regional models depends strongly on the coupling: the CLLJ in coupled simulations is consistently stronger, independently of the resolution and shows little differences between them while the CLLJ in REMO is weaker and shows large different magnitudes, especially in the second half of the year. In winter and spring, the CLLJ in both coupled and uncoupled simulations show values closer (and stronger) to CFSR and ERA-Interim respectively. While in summer and autumn the CLLJ in coupled simulations remain close to CFSR, the uncoupled simulations show a strong weakening with much more dispersion in the values. This could point to the importance of ocean-atmosphere interaction in the region: while the CLLJ remain very close in the coupled models independently of the driving model, the ERA-Interim uncoupled simulation show a stronger CLLJ than the simulation forced by MPI-ESM. The atmospheric grid resolution also seems to influence the CLLJ strength: in the 25 km simulations the CLLJ is stronger than in those with 50 km resolution.
It is seen that the depiction of CLLJ over Plesman field in the CFSR and ROM is reasonable (Fig. 12b). The seasonal peaks in February and June are well represented in the Reanalysis and the simulations forced by ERA-Interim while MPI-ESM and the REMO50M fail to capture the June peak. These results could point to an interplay between the local sea-atmosphere feedbacks and the large-scale processes as variations of the North Atlantic Subtropical High (NASH) and SST in the Pacific and Atlantic. In particular, the REMO simulations show that a correct simulation of the NASH is important, as the biases in the large-scale MPI-ESM fields Fig. 16 a The representation of the Caribbean Low Level Jet (CLLJ) in the driving MPI-ESM, ERA Interim, the coupled CFSR, ROM and REMO. The CLLJ index is calculated as the average of the 925 hPa zonal wind in the 80W-70W, 12.5N-17.5N box. b Validation of the 925 hPa zonal wind against the IGRA radio sonde data at Hato airport leads to a stronger than observed weakening of the CLLJ winds in the second half of the year. Meanwhile, the winds in the MPI-ESM forced ROM simulations are stronger than those of the driving model, and the seasonal cycle is similar to the ERA-Interim forced simulations. In conclusion, the performance of uncoupled and coupled simulations strongly depend on the observational dataset used as reference. While the uncoupled simulations compare with ERA-Interim better than with CFSR, the opposite is true for the coupled simulations.

Summary and discussion
We investigated the role of driving boundary conditions, horizontal resolution, and ocean-atmosphere coupling in coupled (ROM) and uncoupled (REMO) regional simulations of the present-time climate of the Mesoamerican region. The role of the boundary conditions is investigated by forcing REMO and ROM with the global model MPI-ESM and the reanalysis ERA-Interim.
Compared with observations, the MPI-ESM SLP field over the Pacific Ocean is characterized by a weak Aleutian Low in winter and a weak North Pacific High in summer. The uncoupled simulation forced by MPI-ESM reproduces the same large-scale pattern of the driving field (see Fig. 3e-g) including the bias over the Pacific Ocean. Over land, the better representation of the orography over the mountain ranges of North America leads to a different SLP pattern compared to the MPI-ESM. Over the Atlantic Ocean, REMO presents a tendency to simulate a weaker NASH. In ROM, the coupling causes generally colder SSTs, which leads to an overall SLP increase over the ocean, modifying the circulation pattern of the driving model both over the eastern Pacific and the Caribbean. Moreover, it changes the SLP over land, especially in summer, probably due to the advection of colder air from the ocean. This effect of the coupling can be observed both in the ERA-Interim and MPI-ESM driven simulations, suggesting that it is caused by the incorporation of the ocean-atmosphere interactions in ROM. Regarding the effect of higher resolution on the simulated SLP, results provide evidence that a better-resolved Mesoamerican relief improves the simulated SLP over land during the more convectively active summer season. As might be expected, the improvement over ocean is much smaller.
Regarding T2M, both REMO and ROM are generally cooler than MPI-ESM, especially over western North America and Mesoamerica. Over the ocean, REMO reproduces the biases of the driving model, while the coupling cools the surface air temperature, reducing the positive bias that MPI-ESM exhibits along the coasts of Baja California and the Gulf Stream, although T2M cold biases develop in the Caribbean Sea and western North Atlantic over a colder SST, and in the Tropical Pacific due to a displaced ITCZ, especially in JJA. Generally, the ERA-Interim forced simulations show the smallest T2M biases, both for ROM and REMO and the higher resolution seems to add little to the large-scale patterns of surface temperature biases. Contrary to the simulations driven by MPI-ESM, the downscaled surface temperature over land is generally warmer than the driving reanalysis, except over the mountainous regions in western North America and Mexico, which display cooler T2M, much more noticeably in ROM50I and ROM25I. A warm bias is seen over Central America in summer in both REMO50I and the coupled runs, stronger in the latter. Over the ocean, the spatial pattern of the biases resembles these of the MPI-ESM forced runs, albeit with more positive values. We would like to note that part of the cold biases are likely due to the effect of the difference in orography between the model and the reanalysis. The probable dominance of low height stations in the assimilated data might also contribute to it, in line with findings in previous RCM studies (Karmalkar et al. 2011).
The regionalization of the land surface with the help of the Köppen-Geiger climate types allows us to study with more detail the simulated seasonal cycle of temperature and precipitation over land. The relative performance of the models is different across the nine climate types selected for this study, especially for precipitation, as it is influenced by orography. The models' biases influence the representation of the seasonal cycle in precipitation in each region, and the coupled simulations tend to revert the wet bias of REMO over the most common A-type climates in the domain.
For the region of study, the scarcity of data for assimilation makes the climate simulated by the global model (both reanalysis and ESM) very dependent on model details and the advantages of both higher resolution and coupling emerge. Moreover, the computed climate variables show substantial differences among them and with independent observational datasets. This fact plays a very important role for the simulation of the regional climate, especially over land. REMO tends to reproduce the large-scale features of the driving model, albeit showing differences over land, as illustrated by our analysis in the climate type regions. Over the ocean REMO is forced by the global model SST and the simulated fields are close to the global ones. On the other side, ROM adds to the ERA-Interim forcing the air-sea interaction which is also better resolved than the global model.

Conclusions
In this study we have analysed the possible benefits of using coupled regional climate models in the representation of the regional climate for the Mesoamerican region. The main conclusions are as follows: • Air-sea coupling changes the simulated fields over the ocean, especially precipitation. On the other hand, higher resolution leads to improvements mainly over land, particularly over steep orography. However, several largescale variables still depend primarily on the driving fields. • Both higher resolution and coupling lead to an added value compared to the global data, both with ERA-Interim and MPI-ESM forcing. Still, REMO and ROM clearly show different skills in simulating the present time climate. Over the ocean, REMO simulated fields are closer to the global ones, as it is forced by the global model SST. On the other side, ROM adds to the ERA-Interim forcing the air-sea interaction, which is also better resolved than the global model. • The coupling brings a significant improvement in the representation of the Intertropical Convergence Zone and the Caribbean Low Level Jet. The CLLJ is consistently stronger in ROM, regardless of the resolution while it is weaker in REMO and shows large spread in magnitude, especially in the second half of the year. • When forced by ERA-Interim, both REMO and ROM show good skill reproducing the climate of Mexico, Central America and the surrounding water masses of the Eastern tropical Pacific and the Caribbean Sea. Overall, ROM improves the representation of the Mid-Summer Drought, the CADC, and, generally, the seasonal cycle in surface temperature and precipitation, correcting the wet bias of the uncoupled runs over the most extended regions of humid climate and improving the distribution of rainfall in southern Central America. • Driven by MPI-ESM, ROM reduces the dry bias of the global model over the wetter areas in the domain, leading to a seasonal cycle in precipitation comparable to that reproduced by the ERA-Interim-forced experiments.
Future work will focus on a better understanding of the physical processes associated leading to a better representation of the regional distribution of rainfall in the coupled simulations, mainly for the WHWP, the linkages between the CLLJ and the ITCZ, as well as the relevance of surface fluxes. Moreover, possible changes on the regional climate for the Mesoamerican region will be evaluated using climate change projections according to the RCP 8.5 scenario for the XXI-century.