Application and evaluation of a new radiation code under McICA scheme in BCC_AGCM2.0.1

This research incorporates the correlated k distribution BCC-RAD radiation model into the climate model BCC_AGCM2.0.1 and examines the change in climate simulation by implementation of the new radiation algorithm. It is shown that both clear-sky radiation fluxes and cloud radiative forcings (CRFs) are improved. The modeled atmospheric temperature and specific humidity are also improved due to changes in radiative heating rates, which most likely stem from the revised treatment of gaseous absorption. Subgrid cloud variability, including vertical overlap of fractional clouds and horizontal inhomogeneity in cloud condensate, is addressed by using the Monte Carlo Independent Column Approximation (McICA) method. In McICA, a cloud-type-dependent function for cloud fraction decorrelation length, which gives zonal mean results very close to the observations of CloudSat/CALIPSO, is developed. Compared to utilizing a globally constant decorrelation length, the maximum changes in seasonal CRFs by the new scheme can be as large as 10 and 20 W m −2 for longwave (LW) and shortwave (SW) CRFs, respectively, mostly located in the tropics. The inclusion of an observation-based horizontal inhomogeneity of cloud condensate has also a significant impact on CRFs, with global means of ∼ 1.5 W m−2 and∼ 3.7 Wm−2 for LW and SW CRFs at the top of atmosphere (TOA), respectively. Generally, incorporating McICA and horizontal inhomogeneity of cloud condensate in the BCC-RAD model reduces global mean TOA and surface SW and LW flux biases in BCC_AGCM2.0.1. These results demonstrate the feasibility of the new model configuration to be used in BCC_AGCM2.0.1 for climate simulations, and also indicate that more detailed real-world information on cloud structures should be obtained to constrain cloud settings in McICA in the future.


Introduction
Radiation process is crucial for climate simulations.Over the past 2 decades, a lot of progress has been made in atmospheric radiation.For example, in most radiation schemes the traditional band model for gaseous transmittance has been replaced by the correlated k distribution (CKD) method (Fu and Liou, 1992;Mlawer et al., 1997;Zhang et al., 2003;Li and Barker, 2005;Shi et al., 2009, and others).In principle, CKD can be applied to a single absorption line while band models utilize mean values for entire bands.Thus, CKD shows considerable promise in simulating atmospheric radiation accurately and efficiently.The breakthrough in nonspherical scattering makes it possible to accurately calculate ice cloud optical properties (Mishchenko et al., 2002).In an ice cloud optical property parameterization, the full set of single-scattering properties is provided by considering three-dimensional random orientations for multiple ice crystal habits following the observations of CALIPSO (Yang et al., 2005(Yang et al., , 2013;;Baum et al., 2011).All this progress has been included in the Beijing Climate Center Radiation transfer model (BCC-RAD), which is used in the general circulation model (GCM) of Beijing Climate Center (i.e., BCC_AGCM2.0.1).
Published by Copernicus Publications on behalf of the European Geosciences Union.

H. Zhang et al.: Application and evaluation of a new radiation code
How to deal with cloud vertical overlap and cloud internal inhomogeneity has been a very difficult task for atmospheric radiation.This arises mostly from the relatively coarse spatial resolution of GCMs (dozens to hundreds of kilometers), which leaves cloud-relevant processes and inherent subgrid variations of clouds unresolved (Barker and Räisänen, 2005;Zhang et al., 2013).Typically, cloud condensate (water and ice) is treated as horizontally homogeneous (the plane parallel homogeneous, or PPH, assumption) within a GCM grid cell.Additionally, certain predetermined assumptions about the vertical overlap of fractional clouds are required (traditionally maximum-random overlap, or MRO) (Tian and Curry, 1989).Computing on a cloud system resolving model (CSRM) data set, Barker and Räisänen (2005) found that, by a small change to the standard deviation of cloud condensate distribution, zonal mean shortwave (SW) cloud radiative forcing (CRF) could change up to 25 W m −2 at certain latitudes (with a global mean of ∼ 8 W m −2 ).The radiative sensitivity to cloud overlap is of similar magnitude for global averages.Oreopoulos et al. (2012) included a beta distribution function and a latitude-/day-dependent cloud overlap function, both derived from CloudSat/CALIPSO data, in the GEOS-5 atmospheric general circulation model.SW and longwave (LW) CRFs showed significant changes compared to the traditional PPH and maximum-random overlap setup, but the magnitude of changes depends also on the cloud scheme utilized.All of these studies have emphasized the importance of faithfully addressing subgrid cloud variability in GCMs.
To make the representation of subgrid cloud properties flexible and modularized and to maintain computational efficiency, a scheme named the Monte Carlo Independent Column Approximation (McICA) method was developed (Pincus et al., 2003).McICA is a method that does fast spectral integration over given cloudy subcolumns within a domain.The cloud subcolumns required by McICA can be supplied from certain subgrid cloud generators (Räisänen and Barker, 2004) or cloud-resolving models (Hill et al., 2011).The advantages of McICA are that it facilitates adjustment or alteration of both cloud structure and radiative transfer and thus accelerates future development of GCMs.
Though McICA has been extensively studied, there lacks a detailed description on cloud-type-related vertical overlap.An e-folding relationship of cloud overlap has been developed to quantitatively represent cloud overlap for different types of clouds and over different regions (Hogan and Illingworth, 2000;Mace and Benson-Troth, 2002).In current GCMs, a global mean constant value of decorrelation length (hereinafter L cf ), a critical parameter in this overlap algorithm, is often used.Usually, convective cloud has well-organized structures with large vertical correlation due to strong upward motion.Therefore, L cf for convective cloud should be larger than for the other types of clouds.To address the climate impact by distinguishing the cloud-type-related L cf in McICA is one goal of this work.We here incorporate the McICA scheme and a stochastic cloud generator (SCG; Räisänen et al., 2004) into the BCC_AGCM2.0.1 model with the BCC-RAD radiation algorithm.
This work contains two parts: first, we report the improvements on climate simulation by introducing the new BCC-RAD radiation algorithm.Second, we analyze the impact of cloud overlap assumption and cloud condensate inhomogeneity through the McICA scheme.Most attention will be paid to the cloud-type-related decorrelation length by comparing with the results of using a globally constant value.This preliminary work aims to document the impact of the modifications in cloud-radiation process on simulated climate and the model response to these changes and thereby provide suggestions for future development.In Sect.2, the BCC_AGCM2.0.1 model, the BCC-RAD radiation scheme and the McICA scheme are briefly described.The design of experiments is given in Sect.3. Results of the simulations with various model configurations are described in Sect. 4. In Sect.5, we conclude with a brief summary.(Wu et al., 2010).The model runs at T42 spectral resolution (approximately 2.8 • × 2.8 • ) horizontally, and it uses vertical hybrid δ-pressure coordinates including 26 layers with the top located at about 2.9 hPa.An additional layer is added above the topmost layer in the radiative calculation to prevent excessive heating.The default timestep is 20 min, and the radiation code is invoked every three timesteps.
Relative to CAM3, several revisions have been made to improve the physics of the model.These include new reference atmosphere and surface pressures; a revised convection scheme (Zhang and Mu, 2005) that significantly improves the tropical rainfall simulation; a different function for calculating the snow-cover fraction that influences the resulting surface albedo, especially in polar and plateau regions (Wu and Wu, 2004); a new adiabatic adjustment originated by Yan (1987); and new methods for calculating turbulent fluxes over ocean surface that remove the systematic biases in the wind stress and sensible and latent heat fluxes in CAM3.A more detailed description of BCC_AGCM2.0.1 can be found in Wu et al. (2010).In the present research, the interactive Canadian Aerosol Module (CAM) (Gong et al., 2003) with updated aerosol emission sources (Zhou et al., 2012) is used to predict atmospheric aerosol burdens (Zhang et al., 2012).Kiehl and Briegleb, 1991) (SW: Briegleb, 1992) CKD scheme (Zhang et al., 2003(Zhang et al., , 2006a, b) , b) RT solver in LW Absorptivity/emissivity formulations (Ramanathan and Downey, 1986) Two-stream approximation (Nakajima et al., 2000) RT solver in SW δ-Eddington method (Briegleb, 1992) δ-Eddington method (Coakley et al., 1983) Cloud fraction parameterization Diagnostic scheme (Rasch and Kristjansson, 1998) The same as in Old Cloud optics LW: emissivity formulations (Ebert and Curry, 1992); SW: formulas of Slingo (1989) for liquid and of Ebert and Curry (1992) for ice Ice cloud: computed using data from Fu (1996), Yang et al. (2005), and  (Collins, 2001) McICA (Räisänen and Barker, 2004;Barker et al., 2008) Aerosol-radiation coupling scheme BCC_AGCM2.0.1_CAM (Zhang et al., 2012) BCC_AGCM2.0.1_CAM (Zhang et al., 2012) * In the new scheme, contributions from the solar spectrum and terrestrial emission are mixed within 2110-2680 cm −1 .
The BCC-RAD model is substantially different from the previous radiation scheme used in BCC_AGCM2.0.1.To explain the importance of this radiation scheme in modulating climate simulation, it is necessary to describe this revision in advance.A detailed comparison between the old and new schemes is provided in Table 1.
The previous radiation scheme in BCC_AGCM2.0.1 is basically a band model.Although some band models simulated well the broadband fluxes and heating rates, this may have been partly fortuitous because of band overlap effects (Ellington et al., 1991).Another defect of band models is the use of a scaling procedure to account for inhomogeneous atmospheric paths, although these can be made arbitrarily accurate for a homogeneous atmosphere (Kratz, 1995).Therefore, there has been a trend over the past decades to replace band models with CKD methods in GCMs.
The 10-49 000 cm −1 (0.204-1000 µm) spectral range in BCC-RAD is divided into 17 bands (8 LW and 9 SW).Five major greenhouse gases (GHGs) -H 2 O, CO 2 , O 3 , N 2 O, and CH 4 -as well as chlorofluorocarbons (CFCs) are considered.The major absorbers in the solar bands are H 2 O (including continuum absorption), CO 2 , N 2 O, O 3 , and O 2 .The HITRAN2000 database (Rothman et al., 2003) was used to provide line parameters and cross sections.Lu et al. (2012) compared the line parameters in different HITRAN versions and found that the difference in the simulated radiative fluxes between the updated HITRAN2008 and HITRAN2000 is very small, so the use of HITRAN2000 should not affect the final modeled climates in this research.In BCC-RAD, the effective absorption coefficients of CKD are calculated based on the line-by-line radiative transfer model (LBLRTM; Clough and Iacono, 1995) with a spectral interval of onefourth the mean half-width and a 25 cm −1 cutoff for line wings over each band (Clough and Iacono, 1995).The thermal radiation transfer calculation is solved with a two-stream algorithm developed by Nakajima et al. (2000), and the solar radiation transfer is solved with the δ-Eddington method (Coakley et al., 1983).SW radiation model comparisons, including BCC-RAD, are given in Randles et al. (2013).
Cloud and aerosol optical properties in BCC-RAD are also different from those in the original scheme.The optical properties of cloud droplets are from Nakajima et al. (2000), and those of ice crystals are calculated based on several data sets: observational size distribution data from Fu (1996), optical properties of single particles of different shapes from Yang et al. (2005), and the fractional mixing of particles of various shapes suggested by Baum et al. (2005).Aerosol optical properties are from Wei and Zhang (2011) and Zhang et al. (2012).

Description of the McICA scheme
The McICA scheme is based on the ICA algorithm for the computation of domain mean radiation fields.It greatly and effectively reduces computation time while maintaining the accuracy of ICA from a statistical perspective.The basic principles of McICA were first explained in detail by Pincus et al. (2003); Räisänen and Barker (2004) then provided additional ways to diminish the induced noise.For clarity and completeness, we provide a brief summary here.
Conceive a domain R (a GCM grid).The subgrid clouds could be represented by a certain number of subcolumns, which contain individual cells in each layer that are either clear or overcast.Moreover, the domain mean of these subcolumns should hold the cloud profile provided by the GCM.Given these subcolumns, radiative computation can be liberated from the description of partial clouds and their vertical overlap.The required subcolumns could be derived through SCG with consideration of certain overlap and horizontal distribution rules for clouds.For a thorough methodology of SCG, one can refer to Räisänen et al. (2004).
Within the domain R composed of subcolumns, the domain-averaged radiative fluxes can be accurately given by ICA as where x and y are subcolumn counters along the zonal and meridional axis, respectively, S (λ) is the spectral weight at wavelength λ, and F (x, y, λ) denotes the radiative flux at location (x, y) and wavelength λ.
If R is partially cloudy, F ICA can be split into clear F clr and cloudy F cld parts weighted by the cloud fraction Ac: The most time-consuming part of Eq. ( 2) is F cld due to the full spectral integration in all cloudy subcolumns.To diminish the computational burden, Pincus et al. (2003) reduced the two-dimensional integration to a single dimension by introducing a Monte Carlo (random sampling) process: where s rnd is a randomly selected cloudy subcolumn number for radiative calculation at λ. Equation (3) tremendously reduces computation time compared with Eq. ( 2) and represents the kernel of McICA.
It should be noted that the random selection of s rnd in Eq. (3) inevitably introduces random noise.Although this may yield deviated results for a single calculation, averaging over a number of calculations generates almost unbiased results with respect to ICA (Barker et al., 2008).One method for reducing the noise is to increase the number of s rnd for optically critical spectral intervals (Räisänen and Barker, 2004).To date, the McICA scheme has already been operationally utilized in several climate models and numerical weather prediction models (Morcrette et al., 2008;Räisänen and Järvinen, 2010;Neale et al., 2010).

Experimental design
We now have considered two model configurations: the new one with McICA and BCC-RAD to handle the cloudradiative procedure and the old one with the traditional overlap treatment by Collins et al. (2001) and radiation scheme described in Briegleb (1992).The details of these are listed in Table 1.Experiments were designed to reveal (a) the differences in simulated climate between the two configurations and (b) the impact of changing subgrid cloud structures on simulated climate within the new configuration.For all the following experiments, the sea surface temperature (SST) data are from the global Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) data set (Rayner et al., 2003) for years up to 1981 and the Reynolds et al. (2002) for years after 1981; the greenhouse gases are set the same by using the current values; aerosols are produced by a coupled aerosol model (named CAM) that is described in detail in Zhang et al. (2012).All of the experiments are integrated from September 1979 to December 1990, and the results of the last 10 years are used for analysis.

Experiments comparing the new and old model configurations
First, an experiment with the old scheme, denoted OLD, was performed as a control run.Second, a McICA experiment denoted N_MRO, utilizing PPH and the MRO assumptions to be consistent with the OLD run, was done.The comparison between N_MRO and OLD illustrates differences in the climate response due to changes in radiation scheme other than subgrid cloud variation.

Experiments exploring the impacts of subgrid cloud structures
As the McICA scheme is flexible in depicting subgrid cloud structures, three more experiments were implemented to test the model's sensitivity to cloud-structure variations.First, the impact of changing cloud overlap assumption was tested by including the so-called general overlap (hereafter GenO) (Mace and Benson-Troth, 2002).In GenO, the vertically projected cloud fraction of the two cloud layers k and l (C k,l ) is defined as the linear combination of maximum (C max k,l ) and random overlap (C ran k,l ): where and the overlap parameter α k,l is prescribed via an exponential decay function of altitude separation between cloud layers: The Eq. ( 4) is applied to both continuous and discontinuous clouds as in Mace and Benson-Troth (2002).The lapse rate of the decay of α k,l is controlled by the "decorrelation length" (L cf in Eq. 7), which has a global mean value of about 2 km (Barker, 2008).We term the scheme of using a globally constant L cf = 2 km as N_GO2.In reality, L cf is highly related to cloud type and atmospheric dynamics (Naud et al., 2008).Usually, convective cloud has well-organized structure (convective tower) with large vertical correlation due to strong upward motion.Therefore, L cf for convective cloud should be larger than for other types of clouds.Generally, L cf is about 5 km to 10 km for convective cloud and is much smaller (around 1 km) for other types of clouds (H.Barker, personal communication, 2013).In a GCM grid cell, we simply define the grid-scale mean result as where  Barker (2008).We term the scheme of Eq. ( 8) as N_GOF.This can also be called as convective-stratiform contrast scheme, since the stratiform cloud dominates all the other types of cloud in cloud fraction.
Additionally, the impact of breaking the default PPH assumption is addressed by perturbing the horizontal distribution of cloud condensate (water and ice) with an ideal distribution function.The gamma function of cloud condensate applied by Shonk et al. (2010) is used here.In such distribution, the magnitude of inhomogeneity is constrained by the fractional standard deviation (f ), which is defined as where c is the layer mean cloud condensate ignoring the cloud phase, and σ c is the standard deviation of the condensate.In this work, f was set to be 0.75 for both the liquid and ice phases, as was obtained by Shonk et al. (2010) from an extensive collection of observations.This inhomogeneity setup was tested with L cf given as Eq. ( 8), denoted as N_GOF_IH.Because the cloud overlap assumptions are consistent for N_GOF and N_GOF_IH, any discrepancies illustrate the impact of including horizontally inhomogeneous clouds.
In this study, the decorrelation length for cloud condensate in GenO is set to be 1 km in all tests.The variation of the decorrelation length for cloud condensate has smaller influence and is less important than that for the overlap of cloud fractions (Barker and Räisänen, 2005).

Results
This section reports the results in two parts: (i) first, results from OLD and N_MRO are provided to clarify the differences between the new and old model configurations; (ii) second, results from N_GO2, and N_GOF and N_GOF_IH are presented to show the impacts of cloud overlap variations and changing the horizontal distribution of cloud condensate within the McICA scheme.

Radiation budget
We first investigate the difference between the old and new schemes under the same subgrid cloud structure setups.Figure 1 shows the global annual mean radiation fields for various simulations at the top of atmosphere (TOA) and at the surface (SFC) with a comparison against the satellite-derived 11-year (2000-2010) mean CERES_EBAF data sets (http: //ceres.larc.nasa.gov/order_data.php)(Loeb et al., 2009).We focus on the results of OLD and N_MRO in this section.
The central column of Fig. 1 shows that the new scheme obtains much improved net all-sky LW and SW TOA radiative fluxes.This is due to improvements in both the revised cloud optics and the net clear-sky fluxes calculated by the new radiation scheme.
Compared with CERES_EBAF data, the OLD run shows notable discrepancies in TOA LW and SW CRFs (right column of Fig. 1), which are overestimated by ∼ 3 W m −2 and ∼ 7 W m −2 , respectively.The N_MRO run shows large reductions in these biases, with TOA LW and SW CRFs errors reduced to ∼ 1.5 and ∼ 3 W m −2 , respectively.As the same cloud overlap assumptions are used, the improved CRFs in N_MRO should come mainly from the revised cloud optics (see Table 1).
As for the clear-sky net fluxes at TOA (F clr , the left column in Fig. 1 by ∼ 5 and ∼ 1.5 W m −2 , respectively.The biases at the surface are also large, up to ∼ 4 W m −2 for SW F clr .Again, N_MRO produces clear-sky fluxes much closer to the observations (except for LW F clr at surface).The differences between the simulated TOA LW F clr and SW F clr and those from CERES_EBAF observations are reduced to ∼ 1 and ∼ 0.5 W m −2 , respectively, for N_MRO.The improvements in both F clr and CRFs suggest that the implementation of the new radiation scheme fares much better at modeling the inner balance between the radiation components from clear and cloudy regimes.Thus, the new configuration behaves in a more physically coherent manner than the original one in BCC_AGCM2.0.1, and it predictably yields more reasonable all-sky net fluxes (F net , central column of Fig. 1).
Figure 2 displays zonal annual mean F clr , F net , and CRFs at the TOA from the OLD, and N_MRO runs, as well as the CERES_EBAF data set.The simulated zonal distributions of these variables are all in reasonable agreement with observations.However, the N_MRO run gives LW and SW F net much closer to observations, especially at mid-low latitudes (Fig. 2c, d).This occurs mainly because the vast overestimation of LW and SW CRFs by the OLD scheme is reduced overall by the new scheme (Fig. 2e, f).Moreover, N_MRO also shows notable improvement in LW F clr in the subtropics and mid-latitudes (Fig. 2a).The SW F clr is calculated well at most latitudes in all experiments, except at the polar regions where there are noticeable underestimations.This may be linked to the enhanced solar albedo over snow surfaces compared with observations in the Community Land Model version 3 (CLM3) used in the BCC_AGCM2.0.1 model (Oleson et al., 2003), which results in an overestimated solar energy loss to space.
Figure 3 shows the global distribution of errors in annual mean SW and LW CRFs relative to CERES_EBAF, as well as the differences between N_MRO and OLD.The OLD run exhibits negative biases in SW CRF at most low and midlatitudes (see Fig. 3a).The N_MRO run significantly reduces these errors (see Fig. 3c), but the enhanced positive biases appear over subtropical oceans near the west coasts of continents and over East Asia.Figure 3e shows that the differences in SW CRF between the new and old configurations are located mainly in the intertropical convective zone (ITCZ), with maximum values of over 14 W m −2 , as abundant high-level ice clouds exist in all the ITCZ regions.Differences in SW CRF over mid-high latitudes are much smaller.There are only minor differences in areas with large SW CRF along 60 • S, where a large number of low-level clouds (mostly liquid) exist, because the liquid cloud optics in the two configurations are almost equivalent for CRF calculation.Consequently, the changes in ice cloud optical properties are the main cause of the changes in SW in the run with the new configuration when maximum-random overlap of plane-parallel horizontally homogeneous clouds is assumed.
In the OLD run, LW CRF is overestimated over most of the tropical and subtropical oceans with very few exceptions, but it is underestimated over intensively convective tropical regions such as the central Africa, the west Pacific warm pool, and the Amazon forests of South America (Fig. 3b).The N_MRO run produces similar distributions of these biases; however, the positive biases in the tropical and subtropical oceans are reduced, whereas the negative biases are enhanced somewhat (Fig. 3d).The differences in LW CRF between the new and old configurations (see Fig. 3f) show a quite similar geographic distribution to those of SW CRF (see Fig. 3e), with a maximum value of more than −9 W m −2 in the tropical east Pacific.Again, variations in ice cloud optics play a critical role in causing these differences.
The cooling effect by SW CRF and heating effect by LW CRF at the TOA in tropical deep convective regions have been shown to be nearly linearly correlated and generally compensate for each other (Kiehl and Ramanathan, 1990), which means that the (SW CRF) / (LW CRF) ratio is about −1.This ratio is often used as a criterion for showing the performance of modeled CRF.The (SW CRF) / (LW CRF) ratios in the Indonesian region (10 • S-20 • N, 110-160 • E) for the OLD and N_MRO simulations and CERES_EBAF observations are given in Table 2.The table shows that the OLD run overestimates the (SW CRF) / ( LW CRF) ratios for the annual and seasonal means.N_MRO shows a generally noticeable decrease in SW CRF, especially for the annual and summer (JJA) means.This results in decreased (SW CRF) / (LW CRF) ratios (Table 2).However, previous offline comparisons (not shown here) give very similar results for the new and old schemes.Thus, it can be inferred that the two versions of radiation scheme are comparable for diagnosing the (SW CRF) / (LW CRF) ratio, whereas the climate feedback evidently changes the simulated cloud fractions or cloud condensate.Radiative heating/cooling within the atmosphere is a critical driving factor in climate simulations.Figure 4 compares the clear-sky and all-sky LW heating rates (HRs) of N_MRO and OLD, respectively.The differences in SW HR are much smaller than those in LW HR (figure not shown).For the clear-sky condition, N_MRO shows a remarkable (more than 10 %) increased radiative cooling in the lower troposphere within 60 • S-60 • N and a reduced radiative cooling in most of the middle troposphere.These may be related to the different treatments of greenhouse gases, especially O 3 and water vapor.The difference in the all-sky LW HR (see Fig. 4f) is similar to the pattern shown in Fig. 4c, indicating that the differences in the HR of clouds are less important for determining the all-sky HR differences in this case.This pattern tends to increase the stability of the atmosphere below 600 hPa but enhance vertical mixing above 600 hPa.
As shown above, the application of the BCC-RAD radiation scheme, without tuning the subgrid cloud structures, remarkably influences the radiation budget at both boundaries and within the atmosphere.These changes will extensively affect the final simulated climate.(b, d) the differences between N_MRO and observations (red solid lines) and between OLD and observations (blue solid lines).The observation data for temperature are from the ERA-Interim reanalysis, and those for precipitation are from the Xie and Arkin (1997) data set.

Surface climatology
In this subsection, the simulated surface temperature and precipitation are evaluated.Zonal comparisons of surface temperature (ST, for land only) and precipitation rate are shown in Fig. 5.There are substantial differences between the simulations and the ERA-Interim reanalysis (Uppala et al., 2005), which is averaged over the same period as the simulations.For instance, simulated STs are underestimated by about 1.5 K in the mid-latitudes and by about 3 K around the North Pole (see Fig. 5b), but overestimated by about 2 K and 4 K over the tropics and South Pole, respectively.The global distributions of surface temperature biases for the N_MRO and OLD runs are quite similar (figures not shown), with local maximum differences between the N_MRO and OLD runs reaching ±2-4 K.The differences between the simulations and observations are much larger than the differences between the N_MRO and OLD simulations.It should be noted that the SSTs used here are prescribed, which limited the model response.
Similar to Fig. 5a and b, Fig. 5c and d show comparisons of the precipitation rate.Both the N_MRO and OLD simulations capture the main features of the meridional distribution of precipitation, such as the maximum in the tropics and secondary maxima at the mid-latitude storm tracks.However, errors are also clear relative to the observation, especially in the tropics.The two simulations are comparable in the simulation of the zonal mean distribution of precipitation, but there are noticeable local differences in the tropical and subtropical regions (figures not shown).These differences probably stem from the altered atmospheric thermodynamics and dynamics caused by the changes in radiation budget.The increases and decreases in precipitation often coincide with the decreases and increases in surface temperature (figures not shown), respectively; thus, the changes in precipitation also obviously influence the surface energy balance.

Atmospheric states
Simulated atmospheric temperature and specific humidity are analyzed in this subsection.
Figure 6 shows the comparisons of the latitude-height distribution of atmospheric temperature and specific humidity.Notable cold biases relative to ERA-Interim, about 1-2 K in the low-mid troposphere, exist throughout almost the entire troposphere in the OLD case (see Fig. 6a).The N_MRO simulation inherits most of these biases, but the relative warming (up to 0.4-0.8K) within the middle troposphere (800 ∼ 500 hPa) is a desirable change compared with OLD (see Fig. 6c).This is definitely related to the reduced LW cooling rate in the middle troposphere, as shown in Fig. 4.
In addition to the improvements in tropospheric temperature, there are favorable changes in specific humidity.Compared with ERA-Interim, the OLD run is subject to considerable dry biases in the tropical lower troposphere (see Fig. 6d).This is likely caused by LW heating rate biases related to the LW parameterization of water vapor (Collins et al., 2002).Due to changes in heating rate, as shown in Fig. 4 The changes in atmospheric temperature and specific humidity exert influences on the formation and maintenance of cloud water and ice (figures not shown), affecting the modeled local radiation budget, such as altering the (SW CRF) / (LW CRF) ratio mentioned above.
Overall, the incorporation of the new scheme influences radiative fluxes and heating rates remarkably.Due to these changes, the simulated surface and atmospheric climate are comparable or improved relative to the old model configuration.Therefore, the new scheme used here has been demonstrated to be a viable option for long-term climate simulation.
It should also be mentioned that the differences in simulated climate between the two model configurations are relatively smaller than those between the simulations and observations.Nevertheless, the much more flexible cloud structure and internal consistency of the new configuration will benefit further development of model physics.In regard to the convenience of the McICA scheme for modifying subgrid clouds, the impacts of the cloud structure variations are assessed as follows.

The impacts of altering subgrid cloud structures
In this subsection, we discuss the impacts of altering overlapping scheme of fractional clouds and breaking the traditional PPH assumption on simulated radiation and climate.

The impact of altering cloud overlap
In Fig. 7, the top panels show the zonal mean distributions of the fractional contribution by deep convective cloud to total cloud (f con /f tot ) in GCM grid cells.f con /f tot has the maximum value in the upper tropics, then decreases sharply with latitude.There is a seasonal variation as shown in DJF and JJA.In the second row of Fig. 7, the latitude-height distributions of the grid-cell mean L cf calculated using Eq. ( 8) are shown.In general, L cf reaches its maximum value in upper tropics, and decreases with latitude, similar to the distribution of f con /f tot .That the decorrelation length tends to increase upwards has also been alluded by Barker (2008) and Räisänen et al. (2004).
In the bottom panels of Fig. 7, the zonal mean vertically averaged L cf (solid lines) is shown.Also shown is the Shonk et al. (2010) function result based merely on latitude (dashed lines).It is found that L cf not only depends on geographical location but also has seasonal variation as indicated in DJF and JJA.Therefore, it is difficult to parameterize L cf only based on latitude.Though the scheme of N_GOF is not directly from observations, the result shown is very close to that derived from CloudSat/CALIPSO observations (see Fig. 1 in Oreopoulos et al., 2012).Through comparison with the standard scheme of N_GO2, the impact of cloud-typerelated L cf can be explored.
In the upper panels of Fig. 8, the changes in total cloud fraction between N_GOF and N_GO2 are shown for DJF and JJA.Compared to N_GO2, N_GOF yields fewer clouds in low latitudes of the summer hemisphere, and generally produces more clouds in the middle and high latitudes.This is clearly seen in the zonal mean results presented in the lower panels (black lines).In the regions with more strongly convective clouds, the larger L cf causes a larger possibility for maximum overlap, which leads to less cloud fraction.The same argument applies to the regions with fewer convective clouds.In the lower panels, the changes in zonal mean cloud fraction by N_GOF-N_GO2 for low, middle and high clouds are presented separately.For the low and middle clouds, it is found that the cloud fraction generally decreases in the tropics and increases in the middle and high latitudes, which is consistent with the total cloud fraction.However, for high clouds, N_GOF generally produces a larger cloud fraction compared to N_GO2 even in the tropics.As the cloud fractions in individual layers do not change notably, the increased high cloud fraction is mostly because the decorrelation length in the upper troposphere is smaller in N_GOF than in N_GO2.Fewer low/middle clouds in the tropics causes more solar radiation to reach the lower atmosphere and surface, because high clouds have smaller impact on solar radiation compared to the low/middle clouds.The extra solar heating in the lower atmosphere and surface can enhance the surface evaporation and atmospheric convection, which then can produce more cirrus clouds (Emmanuel, 1994).
In the top panels of Fig. 9, the differences in SW CRF between N_GOF and N_GO2 are shown.The differences in SW CRF can exceed 15 W m −2 over a large domain in the tropical regions.In the regions with more convective clouds (for instance, the summer northern Indian Ocean and tropical Atlantic Ocean), the larger chance of maximum overlap leads to less cloud solar reflection, and thereby less SW CRF (positive differences between N_GOF and N_GO2).Similarly, in the regions with fewer convective clouds (for instance, the summer Arctic and winter Southern Hemisphere oceans), the difference between the two schemes becomes negative.The zonal mean results show that SW CRF generally increases in the tropics and decreases in the middle and high latitudes.
The middle panels of Fig. 9 show the differences in LW CRF between N_GOF and N_GO2.In the regions with more convective clouds, the smaller cloud fraction shades less longwave radiation from reaching outer space, thereby causing less LW CRF.Similarly, in the regions of fewer convective clouds, LW CRF is increased.The zonal mean results show that the LW CRF is generally increased by N_GOF in the middle and high latitudes.However the change in the tropics is small.LW CRF is strongly dependent on high clouds.As shown in Fig. 8, the change in high cloud fraction in the tropics is different from that in low/middle clouds.
The changes of SW and LW CRFs are generally opposite.For example, over the summer northern Indian Ocean and tropical Atlantic Ocean, a positive SW CRF difference corresponds to a negative LW CRF difference.However, the solar effect is generally stronger than that of longwave.The net effect in the total CRF does not cancel out.In the lower panels of Fig. 9, the differences in total CRF are shown.From the zonal mean result, a more strongly positive CRF in the summer Arctic is clearly seen.Despite the obvious differences of CRFs at the TOA in different latitudes, the global mean differences in net fluxes are small (Fig. 1b, h) at the TOA, due to compensation between latitudes.
By using diagnostic calculations, Shonk and Hogan (2010) and Oreopoulos et al. (2012) found that more maximumlike cloud overlap leads to decreases in the magnitude of SW and LW CRFs with local maximum shift on an order of 10 W m −2 , which is consistent with our convectivestratiform contrast scheme.One advantage of the clouddependent decorrelation length described in this article is that it can respond to climate change (i.e., if the distribution of convective cloud changes in a climate change experiment, so will the decorrelation length).This is not the case for the latitude-dependent and Julian-day-dependent decorrelation parameterizations of Shonk and Hogan (2010) and Oreopoulos et al. (2012).
In Fig. 10, the top panels show the change in net radiative flux at the TOA, which is equivalent to the change in energy balance at the TOA.It is found that the energy balance at the TOA has been considerably affected by addressing the cloudtype-related L cf .N_GOF leads to an obvious increase of net flux in the tropics (especially for the JJA season) and highlatitude regions.In the subtropical high-pressure regions, the change in net flux is relatively small because the cloud fraction is low there.The change in energy balance from the tropics to the subtropical high regions can influence Hadley circulation.The same as the result at TOA, the largest changes in surface energy balance occur in the tropics and high latitudes.It is worth emphasizing that the change in energy balance in Arctic is large, which should have very strong impact on the evolution of sea ice.Investigating the impact on polar climate by N_GOF will be a subsequent work for us.
Figure 11a and b are the land surface temperature differences between N_GOF and N_GO2.The patterns are similar to those of the differences in net CRF at the surface as shown in Fig. 10c and d over land.For DJF, most land surfaces are simulated warmer by N_GOF than by N_GO2, as net CRF over these areas is more positive for N_GOF.The JJA season sees more negative differences between N_GOF and N_GO2.The range of differences is generally between −3 and 3 K.

The impact of breaking the PPH assumption
In this subsection, we briefly consider the impact of breaking the traditional PPH assumption on the simulated radiation and climate by comparing the N_GOF_IH and N_GOF tests, where N_GOF_IH employs the same vertical overlap scheme as N_GOF but with consideration of the horizontal inhomogeneity in cloud condensate as discussed in Sect.3.
From a global mean perspective, the changes in CRFs and net fluxes caused by including horizontally inhomogeneous cloud condensate (Eq.9) are as large as or even larger than the changes from altering the cloud overlap assumptions (see N_GOF_IH in Fig. 1), which has also been shown by calculations from cloud-resolving models (Barker and Räisänen, 2005).The global mean reductions in LW and SW CRFs at the TOA are about 1.5 and 3.7 W m −2 , respectively.These are a little smaller than the 2.02 and 6.15 W m −2 in Shonk and Hogan (2010) using the ERA-40 data set, but larger than the 0.7 and 2.2 W m −2 in Räisänen et al. (2007)   ECHAM5 data set.Both of these two are offline calculations; thus no climate interactions are included.The consideration of horizontally inhomogeneous clouds here brings the global mean CRFs and F net much closer to observations.In general, N_GOF_IH is the best among all the simulations in terms of the global mean radiative fluxes.However, the spatial distributions of CRF biases of N_GOF_IH are similar to those of other simulations (see Fig. 12), which requires revision of not only the radiation scheme but also the other parts of the model.

using the
Figure 13 shows the differences in LW and SW CRFs between N_GOF_IH and N_GOF for DJF and JJA.It is seen that the inclusion of cloud horizontal inhomogeneity weakens the negative SW CRF nearly everywhere, with local maximum increase of more than 20 W m −2 over the western Pacific warm pool region.This is qualitatively consistent with the well-accepted conclusion that the PPH assumption of Fig. 12.The same as Fig. 3 (c, d), but for the simulation N_GOF_IH.cloud condensate generally overestimates solar reflectance (Carlin et al., 2002).The zonal mean results show that enhancements of SW CRF mostly happen in high latitudes for both the boreal summer and austral summer.Barker and Räisänen (2005) got very similar zonal distribution of SW sensitivity to cloud inhomogeneity by computing on a CSRM data set, with zonal maximum of as large as ∼ 25 W m −2 around 60 • S.
The inclusion of horizontal inhomogeneity of cloud condensate generally reduces LW CRF all over the globe (Fig. 13c, d).This is because the PPH assumption of cloud condensate generally overestimates the LW emissivity upwards (Pomroy and Illingworth, 2000).The patterns of change in SW and LW CRFs are very similar.However, the zonal mean results of LW CRF do not show very special features in the boreal or austral summer.This is also similar to results of Barker and Räisänen (2005).LW CRF shows less dependence on season probably because LW CRF does not increase proportionally with the incoming solar flux at the TOA.
The magnitude of the local change of CRFs by addressing cloud inhomogeneity is about the same as changing L cf (compared with Fig. 10).Thus, it is of great importance to address the cloud water/ice horizontal distribution together with the overlap of fractional clouds in GCMs.
The consideration of cloud horizontal inhomogeneity has a noticeable influence on surface temperature (see Fig. 11c, d).
In the middle and high latitudes of the Northern Hemisphere, there is a remarkable decrease in surface temperature during DJF and an increase during JJA.These changes arise mainly from competition between LW cooling and SW heating.When inhomogeneous clouds succeed homogeneous ones, more LW flux is emitted outward, and more SW flux penetrates to the surface (see Fig. 1).The surface energy budget is then the net effect of the two fluxes.
In general, the modifications of cloud subgrid configurations have distinct impacts on the simulated radiation budget and surface temperature.Considering the improved LW and SW budget for clear-sky and all-sky conditions, the new model configuration can be used in BCC_AGCM2.0.1 to improve physical processes and perform climate simulations.
It should be noted that, in this work, we altered only the subgrid cloud structures used in the radiation calculation, whereas those in precipitation parameterization were not changed.Physically, cloud overlap assumptions in the radiation and precipitation processes should be consistent with each other, but the latter may have a larger effect on the simulated precipitation (Morcrette and Jakob, 2000).However, this is beyond the scope of this study.

Discussion and conclusions
In this work, the BCC-RAD radiation algorithm was incorporated into the BCC_AGCM2.0.1 GCM as a replacement for the original radiation algorithm.The new radiation is entirely distinct from the original one, by including the advanced treatments of gaseous absorption, cloud optics and radiation transfer.The results show that the new scheme markedly improves the representation of the SW and LW radiation budget for both clear-sky and all-sky conditions, whether in the global mean or in geographic distribution.The simulated relationship between SW and LW CRFs in deep convective regions is improved by the new scheme as well.The modeled temperature and specific humidity benefited from changes in the LW heating rate, resulting in a reduction in temperature biases by 0.4-0.8K at the middle troposphere and a reduction in moisture biases by one-third in the tropical lower troposphere relative to the ERA-Interim reanalysis.This shows the superiority of a CKD radiation algorithm to the bandmodel-based CAM3 radiation scheme.All these results indicate that using the BCC-RAD radiation scheme makes the cloud-radiation process more intrinsically coherent and reasonable.
The McICA scheme was applied to the BCC-RAD radiation algorithm to deal with cloud subgrid variability.This new configuration is flexible for treating arbitrarily complex subgrid cloud structures, including the vertical overlap of fractional clouds and the horizontal distribution of cloud condensate.The impact of altering cloud overlap within the McICA scheme was assessed by including a so-called "general overlap" with a global constant L cf of 2 km.In this work we have also considered a cloud-type-related L cf , which especially distinguishes the deep convective cloud from other types of clouds.This scheme yields a similar zonal mean distribution of L cf compared to the CloudSat/CALIPSO result demonstrated by Oreopoulos et al. (2012).One advantage of the cloud-dependent decorrelation length described in this article is that it can respond to climate change.The new cloudtype-related L cf has remarkable effect on local radiation budget compared to the scheme of constant L cf .In the regions of frequent convection such as South Asia monsoon region and ITCZ, CRF is most notably influenced, with local differences of > 20 W m −2 for SW CRF and > 10 W m −2 for LW CRF.
Generally, more net energy reaches surface around the tropics while less reaches the surface in mid-high latitudes.This could enhance the simulated Hadley circulation.Therefore, a constant value of L cf could lead to large biases in climate simulations.
The effect of the horizontal inhomogeneity of cloud condensate was then considered by including an observationbased gamma function in an additional test.The changes compared with its PPH counterpart were strikingly significant, with decreases in the global mean TOA LW and SW CRFs of ∼ 1.5 W m −2 and ∼ 3.7 W m −2 , respectively, making the simulation much closer to observations.This emphasizes the importance of addressing the cloud horizontal distribution in GCMs along with the cloud overlap issue.
For simulated climate, the changes in cloud structures showed a notable effect on surface temperatures in mid-high latitudes, with the largest zonal mean differences being about 1 K, exceeding the differences between the new and old radiation scheme configurations when both use the maximumrandom overlap and PPH assumptions.This will exert an effect on the evolution of sea ice.In this work, we did not include sea-atmosphere interaction, which could enlarge the effect of the signal imposed by the changing cloud structures.
This study shows that N_GOF_IH is the current best choice for climate simulations, and it will be adopted in BCC_AGCM2.0.1.
This study also provided a direction for future improvement of the McICA.More realistic cloud overlap assumption or cloud horizontal distribution, achieved from satellite observations or any other objective sources, could be used to constrain the generation of cloud-type-related L cf and the horizontal distribution of cloud condensate.
developed by the Beijing Climate Center (BCC) at the China Meteorological Administration (CMA) based on the Community Atmosphere Model Version 3 (CAM3) of the National Center for Atmospheric Research (NCAR)

Fig. 1 .
Fig. 1.Global annual mean clear-sky net fluxes (F clr , left panels), all-sky net fluxes (F net , central panels), and CRFs (right panels) for TOA LW (upmost row), surface LW (second row), TOA SW (third row), and surface SW (bottom row) from the various simulations and CERES_EBAF observations.The error bars are the ranges between the maximum and minimum annual mean values for the simulated or observed decades.

Fig. 3 .
Fig. 3.The annual mean differences in SW CRF (left) and LW CRF (right).Differences larger (smaller) than 10 (−10) W m −2 are shaded in yellow (blue).The global mean and root mean square values of each figure are also shown.

Fig. 4 .
Fig. 4. Zonal annual mean clear-sky and all-sky LW heating rates for OLD (a, d), N_MRO (b, e) and the differences between N_MRO and OLD (c, f) (units: K day −1 ).

Fig. 5 .
Fig. 5. Zonal annual mean (a) surface temperature and (c) precipitation from N_MRO, OLD, and observation data, as well as(b, d) the differences between N_MRO and observations (red solid lines) and between OLD and observations (blue solid lines).The observation data for temperature are from the ERA-Interim reanalysis, and those for precipitation are from theXie and Arkin (1997)  data set.

Fig. 6 .
Fig. 6.Biases in zonal annual mean atmospheric temperature and specific humidity compared with the ERA-Interim reanalysis for (a, d) OLD and (b, e) N_MRO simulations and (c, f) the differences between N_MRO and OLD.

Fig. 7 .
Fig. 7. Zonal annual mean (a, b) fractional contribution by deep convective cloud fraction (f con /f tot ) and (c, d) L cf calculated using Eq.(8) from the N_GOF run, as well as (e, f) the vertically averaged L cf for N_GOF (solid lines) and that used by Shonk et al. (2010) (dashed lines), for DJF (left) and JJA (right).

Fig. 10 .
Fig. 10.Differences in net radiative budget (a, b) at the TOA and (c, d) at the surface between N_GOF and N_GO2, for DJF (left) and JJA (right).

Table 1 .
Comparison of the new and old schemes.