The purpose of this paper is to describe all components involved in CNRM‐CM6‐1 including the coupler setup and output server implementation (section 2 ). Section 3 details the calibration method followed, the forcing data set used to meet CMIP6 requirements, and the resulting drift obtained under constant preindustrial forcing. Section 4 assesses the mean state of the simulated climate, while section 5 addresses its variability. Finally, the model climate sensitivity is discussed in section 6 . In sections 4 to 6 , CNRM‐CM6‐1 is compared to CNRM‐CM5.1 when applicable, to provide a comprehensive analysis of progresses made since the CMIP5 exercise. The last section summarizes the main skills and drawbacks of CNRM‐CM6‐1 and provides some perspectives. This paper aims to provide to the scientific community a general picture of CNRM‐CM6‐1, which is required to interpret results obtained within the multi‐model framework encouraged by the CMIP panel. Additionally, CNRM‐CM6‐1 constitutes the core of the CNRM Earth System Model CNRM‐ESM 2‐1 designed to perform Earth system oriented Model Intercomparison Projects (MIPs) as part of CMIP6. It is also the base for a higher resolution version of CNRM‐CM, namely, CNRM‐CM6‐1‐HR.

As a major update of CNRM‐CM5.1, CNRM‐CM6‐1 is expected to better represent the climate system, both in terms of mean state and variability. It is intended that atmospheric processes are better represented by this drastic change in atmospheric parameterizations; however, this may not result in a general improvement of the model as experienced by other models, as shown, for example, by the version 5B of the IPSL‐CM model (Hourdin et al., 2013 )⁠. In this model, the change in atmospheric physics parameterizations led to an improvement of the tropical climate but a worsening of midlatitude and oceanic circulations (Dufresne et al., 2013 )⁠. CNRM‐CM5.1 takes benefit from a long history, which led to a fine tuning of the model, but also to several known or unknown compensations of errors. The main challenge in the course of CNRM‐CM6‐1 development was thus to achieve a satisfying model calibration in a rather short amount of time. This required work was, however, an opportunity to address this process in a more objective and comprehensive way, which will be detailed hereafter.

In comparison to its predecessor CNRM‐CM5.1 (Voldoire et al., 2013 )⁠, which was the basis for the CNRM‐Cerfacs contribution to CMIP5, CNRM‐CM6‐1 (Figure 1 ) corresponds to a major update of its atmospheric and land surface components. It now embarks ARPEGE‐Climat version 6.3 and SURFEX version 8.0, with fully revised physical parameterizations for the atmosphere ( http://www.umr‐cnrm.fr/cmip6/references ) and the land surface (Decharme et al., 2019 ), respectively. All other model components were also updated to their newest respective version but with only limited changes in the choice of parameterizations involved. In CNRM‐CM6‐1, the ocean, sea ice, and river routing components are based respectively on NEMO version 3.6 (Madec et al., 2017 ), GELATO version 6, and a new TRIP version redeveloped at CNRM and called CTRIP (Decharme et al., 2019 ). The coupling software was also updated to OASIS3‐MCT (Craig et al., 2017 ), which now allows for parallel exchanges of fields. A novelty of the system is the inclusion of an output manager, the XIOS server (Meurdesoif, 2018 ), which allows efficient online output processing.

CNRM and Cerfacs collaborate since the 1990s to develop state‐of‐the‐art atmosphere‐ocean coupled general circulation models for climate studies: the CNRM‐CM suite. Since its origin, this model is based on the coupling via OASIS (Valcke et al., 2006 ) of ARPEGE‐Climat (Déqué et al., 1994 ) atmospheric model, which included the ISBA scheme (Noilhan & Planton, 1989 )⁠ for the land surface, with the OPA ocean model (Madec et al., 1998 ). This model was progressively enriched to include other climate components: the sea ice model GELATO (Salas Mélia, 2002 ), the externalized surface interface SURFEX (Masson et al., 2013 )⁠, and the river routing scheme TRIP (Total Runoff Integrating Pathways, Oki & Sud, 1998 ). In preparation of the sixth phase of the Coupled Model Intercomparison Project (CMIP6, Eyring et al., 2016 )⁠, the group has released an updated version, namely CNRM‐CM6‐1.

This approach is highly efficient compared to the classical offline post‐processing approach, especially in the case of CMIP6 Data Request (DR), which combines 248 experiments. The DR involves 1,274 geophysical variables, which can be invoked in 44 so‐called MIP tables, thus leading to 2027 “CMOR variables” covering a variety of spatial shapes and frequencies. An additional tool, called dr2xml ( https://github.com/senesis/dr2pub ), translates the CMIP6 DR for each year of each experiment in a set of XIOS conformant XML definitions. These definitions are then used to activate the outputs in CNRM‐CM6‐1 components, which are all XIOS‐enabled thanks to an alias table mapping model variable names and CMIP6 DR variable names. At run time, CNRM‐CM6‐1 then creates NetCDF output files, which are directly conformant with CMIP6 requirements and ready to publish on the Earth System Federated Grid (ESGF), except for some climatologies. This approach is applicable to any other model as long as it uses XIOS to produce its outputs. It is currently used in the IPSL climate model for its CMIP6 production.

XIOS postprocessing operations include time sampling and averaging, spatial remapping and reduction, vertical interpolation, and simple arithmetics. They are configured at run time using a XML syntax. The processing is actually shared between the model MPI tasks and a number of XIOS additional MPI tasks, called servers. Two levels of servers are used to aggregate and redistribute the parallel output fields so that any given field is gathered into one file written to disk by one single MPI task (thus avoiding parallel writing). Special care was taken in XIOS to overlap model computation and input/output operations using communication buffers and thus to allow for scaling at high number of cores.

Once the distribution of the full model grid among processes is described, the XIOS library has a main interface routine (send_field), which allows each model MPI process to deliver easily its part of the full grid for any field. The routine send_field can be called from any model routine.

CNRM‐CM6‐1 is used in a large number of MIPs, which request numerous diagnostics to be output following specific format norms. In order to ease their production, the model was interfaced with XIOS (Meurdesoif, 2018 ), an input/output parallel server software allowing for a declarative description of output files content and for the realization of online operations on fields, thereby almost fully removing the need for postprocessing.

As the SURFEX tiling approach considers land, ocean, and lake surfaces to calculate fluxes, NEMO only receives the flux component computed over the ocean fraction of the grid cell. Due to ocean coastline mismatch between NEMO and SURFEX, a global conservation procedure is applied at each coupling time step to ensure energy and water conservation in the coupled system.

The energy associated to water fluxes is taken into account in NEMO/GELATO. Liquid precipitation, evaporation, and river discharges enthalpy is calculated considering that they are given at sea surface temperature. Similarly, ice sheets water flux and snow falls are considered to be solid, with a 0 °C temperature; thus, the enthalpy flux corresponding to the melting of ice/snow is considered.

NEMO sends the ocean and sea ice surface properties (sea surface temperature and currents, sea ice fraction and sea ice albedo) that are used in SURFEX to calculate surface turbulent and radiative fluxes, which are then sent back to NEMO. The same is done for the water cycle; SURFEX sends all runoffs and all floodwater relative fluxes to CTRIP. In turn, CTRIP returns back to SURFEX the groundwater table depth and its grid cell coverage and routes continental water runoff to the ocean. When entering NEMO, river discharges are spread vertically depending on the flow intensity to a maximum depth of 10 m. Additionally, SURFEX provides a net water budget over continental ice sheets and lakes to NEMO to close the water budget. This water from the Greenland ice sheet and lakes is incorporated at the ocean surface globally, whereas the Antarctic ice sheet water excess is only spread south of 60S.

The model components are coupled using the OASIS3‐MCT software (Craig et al., 2017 ), which implements exchanges between multiple executables running concurrently using Message Passing Interface (MPI) communication. The OASIS3‐MCT approach has the advantage of requiring a minimal amount of modifications in existing component codes. In CNRM‐CM6‐1, OASIS3‐MCT transfers and interpolates coupling fields between SURFEX, CTRIP, and NEMO at a coupling frequency of 1 hr. In contrast ARPEGE‐Climat and SURFEX models are coupled inline, that is, SURFEX is called as a subroutine of ARPEGE‐Climat at each atmospheric time step, and are thus considered as one executable by OASIS3‐MCT. Similarly, NEMO and GELATO are also coupled inline at a 1‐hr frequency.

Some aspects of the ice‐ocean interface were revised. In the quadratic bulk formula used to calculate the ice‐ocean stress for ice‐ocean momentum exchanges, the drag coefficient is set to 0.001. This value is higher to that used in CNRM‐CM5.1 and other models (e.g., Vancoppenolle et al., 2009 ), in order to take into account that the top ocean level now resides within the relatively thin surface layer (~1‐ to 3‐m thick), following Roy et al. ( 2015 ). The chosen value is a compromise between this rationale and stability constraints. The oceanic heat flux at ice bottom derives from McPhee ( 1992 ).

The elastic‐viscous‐plastic sea ice rheology implemented by Bouillon et al. ( 2009 ) on an Arakawa C‐grid is used, consistently with the formulation of NEMO finite difference scheme, whereas it was solved on an Arakawa B‐grid in Gelato 5. The sea ice transport follows an incremental remapping scheme based on an Arakawa‐B formulation by Hunke et al. ( 1997 ). It was adapted to an Arakawa C‐grid according to previous work from M. Bentsen to ensure a better consistency between the dynamics and transport, which is crucial to represent sea ice transport through one‐grid‐cell‐wide straits. The transported variables are snow density, volume and enthalpy, and ice surface, volume, enthalpy, salinity, and age. Depending on sea ice thickness, rafting, and ridging may take place in case of sea ice convergence, as described in Salas Mélia ( 2002 ).

Within the standard configuration of CNRM‐CM6‐1, Gelato has a time step of 3,600 s. It is used with five ice thickness categories, based on the World Meteorological Organization classification (less than 0.30, 0.3–0.7, 0.7–1.2, 1.2–2, and over 2‐m thick). In each ice category, the snow is treated with a single layer, while ice is simulated with a nine‐layer vertical discretization. The sea ice enthalpy formulation is based on Notz ( 2005 ). As in Hunke and Lipscomb ( 2010 ), an iterative method is used to solve the vertical heat diffusion equation in sea ice. The solar radiative transfert scheme through the snow pack covering sea ice was upgraded following Grenfell and Maykut ( 1977 ). Albedo of dry snow, melting snow, and melting ice are model parameters, set to 0.88, 0.77, and 0.58 in CNRM‐CM6‐1, respectively.

Sea ice within CNRM‐CM6‐1 is represented by Gelato 6. Most upgrades from version 5 of the code (Chevallier et al., 2013 ; Voldoire et al., 2013 ) aimed at improving the overall consistency of the code (in particular, salt, water, and energy conservation) and increasing its computational efficiency. In contrast with Gelato 5, Gelato 6 is a fully parallel code. It is embedded within NEMO, has the same horizontal grid, and inherits the global domain decomposition of the ocean code.

As for the lateral physics, the isoneutral tracer diffusivity and the horizontal viscosity are parameterized as in NOCS‐ORCA1 (Danabasoglu et al., 2014 ). Parameterization of vertical mixing is also similar, with the addition of the internal wave‐induced mixing parameterization of de Lavergne et al. ( 2016 ) which uses the tidal dissipation climatology of de Lavergne et al. ( 2019 ). Finally, the eddy‐induced circulation is parameterized with the Gent and Mcwilliams ( 1990 ) parameterization as in NOCS‐ORCA1, and with the addition of the Fox‐Kemper et al. ( 2011 ) submesoscale mixed layer eddy overturning scheme.

Radiative transfer in the water column is resolved using a chlorophyll‐dependent three‐waveband scheme as described in Lengaigne et al. ( 2007 ) and Mignot et al. ( 2013 ), using a seasonal climatology of surface chlorophyll concentration derived from a former 60‐year‐long simulation run with NEMO‐PISCES (e.g., Lee et al., 2016 ). A vertical profile of chlorophyll concentration is extrapolated from surface concentrations (Morel & Berthon, 1989 ).

CNRM‐CM6‐1 resolves ocean dynamics on 75 vertical levels using a vertical z coordinate with partial step bathymetry formulation (Barnier et al., 2009 ). The model layer thickness increases from 1 m near the surface to 200 m at a depth of 6,000 m. The time step is 30 min. At the surface, the model uses the split‐explicit nonlinear free surface formulation proposed by Shchepetkin and McWilliams ( 2005 ), with a variable volume. Seawater thermodynamics uses the Roquet et al. ( 2015 ) polynomial approximation of TEOS‐10 (IOC et al., 2010 ). Therefore, Conservative Temperature and Absolute Salinity are the model prognostic variables.

In CNRM‐CM6‐1, NEMO is run on eORCA1 horizontal grid, which is an extension of the ORCA1° tripolar grid already used in CNRM‐CM5.1. The eORCA grids family differs from the ORCA grids family by the use of two quasi‐isotropic bipolar grids south of 67°S instead of the former Mercator grid. The eORCA grids family allows for a more realistic representation of the contours of Antarctic ice shelves (Mathiot et al., 2017 ). In eORCA1, a nominal resolution of 1° is chosen with a latitudinal grid refinement of 1/3° in the tropics.

The ocean component of CNRM‐CM6‐1 is based on the version 3.6 of NEMO (Nucleus for European Models of the Ocean; Madec et al., 2017 ). It is based on the NOCS‐ORCA1 configuration described in details in Danabasoglu et al. ( 2014 ). This configuration has been shared among climate modeling group through a dedicated NEMO project called “shaconemo” lead by Institut Pierre Simon Laplace (IPSL).

Over the ocean, SURFEX resolves the exchange of momentum, energy, and water across the air‐sea interface. The radiative properties of the seawater are handled by the ocean surface albedo scheme proposed by Séférian et al. ( 2018 ). The turbulent fluxes of momentum, heat, and water are computed using an improved version of the Exchange Coefficients from Unified Multi‐campaigns Estimates (ECUME) scheme. It is a bulk iterative parameterization developed at CNRM from in situ measurements. This ECUME formulation is considered to better replicate the variability of the observed turbulent fluxes across a wide range of atmospheric and oceanic conditions.

Lakes, including both the Caspian and the Aral seas, are represented using the bulk FLake model, which computes the temporal evolution of the vertical lake temperature profile from the surface mixing layer to the bottom. Modeled lakes can be covered with ice and snow. A skin temperature of a 1‐mm thickness was introduced to simulate a surface temperature representative of the energy budget at the lake surface. The spatial distribution of lakes at the global scale is given by the ECOCLIMAP‐II database, while the lake depth was specified from the 1‐km Global Lake Depth database (Kourzeneva, 2010 ). More details can be found in Le Moigne et al. ( 2016 ). Note that FLake does not simulate the water mass evolution. The associated imbalance between precipitation, runoff, and evaporation over lake is thus artificially spread uniformly over the ocean.

The land surface in one grid cell is tiled into 12 patches in order to account for the variety of soil and vegetation behaviors within a grid point. They aggregate the 500 land cover units at 1‐km resolution present in the ECOCLIMAP‐II database (Faroux et al., 2013 ). Mean seasonal cycles of both the snow‐free albedo and the leaf area index are prescribed from Moderate Resolution Imaging Spectroradiometer products at 1‐km spatial resolution. The soil textural properties (clay, sand, and soil organic carbon content) are given by the Harmonized World Soil Database ( http://webarchive.iiasa.ac.at/Research/LUC/External‐World‐soil‐database/HTML/ ) at a 1‐km resolution. Finally, the mean topography is derived from the 1‐km Global Multi‐resolution Terrain Elevation Data 2010 ( https://topotools.cr.usgs.gov/gmted_viewer/ ). More details can be found in Decharme et al. ( 2019 ).

The land surface is represented using the new Interaction Soil‐Biosphere‐Atmosphere‐CNRM TRIP (ISBA‐CTRIP) coupled system (Decharme et al., 2019 ). ISBA calculates the time evolution of the energy and water budgets at the land surface, while CTRIP simulates river discharges up to the ocean from the total runoff computed by ISBA. CTRIP runs at a half degree resolution with a 30‐min time step and is coupled to ISBA via the OASIS‐MCT coupler (e.g., section 2.5 ). ISBA explicitly solves the one‐dimensional Fourier and Darcy laws throughout the soil using 14 layers down to a 12‐m depth and accounting for the hydraulic and thermal properties of soil organic carbon. The snow is represented by a 12 layer detailed internal process snow model including a simple ice sheet runoff to avoid unrealistic snow accumulation over continental glaciers. A two‐way coupling between ISBA and CTRIP is set up to account for, first, a dynamic river flooding scheme in which floodplains interact with the soil and the atmosphere through free‐water infiltration and evaporation and, second, a two‐dimensional diffusive groundwater scheme to represent unconfined aquifers and upward capillarity fluxes into the superficial soil.

SURFEX is a numerical platform which simulates surface fluxes at the Earth's surface (Masson et al., 2013 ). It is coupled inline to ARPEGE‐Climat and shares the same grid and time step. Three surface types are considered: land (including urban areas treated as rock surface), lakes, and ocean.

Compared to CNRM‐CM5.1, the other atmospheric parameterizations correspond to a major update. The convection scheme follows the work of Piriou et al. ( 2007 ) and Guérémy ( 2011 ) and provides a continuous and prognostic treatment of dry, shallow, and deep convection. Convective microphysical processes are treated outside the convective scheme itself, in a consistent way with those occurring in the convection environment (see below). Entrainment and detrainment of energy, moisture, and microphysical species within the convective updraft depend on the prognostic updraft vertical velocity and follow a buoyancy sorting mechanism (Bretherton et al., 2004 ). The scheme closure is based on a dilute Convective Available Potential Energy (CAPE) relaxation. The stratiform microphysics scheme was designed following the work of Lopez ( 2002 ). It takes into account autoconversion, sedimentation, icing‐melting, precipitation evaporation, and collection processes. Bouteloup et al. ( 2011 ) developed a probabilistic approach for the sedimentation, which allows longer time steps than those associated with the original lagrangian approach. The turbulence scheme follows the approach of Cuxart et al. ( 2000 ), which represents the turbulent kinetic energy with a 1.5‐order scheme prognostic equation. The non‐local mixing length is based on Bougeault and Lacarrère ( 1989 ). A specific treatment of the entrainment at the top of the boundary layer is taken into account based on Grenier and Bretherton ( 2001 ). The turbulence scheme also diagnoses the subgrid‐scale variance of the saturation deficit, which depends on both the total water and liquid‐water potential temperature turbulent fluctuations, and which is then used to compute the stratiform cloudiness and the cloud water content in a consistent way (Ricard & Royer, 1993 ; Sommeria & Deardorff, 1977 ). This diagnosed condensation rate serves as an input to the microphysics scheme. Finally, a non‐orographic gravity wave drag parameterization was introduced in CNRM‐CM6‐1, following Lott et al. ( 2012 ), in addition to the orographic gravity wave drag counterpart which remains similar to that of CNRM‐CM5.1 (Catry et al., 2008 ; Déqué et al., 1994 ).

The same radiation parameterization as the one used in CNRM‐CM5.1 is activated, that is, the Rapid Radiation Transfer Model (Mlawer et al., 1997 ) for the longwave part of the spectrum and a six‐band shortwave scheme originally developed by Fouquart and Bonnel ( 1980 ). Some refinements were undertaken, including a new aerosol climatology (with revised optical properties and first indirect effect) based on an interactive aerosol simulation (paper in preparation). The cloud optical properties were also updated.

The dynamical core is based on a two time level semi‐Lagrangian numerical integration scheme tagged as cycle 37 of the ARPEGE/IFS system. A 15‐min time step is used except for the radiative transfer module, which is called every hour. In addition to the six prognostic variables already considered in the previous version (vorticity, divergence, temperature, specific humidity, surface pressure, and ozone concentration), the model includes 10 new prognostic variables (cloud and precipitating solid and liquid water contents for both stratiform and convective parts, turbulent kinetic energy [TKE], and convective vertical velocity). As the semi‐Lagrangian dynamical core is not fully conservative (Lucarini & Ragone, 2011 ), a global dry air mass and water conservation procedure are activated at each time step after application of a local correction to the water condensates, following Bermejo and Conde ( 2002 ). Note that there is no correction applied on energy yet, such a development is planned for the next model version.

The atmospheric component of CNRM‐CM6‐1 is based on the version 6.3 of the global atmospheric model ARPEGE‐Climat. Its former version 5.2 has been described in Voldoire et al. ( 2013 ). A summary of the main characteristics of the new version is presented in the following with the emphasis on major updates. ARPEGE‐Climat is a spectral model derived from the ARPEGE/IFS (Integrated Forecast System) numerical weather prediction model developed jointly by Météo‐France and the European Center for Medium‐range Weather Forecast. A linear triangular truncation Tl127 is adopted together with a corresponding reduced Gaussian grid (Hortal & Simmons, 1991 ). The model horizontal resolution is about 1.4° at the equator. The CNRM‐CM6‐1 atmospheric component has 91 vertical levels (31 for CNRM‐CM5.1), following a progressive hybrid σ pressure discretization; the highest level is set at 0.01 hPa, while the boundary layer is described with about 15 levels below 1500 m.

Global annual mean surface air temperature (K) evolution for the CNRM‐CM6‐1 piControl experiment over 1,000 years in black and for the 10 members historical experiments in colors; dots indicate the first‐year annual mean for each historical simulations, respectively.

Figure 3 illustrates the two‐step behavior of the CNRM‐CM6‐1 spin‐up: it is particularly strong and rapid during the first 100 years, which are then followed by a weak warming trend in global surface air temperature (GSAT). The warming trend is still present in the piControl experiment but limited to 0.02 K/century. The most noticeable feature is a large multidecadal variability superimposed on this trend. This feature will be discussed further in section 5 , dedicated to variability. Averaged over the 1,000 years of the piControl simulation, the net heat flux over the ocean is 0.15 W/m 2 ; thus, CNRM‐CM6‐1 has not reached equilibrium. This imbalance is consistent with a global sea water potential temperature trend of 0.03 K/century. Compared to CNRM‐CM5.1, the trend in the global mean sea water salinity is significantly reduced. Similarly, the top of the atmosphere imbalance is reduced from 3.6 W/m 2 in CNRM‐CM5.1 to 1.5 W/m 2 in CNRM‐CM6‐1. This confirms that CNRM‐CM6‐1 better conserves energy and mass than its former version.

Tuning is aimed at limiting the spin‐up drift before reaching equilibrium, so that the final state remains close to the observed one. The model spin‐up has been run under fixed preindustrial conditions corresponding to piControl CMIP6 forcings. The ocean component starts from the World Ocean Atlas 2013 (WOA, Locarnini et al., 2013 ). The initial sea ice concentration is set to 100% and its surface temperature at −10 °C where sea surface temperature is at freezing point. The sea ice thickness is initialized to 2 m in the Arctic and 1 m in the Antarctic. The initial thickness of the snow layer upon sea ice is set to zero. The land surface reservoirs are initialized from an equilibrium state obtained by running SURFEX in offline mode forced by atmospheric fields taken from a preindustrial AMIP simulation with ARPEGE‐Climat.

In the development of CNRM‐CM5.1, water conservation was an important focus, which resulted in a weak drift of the model ocean water mass (Table 1 ). However, salt conservation remained an issue, with a relatively high salinity drift of −0.01 psu/century in the ocean (Table 1; Gupta et al., 2013 ). In CNRM‐CM6‐1, a specific effort has been made to ensure the conservation of water and salt in each component and in the coupling procedure. Concerning energy, all components with the noticeable exception of the atmosphere are conservative. In the atmospheric model ARPEGE‐Climat, the effects of physical parameterizations are expressed in a flux‐conservative formalism and therefore the physics is intrinsically conservative for mass and energy. But, as mentioned earlier, the semi‐Lagrangian dynamical core does not ensure energy conservation and no “energy fixer” has been implemented yet in ARPEGE‐Climat. However, under fixed conditions, the model should reach an equilibrium state after a spin‐up phase.

In the following, the preindustrial simulation (known as piControl) of CNRM‐CM6‐1 uses time‐invariant 1850 forcings for GHG, AOD, and ozone coefficients, whereas a present‐day land cover state (taken as an average around the 2000s) is imposed. The simulation of the industrial era climate (as known as historical) uses the observed evolution of these forcings from 1850 to 2014 except for the land cover which is fixed to its present‐day distribution.

CNRM‐CM6‐1 uses a time‐invariant distribution of land‐cover, which is based on ECOCLIMAP‐II (Faroux et al., 2013 ), and thus does not follow the recommendation of Eyring et al. ( 2016 ) to account for land use and land cover change forcings. In absence of a fully resolved carbon cycle, the relevance of land cover change can be questioned. Besides, the impact of land cover changes on surface radiative balance remains small with respect to the other forcings (Myhre et al., 2013 , IPCC AR5).

CNRM‐CM6.1 total ozone column (TOC, Dobson unit, DU) now agrees well with the two reference data sets over the 1980–2000 period, in terms of both absolute values and trend (Figure 2 , top). This was a clear deficiency of CNRM‐CM5.1, as reported in Eyring et al. ( 2013 ). Before 1960, the CNRM‐CM6‐1 TOC is overestimated by about 5 DU compared to the CMIP6 reference ( https://esgf‐node.llnl.gov/search/input4mips/ ). This bias is mainly driven by tropical regions. It is of opposite sign in the northern latitudes (a too weak Brewer‐Dobson circulation could then compensate the positive bias of the tropics).

The ozone treatment has been marginally updated since CNRM‐CM5.1 where the evolution of the ozone mixing ratio is computed following the Cariolle and Teyssèdre ( 2007 ) linear ozone parameterization. In CNRM‐CM6‐1, the parameterization remains linear, but recommendations of Monge‐Sanz et al. ( 2011 ) have been followed so as to include implicitly heterogeneous chemistry in the scheme coefficients. As for AODs, monthly mean coefficients of the linear ozone parameterization have been computed from preliminary AMIP‐type runs of ARPEGE‐Climat 6.3 in which the interactive chemistry scheme REPROBUS (see Morgenstern et al., 2017 , for details) is activated from 1950 to 2014 (a three‐member ensemble has been considered). In CNRM‐CM6‐1, coefficients before 1950 are those of 1950. This approach also ensures a high consistency in the ozone representation between CNRM‐CM6‐1 and its ESM version.

Area‐weighted annual global means of total ozone column (DU, top), stratospheric (middle) and tropospheric (bottom) aerosol optical depth (AOD) at 550 nm used in the historical simulations of CNRM‐CM5.1 (blue) and CNRM‐CM6‐1 (red). For the total ozone column, the NIWA3.3 data set (1980–2010, update of Bodeker et al.,) is in gray line, the CMIP6 ozone data set in black line, and the 90% confidence level has been added in colors based on the spread coming from the ensemble of 10 historical members in CNRM‐CM5.1 and CNRM‐CM6‐1.

As presented in Figure 2 (bottom), the preindustrial tropospheric AOD in CNRM‐CM6‐1 is close to the one in CNRM‐CM5.1, but the increase during the historical period is weaker, albeit in the range of the CMIP5 ensemble AOD (see Figure 9.29 in Flato et al., 2013 ). This is explained by two factors. First, the contribution of natural aerosols, especially sea salt particles, is stronger in CNRM‐CM6‐1 than in CNRM‐CM5.1, while it is the opposite for sulfate aerosols. Second, over the twentieth century, the sulfate AOD increase is weaker in CNRM‐CM6‐1. Note that observations cannot constrain it yet over the whole historical period (Carslaw et al., 2017 ).

For the stratospheric aerosols, mainly emitted by volcanic eruptions, CNRM‐CM6‐1 uses the official CMIP6 data set adapted for the 550‐nm wavelength (Thomason et al., 2018 ). The main differences with the Ammann et al. ( 2007 ) data set used in CNRM‐CM5.1 occur during the 19th century, where volcanic eruption peaks are different. During the twentieth century, the main eruptions are quite similar in both data sets, even if the intensity and the length of the eruptions may change. The CMIP6 data set also includes a significant stratospheric aerosol concentration background that was not present in Ammann et al. ( 2007 ).

As in CNRM‐CM5.1, CNRM‐CM6‐1 uses monthly averages of aerosol optical depth (AOD) for five tropospheric aerosol categories (sulfate, black carbon, organic matter, sea salt, and dust) and one stratospheric aerosol category to represent their interactions with radiation and clouds. These AODs are vertically distributed using a fixed profile for each category. AODs of tropospheric aerosols were precomputed based on an AMIP‐type simulation in which the Tropospheric Aerosols for ClimaTe In CNRM‐CM (TACTIC) interactive aerosol scheme (as detailed in Michou et al., 2015 ) is activated. This simulation is driven by the CMIP6 surface emissions of anthropogenic and biomass burning emissions of short‐lived climate forcers and CMIP6 AMIP SSTs and sea ice cover and covers the period from 1850 to 2014. To remove the effect of model internal variability, simulated monthly AODs of the five aerosol tropospheric species were low‐pass filtered by applying a 11‐year moving average. In CNRM‐CM5.1, the Szopa et al. ( 2013 ) aerosol AODs were used. For CNRM‐CM6‐1, the method employed ensures more consistency first between the aerosol forcing and the atmospheric model climate and second between CNRM‐CM6‐1 and its Earth System Model (ESM) counterpart (Watson et al., 2018 ).

Greenhouse gases (GHG) concentrations are the yearly global averages of Meinshausen et al. ( 2017 ) for CO 2 , N 2 O, CH 4 , CFC12 and a CFC11eq species that includes the effect of all the other GHG of the original data set (39 species); this grouping corresponds to option 2 of Meinshausen et al. ( 2017 ). Yearly global averages of the total solar irradiance forcing were similarly used (Matthes et al., 2017 ), as for CNRM‐CM5.1. Additionally, a geothermal heating flux is applied at the ocean floor (Emile‐Geay & Madec, 2009 ): it varies spatially, is constant in time, and has a global mean value of 86.4 mW/m 2 .

In its CMIP6 modeling setup, CNRM‐CM6‐1 uses, unless otherwise detailed below, the suite of external forcing and boundary conditions recommended by Eyring et al. ( 2016 ), that is, the version 6.2 data collection release.

In the third and last step, coupled ocean‐atmosphere simulations under fixed forcings and starting from observed oceanic conditions were carried out. The final objective was to prevent global‐mean sea surface temperature (SST) of the coupled system from drifting by targeting a net surface energy balance on the ocean surface close to 0 W/m 2 after a spin‐up phase. AMIP‐type simulations were used to reach this zero imbalance at a reduced computational cost. The net energy flux entering the ocean corresponds to the sum of surface turbulent, radiative and enthalpy fluxes associated with water mass exchanges (mainly from runoff). These enthalpy fluxes are not diagnosed in AMIP type simulations. They were estimated from coupled simulations to 0.8 W/m 2 . By accounting for these additional fluxes, the target surface energy balance in a fixed SST experiment must be roughly −0.8 W/m 2 for the Earth to be considered in radiative equilibrium. We focus here on the net surface energy balance rather than the top‐of‐the‐atmosphere energy budget knowing that the atmospheric model is not fully conservative and considering the numerical energy sink of the atmosphere. Furthermore, the Earth is not presently in radiative equilibrium with a global mean imbalance of the order of 0.5 W/m 2 (as estimated by Hansen et al., 2011 , for the recent 2000–2004 period). Therefore, the final AMIP target was fixed to about −0.3 W/m 2 . To perform this final tuning of the radiative imbalance in AMIP and then coupled configurations, only one parameter (shortwave inhomogeneity scaling factor) was modified, in order to prevent overtuning (Dommenget & Rezny, 2018 ). This parameter mainly modifies the SW term and allows to reach global energy balance without modifying much the main properties obtained in the previous steps.

For the land/atmosphere component, simulations using the Atmospheric Model Intercomparison Project (AMIP) protocol were performed with interactive land surface and prescribed observed sea surface temperatures and sea ice cover; the present‐day period, corresponding to years 1979–2014, was used. The zero‐order metrics for these atmosphere‐only simulations were the globally averaged top‐of‐the‐atmosphere and surface energy balances. A second class of targets addressed the mean climate and covered climatological spatial maps of the main climate state variables (e.g., temperature, precipitation, surface pressure, and cloud radiative effects). A focus was also put on the precipitation daily distribution. An empirical iterative method was used to tune the parameters relative to cloud and radiative processes (convective entrainment, autoconversion, and inhomogeneity scaling factor for SW and LW cloud optical thickness).

For the oceanic component, OMIP simulations were performed and the calibration was carried out in strong collaboration with the shaconemo group. For CNRM‐CM6‐1, a specific aspect is the use of the GELATO sea ice model; additional tuning was thus made on sea ice melting albedo and drag under sea ice.

The second phase of tuning was based on stand‐alone simulations of each component, to address the interactions between parameterizations and between parameterizations and the other component elements (e.g., dynamics).

The calibration strategy applied for CNRM‐CM6‐1 followed the three main steps defined by Hourdin et al. ( 2017 ). The first step consists in choosing the model formulation for each of its components, mainly their targeted horizontal and vertical resolutions and set of parameterizations. In the case of the new parameterizations introduced in each component, component working groups (ocean, atmosphere, and land surface) performed calibration using idealized or in situ 1‐D cases (see corresponding papers).

4 Modern Mean Climate Evaluation

In this section, we assess the performance of CNRM‐CM6‐1 to replicate modern observations and compare its skill to that of CNRM‐CM5.1. The four major components of the climate system are assessed: atmosphere, land, ocean, and sea ice. In the following analysis, unless otherwise stated, we have used the 10‐member ensemble historical simulations available for both model versions and focused on the 30 years period 1981–2010 to compare with observational or reanalysis data over the same period. In some cases, the period has been restricted to match with the period of availability of reference data sets. Details on the observationally derived data sets have been grouped in Appendix A.

4.1 Atmosphere As the atmosphere component parameterizations have been fully revised in CNRM‐CM6‐1, a thorough evaluation of the atmospheric representation will be done in a forthcoming paper based on the analysis of the AMIP DECK experiment (http://www.umr‐cnrm.fr/cmip6/references). Here, we summarize mainly the mean surface biases that could be affected by the coupling with the other components. Overall, the surface temperature — air temperature (TAS) over land and SST over the ocean — exhibits rather similar patterns of biases in CNRM‐CM5.1 and CNRM‐CM6‐1 (Figure 4). Figure 4 Open in figure viewer PowerPoint Biases (K) in surface air temperature over land and of sea surface temperature over oceans for CNRM‐CM5.1 (first column) and CNRM‐CM6‐1 (second column) in December–March (DJFM; top row) and June–September (JJAS; bottom row). Biases are estimated as the difference between ensemble means of model historical simulations averaged over 1981–2010 and the average of several data sets over the same period. Details on observation products used can be found in Appendix A . Zonal mean errors over ocean (third column) and land (last column) show the inter‐member range (±1,64 times the inter‐member variance) for simulations in red (resp. blue) color shading for CNRM‐CM6‐1 (resp. CNRM‐CM5.1) and the range of data for observational estimates in gray shading. Over the ocean, the most striking change between the two models occurs in the Southern Ocean where the warm bias of CNRM‐CM5.1 has been significantly reduced in CNRM‐CM6‐1 in both seasons. This improvement is probably due to the increase in the shortwave cloud radiative effect over this region (Figure 6). This warm bias is relatively common in CMIP models (Wang et al., 2014) and has been shown to have an impact on tropical mean climate (Cabré et al., 2017; Mechoso et al., 2016); we may thus expect that this regional improvement may provide benefits on other regions in the model. The warm biases in the eastern tropical ocean basins associated with a deficit in stratiform low clouds is weakly reduced except over the tropical north Pacific during boreal summer. Note that the bias in this region significantly varies among the 10‐member ensemble, which possibly indicates a role of the ocean circulation. This type of SST warm bias is systematic to many climate models (Hourdin et al., 2015; Richter, 2015) and several origins have been pointed out in the literature (Zuidema et al., 2016, and references herein). The lack of stratiform low clouds over these regions is probably instrumental in CNRM‐CM6‐1(Brient et al., 2019). Compared to CNRM‐CM5.1, the CNRM‐CM6‐1 low‐cloud cover is improved (Figure 6) but the associated radiative effect has been tuned down to limit the so‐called “too few‐too bright” compensation error (Nam et al., 2012). Besides, Brient et al. (2019) show that the low‐cloud cover still remains underestimated, due to deficiencies in the cloud parameterization itself and to slower feedbacks involving the large‐scale circulation. Over continents, CNRM‐CM5.1 and CNRM‐CM6‐1 mostly show cold biases during the boreal winter (December–March). These biases are stronger in CNRM‐CM6‐1 over the Himalayas, with a wider pattern of cold bias. Over Greenland, the cold bias is slightly reduced in CNRM‐CM6‐1 (note that the values exceed the axis scale on the zonal mean bias figure, for CNRM‐CM6‐1, the bias is ~8 °C, whereas it reaches 12 °C in CNRM‐CM5.1). CNRM‐CM6‐1 is also warmer over North‐Africa, Europe, Scandinavia, and Russia, which mostly constitutes an improvement compared to CNRM‐CM5.1, and to a certain extent, compared to most CMIP5 models (Cattiaux et al., 2013; Woods et al., 2017). However, this warming leads to the appearance of a warm bias over Eastern Siberia in CNRM‐CM6‐1, which is not present in CNRM‐CM5.1. During the boreal summer (June–September), CNRM‐CM5.1 has warm biases over the northern hemisphere continents, as most CMIP5 models (Cattiaux et al., 2013; Christensen & Boberg, 2012; Mueller & Seneviratne, 2014). This systematic bias has been significantly reduced in CNRM‐CM6‐1, and it has almost disappeared in Eastern Europe, Ukraine, and North America's great plains. This improvement is attributed to the soil moisture scheme and to the aerosol‐cloud scheme changes. Indeed, changes in ISBA physics from a simple bucket force‐restore scheme, as used in CNRM‐CM5.1, to a multilayer diffusive soil combined with the representation of groundwater upward capillarity fluxes and floodplain infiltration, as used in CNRM‐CM6‐1, lead to a larger simulated soil moisture and to increase the seasonal memory of the system that favors summer evaporation (Decharme et al., 2019). Beyond the soil scheme improvements, there is also an atmospheric source of improvement, as pictured on figure 6, there is a net increase in cloud shortwave radiation effect over northern hemisphere continents. On the contrary, results are deteriorated in Central Africa where CNRM‐CM6‐1 presents a warm bias which is consistent with the precipitation dry bias (Figure 5). Figure 5 Open in figure viewer PowerPoint Same as Figure 4 for precipitation but with zonal mean absolute values instead of errors on zonal mean plots. Reference datasets are detailed in Appendix A . DJFM = December–March; JJAS = June–September. As for the surface temperature, the mean precipitation field exhibits biases of similar amplitude and sometimes spatial patterns between the two CNRM‐CM versions (Figure 5). The main differences lie in the western tropical Pacific where CNRM‐CM6‐1 has now a dry bias, and in Central and West Africa where, for instance, CNRM‐CM6‐1 partly missed the summertime West African monsoon. On the opposite, the Indian monsoon rainfall pattern is significantly improved and the June–September precipitation excess in the Southern Ocean is largely reduced. The daily distribution of precipitation is also much improved, most probably in relationship with the major update of the convection scheme (not shown). In CNRM‐CM5.1, the occurrence of light precipitation events was largely overestimated to the detriment of the frequency of dry days. The frequency of moderate to heavy rainfalls was also overestimated, while extreme precipitation were not frequent enough. CNRM‐CM6.1 is much more realistic throughout the spectrum of daily precipitation, both over land and ocean. With regard to the global‐mean radiative budget, the strongest differences between CNRM‐CM5.1 and CNRM‐CM6.1 concern the clear‐sky SW surface radiation and upward SW radiation at the top of the atmosphere (Table 2). The first one is probably due the update of aerosols optical properties in CNRM‐CM6‐1. The second one might also be due to these aerosols properties changes but is more likely related to the cloud cover increase which also increases the planetary albedo. In both cases, these differences make the model closer to the Clouds and the Earth's Radiant Energy System (CERES) reference data set. In addition, the cooling SW cloud radiative effect (CRE) is increased in CNRM‐CM6‐1 in the midlatitudes, where it results in a reduced bias (Figure 6), whereas in the (oceanic and continental) convection and trade cumulus tropical regions, it remains overestimated. Over the eastern part of tropical ocean basins, SW CRE biases are reduced, consistently with reduced SST biases, except over the North Pacific where the new atmospheric model has a persistent warm bias. For LW radiation, CNRM‐CM6‐1 global mean fluxes are slightly lower than in CNRM‐CM5.1 (Table 2). The LW CRE is weakly changed from a zonal mean perspective, but regional biases are generally reduced compared to CNRM‐CM5.1. Table 2. Global Mean Radiative Quantities Estimated From CERES Data (2000–2015) and in Both Model Versions on the Period 1981–2010 Averaged Over 10 Members of Historical Simulations W/m2 CERES CNRM‐CM5.1 CNRM‐CM6‐1 Shortwave Downard surface radiation 187 189 194 Clear sky downward surface radiation 244 232 251 Top of atmosphere outgoing radiation 100 96 102 Top of the atmosphere outgoing clear sky radiation 52 54 54 Longwave Downard surface radiation 345 337 334 Clear sky downward surface radiation 317 312 309 Top of atmosphere outgoing radiation 240 240 237 Top of the atmosphere outgoing clear sky radiation 266 262 259 Figure 6 Open in figure viewer PowerPoint Annual averaged cloud radiative effect bias (W/m2) in shortwave (SW; top) and longwave (LW; bottom) in CNRM‐CM5.1 (left) and CNRM‐CM6‐1 (middle) historical simulations (1981–2010) compared to Clouds and the Earth's Radiant Energy System (CERES; 2000–2014) and (right) zonal mean prole of cloud radiative effect for CERES (black) and model simulations CNRM‐CM6‐1 (red) and CNRM‐CM5.1 (blue). The Hadley circulation plays a major role by redistributing momentum and heat from the tropics to higher latitudes. Figure 7 shows the annual climatology of meridional mass stream function, computed from the zonal‐mean meridional wind for two reanalyses, ERA‐Interim (Dee et al., 2014) and MERRA (Rienecker et al., 2011), and for CNRM‐CM5.1 and CNRM‐CM6‐1 model simulations. In ERA‐Interim, the annual Hadley cell is stronger than in MERRA, with the largest difference being in the ascending branch, which is consistent with the findings of Nguyen et al. (2013), who show that MERRA reanalysis produces one of the weakest meridional circulation. Figure 7 Open in figure viewer PowerPoint Zonal‐mean annual‐mean meridional overturning stream functions in isobaric coordinates in ERA‐Interim, MERRA, CNRM‐CM5.1 and CNRM‐CM6‐1 during 1981–2010. Positive values represent clockwise circulation and negative values represent anticlockwise circulation. The maximum and minimum of the stream functions are shown as black markers for ERA‐Interim and MERRA and as blue dot for CNRM‐CM5.1 and red dot for CNRM‐CM6‐1. The spatial structure of the mass stream function derived from CNRM‐CM6‐1 historical simulations and the two reanalyses are consistent with a southern cell higher than the northern cell. The extension of the southern cell at 5°–10°N is also well captured, and the descending branches of the two cells are located at about 30°S and 30°N. The location and the amplitude of the mass stream function maximum value, commonly used as an index to measure the overturning strength, are close to thoses in the reanalyses in the southern hemisphere. In the northern one, the Hadley cell is stronger, both in CNRM‐CM5.1 and CNRM‐CM6‐1 models. The maximum is located at higher altitudes in CNRM‐CM6‐1. However, other reanalyses (e.g., National Centers for Environmental Prediction (NCEP) and JRA) also produced a higher circulation center (Stachnik & Schumacher, 2011). During boreal winter and summer (not shown), the Hadley cells simulated by CNRM‐CM6‐1 also exhibit a higher localization and a larger amplitude of the stream function maximum value, but the biases remain within the range of uncertainties among reanalyses. Given the major update made on atmospheric parameterizations, the surface mean climate state variable are reasonably simulated.

4.2 Land Many changes have also been done on the hydrological core of the ISBA‐CTRIP land surface system in CNRM‐CM6‐1. Today, the Gravity Recovery and Climate Experiment (GRACE) gravity products (Swenson, 2012) estimate all changes in continental water masses through the time variations in terrestrial water storage (ΔTWS). In CNRM‐CM5.1, ΔTWS is comparable to the sum of time variations in snowpack, canopy water, total soil moisture (including soil ice), and river water mass, while in CNRM‐CM6‐1 floodplains and groundwater storage are added. Figure 8 compares the ΔTWS mean annual cycle simulated by CNRM‐CM5.1 and CNRM‐CM6‐1 to GRACE estimates over four regions: Siberia, Europe, Amazonia, and South Asia. Improvements are clearly seen over the tropics and, to a lesser extent, over the northern high latitudes. As shown in Decharme et al. (2019), over Amazonia, this result is attributable both to the new diffusive soil scheme in ISBA and the new groundwater and floodplain schemes in CTRIP. These increase the memory of the system in CNRM‐CM6‐1 and shift the simulated maximum ΔTWS toward that estimated by GRACE. Over Siberia, the better annual cycle simulated by CNRM‐CM6‐1 is due to the use of a variable river streamflow velocity in CTRIP (Decharme et al., 2010) rather than to an improvement in the representation of snow processes (in CNRM‐CM5.1, this velocity was constant and equal to 0.5 m/s). Figure 8 Open in figure viewer PowerPoint Comparison between observed and simulated monthly mean seasonal cycles of Terrestrial Water Storage time variations over Siberia, Europe, Amazonia, and South Asia. Observations are in black, CNRM‐CM5.1 in blue, and CNRM‐CM6‐1 in red; shading shows the inter‐member range for each model (±1,64 times the inter‐member variance). Details on observation products and processing can be found in Appendix A . CSR = Center for Space Research; GFZ = GeoForschungsZentrum; JPL = Jet Propulsion Laboratory. The snow is an active component of the climate system, especially over northern high‐latitude regions where a drastic decrease in springtime snow cover over recent decades (Derksen & Brown, 2012) contributes to the polar warming amplification (Holland & Bitz, 2003). In addition, the snowpack controls the temperature of the soil and then the possible carbon emissions due to the permafrost thawing (e.g., Schuur et al., 2015). A good simulation of both the northern Arctic snow and surface albedo is thus required in climate models to address scientific questions in these regions. Many efforts have thus been made these last years to improve theirs representations in CNRM‐CM6‐1. The mean annual cycles of snow and surface albedo simulated by CNRM‐CM6‐1 are clearly more realistic than in CNRM‐CM5.1 (Figure 9). These results can be directly related to a better representation of snow processes and snow albedo in the new snow scheme (Decharme et al., 2016) used in CNRM‐CM6‐1, while the appropriate albedo value during summer is more due to the use of the Moderate Resolution Imaging Spectroradiometer albedo for the land surface. Figure 9 Open in figure viewer PowerPoint Comparison between observed and simulated monthly mean seasonal cycles of (a) snow cover extents and (b) surface albedo. Observations are in black, CNRM‐CM5.1 in blue and CNRM‐CM6‐1 in red; shading shows the inter‐member range for each model (±1,64 times the inter‐member variance). Details on observation products and processing can be found in Appendix A . NSIDC = National Snow and Ice Data Center; CERES = Clouds and the Earth's Radiant Energy System. The reader is referred to Decharme et al. (2019) for a more comprehensive evaluation of the land component.

4.3 Ocean Overall, the ocean temperature bias is smaller in CNRM‐CM6‐1 than in CNRM‐CM5.1 in all basins (Figure 10). In the Atlantic, a cold bias persists in CNRM‐CM6‐1 in the northern midlatitudes but the near surface bias is much smaller than in CNRM‐CM5.1, probably in part because of the higher ocean vertical resolution. As already noted on SST, the CNRM‐CM5.1 strong warm bias in the Southern Ocean is largely reduced in CNRM‐CM6‐1 and this improvement is also clear at depth. A warm bias persists in CNRM‐CM6‐1 between 300‐ and 1,000‐m depth from 30°S to 30°N, but it is much smaller than in CNRM‐CM5.1, indicating a better representation of the Antarctic intermediate water masses. Figure 10 Open in figure viewer PowerPoint Zonal mean potential temperature bias (shading) in °C in CNRM‐CM5.1 (left column) CNRM‐CM6‐1 (right column) in the Atlantic (top row), Pacific (middle row), and Indian basins (bottom row). The observed climatology from the World Ocean Atlas (WOA) is shown in contours. A logarithmic scale has been applied to the vertical coordinate (depth) to allow a better visualization. CNRM‐CM6‐1 shows saltier North Atlantic waters than observed at high latitudes, (Figure 11) consistent with the strong deep convection simulated in the Gin Seas (Figure 12). In the Pacific and Indian basins, the temperature bias is quite similar between the two models, albeit smaller at high latitudes in CNRM‐CM6‐1. The salinity bias in the Pacific is more hemispheric in CNRM‐CM6‐1 with water masses that are too salty in the Northern Hemisphere and too fresh in the Southern Hemisphere. The fact that this bias is present from the surface to about 300 m suggests that it could be related to an atmosphere bias in net surface water flux. Figure 11 Open in figure viewer PowerPoint Same as Figure 10 for salinity in psu. Figure 12 Open in figure viewer PowerPoint 2017 ρ = 0.03 kg/m3 with respect to density at 10 m, and in Holte et al. ( 2017 2004 Annual maximum mixed layer depth (a) reconstructed from the 2004–2018 ARGO observations (Holte et al.,), averaged over the 1981–2010 period of the ensemble mean historical for (b) CNRM‐CM5.1 and (c) CNRM‐CM6‐1. The mixed layer maximum is computed from monthly means. In models, it is defined with the density threshold of ∆= 0.03 kg/mwith respect to density at 10 m, and in Holte et al. () climatology it follows the definition from de Boyer Montégut et al. (). Figure 12 displays the annual maximum mixed layer depth as modeled by CNRM‐CM5.1 and CNRM‐CM6‐1 and observed from ARGO measurements (Holte et al., 2017). The main observed deep convection areas are located in the Labrador Sea and the Greenland—Iceland—Nordic Seas. They are reasonably well reproduced in both model versions, although with a predominant deep convection in the latter region for CNRM‐CM6‐1. For the Southern Hemisphere, there is no deep convection observed during the ARGO period, but it has been documented in the past over the Weddell Sea (Gordon, 1982) and Ross Sea (Jacobs et al., 1970). CNRM‐CM5.1 displays weak deep convection in the Weddell Sea, whereas CNRM‐CM6‐1 has a semipermanent polynia in the Indian sector of the Southern Ocean. This explains relatively larger root‐mean‐square errors (RMSEs) in the newer version (151 m compared to 87 m), mean biases being comparable and modest (−9 m compared to +5 m). Opposite to CNRM‐CM5.1, CNRM‐CM6‐1 does not simulate deep convection in the Mediterranean Sea whereas it does in the Japan Sea. Finally, both models simulate intermediate convection at subpolar fronts, with an equatorward bias though. The mean barotropic circulation is very similar in both model versions (Figure 13). It is slightly more intense in the newer version, which causes lower biases (−2 Sv compared to −5 Sv) and RMSEs (26 Sv compared to 31 Sv) with respect to Colin de Verdière and Ollitrault (2016) climatology. The main circulation features are represented by both versions: subtropical and subpolar gyres and the Antarctic Circumpolar Current. However, they are largely underestimated, mostly because of the unresolved mesoscale dynamics along western boundaries. The Antarctic Circumpolar Current reaches respectively 87 Sv in CNRM‐CM5.1 and 108 Sv in CNRM‐CM6‐1, far weaker than the latest estimations of 170−173 Sv (Colin de Verdière & Ollitrault, 2016; Donohue et al., 2016). Finally, consistently with mixed layer biases, the separation between subtropical and subpolar gyres associated with the subpolar front is too equatorward in both models, especially in the Northern Hemisphere. Figure 13 Open in figure viewer PowerPoint 2016 (a) Barotropic stream function reconstructed from the World Ocean Atlas 2009 climatology and mean ARGO float drifts (Colin de Verdière & Ollitrault,) and (b) 1981–2010 mean ocean barotropic stream function from CNRM‐CM6‐1 (shades) and CNRM‐CM5.1 (contours) ensemble mean historical simulations (in Sv). In models, the reference null value is taken over the African‐Eurasian continents. The mean Atlantic Meridional Overturning Circulation (AMOC) in CNRM‐CM6‐1 is much stronger than in CNRM‐CM5.1 with a maximum of about 18 Sv located around 35°N and 1,000‐m depth, while in CNRM‐CM5.1 the maximum was 14 Sv at 26°N and about 800 m (Figure 14). The North Atlantic Deep Water (NADW) cell is deeper in CNRM‐CM6‐1, which contributes to the reduced temperature and salinity biases (Figures 10 and 11). The maximum value of the ensemble mean MOC profile at 26°N (Figure 14c) is very close to the observed value (~17 versus 18 Sv, McCarthy et al., 2017) although it occurs shallower than observed (800 m instead of 1,100 m). Overall, the ensemble mean profile in CNRM‐CM6‐1 is very close to observations in the upper 2,500 m indicating a good representation of the upper NADW cell. The lower NADW cell is too shallow as shown by the return flow that occurs around 3,500 m compared to 4,500 m in observations. The spread in CNRM‐CM6‐1 is larger than in CNRM‐CM5.1. This large spread is due to a marked multidecadal and centennial variability in the Atlantic in CNRM‐CM6‐1 as shown by the fluctuations of the AMOC at 26°N (Figure 14d). The AMOC peak‐to‐peak variations are much larger than in CNRM‐CM5.1 and the period of the oscillations is around 200 years whereas it was closer to 100 years in CNRM‐CM5.1 (Ruprich‐Robert & Cassou, 2015). The mechanisms behind this low‐frequency Atlantic variability in CNRM‐CM6‐1 are still under investigation. Preliminary analysis suggests a possible link with freshwater exchanges in the Arctic. Figure 14 Open in figure viewer PowerPoint Mean Atlantic Meridional Overturning Circulation (AMOC; in Sv) in the CNRM‐CM5.1 (a) and CNRM‐CM6‐1 (b) model averaged over the 10 historical members over the period 1981–2010. Vertical profile (c) of the mean (1981–2010) MOC at 26°N (d) in CNRM‐CM5.1 (blue), CNRM‐CM6‐1 (red) and in RAPID observations (black). The plain line corresponds to the 10 historical mean and the shading to the inter‐member spread. The period used for RAPID is April 2004 to February 2017. All the data have been interpolated on the CNRM‐CM6‐1 vertical grid. AMOC index (d) defined by the annual value of the stream function maximum at 26°N in CNRM‐CM6‐1 (red) and in CNRM‐CM5.1 (blue) preindustrial simulations. The AMOC time series of the 10 historical simulations of CNRM‐CM6‐1 are shown in pink. A 5‐year running mean have been applied to all time series. While the atmosphere is the dominant contributor to the total meridional transport (Trenberth & Caron, 2001), biases in the ocean transport have to be compensated for by the atmosphere and can therefore lead to discrepancies in the whole coupled system. It is thus important to assess the fidelity of the climate simulations in representing the meridional transport by the ocean. The global ocean heat transport shown in Figure 15a indicates that more heat is transported northward in CNRM‐CM6‐1 than in CNRM‐CM5.1 in the Northern Hemisphere. The intensity of the peak around 15°N is closer to observational estimates. In the southern midlatitudes, both models show an equatorward transport, which is in agreement with the observation estimate by Large and Yeager (2009) but in disagreement with most other estimates that suggest a poleward transport. The Southern Hemisphere ocean transport was pointed out to be a problematic aspect of many CMIP5 simulations and reanalysis. This highlights the marked difficulties in representing the large‐scale energy processes in the Southern Ocean as discussed by Trenberth and Fasullo (2010). Figure 15 Open in figure viewer PowerPoint (a) Global and (b) Atlantic ocean heat transport simulated in CNRM‐CM6‐1 (red) and CNRM‐CM5.1 (blue) and in different observational estimates (black and gray, see Appendix A for references). Ten historical simulations are considered for the models with the thick plain line corresponding to the ensemble mean and the shading indicating the spread (estimated by ±1,64 times the inter‐member variance). For the Atlantic, only two members are considered for CNRM‐CM5.1 because the necessary data were not available for other members. The mean Atlantic Ocean heat transport shows a marked improvement in CNRM‐CM6‐1 compared to CNRM‐CM5.1, at all latitudes. The maximum transport is located around 25°–30°N, in agreement with observations. It peaks at a value of about 0.9 PW, which remains below the observational estimates of RAPID (Johns et al., 2011), Ganachaud and Wunsch (2003), and Trenberth and Caron (2001) but is about 20% larger than the transport simulated by CNRM‐CM5.1 around that latitude. This increased meridional transport can be partly explained by the larger overturning heat transport in CNRM‐CM6‐1, which results for the stronger AMOC at the latitude of the maximum transport (Figure 14).