Improved Representation of Clouds in the Atmospheric Component LMDZ6A of the IPSL‐CM6A Earth System Model

The cloud parameterizations of the LMDZ6A climate model (the atmospheric component of the IPSL‐CM6 Earth system model) are entirely described, and the global cloud distribution and cloud radiative effects are evaluated against the CALIPSO‐CloudSat and CERES observations. The cloud parameterizations in recent versions of LMDZ favor an object‐oriented approach for convection, with two distinct parameterizations for shallow and deep convection and a coupling between convection and cloud description through the specification of the subgrid‐scale distribution of water. Compared to the previous version of the model (LMDZ5A), LMDZ6A better represents the low‐level cloud distribution in the tropical belt, and low‐level cloud reflectance and cover are closer to the PARASOL and CALIPSO‐GOCCP observations. Mid‐level clouds, which were mostly missing in LMDZ5A, are now better represented globally. The distribution of cloud liquid and ice in mixed‐phase clouds is also in better agreement with the observations. Among identified deficiencies, low‐level cloud covers are too high in mid‐latitude to high‐latitude regions, and high‐level cloud covers are biased low globally. However, the cloud global distribution is significantly improved, and progress has been made in the tuning of the model, resulting in a radiative balance in close agreement with the CERES observations. Improved tuning also revealed structural biases in LMDZ6A, which are currently being addressed through a series of new physical and radiative parameterizations for the next version of LMDZ.


Introduction
On average, two thirds of the Earth's surface is covered by clouds (King et al., 2013), which are therefore of primary importance in the energy budget of the atmosphere. Similarly, cloud response to global warming is one of the largest sources of uncertainty in climate change simulations Dufresne & Bony, 2008;Vial et al., 2013). From the early stages of climate modeling at the "Laboratoire de Météorologie Dynamique" (LMD) (Sadourny, 1975;Sadourny & Laval, 1984;Laval et al., 1981), efforts were made to develop innovative subgrid-scale parameterizations that correctly represent their effect (Le Treut & Li, 1991;Li, 1999). The current LMD global atmospheric model, called LMDZ for its zooming capability (Hourdin et al., 2006), is the atmospheric component of the Earth system model named after the French climate institute where it is developed: the Institut Pierre-Simon Laplace (IPSL) Climate Model or IPSL-CM. This paper describes the representation of clouds in the latest version of LMDZ, LMDZ6A, which was used for the sixth phase of the Coupled Model Intercomparison project (CMIP6, Eyring et al., 2015). The general descriptions of the IPSL-CM6A model and its atmospheric component, LMDZ6A, can be found in this Special Collection, in two papers by Boucher et al. (2020) and Hourdin et al. (2020), respectively. The two previous versions of LMDZ were LMDZ5A and 5B and are described in Hourdin, Foujols, et al. (2013) and Hourdin, Grandpeix, et al. (2013). Compared to Version 5A, Version 5B was based on a profound rethinking of the parameterization of convection and clouds, on which the new 6A version is built. In the present paper, we will focus on comparing LMDZ5A with LMDZ6A directly, because Version 5B was in many respects the prototype of Version 6A. CMIP5 revealed a number of biases in LMDZ5A cloud properties: • Despite major development efforts, tropical and subtropical low-level cloud fractions were underestimated; • Mid-level clouds were almost inexistent; • Large biases were found in the total cloud radiative effect (CRE) over the Southern Ocean; • Low-level cloud cover was underestimated, and cloud reflectance was overestimated (Konsta et al., 2016); • The altitude of low-level clouds was too low (Konsta et al., 2016).
Our goal in this paper is twofold: to review the entire set of cloud parameterizations developed for LMDZ, and to present the main improvements of the newest version, LMDZ6A. Particular care has been given in the LMDZ parameterizations to the representation of convection, for which a deliberate choice was made to separate deep and shallow convection and which is coupled to cloud description through the specification of subgrid-scale distribution of total water or saturation deficit. These developments have been described through a series of publications but always focusing on one particular aspect. The present paper provides a full description of the parameterizations that control clouds in LMDZ as well as their interactions. In terms of evaluation, particular attention will be paid to the global cloud distribution and its role in maintaining the global radiative balance in the model. The discussion and conclusion will highlight the remaining biases and present the current development efforts to address them.

Parameterization of Clouds in LMDZ6A
The challenge in modeling clouds resides in the various scales of atmospheric processes controlling their macrophysical and microphysical properties. They depend on both km-scale and μm-scale processes evolving on time scales ranging from minutes to seconds. In the last two decades, the LMDZ team worked on a set of innovative parameterizations that describe the subgrid-scale vertical motions and their connections to cloud properties. Clouds in LMDZ depend on (1) turbulent mixing, shallow convection, deep convection, and large-scale horizontal advection, and (2) cloud statistical schemes that use the physical information provided by these processes to compute their opacity and the fraction of the gridbox they cover. To do so, atmospheric properties such as the area covered by thermal plumes in the boundary layer or mass fluxes in deep convective clouds are used to shape the subgrid-scale distributions of water vapor. The general approach is to represent these distributions by probability density functions (PDFs) that can be unimodal or bimodal and whose variance and asymmetry toward high humidity values increases when convective plumes bring near-surface moist air toward the drier free troposphere. Since the temperature of the gridbox is known, it is possible to derive, from these distributions, the populations of air parcels that are supersaturated and to deduce the cloud fraction and water content.
All the processes occurring in a gridbox (turbulent mixing, shallow and deep convection) are called sequentially in LMDZ, as represented in Table 1. In the following sections, the different steps of this diagram will be described, from the main model prognostic variables to the final cloud fraction α c and water content q in c , which are the two information used by the radiative transfer scheme to compute cloud radiative heating rates.

Evaporation
The first procedure of the LMDZ physical package is the evaporation of all condensates, because most parameterizations of convection work with the total water mass mixing ratio q t (see the early work of Betts, 1973). This does not mean that clouds are purely diagnostic. The cloud liquid and ice mixing ratios are "semi-prognostic" variables in the sense that they are advected by the dynamical core, but they are evaporated/sublimated at each time step at the beginning of the physical package. This assumption may hold for liquid droplets whose lifetime is often smaller than the physics time step of ∼15 min but can be a limitation for ice or mixed-phase clouds. This first procedure is represented in Table 1 and affects the three water phases (the water vapor, liquid water, and ice mass mixing ratios, noted q v , q l , and q i ) as well as the potential temperature θ, through evaporative cooling. It returns the total water content q t which is then used and updated by all the cloud parameterizations. The only other procedure affecting the prognostic variables q l and q i is the so-called large-scale condensation scheme, which condenses, before calling the radiative transfer scheme, all the water vapor in excess of saturation coming from the different parameterizations.

Local Turbulent Mixing
The first process that is accounted for is the local turbulent mixing in the boundary layer, which was revisited in LMDZ6A. It now includes a 1.5 order closure K-gradient scheme and a prognostic equation for the TKE (turbulent kinetic energy). The K-gradient scheme is based on the work of Yamada (1983) and was improved for stable boundary layers Vignon et al., 2017). The total water vapor mass mixing ratio q t is vertically mixed assuming a down-gradient Fick's type diffusion whose intensity Table 1 Architecture of the Physical Package, Showing All Cloud-Related Variables Note. The first column gives the names of the different procedures, that are also used as subsection titles in section 2. The second column indicates the main variables used by the procedure on the left, and the prognostic variables that are updated at the end of the procedure on the right, in gray. The other useful variables computed by each procedure are given in the last column. Variables colored in blue are related to cloud properties, and are those used by the radiative transfer scheme to compute the cloud radiative effect. All the notations are given in the text and summarized in Appendix A. depends on the TKE. As is classical in climate models, the turbulence scheme includes the representation of exchanges with the surface, including the evaporation and sensible heat fluxes, which are essential to cloud formation.

Deep Convection
The deep convection scheme of LMDZ computes heating, moistening, and momentum changes using a modified version of the Emanuel mass flux scheme (Emanuel, 1991) to which a parameterization of cold pools was added Grandpeix, Lafore & Cheruy 2010). Version 6A differs significantly from Version 5A which was using the Emanuel scheme without the improved mixing representation (Grandpeix et al., 2004) and the various improvements described in the present section.
Once the turbulent mixing in the boundary layer has been computed, deep convection can be initiated, depending on the ALE (available lifting energy) inherited from the previous time step. The ALE can be provided by frontal lifting at the edge of cold pools or by boundary layer thermals, which are noted ALE wk and ALE th in Table 1, respectively. The ALE finally used by the deep convection scheme is the largest of the two energies. Deep convection is triggered if the ALE exceeds the CIN (Convective INhibition) and if at least one of the cumulus of the domain reaches a given threshold size and evolves into a congestus or cumulonimbus cloud. This latter process is represented by a stochastic triggering scheme (Rochetin et al., 2014) and is also a new feature of LMDZ6A. Another important new feature of Version 6A is the inclusion of the latent heat exchange due to the liquid ↔ ice phase change in the deep convection scheme.
Once deep convection has been triggered and the cold pools have been initiated, the column is split into two separate fields: the cold pool area and its environment, each having their own temperature and humidity. Deep convection then "sees" the environment of cold pools, rather than the mean grid cell, while downdrafts fall inside the cold pool region. This so-called splitting technique is essential to maintain deep convection within the grid cell. The deep convection closure is based on the ALP (available lifting power; see Grandpeix, Lafore & Cheruy 2010), which is inherited from the previous time step and is the sum of the ALP provided by the cold pools and by the thermal plumes of the boundary layer.
The deep convection scheme then computes the in-cloud water mass mixing ratio q in; cv c , which is the ratio of condensed water mass to in-cloud air mass. Note that this quantity is different from the liquid or ice mass mixing ratios within a gridbox q l and q i which correspond to the ratio of condensed water mass to gridbox air mass. It also computes the convective rainfall and snowfall P cv l; i . The precipitation mechanism follows Emanuel and Ivkovi-Rothman (1999): All the condensate in excess of a temperature-dependent conversion threshold is converted into large hydrometeors that will eventually fall. The precipitation efficiency (i.e., the fraction of large hydrometeors in the total condensate) is bounded by a maximum value ep max , which is usually slightly lower than 1 (see Table 3) to always keep some cloud water in the atmosphere (Bony & Emanuel, 2001). All the condensate is carried up in the updrafts till their ends, at which point the large hydrometeors fall as precipitation with a prescribed terminal velocity.
In our scheme, both the undiluted updrafts and the mixed drafts contribute to the in-cloud water content of deep convective clouds. The deep convective cloud fraction α cv c is computed (as explained in section 2.4) from the in-cloud water content of deep convective clouds q in; cv c , which is itself deduced from the different mass fluxes and coverage fraction of undiluted and mixed updrafts. In the case of undiluted updrafts, the coverage fraction α a is given by α a ¼ M a =ðρw a Þ where M a is the mass flux density and w a the vertical velocity. In the case of the mixed drafts, the entrained air at each level feeds cloud formation, and these clouds dissipate with a time constant τ m . Therefore, the time evolution of the cloud water mass in a layer of thickness δz can be written as where M t is the mass flux density of the mixed drafts and q m its condensed water mixing ratio. The coverage fraction of mixed drafts can then be deduced from Equation 1 by assuming a steady state, which gives α m ¼ M t τ m =ðρδzÞ. The in-cloud water content is finally calculated as a linear combination of the cloud water of the undiluted updraft and mixed drafts: where q a is the condensed water mixing ratio of the undiluted updraft. In Equation 1, the saturated draft dissipates with a time constant τ m of the order of 100 s.
This in-cloud water content q in; cv c is computed for use in the radiative transfer scheme and in the deep convective cloud statistical scheme (see the next section), but it is not removed from the vapor phase or used to derive the prognostic variables q l and q i . At the end of the deep convection scheme, the vertical profiles of convective rainfall and snowfall P cv l; i are returned and removed from the vapor phase, and only θ and the total water mass mixing ratio q t are changed accordingly. The deep convection scheme also returns the change in both temperature and water content due to downdrafts dθ cv dw and dq cv t; dw , which are later used by the cold pool scheme (see section 2.5).

Deep Convection PDF
As briefly mentioned at the beginning of this section, the cloud statistical schemes used in LMDZ are tightly connected to the information provided by the shallow and deep convection schemes. Such statistical schemes rely on a PDF describing the subgrid-scale distribution of water vapor or saturation deficit. In the case of deep convection, the total mass mixing ratio of water q within the gridbox is assumed to be a random variable of mean value q t . The latter can be written as The cloud mixing ratio q l or q i and cloud fraction α c can then be computed as ðq − q sat Þ PðqÞ dq; and (4) where q sat is the water vapor saturation mixing ratio at the gridbox mean temperature and pressure, that is, q sat ðT ; pÞ. We neglect in this case the effect of temperature heterogeneities on q sat . The gridbox mean amount of both condensates and in-cloud vapor, q tc , can be written as In this context, the in-cloud water content q in c is given by The deep convection scheme provides the in-cloud water content q in; cv c , as described in section 2.3. Therefore, the three free parameters of a lognormal PDF are then deduced from Equations 3 and 8 by an inverse procedure, assuming that the PDF equals 0 for q ¼ 0 (Bony & Emanuel, 2001, Appendix A). The PDF is then used to compute α cv c , which is later used, together with q in; cv favor new convective zones at their edges on the other hand. Therefore, they play an important role in the lifecycle of convective systems. Their representation is a new feature of LMDZ6A. To account for this process, the deep convection scheme assumes that a fraction of precipitation (15% above cloud base and 100% below) falls outside the cloud and evaporate to form precipitating downdrafts. The cold pool scheme then uses the change in both temperature and water content due to these downdrafts dθ cv dw and dq cv t; dw . As explained earlier, we use a splitting technique so that cold pools can have their own temperature and humidity. The cold pool scheme also derives its own ALE and ALP quantities that will be later used, at the next time step, by the deep convection scheme for its triggering and closure (Grandpeix, Lafore & Cheruy 2010). Density currents affect clouds indirectly in two ways. First, they redistribute heat and water vapor vertically. Second, they play a role, via the term ALE wk , in triggering deep convection.
2.6. Shallow Convection 2.6.1. Thermal Plume Model and Shallow Cumulus Convection Version 6A uses a mass flux parameterization of thermals (Hourdin et al., 2002) instead of using a countergradient term in the vertical derivative of potential temperature and a dry convective adjustment as was the case in Version 5A (Hourdin, Foujols, et al., 2013). This thermal plume model was extended to the representation of shallow cumulus convection by Rio and Hourdin (2008). Conceptually, this model represents two subgrid-scale objects: a given coverage fraction of thermals, and their environment. The splitting technique mentioned in the previous section is also applied to the shallow convection scheme and thermals develop outside the cold pool region and in the same environment as the convective updrafts, that is, in a more unstable environment than that of the mean atmospheric grid cell. To do so, the potential temperature and total water content outside the cold pool region ( θ wk env and q wk t; env in Table 1) are used as inputs of the shallow convection scheme, thereby improving the buoyancy calculations and thermals development. In LMDZ6A, the thermal plume model was also improved by changing the detrainment formulation to better represent the transition from stratocumulus to cumulus clouds. This was done by using in the buoyancy formulation the difference in virtual potential temperature between the updraft and the environment at two different vertical levels, instead of computing the temperature difference on a same level. This method significantly improved the representation of clouds in regions of subsidence (for more details, see Hourdin et al., 2019).

Statistical Cloud Scheme
The shallow convection scheme is tightly connected to a statistical cloud scheme that uses a bi-Gaussian distribution Q of the saturation deficit s . The parameters required to compute the bi-Gaussian distribution are given by the thermal plume scheme and provided to the so-called large-scale condensation scheme described in the next section. In order to partly account for subgrid-scale temperature fluctuations, each Gaussian distribution is characterized by the mean saturation deficit and standard deviation of the thermal plume (s th and σ th ) and its environment (s env and σ env ), where the environment corresponds to the main mode of the bimodal distribution. The bi-Gaussian PDF can therefore be written as where α th is the coverage fraction of thermals and f is the classical Gaussian PDF: The in-cloud water content and cloud fraction can then be expressed as QðsÞ ds: The two mean saturation deficits s th and s env are computed automatically by the thermal plume model, and the variances are parameterized based on the coverage fraction of thermals α th (see Equations 7 and 8 of Jam et al., 2013). The shallow convection scheme does not remove the condensates from the prognostic total water variable at this stage, and only contributes to the mixing of q t (see Table 1). Shallow convective cloud formation and conversion to precipitation is computed afterward by the large-scale condensation scheme.

Large-Scale Condensation
The role of the large-scale condensation scheme is to condense the water vapor in excess of saturation coming from all the other procedures, as well as the water vapor brought to saturation by the large-scale horizontal circulation (which obviously affects q t and θ as well). It is in charge of the final calculation of the prognostic variables q l and q i and rebuilds the cloud macrophysical properties q in c and α c for further use in radiative transfer computations. It also computes the large-scale rainfall and snowfall rates P lsc l and P lsc i . Note that the term "large-scale" is a bit abusive in the sense that the cloud amounts and rainfall/snowfall rates computed by the large-scale condensation scheme include both large-scale clouds and shallow cumulus and stratocumulus clouds associated with the thermal plume model.
In practice, the large-scale condensation scheme computes, for each atmospheric column, the different processes using a vertical top-tobottom loop. In this section, the current layer will be referred to as z k , with z k + 1 the overlying layer and z k − 1 the underlying layer. The procedure computes all the condensed water contents in three steps: (1) It computes the reevaporation/sublimation of rain/snow coming from the overlying level z k + 1 (simply called reevaporation hereinafter), (2) it computes the amount of clouds that form in the gridbox at level z k using a subgrid-scale PDF, and (3) it converts part of the cloud into rain/snow. These three tasks are performed sequentially in this order. No structural changes were made to this scheme between Versions 5A and 6A, but many existing parameterizations were improved, and these adjustments will be noticed in the following subsections. 2.7.1.
Step 1: Reevaporation The loop starts with the reevaporation at level z k of the rain or snow coming from level z k + 1 . This reevaporation is based on the work by Klemp and Wilhelmson (1978) and Schlesinger et al. (1988) and can be written as where P l, i is the liquid or solid precipitation mass flux density in kg m −2 s −1 . It depends on the relative humidity q t /q sat and on a parameter called β, which is the same for rain and snow in LMDZ. Reevaporation is such that water vapor in the fraction of the gridbox below clouds does not exceed the saturation mixing ratio. In LMDZ5A, the reevaporation at level z k is limited to α ev c ðq sat − q t Þ, where α ev c ðz k Þ ¼ α c ðz k þ 1 Þ, with α c the actual cloud fraction simulated by the model (see the dashed line in Figure 1). This means that at two levels below cloud base, α ev c is set to 0 and reevaporation is no more possible. In LMDZ6A, α ev c was changed to the maximum cloud fraction found in the overlying layers and is reset back to 0 only if precipitation at level z k + 1 stops (see the solid line in Figure 1). This method implies that reevaporation is more efficient in Version 6A than in Version 5A (see the shaded gray area in Figure 1), if of course the value of the β coefficient in Equation 12 is unchanged. 2.7.2.
Step 2: Cloud Formation Cloud formation comes next, and the computation of the amount of condensates differs whether shallow convection is active in the gridbox or not. If shallow convection is active, cloud amount and fraction are computed using the bi-Gaussian PDF described in section 2.6. To do so, it uses the mean saturation deficits (s th , s env ) and standard deviations (σ th , σ env ) computed by the shallow convection scheme (Table 1). Otherwise, outside the grid cells where shallow convection is active, q in; lsc c and α lsc c are computed using a generalized lognormal PDF whose standard deviation σ is computed as σ ¼ ξq t . ξ is a function of pressure that has changed through the different versions of the model, as shown in Figure 2. In all versions, ξ is chosen so as to increase from the bottom of the troposphere to the top. Indeed, in the low and middle troposphere, the shallow convection scheme already computes the subgrid-scale water distributions and the large-scale standard deviation σ is therefore kept close to 0. In the case where the shallow convection scheme is not active, the standard deviation σ being close to 0, the scheme is almost equivalent to an "all-or-nothing" cloud scheme. The variance of the lognormal PDF in the lower and middle troposphere was set to a higher value in LMDZ5A than in LMDZ6A (using the ξ parameter represented in Figure 2) because the bi-Gaussian PDF was not implemented at the time and shallow convective clouds had to be represented by the lognormal PDF. In LMDZ6A, this becomes useless and the variance of the lognormal PDF is strongly reduced in the lower and middle troposphere to let the bi-Gaussian PDF of the shallow convection scheme do the calculation. In the high troposphere, ξ increases to reach a maximum value ξ 300 , which is used as a tuning coefficient. It exerts a strong control on the upper troposphere relative humidity and cloud cover (see section 3 of Hourdin, Grandpeix, et al., 2013).
Once q in; lsc c and α lsc c are computed, the cloud phase is distributed among liquid droplets and ice crystals according to temperature, resulting in some of the liquid droplets to be supercooled. The fraction of cloud water in the liquid-phase x liq is computed as where T min , T max , and n were set in Version 6A to −30°C, 0°C, and 0.5 respectively. As can be seen in Figure 3, the proportion of supercooled droplets was increased in LMDZ6A to be more consistent with the most recent satellite observations (Cesana & Chepfer, 2013;Cesana et al., 2015;Choi et al., 2014;Doutriaux-Boucher & Quaas, 2004). 2.7.3.
Step 3: Autoconversion Part of the cloud water is converted to precipitation, depending on cloud phase. For liquid clouds, this corresponds to a sink term that can be written as where τ conv is an autoconversion time constant and q clw is a threshold condensed water amount above which autoconversion sharply increases. Note that in Equation 14, q l is the liquid water mass mixing ratio within the gridbox, and that q l /α c is therefore the in-cloud liquid water content q in c . For ice clouds, the corresponding sink term follows: where q i is the water ice mass mixing ratio within the gridbox and where w iw ¼ γ iw w 0 . The fall velocity w iw depends on γ iw which is widely used as a tuning parameter of climate models Mauritsen et al., 2012). The terminal fall velocity is computed according (Heymsfield, 1977;Heymsfield & Donner, 1990), and depends on the mass of cloud ice without taking into account any actual size or shape of the particles. The conversion from cloud water to liquid or solid precipitation is done using a subtime step 5 times smaller than the physics time step. It is worth adding that in both

10.1029/2020MS002046
Journal of Advances in Modeling Earth Systems LMDZ5A and LMDZ6A, the cloud water content provided by the large-scale condensation scheme to the radiative transfer scheme is not what remains in the cloud at the end of the time step, but a mean cloud water content over the duration of the physics time step. Therefore, part of the cloud water that is converted to precipitation during the physics time step is "seen" by the radiative transfer scheme.
In Version 6A, the latent heat exchange due to the liquid ↔ ice phase change is not only implemented in the deep convection scheme (see section 2.3) but also in the large-scale condensation scheme. Moreover, when supercooled droplets are converted to precipitation, they freeze instantly, which was not the case in Version 5A. When freezing, rain releases latent heat, which can potentially bring the temperature back to above freezing. If this is the case, a small amount of rain remains liquid to stay below freezing. At the end of the large-scale condensation scheme, both the water vapor content q v and amount of condensates q l, i are known, as well as the in-cloud water content q in; lsc c and cloud fraction α lsc c provided by either the bi-Gaussian PDF used for shallow convection or generalized lognormal PDF used for large-scale condensation. The prognostic variables are ready for advection by the dynamical core, and the cloud water contents and fractions can be used by the radiative transfer scheme for heating rate calculations.

Radiative Transfer
Once the two cloud fractions α cv c and α lsc c are known, the total cloud fraction is estimated using where α lsc c includes both the cloud fraction coming from shallow convective clouds (bi-Gaussian PDF) and large-scale clouds (lognormal PDF), and where α cv c is the cloud fraction computed by the deep convection scheme. Similarly, the mean gridbox-averaged condensed water can be written as q rad is used in the radiative transfer scheme to compute the optical depth and α c is used to weight clear-sky and cloudy heating rates (precipitation is not radiatively active in LMDZ). The radiative transfer scheme uses the maximum random overlap assumption (Hogan & Illingworth, 2000;Morcrette & Fouquart, 1986). Cloud phase is determined using Equation 13. For liquid droplets, number concentration is parameterized using a modified version of Boucher and Lohmann (1995): where CDNC is the cloud droplet number concentration and m aer the soluble aerosol mass (instead of the sulfate aerosol mass used in Boucher & Lohmann, 1995, Equation D). Droplet sizes are then computed following Equations 2 and 4 of Boucher and Lohmann (1995). For ice crystals, particle sizes are parameterized following Equation 6 of Iacobellis and Somerville (2000) and vary in radius from r min at T < −81.4°C to 61 μm at 0°C (Heymsfield, 1986), where r min is a tuning parameter that varies between 3.5 and 20 μm. Note that aerosols have an impact on the size of the droplets but not on the size of the ice crystals. The first indirect effect of aerosols is therefore represented through the aerosol-dependent size of the droplets only. Liquid cloud radiative properties follow Fouquart (1988) and Smith and Shi (1992) in the SW and LW domain, respectively. Ice cloud radiative properties both in SW and LW domains are computed according to Ebert and Curry (1992). Aerosol radiative properties are computed as described in Lurton et al. (2020). LMDZ5A uses the Fouquart and Bonnel (1980) radiative transfer scheme in the SW (two bands) and LW domains (Morcrette, 1991), whereas LMDZ6A uses Fouquart and Bonnel (1980) only in the SW domain (and with six bands) and RRTM in the LW domain (Mlawer et al., 1997).

Summary of the Main Improvements
The changes affecting clouds made in Version 6A compared to Version 5A are therefore abundant and can be summarized as follows: • New scheme for local turbulent mixing (section 2.2); • New shallow convection scheme based on the so-called eddy-diffusivity-mass flux (EDMF) approach, coupled with the deep convection scheme; use of an improved statistical cloud scheme and bi-Gaussian PDF of the subgrid-scale distribution of the saturation deficit; new detrainment formulation (section 2.6); • New deep convection scheme that includes an improved mixing representation, new closure and a stochastic formulation of deep convection triggering (section 2.3); • New parameterization of cold pools coupled with the deep convection scheme; splitting technique applied to the grid cell to distinguish the cold pool region from its environment, and to allow both the shallow and deep convections to develop outside the cold pool region (section 2.5); • New vertical profile of the lognormal distribution's variance used for large-scale clouds ( Figure 2); • Inclusion of the latent heat exchange due to the liquid ↔ ice phase change in both the deep convection and large-scale condensation schemes; • New formulation of the subgrid-scale rain reevaporation rate ( Figure 1); • New phase-partitioning in mixed-phase clouds ( Figure 3); • New radiative transfer scheme (section 2.8).
2.10. Lessons Learned From the Development of LMDZ6A 2.10.1. Pros and Cons of a Multiobject Framework One of the most important aspects of LMDZ6A is the interplay between the different cloud parameterizations, that is, the shallow convection scheme, the deep convection scheme, and the so-called large-scale condensation scheme. The deep convection scheme forms a set of interconnected parameterizations that includes mixing, microphysics and thermodynamics. The representation of shallow convective clouds comes from two parameterizations, the thermal plume model and the large-scale condensation scheme. The thermal plume model transfers water from the surface to the cloud layer and provides the parameters of the subgrid-scale bi-Gaussian water distribution. The large-scale condensation scheme handles cloud formation and computes the reevaporation and autoconversion processes inside this newly formed cloud. This whole framework allows us to split into pieces complex processes and gradually link them together. It also enables the coupling between the different parameterizations, for example the deep convection triggering by thermals and cold pools (for more context on the state of the art of deep convection schemes, see Rio et al., 2019). However, each scheme provides its own cloud PDF, and ensuring a smooth transition between the different cloud PDFs is sometimes a difficult task.

Importance of Splitting the Grid Cell Into Two Regions
One key technical step was also distinguishing temperature and humidity inside and outside the cold pool region in both the shallow and deep convection schemes, so that both schemes run outside the cold pool region, in a more unstable environment than that of the mean atmospheric column. In Version 6A, both the thermal plumes and the deep convective updrafts thus develop in a same environment of given temperature and humidity, instead of using the mean grid cell values. Applying this splitting technique not only to the deep convection scheme (as was the case in some intermediate versions of the model) but also to the thermal plume model led to a strengthening of shallow convection relative to deep convection, and resulted in a major improvement in rainfall variability over tropical oceans. It also prevented the inhibition of shallow convection by deep convection, and that of deep convection by downdrafts and cold pools. This concept of splitting the atmospheric column in different subcolumns might be extended, in the future, to the boundary layer turbulence scheme and large-scale condensation scheme. It would allow the processes to affect temperature and humidity differently in the cloudy and clear portions of the cells. Adjusting the reevaporation rate was an essential part of the development of LMDZ6A. This rate is based on the fraction of overlying clouds (see Step 1 of section 2.7) but still affects the humidity of the whole gridbox. This splitting technique would make it possible to reevaporate rain only in the cloudy portion of the cell.

Revisiting Basic Thermodynamics
The development of LMDZ6A also revealed the importance of a consistent thermodynamics by the implementation of the heat exchange due to the liquid ↔ ice phase change and resulting changes in the entire cloud distribution. A disadvantage of a multiobject framework is the difficulty in ensuring thermodynamical consistency and energy conservation in the three different schemes.

Tuning as a Tool for Identifying Model Weaknesses
Finally, one essential lesson learned during the development of Version 6A is the need to tune the free parameters of the cloud schemes using well identified radiative targets. Beyond the technical need to tune climate models, tuning helps improve the physical formulations and identify model deficiencies "if parameter values 10.1029/2020MS002046 Journal of Advances in Modeling Earth Systems needed to satisfy a given metric are outside the acceptable range, or if different values are needed for different regions or climate regimes" . We will later see, for example, that the tuning of Version 6A revealed a probable deficiency in the computation of high-level cloud cover and associated overlap assumptions (see section 5). The tuning process is also a good way to reveal compensating errors.

Model Setup and Evaluation
The impact of the physics improvements described in section 2 on the cloud structure and properties is analyzed using two 20 year AMIP-typed simulations that are described in Table 2. We focus on the differences between Versions 5A and 6A of LMDZ, or more specifically between the atmospheric components of the IPSL-CM5A-MR and IPSL-CM6A-LR models, which share the same horizontal grid (144 × 142). However, we don't compare the two versions on the same vertical grid because the vertical resolution is strongly tied to the physical parameterizations of each version (39 levels for Version 5A and 79 levels for Version 6A). Thanks to the backward compatibility of LMDZ , the two simulations are run using the same source code, but the simulation is configured with LMDZ5A parameterizations in one case, and LMDZ6A parameterizations in the other. Version 5B is not analyzed in this paper because it was in many respects a prototype of Version 6A, as mentioned in section 1. The same aerosol concentration is used in both simulations, and is the one used for the CMIP6 project (Lurton et al., 2020). Both simulations are run using the most recent version of the ORCHIDEE soil and vegetation scheme. This scheme computes the vertical water transport in the soil using the Richard's equation (de Rosnay et al., 2002;D'Orgeval et al., 2008) discretized with 11 layers (see Cheruy et al., 2020, for more details on the scheme and its impact on the results of IPSL-CM6A).
The two simulations are tuned, meaning that some cloud parameters are adjusted . The tuning of LMDZ5A is described in section 3.4 of Hourdin, Grandpeix, et al. (2013), and the tuning of LMDZ6A is presented in Hourdin et al. (2020). When comparing the two simulations of the present paper, it is therefore important to keep in mind that the two simulations are tuned by targeting in particular a good TOA (top of atmosphere) global net flux. Some terms introduced in section 2 differ between LMDZ5A and LMDZ6A: ξ 300 in Figure 2, β in Equation 12, τ conv and q clw in Equation 14, γ iw in Equation 15 and r min (the smallest ice particle size) in section 2.8. The maximum precipitation efficiency for deep convection ep max is the same in the two simulations. The different values used for these parameters are summarized in Table 3. The role of each parameter in the tuning process is described in detail in Hourdin, Grandpeix, et al. (2013) and can be summarized as follows. Increasing β, τ conv or q clw tends to increase the amount of low-level clouds but impacts differently on their vertical profile. Increasing r min decreases the emissivity of high-level clouds. Increasing the γ iw coefficient increases the conversion to precipitation in ice clouds and decreases their water content. Increasing ep max decreases the amount of detrained water and high-level clouds in convective regions. As mentioned in section 2.7, the ξ 300 parameter has a strong impact on the relative humidity in the tropical upper troposphere and controls the variance of the lognormal PDF used in the cloud statistical scheme of high-level clouds. The latter three parameters (γ iw , ep max , and ξ 300 ) all affect the relative humidity of the tropical upper troposphere as they impact on the sources (ep max ) and sinks (γ iw and ξ 300 ) of water vapor.
Since the two simulations are tuned, both simulations correspond to the same mean climate state. Therefore, differences between the two simulations mainly arise from changes in the model parameterizations, and to a lesser extent from slight changes in the values of the tuning parameters themselves. The impact of the physics time step and the vertical resolution were also assessed using sensitivity experiments. Changing the physics time step from 30 to 15 min in Version 5A has almost no impact on the results. The vertical resolution, however, has a noticeable impact on the results (as also noticed in other models, e.g., Xie et al., 2018), and changing the number of vertical levels from 39 to 79 levels in Version 6A increases   Table 2 Tuning parameter LMDZ5A LMDZ6A ξ 300 (see Figure 2) (see Figure 2) β ((kg m −2 s −1 )

Cloud Spatial Distribution
We first compare the simulated cloud distribution to the lidar-based GOCCP data set (GCM Oriented CALIPSO Cloud Product, Chepfer et al., 2010). To do so, the cloud water contents and fractions predicted by LMDZ are processed by the CALIPSO-COSP simulator (Bodas-Salcedo et al., 2011) to derive the cloud fractions and covers the instrument would see if it was observing the model. To do so, the simulator uses the same overlap assumption as the LMDZ radiative transfer. Note that in the present paper, the term "cloud fraction" refers to the 3-D cloud fraction at each level and in each gridbox, whereas the term "cloud cover" refers to the total cloud cover seen from above, computed by integrating the 3-D cloud fractions vertically assuming a given overlap of clouds within the vertical column of the model gridboxes. This integral can be over the entire column or over a given pressure interval. In our case, we use three cloud covers that correspond to three cloud categories: low-level clouds (below 680 hPa or ∼3 km), mid-level clouds (between 680 and 440 hPa, that is, 3 and 6.5 km) and high-level clouds (above 440 hPa or ∼6.5 km). Figure 4 shows the cloud cover maps and bias maps of the three cloud categories, whereas Figure 5 shows the 3-D cloud fractions. Table 4 also summarizes the mean bias between the model and the observations, the RMSE and the correlation coefficient.
Starting with low-level clouds, comparing Figures 4a, 4d, and 4g reveals a significant improvement in the low-level cloud covers over the tropical oceans in LMDZ6A. On the west side of ocean basins, trade wind cumulus clouds were underestimated in Version 5A, as can be seen in Figures 4a and 4g. In LMDZ6A, they reach a better agreement with the observations (see Figures 4a and 4d). On the east side of ocean basins, stratocumulus clouds are improved in LMDZ6A due to the new statistical cloud scheme and change in the detrainment formulation of the thermal plume model (see section 2.6). Low-level clouds were underestimated over the Indo-Pacific warm pool in LMDZ5A and are now better represented as well. As can be seen in the bias plots Figures 4j and 4m, the overall bias is reduced in Version 6A over the tropical oceans but stratocumulus cloud cover maxima are slightly shifted away from the coasts. As described in Hourdin et al. (2019), this shift might be due to the tendency of the LMDZ6A model to maintain a 100% cloud deck for too long during the transition from stratocumulus to cumulus clouds. Outside the tropical belt, low-level clouds are overestimated over the Arctic and Southern Oceans. As evidenced in  Note. See Figure 4 for context.

Journal of Advances in Modeling Earth Systems
Figures 4j and 4m, this bias is stronger in Version 6A than in Version 5A. The overall RMSE for low-level clouds is reduced in Version 6A (see Table 5), mostly thanks to the improvements seen in the tropical regions.
The mid-level cloud distribution is one of the most striking improvements of LMDZ6A. A comparison of Figure 4b and 4e shows a reasonable agreement between the model and the observations, whereas previous versions of the model were systematically underestimating mid-level clouds. This is due to the improvement of the deep and shallow convection schemes in the tropical and mid-latitude regions (see sections 2.3 and 2.6), and to the new phase-partitioning of clouds in the mid-latitude to high-latitude regions (see Figure 3). As mentioned in section 3, the increase in vertical resolution from 39 levels in Version 5A to 79 levels in Version 6A also improved mid-level cloud covers in the ITCZ. High-level clouds are however underestimated in LMDZ6A, which was not the case before (Figure 4, right column). We had to reach a compromise in the tuning of the fall velocity parameter γ iw , which is relatively high in Version 6A (see Table 3). This tends to reduce the amount of high-level clouds globally to meet the LW CRE tuning target. Figure 5 shows the zonal mean cloud fractions averaged over 20 years of simulation in the two versions of the model and in the CALIPSO-GOCCP data set. As already noticed in Figure 4, outside the tropical belt, low-level clouds are overestimated in both LMDZ5A and 6A, but their altitude and fraction are improved in LMDZ6A. Their altitude of around 2 km is now slightly too high compared to the observations where low-level clouds are mostly below 1.5 km. Interestingly, comparing Figure 5e and 5h reveals that in Version 6A, we actually decreased the 3-D cloud fraction, but increased the geometrical thickness of low-level clouds, thereby increasing the low-level cloud cover (see Figure 4d). Mid-level clouds were mostly absent in LMDZ5A and are now better represented (see Figure 5e and 5h), especially over mid-latitude to high-latitude regions. This is also evidenced by the mean bias, RMSE and correlation coefficient shown in Table 4. In the tropics, LMDZ6A shows a local maximum in mid-level cloud cover slightly below 5 km altitude. The same maximum is located more than a 1,000 m higher in the observations, at elevations devoid of any cloud in the model. Another striking improvement of Version 6A is the water phase-partitioning in mid-level to high-level clouds. In LMDZ5A, the ice phase cloud fraction was clearly overestimated (see Figure 5i) and not consistent with the observations (Cesana et al., 2015). Changing the phase-partitioning in mixed-phase clouds (as shown in Figure 3) significantly improved the ice phase cloud fractions in LMDZ6A ( Figure 5, right column), as well as the liquid-phase cloud fractions in mid-level clouds (middle column). As previously mentioned, high-level cloud cover remains underestimated due to a compromise in the tuning of the model, but their spatial distribution is improved (see Figure 4f and correlation coefficients in Table 4). High-level 3-D cloud fractions are overestimated in the tropical regions if we compare Figures 5c and 5f, but their total column cloud cover is underestimated in this same region if we look at Figure 4f and upper left panel of Figure 8. This suggests, as will be discussed in section 5, that the cloud cover computed by the model for high-level clouds is too low and compensated by a too high 3-D cloud fraction. Figure 6 focuses on the cloud fraction in the tropical regions, more exactly on the GPCI transect, which spans from San Francisco to Honolulu (see Teixeira et al., 2011, for more details). This transect is especially useful Note. See Figure 7 for context. to evaluate the representation of the stratocumulus to cumulus (Sc-to-Cu) and shallow to deep convection transitions in climate models. In LMDZ5A, the Sc-to-Cu transition was visible, but stratocumulus clouds were too close to the surface and high-level cloud fractions were overestimated. Version 6A nicely represents the Sc-to-Cu transition and shows a better evolution of the cloudy boundary layer, but clouds tend to extend beyond the 2 km height seen in the observations. Over the warmer waters of the trade wind boundary layer (around 5°N), the model cloud fractions remain too low compared to the observations. Mid-level cloud fractions are also underestimated in deep convective regimes, as is also noticed in Figure 5d at around 8 km altitude. This altitude range is where the ξ(p) function sharply increases (see Figure 2). It is therefore in the transition zone between the PDFs of the shallow convection, deep convection and large-scale condensation schemes, and suggests that the interplay between the schemes need to be improved in this region. The high-level cloud fraction is better represented in LMDZ6A but clouds remain too geometrically thin compared to the observations.

Cloud Radiative Effect
Clouds play a crucial role in the radiative budget of the atmosphere, and a compromise has often to be found between a good representation of their properties and a good TOA energy budget of the model. The tuning method of LMDZ6A is described in Hourdin et al. (2020), and we focus here on the role of clouds in the radiative budget. Figure 7 shows the observed and simulated Cloud Radiative Effect (CRE) in the SW and LW domains, as well as the bias maps. The left column of this figure shows a clear improvement of the SW CRE, especially in mid-latitude to high-latitude regions where reflection by low-level clouds was too high in Version 5A. This improvement results in a 5 W m −2 reduction of the SW CRE mean bias and RMSE in LMDZ6A, as shown in Table 5. An improvement of the same magnitude is seen in the LW CRE, but in this latter case, the spatial distribution is also improved (see the increase in the correlation coefficient in Table 5), which is less the case of the SW CRE, especially in the tropical regions. Indeed, despite a clear improvement of the SW CRE in the ITCZ (see Figure 7c), the SW radiative effect of stratocumulus clouds is shifted away from the coast over the eastern part of tropical ocean basins, and trade wind cumulus clouds reflect less sunlight than in the observations (see Figure 7g). These biases are consistent with those of the low-level cloud cover described in section 4.1. The left column of Figure 8 summarizes the zonal mean cloud cover of the three cloud categories and the corresponding radiative forcings in the right column. In the tropics, cloud covers are improved in LMDZ6A at all levels, but remain slightly lower than in the observations. The right column of Figure 8 shows that in this region a realistic CRE is reached even though the cloud covers are slightly biased low. For high-level clouds, this suggests that the underestimated cloud covers are probably compensated by a too high 3-D cloud fraction. For low-level clouds, it suggests that the underestimated cloud cover is compensated by overly bright low-level clouds, as will be discussed in section 5. The situation is different over the Arctic and Southern oceans, where a realistic CRE is reached even though the low-level cloud covers are biased high (see Figures 7c,7d,and 8,lower left panel). In these regions, the LMDZ6A SW CRE is in better agreement with the observations than that of LMDZ5A, and this has to do with cloud phase and opacity, as we will see in the next paragraph. It is worth noting that this difference in low-level cloud covers between the two versions could have come from the results of the simulator because of the possible screening of low-level clouds by high-level clouds. In our case, high-level cloud covers are biased low relative to the observations in Version 6A (see Figure 4l) and could increase the signal coming from low-level clouds and partly explain the positive cloud cover bias seen in Figure 4d. But this is not the case. We find the same difference between the low-level cloud covers of Versions 5A and 6A using the results of the model radiative transfer itself (not the simulator).
Let's now return to the good total CRE simulated in LMDZ6A in mid-latitude to high-latitude regions despite the biases seen in the various cloud covers (Figure 8). In mid-latitude to high-latitude regions, phase-partitioning has been found to be strongly connected to the SW CRE in many models (McCoy et al., 2016). In our case, sensitivity experiments show that increasing the temperature range of supercooled droplets leads to a greater vertical extension of liquid clouds, which are otherwise confined to lower layers. This results in a higher concentration, in LMDZ6A, of liquid droplets in mid-level clouds, where droplets are more reflective than ice (Liou, 2003), but more importantly in a lower concentration of droplets in low-level clouds. This decrease in the concentration of liquid droplets in low-level clouds explains why the SW CRE is in better agreement with the observations in LMDZ6A, despite the overestimation of the low-level cloud cover. The LW CRE is also sensitive to phase-partitioning in mixed-phase clouds. The left column of Figure 8 shows that LMDZ6A has less high-level clouds and more mid-level clouds in mid-latitude to high-latitude regions. Decreasing the high-level cloud cover decreases the LW CRE, but on the other hand, the increase in mid-level cloud covers of high liquid content strongly increases it. In the end, the LW CRE in LMDZ6A is reduced by the right amount compared to that of LMDZ5A and is in good agreement with the

10.1029/2020MS002046
Journal of Advances in Modeling Earth Systems observations. The overall cloud liquid water path in mid-latitude to high-latitude regions is increased, as illustrated in Figure 9, whereas the ice water path is strongly decreased. Satellite retrieval of the LWP and IWP is not an easy task, but a comparaison of the simulated LWP with the work of O'Dell et al. (2008) suggests that it is in good agreement with the observations in the tropical regions and slightly too high in mid-latitude to high-latitude regions. It is more difficult to compare the simulated cloud IWP (i.e., the nonprecipitating ice) to existing observations, but the sharp decrease in the cloud IWP of Version 6A is more in line with the cloud IWP found in other models, including the ERA5 and MERRA-2 reanalyses (see Figure 3 of Duncan & Eriksson, 2018). Figure 10 focuses on the tropics and shows the simulated CRE as a function of the dynamical regimes (through the vertical velocity ω at 500 hPa). This type of analysis, introduced by Bony et al. (2004), shows how well the CRE is represented in regions of subsidence (ω > 0) and updraft (ω < 0). Both the SW and LW CREs show a gradual decrease (in terms of absolute value for the SW CRE) from regions of strong updrafts where clouds are abundant to regions of strong subsidence where clouds dissipate. The lower panel of Figure 10 shows a clear improvement of the total CRE in LMDZ6A in convective regions (ω < 0). This is mostly due to an improvement of the SW CRE (upper panel), and reflects the changes applied to the thermal plume parameterization, which improved both the stratocumulus clouds over the eastern part of tropical ocean basins and trade wind cumulus clouds (see section 4.1). However, the SW and LW CREs are still too weak in magnitude in strongly convective regions (ω < −40 hPa/day) and the SW CRE is higher than observed in regions of strong subsidence (ω > 20 hPa/day).

Discussion
Thanks to the improvements of the physical parameterizations and to an experience gained in the tuning of the model, the cloud distribution and radiative effects have been significantly improved in LMDZ6A. But the refined tuning of the model has also underlined structural problems, especially in the detailed cloud radiative properties. In particular, the difficulty to tune high-level clouds points to an inappropriate representation of their radiative properties, which impacts on all clouds. Difficulties in modeling the properties of high-level clouds were already found in the early versions of LMDZ (Webb et al., 2001). Figure 11 shows the PDF of the high-level cloud cover over the tropical oceans based on the daily outputs of the CALIPSO-GOCCP observations (left panel) and results of the LMDZ model simulator (middle and right panels). The observed PDF is a highly skewed-right distribution with a peak at 0-5% cloud cover and an outlier at 97.5-100%. The LMDZ5A PDF shows a lower peak at 0-5% but an otherwise similar distribution, with a smaller outlier at 97.5-100%. LMDZ6A shows a skewed-right distribution similar to the observations for cloud cover lower than 20%, but its PDF differs significantly for higher cloud covers, with a strong decrease above 50% and no outlier at 97.5-100%. This difficulty of LMDZ6A to attain complete coverage for high-level clouds might explain why these clouds are hard to tune in this version. Therefore, work is underway to improve the ξ function (see section 2.7) using a more physical parameterization, as well as the overlap assumptions and subgrid-scale heterogeneities of high-level clouds.
Regarding tropical low-level clouds, Figure 12 shows the density of points of a given cloud reflectance and cloud cover in the observations (left panel) and in the model (middle and right panels, see Konsta Figure 10. Regime sorted plots of the SW (top), LW (middle), and net (bottom) CRE as a function of ω500 in hPa/day between 30°S and 30°N and over the oceans. For comparison, the black line shows the same diagnostics obtained using ERA reanalysis and the CERES data (EBAF data set Loeb et al., 2009). et al., 2016, for more details on the method). Cloud reflectance in Figure 12 is a function of the vertically integrated cloud optical depth, whereas cloud cover will be more dependent on the cloud fraction vertical profiles and overlap assumption. Two populations can be identified in the observations (Figure 12, left panel): Trade wind cumulus clouds have low reflectance and cover values, whereas stratocumulus clouds have medium reflectance and high cover values. The observations also show an increase in cloud reflectance with increasing cloud cover. LMDZ5A was showing the opposite tendency ( Figure 12, middle panel) and trade wind cumulus clouds were too bright in this version of the model, a problem commonly referred to as the "too few, too bright" problem (Nam et al., 2012). As explained in Konsta et al. (2016), this increase in reflectance with decreasing cloud fraction in LMDZ5A was due to the activation of the deep convection scheme in trade wind regions, which affected the low-level cloud PDFs. The implementation of the thermal plume model in LMDZ6A clearly improved the distribution, which is now closer to the observations (Figure 12, right panel). However, in LMDZ6A, trade wind cumulus clouds are still too reflective and their cover is too low. Stratocumulus clouds are well represented and show medium reflectance and high cover values, in agreement with the observations. Between these two populations, a third population appears in the model, and is characterized by cloud reflectance values of around 0.2 and cover values between 0.6 and 0.9. The too few, too bright bias was thus reduced but not fully solved. Despite the high number in LMDZ6A of low cloud cover values compared to the observations (Figure 12, right panel), the SW CRE is still in good agreement with the observations. This suggests that this too low cover is compensated by an excessive brightness in the tuning process, which targets the CRE as a  priority. We thus see cloud reflectances of around 0.3 in LMDZ6A, compared to less than 0.1 in the observations (see Figure 12, left and right panels). This shows the limit of the maximum random overlap assumption used in LMDZ6A. Preliminary sensitivity experiments performed with LMDZ6A shows that using the exponential-random overlap assumption (Hogan & Illingworth, 2000) instead of the maximum random overlap assumption may improve the distribution shown in Figure 12 by increasing cumulus cloud cover. Another way to increase low-level cloud covers is to represent subgrid-scale vertical heterogeneities by distinguishing the cloud fraction by volume from the cloud fraction by surface. The latter was found to be 20% greater on average than the cloud fraction by volume (Brooks et al., 2005). The cloud fraction by surface is more appropriate for coupling with radiative transfer schemes but most climate models do not yet distinguish between the two quantities and by doing so, assume that the cloudy area of a gridbox fills the entire gridbox in the vertical. The difference between the cloud fraction by volume and the cloud fraction by surface can be computed by a parameterization of subgrid-scale heterogeneities that will depend on the vertical resolution and various physical information, such as wind shear (Sulak et al., 2020). Work is underway to implement such parameterization in LMDZ (Jouhaud et al., 2018). This could improve the CRE of low-level clouds but also high-level clouds.

Conclusion
After a series of parameterization changes (summarized in section 2.9) and a finer tuning of the radiative budget, several cloud features were improved in Version 6A of the LMDZ climate model: • Low-level (below 3 km) cloud covers are improved both in trade wind regions and in the east side of ocean basins (see Figure 4d), due to the new shallow convection scheme; • Mid-level clouds, which were almost inexistent in LMDZ5A, are much better represented in LMDZ6A (see Figure 4e). Mid-level to high-level cloud phase is also more realistic and now includes a more realistic fraction of supercooled liquid droplets (see Figure 5e). These improvements mostly come from the changes made in the deep convection scheme in the tropical regions, and in the new phase-partitioning in the midlatitude to high-latitude regions; • CREs are improved (see Figure 8, right column) and LMDZ6A shows a 5 W m −2 improvement in both the SW and LW CRE compared to LMDZ5A (see Table 5), due to the combined effect of the new shallow convection scheme and new phase-partitioning; • A 20 W m −2 bias in the SW CRE of the convective regions is corrected (see Figure 10, upper panel), thanks mostly to the new shallow convection scheme; • Tropical low-level cloud reflectance and cover are significantly improved (see Figure 12, right panel) due to the shallow convection scheme and its new statistical cloud scheme based on a Bi-gaussian PDF.
The finer model tuning performed for LMDZ6A also revealed structural errors and inconsistencies that call for a revisit of some existing parameterizations. Indeed, the model reaches a good radiative balance for cloud covers that are sometimes strongly biased. This is true for low-level clouds but more importantly for high-level clouds, whose covers need to be lower than observed to restore the radiative balance. For clouds of all levels, work is underway to improve the overlap assumptions of the radiative transfer scheme and to better account for the cloud subgrid-scale heterogeneities (see, e.g., Jouhaud et al., 2018). High-level clouds also rely on a fixed value of the lognormal PDF variance (ξ 300 ) which must be improved and more physically based. Mid-level clouds are also the focus of current development efforts, in order to better represent the deepening of shallow cumulus clouds into congestus clouds (see Figure 6). Improvement of the cloud microphysical scheme is also underway, with a particular focus on cold and mixed-phase clouds. Priorities include the improvement of the conversion of ice clouds to solid precipitation (Lemonnier et al., 2020), the implementation of supersaturation with respect to ice , and the representation of subgrid-scale processes in mixed-phase clouds.
Appendix A: Notations ρ Atmospheric density (kg m −3 ) ω500 Large-scale vertical velocity at 500 hPa (hPa day −1 ) θ Potential temperature (K) q v Water vapor mass mixing ratio (kg kg −1 ) q l Liquid water mass mixing ratio (kg kg −1 ) q i Ice mass mixing ratio (kg kg −1 ) q t Total water mass mixing ratio (kg kg −1 ) q tc Gridbox mean amount of condensate and in-cloud vapor (kg kg −1 ) q sat Saturation mass mixing ratio (kg kg −1 ) s Saturation deficit (see Equation 3 of Jam et al., 2013) (kg kg −1 ) P(q) Probability Density Function (PDF) of water vapor q Q(s) Probability Density Function (PDF) of the saturation deficit s ALE Available Lifting Energy (J kg −1 ) ALP Available Lifting Power (W m −2 ) w iw Fall velocity of ice crystals (m s −1 ) w 0 Terminal fall velocity of ice crystals (m s −1 ) P l, i Liquid/Ice precipitation flux density (kg m −2 s −1 ) dθ cv dw Temperature tendency due to downdrafts (K s −1 ) dq cv t; dw Total water tendency due to downdrafts (kg kg −1 s −1 ) α th Coverage fraction of thermals θ env θ in the environment of the plume (K) q t, env Mean q t in the environment of the plumes (kg kg −1 ) s env Saturation deficit in the environment of the plumes (kg kg −1 ) σ env σ of the PDF related to the environment of the plumes (kg kg −1 ) s th Saturation deficit inside the plumes (kg kg −1 ) σ th σ of the PDF related to the plumes (kg kg −1 ) q in c In-cloud water mass mixing ratio (kg kg −1 ) α c Cloud fraction q m Condensed water mixing ratio in the mixed drafts (kg kg −1 ) M t Mass flux density of the mixed drafts (kg m −2 s −1 ) α m Coverage fraction of mixed drafts τ m Dissipation time constant of the saturated drafts (s) δz Vertical spacing of gridboxes (m) M a Mass flux density of the undiluted updrafts (kg m −2 s −1 ) α a Coverage fraction of undiluted updrafts w a Vertical velocity of the undiluted updrafts (m s −1 ) When written in superscript, th, wk, cv and lsc indicates variables related to thermal plumes, wakes, deep convection, and large-scale condensation, respectively.
For a list of the tuning parameters and their notations, see Table 3.

Data Availability Statement
The last version of the LMDZ source code can be downloaded freely on the LMDZ website (https://lmdz. lmd.jussieu.fr). The version used for the specific simulation runs of this paper is the svn release 3404 from 16 October 2018