Implementation of the CMIP6 Forcing Data in the IPSL‐CM6A‐LR Model

The implementation of boundary conditions is a key aspect of climate simulations. We describe here how the Climate Model Intercomparison Project Phase 6 (CMIP6) forcing data sets have been processed and implemented in Version 6 of the Institut Pierre‐Simon Laplace (IPSL) climate model (IPSL‐CM6A‐LR) as used for CMIP6. Details peculiar to some of the Model Intercomparison Projects are also described. IPSL‐CM6A‐LR is run without interactive chemistry; thus, tropospheric and stratospheric aerosols as well as ozone have to be prescribed. We improved the aerosol interpolation procedure and highlight a new methodology to adjust the ozone vertical profile in a way that is consistent with the model dynamical state at the time step level. The corresponding instantaneous and effective radiative forcings have been estimated and are being presented where possible.


Model and Input Data Description
The Institut Pierre-Simon Laplace (IPSL) Climate Modelling Centre (see https://cmc.ipsl.fr) has set up a new version of its climate model for Phase 6 of the Coupled Model Intercomparison Project (Eyring et al., 2016). The model components and the coupled model itself (known as IPSL-CM6A-LR) are described in companion papers (Boucher et al., 2020;Hourdin et al., 2020). Here we describe how the various boundary conditions that have been prepared for the Climate Model Intercomparison Project Phase 6 (CMIP6) have been processed and implemented in IPSL-CM6A-LR. While CMIP6 is more prescriptive than CMIP5 in the way boundary conditions should be implemented, there are still choices to be made during model implementation with possible implications for the simulated climate forcings and model performance. This is

Long-Lived Greenhouse Gases
The atmospheric (volume) mixing ratios of long-lived greenhouse gases (carbon dioxide CO 2 , methane CH 4 , nitrous oxide N 2 O, and halocarbons) are prescribed from Meinshausen et al. (2017) as per the CMIP6 protocol. We use Version 1.2.0 of the data. For simplicity we prescribe global annual mean mixing ratios without any temporal interpolation (i.e., only the global mean mixing ratio is used and the value is updated on 1 January of each year in experiments where greenhouse gases vary). The differences between CMIP5 and CMIP6 for CO 2 , CH 4 , and N 2 O are negligible for the historical period but are noticeable for the scenario levels that are common, for example, the RCP85 (Pathway 8.5) and SSP585 (Pathway 5-8.5) scenarios (see Figures 1a-1c) (Meinshausen et al., 2020). It should be noted that methane presents latitudinal and vertical variations that may be interesting to represent in the model in the future. Meinshausen et al. (2017) propose to prescribe mixing ratios of halocarbons using one of the two following groupings: either CFC11-eq and CFC12 (chlorofluorocarbons) or HFC134-eq (hydrofluorocarbons) and CFC12-eq. We choose the second option in the IPSL model as it allows to identify the contribution of the gases controlled under the Montreal protocol, represented in terms of "equivalent CFC-12," from  (ppmv, ppbv, or pptv according to the gas considered) as prescribed in IPSL-CM5A-LR (CMIP5) and IPSL-CM6A-LR (CMIP6) for CO 2 , CH 4 , N 2 O, CFC11-eq, and CFC12-eq in the historical experiment, in CMIP5's RCP scenario experiments, and in CMIP6's four corresponding SSP scenario experiments. Additional points show the values for the PMIP4 (paleoclimate) experiments. The CMIP6 data are from Meinshausen et al. (2017). the others represented in terms of "equivalent HFC-134a." Because CFC11 and CFC12 are the only two halocarbons that can currently be fed to our radiative code, we first convert HFC134a-eq into a CFC11-eq mixing ratio using the Intergovernmental Panel on Climate Change Fifth Assessment Report values for their respective radiative efficiencies (Myhre et al., 2013, their Table 8.2; CFC11, RF = 0.062 W m −2 for 238 pptv and HFC134-eq, RF = 0.0100 W m −2 for 62.7 pptv), which gives a conversion factor of (0.0100∕62.7)∕(0.062∕238.0) to go from HFC134a-eq to CFC11-eq mixing ratios (fed to the model CFC11). The CMIP6 CFC12-eq mixing ratio is fed to the model CFC12 without any change. It should be noted that the CFC11-eq and CFC12-eq mixing ratios are nonzero for preindustrial conditions because they include some natural halocarbons. The time series of the mixing ratios for CFC11-eq and CFC12-eq are shown in Figures 1d and 1e, which shows some difference between CMIP5 and CMIP6 because of the different lumping of gases.
In the 1pctco2 experiment and its variants, the CO 2 mixing ratio is incremented by 1% from its preindustrial value of 284.32 ppm every 1 January starting from Year 2.
In all experiments aiming at representing past equilibrium climate states, that is, for the mid-Holocene (midHolocene), Last Interglacial (lig127k), mid-Pliocene Warm Pariod (midPliocene-eoi400), and Last Glacial Maximum (lgm), long-lived greenhouse gases are set to a constant value for the duration of the simulation, as recommended by the PMIP4-CMIP6 protocols (Kageyama et al., 2018(Kageyama et al., , 2017Haywood et al., 2016;Otto-Bliesner et al., 2017). Those values are particularly low for the lgm experiment, representing a glacial climate, and CO 2 is significantly higher than the piControl value for the midPliocene-eoi400 experiment, which represents a climate warmer than present (Figure 1).
We have computed the instantaneous radiative forcings (IRFs) associated with a doubling and quadrupling of the CO 2 concentration from double radiation calls. They amount to 2.50 and 5.75 W m −2 , respectively, at the top of atmosphere, and 4.06 and 8.50 W m −2 , respectively at the tropopause level, using the preindustrial climate as a reference (see Table 2). It should be noted that the evolution of IRF is not strictly logarithmic with respect to the increase in the CO 2 concentration, with the IRF at the tropopause for 4xCO 2 being 5% larger than twice the IRF for 2×CO 2 . This supralogarithmicity of the CO 2 radiative forcing was already flagged by Zhong and Haigh (2013) at least for very large concentrations. A similar result was also found by Etminan et al. (2016) for the stratospherically adjusted radiative forcing.
We have also computed the effective radiative forcings (ERF) for both 2×CO 2 and 4×CO 2 experiments. A variety of methods have been used. First, fixed-SST (sea surface temperature) runs have been performed to diagnose the ERF (see Table 2, fourth column), following Forster et al. (2016). Second, we use the method by Gregory (2004), with linear regressions of the top-of-atmosphere net radiative flux against surface temperature change (Figure 2). In this case the linear regressions were carried out over the first 20 years of the abrupt experiments, with a subtraction of the piControl variables (for temperature and radiative flux) Figure 2. Computation of ERF for abrupt simulations using the Gregory (2004) method for 2 × CO 2 (a) and 4 × CO 2 (b, c). For (a) and (b), the linear regressions are carried out on the first 20 points (i.e., first 20 years) of the diagrams. On (c), 12 members are considered, each of them with a different starting month, over 5 years; the regression is carried out on the corresponding 60 monthly points. The corresponding forcings are listed in Table 2. both on a year-to-year basis (to remove the effect of low-frequency variability) or using 20-year climatologies (see Table 2, fifth and sixth columns, and Figures 2a and 2b-the climatology subtraction method is not illustrated in Figure 2). For the 4×CO 2 experiment, we also used 12 different members with different starting months and carried out the linear regression over 5 years ( Figure 2c). The values for all IRF and ERF discussed above are synthetized in Table 2, and ERF are displayed in Figure 3. Table 3 presents IRF and ERF values from fixed-SST sensitivity simulations using GHG (greenhouse gases) concentrations from Year 2014, considered as "present-day" in the RFMIP (Radiative Forcing Model Intercomparison Project)-ERF protocol (Pincus et al., 2016), and three paleoclimate periods considered in PMIP4-CMIP6: Pliocene, Last Interglacial, and Mid-Holocene. Note that the Pliocene CO 2 value is almost identical to the one for Year 2014.

Original and Hybrid Data Sets
With the two exceptions discussed below, we have used three-dimensional ozone fields (Version 1.0) from University of Reading (Checa-Garcia et al., 2018;Hegglin et al., 2016) for the piControl, historical, and scenario experiments. For the Detection and Attribution MIP experiments, we have used the ozone fields provided by Chemistry-Climate Model Initiative (v1.0 for hist-stratO3 and ssp245-nat and v1.1 for hist-nat, hist-sol, and hist-volc) through the input4MIPs repository. These ozone fields have been corrected to reflect the impact of the forcings that are being considered in the individual Detection and Attribution MIP experiments. For instance, the ozone fields for the hist-nat experiments include the effects of natural forcings on the ozone profiles but not the effects from anthropogenic forcings. Some ozone fields for Tier 2 ScenarioMIP experiments were not available at the time of running the experiments. We have thus performed a few scenario experiments with ozone fields that were scaled using monthly outputs from the CNRM-CM6 model (Michou et al., 2020;Voldoire et al., 2019). This followed three steps. First, since the CNRM-CM6 output data and the official ozone fields have different vertical levels, an interpolation of the CNRM-CM6 sets of data was performed onto the vertical levels of the official UReading data. Then, a grid-point-to-grid-point scaling was operated. Considering a sspN scenario, let Hybrid N denote the Tier 2 ozone field that was to be computed, CNRM N the ozone fields output by the CNRM-CM6 model, and UReading N the official ozone field (only available for Tier 1 SSP at the time of computation). For each of the Tier 2 SSP to be processed, we determined the two neighboring Tier 1 SSP, lo(N) and hi(N), in terms of radiative forcing reached in 2100-with the exception of SSP1-1.9 having only one neighbor-such as   Etminan et al. (2016) and in the IPCC Third Assessment Report (Ramaswamy et al., 2001 We performed two scalings for each Tier 2 SSP required, using the following equations: where the multiplications and divisions are meant point-to-point over the grids. Every Hybrid {hi(N),lo(N)} pair was then averaged into a single hybrid file Hybrid N , with the exception of ssp 119 being simply scaled from ssp 126 . As for the lower troposphere region, CNRM fields below 560 hPa are relaxed to the 560-hPa level (Michou et al., 2020). Our method is thus imperfect. However, we could assess that the ratios between ozone fields of a CNRM Tier 1 ssp and a CNRM Tier 2 ssp are close to unity in the troposphere. Therefore, in the tropospheric region, the final averaging between the hybrid files becomes close to averaging the two hi and lo UReading ozone fields, for which the troposphere is duly represented. Since the CNRM-CM6 output data do not reach as high in altitude as the official forcing files, our last step was a "hatting" process, designed to complement the upper part of the ozone field. For this purpose, we simply averaged the high tops of the two neighboring Tier 1 official fields and appended it to the synthetized hybrid Tier 2 ozone fields. A second set of experiments was performed with the official CMIP6 fields when these became available. These experiments are labeled as r?i1p1f1 when using the official CMIP6 ozone fields, and r?i1p1f2 when using the hybrid CNRM ozone fields. As an example, Figure 4 shows the time evolution of the zonal mean ozone column for the 1960-2014 period of the historical and the scenario SSP4-6.0 as per the original input4MIPs data.

Interpolation Procedure
The ozone volume mixing ratios forcing fields are interpolated in space and time for insertion into the atmospheric model of IPSL-CM6A-LR in three steps: • a horizontal interpolation onto the IPSL-CM6A-LR atmospheric grid (performed offline); • a linear interpolation in time from monthly means to daily values (performed online, to avoid heavy offline manipulation of forcing files); and • a vertical linear interpolation on log-pressure levels (performed online as this depends on the model state).
The third step concerns the interpolation in the tropopause region. It requires a particular attention because the tropopauses of the original CMIP6 ozone field and the IPSL atmospheric field are not identical. The  ) and, as an example, the SSP4-6.0 scenario (2015-2100). Both the official CMIP6 data and the hybrid CNRM data are displayed. The data have been annually averaged. The internnual variability in the two data sets, which have been ratioed results a higher variability in the processed hybrid data. model used to prepare the ozone forcing fields and the IPSL climate model follow different dynamical trajectories and have different dynamical biases, so their tropopause levels differ not only day to day but also in their climatological mean. This means that a simple vertical interpolation of the ozone on the IPSL model levels could lead to a net transfer of ozone from the lower stratosphere to the upper troposphere or vice versa. For a given ozone perturbation, the radiative forcing is largest around the tropopause region (see Forster & Shine, 1997, their Figure 9), so a spurious net ozone transfer above or under the model tropopause may affect not only the day-to-day variations in the ozone radiative forcing but also, and more significantly, its climatological mean. Indeed, Figure 5 shows a significant mean height difference: the chemical tropopause of the climatological input fields is about 400 m higher than the dynamical tropopause of the IPSL-CM6A-LR model on average, with significant latitudinal variations: about +1 km below 50 • S, −300 m above 40 • N and +300 m in between these latitudes. To remove these effects, each time an ozone field is used, its vertical profiles are deformed in order to match the IPSL-CM6A-LR model tropopauses. The deformation is confined to a thin region (of a few kilometers) around the tropopause, so most of the vertical profile shape remains untouched.

Deformation Procedure in the Tropopause Region
Two definitions of the tropopause height at a particular geographical location are used depending on whether the forcing field or the model field is considered: • the dynamical tropopause is identified as the lowest level between the level with an Ertel potential vorticity of |PV| = 2 PVU and the level with a potential temperature of = 380 K; and • the chemical tropopause is identified as the level with an ozone concentration of [O 3 ] = 100 ppbv.
The original CMIP6 ozone forcing data were not provided with the fields needed to compute |PV| and , so the dynamical definition is not an option and the chemical definition is used. Conversely, the IPSL model has no interactive ozone, so the tropopause location is determined using the dynamical definition. The purpose of the algorithm is to deform the CMIP6 ozone forcing field to force both tropopauses to coincide. Our method to perform this is new, and similar to that of Hardiman et al. (2019), which was developed separately. The different steps of the algorithm are explained hereafter: • Locate the chemical tropopause height, f , in the ozone forcing file. The tropopause is identified in the forcing file as the chemical tropopause height, defined here as the level of [O 3 ] = 100 ppbv. We use a simple search algorithm because forcing fields are relatively smooth and their vertical resolution is limited. After preliminary testing, a pressure level of 90 hPa was found to be an adequate starting level for searching down the chemical tropopause. Starting from a much higher level (lower pressure) runs the risk of identifying as chemical tropopause the very low spring ozone minimum that forms in the lower stratosphere under Antarctic ozone hole conditions. • Locate the IPSL model dynamical tropopause, m . The rather noisy PV vertical profiles (built from instantaneous fields) are smoothed with a Gaussian-shaped function extending on five levels. Two iso-PV points are detected, respectively from top to ground ( = 0 to 1) and from ground to top ( = 0.75 to 0). If the PV between them is greater than 2 PVU, the first point is retained, the second one otherwise. The height of the model dynamical tropopause is extremely variable. For example, the very low tropopause heights typical of stratospheric ozone filaments entering the upper troposphere, would lead the tropopause matching algorithm to inject ozone at low altitudes. This is not desirable, because such fine structures are not accurately enough captured (a larger time-space resolution of the model and refreshing rate of the forcing fields one day would be required) to reliably trigger net ozone transfers. Therefore, tropopause pressures outside the 80-400 hPa interval are clipped, limiting the amplitude of the stretching of the ozone field in extreme situations. Stretching the forcing ozone profile and interpolating it on actual model levels would be inaccurate, considering the low vertical resolution of the input data and the two stages of interpolation needed. We chose instead to strain the model levels with the reciprocal transformation and interpolate the original ozone vertical profile directly on these adjusted levels. The following vertical deformation transformation, limited to a bounded region around m , is used:̃= where • = ln ( ∕ m )∕ ln m is the exponent required for a power law to move the m point to f • (see Figure 6) is a pseudo-hat function restricting the effect of the power law to a bounded region, defined in two parts: ln( m ∕ 0 ) . In every column and at each time, the deformed layer is bounded with prescribed values b (bottom) and t (top); it must contain the region bounded by m and f and have sufficient thickness to ensure that (i) the deformed region is discretized on several layers of the atmospheric model, (ii) the deformed pressure profile variations stay monotonic, and (iii) the deformed ozone profile stays smooth. In practice, 40% of the minimum thickness (e.g., the distance between forcing and model tropopauses) on each side is a sufficient margin. A typical result of our interpolation scheme is depicted in Figure 7. Comparison of two identical simulations except for the tropopause deformation procedure indicates a radiative budget difference of about −0.25 W m −2 at the top of the atmosphere. Our tropopause adjustment procedure largely prevents the occurrence of this artificial radiative forcing linked to the difference in tropopause height between the model and the initial ozone forcing field. The negative sign is compatible with the fact that the tropopause inferred from the forcing files is on average slightly higher (about 400 m) than its atmospheric model counterpart (see Bekki et al., 2013, and Figure 5).
It should be noted that the tropopause adjustment was used in all experiments including those of the DECK (Diagnosis, Evaluation, and Characterization of Klima). The resulting tropospheric and stratospheric ozone columns in the IPSL-CM6A-LR model are shown in Figure 8 for the historical and some of the scenario experiments.

Tropospheric Aerosols
Tropospheric aerosol loads were computed using LMDZOR-INCA historical runs, in order to then be prescribed as climatological data sets in the IPSL-CM6A-LR models. The model was run in atmosphere-only mode; hence, SST and sea ice cover were required as boundary conditions. SST and sea ice from input4MIPs were used for the historical period 1850-2014 (these data have an only version, though augmented through time). For the scenario period (2015-2100), SST and sea ice cover were derived from CMIP5 outputs of the IPSL-CM5A-MR model driven by the RCP scenarios. They were bias corrected for the period 2005-2017 based on the input4MIPS historical data set, with the sea ice cover bounded between 0 and 1. The anomalies for SST and sea ice cover were then interpolated from the four RCP scenarios to the CMIP6 scenarios based on the 2100 radiative forcing target values. LMDZOR-INCA was used in its aerosol-only configuration (i.e.. without the full gas phase chemistry switched on). Oxidants (i.e., OH, O 3 , and H 2 O 2 ) were prescribed from the CMIP5 outputs of the IPSL-CM5A-LR model coupled with the INCA chemistry-aerosol model (Szopa et al., 2013). The same mapping as for the SST is used from the RCP to the SSP scenarios.
For the historical period, five different time segments were defined, each of them tallying with 33 years of the complete 1850-2014 period. For each of the five segments, a prior 3-year spin-up period was run, starting with an atmosphere free of SO 2 and aerosols, and with forcings corresponding to the first year of the segment considered. This made the five runs actually cover 36 years each, from which the 3-year spin-up periods were eventually withdrawn and not considered in the final synthesis of climatologies. This procedure was adopted to speed up the production of the aerosol forcing fields in order to start the historical simulations as early as possible in the CMIP6 calendar. For future scenarios, the whole 2015-2100 period was computed in a single LMDZOR-INCA run; the scenario runs directly restarted from the end of the fifth historical segment (i.e., 31 December 2014).
Only surface emissions as prescribed by the CMIP6 protocol were considered. Three-dimensional emissions from aviation (AIR) were neglected. All aerosol emissions were prescribed using the CMIP6 input data (Hoesly et al., 2018;van Marle et al., 2017), with Version 1.2 for biomass burning emissions, and v2017-05-18 for anthropogenic emissions. In the absence of further information in the CMIP6 data set, biomass burning emissions were emitted at the surface. For the three spin-up years preceding 1850, emission rates were duplicated from the 1850 data set.  . An example of the interpolation correction as derived from the algorithm exposed in Sheng and Zwiers (1998). The CMIP6 original data (light blue line), when input as midmonth values and then interpolated over time, lead to erroneous monthly averages (red lines); this is corrected with the use of the Sheng and Zwiers algorithm (black lines).
For nitrate aerosols, emissions of nitrogen containing species (i.e., nitrogen dioxide NO 2 and ammonia NH 3 ) are considered for both anthropogenic and openburning emissions and for both the historical period and future scenarios. The nitrate aerosol mechanism described by Hauglustaine et al. (2014) has been used and modified for this simplified version of the INCA model without interactive gas phase chemistry. Several new gas phase species (NH 3 , NO 2 , and nitric acid HNO 3 ) have been introduced in this INCA aerosol only version and are now calculated interactively with the main oxidants prescribed as described before. According to the simplified gas phase mechanism considered in this version, NO 2 is oxidized by OH to form HNO 3 and nitric acid is photolyzed back to NO 2 . Wet and dry deposition of the gas phase species and the formation of ammonium sulfate and ammonium nitrate particles are calculated as described by Hauglustaine et al. (2014).
The aerosol loads output from the LMDZOR-INCA runs underwent a smoothing process, using a sliding, weighted average over a 3-year window centered on the year of interest. The chosen weights were [0.25; 0.5; 0.25]. For boundary years of the historical period 1850 and 2014, the previous (resp. following) year data were considered as equal as those of the year of interest. For future scenarios, this concerned only 2100, since 2015 could be readily averaged with 2014 values from the historical run. This procedure was carried out on a month-to-month basis, so as to keep a monthly evolution in the final forcing climatology files.
Since CMIP6 emissions rates are provided at a monthly frequency, and to avoid step-like evolution of emission rates throughout the year, the input aerosol emissions were linearly interpolated in time. Though easy to implement, this operation obviously entails that monthly averages of the interpolated values shall differ from the original monthly values. To maintain consistency, we relied on the algorithm exposed in Sheng and Zwiers (1998): The CMIP6 original monthly values are first considered as midmonth values and are then corrected in a manner such that when the linear time interpolation is effectively carried out, monthly averages calculated over one calendar month should yield the actual CMIP6 original value. The linear time interpolation combined with the Sheng and Zwiers algorithm ensure both (i) continuity in the emission rates of aerosols in time and (ii) consistency with the monthly averaged CMIP6 values.
The Sheng and Zwiers algorithm in its original form (Sheng & Zwiers, 1998) can correct some small numbers into negative numbers when there is a large change from one month to the next, which of course is non-physical for emission rates. This can happen in particular when dealing with biomass burning emissions which exhibit a large seasonal cycle and small emissions during the nonburning season. To avoid this, we improved the Sheng and Zwiers algorithm with a "freezing" process inspired by Taylor et al. (2000). Considering a single point of the grid, any month that would yield a negative emission value when corrected through the Sheng and Zwiers algorithm is declared as "frozen." In this case, the whole of the month duration will be associated to a constant emission value, corresponding to the prescribed original input value, which avoids the correction to a negative number. The original Sheng and Zwiers algorithm is then slightly adapted: For a month preceding a frozen month, the linear interpolation now occurs from midmonth to the beginning of the frozen month, and likewise (mutatis mutandis) for a month following a frozen month. This  way we keep the continuity of the emission rates in time and still have monthly averages equal to the prescribed monthly value, avoiding the issues arising from possible correction to negative numbers. By simple linearity, global total emissions or averages are kept consistent with the input CMIP6 values, as illustrated in Figure 9. Figure 10 shows the burdens of the main aerosol types for the two decades 1850-1859 and 2000-2009. The concentrations are prescribed in IPSL-CM6A-LR with a daily time step at the same horizontal and vertical resolutions using a linear time interpolation scheme between the monthly means. The aerosols are then hydrated according to the local conditions and the aerosol properties (extinction coefficient, single scattering albedo, and asymetry parameter). Figure 11 shows the time evolution of the global mean aerosol optical depth in the historical and scenario experiments.
As mentioned above, tropospheric aerosols are prescribed as monthly three-dimensional fields in IPSL-CM6A-LR. The fields are interpolated daily between the current month and either the previous or the next month depending on the day in the month. Aerosols are treated as an external mixture for aerosol-radiation interactions. Their optical properties have been recomputed for IPSL-CM6A-LR for the six wavebands of the shortwave radiative transfer model using wavelength-dependent refractive indices, prescribed size distribution and Mie theory. Only the cloud albedo effect is considered for aerosol-cloud interactions. The cloud droplet number concentration (CDNC, in cm −3 ) depends on the total sum of accumulation-mode aerosols (in g m −3 ) based on a revised formulation of Boucher and Lohmann (1995): where b 0 = 1.3 and b 1 = 0.2 have been slightly adjusted as compared to the previous version in order to obtain a better fit to the radiation budget. It should be noted that CDNC is bound between 10 and 1,000 cm −3 . Furthermore, the effective cloud droplet radius computed from the cloud liquid water content and the CDNC has a minimum value of 5 μm. Aerosol radiative forcing estimates are provided in Table 4 using double radiation calls to compute the IRF and fixed-SST experiments to compute ERF. The total aerosol radiative forcing is relatively small (in magnitude) at −0.62 and −0.73 W m −2 for the ERF and IRF estimates, respectively.

Stratospheric (Volcanic) Aerosols
Stratospheric aerosols were provided by CMIP6 as a two-dimensional (latitude-height) time-varying climatology of the aerosol extinction coefficient, single scattering albedo, and asymmetry parameter (Thomason et al., 2018) averaged over the wavebands of each model radiative code. The radiative code of the IPSL-CM6A-LR comprises six wavebands in the shortwave and 16 wavebands in the longwave parts of the spectrum. For CMIP6, we used Version 3 of the data set, keeping in mind that Version 4 has become available after we started our CMIP6 experiments. Version 4 corrects the gap filling process applied after the Pinatubo eruption, which results in extinction coefficient that is lower around 22-to 23-km altitude in the tropics (i.e., the region of largest mass loading) and larger around 15 km, with a decrease in aerosol heating rate. The resolution of the original data set is 0.5 km on the vertical, 5 • in latitude and monthly in time, and has to be mapped onto the resolution of the atmospheric model. In latitude we simply map the coarse resolution of the original data set on the more refined resolution of the atmospheric model by choosing the closest latitude. On the vertical two issues had to be solved to interpolate the original CMIP6 data set onto the IPSL-CM6A-LR atmospheric grid. First, the vertical coordinate is altitude while the atmospheric model is on a coordinate. Second, the original data set extends below the tropopause and the tropopause is not specified (and its height cannot be diagnosed from the aerosol vertical profile). To address these issues, we interpolate offline the original data set onto the vertical grid of the model by using a 10-year average monthly climatology of the model geopotential field. In this way seasonal and latitudinal variations in the relationship between pressure and altitude due to temperature variations be accounted for. We then diagnose the pressure of the tropopause using the algorithm from Reichler et al. (2003) and only consider the aerosol above the tropopause. This is justified by the fact that most of the aerosol under the tropopause is not of volcanic origin and it is therefore already accounted for by our tropospheric aerosol scheme (see section 4.1). For some years and latitudes we may nevertheless miss volcanic aerosols sedimenting into the troposphere (e.g., as it is the case in the polar troposphere a year after the Pinatubo eruption). Finally, we prescribe the aerosol as monthly mean without any temporal interpolation as we do not want to smooth the large peaks in aerosol optical depth introduced by large volcanic eruptions. Figure 12 shows the stratospheric aerosol optical depth after interpolation and removal of the tropospheric contribution in the IPSL-CM6A-LR model.
For the future (scenario) period, we use yearly data for the Years 2015 to 2023 (included), after which an average input is used with no yearly dependence. This follows the recommendations of the CMIP6 protocol. It should be noted that a minor decrease in stratospheric aerosol optical depth during scenarios 2025 onward because there is a slight upward shift in the tropopause level with global warming.

Oceanic External Supply of Nutrients
NEMO-PISCES is the ocean biogeochemical model used for IPSL-CM6A-LR (Aumont et al., 2015). As the biogeochemical cycles are not fully coupled between the oceanic, atmospheric and land components, fluxes of nutrients are required for PISCES. These are supplied to the ocean from different sources and more particularly from atmospheric deposition and rivers.

Atmospheric Deposition
The model considers the atmospheric supply of iron Fe, silicon Si, phosphorus P, and nitrogen N elements. The atmospheric cycles of these elements have all been perturbed by anthropogenic activities but to different extents (e.g., Wang, Balkanski, Boucher, Bopp, et al., 2015;. Here we only consider the natural sources of Fe, Si, and P through dust deposition, while we also consider N from anthropogenic sources (reactive nitrogen emitted by combustion).
For the atmospheric supply of soluble iron, phosphorus, and silicate, we use modeled deposition rate from outputs of the LMDZOR-INCA model. These are from the same simulations, and we apply the same smoothing process as for tropospheric aerosols (see section 4.1). We computed the total dust flux by summing the monthly mean dry deposition, wet deposition and sedimentation fluxes. Figure 13 shows the time evolution of the global mean of the total dust over the ocean for the historical and scenario experiments.
These natural dust forcings are applied in PISCES under the assumption that mineral soluble fraction of dust in P, Si, and Fe is fixed regardless of the environmental conditions. Once it has left the surface layer, particulate inorganic iron from dust is still assumed to experience dissolution. The dissolution rate is computed assuming that mineral particles sink at a constant speed and that about 0.01% of the particulate iron dissolves in a day. The iron deposition at surface computed by the model is illustrated in Figure 14.
Atmospheric deposition of Si is restricted to the first level of the ocean model. Atmospheric deposition of P is computed from dust deposition assuming that the total phosphorus content of dust is 750 ppm and that the solubility in surface sea water is 10%. As for Si, deposition is restricted to the first level of the ocean model.
The nitrogen deposition forcing is coming from the wet and dry deposition of inorganic nitrate (drynoy and wetnoy variables) and ammonium (drynhx and wetnhx variables). These data are coming from the National  Center for Atmospheric Research-Chemistry-Climate Model Initiative nitrogen deposition rate prepared for input4MIPs. The nitrogen deposition fields read in PISCES consist of the sum of NO 3 and NH 4 fields which are both deposited at the ocean surface and are supposed to be 100% soluble independently of atmospheric conditions. As the nitrogen deposition forcing is unchanged from the input4MIPs data, it is not shown on a figure.

River Discharge
River supply of biogeochemical elements to the ocean is taken from the GLOBAL-NEWS2 data set (Mayorga et al., 2010) for dissolved inorganic nitrogen, dissolved inorganic phosphorus, dissolved organic nitrogen, dissolved organic phosphorus, and Si. For dissolved inorganic carbon and alkalinity, we use results from the Global Erosion Model of Ludwig et al. (1996), neglecting the particulate organic carbon delivery as most of it is lost in the estuaries and in the coastal zone (Smith & Hollibaugh, 1993). All fields are interpolated onto the ocean model grid and colocalized with the river runoff prescribed in the physical model. Iron is also delivered to the ocean by rivers. The amount of supplied iron is computed from the river supply of inorganic carbon, assuming a constant Fe/dissolved inorganic carbon ratio. This ratio is determined so that the total Fe supply equals 1.45 Tg Fe year −1 as estimated by Chester (1990). In the absence of specific historical reconstructions and future scenarios for the evolution of riverine inputs of carbon and nutrient, these inputs are fixed at constant present-day values for the entire duration of all simulations performed with IPSL-CM6A-LR.

Solar Forcing
We apply Version 3.2 of the solar forcing data set prepared by Matthes et al. (2017). We have used the MATLAB routine provided by the authors to compute the fraction of the solar constant in each of the six shortwave wavebands of the radiative transfer model of IPSL-CM6A-LR. The model follows a Gregorian calendar and data are prescribed at the daily resolution. We do not make use of cosmic rays, solar protons, and medium-energy electrons. As for the future period, we use the standard (higher) forcing option, though the lower option could be used in experiments to come. Figure 15 shows the variations of total solar irradiance through time, over both historical and scenario periods.
When going back in time, changes in Earth's orbit have to be accounted for. The five time periods considered in CMIP6 encompass various orbital configurations (Kageyama et al., 2018). For these periods, changes  in the eccentricity, obliquity, and precession of the equinoxes have little impact on the global annual mean incoming solar radiation at the top of the atmosphere (insolation). They mostly modulated the annual mean solar forcing between the equator and the poles and the seasonal distribution ( Figure 16). For each simulation we prescribed the orbital parameters as requested in the corresponding protocol (i.e., Haywood et al., 2016, for midPliocene-eoi400; Otto-Bliesner et al., 2017, for lig127k and midHolocene; Kageyama et al., 2017, for lgm;and ;Jungclaus et al., 2016. The changes in the Earth's orbit are associated with different lengths of the seasons, or time spent between equinoxes and solstices (Berger & Loutre, 1991, 1994Joussaume & Braconnot, 1997). It is thus important for model comparisons that all simulations have the same calendar date assigned on a reference point on the ellipse. For all time periods, the vernal equinox is prescribed to occur on 21 March each year. The use of the Gregorian calendar is such that the bissextile years are also accounted for in past climates, and correspond to the artificial dates used to run the simulations. The Pliocene is run with the same insolation as today. For the other period the zonal mean changes plotted as a function of months highlight the large seasonal changes induced in insolation for the last interglacial (lig127k) and the mid-Holocene (6,000 years BP). These two climates are characterized with increased (decreased) insolation seasonality in the Northern Hemisphere (Southern Hemisphere), with larger differences with piControl found for lig127k (Otto-Bliesner et al., 2017). These changes reach 60 W m −2 during norther hemisphere summer in high latitude for lig127k and 30 W m −2 for midHolocene.
The changes imposed for the last glacial maximum are very small in comparison and small compared to the ice sheet forcing. As expected, these figures realized from the output or the different simulations are similar to the one of the corresponding protocol papers (see Otto-Bliesner et al., 2017, Figure 3; Kageyama et al., 2017, Figure 2), when the modern calendar is used to compute monthly values. Note however that the orbital parameters in piControl are set to 1990 and not 1850 as recommended in CMIP6 (Eyring et al., 2016). This introduces a difference smaller than 0.5 W m −2 in the insolation forcing (Otto-Bliesner et al., 2017). The changes in the Earth's orbit are also imposed for the past1000 simulations. The initial state in Year 850 is still affected by seasonal changes induced by the Earth's orbit. Even though it is reduced compared to the mid-Holocene, orbital parameters still impose changes in seasonality that combines with the other trace gazes and land use forcing at the global and regional scales. , whose partitioning varies at the grid cell level. The model uses yearly land use maps prescribing the spatial distribution of the different vegetation types and the had biomass (sections 6.1. and 6.2). Additionally, a single soil texture map for all years is used, describing the global distribution of different soil types used to compute the thermal and hydrological soil properties. Each model grid cell is assigned to a single dominant soil type from three possible types, coarse (sandy loam), medium (loam), and fine (clay loam), which are derived from the Zobler soil texture map (Zobler, 1986), where the initial five soil classes are reduced to three classes. Note also that soil albedo is prescribed using information from the MODerate-resolution Imaging Spectro-radiometer observations. We used the bihemispherical reflectances (white-sky albedo) for the Years 2001-2010 at 0.5 • resolution for the visible and near infrared spectral domains to define spatially varying but fixed in time soil albedo.

Land Cover
We describe here the input data and algorithms used to create the land cover maps specific for our CMIP6 simulations using the historical/future reconstruction of land use states provided as reference data sets for CMIP6 within the land use harmonization database (LUH2v2h, Hurtt et al., 2011). More details are provided on the devoted web page (https://orchidas.lsce.ipsl.fr/dev/lccci/), which shows further tabular, graphical, and statistical data. The overall approach relies on the combination of the LUH2v2 data with present-day land cover distribution derived from satellite observations for the past decades. The main task consists in allocating the land use types from LUH2v2 in the different PFTs for the historical period and the future scenarios. The C3 crop fraction in IPSL-CM6A-LR is directly set as the sum of the C3 annual crop, C3 perennial crop and C3 N fixing crop from LUH2v2, and similarly C4 crop fraction as the sum of annual C4 crop and perennial C4 crop. The pasture fraction from LUH2v2 feed the C3 and C4 PFTs, using the partitioning from (Still et al., 2003) at the grid cell level. The total fraction-within each grid cell-occupied by the other land use states from LUH2v2 (i.e. primary forested land, secondary forested land, primary nonforested land, secondary nonforested land and rangeland) is assumed to correspond to the natural vegetation. To characterize the natural vegetation in terms of PFT distribution, one makes use of the land cover (LC) product, v2.0.7, developed within the land cover project of the Climate Change Initiative (CCI) program of the European Space Agency (ESA). Consistent time series of global LC maps are available at 300-m spatial resolution on an annual basis from 1992 to 2015. The data are classified into 38 land cover classes from the Food and Agricultural Organisation of the United Nations Land Cover Classification System (LCCS) (Di Gregorio & Jansen, 2005). To be usable within the climate model, the data have been remapped onto 15 generic PFTs at 0.25 • resolution with the help of the software provided with the ESA-CCI land cover project (user tool, v3.13) and a cross-walking table defining the correspondence between the LCCS classes and generic PFTs (Poulter et al., 2015). When allocating LCCS classes in generic PFTs, additional grid cell information on bioclimatic zones and the partitioning between grasslands with C3 and C4 photosynthetic pathway are also used following the Köppen-Geiger climate classification map and the data set from Still et al. (2003), respectively, priorly enlarged to five major climate groups. Once disaggregated to generic classes, the data are then reaggregated into the model-specific PFTs. Note that in our model configuration, shrub, and tree generic PFTs are grouped into unique tree PFTs.
The natural vegetation in each grid cell is defined as the PFT distribution derived from the ESA-CCI land cover product for the year 2010 to which pasture fraction and crop fraction from LUH2v2 (for the Year 2010) have been subtracted from grass and crop PFTs. This characterization of the natural vegetation in terms of PFT distribution is assumed invariant in time and is used for both the historical period and the different future scenarios. Note finally that cities, lakes, and inland water are not treated separately but assigned to the bare soil fraction.
Resulting historical and future (for a set of different scenarios) time series for the major model PFT groups are shown in Figure 17. For the historical period, the crops cover increased by 10.7 Mkm 2 from 1860 to 2015 (148.7%), while the forest cover decreased by approximately the same amount (10.1 Mkm 2 or 15.3%). Grasses and bare soil PFT covers are relatively stable, although at regional scale the temporal dynamic of the different PFT groups is more variable (see PFT maps for different time periods at the reference web page). For the future scenarios, the selected socioeconomic scenarios lead to substantial differences in the evolution of crop, grass and forest covers, with potentially large impacts on the land surface properties. For the lgm, midHolocene, and lig127k simulations, land cover is not modified (apart from the new ice sheets for the

Wood Harvest
Removal of harvested wood biomass is accounted for in our modeling framework, using the data from LUH2v2 (Hurtt et al., 2011). Harvested wood biomass is imposed on a carbon (C) mass basis. A low estimate of the annual harvested C wood biomass per unit of forest area in each grid cell is computed as the ratio of the total harvested wood C biomass (sum of the transition variables primf_bioh, primn_bioh, secmf_bioh, secyf_bioh, and secnf_bioh of LUH2v2) divided by the related land cell area (sum of the states variables  For the paleoclimate conditions, the altitude is shown with two different scales: for land covered by ice sheets in the lgm case or by the piControl ice sheets in the midPliocene-eoi400 case, the red to blue scale is used; in other cases, the green to brown scale is used. The contour of the continents is shown in white for the piControl case, in purple on the other maps. The ice sheet extent is shown with a magenta contour for the piControl case, a blue contour for the lgm case, and a red contour in the midPliocene-eoi400 case. Bering Strait and Canadian Archipelago passage are closed in both the lgm and midPliocene-eoi400 cases. primf, primn, secdf, and secdn multiplied by the area). On an annual time step, this amount of C is removed from the aboveground biomass of tree PFTs.

SST and Sea Ice Cover for piClim Experiments
The climatic SSTs and sea ice cover used as boundary conditions in the piClim experiments of RFMIP were obtained by averaging 200 years of model output from the CMIP6 DECK piControl experiment performed with IPSL-CM6A-LR. More specifically, the average was performed on simulation years 51 to 250 (i.e., model years 1900 to 2099) of the piControl simulation.

Other Specific Forcings for Paleoclimate Simulations: Land-Sea Mask, Land Ice Extent, and Topography
The forcings used for the paleoclimate simulations defined for CMIP6 via the PMIP4 project include changes in long-lived greenhouse gases (section 2 and Figure 1), insolation (section 5 and Figure 16) and, in the case of the midPliocene-eoi400 simulation, in land cover (section 6 and Figure 18). Here, we further describe the changes in boundary conditions required for the midPliocene-eoi400 and lgm simulations.
For both these simulations, the coastlines have been modified to account for sea level and topographic changes ( Figure 19). In both cases, the Bering Strait is closed, the Canadian Archipelago belongs to the American continent and is linked to Greenland. Hudson Bay is replaced by land in the midPliocene-eoi400 boundary conditions, and by the land ice surface type for the lgm one. For the lgmsimulation, other large changes include the Sunda and Sahul shelves, between South East Asia and Australia, being exposed due to the lowering of the sea level results from the massive ice sheets. The Fennoscandian ice sheet also covers the Barents-Kara Sea, which results in a shrinking of the passage between the Nordic Seas and the Arctic Ocean. The coastlines are first modified for the definition of the ocean model domain-in practice, new land points are simply assigned a 0-m bathymetry. For the LGM, the ocean bathymetry is adjusted according to the reconstruction (see below) in crucial areas such as west and east of Iceland, where the depth of the sills could have an impact of deep water formation and transfer to the North Atlantic.
For both the midPliocene-eoi400 and lgm simulations, the ice sheet extent has been modified ( Figure 19), with reduced ice sheets for the warmer mid-Pliocene Warm Period climate and new, extensive ice sheets over northern North America and Fennoscandinavia for the LGM climate. Ice sheet changes also result in changes in altitude, over and outside ice sheet locations ( Figure 19). For the midPliocene-eoi400 simulation, both the Greenland and Antarctic ice sheets are less extensive. For the lgm simulation, the North American and European ice sheets represent changes in altitude up to 3 km higher than the present configuration.
Two reconstructions are being tested, from W. R. Peltier (ICE-6G_C) and L. Tarasov (GLAC-1D), following the PMIP4 protocol (Kageyama et al., 2017). In practice, we compute a topography at the same resolution as used to define the piControl altitude-related parameters, that is, the mean altitude, and also the second moment characteristics needed for the orographic drag parameterization (Lott & Miller, 1997), by adding the reconstructed LGM-present anomaly to the present altitude. We then smooth the altitude over the ice sheets to avoid spurious spatial variability arising from the reference (present-day) data set and use this new topography to compute updated parameters for the orographic drag scheme.

Conclusions
We have described the details of how the IPSL Climate Modelling Centre has implemented the boundary conditions of the CMIP6 experimental protocol in the IPSL-CM6A-LR model. Although the CMIP6 protocol is more precise than that of CMIP5, there are choices to be made on how to implement some of the forcing terms. As IPSL-CM6A-LR is run without interactive chemistry, tropospheric, and stratospheric aerosols as well as ozone have to be prescribed. Two major improvements have been implemented relative to the IPSL-CM5A-LR model. First, we have improved the interpolation of aerosol emission with a procedure that conserves monthly emissions and ensure positive values by "freezing" daily values in months that would have seen negative values through the original implementation of the Sheng and Zwiers algorithm. Second, we highlight a new methodology to adjust the ozone vertical profile in a way that is consistent with the model dynamical state of the tropopause at the time step level.
CMIP6 simulations were run successfully following this implementation of the forcing terms, which allowed the rapid publication of our climate data sets on the ESGF. A companion paper will present and assess the results obtained in the historical and scenario runs (Boucher et al., 2020). A comparison between the performance of the IPSL-CM5 and IPSL-CM6A-LR models will be shown. Particular aspects of the forcing terms not covered here (e.g., past millennium) will be covered in a separate publication.
Implementing and testing the CMIP6 forcing terms in the IPSL-CM6A-LR climate model has represented an important piece of work. While it is crucial that various aspects of climate simulations keep improving, including the representation of boundary conditions or forcing terms, it is important that this is done in a way that maintain current format and publication procedure to make such technical tasks as easy as possible.

Author Contribution
T. L. performed most of the preparation of the IPSL-CM6A-LR forcings, in particular the climatologies for tropospheric aerosols and the CO 2 ERF-abrupt calculations, and coordinated the writing of the manuscript. D. C. designed, developed, and evaluated the matching tropopause scheme for O 3 . R. M. H. prepared the solar forcing. A. S. performed the IRF and ERF fixed-SST calculations. N. L. and M. K. tested the climatologies for stratospheric aerosols. C. É. prepared the nutrient forcing. N. V., V. B., and P. P. designed and prepared the land use forcing. O. B. supervised the IPSL-CM6A-LR model development, the preparation of forcing terms, and the writing of this manuscript. All other authors contributed to the IPSL-CM6A-LR model development, the forcing preparation, and the writing the manuscript.