4.8 Parameterization of Shortwave Radiation

4.8.1 Diurnal cycle

With standard name-list settings, both the longwave and shortwave heating rates are evaluated every model hour. Between hourly evaluations, the longwave and shortwave fluxes and flux divergences are held constant.

In CAM 3.0, insolation is computed using the method of Berger [17]. Using this formulation, the insolation can be determined for any time within $10^6$ years of 1950 AD. This facilitates using CAM 3.0 for paleoclimate simulations. The insolation at the top of the model atmosphere is given by

$$S_I = {S_0 {\rho^{-2}} \cos\mu} ,$$ (4.174)

where $S_0$ is the solar constant, $\mu$ is the solar zenith angle, and $\rho^{-2}$ is the distance factor (square of the ratio of mean to actual distance that depends on the time of year). In the standard configuration, $S_0 = 1367.0$ W/m${}^2$. CAM 3.0 includes a mechanism for treating the slow variations in the solar constant over the 11-year cycle and during longer secular trends. A time series of $S_0$ for 1870-2100 based upon Lean et al. [104] is included with the standard model.

We represent the annual and diurnal cycle of solar insolation with a repeatable solar year of exactly 365 days and with a mean solar day of exactly 24 hours, respectively. The repeatable solar year does not allow for leap years. The expressions defining the annual and diurnal variation of solar insolation are:

$$\cos\mu = \sin\phi \sin\delta - \cos\phi \cos\delta \cos(H)$$ (4.175)

$$\delta = \arcsin(\sin\epsilon\sin\lambda)$$ (4.176)

$$\rho = \frac{1-e^2}{1+e \cos(\lambda - \tilde\omega )}$$ (4.177)

$$\tilde\omega = \Pi + \psi$$ (4.178)
where

$$\phi = {\rm latitude in radians} \nonumber$$

$$\delta = {\rm solar declination in radians}$$

$$H = {\rm hour angle of sun during the day}$$

$$\epsilon = {\rm obliquity}$$

$$\lambda = {\rm true longitude of the earth relative to vernal equinox}$$ (4.179)

$$e = {\rm eccentricity factor}$$

$$\tilde\omega = {\rm longitude of the perihelion}+ 180^\circ$$

$$\Pi = {\rm longitude of perihelion based on the fixed equinox}$$

$$\psi = {\rm general precession} .$$

Note that $\Pi$ is denoted by $\pi$ in Berger [17].

The hour angle $H$ in the expression for $\cos\mu$ depends on the calendar day $d$ as well as model longitude:

$$H = 2 \pi\left(d + \frac{\theta}{360^\circ}\right) ,$$ (4.180)

where $\theta$ = model longitude in degrees starting from Greenwich running eastward. Note that the calendar day $d$ varies continuously throughout the repeatable year and is updated every model time step. The values of $d$ at 0 GMT for January 1 and December 31 are 0 and 364, respectively. This would mean, for example, that a model calendar day $d$ having no fraction (such as 182.00) would refer to local midnight at Greenwich, and to local noon at the date line ($180^\circ$ longitude).

The obliquity $\epsilon$ may be approximated by an empirical series expansion of solutions for the Earth's orbit

$$\epsilon = \epsilon^* + \sum_{j=1}^{47} A_j \cos\left(f_j t + \delta_j\right) \$$ (4.181)

where $A_j$, $f_j$, and $\delta_j$ are determined by numerical fitting. The term $\epsilon^* = 23.320556^\circ$, and $t$ is the time (in years) relative to 1950 AD.

Since the series expansion for the eccentricity $e$ is slowly convergent, it is computed using

$$e = \sqrt{\left(e\cos\Pi\right)^2+\left(e\sin\Pi\right)^2}$$ (4.182)

The terms on the right-hand side may also be written as empirical series expansions:

$$e\left\lbrace\begin{array}{c}\cos\ \sin\end{array}\right\rbrace ... ...\begin{array}{c}\cos\ \sin\end{array}\right\rbrace \left(g_j t+\beta_j\right)$$ (4.183)

where $M_j$, $g_j$, and $\beta_j$ are estimated from numerical fitting. Once these series have been computed, the longitude of perihelion $\Pi$ is calculated using

$$\Pi = \arctan\left(\frac{e \sin\Pi}{e \cos\Pi}\right)$$ (4.184)

The general precession is given by another empirical series expansion

$$\psi = \tilde\psi t + \zeta + \sum_{j=1}^{78} F_j \sin\left(f'_j t + \delta'_j\right)$$ (4.185)

where $\tilde\psi = 50.439273''$, $\zeta = 3.392506^\circ$, and $F_j$, $f'_j$, and $\delta'_j$ are estimated from the numerical solution for the Earth's orbit.

The calculation of $\lambda$ requires first determining two mean longitudes for the orbit. The mean longitude $\lambda_{m0}$ at the time of the vernal equinox is :

$$\lambda_{m0} = 2\left\lbrace \left(\frac{e}{2} + \frac{e^3}{8}\right) (1+\beta)\sin(\tilde\omega ) \right.$$

$$\phantom{-2\left\lbrace\right.} -\frac{e^2}{4} \left(\frac{1}{2}+\beta\right) \sin(2 \tilde\omega )$$ (4.186)

$$\phantom{-2\left\lbrace\right.} +\left.\frac{e^3}{8} \left(\frac{1}{3}+\beta\right) \sin(3 \tilde\omega ) \right\rbrace$$

where $\beta = \sqrt{1-e^2}$. The mean longitude is

$$\lambda_m = \lambda_{m0} + \frac{2 \pi (d-d_{ve})}{365}$$ (4.187)

where $d_{ve}=80.5$ is the calendar day for the vernal equinox at noon on March 21. The true longitude $\lambda$ is then given by:

$$\lambda = \lambda_m + \left(2 e - \frac{e^3}{4}\right)\sin(\lambda_m-\tilde\omega )$$

$$+ \frac{5 e^2}{4} \sin\left[2(\lambda_m-\tilde\omega )\right]$$ (4.188)

$$+ \frac{13 e^3}{12}\sin\left[3(\lambda_m-\tilde\omega )\right]$$

The orbital state used to calculate the insolation is held fixed over the length of the model integration. This state may be specified in one of two ways. The first method is to specify a year for computing $t$. The value of the year is held constant for the entire length of the integration. The year must fall within the range of $1950 \pm 10^6$. The second method is to specify the eccentricity factor $e$, longitude of perihelion $\tilde\omega - 180^\circ$, and obliquity $\epsilon$. This set of values is sufficient to specify the complete orbital state. Settings for AMIP II style integrations under 1995 AD conditions are $\epsilon = 23.4441$, $e = 0.016715$, and $\tilde\omega - 180 = 102.7$.

4.8.2 Formulation of shortwave solution

The $\delta$-Eddington approximation of Joseph et al. [80] and Coakley et al. [35] has been adopted and is described in Briegleb [27]. This approximation has been shown to simulate quite well the effects of multiple scattering. The major differences between the shortwave parameterizations in CCM3 and CAM 3.0 are

1. the new treatment of cloud vertical overlap [38];
2. updated parameterization for near-infrared absorption by water vapor; and
3. inclusion of prescribed aerosol data sets for computing shortwave aerosol radiative forcing.

The solar spectrum is divided into 19 discrete spectral and pseudo-spectral intervals (7 for $\mathrm{O_{3}}$, 1 for the visible, 7 for $\mathrm{H_{2}O}$, 3 for $\mathrm{CO_{2}}$, and 1 for the near-infrared following Collins [37]). The CAM 3.0 model atmosphere consists of a discrete vertical set of horizontally homogeneous layers within which radiative heating rates are to be specified (see Figure 3.1). Each of these layers is considered to be a homogeneous combination of several radiatively active constituents. Solar irradiance, surface reflectivity for direct and diffuse radiation in each spectral interval, and the cosine of the solar zenith angle are specified. The surface albedo is specified in two wavebands (0.2-0.7 $\mu$m, and 0.7-5.0 $\mu$m) and distinguishes albedos for direct and diffuse incident radiation. Albedos for ocean surfaces, geographically varying land surfaces, and sea ice surfaces are distinguished.

The method involves evaluating the $\delta$-Eddington solution for the reflectivity and transmissivity for each layer in the vertical under clear and overcast conditions. The layers are then combined together, accounting for multiple scattering between layers, which allows evaluation of upward and downward spectral fluxes at each interface boundary between layers. This procedure is repeated for each spectral or pseudo-spectral interval and binary cloud configuration (see ``Cloud vertical overlap'' below) to accumulate broad band fluxes, from which the heating rate can be evaluated from flux differences across each layer. The $\delta$-Eddington scheme is implemented so that the solar radiation is evaluated once every model hour (in the standard configuration) over the sunlit portions of the model earth.

The $\delta$-Eddington approximation allows for gaseous absorption by $\mathrm{O_{3}}$, $\mathrm{CO_{2}}$, $\mathrm{O_{2}}$, and $\mathrm{H_{2}O}$. Molecular scattering and scattering/absorption by cloud droplets and aerosols are included. With the exception of $\mathrm{H_{2}O}$, a summary of the spectral intervals and the absorption/scattering data used in the formulation are given in Briegleb [27] and Collins [37]. Diagnostic cloud amount is evaluated every model hour just prior to the solar radiation calculation.

The absorption by water vapor of sunlight between 1000 and 18000 cm${}^{-1}$ is treated using seven pseudo-spectral intervals. A constant specific extinction is specified for each interval. These extinctions have been adjusted to minimize errors in heating rates and flux divergences relative to line-by-line (LBL) calculations for reference atmospheres [3] using GENLN3 [57] combined with the radiative transfer solver DISORT2 [163]. The coefficients and weights have the same properties as a k-distribution method [98], but this parameterization is essentially an exponential sum fit (e.g., Wiscombe and Evans [191]). LBL calculations are performed with the HITRAN2k line database [153] and the Clough, Kneizys, and Davies (CKD) model version 2.4.1 [33]. The Rayleigh scattering optical depths in the seven pseudo-spectral intervals have been changed for consistency with LBL calculations of the variation of water-vapor absorption with wavelength. The updated parameterization increases the absorption of solar radiation by water vapor relative to the treatment used in CCM and CAM since its introduction by Briegleb [27].

For some diagnostic purposes, such as estimating cloud radiative forcing [94] a clear-sky absorbed solar flux is required. In CAM 3.0, the clear-sky fluxes and heating rates are computed using the same vertical grid as the all-sky fluxes. This replaces the 2-layer diagnostic grid used in CCM3.

4.8.3 Aerosol properties and optics

4.8.3.1 Introduction

The treatment of aerosols in CAM 3.0 replaces the uniform background boundary-layer aerosol used in previous versions of CAM and CCM. The optics for the globally uniform aerosol were identical to the sulfate aerosols described by Kiehl and Briegleb [88]. In the visible, the uniform aerosol was essentially a conservative scatterer. The new treatment introduces five chemical species of aerosol, including sea salt, soil dust, black and organic carbonaceous aerosols, sulfate, and volcanic sulfuric acid. The new aerosols include two species, the soil dust and carbonaceous types, which are strongly absorbing in visible wavelengths and hence increase the shortwave diabatic heating of the atmosphere.

The three-dimensional time-dependent distributions of the five aerosol species and the optics for each species are loaded into CAM 3.0 during the initialization process. This provides considerable flexibility to:

• Change the speciated aerosol climatology / time-series as aerosol modeling improves;
• Vary the aerosol distributions for climates different from present-day conditions;
• Examine the effects of individual aerosol species and arbitrary combinations of aerosol species; and
• Change aerosol optical properties.

In its present configuration, CAM includes the direct and semi-direct effects of tropospheric aerosols on shortwave fluxes and heating rates. The first indirect effect, or Twomey et al. [175] effect, is not included in the standard version of CAM 3.0.

4.8.3.2 Description of aerosol climatologies and data sets

The data sets for the tropospheric and stratospheric aerosols are treated separately in the model.

The annually-cyclic tropospheric aerosol climatology consists of three-dimensional, monthly-mean distributions of aerosol mass for:

• sulfate from natural and anthropogenic sources;
• sea salt;
• black and organic carbon derived from natural and anthropogenic sources; and
• soil dust

There are four size categories of dust spanning diameters from 0.01 to 10 $\mu$m, and the black and organic carbon are represented by two tracers each for the hydrophobic (new) and hydrophilic (aged) components. The climatology therefore contains ten types of aerosol: sea salt, four size bins of soil dust, sulfate, new and aged black carbon, and hydrophobic and hydrophilic organic carbon.

The climatology is produced using an aerosol assimilation system [39,41] integrated for present-day conditions. The system consists of the Model for Atmospheric Chemistry and Transport (MATCH) [142] and an assimilation of satellite retrievals of aerosol optical depth. MATCH version 4 is integrated using the National Centers for Environmental Prediction (NCEP) meteorological reanalysis at T63 triangular truncation [83]. The satellite estimates of aerosol optical depth are from the NOAA Pathfinder II data set [167].

The formulation of the sulfur cycle is described in Barth et al. [12] and Rasch et al. [143]. The emissions inventory for SO${}_2$ is from Smith et al. [162]. The sources for mineral dust are based upon the approach of Zender et al. [196] and Mahowald et al. [120]. The emissions of carbonaceous aerosols include contributions from biomass burning [113], fossil fuel burning [42], and a source of natural organic aerosols resulting from terpene emissions. The vertical profiles of sea salt are computed from the 10m wind speed [21].

The monthly-mean mass path for each aerosol species in each layer is computed in units of kg/m${}^2$. During the initialization of CAM 3.0, the climatology is temporally interpolated from monthly-mean to mid-month values. At each CAM 3.0 time step, the mid-month values bounding the current time step are vertically interpolated onto the pressure grid of CAM 3.0 and then time interpolated to the current time step. The interpolation scheme in CAM 3.0 preserves the aerosol masses for each species to 1 part in 10${}^7$ relative to the climatology, and it is guaranteed to yield positive definite mass-mixing ratios for all aerosols.

The stratospheric volcanic aerosols are treated using a single species in the standard model. Zonal variations in the stratospheric mass loading are omitted. The volcanic input consists of the monthly-mean masses in units of kg/m${}^2$ on an arbitrary meridional and vertical grid. The time series for the recent past is based upon Ammann et al. [2] following Stenchikov et al. [166].

4.8.3.3 Calculation of aerosol optical properties

The three intrinsic optical properties stored for each of the eleven aerosol types are specific extinction, single scattering albedo, and asymmetry parameter. These properties are computed on the band structure of CAM 3.0 using Chandrasekhar weighting with spectral solar insolation. The aerosol types affected by hygroscopic growth are sulfate, sea salt, and hydrophilic organic carbon. In previous versions of CCM and CAM 3.0, the relative humidity was held constant in calculations of hygroscopic growth at 80%. In CAM 3.0, the actual profiles of relative humidity computed from the model state each radiation time step are used in the calculation.

The optics for black and organic carbon are identical to the optics for soot and water-soluble aerosols in the Optical Properties of Aerosols and Clouds (OPAC) data set [69]. The optics for dust are derived from Mie calculations for the size distribution represented by each size bin [196]. The Mie calculations for sulfate assume that it is comprised of ammonium sulfate with a log-normal size distribution. The dry size parameters are a median radius of 0.05 $\mu$m and a geometric standard deviation of 2.0. The optical properties in the seven H${}_2$O pseudo-spectral intervals are averaged consistently with LBL calculations of the variation of water-vapor absorption with wavelength. This averaging technique preserves the cross correlations among the spectral variation of solar insolation, water vapor absorption, and the aerosol optical properties. The volcanic stratospheric aerosols are assumed to be comprised of 75% sulfuric acid and 25% water. The log-normal size distribution has an effective radius of 0.426 $\mu$m and a standard deviation of 1.25.

The bulk formulae of Cess [30] are used to combine the optical properties of the individual aerosol species into a single set of bulk aerosol extinctions, single-scattering albedos, and asymmetry parameters for each layer.

4.8.3.4 Calculation of aerosol shortwave effects and radiative forcing

CAM 3.0 includes a mechanism to scale the masses of each aerosol species by user-selectable factors at runtime. These factors are global, time-independent constants. This provides the flexibility to consider the climate effects of an arbitrary combination of the aerosol species in the climatology. It also facilitates simulation of climates different from present-day conditions for which the only information available is the ratio of globally averaged aerosol emissions or atmospheric loadings. A mechanism to scale the carbonaceous aerosols with a time-dependent unitless factor has been included to facilitate realistic simulations of the recent past.

CAM 3.0 also includes a run-time option for computing a diagnostic set of shortwave fluxes with an arbitrary combination of aerosols multiplied with a separate set of user-selectable scale factors. This option can be used to compute, for example, the aerosol radiative forcing relative to an atmosphere containing no aerosols.

The diagnostic fields produced the aerosol calculation include the column-integrated optical depth and column-averaged single-scattering albedo, asymmetry parameter, and forward scattering parameter (in the $\delta$-Eddington approximation) for each aerosol species and spectral interval. These fields are only computed for illuminated grid points, and for non-illuminated points the fields are set to zero. The fraction of the time that a given grid point is illuminated is also recorded. Time averages of, for example, the optical depth can be obtained by dividing the time-averaged optical depths in the history files by the corresponding daylit fractions.

4.8.3.5 Globally uniform background sulfate aerosol

The option of introducing a globally uniform background sulfate aerosol is retained, although by default the optical depth of this aerosol is set to zero. Its optical properties are computed using the same sulfate optics as are used for the aerosol climatology. However, for consistency with the uniform aerosol in previous versions of CAM 3.0 and CCM3, the relative humidity used to compute hygroscopic growth is set to 80%.

4.8.4 Cloud Optical Properties

4.8.4.1 Parameterization of effective radius

Observational studies have shown a distinct difference between maritime, polar, and continental effective cloud drop size, $r_e$, for warm clouds. For this reason, CAM 3.0 differentiates between the cloud drop effective radius for clouds diagnosed over maritime and continental regimes [89], and over pristine surfaces (sea ice, snow covered land). Over the ocean, the cloud drop effective radius for liquid water clouds, $r_{el}$, is specified to be 14$\mu$m. Over sea ice, where we presume pristine conditions, $r_{el}$ is also specified to be 14$\mu$m. Over land masses $r_{el}$ is determined using

$$r_{el} = \begin{cases}8 \; \mu{\rm m} & -10^\circ C < T \ 8 - 6(... ...^\circ C \le T \le -10^\circ C\ 14 \; \mu{\rm m} & -30^\circ C > T \end{cases}$$ (4.189)

This does not necessarily correspond to the range over which the cloud ice fraction increases from 0 to 1. In addition, $r_{el}$ ramps linearly toward the pristine value of 14$\mu$m as water equivalent snow depth over land goes from 0 to 0.1 m.

An ice particle effective radius, $r_{ei}$, is also diagnosed by CAM 3.0. Following Kristjánsson and Kristiansen [97], the effective radius for ice clouds is now a function only of temperature, as shown in Figure 4.2.

Figure 4.2: Ice effective radius and terminal velocity. Top, ice effective radius versus temperature. Bottom, ice velocity versus radius (left) and temperature (right); the Stokes terminal velocity is solid and the actual velocity is dashed.

4.8.4.2 Dependencies involving effective radius

For cloud scattering and absorption, the radiative parameterization of Slingo [160] for liquid water droplet clouds is employed. In this parameterization, the optical properties of the cloud droplets are represented in terms of the prognosed cloud water path (CWP, in units of kg m$^{-2}$) and effective radius $r_e = \int r^3 n (r) dr / \int r^2 n (r) dr$, where $n(r)$ is the cloud drop size distribution as a function of radius $r$.

Cloud radiative properties explicitly account for the phase of water. For shortwave radiation we use the following generalization of the expression used by Slingo [160] for liquid water clouds. The cloud liquid optical properties (extinction optical depth, single scattering albedo, asymmetry parameter and forward scattering parameter) for each spectral interval are defined as

$$\tau _l^c = CWP\left[ a_l^{i}+\frac{b_l^i}{r_{el}} \right] \left( 1-f_{ice} \right)$$ (4.190)

$$\omega _l^c = 1-c_l^i-d_l^ir_{el}$$ (4.191)

$$g_l^c = e_l^i+f_l^ir_{el}$$ (4.192)

$$f_l^c = \left( {g_l^c} \right)^2$$ (4.193)

where superscript $i$ denotes spectral interval. The spectral intervals and coefficients for liquid water are defined in Slingo [160].

The radiative properties of ice cloud are defined by

$$\tau_i^c = CWP\left[ a_i^i+ \frac{b_i^i}{r_{ei}} \right] f_{ice}$$ (4.194)

$$\omega _i^c = 1-c_i^i-d_i^ir_{ei}$$ (4.195)

$$g_i^c = e_i^i+f_i^ir_{ei}$$ (4.196)

$$f_i^c = \left( {g_i^c} \right)^2$$ (4.197)

where the subscript $i$ denotes ice radiative properties. The values for the coefficients $a_i-f_i$ are based on the results of Ebert and Curry [55] for the four pseudo-spectral intervals (.25-.69 $\mu$m, .69-1.19 $\mu$m, 1.19-2.38 $\mu$m, and 2.38-4.00 $\mu$m) employed in the CAM 3.0 shortwave radiation model. Note that when $0<f_{ice}<1$, then the combination of these expressions in (4.204 - 4.207) represent the radiative properties for a mixed phase cloud.

4.8.5 Cloud vertical overlap

The treatment of cloud vertical overlap follows Collins [38]. The overlap parameterization is designed to reproduce calculations based upon the independent column approximation (ICA). The differences between the results from the new parameterization and ICA are governed by a set of parameters in the shortwave code (Table 4.1 on page and section 4.8.9). The differences can be made arbitrarily small with appropriate settings of these parameters. The current parameter settings represent a compromise between computational cost and accuracy.

The new parameterizations can treat random, maximum, or an arbitrary combination of maximum and random overlap between clouds. The type of overlap is specified with the same two variables for the longwave and shortwave calculations. These variables are the number of random-overlap interfaces between adjacent groups of maximally-overlapped layers and a vector of the pressures at each of the interfaces. The specification of the overlap is completely separated from the radiative calculations, and if necessary the type of overlap can change at each grid cell or time step.

4.8.5.1 Conversion of cloud amounts to binary cloud profiles

The algorithm for cloud overlap first converts the vertical profile of partial cloudiness into an equivalent collection of binary cloud configurations. Let $\cfrac (i)$ be the fractional amount of cloud in layer $i$ in a profile with $K$ layers. The index $i = 1$ corresponds to the top of the model atmosphere and $i = K$ corresponds to the layer adjacent to the surface. Let $N_m$ be the number of maximally-overlapped regions in the column separated by random-overlap boundaries. If the entire column is maximally overlapped, then $N_m = 1$, and if the entire column is randomly overlapped, then $N_m = K$. Each region $j$ includes all layers $i$ between $i_{{j},\min}$ and $i_{{j},\max}$. Within each region, identify the $n_j$ unique, non-zero cloud amounts and sort them into a descending list $\cfrac _{{j},{k_{j}}{}}$ with $1 \le k_j \le n_j$. Note than in CAM 3.0, cloud amounts are not allowed to be identically equal to 1. It is convenient to define $\cfrac _{j,0} = 1$ and $\cfrac _{j,n_j+1} = 0$. By construction $\cfrac _{{j},{k_{j}}{-1}}> \cfrac _{{j},{k_{j}}{}}$ for $1 \le k_j \le n_j+1$.

The binary cloud configurations are defined in terms of the sorted cloud amounts. The number of unique cloud binary configurations in region $j$ is $n_j+1$. The $k_j^{\hbox{th}}$ binary cloud configuration $\tilde{\cfrac }_{{j},k_{j}}$ in region $j$ is given by

$$\tilde{\cfrac }_{{j},k_{j}}(i) = \left\{ \begin{array}{ll} 1 & \m... ...(i) \ge \cfrac _{{j},{k_{j}}{-1}}$} \ 0 & \mbox{otherwise} \end{array} \right.$$ (4.198)

with $1 \le k_j \le n_j+1$. The fractional area of this configuration is

$$\tilde{A}_{{j},k_{j}}= \cfrac _{{j},{k_{j}}{-1}}- \cfrac _{{j},{k_{j}}{}}$$ (4.199)

The binary cloud configurations for each maximum-overlap region can be combined into cloud configurations for the entire column. Because of the random overlap boundaries between regions, the number of column configurations is

$$N_c = \prod_{j' = 1}^{N_m} (n_{j'}+1)$$ (4.200)

Let $\tilde{\cfrac }[k_1,\ldots,k_{N_m}]$ represent the column configuration with $\tilde{\cfrac }_{1,k_1}$ in region 1, $\tilde{\cfrac }_{2,k_2}$ in region 2, etc. The vertical profile of binary cloud elements is given by:

$$\tilde{\cfrac }[k_1,\ldots,k_{N_m}](i) = \sum_{j' = 1}^{N_m} \tilde{\cfrac }_{{j'},k_{j'}}(i)$$ (4.201)

The area of this configuration is

$$\tilde{A}[k_1,\ldots,k_{N_m}] = \prod_{j' = 1}^{N_m} \tilde{A}_{{j'},k_{j'}}$$ (4.202)

4.8.5.2 Maximum-random overlap assumption

The cloud overlap for radiative calculations in CAM 3.0 is maximum-random (M/R). Clouds in adjacent layers are maximally overlapped, and groups of clouds separated by one or more clear layers are randomly overlapped. The two overlap parameters input to the radiative calculations are the number of random-overlap interfaces, which equals $N_m$, and a vector of pressures $\vec p$ at each random-overlap interface. These parameters are determined for each grid cell at each radiation time step. Suppose there are $M \ge 0$ groups of vertically contiguous clouds in a given grid cell. The first parameter $N_m = \max(M,1)$. Let $p_j$ represent the pressure at the bottom interface of each group of contiguous clouds, and let $p_s$ denote the surface pressure. Both $j$ and $p_j$ increase from the top of the model downward. Then

$$\vec p = \left\{\begin{array}{ll} \left[ p_s \right] & \hbox{if $... ...p_1, \>\ldots\>, p_{M-1}, p_s \right] & \hbox{if $M \ge 2$} \end{array} \right.$$ (4.203)

4.8.5.3 Low, medium and high cloud overlap assumptions (diagnostics)

For diagnostic purposes, the CAM 3.0 calculates three levels of cloud fraction assuming the same maximum-random overlap as in the radiative calculations. These diagnostics, denoted as low, middle, and high cloud, are bounded by the pressure levels $p_s$ to 700 mb, 700 mb to 400 mb, and 400 mb to the model top.

4.8.5.4 Computation of fluxes and heating rates with overlap

The solution for the shortwave fluxes is calculated by determining all possible arrangements of binary clouds which are consistent with the vertical profile of partial cloudiness, the overlap assumption, and the parameters for accelerating the solution (Table 4.1 and section 4.8.9). The shortwave radiation within each of these configurations is calculated using the same $\delta$-Eddington solver introduced in CCM3 [27]. The all-sky fluxes and heating rates for the original profile of partial cloudiness are calculated as weighted sums of the corresponding quantities from each configuration. The weights are equal to the horizontal fractional area occupied by each configuration. The number of configurations is given by eqn. (4.200), and the area of each configuration is given by eqn. (4.202). There are two steps in the calculations: first, the calculation of the cloud-free and overcast radiative properties for each layer, and second the combination of these properties using the adding method to calculate fluxes. These two processes are described below.

4.8.6 $\delta$-Eddington solution for a single layer

Details of the implementation are as follows. The CAM 3.0 model atmosphere is divided into $K+1$ layers in the vertical; an extra top layer (with index 0, above the $K$ layers specified by CAM 3.0) is added. This extra layer prevents excessive heating in the top layer when the top pressure is not very low; also, as the model does not specify absorber properties above its top layer, the optical properties of the top layer must be used for the extra layer. In CAM 3.0, clear-sky and all-sky solar fluxes are calculated and output for the top of model (TOM) at layer 1 and the top of atmosphere (TOA) corresponding to layer 0. The TOM fluxes are used to compute the model energetic balance, and the TOA fluxes are output for diagnostic comparison against satellite measurements. The provision of both sets of fluxes is new in CAM 3.0. Layers are assumed to be horizontally and vertically homogeneous for each model grid point and are bounded vertically by layer interfaces. For each spectral band, upward and downward fluxes are computed on the layer interfaces (which include the surface and top interface). The spectral fluxes are summed and differenced across layers to evaluate the solar heating rate. The following discussion refers to each of the spectral intervals.

In general, several constituents absorb and/or scatter in each homogeneous layer (e.g. cloud, aerosol, gases...). Every constituent is defined in terms of a layer extinction optical depth $\tau$, single scattering albedo $\omega$, asymmetry parameter $g$, and the forward scattering fraction $f$. To define bulk layer properties, the combination formulas of Cess [30] are used:

$$\tau = \sum_i \tau_{i} ,$$ (4.204)

$$\omega = {\sum_i \omega_{i} \tau_{i}}{\tau} ,$$ (4.205)

$$g = \frac{\sum_i g_{i} \omega_{i} \tau_{i}}{\omega\tau} ,$$ (4.206)

$$f = \frac{\sum_i f_{i} \omega_{i} \tau_{i}}{\omega \tau} ,$$ (4.207)

where the sums are over all constituents.

The $\delta$-Eddington solution for each layer requires scaled properties for $\tau$, $\omega$, $g$, given by the expressions:

$$\tau^\ast = \tau (1-\omega f) ,$$ (4.208)

$$\omega^\ast = \omega \left(\frac{1 - f}{1 - \omega f}\right) ,$$ (4.209)

$$g^\ast = \frac{g - f}{1 - f} .$$ (4.210)

The scaling accounts for the scattering effects of the strong forward peak in particle scattering. The $\delta$-Eddington nonconservative ($\omega<1$) solutions for each layer for direct radiation at cosine zenith angle $\mu_0$ are (following the notation of Coakley et al. [35]:

$$R(\mu_0) = (\alpha - \gamma) {\overline T}e^{{-\tau^\ast/\mu_0}} + (\alpha + \gamma) {\overline R} - (\alpha-\gamma) ,$$ (4.211)

$$T(\mu_0) = (\alpha + \gamma) {\overline T} + (\alpha - \gamma){\overline R}e^{{-\tau^\ast/\mu_0}} - (\alpha + \gamma - 1) e^{{-\tau^\ast/\mu_0}} ,$$ (4.212)

$${\overline R} = (u+1)(u-1)(e^{\lambda\tau^\ast} - e^{-\lambda\tau^\ast}) N^{-1} ,$$ (4.213)

$${\overline T} = 4u N^{-1} ,$$ (4.214)
where

$$\alpha = \frac{3}{4} \omega^\ast \mu_0 \left( \frac{1+g^\ast(1-\omega^\ast)} {1 - \lambda^2\mu_0^2} \right) ,$$ (4.215)

$$\gamma = \frac{1}{2} \omega^\ast \left(\frac{1+3g^\ast(1-\omega^\ast)\mu_0^2} {1 - \lambda^2\mu_0^2} \right) ,$$ (4.216)

$$N = (u+1)^2 e^{\lambda\tau^\ast} - (u-1)^2 e^{-\lambda\tau^\ast} ,$$ (4.217)

$$u = \frac{3}{2} \left({1-\omega^\ast g^\ast}{\lambda}\right),$$ (4.218)

$$\lambda = \sqrt{3(1-\omega^\ast)(1-\omega^\ast g^\ast)} ,$$ (4.219)

where $R(\mu_0)$, $T(\mu_0)$ are the layer reflectivity and transmissivity to direct radiation respectively, and $\overline R$, $\overline T$ are the layer reflectivity and transmissivity to diffuse radiation respectively. It should be noted that in some cases of small but nonzero $\omega$, the diffuse reflectivity can be negative. For these cases, $\overline R$ is set to 0, which produces negligible impact on fluxes and the heating rate. Note that in the new overlap scheme, these properties are computed separately for the clear and cloud-filled portions of each layer [38].

4.8.7 Combination of layers

To combine layers, it is assumed that radiation, once scattered, is diffuse and isotropic (including from the surface). For an arbitrary layer 1 (or combination of layers with radiative properties $R_1(\mu_0)$, $T_1(\mu_0)$, $\overline{R}_1$, $\overline{T}_1$) overlaying layer 2 (or combination of layers with radiative properties $R_2(\mu_0)$, $T_2(\mu_0)$, and $\overline{R}_2,\overline{T}_2$), the combination formulas for direct and diffuse radiation incident from above are:

$$R_{12}(\mu_0) = R_1(\mu_0) + \frac{\overline{T}_1 \{ (T_1(\mu_0) - e^{-\tau_1^\as... ..._2 + e^{-\tau_1^\ast/\mu_0}R_2(\mu_0) \}} {1 - \overline{R}_1 \overline{R}_2} ,$$ (4.220)

$$T_{12}(\mu_0) = e^{-\tau_1^\ast/\mu_0} T_2(\mu_0) + \frac{\overline{T}_2 \{ (T_1(... ...1^\ast/\mu_0}R_2(\mu_0) \overline{R}_1 \}}{1 - \overline{R}_1 \overline{R}_2} ,$$ (4.221)

$$\overline{R}_{12} = \overline{R}_1 + \frac{\overline{T}_1 \overline{R}_2 \overline{T}_1} {1 - \overline{R}_1 \overline{R}_2} ,$$ (4.222)

$$\overline{T}_{12} = \frac{\overline{T}_1 \overline{T}_2} {1 - \overline{R}_1 \overline{R}_2} .$$ (4.223)

Note that the transmissions for each layer ( $T_1(\mu_0),T_2(\mu_0)$) and for the combined layers $(T_{12}(\mu_0))$ are total transmissions, containing both direct and diffuse transmission. Note also that the two layers (or combination of layers), once combined, are no longer a homogeneous system.

To combine the layers over the entire column, two passes are made through the layers, one starting from the top and proceeding downward, the other starting from the surface and proceeding upward. The result is that for every interface, the following combined reflectivities and transmissivities are available:

$$e^{-\tau^\ast/\mu_0} = \;$$ direct beam transmission from top-of-atmosphere to the

$$\tau^\ast$$ to the interface),

$$R_{up}(\mu_0) = \;$$ reflectivity to direct solar radiation of entire atmosphere below the interface,

$$T_{dn}(\mu_0) = \;$$ total transmission to direct solar radiation incident from above to entire atmosphere above the interface,

$${\overline R}_{up} = \;$$ reflectivity of atmosphere below the interface to diffuse radiation from above,

$${\overline R}_{dn} = \;$$ reflectivity of atmosphere above the interface to diffuse radiation from below.

With these quantities, the upward and downward fluxes at every interface can be computed. For example, the upward flux would be the directly transmitted flux ( $e^{-\tau\ast/\mu_0}$) times the reflection of the entire column below the interface to direct radiation ( $R_{up}(\mu_0)$), plus the diffusely transmitted radiation from above that reaches the interface ( $T_{dn}(\mu_0) - e^{-\tau\ast/\mu_0}$) times the reflectivity of the entire atmosphere below the interface to diffuse radiation from above ( $\overline{R}_{up}$), all times a factor that accounts for multiple reflections at the interface. A similar derivation of the downward flux is straightforward. The resulting expressions for the upward and downward flux are:

$$F_{up} = \frac{e^{-\tau^\ast/\mu_0} R_{up}(\mu_0) + (T_{dn}(\mu_0) - e^{-\tau^\ast/\mu_0})\overline{R}_{up}} {1 - \overline{R}_{dn}\overline{R}_{up}},$$ (4.224)

$$F_{dn} = e^{-\tau^\ast/\mu_0} + \frac{(T_{dn}(\mu_0) - e^{-\tau^\ast/\mu... ...\mu_0}R_{up}(\mu_0)\overline{R}_{dn}} {1 - \overline{R}_{dn}\overline{R}_{up}}.$$ (4.225)

Note that in the new overlap scheme, the calculation of the combined reflectivities, transmissions, and fluxes at layer interfaces are computed for each binary cloud configuration, subject to techniques for significantly accelerating these calculations (below) [38].

4.8.8 Acceleration of the adding method in all-sky calculations

If two or more configurations of binary clouds are identical between TOA and a particular interface, then $T_{dir}= e^{-\tau^*/\mu_0}$, $T_{dn}$, and $\bar{R}_{dn}$ are also identical at that interface. The adding method is applied once and the three radiative quantities are copied to all the identical configurations. This process is applied at each interface by constructing a binary tree of identical cloud configurations starting at TOA down to the surface. A similar method is used for $R_{up}$ and $\bar{R}_{up}$, which are calculated using the adding method starting the surface and continuing up to a particular interface. The copying of identical radiative properties reduces the number of calculations of $T_{dir}$, $T_{dn}$, and $\bar{R}_{dn}$ by 62% and the number of calculations of $R_{up}$ and $\bar{R}_{up}$ by 21% in CAM 3.0 with M/R overlap.

4.8.9 Methods for reducing the number of binary cloud configurations

The computational cost of the shortwave code has two components: a fixed cost for computing the radiative properties of each layer under clear and overcast conditions, and a variable cost for applying the adding method for each column configuration $\tilde{\cfrac }[k_1,\ldots,k_{N_m}]$. The variable component can be reduced by omitting configurations which contribute small terms in the shortwave fluxes. Several mechanisms for selecting configurations for omission have been included in the parameterization. The parameters that govern the selection process are described in Table 4.1.

Table 4.1: Parameters for Decreasing Number of SW Calculations.

Parameter	Symbol	Definition	Value in CAM 3.0
cldmin	$\cfrac _{\min}$	Minimum cloud area	0	.
cldeps	$\cfrac _{eps}$	Minimum cloud area difference	0	.
areamin	$\tilde{A}_{\min}$	Minimum configuration area	0	.	01
nconfgmax	$N_{c,\max}$	Maximum # of configurations	15	.

Any combination of the selection conditions may be imposed. If the parameter $\cfrac _{\min}> 0$, cloud layers with $\cfrac (i) \le \cfrac _{\min}$ are identified as cloud-free layers. The configurations including these clouds are excluded from the flux calculations. If the parameter $\cfrac _{eps}> 0$, the cloud amounts are discretized by

$$\cfrac (i) \rightarrow \left[{\cfrac (i) \over \cfrac _{eps}}\right] \cfrac _{eps}$$ (4.226)

where $[x]$ represents rounding to the nearest integer less than $x$. This reduces the number of unique, non-zero cloud amounts $n_j$ in each maximum-overlap region $j$. For example, if $\cfrac _{eps}= 0.01$, then two cloud amounts are distinguished only if they differ by more than 0.01. If the parameter $\tilde{A}_{\min}> 0$, only configurations with $\tilde{A}[k_1,\ldots,k_{N_m}] \ge \tilde{A}_{\min}$ are retained in the calculation. The fluxes and heating rates are normalized by the area of these configurations:

$$\tilde{A}_{tot}= \sum_{k_1 = 1}^{n_1 + 1} \cdots \sum_{k_{N_m} = ... ..._{N_m}] \;\;\theta\left(\tilde{A}[k_1,\ldots,k_{N_m}] - \tilde{A}_{\min}\right)$$ (4.227)

where $\theta$ is the Heaviside function. In CAM 3.0, $\tilde{A}_{\min}= 0.01$. Finally, if the number of configurations $N_c > N_{c,\max}$, then only the $N_{c,\max}$ configurations with the largest values of $\tilde{A}[k_1,\ldots,k_{N_m}]$ are retained. This is equivalent to setting $\tilde{A}_{\min}$ so that the largest $N_{c,\max}$ configurations are selected. The fluxes and heating rates are normalized by $\tilde{A}_{tot}$ calculated with this value of $\tilde{A}_{\min}$. With the current cloud parameterizations in CAM 3.0 and with $\tilde{A}_{\min}= 0.01$, the mean and RMS $N_c$ are approximately 5. $N_{c,\max}$ is set to 15, or 2 standard deviations above the mean $N_c$. Only 5% of cloud configurations in CAM 3.0 have $N_c \ge N_{c,\max}$. The errors of the solutions relative to ICA are relatively insensitive to $\tilde{A}_{tot}$ [38].

4.8.10 Computation of shortwave fluxes and heating rates

The upward and downward spectral fluxes at each interface are summed to evaluate the spectrally integrated fluxes, then differenced to produce the solar heating rate,

$$Q_{\rm sol} = \frac{g}{c_p} \frac{F_{dn}(p_{k+1}) - F_{up} (p_{k+1}) - F_{dn}(p_k) + F_{up}(p_k)} {p_{k+1} - p_k}$$ (4.228)

which is added to the nonlinear term $(Q)$ in the thermodynamic equation.