HighResMIP versions of EC-Earth: EC-Earth3P and EC-Earth3P-HR. Description, model performance, data handling and validation

A new global high-resolution coupled climate model, EC-Earth3P-HR has been developed by the ECEarth consortium, with a resolution of approximately 40 km for the atmosphere and 0.25 degree for the ocean, 25 alongside with a standard resolution version of the model, EC-Earth3P (80 km atmosphere, 1.0 degree ocean). The model forcing and simulations follow the HighResMIP protocol. According to this protocol all simulations are made with both high and standard resolutions. The model has been optimized with respect to scalability, performance, data-storage and post-processing. In accordance with the HighResMIP protocol no specific tuning for the high resolution version has been applied. 30 Increasing horizontal resolution does not result in a general reduction of biases and overall improvement of the variability, and deteriorating impacts can be detected for specific regions and phenomena such as some EuroAtlantic weather regimes, whereas others such as El Niño-Southern Oscillation show a clear improvement in their spatial structure. The omission of specific tuning might be responsible for this. The shortness of the spin-up, as prescribed by the HighResMIP protocol, prevented the model to reach 35 equilibrium. The trend in the control and historical simulations, however, appeared to be similar, resulting in a warming trend, obtained by subtracting the control from the historical simulation, close to the observational one. 40 https://doi.org/10.5194/gmd-2019-350 Preprint. Discussion started: 2 March 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The CMIP6 protocol requests modeling groups to use specific forcing datasets that are common for all participating models. Table 1 lists the forcings that have been implemented in EC-Earth3P(-HR). Because of the HighResMIP protocol, EC-Earth3P(-HR) distinguish themselves in several aspects from the model configurations 125 used for the CMIP6 experiments (Doescher et al., 2019): 1. The stratospheric aerosol forcing in EC-Earth3P(-HR) is handled in a simplified way that neglects the details of the vertical distribution and only takes into account the total aerosol optical depth in the stratosphere which is then evenly distributed across the stratosphere. This approach follows the treatment 130 of stratospheric aerosols as it was used by EC-Earth2 for the CMIP5 experiments yet with the stratospheric aerosol optical depth (AOD) at 500 nm updated to the CMIP6 data set.
2. A sea surface temperature (SST) and sea-ice forcing data set specially developed for HighResMIP is used for AMIP experiments (Kennedy et al., 2017). The major differences compared to the standard SST forcing data sets for CMIP6 are the higher spatial (0.25 deg vs. 1 deg) and temporal (daily vs. monthly) The aim of the performance activities for EC-Earth3P-HR is to adapt the configuration to be more parallel, scalable and robust, and to optimize part of the execution when this high-resolution configuration is used. The

165
performance activities are focused on three main challenges: (1) scaling of EC-Earth3P-HR to evaluate the ideal number of processes for this configuration, (2) analyses of the main bottlenecks of EC-Earth3P-HR and (3) new optimizations for EC-Earth3P-HR.

170
The results of the scalability analyses of the atmosphere (IFS) and ocean (NEMO) components of EC-Earth3P-HR are shown in Fig. 1, and for the fully coupled model in Fig. 2. Acosta et al. (2016) showed that, while for coupled application the load balance between components has to be taken into account in the scalability process, the process needs to start with a scalability analysis of each individual component. However, if the optimization 175 of one component (e.g. the reduction of the execution time of IFS) does not reduce the execution time of the coupled application, because of other slower components, a load balance analysis is required. The final choice depends on the specific problem, where either time or energy can be minimized. In section 3.2, we describe how the optimal load balance between the two components, where NEMO is the slowest component, was achieved (Acosta et al., 2016).

Bottlenecks
For the performance analysis, the individual model components (IFS, NEMO and OASIS) are benchmarked and 185 analyzed using a methodology based on extracting traces from real executions. These traces are displayed using the Paraver software and processed to discover possible bottlenecks (Acosta et al., 2016). Eliminating these bottlenecks not only involves an adjustment of the model configuration and a balance of the number of cores devoted to each one of its components, but also modifications of the code itself and work on the parallel programming model adopted in the different components.

190
The first step of a performance analysis consists in analyzing parallel programming model codes using targeted performance tools. Figure 3a illustrates an example of the performance tool's output from one single EC-Earth3P-HR model execution as provided by the Paraver tool, focusing only on its two main components: NEMO and IFS.  2) From other parts of the application, we also notice the expensive cost of the IFS output process for each time step. A master process gathers the data from all MPI sub-domains and prints the complete outputs at a regular time interval of three and six hours. During this process, the rest of processes are waiting for this step to be 225 completed. Due to the large data volumes, this sequential process is very costly, increasing the execution time of IFS by about 30% when outputs are required.
3) The bottom part of Fig. 3a shows that the communication in NEMO is not very effective and that a large part of it is devoted to global communications, which appear in purple. Those communications belong to the horizontal 4) Due to the domain decomposition used by NEMO some of the MPI processes, which are used to run part of 235 the ocean domain in parallel, were computing without use. This is because domain decomposition is done on a regular grid and a mask is used to discriminate between land and sea points. The mask creates subdomains of land points whose calculations are not used. This is illustrated in Fig. 4 showing a particular case in which 12 % of the depicted subdomains do not contain any sea-point. This is because all-to-one/one-to-all MPI communications are replaced by global communications (gather/scatter and reduction) and the coupling calculations are done by all the IFS processes instead of only the IFS master process.

255
Another functionality of OASIS consists in gathering all fields sent from IFS to NEMO in a single group (Acosta et al., 2016). Coupling field gathering, an option offered by OASIS3-MCT, can be used to optimize coupling exchanges between components. The results show that gathering all the fields that use similar coupling transformations reduces the coupling overhead. This happens because OASIS3-MCT is able to communicate and interpolate all of the fields gathered at the same time.
260 Figure 3b shows the execution when "opt" and "gathering" options are used, with the 90% reduction in coupling time clearly visible (large green section). In the case of the first time step in the trace, the coupling time is replaced by waiting time, since NEMO is finishing its time step and both components have to exchange fields at the end of the time step.

265
2) For the output problem, the integration of XIOS as the I/O server for all components of EC-Earth can increase performance dramatically. XIOS is already used for the ocean component NEMO and the I/O server receiving also all the data from IFS processes and doing the output work in parallel and in an asynchronous way is the best solution to remove the sequential process when an IFS master process is required to do this work. This is being developed and will be included in the next version of EC-Earth.  processor cores 41% faster using 25% less resources). The increase in throughput was due to less computations and related to that less communications. In addition, ELPiN allows for the optimal use of the available resources

305
CMOR-library for the production of netCDF files with the appropriate format and metadata. The latter is the only supported step for the ocean output.
To speed up the atmosphere post-processing, the tool can run multiple CDO commands in parallel for various requested variables. Furthermore, we optimized the ordering of operations, performing the expensive spectral 310 transforms on time-averaged fields wherever possible. We also point out that the entire procedure is driven by the data request, i.e. all post-processing operations are set up by parsing the CMOR tables and a single dictionary relating EC-Earth variables and CMOR variables. This should make the software easy to maintain with respect to changes in the data request and hence useful for future CMIP6 experiments. including O3 and aerosol loading for a 1950s (~10 year mean) climatology. Output from the initial 50 year spinup is saved to enable analysis of multi-model drift and bias, something that was not possible in previous CMIP exercises, with the potential to better understand the processes causing drift in different models.
-Control simulation; control-1950 geopotential height at 500 hPa (Z500) during DJF is well represented with the exception of the region south of Greenland, which is consistent with the MSLP bias (Fig. 9).
Doubling of the atmospheric horizontal resolution has only modest impact on the large-scale structures of the main meteorological variables, as illustrated by the global MSLP, SAT, and precipitation (Fig. 10). For SAT the

spinup-1950
As discussed in the outline of the HighResMIP protocol, the spin-up was started from an initial state that is based For the ocean no data assimilation has been performed, which can result in imbalances between the density and velocity fields giving rise to initial shocks and waves.

430
This drift also has an impact on the fast components of the climate system, which therefore still might reveal trends.

control-1950
435 After the spin-up the SAT each of the three members of EC-Earth3P-HR is in quasi-equilibrium and the global mean temperature oscillates around 13.9 °C (Fig. 11-left, black). The ocean is still warming as expressed by a negative net surface heat flux in the order of -1.5 Wm -2 (positive is upward) (Fig. 11-right, black). This imbalance is reduced during the simulation, but without an indication that the model is getting close to its equilibrium state.

440
Contrary to EC-Earth3P-HR, the global annual mean SAT of EC-Earth3P displays a significant upward trend, with an indication of stabilizing after about 35 years ( Fig. 11-left, red). This warming trend is caused by a large warming of the North Atlantic as revealed by Fig. 12 showing the difference between the first and last 10 years of the control-1950 run. This warming is caused by the activation of the deep convection in the Labrador Sea (not shown), that started about 10 years after the beginning of the control simulation, which was absent in the spin-up state does not strongly affect the slow warming of the deeper ocean, which is reflected in a similar behavior of the net surface heat flux as for EC-Earth3P-HR (Fig. 11-right).
The control-1950 experiment is also analyzed to evaluate model performance of internally-generated variability 450 in the coupled system; the targets are: El Niño-Southern Oscillation (ENSO), the North Atlantic Oscillation (NAO), sudden stratospheric warmings (SSWs) and the Atlantic Meridional Overturning (AMOC). Figure 13 depicts the seasonal cycle of the Niño3.4 index (SST anomalies averaged over 5ºS-5ºN/170ºW-120ºW).

455
As it was also shown for EC-Earth3.1 (Yang et al 2019), both EC-Earth3P and EC-Earth3P-HR still have a systematic underestimation of the ENSO amplitude from late-autumn to mid-winter and yield the minimum in July, 1-2 months later than in observations. Increasing model resolution reduces the bias in early-summer (May-June) but worsens it in late-summer (July-August). Overall, EC-Earth3P-HR shows lower ENSO variability than EC-Earth3P, which following Yang et al.'s (2019) arguments suggests that the ocean-atmosphere coupling 460 strength over the tropical Pacific is stronger in the high-resolution version of the model. On the other hand, Fig. 14 displays the spatial distribution of winter SST variability and the canonical ENSO pattern, computed as linear regression onto the Niño3.4 index. Increasing model resolution leads to a reduction in the unrealistic zonal extension of the cold tongue towards the western tropical Pacific, which was also present in EC-Earth3.1 (Yang et al., 2019) and is a common bias in climate models (e.g. Guilyardi et al., 2009): EC-Earth3P reaches longitudes larger the optimal ratio, the more clustered are the data. The sharpness is an indicator of the statistical significance 570 of the regime structure in the dataset in comparison with a randomly sampled multinormal distribution (Straus et al., 2007). The closer the value is to 100, the more significant is the multimodality of the distribution. The sharpness tends to saturate at 100 for very long simulations, so the values reported in Table 3 are obtained from a bootstrap on 30 years randomly chosen. Both the optimal ratio and the sharpness are too low in the EC-Earth3P simulations, as is usually seen for all models. A significant increase with EC-Earth3P-HR is seen for the optimal 575 ratio, and a smaller (non-significant) one is seen for the sharpness. greenhouse-forced warming from 1950 onward in EC-Earth3P can be inferred by subtracting both simulations.

610
The resulting warming pattern compares well with the observed one and is similar to the warming pattern simulated by EC-Earth3P-HR. Due to the transition, the control-1950 does not provide a near-equilibrium state.
It was therefore decided to extend the control-1950 run for another 100 year to allow process studies, that will be documented elsewhere.

615
Analysis of the kinetic energy spectrum indicates that the sub-synoptic scales are better resolved at higher resolution (Klaver et al., 2019) in EC-Earth. Despite the lack of a clear improvement with respect to biases and synoptic scale variability for the high resolution version of EC-Earth, the better representation of sub-synoptic scales results in better representation of phenomena and processes on these scales such as tropical cyclones (Roberts et al., 2019) and ocean-atmosphere interaction along western boundary currents (Tsartsali et al. in 620 preparation).