IPSL-CM5A2 – an Earth system model designed for multi-millennial climate simulations

Abstract. Based on the fifth phase of the Coupled Model Intercomparison Project (CMIP5)-generation previous Institut Pierre Simon Laplace (IPSL)
Earth system model, we designed a new version, IPSL-CM5A2, aiming at running multi-millennial simulations typical of deep-time paleoclimate studies. Three priorities were followed during the setup of the model: (1) improving the overall model computing performance, (2) overcoming a persistent cold bias depicted in the previous model generation and (3) making the model able to handle the specific continental configurations of the geological past. These developments include the integration of hybrid parallelization Message Passing Interface – Open Multi-Processing (MPI-OpenMP) in the atmospheric model of the Laboratoire de Météorologie Dynamique (LMDZ), the use of a new library to perform parallel asynchronous input/output by using computing cores as “I/O servers” and the use of a parallel coupling library between the ocean and the atmospheric components. The model, which runs with an atmospheric resolution of 3.75∘×1.875∘ and 2 to 0.5∘ in the ocean, can now simulate ∼100 years per day, opening new possibilities towards the production of multi-millennial simulations with a full Earth system model. The tuning strategy employed to overcome a persistent cold bias is detailed. The confrontation of a historical simulation to climatological observations shows overall improved ocean meridional overturning circulation, marine productivity and latitudinal position of zonal wind patterns. We also present the numerous steps required to run IPSL-CM5A2 for deep-time paleoclimates through a preliminary case study for the Cretaceous. Namely, specific work on the ocean model grid was required to run the model for specific continental configurations in which continents are relocated according to past paleogeographic reconstructions. By briefly discussing the spin-up of such a simulation, we elaborate on the requirements and challenges awaiting paleoclimate modeling in the next years, namely finding the best trade-off between the level of description of the processes and the computing cost on supercomputers.


boundary layer, cumulus clouds and the diurnal cycle of continental convective rainfall is driven both by (i) the fact that no coupled version of IPSL including the new parameterizations was available when IPSL-CM5A2 project was initialized and (ii) 110 computing cost. The NP parameterizations in LMDZ have been validated for higher vertical resolution (79 levels) and require higher-frequency time coupling between the physics and the dynamics. SP indeed requires a time step for the dynamics of 3 minutes and a coupling every 30 minutes with the physics, whereas dynamics is computed every 2.15 minute and coupled every 15 minutes in NP, leading to significant increase in computation time (see next section). In terms of revision nomenclature, IPSL-CM5A included LMDZ5, revision 2063, whereas IPSL-CM5A2 includes revision 3342. This represents a >2-year gap 115 between both versions.The differences mostly include bug corrections, various improvements in the physics/dynamics interface, and energy conservation. Changes also concern the Emanuel convection scheme, with bug corrections and changes in the upper bound of the deep convection loop.

ORCHIDEE
ORCHIDEE (Krinner et al., 2005) is a dynamic global vegetation model (DGVM) that can run either in stand-alone mode, 120 coupled to LMDZ, or included as the land surface component in the IPSL climate model. ORCHIDEE is made of three parts : (i) the "hydrology module" SECHIBA (Sche´matisation des Echanges Hydriques a`l'Interface entre la Biosphe`re et l'Atmosphère (De Rosnay and Polcher, 1998;Ducoudré et al., 1993), (ii) parameterizations regarding vegetation dynamics based on plant functional types (PFTs) and derived from the LPJ DGVM (Sitch et al., 2003), and (iii) a "carbon module" called STOMATE (Saclay Toulouse Orsay Model for the Analysis of Terrestrial Ecosystems), that simulates vegetation phenology 125 and carbon dynamics. SECHIBA "describes exchanges of energy and water between the atmosphere and the biosphere, and the soil water budget" (Krinner et al., 2005) with a time step of 30 minutes. In IPSL-CM5A and IPSL-CM5A2, SECHIBA includes a two-layer soil hydrology scheme, in which maximum water storage capacity is set globally to 300 kg m-2 over a two-meter parallel coupling library used in both NEMO and LMDz components. It ensures fully parallel interpolation and exchange of 160 the coupling fields. As in IPSL-CM5A, 24 coupling fields are exchanged, 7 from the ocean to the atmosphere and 17 from the atmosphere to the ocean (Valcke, 2013).

XIOS library
The use of XIOS as input-output library allows performing parallel asynchronous input/output by using computing cores as "IO servers". These servers are dedicated to the reading and the writing of the data and permit model computation to continue 165 while data is written or read.

Integration of hybrid parallelization MPI-OpenMP in LMDz atmospheric component
In LMDz, longitudinal filtering is applied on the dynamical equations for latitude higher than 60°in both hemispheres (Hourdin et al., 2013). Thus the choice of MPI domain decomposition is optimal only along latitudinal bands. This decomposition has been initiated in the LMDz version of IPSL-CM5A with a minimum of three latitude bands per MPI process, reduced to two 170 in IPSL-CM5A2. In addition, the use of multi-core processors in HPC led us to add a shared memory parallel processing (OpenMP) in LMDz. The use of this shared memory parallelism on vertical levels in the dynamics allows to increase the number of cores used on the supercomputer and consequently to reduce the elapsed time of the simulation (see next section).

Computing performances
2.3.1 IPSL-CM5A2 vs IPSL-CM5A: computing performances on supercomputer CURIE Table 1. Computing performances metrics collected on CURIE and JOLIOT-CURIE platforms. Parallelization is defined as the total number of computing cores (Number of MPI (x Number of OMP) per component). One core is always used for XIOS process.
(IPSL-CM5A2-S) running on 302 cores. IPSL-CM5A2-T configuration allows to reach 38 SYPD (3.8 times faster than IPSL-CM5A) with a 10.8 % increase in computing cost (CHSY, Table 1 and Fig. 2). IPSL-CM5A2-S reaches 56 SYPD, for a CHSY 42 % larger than for IPSL-CM5A (129,000 hours for 1,000 simulated years). Besides, CHSY is greater for IPSL-CM5A2-S, 190 which means there is a better use of a resource allocation than with IPSL-CM5A. The replacement of the sequential OASIS3.3 coupler by the parallel OASIS3-MCT library in IPSL-CM5A2 also allows reducing the cost of the coupling from 12 to 1 % of the total CPU time of a simulation. Regarding inputs/outputs aspects, IPSL-CM5A used the IOIPSL library, which performed sequential writing and imposed a rebuild post-processing step to merge multiple files into one single file. Instead, IPSL-CM5A2 uses the XIOS library that writes out data asynchronously using dedicated I/O servers. Parallel writing allows users to obtain 195 directly one single file without any time-consuming rebuild step. Still, although the writing is performed asynchronously, the I/O cost remains important in IPSL-CM5A2 (10 %) due to operations (temporal operations, variables combinations) that are now performed by the library (on the client side, i.e. model process) rather than previously into the model.

Performances of IPSL-CM5A2 on the supercomputer JOLIOT-CURIE
Switching from CURIE to JOLIOT-CURIE supercomputer opened opportunities to increase performances of IPSL-CM5A2 200 through several aspects: first, replacing Intel Sandy Bridge by more recent Intel Skylake processors was expected to reduce computation time. Second, we decided to test higher parallelization levels, by combining (i) two-latitude band parallelization (instead of three), (ii) an increased number of OMP processes for LMDz, and (iii) more MPI processes for NEMO. Figure   2 depicts the scaling of IPSL-CM5A2 on JOLIOT-CURIE and CURIE. It shows that for the same parallelization set at 302 cores, IPSL-CM5A2 computes about 43% faster on JOLIOT-CURIE than on CURIE (80 SYPD compared to 56 SYPD), with 205 a cost reduced by ∼30 % (91,000 hours for 1000 year). Finding the most valuable configuration requires both to find a tradeoff between maximizing SYPD and minimizing CHSY for each component and to minimize each component waiting time in coupled configuration. To achieve the latter goal, we used LUCIA (Maisonnave and Caubel, 2014) to determine LMDZ and NEMO computing and waiting time. From the 302-cores configuration, we made a step-by-step increase in LMDZ and NEMO resources that shows that both components waiting times are minimized when LMDZ runs on 47 MPI processes and 8 OMP 210 threads together with NEMO using 60 MPI processes (Fig. A2). This configuration uses 437 cores and allows simulating 98 years per day, leading to a cost of 107,000 hours for a thousand-year simulation. Scaling up does not provide significant improvements in SYPD numbers but the cost of the simulation is dramatically enhanced (39 % increase when we dedicate more than 600 cores to the model (Fig. 2 and Table 1). Consequently, the 437-core configuration has been chosen as the distributed version of reference.

Tuning: target and strategy
We focus here on the general tuning strategy adopted by defining the target and how we reached it, rather than giving a comprehensive description of the numerous simulations that have been designed to obtain the final pre-industrial simulation.
As depicted in Hourdin et al. (2017), several choices need be to made when tuning a model, namely defining a relevant 225 target and associated metrics, defining in what mode the model is to be tuned (e.g. atmospheric only, fully coupled, etc.), and choosing which parameters will be used to reach the target. As most of modeling groups (Hourdin et al., 2017) our choice was to tune IPSL-CM5A2 through pre-industrial runs, using alternatively atmospheric and fully-coupled experiments. A first untuned 1,000-year spin-up of IPSL-CM5A2 forced by CMIP5 pre-industrial boundary conditions provided an adjusted global surface air temperature of 11.3°C. Such a cold bias had been depicted in the previous IPSL-CM5A . We 230 targeted to increase global-mean surface temperature by 2.2°C to reach 13.5°C in pre-industrial conditions with IPSL-CM5A2, expecting this value to translate into 15.5°C in smulations with present-day conditions.
Regarding the free parameters, the choice was to use a parameter which controls the conversion of cloud water to rainfall to alter the cloud radiative forcing (CRF), and thereby the TOA net radiation (Q TOA , counted positive downward) and the globalmean surface temperature. Using clouds parameters for tuning, the most uncertain parameters that affect the most radiation, is 235 also well shared practice among modeling groups. Previous results obtained with LMDZ and IPSL-CM5A show that changing TOA balance by 1 W.m −2 results in a change in temperature of 1 K. It is also the typical value of the sensitivity of global temperature to greenhouse gas concentration. Our +2.2°C target thus translates into an increase of Q TOA by 2.2 W.m −2 . In LMDZ, the conversion of cloud water to rainfall is computed following Sundqvist (1978), as where q lw is the mixing ratio, c lw is the in-cloud water threshold for autoconversion, set at 0.418 g.kg-1 in IPSL-CM5A, For q lw » c lw the time constant for auto-conversion is τ conv (here set at 1800 s) while it goes to infinity (no auto-conversion) for q lw « c lw .

245
Decreasing c lw is expected to lower cloud water content and reduce the CRF. First, we carried out atmospheric-only simulations with varying c lw values to define the sensitivity of Q TOA to c lw . The resulting linear relationship provided a value of 0.325 g.kg −1 for c lw to obtain an increase of 2.2 W.m −2 in Q TOA . Thus we ran a second pre-industrial fully-coupled simulation of 500 years with this tuning. This experiment was interrupted when two of the bugs mentioned above were identified and fixed.
A second set of tuning was required as the bug correction altered the global 2-meter temperature. After several adjustments, the 250 target of +2.2 W.m −2 at TOA was expected to be obtained by setting c lw at 0.343 g.kg −1 . A ultimate fully-coupled experiment (called here PREIND) including this tuning was then branched on the previous simulation and run for 2,800 years.
3.2 Energetic and freshwater balances, equilibrium and drifts  (Hobbs et al., 2016), is small (0.094 W.m −2 over the last 1000 years) compared to IPSL-CM5A (0.18 W.m −2 ), which might be an improvement linked to correc-260 tions of energy conservation between IPSL-CM5A and IPSL-CM5A2 versions of LMDZ (see section 2.1). Figure 3a indicates that this imbalance is constant after ∼300 years of preindustrial simulation, and that OHC stabilizes after ∼2000 years of simulation. This ensures that the model drift is small, therefore validating the use of IPSL-CM5A2 for multi-millennial simulations.
It takes about 500 years for sea surface temperatures and 2-meter air to stabilize after PREIND initialization (Fig. 3b). They 265 converge to 17.5°C and 13.2°C, when globally averaged, respectively. We thus missed the 13.5°C target by 0.3°C which is reasonable when considering the very simple procedure used. Such a difference is also small enough compared to the regional biases in sea surface temperature. An additional 1,500 years is required for the global average ocean temperature to stabilize, as deeper ocean temperatures still cool by ∼0.1°C to 0.2°C during this period (Fig. 3c).   (Hobbs et al., 2016), dashed black lines show the heat content implied by a QTOA of ± 0.5 W.m −2 , the observed late-twentieth-century QTOA (Roemmich et al., 2015). b) Sea-surface temperature (orange) and 2-meter air temperature (blue). c) 1000-m (orange) and 500-m (blue) ocean potential temperature. are characterized by a drift in SSH reduced to 0.01 m / century and a corresponding weaker freshening of the global ocean (4.7x10 −4 psu / century, not shown). Such a drift might be problematic in the future if IPSL-CM5A2 is to be used for very long (transient-like) integrations. Improvements on this concern are expected from the future replacement of LIM2 by the more recent LIM3 version (Rousset et al., 2015), for which conservation is insured. Yearly evolution of long-lived greenhouse gases concentrations and total solar irradiance (TSI) correspond to the historical requirements of the CMIP5 project, with notably reduced TSI to consider the volcanic radiative forcing based on an updated version of Sato et al. (1993). Tropospheric aerosols, tropospheric and stratospheric ozone were prepared as described in (Szopa 285 et al., 2013) and land use changes processed from crop and pasture datasets developed by Hurtt et al. (2011), as described in . Unless otherwise expressed, comparison between the historical run and observations is made by using  LMDZ (and hereby IPSL-CM5A) has been shown to underestimate low and mid-level cloud cover and overestimate high-level clouds at mid and high latitudes (Hourdin et al., 2013), biases shared with many climate models (Zhang et al., 2005). LMDZ cloud representation has been since strongly improved with the new set of physical parameterizations and associated increase 295 in vertical resolution (Hourdin et al., 2013) that is currently being used for CMIP6 exercise. In IPSL-CM5A2, the biases of the standard parameterization translate into CRF, with an overestimation of shortwave CRF at mid and high latitudes (+15 to +25 W.m −2 ), and an underestimation of outgoing longwave trapping in some convective regions of the tropics (-5 to -15 W.m −2 ), when compared to the CERES dataset. The tuning depicted earlier leads to a global decrease in total CRF at all latitudes compared to IPSL-CM5A (Fig. 4).

300
The CRF strong overestimation in the southern hemisphere mid-latitudes is only slightly corrected, as less negative LW CRF 0. 10.
-20.  is likely related to a bad representation of high cloud decks over these regions (Fig. 4d). These biases translate into surface air temperature biases (Fig. 4e). Regions with an underestimated CRF (e.g. the eastern boundaries of the Pacific and the south Atlantic) are associated with up to 2.5°C warm biases. Conversely, regions where CRF is overestimated, like the northern and southern Atlantic and several mountain ranges (the Andes, the Rockies and the Tibetan Plateau) show cold biases (down to -6°C for the Tibetan Plateau).

Temperature and general circulation
The tuning of IPSL-CM5A2 induces a low latitude warming which amplitude varies with height, from +1°C at 850 hPa to +2.5°C at 200 hPa, when compared to IPSL-CM5A, as expected from the adjustment of the moist adiabatic lapse rate (Fig.   5a). As a consequence of the increased the latitudinal temperature gradient between 200 and 300 hPa, the zonal jets expand towards high latitudes, especially in the Northern Hemisphere (Fig. 5b). This is an improvement, as IPSL-CM5A had been 315 shown to underestimate the jets latitudinal extension (Barnes and Polvani, 2013). Despite this improvement, IPSL-CM5A2 still underestimates the jet extension in the Southern Hemisphere, as shown by a negative anomaly in zonal wind in the high latitudes and an overestimation of the zonal winds in the low latitudes. At 850 hPa and surface, IPSL-CM5A2 provides a better representation of the zonal winds in the high latitudes: the negative anomaly with ERA-Interim is reduced along the Antarctic, showing that the westerlies have moved southward. The 850 hPa temperature anomalies also show that tuning has 320 been effective in IPSL-CM5A2 to reduce the cold bias of the model, but also that the warm bias over the eastern tropical Pacific has strengthened. In the tropics, IPSL-CM5A2 depicts biases similar to IPSL-CM5A, i.e. over the Pacific Ocean, surface easterlies bad seasonal cycle lead to an overall yearly misposition, in relation with the double ITCZ, and over the Atlantic, easterlies strength is systematically underestimated. IPSL-CM5A2 also shows overestimated winter westerlies in the northern mid-latitudes, that translate into an anomalous wind stress over the northern Atlantic.

325
Sea-level pressure (slp) and surface wind patterns simulated by IPSL-CM5A2 are compared to ERA-Interim reanalysis in This is particularly visible over the Sahara region in DJF (Fig. 6ab) and the northern part of South America in JJA (Fig. 6de).

Tropical precipitation
Regarding precipitation, IPSL-CM5A2 and IPSLCM5A essentially share the same qualities and biases, i.e. an underestimation 340 of rainfall in the mid-latitudes (20°-40°) of both hemispheres, an overestimation of rainfall in the southern tropics, characterized on the yearly zonal average by a peak at ∼10°S corresponding to a "double ITCZ" 7e. Dry biases (-0.5 to -2 mm.d −1 ) are marked over the western edges of the northern and southern Atlantic at mid-latitudes, as well as the northwestern Pacific in JJA. The regional anomaly between IPSL-CM5A2 and GPCP data also show that the tropical biases vary with longitude: namely rainfall is overestimated over the equatorial western Indian, equatorial Africa and south equatorial Atlantic (+1 to +3 345 mm.d −1 ) in DJF ( Fig. 7ac), whereas there is a strong underestimation of tropical rainfall over South America, specifically over the Amazon basin (down to -5 mm.d −1 ).
South America. The dry bias over South America and every regions of the Amazon basin (Fig. 8) is a long-standing feature of the IPSL coupled models that likely has multiple origins, as precipitation over this region depends both on moisture 350 advection from the Atlantic (and thus from the capability of the model to reproduce correct atmospheric dynamics and SSTs) and continent recycling of water through evapotranspiration.  showed that soil depth had to be increased from 2 to 4 meters in IPSL-CM5A to favor seasonal water retention in the soil, allow vegetation to grow, and thereby reduce a too strong dry bias over the Amazon. We kept this tuning in IPSL-CM5A2. Still, the seasonal 850 hPa wind divergence shows that a strong dipole between restricted convective zones (Fig. A4) and expanded subsiding zones is simulated over tropical 355 South America. This pattern induces restricted areas with strong (>8 mm.d −1 ) convective precipitation (northwestern Amazon in JJA, south-central Amazon in DJF) connected to larger zones where rainfall is limited by subsidence. The increase in the dry bias over South America between IPSL-CM5A and IPSL-CM5A2 results from a feedback loop linked to vegetation change. The tuning of the clouds produces excessive warming (Fig. 4), especially over the northwestern Amazon basin, that Top is DJF, bottom is JJA. e) Zonal average of absolute annually-averaged precipitation (mm/day) for IPSL-CM5A2, IPSL-CM5A and GPCP data.
combines with DJF subsidence to strongly reduce rainfall amounts and thereby tropical rainforest cover (not shown). In line 360 with earlier deforestation experiments carried out with the IPSL model (Davin and de Noblet-Ducoudré, 2010), such a decrease in tree cover has a strong feedback on temperature, as it leads to a decrease of latent heat flux and increase in sensible heat fluxes (overpassing the albedo increase, that tends to cool surface air temperature). As an illustration, the 2-meter temperature anomaly over northwestern Amazonia between IPSL-CM5A2 and IPSL-CM5A reaches + 3.5°C, with IPSL-CM5A2 annual mean two-meter temperature reaching 29.9°C (not shown). Vegetation decrease also further reduces water recycling through 365 lower evapotranspiration as confirmed by a supplemental historical experiment with IPSL-CM5A2 initialized with observed values of vegetation rather than the results from PREIND that shows that prescribing rainforest over the Amazon basin cools two-meter temperature by 1 to 3°C and increase annual rainfall amount by 2 mm.d −1 (not shown). The seasonal cycles of precipitation computed over the northern and southern parts of the Amazon basin show the simulated rainfall collapse during the respective hemispheric winters, while the year-round underestimation of rainfall in the western part of the basin seems to 370 related to a mislocation of rainfall, that occur on top of the Andes instead of the foothills.
Africa. Conversely, precipitation rates are overestimated over equatorial Africa (Fig. 7). The divergence of 850 hPa wind shows that the model overestimates convergence of moist air masses south of the equator in DJF (Fig. A4). Moreover the +0.5 to +2°C anomalies in the eastern tropical Atlantic SST likely favor enhanced evaporation and moisture advection over the continent.
Over central Africa (Fig. 15), rainfall seasonal cycle is bimodal with two humid seasons, but the maxima lag by one month 375 when compared to data. More importantly, rainfall amount are largely overestimated by ∼ + 66 % during the humid seasons.
Over the western African monsoon region, IPSL-CM5A and CM5A2 both lag the start of the rainy season, leading to an overall  underestimation of rainfall before the monsoon peak. However this rainfall maximum occurs in August as for observations, and the subsequent decrease in simulated correctly.
Indian monsoon and Indonesia . Compared to IPSL-CM5A, the Indian monsoon has been improved in IPSL-CM5A2 ( Fig.   380 10ab). The 1-month lag that characterized IPSL-CM5A is still present, but the amplitude of the seasonal cycle of precipitation, although still underestimated, has been enhanced, peaking at 7.3 mm.d −1 in August (5.6 mm.d −1 in August in IPSL-CM5A), closer to observations (9.3 mm.d −1 in July in GPCP data). This improvement is linked both to higher large-scale and convective precipitation during summer, consecutive of the enhanced hydrological cycle linked to the tuning-induced warming. However, the artificially high precipitation rates on the Himalayan foothills have not been corrected. They are made mostly of large-scale 385 precipitation that are overestimated due to a too strong cooling over the highest altitude grid points. Fig. 10c shows Hovmöller diagrams for GPCP dataset compared to IPSL-CM5A2. Although at the first order, the seasonal cycle of precipitation is correctly represented over the "Maritime Continent", IPSL-CM5A2 depicts a strong humid bias over the region (up to +5 mm.d −1 during the humid season). This bias has been described earlier with LMDZ pre-industrial simulations run either in atmospheric-only or fully-coupled mode (Toh et al., 2018;Sarr et al., 2019) and it was shown that LMDZ system-390 atically overestimates rainfall over both land and sea surfaces of the Maritime Continent, together with western Indian Ocean and Bengal Bay. Given the complexity of the Indonesian region, with a myriad of islands that the model spatial resolution cannot capture, it is no surprise that rainfall are misrepresented in this model (Qian, 2008). 12.0m/s Figure 10. Hovmöller diagrams of precipitation (mm.d −1 ) averaged over a 95°E-115°E range. Left, GPCP, right, IPSL-CM5A2.

Ocean mean state and identified biases
The annual SST averaged between 50°S-50°N simulated by IPSL-CM5A2 for the last decades of the historical run is 22.3°C 395 (1°C warmer than IPSL-CM5A), a value in the higher range of CMIP5 ESMs and consistent with observations (Fig. A3).
The warm bias over eastern subtropical Pacific and Atlantic oceans is more marked in IPSL-CM5A2 than in IPSL-CM5A as a response to CRF tuning (Fig. 11ab). This bias has already been reported in several coupled models using NEMO, namely CNRM-ESM1 (Séférian et al., 2016), CNRM-CM5.1 (Voldoire et al., 2013), IPSL-CM4 (Marti et al., 2010) and IPSL-CM5A . As stated in Voldoire et al. (2013), ocean-only experiments (Griffies et al., 2009) have suggested that 400 such a warm bias was likely linked to "poorly resolved coastal upwellings and underestimated associated westward mass transport due to the coarse model grid resolution". The 0.5-2°C warm bias over the Antarctic circumpolar current (ACC) is comparable to IPSL-CM5A and can be related to the poor simulation of the CRF over this region (Hyder et al., 2018) and the associated mispositioning of mid-latitude westerlies depicted above. Conversely, the tuning has reduced the cold bias in the Pacific and Atlantic gyres. Surface salinity biases are rather similar between both models. The surface salinity anomaly ocean freshwater balance. Both models depict too fresh surface waters in the northern Atlantic and locations of deep water formation (Fig. 11cd). The ocean also depicts temperature biases at depth both in the Atlantic and Pacific oceans (Fig. 12).
Zonally-averaged potential temperature anomalies show that the warm bias in the northern Atlantic between 1000 and 3500 m has been amplified, as well as the warm bias in the northern Pacific, that concerns the entire water column in IPSL-CM5A2.

410
The strong (> + 2.5°C) overestimation of temperatures in the subsurface waters of the southern hemisphere is almost identical in both versions of the model. The cold (∼ − 1°C) bias in the high latitudes of the southern hemisphere has been slightly reduced, but is still present due to the lack of strong enough overturning in this region. Sea-ice extent (Fig. 13)

Ocean Heat Transport and Meridional Overturning Circulation
IPSL-CM5A2 ocean meridional heat transport is higher (in absolute value) than the one simulated with IPSL-CM5A in both hemispheres (Fig. 14). This increase in the Northern Hemisphere is mainly related to the increase in the Atlantic meridional heat 420 transport in the mid-latitudes (not shown). When compared to observations, the two model versions have asymmetric biases, with the northward transport being underestimated in the Northern Hemisphere and the southward transport overestimated in the Southern Hemisphere low latitudes (between the equator and 35°S). Southward heat transport is too strong between the equator and 40°S and is in better agreement with data than IPSL-CM5A between 40°S and 60°S. The Atlantic Meridional Overturning Circulation (AMOC) is more vigorous in IPSL-CM5A2, with a maximum in depth (below 500 m) and latitude 425 (from 30°S to 60°N) in preindustrial conditions moving from 10.3 ± 1.2 Sv (Escudier et al., 2013) to 12.02 ± 1.1 Sv (taken for the last 1000 years of simulation, Fig. 15a). The last 20 years of the 20th century in the historical run depicts a less vigorous circulation in the Atlantic, with a mean value for AMOC index of 10.5 ± 0.9 Sv. The convection sites have been slightly modified between the two versions of the model. While the strongest convection was located in the North Atlantic subpolar gyre (SPG), south of Iceland, in the IPSL-CM5A version (Escudier et al., 2013), in IPSL-CM5A2, the Greenland sea shows the 430 deepest mixed layer depth, which is larger than 1500 meters in winter (Fig. 13cd). Nevertheless, there is still a very active deep convection site south of Iceland, which is not realistic. The strengthening of the Nordic Seas convection can be related with the decrease of its sea-ice cover in winter (Fig. 13). As a consequence, denser water is now produced in this site, which improves the overflow through Greenland-Iceland-Scotland straits. Indeed, the water denser than 27.8 kg.m −3 flowing southward reaches 1.5 Sv in IPSL-CM5A2, while it was only 0.2 Sv in IPSL-CM5A. It remains largely underestimated as the observations from 435 (Olsen et al., 2008) indicates around 6 Sv of overflow in total from observation-based estimates over the last few decades.
Besides, the poleward shift of the westerlies in the North Atlantic stimulate a northward shift of the boundary between the Atlantic subtropical and subpolar gyres, as shown by the barotropic streamfunction anomaly (Fig. 15b). This may contribute to bringing more salt into the subpolar gyre, which may also act to intensify the AMOC in IPSL-CM5A2.   440 We computed oceanic mass transport through selected longitudinal and latitudinal cross-sections corresponding to the major gateways (Table 3) that uses NEMO as well (Voldoire et al., 2013). These diagnostics highlight that the Antarctica Circumpolar Current (ACC) is underestimated by IPSL-CM5A2 by around 25%, as for the other configuration considered here, whereas despite the absence of high resolution and explicit resolution of Indonesian gateways, the overall Indonesian throughflow tends to be close to 445 observation-based estimate (underestimation of 16%). Compared to observations and CNRM-CM5.1, the flow in the Florida strait is quite weak (weaker than observation by 37%), in line with the simulated "sluggish" AMOC. The increase of the AMOC from IPSL-CM5A to IPSL-CM5A2 has not been enough to strengthen this current, also largely resulting from wind forcing.

Fluxes in the major ocean gateways
Also, the absence of resolved eddies can also play a significant role in this underestimation.

450
To illustrate some of the capabilities of IPSL-CM5A2, we also evaluate the simulated oceanic NPP (Fig. 16). The global integrated NPP over the historical period (here averaged over 1980-1999 from the historical simulation) is 47.5 PgC.y −1 , which is a strong increase when compared to 30.9 PgC.y −1 simulated for IPSL-CM5A , making IPSL-CM5A2 being closer to the 52.1 PgC.y −1 global estimate based on SeaWifs remote-sensing observations (Behrenfeld and Falkowski, 1997). The representation of NPP at the spatial scale compares well to observation-based estimates for the tropical 455 and Southern Hemisphere, but depicts large biases in the Northern Hemisphere, especially in the North Atlantic where NPP is largely under-estimated (Fig. 16c).

Ocean variability
The main modes of variability of the AMOC are analyzed with an empirical orthogonal function (EOF) of the yearly Atlantic meridional overturning streamfunction (Fig. 17, top). The first EOF is a monopole in IPSL-CM5A and IPSL-CM5A2, but 460 IPSL-CM5A2 shows less loading in the subpolar gyre between 40°N and 60°N. The associated principal component shows a significant 20-yr variability in the case of IPSL-CM5A. Such variability was previously linked to advection of surface salinity anomalies from the western boundary to the Eastern Atlantic (Escudier et al., 2013), and westward propagating planetary waves in subsurface (Ortega et al., 2015). However, the peak of variability at 20-yr timescale is not significant any more in IPSL-CM5A2 (see spectra in Fig. 17, center). This difference can be explained by the different tuning of the two model ver- in IPSL-CM5A-MR. Such a warmer mean state modifies the subpolar gyre currents and the subsurface stratification in the Eastern Atlantic Ocean where the mixed layer is deepest in IPSL-CM5A models. These changes were suggested to explain 470 a smaller growth of baroclinic instabilities triggering the propagation of westward propagating anomalies (Gastineau et al., 2018). It is therefore likely that the different tuning of IPSL-CM5A2 led to similar changes. The SST associated with the first principal component of the AMOC also shows weaker SST anomalies in the subpolar gyre when compared to IPSL-CM5A (not shown). The Atlantic multidecadal variability (AMV), was calculated from the low-pass filtered (61-month running mean) North Atlantic SST in 0°N-60°N 80°W-0°E. The AMV pattern is given by the regression of the SST onto the AMV index.

475
The AMV patterns are simimar in IPSL-CM5A and IPSL-CM5A2, with a slightly more intense subpolar gyre SST (6.8°C vs 6.2°C for the maximum value in the subpolar gyre) in IPSL-CM5A2. Larger positive SST anomalies in IPSL-CM5A2 are also simulated in the Nordic Seas and Irminger sea associated with the AMV. The larger AMV in IPSL-CM5A2 is consistent with the smaller sea-ice cover in the North Atlantic, and a more intense signature of the atmospheric forcing. American pattern, as in other CMIP5 models (Weare, 2013). The Pacific Decadal Oscillation (PDO) was also calculated as the first EOF of the monthly Pacific SST anomalies north of 20°N (Fig. A7), which illustrate a pattern with positive SST anomalies over the equatorial Pacific, negative anomalies in the northwest Pacific, and horseshoe pattern of positive SST anomalies in the Northern Hemisphere. The IPSL-CM5A model generally simulates the characteristics of the observed PDO, as found by (Fleming and Anchukaitis, 2016) or (Nidheesh et al., 2017), as well as the associated anomalies in the Atlantic and ocean 495 basins. The IPSL-CM5A2 version simulates a PDO nearly identical to that of IPSL-CM5A, so that the mean state changes have no impact on the PDO state. simulation of the Cretaceous as a case study. We present a very brief analysis of the model timeseries to illustrate the long adjustment required to reach equilibrium in deep-time simulations. As the focus of the paper is the description of IPSL-CM5A2, we only provide a very basic description of the simulated climate, whose characterization will be the focus of a subsequent study.

Generating a new grid for NEMO
In addition to the aforementioned computational performances, technical issues related to the oceanic component of the model also have long prevented the use of any version of the IPSL ESM for deep-time climate modelling. The main lock was that NEMO runs on a tripolar grid (ORCA2, with two poles in the northern hemisphere and one over Antarctica) that was designed, amongst other objectives, to avoid any singularity points on the computational domain (i.e. the ocean domain) (Madec, 2015; 510 Madec and Imbard, 1996). Specifically, north of 20°N the mesh is made of a series of embedded ellipses, the foci of which are the two north poles of the grid that are by design placed over continents (Fig. 18a). Southward the grid is a Mercator (i.e. identical longitudinal ant latitudinal grid spacing) grid, defined down to 78°S, except from local refinements applied in the tropics, so that the meridional resolution reaches 0.5°at the Equator and over the Red, Black, Caspian and Mediterranean Seas where resolution is about 1°.

515
Because the grid poles are located over North America, Asia and Antarctica, changing the land-sea mask to account for continental drift can shift the singularity points into the computation domain. In addition, a number of ocean routines include hard-coded specifications that apply to some narrow modern straits, such as the Gibraltar strait. In these regions, the resolution of the grid is much larger than the width of the straits, thereby requiring the outflow of water to be prescribed. Although needed to better simulate modern climates, these features are irrelevant for deep-time simulations. To overcome these issues, we design a new PALEORCA grid based on techniques and code designed by G. Madec and A. Sellar. The southern hemisphere (SH) pole is conserved, which restricts the application of IPSL-CM5A2 to the last 100 million-years (Antarctica has been located over the South Pole since ca. mid-Cretaceous). In the northern hemisphere (NH), new ellipses are computed, changing the location of the two singularity points to eventually obtain a grid that can effectively be used for the last 100 million-years. To this end, and because the grid structure requires both NH poles to be diametrically opposed and limits the range of latitude 525 in which the location of the poles can be chosen, the two poles are placed respectively at 106°E, 60°N and 74°W, 60°N (Fig.   18b). Above-mentioned refinements over the Mediterranean and Eurasian Seas are removed, but the tropical refinement is kept. We also take advantage of the well-defined coding guidelines of NEMO to create a "PALEORCA" configuration within the available oceanic configurations of the model. In this configuration, no hard-coded water exchanges across the globe are prescribed.

Application: Generating boundary conditions for deep-time climate modelling
The Cretaceous is historically known for its warmer than present climatic state without permanent ice sheets at the poles. Here, we present a simulation of the Turonian stage (circa 90 Ma), during which the apex of the Cretaceous warmth was reached.
CO2 estimates reconstructed from proxy-based studies vary from 500 ppm to > 2000 ppm with large uncertainties (Wang et al., 2014). In our simulation, we impose an atmospheric CO2 concentration of 1120 ppm (4x preindustrial values) and we remove 535 polar ice sheets. In addition, we keep the modern orbital configuration as no numerical solution for orbital parameters exists for the Cretaceous (Laskar et al., 2004). The basis of our Turonian configuration is the Cenomanian/Turonian paleogeography of (Sewall et al., 2007).

Oceanic boundary conditions
Bathymetry. The Turonian bathymetry (Fig. 19b) is obtained from the merging of two datasets. First, the regular longitude-540 latitude Cenomanian/Turonian paleogeography from (Sewall et al., 2007) is masked over continental areas and remapped onto the curvilinear PALEORCA grid. Given that the bathymetry of the Pacific and Indian oceans is mostly flat (without ridges) and deep (∼5000 m) in Sewall et al. (2007), we have decided to merge the more recent bathymetric dataset of Muller et al. (2008) with the land/sea mask and continental margins of Sewall et al. (2007) to create a more realistic deep bathymetry. Handmade corrections are then applied to the bathymetry to ensure that narrow straits (e.g. West Africa) are at least two grid points wide.

545
This avoids the implementation of specific hard-coded schemes in the model to represent water advection across those straits.
Geothermal heating. The model modern boundary conditions include a realistic map of geothermal heating, whose impact on the ocean abyssal dynamics is non-negligible (Emile-Geay and . As the geothermal heat flow q(t) can be related to the ocean lithosphere age t, we construct a geothermal heating map for the Turonian based on the reconstruction of the Turonian ocean crustal age of Müller et al. (2008). We derive the geothermal heatflux using the age-heatflux relationship of the For t ≤ 55 Ma: For t > 55 Ma: In regions where crust age is missing in (Muller et al., 2008) dataset, we impose the minimal geothermal heat flux obtained with the above equations. The geothermal heating map created is then remapped on the PALEORCA grid and we limit the maximal heat flux to 400 mW.m −2 because the Stein and Stein (1992)  dissipation values for the Cretaceous, we neglect this forcing and set these fields to zero.
Initial conditions.Salinity is initialized as a uniform value of 34.7 per mil equal to the modern mean ocean value . The ocean initial temperature follows a slightly altered formulation of the zonally symmetric temperature distribution These initialization temperatures are quite high, as the ice-free Cretaceous periods have been shown to have far much warmer oceans than present, associated with high atmospheric pCO2. No sea ice is prescribed at the beginning of the simulation.

Continental boundary conditions
Topography We use the Cenomanian/Turonian topography from Sewall et al. (2007) but modify the Antarctic continent because the PALEORCA mesh is not defined beyond 78°S, which corresponds to the maximal extent of the modern Antarctic 575 coastline (including ice shelves). This requires that the Antarctic coastline reaches at least the same latitude for paleoclimate simulations, which is not the case in the Cenomanian/Turonian continental configuration of Sewall et al. (2007). We therefore impose the maximal Antarctic topography proposed by ? for the Eocene-Oligocene, which has the advantages of meeting the required Antarctic coastline criterion and is a relatively recent estimate of a deep-time Antarctic topography, which remains generally poorly constrained. The final topography/bathymetry is shown on Fig. (19b). LMDZ requires a 10-minute global 580 elevation dataset to compute subgrid-scale parameters (elevation standard deviation, slope, peaks and valleys) that are needed to compute surface roughness length as well as to activate orographic drag parameterizations (Lott and Miller, 1997). As the paleotopography provided by Sewall et al. (2007) is 2.8°x1.4°in resolution, we simply artificially upscaled it to a 10-minute resolution through the use of cdo remapping tools (Schulzweida, 2019).
River routing The river routing system (i.e. length of rivers, direction of flows, locations of river outflows) is constructed by 585 first upscaling the low-resolution topography of Sewall et al. (2007) to a 1°x1°resolution and second using the (Wang and Liu, 2006) algorithm within the QGIS software (QGIS Geographic Information System. Open Source Geospatial Foundation Project. ). This algorithm ensures that inland depressions are filled and therefore that all continental freshwater is routed to the ocean in the coupled model.

Vegetation
In the absence of globally consistent vegetation reconstructions for the Turonian, the imposed vegetation follows a 590 very simple latitudinal distribution of Plant Functional Types (4).

5.4
Results of the case-study experiment

Equilibrium
Time series of SST and T2M shows the strong model departure from the initial conditions in the first centuries of the simulation.
After 3,000 years, SST and T2M are adjusted at 28.3°C and 26°C, respectively (Fig. 20e). OHC of the first 300 meters is 595 stabilized as well. The long-term adjustment of the ocean is, as for our initial pre-industrial experiment, driven by the initial radiative unbalance resulting from the initial ocean temperatures and the greenhouse forcing. The vertical cross sections at initial and final states of the experiment shows the spatial patterns of adjustment of ocean temperature and salinity (Fig.   20abcd). The constant salinity together with the latitudinal profile of temperature imposed as initial conditions also induce a latitudinal gradient in seawater density, and thereby an adjustment of ocean dynamics. Surface salinity and temperature are 600 also adjusting driven by the precipitation minus evaporation balance (PE) and freshwater inputs from river routing (see next section) + radiative balance. Contrary to surface signals, the total OHC shows that the deep ocean is still adjusting after 3,000 years ( Fig. 20e), showing the need for very long integrations when one model deeptime paleoclimates.

Mean surface climate
The elevated values of SST and T2M at equilibrium are primarily explained by the high atmospheric CO2 concentration 605 prescribed. In the ocean, warmest temperatures (> 36°C) are found in the equatorial band in the Atlantic and Indian oceans as well as in the western Pacific (Fig. 21ab), whereas the eastern equatorial Pacific is cooler by a few degrees. The Southern Ocean is relatively warm (> 10°C) as is the North Pacific while the Arctic Ocean exhibits the coldest surface temperatures (reaching less than 8°C at the highest latitudes). There is no sea ice formed in the polar oceans. In continental areas, the radiative impacts of the high atmospheric CO2 concentration and the absence of ice sheets also substantially affect the mean annual t2m. Values 610 > 40°C are found in equatorial Africa, South America and southeastern Asia, whereas below 0°C mean annual continental t2m are only located in Antarctica (Fig. 21b). The winter season yields well below freezing t2m over large parts of the NH high latitudes, with lowest values in northern Siberia. Similarly, the whole Antarctic continent suffers below freezing t2m, with values below -16°C over large parts of the continental interior). In the summer however, Antarctic t2m are mostly above 20°C and reach up to more than 28°C close to the modern Weddell Sea (not shown). Tropical continents are affected by 615 extremely high summer t2m reaching more then 44°C whereas NH high latitudes exhibit the lowest summer values with < 16°C in northern Siberia.
The mean annual precipitation-rainfall (PE) shows expected patterns across the globe (Fig. 21d) This simulation, and by extant the overall response of IPSL-CM5A2 to a Cretaceous paleogeography without ice sheets and elevated CO2 levels, would need to be evaluated in more details against previous work and proxy data. Still, the basic analysis presented here demonstrates that the mean climate simulated by our model shares many similarities with that simulated by other investigations of Cretaceous warmth with coupled models (Niezgodzki et al., 2017;Tabor et al., 2016;Upchurch et al., 2015). It therefore suggests that our goal of adapting the IPSL ESM to deep-time boundary conditions has been reached with 640 success.  large scale biases when forced by four-time atmospheric CO2 concentration. Answering the same question for paleoclimate contexts in which the boundary conditions, such as land-sea mask, topography and bathymetry, insolation and greenhouse gas forcing, are so different than present that they impose very strong forcings on the climate system, remains to be done. Until this challenge is tackled, the use of a computing-efficient ESM to run sensitivity experiments to individual forcing factors, for which the geologic record provides good constraints, appears to be the key.

660
On the technical point of view, the ongoing race to higher resolution in ESMs asks for consistency between code architecture and supercomputer design. For example, a simple comparison of computing performances of IPSL-CM5A2 and IPSL-CM6-LR, the higher resolution IPSL model designed for CMIP6, shows that the latter requires ∼8 times more computing resource than IPSL-CM5A2 to run a simulation of a given length (Fig. A10). Specifically, for a 1000-year experiment, IPSL-CM6 require more than 1.6 millions of computing hours, simulating 16 SYPD on more than 1000 cores. This model is being used 665 for paleoclimate snapshots experiments, like in the PMIP project, but its use for long, multi-millennial simulations that are typical of deep-time paleoclimate studies is compromised by these performances and requirements. Similar constraints have led colleagues from other climate modelling group to keep updating and using climate models at the resolution they were originally designed more than 15 years ago (Valdes et al., 2017). Given the apparent contradiction between ( of IPSL ESMs, with the inclusion of passive tracers in the ocean (Arsouze et al., 2007), the water isotopes in the atmosphere (Risi et al., 2010), and the explicit resolution of chemical photoreactions in the stratosphere (Jourdain et al., 2008). All of them appear to be particularly relevant for paleoclimate studies (Botsyun et al., 2019;Sepulchre et al., 2014;Szopa et al., 2019), but most drastically reduce the computing efficiency of the model. A trade-off between the number of tracers and/or processes to be included in the model, its spatial resolution, and the length of simulations required for paleoclimate simulations will thus 680 have to be found in the next years to setup relevant experiments designs in the deep-time climate modelling field.  (Crameri, 2018a, b), to prevent visual distortion of the data. We also   Competing interests. The authors declare no competing interests.
Disclaimer. The simulations analyzed in this paper required approximately 500,000 CPU hours, emitting approximately 1.7 teqCO2.