## Ld2fT2X Mrx SXTX695

with the angle-independent source function S\ given by

where ux is the single-scattering albedo and rx is the extinction optical depth. When the expression for the source function is inserted into the above diffusion equation we get l^f1 = C1 " ^)(Mrx) ~ Bx(rx)). (6.97)

Given appropriate boundary conditions, one can rapidly solve this boundary-valued problem using a numerical technique such as the Thomas algorithm, described below, if one replaces the second derivative by a three-term finite-differences expression.

### 6.7.2 Thomas algorithm

The Thomas algorithm is Gauss elimination applied to the inversion of a tridi-agonal matrix (see Williams 1972). If the first and second derivatives involved in the radiation diffusion equation are replaced by finite differences and keeping only three terms then the Thomas algorithm can be applied and is very rapid. We may replace first and second derivatives by expressions of the form

so that the radiation diffusion equation at atmospheric level i takes the form

We then solve the linear system of n equations

where J and V are vectors and A is a matrix of form

with imposed upper and lower boundary conditions to provide information on J0 and Jn+i, required by the finite-difference expressions for the first and second derivatives.

6.7.2.1 Solution procedure For i = 1 we set ai = ai Yi = ci/ai ui = vi/ai , for i = 2, 3, ...n we set ai = ai — biYi-i, Ui = (vi — biUi-i)/a.i whilst for i = 2, 3, ...n — 1 we have

with the back-substitution procedure for i = n — 1,n — 2, ...1

Note that matrix inversion operation counts (and hence computer time) depend on n3, whilst the Thomas algorithm operation counts vary as n. This makes the Thomas algorithm relatively fast, especially if one needs to solve for many atmospheric layers. Furthermore, one can apply the method to a tridiagonal matrix where the matrix elements are themselves matrices, an example is the solution of the continuity equation for many species concentrations in atmospheric photochemistry (see Chapter 7, also Lavvas et al. 2007). In this case divisions by scalar factors are replaced by matrix inversions in the Thomas solution procedure. Generally, one can speed up matrix inversions by partitioning the matrix into four submatrices and so reduce the size by 2 and the time for each inversion by 8, at the expense of more matrix multiplications.

### 6.7.3 Anisotropic scattering solution

6.7.3.1 Truncated phase-function approximation Shettle and Weinman (1970) expressed the scattering phase function, p(\$), as a two-term truncated series of associated Legendre polynomials in the manner of Chandrasekhar (1960) to obtain p(\$) = 1 + 3g cos ■&, (6.106)

where \$ is the scattering angle. This form arises by applying two conditions, that this phase function is normalized and that it satisfies eqn (3.58). For an azimuthally averaged radiation field, we can average also the phase function

over the incoming azimuthal direction by replacing cos ß given by eqn (3.69) to obtain p(p, p ) = 1 + 3gpp , (6.107)

where p is the cosine of the incoming zenith direction and p that of the outgoing zenith direction. The angle-dependent source function can now be expressed as

where the thermal emission source function is

SeX = (1 - ux)Bx, (6.109) the scattering source function is ux f1 ' ' '

and the source function for the incoming direct solar radiation is