A unified parameterization of clouds and turbulence using CLUBB and subcolumns in the Community Atmosphere Model

Most global climate models parameterize separate cloud types using separate parameterizations. This approach has several disadvantages, including obscure interactions between parameterizations and inaccurate triggering of cumulus parameterizations. Alternatively, a unified cloud parameterization uses one equation set to represent all cloud types. Such cloud types include stratiform liquid and ice cloud, shallow convective cloud, and deep convective cloud. Vital to the success of a unified parameterization is a general interface between clouds and microphysics. One such interface involves drawing Monte Carlo samples of subgrid variability of temperature, water vapor, cloud liquid, and cloud ice, and feeding the sample points into a microphysics scheme. This study evaluates a unified cloud parameterization and a Monte Carlo microphysics interface that has been implemented in the Community Atmosphere Model (CAM) version 5.3. Model computational expense is estimated, and sensitivity to the number of subcolumns is investigated. Results describing the mean climate and tropical variability from global simulations are presented. The new model shows a degradation in precipitation skill but improvements in shortwave cloud forcing, liquid water path, long-wave cloud forcing, precipitable water, and tropical wave simulation.


Introduction
Most climate models today use separate parameterizations to model separate cloud types, such as stratiform clouds, shallow cumuli, and deep cumuli.Each parameterization uses its own separate equation set.The resulting suite of parameterizations is intended, collectively, to represent the full range of subgrid-scale clouds included in the climate model.
While the use of separate parameterizations for separate cloud regimes offers several advantages, it also suffers disadvantages.First, the use of multiple, separate cloud parameterizations leads to unnecessary complexity.Some of the complexity is of a practical sort: it is hard to understand a suite of parameterizations written by different authors that use differing coding conventions and assumptions.Some of the complexity is more conceptual in nature: even if each parameterization is simple, the interactions among the parameterizations might be complex (Zhang and Bretherton, 2008;Bretherton, 2007).Second, it is difficult to formulate, in a realistic way, the triggers that are used to activate cumulus parameterizations.For instance, deep convection does not appear instantaneously; rather, in many instances, deep clouds are initiated by the gradual and continuous growth of shallow clouds (Grabowski et al., 2006;Wu et al., 2009).Accurately parameterizing the gradual onset of deep convection is important for modeling tropical phenomena such as the Published by Copernicus Publications on behalf of the European Geosciences Union.
To avoid such difficulties, some past researchers have parameterized two or more cloud types using a single equation set, thereby partly unifying the description of clouds.The greater the degree of unification, the greater the reduction in the number of interacting parameterizations and trigger functions.
For instance, to avoid the difficulties of coupling shallow and deep cumulus parameterizations, some researchers have represented both cloud types using a single parameterization (Kain, 2004;Park, 2014a, b).However, the aforementioned parameterizations are only partly unified because they do not include stratiform clouds; instead, those clouds must be handled by a separate parameterization.
To avoid the difficulties of coupling stratocumulus and shallow cumulus parameterizations, some researchers have parameterized both cloud types with a single equation set (Lappen and Randall, 2001;Golaz et al., 2002;Larson and Golaz, 2005;Cheng andXu, 2006, 2008;Firl, 2009;Bogenschutz and Krueger, 2013).To close some higher-order terms in the equation set, these parameterizations make an assumption about the shape of the probability density function (PDF) of subgrid variability.Assumed PDF parameterizations have a long history in atmospheric science (e.g., Manton and Cotton, 1977;Sommeria and Deardorff, 1977;Mellor, 1977;Bougeault, 1981a, b;Lewellen and Yoh, 1993).For several decades, PDF parameterizations have been implemented in regional or global models (e.g., Smith, 1990;Tompkins, 2002;Nakanishi and Niino, 2004).Recently, the Cloud Layers Unified By Binormals (CLUBB) parameterization has been implemented and evaluated in two global climate models (Bogenschutz et al., 2013;Guo et al., 2014).In these implementations, CLUBB does unify the representation of boundary layer clouds, but both implementations parameterize deep convection separately.Guo et al. (2015) used a similar configuration to Guo et al. (2014), but also used CLUBB to parameterize deep clouds.However, this configuration does not parameterize, in a unified way, subgrid variability in ice clouds.
The configurations used by Bogenschutz et al. (2013), Guo et al. (2014), and Guo et al. (2015) share three drawbacks.First, none of those three configurations fully unifies the description of all cloud variability because in all three configurations, cloud ice is not seen by CLUBB.Specifically, cloud ice is not included in CLUBB's subgrid PDF.Second, even for liquid clouds, the description is, in certain respects, internally inconsistent.For instance, a different marginal PDF shape of cloud water is assumed by CLUBB in order to diagnose cloud liquid water content (namely, a truncated normal mixture) than is assumed by the microphysics in order to compute autoconversion and accretion (namely, a gamma function) (Morrison and Gettelman, 2008).(A uni-variate marginal PDF is the PDF that remains when a multivariate PDF is integrated over all variates but one.)Third, certain aspects of the subgrid variability, such as the precipitation fraction, are treated by a microphysics scheme that is designed to parameterize stratiform cloud (Morrison and Gettelman, 2008) and whose assumptions about subgrid variability may not be well-suited to cumulus clouds.These three drawbacks might be related to certain errors seen in the simulations, such as the overestimate of precipitable water and underestimate of cloud ice noted by Guo et al. (2015).
One key to parameterizing deep convection is accurately parameterizing the subgrid coupling between clouds and microphysics (Emanuel, 1991;Donner, 1993).The reason is that interactions among condensed water content, clear-air relative humidity, and precipitation evolution are strong.In fact, Hohenegger and Bretherton (2011) stated that "the main difference between shallow and deep convection is precipitation (both rain and snow) and its effects."If true, this hints that a PDF parameterization that can accurately parameterize shallow convection can, in conjunction with a suitable coupling to the microphysics, also parameterize deep convection.
Here, in order to interface clouds and microphysics, we use a Monte Carlo integration technique named the Subgrid Importance Latin Hypercube Sampler (SILHS) (Larson et al., 2005;Larson and Schanen, 2013).SILHS samples the subgrid PDFs predicted by CLUBB, thereby providing a set of vertical profiles, or subcolumns, of sample points.The subcolumns are then fed into a single microphysics scheme, thereby allowing the microphysics to respond to subgrid variability in clouds (Jakob and Klein, 1999;Jess et al., 2011;Tonttila et al., 2013Tonttila et al., , 2015)).Within an individual subcolumn, each grid level has uniform properties (e.g., all cloudy or all clear), but collectively, a set of subcolumns represents the subgrid variability within a grid column.This may improve the representation of nonlinear microphysical process rates (Pincus and Klein, 2000;Larson et al., 2001;Jess et al., 2011).Subcolumn approaches have long been used for radiative transfer applications in large-scale models (e.g., Barker et al., 2002Barker et al., , 2008;;Pincus et al., 2003Pincus et al., , 2006;;Räisänen et al., 2004Räisänen et al., , 2005Räisänen et al., , 2007Räisänen et al., , 2008;;Räisänen and Barker, 2004).
The use of SILHS helps mitigate the three aforementioned drawbacks of the configurations of Bogenschutz et al. (2013) and Guo et al. (2014Guo et al. ( , 2015)).First, cloud ice is included in CLUBB's subgrid PDF and is sampled by SILHS, thereby driving ice microphysics with subgrid variability.Second, SILHS feeds within-cloud variability directly and consistently into microphysics, ensuring that the same marginal PDF that is used to diagnose cloud water content is also used to diagnose autoconversion.Third, assumptions about subgrid variability, such as those regarding vertical overlap of condensate and vapor, are removed from the microphysics scheme and instead embedded in SILHS (Larson and Schanen, 2013;Storer et al., 2015).This facil-itates the implementation of subgrid assumptions that are more general.
Here, we evaluate a new configuration of the CAM climate model that we call CAM-CLUBB-SILHS.It shuts off the Zhang and McFarlane (1995) parameterization of deep convection and instead uses CLUBB to parameterize deep cumulus, shallow cumulus, stratiform liquid clouds, and stratiform ice clouds.SILHS is used in order to feed samples of the subgrid variability into a microphysics scheme, following the approach of Storer et al. (2015).A single microphysics scheme is used in all cloud types.This model configuration provides a more fully unified parameterization of clouds.The purpose of the present paper is 2-fold.First, it outlines the subcolumn software framework in CAM.This software framework contains SILHS.Second, unlike Storer et al. (2015), this paper evaluates the behavior of CLUBB-SILHS in a global context, including climatologies of cloud-related fields and some aspects of tropical variability.
This paper is organized as follows.Section 2 describes the CLUBB-SILHS methodology and its implementation in CAM.Section 3 estimates the computational cost of CAM-CLUBB-SILHS.Section 4 evaluates the mean climate versus satellite observations.Section 5 evaluates CAM-CLUBB-SILHS' simulation of tropical variability.Section 6 illustrates the sensitivity to the number of subcolumns.Section 7 summarizes the evaluation and concludes.

Description of the CLUBB moist turbulence parameterization
CLUBB's methodology is described in Golaz et al. (2002), and an up-to-date listing of CLUBB's equations is contained in Storer et al. (2015).CLUBB parameterizes subgrid turbulence in both clear and cloudy air, and subgrid variability in all cloud types, including stratiform, shallow cumulus, and deep cumulus.If CLUBB's single equation set is to represent turbulence and all cloud types, the equation set must be sufficiently rich and general.CLUBB's equation set includes prognostic equations for various moments of the vertical air velocity w, the liquid water potential temperature θ l , and total water mixing ratio (vapor plus liquid cloud water) r t .The grid-averaged means of these variables are prognosed by the host model, CAM.CLUBB adds prognostic Reynoldsaveraged equations for the following moments: w θ l , w r t , w 2 , w 3 , r t 2 , θ l 2 , r t θ l (Golaz et al., 2002;Larson and Golaz, 2005).CLUBB parameterizes momentum fluxes using down-gradient diffusion, but CLUBB does not explicitly parameterize subgrid-scale mesoscale convective organization (e.g., Moncrieff, 1992;Donner, 1993;Moncrieff and Liu, 2006).
These prognostic equations include several higher-order moments that are unclosed.To close them, CLUBB integrates them over a PDF of subgrid variability.CLUBB contains a multivariate subgrid PDF for r t , θ l , w, cloud ice (mass) mixing ratio r i , and cloud ice number mixing ratio N i .The inclusion of r t and θ l allows both moisture and temperature fluctuation to enter the diagnoses of cloud fraction and cloud water mixing ratio.The inclusion of w allows the buoyancy flux, w θ v , to be computed consistently with cloud fraction and cloud water.The inclusion of ice in the PDF allows ice processes to be coupled to the drafts and thermodynamics on the subgrid scale.The marginals of w, r t , and θ l are normal mixtures, that is, the sum of two Gaussians.This PDF shape has been shown to compare favorably with aircraft observations and large-eddy simulations of stratiform, shallow cumulus, and deep cumulus clouds (Larson et al., 2002;Bogenschutz et al., 2010).The marginal PDF for r i and N i is a delta double lognormal.That is, the PDF shape for ice is the sum of a delta function representing the ice-free area and the sum of two lognormal distributions.This PDF shape has recently been evaluated against large-eddy simulations (Griffin and Larson, in preparation).The within-ice standard deviation of r i is assumed to be proportional to the within-cloud mean (Lebo et al., 2015).The same is true for N i .The correlations among hydrometeors -including mass and number mixing ratios of liquid and ice -are prescribed as in Storer et al. (2015).

The interface between clouds and microphysics: SILHS
CLUBB computes the transport of hydrometeors and production of cloud water via saturation adjustment, but CLUBB must be coupled to a microphysics scheme in order for other microphysical process rates to be computed.The coupling between clouds and microphysics is accomplished by use of a Monte Carlo sampler called SILHS.SILHS' methodology is described in Larson et al. (2005) and Larson and Schanen (2013).SILHS draws n samples from the subgrid PDF at each grid level.When the liquid cloud fraction is moderate, half the samples are drawn from liquid cloud and half are drawn from the remainder of the grid box, with appropriate weighting, using the method described in Larson and Schanen (2013).The n samples at each grid level are used to construct n vertical profiles of sample points, or subcolumns.
In order to parameterize cloud overlap, non-zero vertical correlation between vertical grid levels is allowed.The vertical correlation between samples is assumed to drop off exponentially with vertical distance (Larson and Schanen, 2013).Each subcolumn is fed into Version 1.0 of the Morrison-Gettelman (MG1) microphysics scheme (Morrison and Gettelman, 2008).MG1 provides a simplified initial test for the subcolumn methodology because MG1 diagnoses rain and snow.Therefore, rain and snow are not inputs to MG1, and hence the subcolumns need not contain rain or snow vari- ates.In the future, we hope to use SILHS with version 2.0 of Morrison-Gettelman (MG2) microphysics scheme (Gettelman and Morrison, 2015;Gettelman et al., 2015).MG2 prognoses rain and snow, and hence using it will require us to draw subcolumns with rain and snow.Although this will add complexity and expense, the higher-dimensional PDF will offer greater control over processes that involve two or more hydrometeor species, such as accretion of cloud water by rain water.When subcolumns are used, MG1's native assumptions about subgrid variability, including a gamma distribution of cloud water, are shut off, and MG1 is made to assume that each grid level has uniform properties, e.g., is overcast or clear.MG1 calculates time tendencies for cloud ice, cloud liquid water, water vapor, and other relevant microphysical variables.One set of microphysical tendencies is calculated per each subcolumn.The tendencies are then averaged in order to produce a grid-mean tendency.The grid-mean tendencies are then fed into the host model's grid-mean equations for microphysical species, temperature, and moisture.The averaging is weighted appropriately to account for the fact that different subcolumns may represent different-sized areas of a grid column, as described in Larson and Schanen (2013).
Ice processes are coupled to CLUBB's grid-mean thermodynamical variables, θ l and r t , through the microphysics.Subcolumns that include subgrid variability in vapor, liquid, and ice are fed into the microphysics, and the effects of ice, such as the Bergeron effect, are computed by the microphysics at the subgrid scale.These effects of ice are expressed in terms of microphysical tendencies of vapor, liquid, and ice.These tendencies are used to update θ l , r t , and r i .These updated values influence ice during the subsequent time step.In this sense, ice and liquid processes interact on the subgrid scale.Although information about the subgrid PDF of ice is contained within CLUBB, SILHS is needed in order to carry out the subgrid (Monte Carlo) integration of complex, nonlinear ice microphysical processes.
Although CLUBB is substepped with a 5 min time step, MG1 is called with a 30 min (physics) time step.At each physics time step, new SILHS sample points are drawn from CLUBB's PDF from CLUBB's most recent substep.The subcolumn-averaged microphysical tendencies are fed back into the host model at the end of the physics time step.SILHS retains no memory of sample points from one time step to the next.Rather, the memory is retained within CLUBB's prognosed moments.

Comparison of CLUBB-SILHS with other modeling techniques
Now that CLUBB-SILHS' methodology has been described, we pause and briefly contrast CLUBB-SILHS with other methods.
First, we compare and contrast CLUBB-SILHS with the eddy-diffusivity mass-flux (EDMF) approach (e.g., Soares et al., 2004;Siebesma et al., 2007;Neggers et al., 2009;Neggers, 2009;Sušelj et al., 2012Sušelj et al., , 2013Sušelj et al., , 2014)).Broadly speaking, two types of grid-box averaging ought to be performed, explicitly or implicitly, in large-scale models: (1) grid averaging of subgrid turbulent fluxes, and (2) grid averaging of source terms, such as microphysical tendencies.Whereas CLUBB prognoses the turbulent fluxes of moisture and heat content based on the parameterization of each individual term in the flux budget, EDMF diagnoses those turbulent fluxes based on physical considerations.Whereas CLUBB-SILHS averages microphysical tendencies by Monte Carlo integration, EDMF per se delegates the averaging of those tendencies to other parameterizations.CLUBB-SILHS is more expensive than EDMF, but CLUBB-SILHS' foundation in PDFs facilitates the consistent calculation of, e.g., cloud fraction and virtual potential temperature, and allows the global use of a single microphysics scheme for all clouds.
Second, we distinguish CLUBB-SILHS from methods that alter the grids on which the equations are solved.We consider two such example methods.One is the multiscale modeling framework (MMF; e.g., Grabowski, 2001).It embeds a convection-permitting model within each grid column of a climate model, thereby unifying the description of cloud features larger than about 4 km in the horizontal extent.Another is the method of Yano et al. (2005), which spectrally decomposes the equations into wavelet modes, and thereby unifies the description of those cloud features that are resolved by the wavelet models.These two methods are more akin to nested gridding or variable-resolution gridding techniques than to parameterizations such as CLUBB.These two methods have the advantage of containing information about the horizontal spatial arrangement of cloud parcels, but they are computationally expensive.For instance, a standard MMF configuration is on the order of 180 times slower than conventional climate models (Khairoutdinov and Randall, 2001).
Finally, we note that CAM-CLUBB-SILHS deviates from common practice in microphysical parameterization.Namely, climate models typically use separate microphysics schemes for separate cloud types, such as stratiform and cumulus clouds.For instance, a relatively sophisticated microphysics scheme might be used in stratiform cloud, and a simpler microphysics scheme might be used in a massflux parameterization (e.g., Donner et al., 2011;Neale et al., 2012).In contrast, CAM-CLUBB-SILHS uses a single microphysics scheme, MG1, in all cloud types.Although we have previously mentioned some advantages of using a single, unified parameterization for clouds and turbulence, there are also advantages to using a single, unified scheme for microphysics.For instance, use of a single microphysics scheme avoids complexity and allows aerosol effects on clouds to be parameterized in all cloud types.

The subcolumn software framework in CAM
The subcolumn software framework in CAM is a newly developed piece of infrastructure that allows subcolumn samplers, such as SILHS, to feed subcolumn values from clouds to microphysics.The subcolumn framework will be available publicly in the release of CAM 5.4 and later versions, and is described in Appendix A.
The call sequence involving subcolumns is as follows: 1. CLUBB calculates a multivariate PDF that contains information about the subgrid variability of temperature, vapor, cloud liquid (mass) mixing ratio, cloud droplet number mixing ratio, cloud ice (mass) mixing ratio, cloud ice number mixing ratio, and vertical velocity.
2. The subcolumn software framework passes information about CLUBB's PDF to the SILHS sampler.
Each subcolumn includes all the aforementioned variates in CLUBB's PDF.The subcolumn framework creates a new model state data structure with these profiles.
4. Microphysics computes tendencies for all microphysical variates for each subcolumn, on the assumption that each subcolumn is horizontally uniform (e.g., overcast or cloud free).Aerosol tendencies are not computed on subcolumns.
5. The subcolumn tendencies are averaged together to obtain a grid-mean tendency.This averaging is done by the subcolumn framework using weights provided by SILHS.
6.The grid-mean tendency is applied to the grid-scale values in each column.Energy and water conservation checks are performed.
In order to visualize the flow of the calculations in CAM-CLUBB-SILHS, a schematic is provided in Fig. 1.
In order to ensure conservation of water and energy, the version of CAM-CLUBB-SILHS presented here modifies the sample values such that the weighted mean of all samples is constrained to be the same as the grid-mean value.In the limit of many sample points, the sample mean of the subcolumns converges to the grid mean seen by CLUBB.With a finite number of samples, however, the sample mean will in general differ from the grid mean.If, hypothetically, microphysics was evaluated on a set of samples whose mean exceeded the grid mean, the averaging could produce a mean drying tendency that is larger than the amount of water actually present in the grid column, even though the microphysics guarantees that each subcolumn individually returns non-negative values of water.If this excessive tendency were applied to the grid mean, the resulting negative water would be reset to zero by the energy checker, and a spurious source of water would be created.We prevent this from occurring by scaling the subcolumn values at each level and each time step by a constant factor, so that the weighted mean of the subcolumns exactly matches the grid-mean value at that point.The scaling occurs after SILHS has drawn sample values but before those values have been fed into the microphysics.This scaling has the undesirable side effect of effectively reducing the standard deviation of the subgrid PDFs.However, CLUBB's assumption that the standard deviation is proportional to the mean has uncertainty regardless of whether any scaling is done.Other than this scaling, no upper limit is placed on the values of the samples.We constrain the means of water vapor, liquid and ice mass mixing ratio, and liquid and ice number mixing ratio, but not temperature and vertical velocity.It is currently unknown how much the cost per subcolumn can be reduced by optimization.Another way to reduce the cost is to draw more representative subcolumns, so that fewer subcolumns are needed.In the future, we will evaluate a new sampling method that produces equal accuracy with about half as many subcolumns (Raut and Larson, 2015).
These test runs for timing do not attempt to vectorize subcolumn calculations.Since subcolumns do not communicate with each other, they can be efficiently parallelized.For this reason, subcolumn-based methodologies are well-suited to take advantage of vector processing and the next generation of high-performance computers.

Mean climate
This section evaluates the time-averaged climatology simulated by CAM-CLUBB-SILHS.We compare three versions of CAM -CAM-CLUBB-SILHS with 2 • horizontal resolution, CAM-CLUBB-SILHS with 1 • horizontal resolution, and CAM 5.3 with 2 • horizontal resolution -to a range of observational data sets that are summarized in Table 3.More information on each observational field, including specific references and discussion of observational uncertainties, can be found online with the NCAR Climate Data Guide at https://climatedataguide.ucar.edu/.In all figures in this section, the first row of plots shows the total field, and the second row shows differences from observations (model − obs).
Total surface precipitation rates for the three model versions and the Global Precipitation Climatology Project (GPCP) observations are presented in Fig. 2. CAM 5.3 exhibits a moderate, spurious double Inter-Tropical Convergence Zone (ITCZ), that is, a double band of precipitation in the Indian Ocean, and, to a lesser extent, in the equatorial Pacific.Both versions of CAM-CLUBB-SILHS produce a single band of rain through the tropics, thereby reducing the double-ITCZ bias.However, CAM-CLUBB-SILHS' precipitation is too intense and its ITCZ is too narrow, as compared to GPCP observations.The overall pattern of precipitation is similar between the 2 • and 1 • simulations, but the root mean  CAM-CLUBB-SILHS slightly improves the mean climatological column-integrated water vapor (Fig. 3), as compared to National Aeronautics and Space Administration (NASA) Water Vapor Project (NVAP) observations.CAM 5.3's overestimate of precipitable water is reduced in both the CAM-CLUBB-SILHS 2 • and 1 • simulations.The 2 • and 1 • simulations resemble each other, with the 1 • simulation providing a closer match to observations.Furthermore, the bias in precipitable water for the 1 • simulation is reduced by a factor of 4 as compared to the results of Guo et al. (2015).The improvement may be related to the fact that SILHS contains a detailed representation of hydrometeor/vapor overlap (Larson and Schanen, 2013;Storer et al., 2015), which influences the evaporation or accretional growth of precipitation as it falls to the ground or ocean (Jakob and Klein, 1999).
The top-of-the-atmosphere (TOA) long-wave cloud forcing (LWCF) for the three models is compared to observations in Fig. 4.Both versions of CAM-CLUBB-SILHS have smaller bias and lower RMSE in LWCF than does CAM 5.3.Furthermore, CAM-CLUBB-SILHS' bias is about a factor of 4 less than that of the simulation of Guo et al. (2015).The representation of LWCF in CAM-CLUBB-SILHS is aided by the fact that SILHS samples within-cloud variability of ice and feeds it into the microphysics scheme.Within-cloud subgrid-scale variability in ice is important because several ice processes are nonlinear (Morrison and Gettelman, 2008).The use of CLUBB-SILHS improves the TOA shortwave cloud forcing (SWCF) (Fig. 5).CAM 5.3 produces excessively reflective clouds over tropical land masses, probably because the deep convective microphysics does not precipitate out sufficient liquid cloud water.Use of CAM-CLUBB-SILHS, however, mitigates the excessive reflectivity of deep cumuli.The improvement may be related to the fact that accurate parameterization of the SWCF of deep cumuli requires accurate coupling of subgrid variability of clouds and precipitation (which in CAM-CLUBB-SILHS is handled by SILHS) and also requires accurate parameterization of deep convective microphysics itself (which in CAM-CLUBB-SILHS is handled by MG1).
The total grid-box cloud fraction for the three models and observations is presented in Fig. 6.Both versions of CAM-CLUBB-SILHS have a slightly lower cloud fraction (by about 5 %) than do CAM 5.3 or the observations.This is largely due to a lower cloud fraction throughout the tropics and subtropics in CAM-CLUBB-SILHS.
CAM-CLUBB-SILHS has about 35 % more total gridmean liquid water path (LWP) than does CAM 5.3, improving the agreement with observations (Fig. 7).It is notable that CAM-CLUBB-SILHS improves (increases) LWP without degrading (increasing the magnitude of) SWCF.How do the clouds in CAM-CLUBB-SILHS increase in water mass without increasing in reflectivity?A first reason is that CAM-CLUBB-SILHS' cloud fraction is slightly decreased in the tropics, as noted earlier.The decrease in cloud fraction, coupled with the increase in LWP, indicates that within-cloud cloud liquid water is increased in CAM-CLUBB-SILHS, ei- ther because the cloud liquid water has a more adiabatic profile, is more vertically stacked, or is more temporally intermittent.This piled-up vertical structure of LWP allows more solar radiation to reach the ocean or land surface (not shown) and thereby leads to reduced cloud reflectivity per unit of LWP.A second reason is that CAM-CLUBB-SILHS' cloud droplet effective radius is increased (not shown), thereby decreasing the reflectivity per unit of within-cloud LWP.Accurate simulation of droplet radius in deep convection requires accurate formulation of microphysics, which in CAM-CLUBB-SILHS is handled by the MG1 microphysics.
Figure 8 shows the Taylor score diagram for CAM 3.5, CAM 5.3, and 2 • CAM-CLUBB-SILHS (Taylor, 2001).CAM-CLUBB-SILHS is competitive with CAM 5.3 on most metrics, but has a higher RMSE in land rainfall, ocean rainfall, and the Pacific surface stress.The fact that Pacific sur-face stress is degraded suggests that CLUBB's formulation of vertical momentum flux, which is based on downgradient diffusion, needs to be modified in future work.
Table 4 shows the TOA global mean values for several radiation and energy balance terms, with mean values calculated from the observational data sets described in Table 3 and estimated uncertainties from Stephens et al. (2012).Unlike CAM 5.3, CAM-CLUBB-SILHS has not yet been tuned for top-of-model (TOM) radiative balance.Such tuning will be necessary before coupled simulations are attempted.
The differences between CAM-CLUBB-SILHS at 2 • or 1 • horizontal resolution are minor in both globally averaged radiation (Table 4) and in the spatial patterns of radiation and cloud fields (Figs. 3 to 8).This suggests that CAM-CLUBB-SILHS is relatively insensitive to small changes in horizontal resolution, aside from localized phe-  nomena such as near-coastal marine stratocumulus clouds.
In CAM-CLUBB-SILHS, the treatment of subgrid variability is removed from the microphysics and handled instead by CLUBB and SILHS.This removes one potential source of sensitivity to grid scale.

Madden-Julian Oscillation and tropical variability
The outgoing long-wave radiation (OLR) power divided by the background spectrum for various zonal wavenumbers and frequency (as in Wheeler and Kiladis, 1999, Fig. 3b) is shown in Fig. 9 for CAM 5.3, CAM-CLUBB-SILHS, and observations.This figure shows that, as compared to CAM 5.3, CAM-CLUBB-SILHS has increased power in the lowwavenumber, low-frequency, eastward-propagating region of the spectrum associated with the Madden-Julian Oscillation (MJO).The MJO power is not as strong as in the observations, and the frequency is slightly too high.The power associated with Kelvin waves is also increased in CAM-CLUBB-SILHS as compared to CAM 5.3, and compares well to the observations.However, CAM-CLUBB-SILHS has too much power in the high-frequency, westward-propagating side of the spectrum often associated with large convective systems advected westward by the mean flow (Wheeler and Kiladis, 1999).
Figure 10 shows the 20-80 day bandpass filtered precipitation and U 850 hPa winds at a given lag relative to a composite MJO passage and at a given longitude (top) and latitude (bottom).The MJO precipitation for CAM-CLUBB-SILHS is weaker than both the observations and CAM 5.3, but shows eastward propagation at the correct phase and speed.The overall coherence and structure of the MJO is much better in CAM-CLUBB-SILHS than in CAM 5.3.Figure 10 indicates that CAM 5.3 has primarily westward propagation of disturbances at this scale and has westerly winds nearly in phase with the maximum in precipitation.In contrast, CAM-CLUBB-SILHS simulates eastward propagation, with east- ward winds leading the precipitation and westward winds following, as seen in the observations.
In order to investigate differences in tropical convective processes between CAM 5.3 and CAM-CLUBB-SILHS, Fig. 11 shows average profiles of relative humidity, total physics moisture tendency and total physics temperature tendencies per value of rain rate for latitudes between 15 • N and 15 • S and longitudes between 60 and 180 • E (the Indian Ocean and west Pacific warm pool).These are similar to diagnostics used to evaluate MJO fidelity in Thayer-Calder and Randall (2009), Kim et al. (2009), Xavier (2012), andKim et al. (2014).All of these studies stress the importance of a smooth, gradual build-up in moisture from shallow convection (and light precipitation) to deep convection (and intense precipitation).
In observations, and in most models with a realistic MJO simulation, deep convection occurs in a nearly saturated column (Bretherton et al., 2004;Kim et al., 2009;Halloway and Neelin, 2009).Figure 11 shows that CAM 5.3 does not produce rain rates as intense as those simulated in CAM-CLUBB-SILHS; thus, the right-most profiles are missing.However, both models have nearly saturated profiles for the most intense rain rates that do occur.CAM-CLUBB-SILHS has a deeper boundary layer with higher relative humidity for mid-range precipitation values (between 0.5 and 10 mm day −1 ) than CAM 5.3.The relative humidity contours also show a smoother transition between light and intense precipitation than CAM 5.3.The transition from 80 % relative humidity in the boundary layer to near saturation around 11 mm day −1 in CAM 5.3 is more abrupt than reanalysis shown in similar results from Kim et al. (2009, Fig. 13) and Xavier (2012, Fig. 3).This abrupt transition may be an ill effect of poor deep convection triggering function.In contrast, the unified convection in CAM-CLUBB-SILHS produces a smooth deepening of the boundary layer into a fully saturated column at high rain rates.
Figure 11 also shows the total physics moisture and temperature tendencies for both models.CAM-CLUBB-SILHS shows strong moistening in shallow convective layers that transitions smoothly to intense drying through the entire column for deep convection.Similarly, the temperature tendencies smoothly change from low level heating, to convection rising in depth, to intense heating through nearly the entire column.These profiles resemble results for the SP-CAM presented in Thayer-Calder and Randall (2009, Figs. 4 and 9).The SP-CAM has been shown to simulate a realistic MJO (Khairoutdinov et al., 2008;Benedict and Randall, 2009), and so producing similar results in these diagnostics is promising.
In contrast, CAM 5.3 seems to have two main regimes.In the first, shallow convection produces light moistening tendencies above and below a layer of cloud-related drying around 900 hPa.This cloud layer produces a positive temperature tendency above a layer of cooling for all precipitation rates between 0.0003 and 2.5 mm day −1 .Past this point, there is an abrupt transition to convective drying and warming below 700 hPa, and then to a full column of drying above about 30 mm day −1 .However, unlike CAM-CLUBB-SILHS, the most intense precipitation in CAM 5.3 has strong heating only above 600 hPa.Again, the transition in moistening and heating rates is more abrupt than that seen in similar plots by Thayer-Calder and Randall (2009).There is a clear signal in CAM 5.3 of an unrealistic transition from convection handled by the shallow/stratiform parameterizations to convection produced by the Zhang and McFarlane (1995) deep convection parameterization.
There are still deficiencies in the simulation of the MJO by CAM-CLUBB-SILHS, but our unified parameterization of clouds produces promising improvements in the build-up of tropical moisture and the transition from shallow to deep convection.Boyle et al. (2015) show that an acceptable MJO in CAM 5.3 can be produced with tuning changes, but only at the expense of the mean climate.Our structural changes to CAM 5.3 have, in one and the same simulation, produced a realistic mean climate and improved tropical variability.

Subcolumn impact
In order to evaluate the impact of the number of subcolumns on these simulations, we performed four sensitivity experiments.All four simulations use the exact same settings and tuning parameters as the main 20-year, 2 • simulation described in Sect.assumptions about subgrid variability that are operative in CAM5.3, including a subgrid integration over cloud liquid water.The deep convection parameterization remains turned off here.This simulation indicates how CAM5.3 would behave if it used CLUBB as a unified parameterization and it used MG1's subgrid assumptions, developed for stratiform clouds.The three other simulations varied the number of subcolumns from 4 to 10 to 50.Because of restrictions in the SILHS importance sampling algorithm (Larson and Schanen, 2013), the number of subcolumns must always be divisible by two.
As expected, the simulation without subcolumns produces an unrealistic climate.Figure 12 shows that the No Subcolumns simulation has very low long-wave cloud forcing, and Table 4 shows this simulation has the highest OLR, largest radiative imbalance, and greatest error in SWCF.This is likely because the convection is not penetrating as deeply into the atmosphere, and the clouds are not cold and icy enough.This is supported by the large SWCF for the simulation (Table 4), which has a high bias and RMSE in Fig. 13. Figure 14 shows that this simulation has an even lower LWP than that of CAM 5.3.
The simulation with only four subcolumns shows marked improvement over the No Subcolumns simulation.Table 4 shows a large decrease in both net solar TOA flux and OLR, with reasonable values of LWCF, but a lower SWCF corresponding to brighter clouds.This is also seen in Fig. 13, where the low bias in SWCF is distributed over all oceans.This low bias in SWCF is tied to the higher cloud LWP for this simulation (Fig. 14).The 4-subcolumn simulation appears to have a lower precipitation efficiency than the 10subcolumn simulation.The reason, we speculate, is that use of a limited number of subcolumns leads to poor sampling of the tails of the distribution, which is where precipitation forms and grows.
The 10-and 50-subcolumn simulations are similar, suggesting that climatological averages are fairly close to converged even when only 10 subcolumns are used.Table 4 shows that increasing to 50 subcolumns decreases the OLR by 1 W m −2 and increases the net solar flux by 0.5 W m −2 .Figure 12 shows that the LWCF is very similar between the 10-and 50-subcolumn runs.Both simulations have similar SWCF (Fig. 13) and LWP (Fig. 14).The fact that LWP decreases when the number of subcolumns is increased to 50 supports the hypothesis that increasing subcolumns increases precipitation efficiency, although there is diminishing effect after 10 subcolumns.

Summary and conclusions
This paper evaluates a version of CAM, "CAM-CLUBB-SILHS", that uses a single equation set to parameterize all cloud types, including shallow convective, deep convective, and stratiform liquid and ice clouds.The equation set is CLUBB's set of equations for higher-order moments.CLUBB uses the higher-order moments to construct a multivariate subgrid PDF, which, in turn, is sampled by SILHS.The samples are then used to drive a single microphysics scheme, MG1, that acts on all cloud types.In CAM-CLUBB-SILHS, clouds are parameterized in a more fully unified way, and so are microphysical processes.
The use of a single, multivariate subgrid PDF fosters consistency in the sense that all cloud and microphysical processes see the same subgrid PDF.In this paper, the PDF has been extended to include cloud ice mass and number, thereby incorporating subgrid variability in ice processes.
As compared to CAM5, the most important degradation in the CAM-CLUBB-SILHS simulations is the RMSE in surface precipitation rate.In particular, the surface precipitation field is stronger in the precipitating regions than that observed by satellite.However, several aspects of the simulations have been improved.We list the improvements here, even though it is difficult to pinpoint their causes.
First, CLUBB-SILHS slightly reduces CAM5's overestimate of precipitable water.This may be related to the fact that CLUBB-SILHS contains a detailed representation of vertical overlap, which affects the relative rates of evaporation and accretional growth of precipitation.
Second, CLUBB-SILHS improves LWCF.In general, CLUBB-SILHS offers a more detailed representation of subgrid variability in ice because cloud ice mass and number mixing ratio are included in the subgrid PDF.The inclusion of ice in the PDF, in turn, allows for subgrid variability in ice to drive ice-related microphysical processes.
Third, CLUBB-SILHS simultaneously improves the simulation of both LWP and SWCF.In CAM5, LWP is underestimated by almost a factor of 2, and deep convective clouds are too reflective over the tropical continents.In CAM-CLUBB-SILHS, LWP is increased without unduly increasing the magnitude of SWCF.In part, this is related to the fact that CAM-CLUBB-SILHS predicts smaller cloud fraction.That is, CAM-CLUBB-SILHS' liquid water content is more vertically and/or temporally correlated and less horizontally extended, allowing more LWP to be present without causing excessive cloud albedo.In addition, in CAM-CLUBB-SILHS, the cloud liquid droplet radius is increased, thereby reducing reflectivity of clouds.Fourth, although CLUBB-SILHS underestimates MJO wave activity, it improves (strengthens) the spectral power associated with the MJO and convectively coupled Kelvin waves.The improvement may be related to the fact that CLUBB-SILHS is a unified parameterization in which there is no categorization of clouds nor a cumulus trigger function.This allows for a smoother, more realistic transition between shallow and deep convection in the tropics.
The simultaneous improvement of LWP, SWCF, and tropical power spectrum is significant.Use of automated parameter estimation reveals that although CAM5's MJO can be improved by changes in parameter values, the improvement comes at the expense of the simulated climatology, including the absorption of shortwave radiation (Boyle et al., 2015).This suggests that, in order to simultaneously improve CAM 5.3's MJO and mean state, structural modifications to the parameterization suite are required.The use of CLUBB-SILHS is one possible structural modification.
The results are relatively insensitive to an increase in resolution from 2 to 1 • .Avoiding undesirable grid-scale sensitivity is aided by the fact that CAM-CLUBB-SILHS does not require the microphysics scheme to internally account for resolution changes.Instead, any model awareness of horizontal resolution is contained in CLUBB and is communicated to the microphysics via SILHS.As cloud-resolving resolutions are approached, CLUBB is designed to gradually shut itself off by reducing its turbulent dissipation timescale (Larson et al., 2012).Whether in practice the output of CAM-CLUBB-SILHS proves to be sensitive to significant changes in resolution is left for future work.
Although acceptable results can be found with as few as four sample points per grid box and physics time step, the results are moderately sensitive to the number of sample points.This suggests that climate simulations are sensitive to the details of subgrid variability within clouds and how such variability is communicated to the microphysics.Therefore, it is worth investigating subgrid integration methods, whether they be Monte Carlo methods or alternative methods.One alternative method is analytic integration, which is computationally inexpensive but is restricted to simple microphysical formulations (Morrison and Gettelman, 2008;Larson and Griffin, 2013;Griffin and Larson, 2013).Another alternative method is deterministic quadrature, which requires somewhat intrusive software changes but is more generally applicable than analytic integration (Golaz et al., 2011;Chowdhary et al., 2015).
Each subcolumn that is added increases the total model computational cost by about 6 %.This cost is reasonable, considering the wealth of detail that is output by subcolumns.
Much further unification of parameterizations of subgrid variability is possible in the future.Although CLUBB-SILHS unifies the parameterization of subgrid-scale variability in clouds and feeds that information into a microphysics scheme, that information is not fed consistently into aerosol, radiative, or land surface processes.That extension is left for future work.

Appendix A: CAM subcolumn implementation A1 Description of subcolumn implementation
Subcolumns were implemented in CAM to assist in the study of subgrid-scale physics.The implementation supports both studies based on spatial subdivision of a physics column and studies based on statistical sampling of subgrid variability (e.g., SILHS).Other features of the CAM implementation of subcolumns are -Subcolumns are used to study subgrid-scale physics in a select subset of physics parameterizations.
-Subcolumn data may be shared between parameterizations (for example, passing microphysics subcolumns to the radiative transfer scheme).
-The subcolumn scheme (see below) may specify a different number of subcolumns per grid column (e.g., 15 subcolumns per grid column in the tropics, 2 elsewhere).
-The memory layout provides for efficient, threaded performance and seamless use in current, portable code layers (see Fig. A1).
-Subcolumn data may persist across time steps.
-If subcolumns are not invoked, the basic model state is not altered.
-Parameterizations themselves do not need to know about subcolumns because information is passed at the interface and driver levels.
-Subcolumn information can be output for analysis.
Subcolumns in CAM are considered static: once the number of subcolumns in any grid column is set at the beginning of simulation, this number should not be changed.The subcolumn framework supports only instantaneous history output of subcolumn fields.Currently, the only CAM physics parameterization that accepts subcolumn input is the Morrison-Gettelman microphysics (Morrison and Gettelman, 2008).However, the software framework allows subcolumns to be applied to other parameterizations.A key goal is to apply subcolumns uniformly across the column physics: for example, currently there are separate subcolumn generators for radiation and satellite simulators in CAM, these could be made consistent with this framework.
Use of subcolumns begins with sampling or generation of subgrid fields based on the current physics state.In this way, a complete state on subcolumns is passed to the parameterization.The sampling can occur by any method (in this case SILHS) and for arbitrary fields.Parameterizations then use these fields to produce subgrid tendencies.Finally, the subgrid state and tendencies are averaged back to the grid scale.The number of subcolumns varies by grid column, as shown in the conceptual layout (left).Internally, the subcolumns are stored in a compressed layout (right).The information on per-chunk linkages is stored in a series of software parameters, some of which can be set by the user (black) and others of which are internally calculated (blue).pcols is the maximum number of grid columns, psubcols is the maximum number of subcolumns, psetcols is the maximum number of all columns (i.e., pcols • psubcols), ngrdcol is the actual number of grid columns that contain data, nsubcol(pcols) is the number of subcolumns in each grid column with data, ncol is the total number of all subcolumns with data, and indcol(psetcols) is the grid index for each subcolumn that is used for mapping subcolumns back to grid columns.
The subcolumn "gather" or averaging routines can be customized so that averaging can be performed using weights or masking if desired.Organization of different methods for drawing or generating subcolumns and averaging them back to the grid is described below.For more details or for documentation on making a parameterization subcolumn aware, see the CAM reference manual (Eaton et al., 2015).

A2 Implementing a new subcolumn scheme within CAM
Different methods or "schemes" for generating and averaging subcolumn fields can be invoked.SILHS is one subcolumn scheme.The use of a specific subcolumn scheme is controlled by the CAM subcol_scheme namelist variable.Each of the generic subcolumn interfaces listed below call the scheme-specific version based on the value of subcol_scheme.Scheme-specific versions of the following routines will need to be supplied, even if they contain no executable code.Typically the schemespecific versions are designated by the generic name followed by "_schemeName" , noted below by XXX (e.g., subcol_register_SILHS).The routines are -subcol_register_XXX: register any subcolumn-specific physics buffer fields using pbuf_add_field; -subcol_readnl_XXX: read any subcolumnscheme-specific namelist parameters; -subcol_init_XXX: perform subcolumn-specific initialization, set up any output calls for subcolumn diagnostics (via addfld), and initialize any subcolumn physics buffer fields, if required; -subcol_gen_XXX: contains the details of mapping state, physics tendencies, and physics buffer fields from the grid to subcolumns.Typically, this routine will be the interface between CAM and the unique code for generating the subcolumns or drawing them from PDFs.
Once physics tendencies and/or updates are computed for each subcolumn, the subcolumn values need to be averaged back onto the CAM grid.This is accomplished via calls to averaging routines.The default behavior of these routines is to perform a simple average, applying optionally supplied scheme-specific weights and/or filters, such as a cloud mask or conditional sampler.If a more sophisticated method is required, scheme-specific routines may be supplied for these two routines.
-subcol_ptend_avg_XXX: average the subcolumn physics tendency values back to the grid so that these values can be applied to the grid-resolved state.
-subcol_field_avg_XXX: average the physics buffer fields from subcolumn values back to the grid.This function only needs to be called for physics buffer fields that are used in other parameterizations on the grid.
The data layout for subcolumns is illustrated in Fig. A1.The number of subcolumns varies by grid column, as shown in the conceptual layout (Fig. A1, left).Internally, the subcolumns are stored in a compressed layout (Fig. A1 right).The information on organization is stored in a series of parameters some of which can be set by the user (black) and others which are internally calculated (blue).For further details, see Eaton et al. (2015).

Figure 1 .
Figure 1.The sequence of calculations in CAM-CLUBB-SILHS.Red lines represent temperature profiles and dark blue lines represent moisture profiles, as an example.Light blue lines represent figurative microphysical tendencies, for both temperature and moisture.For details on the SILHS and subcolumn methodology, see Sect. 2.

Figure 9 .
Figure 9. OLR power divided by the background spectra for various wavenumbers and frequencies in the tropics for CAM 5.3 (left), CAM-CLUBB-SILHS (center), and NOAA OLR observations (right).CAM-CLUBB-SILHS has stronger MJO signal, and a stronger signal for Kelvin waves, than CAM 5.3

Figure 11 .
Figure 11.Daily average profiles of fields per daily average value of precipitation for the region between latitudes 15 • N and 15 • S and longitudes 60 to 180 • E over 1 year.Relative humidity for CAM-CLUBB-SILHS (top left) and CAM5 (top right), total physics moisture tendencies for CAM-CLUBB-SILHS (middle left) and CAM5 (middle right), and total physics temperature tendencies for CAM-CLUBB-SILHS (bottom left) and CAM5 (bottom right) are shown.Because CAM-CLUBB-SILHS is a unified parameterization, there is a smoother transition from light to intense precipitation values for all fields.

Figure 12 .Figure 13 .
Figure 12.LWCF difference from CERES-EBAF observations for 5 years of simulation without subcolumns at all (top left), 4 subcolumns (top right), 10 subcolumns (bottom left), and 50 subcolumns (bottom right).Without subcolumns, the LWCF has a severe low bias.The simulation with 4 subcolumns has the lowest global error, and very little changes between the simulations with 10 and 50 subcolumns.

Figure 14 .
Figure 14.LWP for 5 years of simulation without subcolumns at all (top left), 4 subcolumns (top right), 10 subcolumns (bottom left), and 50 subcolumns (bottom right).Without subcolumns, the model has little cloud liquid water.The simulation with 4 subcolumns has a large amount of cloud liquid, which moderates in the 10 and 50 subcolumn simulations.

Figure A1 .
Figure A1.Schematic for CAM's subcolumn software framework.The number of subcolumns varies by grid column, as shown in the conceptual layout (left).Internally, the subcolumns are stored in a compressed layout (right).The information on per-chunk linkages is stored in a series of software parameters, some of which can be set by the user (black) and others of which are internally calculated (blue).pcols is the maximum number of grid columns, psubcols is the maximum number of subcolumns, psetcols is the maximum number of all columns (i.e., pcols • psubcols), ngrdcol is the actual number of grid columns that contain data, nsubcol(pcols) is the number of subcolumns in each grid column with data, ncol is the total number of all subcolumns with data, and indcol(psetcols) is the grid index for each subcolumn that is used for mapping subcolumns back to grid columns.

Table 3 .
Information about observational data sets used for comparison in this paper.More information about each of these can be found on the website for the National Center for Atmospheric Research (NCAR) Climate Data Guide at https://climatedataguide.ucar.edu/.

Table 4 .
Globally averaged top-of-the-atmosphere (TOA) radiation fields, and the top-of-model (TOM) radiation imbalance for various configurations of CAM-CLUBB-SILHS.Estimates of observational uncertainty are from Stephens et al. (2012).Values are in units of W m −2 .