The parameterization of radiative processes has been prepared for NMC by Stephen B. Fels and M. Daniel Schwarzkopf of GFDL and their generous contribution of time and effort is gratefully acknowledged. This section documents the radiation scheme as implemented at NMC in August 1987. Modeling of both shortwave and longwave processes include the effects of the major radiatively active atmospheric constituents - water vapor (H2O), carbon dioxide (CO2), ozone (O3), and clouds. Outputs from the parameterization are forecast model layer radiative heating rates (_K/sec), as well as downward longwave and net shortwave fluxes (cal cm-2 min-1) at the earth's surface.

The radiation calculations are computationally quite expensive. Using compiler optimized CYBER 205 coding, the average radiation computation time for one 18 layer column of atmospheric data is .0031 seconds. This requires approximately 40 seconds of computation time for a rhomboidal-40 grid array. The faster shortwave calculations account for only 18% of this time. Because of their computational expense, radiation calculations are made only once every 12 hours. Longwave heating rates are held fixed for the entire 12 hours; however, provision is made within the forecast model for an approximate diurnal cycle, so shortwave heating rates and surface radiative fluxes are allowed to change during this interval.

The first part of this report describes the necessary alteration of model temperature and moisture data for the radiation parameterization. Then, climatological surface albedo, ozone concentration, and cloudiness are discussed. Finally, basic techniques which are used to model both shortwave and longwave processes are described. A schematic of the radiation calculations, including basic subroutine names, is shown in Fig 3.1.

  1. Preprocessing Input Data
    1. Atmospheric Temperature
Layer virtual temperatures, Tv, which are received from the forecast model are converted to thermodynamic temperatures, T, as in eqn (3.1):

where q is layer specific humidity. This calculation is made only for the forecast model's lower 12 moisture bearing layers. The extrapolated moisture in the upper model layers (described below) is small and is neglected in eqn (3.1).

  1. Moisture Extrapolation
Specific humidities for the dry upper layers of the forecast model are obtained for radiation computations by extrapolation. First, q in moisture bearing layers of the model are adjusted in order to remove negative values or supersaturations. The constraint is that layer relative humidities must lie between .15 and 1.0. Then extrapolation of q into the dry layers is made using a power law which reasonably describes the mean vertical moisture distribution for the entire atmosphere:

where qk,pk are dry layer values and qm,pm are uppermost moist layer (currently layer 12) values of specific humidity and pressure. The exponent, l, is computed at each model grid point by inverting eqn (3.2) , assuming qk=3.x10-6 g/g at pk=5 cb. Moisture is assumed constant for pressures less than 5 cb. Relative humidity is constrained to 100% to avoid excessive extrapolated q in tropical regions.

  1. Initial Land Surface Temperature
There is no available analysis of land surface temperature. There are two options:

(1) Assuming an initial value of 290K at all land points, compute the radiative fluxes and, then, solve a surface energy balance equation for surface temperature. With these new land temperatures, make another pass through the radiation scheme.

(2) Extrapolate (linear in pressure) temperature down to the earth's surface from the lowest two forecast model layers.

There are several problems with the first option. The surface energy balance algorithm is not the same as being used in the forecast model and produces too large initial surface temperatures in regions where the sun is "high in the sky". Additionally, the process of solving for an energy balanced temperature and recomputing radiative heating rates should be iterated more than once (at least 2 or 3 iterations seem necessary). Shortcomings of both methods can be neglected in most regions, because the forecast model's "memory" of the initial surface temperature is quite short once the sun rises.

  1. Climatological Input Data
    1. Surface Albedo
Monthly values of background land surface albedo, obtained from Matthews (l985), are available on the forecast model's computational grid. Because these arrays change very slowly during the year, there is no interpolation from monthly data to the day of the forecast. The albedo is adjusted for snow depth (see section 3.3.5) prior to the computation of shortwave heating rates.

  1. Ozone
Seasonal zonal mean ozone amounts are available in block data (BLCKFS) form for all model layers. The 5_ latitude climatology (Fig 3.2) has been constructed from Hering and Borden (l965) and London (l962) data. Zonal mean data is horizontally interpolated to forecast model latitudes, and it changes for each radiation computation because the seasonal values are interpolated, in time, to the specific hour of the forecast.

  1. Clouds
Three stratiform cloud types (high, middle, low) are specified (Telegadas and London, l954). Seasonal zonal mean cloud amount (Fig 3.3) and cloud top/bottom (Fig 3.4) are available at 5_ latitude intervals in block data (BLCKFS) form. Cloud amount is not horizontally interpolated to forecast model latitudes, rather the nearest 5_ value is used. However, it changes for each radiation computation because there is a time-interpolation from seasonal data to the specific hour of the forecast. Cloud top and bottom are also obtained from the nearest 5_ latitude data, but they are not time-interpolated from seasonal data to the hour of the forecast.

  1. Shortwave Processes
Incoming solar radiation is the ultimate driving force for atmospheric circulations. The major atmospheric gases which absorb shortwave solar radiation are ozone, water vapor, and carbon dioxide. Ozone absorption is strong in ultraviolet wavelengths and weaker in the visible wavelengths, whereas water vapor absorbs in the visible and near infrared wavelengths. Carbon dioxide absorption in the near infrared region of the solar spectrum is the smallest of the three gases. Multiplication of transmission functions for the gases accounts for the overlapping of ozone and water vapor absorption bands and the overlap of water vapor and carbon dioxide bands.

In order to model the wavelength dependent absorption coefficients for water vapor across the solar spectrum, a probability distribution, p(k), is assumed for the absorption coefficients in a manner virtually identical to the work of Lacis and Hansen (l974). As in eqns (28) and (29) of Lacis and Hansen the absorption, A, by water vapor is:

where is water vapor amount, k is the absorption coefficient for water vapor, p(k) dk is the fraction of incident solar flux associated with an absorption coefficient between k and k+dk, p(kn) is the discrete probability distribution for N subintervals over the solar spectrum, and kn, is the mean absorption coefficient for each subinterval. In order to overlap ozone absorption with water vapor processes in the ultraviolet part of the solar spectrum, Lacis and Hansen's first subinterval is split in two, giving N=9 (table 3.1). Thus O3 and H2O overlap is assigned to the n=1 subinterval, and CO2 and H2O overlap is assigned to n=2,3,...,9. Since O3 and CO2 are each relegated to intervals covering 1/2 the solar spectrum, their respective absorptions will be weighted by 1/2 (the p(kn)'s), so a compensating multiplication by 2 is necessary.

Absorber amounts (cm) for each of the three gases are computed as total amounts from top of the atmosphere down to each model level (using model layer values of H2O,O3,CO2), plus, for reflected paths from the earth's surface, additional amounts up to each model level. Thus for a model with j levels, there are (2j-1) total amounts of each absorbing gas. For the downward radiation path, total absorber amount, xl, between the top of the atmosphere and the lth model level, traversed by the direct solar beam and by diffuse solar radiation below the top of the highest cloud (if any), lc, is:

where ak is amount of an absorber in the kth model layer, M is a magnification factor accounting for slant path and refraction, and m is an effective magnification (diffusivity) factor for diffuse radiation. For the three absorbing gases, M and m are defined as:

where is the cosine of the solar zenith angle, M is from Rodgers (l967), and m is from Lacis and Hansen (l974).

For the upward returning path, where ls is the level of the earth's surface, total absorber amount is:

Transmission functions for each of the three gases are computed from the top of the atmosphere to a particular model level. Shortwave flux at each level is obtained by weighting solar flux at the top of the atmosphere by p(kn) and multiplying it by the transmission functions. The solar flux at the top of the atmosphere is , where S=1367.4 Watts/m2, and , is the cosine of the solar zenith angle. Net fluxes (upward-downward) are computed at each model level and the layer heating rates, ,are computed from the vertical component of radiative flux divergence:

where F and p are level values of net radiative flux and pressure respectively.

  1. Ozone
The ozone absorption of solar radiation by the weak visible (Chappius) band and the strong ultraviolet (Hartley and Huggins) bands are estimated by analytic expressions of precise frequency integrated absorption curves taken from Lacis and Hansen (see their eqns 8, 9, 10 and Fig. 6). Ozone amount, , is calculated from top of the atmosphere down to each model level as in eqns (3.4) - (3.6). Defining A, as the fraction of absorbed incident solar flux, the absorption in the weak visible band is:

For the ultraviolet bands:

Total ozone absorption ,, is . Ozone transmission down to each model level, , is computed and later multiplied by H2O transmission in the n=1 subinterval (table 3.1) to account for absorption band overlap.

Ozone absorption occurs primarily in the high atmospheric layers, where molecular concentration is low enough to neglect scattering. Most of the scattering occurs in the lower atmosphere where there is little ozone absorption. Thus Lacis and Hansen model the scattering process by assuming a simple 2 layer atmosphere; a purely absorbing layer on top of a purely reflecting layer. The albedo of the earth surface region, , is a combination of the reflectivity of the earth's surface, Rg, and the effective albedo of the non-absorbing lower atmosphere, , while accounting for multiple reflections between ground and the atmosphere:

where is the average of over all solar zenith angles. Simple expressions for the various reflectivities are obtained, in a least squares sense, for clear skies by Lacis and Hansen (their cloudy sky formulation is not used):

  1. Water Vapor
The parameterization of water vapor absorption is more complicated because there are many absorbing bands as well as significant pressure dependence of the absorption coefficients. Scattering and absorption can occur in the same atmospheric layer, thus precluding use of the simple ozone reflecting model noted in the preceding section. Water vapor amount, , is calculated from top of the atmosphere down to each model level as in eqns (3.4) - (3.6). A scaling approximation is used as a correction for the pressure dependence of the absorption coefficients, so the layer absorber amount is pressure weighted by pk/psfc where pk is layer pressure and psfc is surface pressure. Transmission functions, Tr, for each of the 9 subintervals in the solar spectrum are computed down to each model level:

where kn is the mean absorption coefficient for each subinterval in Table 3.1.

  1. Carbon Dioxide
Parameterization of CO2 absorption follows the work of Sasamori, et al. (l972). CO2 mixing ratio is 330 ppmv and, similar to H2O, layer absorber amounts are scaled by pk/psfc to account for the pressure dependence of absorption coefficients along the optical path. Total carbon dioxide amount,, from the top of the atmosphere to each model level is calculated as in eqns (3.4) - (3.6). Absorption, A, is :

Transmission functions to each model level, , are computed and later multiplied by H2O transmission functions in subintervals 2 through 9 to account for absorption band overlap.

  1. Clouds
The previous sections discussed the calculation of clear sky transmission functions from top of the atmosphere to each level, however, the occurrence of cloud layers complicates the computations. Clear sky transmission is weighted by cloud effects, which are functions of cloud fraction, cloud albedo, and cloud (water) absorption. High, middle, and low clouds are assumed randomly located with respect to each other within each grid box. Universal cloud radiative properties for each cloud type are listed in Table 3.2. No provision is made to vary cloud albedo or absorption with (estimated) cloud water content or solar zenith angle. There is no cloud absorption in the ultraviolet part of the solar spectrum (1st subinterval of Table 3.1).

For convenience, the top of the atmosphere and earth's surface are considered cloud layers. If one makes the additional assumption that the reflectivity of both top and bottom of the 3 physical clouds are equal, then 5 values of cloud reflectivity and transmissivity define the entire cloud system. Reflectivity, (CR)i, and transmissivity, (CT)i, weighted for fractional cloudiness are defined:

where Ci is fractional cloud amount, (alb)i and (abs)i are cloud albedo and absorptivity (Table 3.2), and i is cloud type. At the earth's surface (i=5), (CR)i is calculated as in eqn (3.14); with C5=1., (alb)5=Rg for the second through ninth subinterval of the solar spectrum, and (alb)5= for the first subinterval (see eqn (3.10)).

The calculation of shortwave fluxes in the presence of clouds entails a computation of both upward and downward radiative fluxes at, potentially, 8 reflecting surfaces (defined above). There is provision for multiple reflections, and, of course, the downward flux at the top of the atmosphere is known. Consider two clouds of index i and i+1 as depicted in Fig 3.5. Assuming the downward flux at the top of ith cloud ( ) is known, the following fluxes are defined:

where ti and bi represent model levels of top and bottom for cloud index i; and represents clear sky transmission function from level ti to level bi. Manipulation of eqns (3.15-3.18) yields equations for upward flux at the top of cloud index i and downward flux at the top of the next lower cloud:

Now assume that the transmissivities between two surfaces are independent of direction (i.e. ) and that (i.e. the transmissions are exponential). Both assumptions are exactly true for bands 2-9 (disregarding the CO2 transmissivity) and are assumed true in the O3 band 1. Bearing this in mind, define

where (TCLU)i represents the transmission between tops of adjacent reflecting surfaces and (TCLD)i is clear sky transmission between the bottom of cloud index i and the top of the next lower reflecting surface, i+1 (Fig 3.5). Using the above definitions as well as the assumption about the exponential nature of the transmission functions, eqns ( 3.19, 3.20) may be rewritten as:

Recalling that is known, and that for the earth's surface (i=5) the relation between and is

one obtains (after suitably manipulating eqns (3.23-3.25) for i=4):

Simplifying (3.26):

The quantity (ALFA4+(CR)4) thus represents the effective albedo of the low cloud-ground system accounting for multiple reflections from below the cloud. Perhaps this can be visualized by looking at the low cloud/earth surface region below level t4 in Fig 3.5. For one reflection at the ground

With an additional reflection off the bottom of the low cloud

where (CR)4 (CR)5 (TCLD4)2 accounts for clear sky transmission to the base of the low cloud, a reflection off it, then transmission back and a subsequent reflection off the earth's surface. With "n" reflections

Since (CR)4 (CR)5 (TCLD4)2 < l, the series in parentheses converges in the limit as n approaches infinity and can be written as:

The term in parentheses is a multiple reflection factor and is seen in eqn (3.26). One may now proceed to solve eqns (3.23,3.24) for the next surface, i=3, using the lower boundary condition (3.27) in place of eqn (3.25). The resulting solution for is the same form as (3.26), with all indices decremented by l and (CR)5 replaced by (ALFA4+(CR)4) :

For the nth cloud, the ALFA term may be generalized (n=l-4):

where ALFA5=0.

The above procedure is repeated until the upward flux has been calculated for the top of the highest physical cloud (i=2). Since is a clear sky flux, and is already known, eqns 3.30, 3.23, 3.24, 3.15-3.18, may be employed to obtain upward and downward fluxes for all reflecting surfaces. Fluxes for a model level between the clouds are easily obtained using clear-sky transmittances between the appropriate cloud top and bottom and the desired level. Finally, for cases where the cloud is more than one model layer thick (currently low cloud), linear interpolation of transmission functions between cloud top and bottom produces values at model levels inside the cloud. Application of eqn (3.7) results in identical shortwave heating rates for these internal model layers.

  1. Surface Albedo Adjustment
Surface albedo is adjusted for the presence of snow at all land and sea-ice points. Poleward of approximately 70_ latitude, surface albedo for snow covered land or sea ice is set to 0.75. Elsewhere the albedo, Rg is a function of the background surface albedo, A0 (A0 = 0.5 for sea ice), and the snow depth, d (cm), as in eqn (3.31):

A zenith angle dependent surface albedo over open water is obtained from Payne's (l972) tabulated data. Payne accounts for the effect of atmospheric scattering and absorption on the surface albedo by tabulating the data as a function of atmospheric transmittance, TRANS, as well as zenith angle. The radiation parameterization uses a value of TRANS = .537 for all grid points, where TRANS=1. represents clear atmosphere.

  1. Approximate Diurnal Cycle
Since the sun governs the diurnal cycle, shortwave (SW) processes are the most important to consider. Letting be the cosine of solar zenith angle at a point on the earth's surface, the "exact" shortwave flux, Fsw , at a level in the atmosphere can be represented as:

where S is the solar constant, t is time, and Tr is the atmospheric transmissivity. Tr is a function of H2O, CO2, O3, other trace gases and aerosols (not considered here), and clouds. Holding these atmospheric profiles fixed (), the SW flux integrated over a complete day is:

where the final equality exists because at night.

Ideally, one could simulate the diurnal cycle at a grid point by making radiation calculations each model time step, thereby accounting for changes in zenith angle, H2O, temperature (for longwave computations) and clouds. This is computationally prohibitive however, so some approximations must be made. For example, one could either calculate the "less expensive" SW effects more frequently than longwave (LW) processes, or calculate both SW and LW effects more frequently on a coarse subset of the grid and interpolate to the full grid. An alternative approximation, developed by Ellingson (U of Md) and Campana, for SW processes adds little computational overhead (Campana, l986). A mean cosine of the solar zenith angle, , computed for each model latitude as in eqn (3.34) of Manabe and Strickler (1964), and is used to calculate the SW fluxes at each model "grid point".

Shortwave radiative fluxes are computed every 12 hours at all grid points using this mean cosine of the zenith angle, , as well as atmospheric transmissivities valid at those times (eqn (3.36)).

For the nondiurnal option these fluxes are fixed, but for the diurnal cycle option, grid point values of SW fluxes (and thus heating rates) are weighted by the exact cosine of the zenith angle, , at each forecast model time step, as in eqn (3.35):


and whererepresents grid point transmissivity at the beginning of a 12-hour interval. Though changes in H2O will occur during the 12-hour interval, a (reasonable) assumption is made that changes in are much more important than changes in water vapor profiles in determining SW fluxes. Integration of this SW flux over 24 hours, again holding fixed, yields the same result as eqn (3.33)

  1. Longwave Processes
Longwave radiative processes occur in a region of the electromagnetic spectrum distinct from shortwave influences; at wavelengths greater than 4 (=10-4 cm). In general, longwave effects are more difficult to model than those in the shortwave part of the spectrum since both absorption and emission occur along optical paths and effects from all atmospheric layers must be included. The parameterization scheme, which has been developed by Fels and Schwarzkopf at GFDL, includes radiative effects of the ozone 9.6 band, the carbon dioxide 15 band, the water vapor band at 6.3 , the H2O "rotational " bands at wavelengths greater than 12, and weak continuous H2O absorption in the 8- 18range. The 8- 12region contains virtually no specific H2O absorption bands and the weak continuous absorption is thought to be due, in part, to effects from strong H2O absorption outside this "window". Radiative effects from this "continuum" region occur primarily in the lower moist tropical atmosphere.

As in the parameterization of shortwave processes, the desired output from the longwave scheme are net radiative fluxes at all model levels which then are used to compute model layer heating rates (eqn (3.7)). Appropriate radiative transfer equations for the calculation of upward (F" ) and downward (F# ) longwave flux at a specified atmospheric level, z, are (see Stephens, l984, p. 828):

where n is wavenumber (i. e., wavelength-1), Bn is the Planck function, zT is top of the atmosphere, and is the longwave transmission function between level and z . is defined over all zenith angles along an optical path, du, in eqn (3.40):

where is cosine zenith angle and is the absorption coefficient (for an absorber concentration, u) which varies with pressure, p, and temperature, T, along the path from to z. The four integrals imbedded in eqns (3.38) - (3.40) make the numerical calculation of longwave effects quite expensive.

One can rewrite the radiative flux equations in a more convenient form by integrating eqns (3.38) - (3.39) by parts to obtain:

The net flux at level z, using eqns (3.41) - (3.42) is:

Changing vertical coordinate from z to pressure, p, noting that p=0 at z=zT and p=psfc at z=0, the net longwave flux equation (3.44) is virtually identical to eqn (1) in Fels and Schwarzkopf (l975):

Computation of net flux at each model level requires the evaluation of integrals in eqn (3.40) and (3.44) for

1. zenith angle (),
2. absorber optical path (),
3. spectral interval (), and
4. transmission from all levels ().

Fels and Schwarzkopf parameterize longwave processes using techniques designed to reduce the computational overhead inherent in these integrals, while, at the same time, retaining accuracy in the calculations. The integral over zenith angle ( ) is removed by scaling the gaseous absorber path (i.e., amount) by a diffusivity factor of 1.66. In this manner diffuse transmission from to p is represented, with little loss of accuracy (Rodgers and Walshaw, l966), as a beam of radiation along a path whose zenith angle,, is approximately 53_.

The integral over optical path (du in eqn (3.40)) is difficult to compute because the absorption coefficient depends upon pressure and temperature, both of which can change rapidly along du. There are techniques which account for these inhomogeneous paths by adjusting pressure and path length to create a homogeneous path having a mean absorption coefficient. Fels and Schwarzkopf use two methods for their water vapor calculations. A simple scaling of the optical path by is used in their emissivity heating rate computations (discussed below), where =1013.25 mb. A more accurate scaling is the 2 parameter Curtis-Godson approximation where both pressure and absorber amount are adjusted to create the homogeneous path. Fels and Schwarzkopf employ this technique with the H2O random band model used for their exact heating rate calculations (discussed below).

The integral over spectral interval ( in eqn. (3.44)) is quite difficult to calculate because the fine scale of the absorption lines must be multiplied by the vertical derivative of the much smoother Planck function, Bn. Typically, methods are used which separate the longwave spectrum into spectral intervals small enough both to consider the Planck function a constant and to produce transmission functions which are exponential in nature (Stephens, l984). Though O3 and CO2 are evaluated as one band, H2O processes are complex enough to require the smaller intervals. Using statistical properties for groups of absorption lines, an absorption coefficient can be defined as a function of the strengths, separations, and positions of the detailed line structure within each spectral interval (or band). In "exact" calculations for water vapor, Fels and Schwarzkopf use the Rodgers and Walshaw (1966) random band model distribution of absorption lines and spectral intervals. In order to simplify the overlapping region of H2O and CO2 absorption (15), the Rodgers and Walshaw bands 6,7,and 8 have been restructured to place all of the CO2 effects in band 7 (see table 3.3, where, is mean line strength, is mean distance between absorption lines, andis (Lorenz) line width).

A less accurate but computationally more efficient alternative to the band model approximation for is the emissivity approximation (see Stephens, l984; Paltridge and Platt, l976, p. 173). This treats the absorption as a constant for all wavelengths within each spectral band and allows the integral over the entire spectrum to be performed once, thus saving computer time. Fels and Schwarzkopf apply the emissivity approximation to eqn (3.44) and obtain the following flux equation (for H2O):

where are pre-calculated tables for temperature between 100K - 370K (every 10K) and for water vapor amount between 10-16 and 102 g/cm2 (180 values, equal spacing in ln u):

where n is the number of spectral intervals for H2O and is the transmission function for each band. The integral in eqn (3.45) is evaluated as

where j, j+1 refer to model levels and j+1/2, i+1/2 refer to model midlayers. For nearby layers, where , a precomputed array, G3, is used to evaluate the integral:

Further details may be found in Fels and Schwarzkopf (l975).

One integral remains to be discussed - the integration over all atmospheric layers ( in eqn (3.44)). Fels and Schwarzkopf simplify this integration by separating the net radiative flux equation (3.44) into two terms (Green, l967):

1. Cooling to space (CTS), which is longwave emission from level p directly to space , and

2. Internal EXCHANGE between atmospheric levels, accounting for absorption and re-emission. Stephens (1984) notes that this process primarily occurs in "nearby" layers, because at greater distances its importance drops with the exponential decrease of the transmission function.

In an isothermal atmosphere , and eqn (3.44) can be written as:

where Bn is the Planck function for each spectral interval, n, in the band model representation of the longwave spectrum. The CTS approximation assumes that (3.48) can be used to calculate flux at all model levels even though the atmosphere is not isothermal. The simplicity of this scheme can be seen in the representation of transmission functions in (3.44) and (3.48); the matrix calculation implied by in (3.44) is replaced by a function in (3.48) which depends only on absorber amount above level p. Of course, the CTS approximation is not accurate enough by itself and the EXCHANGE term must be added.

Fels and Schwarzkopf treat the cooling-to-space and exchange calculations in different manners in order to obtain quite accurate results with a computationally efficient method. Recall from eqn (3.7) that the heating rate, Q, is the vertical component of the radiative flux divergence. Computation of F(p) from eqn (3.44) produces a total heating rate, QTOTAL. Separation of the flux calculation into two terms can be used to obtain a CTS heating rate, QCTS (using eqn 3.48) and an exchange rate, QEXCHANGE (ie. the integral term in eqn 3.44). Thus

Noting that the latter term is simply

Fels and Schwarzkopf find that the EXCHANGE is relatively insensitive to the method of calculation. Thus they compute QEXCHANGE using the computationally more efficient emissivity () method (3.45). The calculation of QCTS (in eqn (3.49) however is done more accurately - via either the Rodgers and Walshaw random band (RB) model for H2O or accurate pre-calculated transmission functions for CO2. The "exact" CTS calculation adds relatively little overhead to the efficient computation for the EXCHANGE term because the former is a vector, not matrix, operation. The approximation to total longwave heating rate for each model layer, using (3.49, 3.50), now becomes:

Heating rate calculations for both H2O and CO2 employ equation (3.51), however O3 processes are calculated exactly, using eqn (3.44), for a one band model.

  1. Water Vapor
For the "exact" CTS calculations (i.e., in eqn 3.51), the Rodgers and Walshaw random band model (Table 3.3) is used. The transmission function, in eqn (3.48), is computed for each spectral interval as:

where , is absorber amount and (the diffusivity factor). More details can be found in Paltridge and Platt (l976), chapter 7. For the water vapor rotation bands (1-10), the Curtis-Godson approximation is employed to remove the effects of pressure on the transmission function along the optical path, du (see eqns (7.47) - (7.51) Paltridge and Platt, l976). An additional term is added to eqn (3.52) for bands 7-11, which partially covers the region of the water vapor continuum (band 21 completes the continuum region but is not included here, see section 3.4.2),

Note that the first term in (3.53) is zero for n=11, because A11=0. Eqn (3.53) is equivalent to a multiplication of exponential transmission functions, where the pressure weighting accounts for pressure effects on the continuum absorption coefficient along the optical path and

where (SK)n = 38.483, l6.288, 11.312, 7.428, 4.913 , for n=7,...,11. Equation (3.54) is from Roberts, et al. (l976) and accounts for the temperature, T, dependence of the absorption coefficient.

For the emissivity calculations (in eqn 3.51), equations (3.45) - (3.47) are used, where G1, G2, G3 are precomputed in tabular form using the Rogers and Walshaw parameters in Table 3.3 for strong absorption line transmission (see Paltridge and Platt, l976, eqn 7.38):

Water vapor amount scaled by is used during the table "look-up" for G1, G2, and G3. Water vapor rotational band 7 for these emissivity calculations is handled in the CO2 overlap (section 3.4.3). The contribution due to lines in bands 8-11 and 21 is neglected in eqn (3.55), however the effect of lines in bands 8-11 is included in the "exact" CTS computation (eqn 3.53). For the emissivity calculations, the continuum effects in bands 8-11 are included as a one broad band calculation and added to and , after the G1, G2, G3 table "look-up". Transmission in this one band continuum is:

where 8.6658 is a frequency weighted mean of the (SK)n (n=8-11) in eqn (3.54). Emissivity calculations for the continuum in band 7 are included in the CO2 overlap (section 3.4.3).

  1. Ozone
Ozone processes are not split into CTS and EXCHANGE terms. Parameterization of O3 effects is "exact" using the one interval random band model of Rodgers (l968). The transmission function in eqn (3.44) is computed for band 21 (table 3.3) as:

where p is pressure in units of atmospheres, =1.66, and and are absorber amounts. Using Rodgers notation

with =.28 cm-1, =.1 cm-1, k=208. cm g-1. The overlap with water vapor continuum (thus treated "exactly") is included as the last term in eqn (3.56) where

  1. Carbon Dioxide
Transmission functions for carbon dioxide are precomputed using techniques described by Fels and Schwarzkopf (l981), as updated in Schwarzkopf and Fels (l985). Detailed "line by line" procedures are used to obtain CO2 (330ppmv) transmission functions and their first and second temperature derivatives for the forecast model's vertical coordinate. Standard atmosphere temperature profiles and two specified surface pressures (1013.25 mb, 810.6 mb) are used. The CO2 transmission and derivatives are available in block data, BD2. Interpolation of the transmission, , and its derivatives to the vertical structure of each model grid point is simply:

where and represent the two model profiles available in block data. The interpolation of transmission functions to the temperature profile of the particular grid point is made from a second order expansion below:

where is from (3.57) and (eqn (33) in Fels and Schwarzkopf (1981)) is

Numerical experimentation shows that

where po=30 mb and p is model level pressure using psfc=1013.25 mb. The function G(p) is precomputed and resides in block data, BD2 (GTEMP).

A temperature correction to in (3.58) is necessary because the spectral interval (band 7) is not really narrow enough to consider the Planck function a constant. The final transmission function, is (eqn (7) in Schwarzkopf and Fels (l985)):


For the "exact" CTS heating rate calculation, in eqn (3.51), CO2 transmission is multiplied by the H2O band transmission (eqn 3.53)) to account for absorption overlap. For the emissivity EXCHANGE heating rate calculation,, in eqn (3.51), a one band version of eqn (3.44) is used. Overlapping H2O transmission is computed as (see eqn (3.54) and (3.55)):

where C7 is defined in eqn (3.54), and the strong line approximation is used for H2O band 7.

  1. Clouds
Three types of stratiform cloud are allowed (section 3.2.3) and for longwave processes all are considered blackbodies. Model levels are vertical boundaries for the clouds; however, unlike the shortwave parameterization, the bounding surfaces are considered to be in clear sky (Fig 3.6a). The longwave radiative flux at a particular level, p, is weighted by the percentage of clear sky in the optical paths from all atmospheric levels to level p. When more than one cloud type exists in a path, the clouds are randomly overlapped (Fig 3.6b).

In the case of thick (i.e., greater than one model layer) low cloud, longwave processes result in strong cooling in the upper cloud layer and much less cooling in the lower cloud layer. This produces a strong destabilizing effect, which could be harmful when used with the fixed zonal mean cloud climatology. An adjustment is made to the vertical gradient of longwave flux in the thick cloud (a smoothing!) so that it is constant for the cloudy part of the column. For completely overcast low cloud situations (i.e., cloud fraction =1), the resulting model layer longwave heating rates within the cloud will be identical.

  1. Adjustment to the Surface Downward Flux
When the approximate diurnal cycle is used for SW processes (section 3.3.6), changes to model temperature and H2O profiles in the lower atmosphere can be large. Since longwave heating rates are held fixed for each 12-hour interval, they could be seriously in error during parts of the model-day. There is, however, an adjustment to downward LW flux at the earth's surface ( ), for these large diurnal changes, which reduces the potential "12-hour shock" to the land surface temperature. Changes to are computed by assuming that the ratio is constant at each forecast model time step, where T1 is lowest model layer temperature. Thus when T1 increases (decreases) during the forecast, an estimate of the increase (decrease) of is obtained. The downward longwave flux is computed as in eqn (3.62), where ( )o denotes values at the beginning of each 12-hour radiation computation interval and ( )n are values at each model time-step.

Table 3.1

Discrete probability distribution,, of water vapor absorption coefficients, kn, for N=9 subintervals of the solar spectrum (see Lacis and Hansen (l974), table 1)

Table 3.2

Shortwave radiative properties of clouds for the 9 subintervals of the solar spectrum

Table 3.3

Water vapor band parameters used for exact CTS calculations. Bands (7-11,21) are H2O continuum and band 21 is for exact O3 computation, but neither is part of Rodgers and Walshaw (1966) model. H2O rotation bands (1-10) and 6.3 bands (12-20) from Rodgers and Walshaw. Adjustments made to bands 6-8 so that CO2 effects lie entirely within band 7.


Campana, K. A., l986: "An Approximation to the Diurnal Cycle for use in NMC's Medium-Range Forecast Model", Medium-Range Modeling Br. Note, unpublished.

Fels, S. B. and M. D. Schwarzkopf, l975: "The Simplified Exchange Approximation - A New Method for Radiative Transfer Calculations", J. of Atmos. Sci., pp.1475-1488.

Fels, S. B. and M. D. Schwarzkopf, l981: "An Efficient, Accurate Algorithm for Calculating CO2 15Band Cooling Rates", J. of Geophys. Res., pp. 1205-1232.

Green, J. S. A, l967: "Division of Radiative Streams into Internal Transfer and Cooling to Space", Quart. J. of Roy. Met. Soc., pp. 371-372.

Hering, W. S. and T. R. Borden, Jr., l965: "Mean Distribution of Ozone Density over North America, l963-l964", Environmental Research Papers, Report 162 USAF Cambridge Research Laboratory.

Lacis, A. A. and J. E. Hansen, l974: "A Parameterization for the Absorption of Solar Radiation in the Earth's Atmosphere", J. of Atmos. Sci., pp. 118-133.

London, J., l962: "Mesospheric Dynamics, 3, the Distribution of Total Ozone in the Northern Hemisphere, Final Report", Dept. of Meteorology/Oceanography, New York University.

Manabe, S. and R. F. Strickler, l964: "Thermal Equilibrium of an Atmosphere with a Convective Adjustment", J. of Atmos. Sci., pp. 361-385.

Matthews, E., l985: "Atlas of Archived Vegetation, Landuse, and Seasonal Albedo Data Sets", NASA Technical Memorandum 86199, Goddard Institute for Space Studies, New York.

Paltridge, G. W. and C. M. R. Platt, l976: "Radiative Processes in Meteorology and Climatology, Elsevier Scientific Publishing Company.

Payne, R. E., l972: "Albedo of the Sea Surface", J. of Atmos. Sci., pp. 959-970.

Roberts, R., J. Selby, and L. Biberman, l976: "Infrared Continuum Absorption by Atmospheric Water Vapor in the 8-12Window", Applied Optics, pp. 2085-2090.

Rodgers, C. D. and C. D. Walshaw, l966: "The Computation of Infrared Cooling Rate in Planetary Atmospheres", Quart. J. of Roy. Met. Soc., pp. 67-92.

Rodgers, C. D., l967: "The Radiative Heat Budget of the Troposphere and Lower Stratosphere", Report No. A2, Planetary Circulations Project, Dept. of Meteorology, M. I. T.

Rodgers, C. D., l968: "Some Extensions and Applications of the New Random Model for Molecular Band Transmission", Quart. J. of Roy. Met. Soc., pp. 99-102.

Sasamori, T., J. London and D. Hoyt, l972: "Radiation Budget of the Southern Hemisphere", Meteorological Monagraphs, vol. 13, number 35, pp. 9-23.

Schwarzkopf, M. D. and S. B. Fels, l985: "Improvements to the Algorithm for Computing CO2 Transmissivities and Cooling Rates", J. of Geophys. Res., pp. 10541-10550.

Stephens, G. L, l984: "The Parameterization of Radiation for Numerical Weather Prediction and Climate Models," Mon. Wea. Rev., pp. 826-867.

Telegadas, K and J. London, l954: "A Physical Model of the Northern Hemisphere Troposphere for Winter and Summer," Sci. Rpt. 1, Research Division, College of Engineering New York University, 55 pp.

Fig 3.1 Schematic of radiation calculations, with subroutine names.

Fig 3.2 Seasonal zonal mean ozone (parts per million) in upper nine layers of the 18 layer model. Mean layer pressure (mb) for surface pressure=1000 mb. Contour interval = 0.5 ppm.

Fig 3.3 Seasonal zonal mean cloud fraction for high, middle, and low cloud climatology. (a) Northern hemisphere winter. (b) Northern hemisphere spring. Values hemispherically flipped for opposite season.

Fig 3.4 Seasonal cloud thickness in model sigma layers (surface pressure = 1000 mb). High and middle clouds are one layer thick. (a) Northern hemisphere winter. (b) Northern hemisphere spring. Data hemispherically flipped for opposite season.

Fig 3.5 Multiple reflections for the short wave cloud calculation. Subscripts for TCLU and TCLD are inverted in subroutine SWR983 ; i.e. TCLU4 becomes TCLU1 , etc...

Fig 3.6 Model clouds. (a) Schematic representation of cloud bounding surfaces for short wave (SW) and long wave (LW) processes. (b) Cloud overlap, where H, M, L are cloud fraction for high, middle, and low cloud respectively.