## Abstract

We summarize an advanced, thoroughly documented, and quite general purpose discrete ordinate algorithm for time-independent transfer calculations in vertically inhomogeneous, nonisothermal, plane-parallel media. Atmospheric applications ranging from the UV to the radar region of the electromagnetic spectrum are possible. The physical processes included are thermal emission, scattering, absorption, and bidirectional reflection and emissionat the lower boundary. The medium may be forced at the top boundary by parallel or diffuse radiation and by internal and boundary thermal sources as well. We provide a brief account of the theoretical basis as well as a discussion of the numerical implementation of the theory. The recent advances made by ourselves and our collaborators—advances in both formulation and numerical solution—are all incorporated in the algorithm. Prominent among these advances are the complete conquest of two ill-conditioning problems which afflicted all previous discrete ordinate implementations: (1) the computation of eigenvalues and eigenvectors and (2) the inversion of the matrix determining the constants of integration. Copies of the fortran program on microcomputer diskettes are available for interested users.

© 1988 Optical Society of America

## I. Introduction

The discrete ordinate method for radiative transfer is commonly ascribed to Chandrasekhar.[1] Computer implementations of that method were, however, plagued by numerical difficulties[2] to such an extent that researchers made little use of it. The purpose of this paper is to alert the community to a new numerical implementation of the discrete ordinate method for vertically inhomogeneous layered media which is free of these difficulties and to give a summary of its equations and its various advanced features. The resulting computer code represents the culmination of years of effort[3]–[9] on the part of ourselves and our collaborators to make it the finest algorithm available. Our intent is that the code be so well documented, so versatile, and so error-free that other researchers can easily and safely use it both in data analysis and as a component of large models.

The problem to be solved is the transfer of monochromatic radiation in a scattering, absorbing, and emitting plane-parallel medium with a specified bidirectional reflectivity at the lower boundary. Section II summarizes the equations and boundary conditions. Section III discusses the numerical implementation of the theory, in particular, (a) how to compute eigenvalues and eigenvectors reliably and efficiently and (b) how to avoid fatal overflows and ill-conditioning in the matrix inversion needed to determine the constants of integration. Reference [9] provides a more detailed account as well as documentation, test problems, and a listing of the fortran code.

## II. Theory

#### A. Basic Equations

The purpose of this section is to present the basic radiative transfer formulas with a minimum of definition and explanation to establish notation and conventions. (For more comprehensive discussions of the subject, the reader is referred to textbooks[1],[10] and review articles.[11]–[15])

The equation describing the transfer of monochromatic radiation at frequency *ν* through a plane-parallel medium is given by[1]

where *u** _{ν}*(

*τ*

*,*

_{ν}*μ*,

*ϕ*) is the specific intensity along direction

*μ*,

*ϕ*at optical depth

*τ*

*measured perpendicular to the surface of the medium (*

_{ν}*ϕ*is the azimuthal angle, and

*μ*is the cosine of the polar angle).

*S*

*is the source function given by*

_{ν}where *ω** _{ν}*(

*π*

*) is the single-scattering albedo, and*

_{ν}*P*

*(*

_{ν}*τ*

*,*

_{ν}*μ*,

*ϕ;μ*′,

*ϕ*′) is the phase function. For thermal emission in local thermodynamic equilibrium (LTE), the source term

*Q*

*is*

_{ν}where *B** _{ν}*(

*T*) is the Planck function at frequency

*ν*and temperature

*T*. If the usual diffuse–direct distinction is made (Ref. [1], p. 22), so that

*u*

*in Eqs. (1) and (2) describes the diffuse radiation only, then for a parallel beam incident in direction*

_{ν}*μ*

_{0},

*ϕ*

_{0}on a nonemitting medium,

where *μ*_{0}*I*_{0} is the incident flux. We take

(Many radiative transfer models consider one or the other source term but not both.)

To simplify the notation we omit the *ν* subscript in the following. By expanding the phase function *P*(*τ*, cos*θ*) in a series of 2*N* Legendre polynomials and the intensity in a Fourier cosine series,[1],[9]

we find that Eqs. (1) and (2) are replaced by 2*N* independent equations (one for each Fourier component)

where

*P** _{l}*(

*cosθ*) is the Legendre polynomial, ${P}_{l}^{m}(\mu )$ is the associated Legendre polynomial, and

*θ*is the angle between the direction vectors before and after scattering. Solutions to Eq. (6) give the azimuthal components, and then Eq. (5) gives the complete azimuthal dependence of the intensity. [When there is no beam source, the sum in Eq. (5) reduces to the

*m*= 0 term, and

*μ*

_{0},

*ϕ*

_{0}are irrelevant.]

#### B. Discrete Ordinate Approximation—Matrix Formulation

The discrete ordinate approximation to (6) can be written[3],[4]

where *μ** _{i}* and

*w*

*are quadrature points and weights.*

_{i}Since the single-scattering albedo *ω*(*τ*) and the phase function *P*(*τ*,*μ*,*ϕ;μ*′,*ϕ*′) are functions of *τ* in a vertically inhomogeneous medium, Eqs. (7) constitute (for each *m*) a system of 2*N* coupled differential equations with nonconstant coefficients for which analytic solutions do not exist. To obtain analytic solutions, the medium is assumed to consist of *L* adjacent homogeneous layers in which the single-scattering albedo and the phase function are taken to be constant within each layer (but allowed to vary from layer to layer, see Fig. 1), and the thermal source term is approximated by a polynomial in *τ*. For the time being, we shall consider a single homogeneous layer for which *τ*_{p}_{−1} ≤ *τ* ≤ *τ** _{p}*. The

*τ*arguments of

*D*and

*X*will be omitted since

*D*and

*X*are by assumption independent of

*τ*in any one layer. We will also omit the

*p*subscripts, which are implicit on all quantities, until later.

For any layer in Fig. 1, then, we write Eq. (7) in matrix form as

#### C. Quadrature Rule

The use of Gaussian quadrature is essential because it makes phase function renormalization unnecessary, i.e.,

implying that energy is conserved in the computation.[17]

The quadrature points and weights of the double-Gauss scheme adopted here satisfy *μ*_{−}* _{j}* = −

*μ*

*and*

_{j}*w*

_{−}

*=*

_{j}*w*

*. Double-Gauss simply refers to a quadrature rule suggested by Sykes[18] in which the Gaussian formula is applied separately to the half-ranges −1 <*

_{j}*μ*< 0 and 0 <

*μ*< 1. The main advantage of this double-Gauss scheme is that the quadrature points (in even orders) are distributed symmetrically around |

*μ*| = 0.5 and clustered both toward |

*μ*| = 1 and

*μ*= 0, whereas in the Gaussian scheme for the complete range, −1 <

*μ*< 1, they are clustered toward

*μ*= ±1. The clustering toward

*μ*= 0 will give superior results near the boundaries where the intensity varies rapidly around

*μ*= 0. A half-range scheme is also preferred since the intensity is discontinuous at the boundaries. Another advantage is that upward and downward fluxes and average intensities are obtained immediately without further approximations.

#### D. Basic Solution

Equation (7b) is a system of 2*N*-coupled ordinary differential equations with constant coefficients. These coupled equations are linear, and our goal is to uncouple them by using well-known methods of linear algebra.

Seeking solutions to the homogeneous version (
$\tilde{\mathbf{Q}}$ = **O**) of Eq. (7b) of the form

we find

which is a standard algebraic eigenvalue problem of order 2*N* × 2*N* determining the eigenvalues *k* and the eigenvectors **G**^{±}.

Because of the special structure of the matrix in Eq. (8b), the eigenvalues occur in pairs (±*k*), and the order of the algebraic eigenvalue problem [Eq. (8b)] may be reduced by a factor of 2 as follows.[3] Breaking Eq. (8b) into

and then adding and subtracing these two equations, we find

Combining Eqs. (8c) and (8d), we obtain

which completes the reduction of the order. Following Stamnes and Swanson[3] we solve Eq. (8e) to obtain eigenvalues and eigenvectors (**G**^{+} + **G**^{−}). We then use Eq. (8d) to determine (**G**^{+} − **G**^{−}).

For beam sources *Q*^{(beam)}(*τ*,*μ*) = *X*_{0}(*μ*) exp(−*τ*/*μ*_{0}) [Eq. (6d)], it is easily verified that a particular solution of Eq. (7) is (omitting *m* superscript)

For thermal sources the emitted radiation is isotropic, so that *Q** ^{m}* = 0,

*m*> 0,

*Q*(

*τ*) ≡

*Q*

^{0}(

*τ*) = (1 −

*ω*)

*B*(

*τ*). By approximating the Planck function for each layer by a polynomial in

*τ*,[3],[19] i.e.,

we find that the particular solution may also be expressed as a polynomial in *τ*:

The coefficients *Z*_{0}(*μ** _{i}*) and

*Y*

*(*

_{l}*μ*

*) are readily determined by solving systems of linear algebraic equations (see Refs. [3] and [9] for details).*

_{i}The general solution to Eq. (7a) consists of a linear combination, with coefficients *C** _{j}*, of all the homogeneous solutions plus the particular solutions for beam and thermal emission sources (omitting the

*m*superscript):

The *k** _{j}* and the

*G*

*(*

_{j}*μ*

*) are the eigenvalues and eigenvectors obtained as described in Sec. III.A, the*

_{i}*μ*

*the quadrature angles, and the*

_{i}*C*

*the constants of integration. We have for convenience included the beam solution in the sum by defining for*

_{j}*j*= 0

#### E. Boundary Conditions

Below we shall limit ourselves to a linear-in-optical-depth dependence of the Planck function within each layer, since this will give sufficient accuracy for most cases of practical interest. We shall discuss the inclusion of boundary conditions in some detail, since a proper formulation of this aspect of the problem (in terms of the bidirectional reflectivity) is rarely discussed in the open literature.

In general, Eq. (1) must be solved subject to the boundary conditions

where *u*_{∞} and *u** _{g}* are the intensities incident at the top and bottom boundaries, respectively, and

*τ*

*is the total optical depth.*

_{L}Quite general boundary conditions can be accommodated in the discrete ordinate method. To be specific, we shall assume that the medium is illuminated from above by known diffuse radiation (incident parallel beams are treated as pseudosources) and that the bottom boundary has a known bidirectional reflectivity. Equation (11b) may now be rewritten as

where *∊*(*μ*) is the directional emissivity and *T** _{g}* is the temperature of the bottom boundary;

*ρ*

_{d}(

*μ*,

*ϕ*;−

*μ*′,

*ϕ*′) is the bidirectional reflectivity, and

*I*

_{0}is the incident beam intensity at the upper boundary. Kirchhoff’s law

is applied to calculate the emissivity from the reflectivity. Note that some authors include the (1/*π*) factor in *ρ** _{d}*; we find our definition more useful, however, since when

*ρ*

*is independent of angle it reduces to an albedo.*

_{d}Below we will assume that the bidirectional reflectivity is a function only of the angle *θ* between the incident and the reflected radiation (i.e., there are no preferred directions at the lower boundary). Then it can be expanded in a series of 2*N* Legendre polynomials:

where

Using the addition theorem for spherical harmonics, we find

where

Substituting Eq. (13c) into Eq. (12) and using Eq. (5), we find that each Fourier component must satisfy the bottom boundary condition

where

In a multilayered medium we must also require the intensity to be continuous across layer interfaces.[7] Thus Eq. (7) must satisfy boundary and continuity conditions as follows:

where

Evaluation of the integral in Eq. (16b) by numerical quadrature as in Eq. (17d) requires a quadrature rule which integrates separately over the downward and upward directions. As noted previously, the double-Gauss rule adopted here satisfies this requirement.

When discussing boundary conditions, it is convenient to write the discrete ordinate solution in the following form (*k** _{jp}* > 0 and

*k*

_{−}

*= −*

_{jp}*k*

*):*

_{jp}where the sum contains the homogeneous solution involving the unknown coefficients (the *C** _{jp}*) to be determined, and

*R*

*(*

_{p}*τ*,

*μ*

*) is the particular solution given by [Eq. (10a)]*

_{i}Insertion of Eq. (18) into Eqs. (17a)–(17c) yields (omitting the *m* superscript)

where

Equations (20a)–(20c) constitute a (2*N* × *L*) × (2*N* × *L*) system of linear algebraic equations from which the 2*N* × *L* unknown coefficients, the *C** _{jp}* (

*j*= ±1,…,±

*N; p*= 1,…,

*L*) are determined. To solve this system of equations, we take advantage of the fact that the coefficient matrix is a (6

*N*− 1) diagonal band matrix.[7]

## III. Numerical Implementation

#### A. Computation of Eigenvalues and Eigenvectors

Stamnes and Swanson[3] obtained the eigenvalues and eigenvectors of Eq. (8e) using subroutines in eispack[20] for a general real asymmetric matrix. Since such matrices may have complex eigenvalues/eigenvectors, the eispack routines use complex arithmetic. Yet invariably the eigenvalues returned were purely real as they should be.

As discussed by Stamnes *et al*.,[8] the advantage of solving the eigenvalue problem involving the asymmetric matrix (*α* − *β*) (*α* + *β*) in Eq. (8e) is that only one matrix multiplication is necessary. Nakajima and Tanaka[16] showed that the eigenvalue problem can be reduced to finding eigenvalues and eigenvectors of the product of two symmetric matrices, one of which is positive definite. Stamnes *et al*.[8] used the Cholesky decomposition of this positive definite matrix to convert the eigenvalue problem into one involving only one symmetric matrix. But the transformation to symmetric matrices (whose eigensolutions can be found faster and by using real arithmetic) and subsequent solution procedures[8],[16] introduce matrix operations in which the effect of rounding error can seriously degrade the results. Stamnes *et al*.[8] compared the original Stamnes/Swanson method with the new procedure based on the Cholesky decomposition and with Nakajima and Tanaka’s symmetrizing procedure. They summarize the outcome of these comparisons as follows:

- In single precision the Stamnes/Swanson procedure is more accurate than the Cholesky procedure, which in turn is slightly more accurate than the Nakajima/Tanaka procedure.
- The Cholesky procedure is slightly faster than Nakajima/Tanaka’s, which is slightly faster than Stamnes/Swanson’s.

Based on the findings, they decided to speed up the Stamnes/Swanson procedure since it is the most accurate. By eliminating complex arithmetic in the eispack subroutines, the speed became comparable with that of the Cholesky approach. We have adopted this faster algorithm in our code.

#### B. Scaling Transformation

As discussed by Stamnes and Conklin,[7] to avoid numerical ill-conditioning, it is necessary to remove the positive exponentials in Eq. (20) (remember *k** _{jp}* > 0 by convention). This is achieved by the scaling transformation[7]

Inserting Eqs. (22) in Eqs. (20) and solving for the
${\tilde{C}}_{jp}$ instead of the *C** _{jp}*, we find that all the exponential terms in the coefficient matrix have the form exp[−

*k*

*(*

_{jp}*τ*

*−*

_{p}*τ*

_{p}_{−1})];i.e., all the exponentials in the coefficient matrix have negative arguments (

*k*

*> 0,*

_{jp}*τ*

*>*

_{p}*τ*

_{p}_{−1}). Consequently, numerical ill-conditioning is avoided implying that our numerical scheme for finding the ${\tilde{C}}_{jp}$ is unconditionally stable for arbitrary total optical depths and arbitrary individual layer thicknesses. [To solve Eqs. (20) for the ${\tilde{C}}_{jp}$ we use the linpack[21] routine sgbsl designed to solve a banded system of linear algebraic equations.]

#### C. Scaled Solutions—Angular Distributions

Since only the homogeneous solution is affected, we start by introducing the scaling transformation (22) into Eq. (18), i.e.,

Since *k** _{jp}* > 0 and

*τ*

_{p}_{−1}≤

*τ*≤

*τ*

*, all exponentials in Eq. (23) have negative arguments as they should to avoid fatal overflows.*

_{p}For a slab of thickness *τ** _{L}*, we may solve Eq. (1) formally to obtain (

*μ*> 0)

where we have again omitted the *m* superscript. These two equations show that if we know the source function *S*(*t*,*±μ*), we can find the intensity at arbitrary angles by integrating the source function. The discrete ordinate solutions are used to derive explicit expressions for the source function, which can be integrated analytically.[3]–[5],[9] This procedure is sometimes referred to as the iteration of the source-function technique, although as pointed out by Kourganoff[22] and discussed in some detail by Stamnes,[5] this procedure essentially amounts to an interpolation.

In a multilayered medium we may now integrate Eq. (24a) layer by layer[5] to obtain the following expression for the upward intensity at arbitrary *τ* and *μ*:

where the scaling given by Eqs. (22) has been incorporated and

with

for *n* > *p* and

Similarly, we find that the downward intensity becomes

where

with

for *n* < *p* and

We see that all exponentials in Eqs. (25b)–(25f) and (26b)–(26f) have negative arguments as they should. This ensures that fatal overflow errors are avoided in the computations.

The *G** _{jn}* terms and
${\overline{V}}_{0n}(\mu )$ in Eqs. (25a) and (26a) are given by the following explicit expressions (omitting the

*n*subscript):

where *b*_{0} is defined in Eq. (9b). Formulas (27a)–(27c) clearly reveal the interpolatory nature of Eqs. (25) and (26). Although we have generally omitted the *m* superscript above, we have explicitly written *D*^{0} in Eq. (27c) to indicate that thermal radiation contributes only to the azimuth-independent part of the intensity.

#### D. Computational Shortcut for Highly Absorbing Media

For highly absorbing media without thermal emission, very little radiation will reach levels for which the absorbing optical depth is greater than, say, 10. In such cases we ignore the medium below this depth and set the surface reflection to zero. This may result in enormous computation savings. For example, in an atmosphere with clouds it prevents unnecessary scattering computations at wavelengths where sunlight never reaches the clouds.

#### E. Simplified Albedo and Transmissivity Computations

By appealing to the reciprocity principle we may derive simple expressions for the plane albedo *a*(*μ*) and the transmissivity *t*(*μ*) of a vertically inhomogeneous plane-parallel medium lacking thermal sources. Thus

- the plane albedo for a given angle of parallel incidence equals the azimuthally averaged reflected intensity
*u*^{0}(0,*μ*) for isotropic illumination of unit intensity incident at the top boundary; - the transmissivity equals the azimuthally averaged transmitted intensity
*u*^{0}*(0,*μ*) for isotropic illumination of unit intensity incident at the bottom boundary.

Consequently, for an inhomogeneous medium the following duality relations apply[6]:

where *τ** _{L}* is to the optical thickness of the inhomogeneous medium, and the asterisk refers to illumination from below.

Equations (28) and (29) offer substantial computational advantages when only albedo and transmissivity are required. The conventional procedure requires computation of the angular distribution of the diffuse intensity, which must then be integrated over angle to obtain the desired result. This necessitates the computation of a particular solution, which can be quite costly in an inhomogeneous (multilayered) medium, since such a solution is required for each layer and every direction of incidence considered. In contrast, by applying an isotropic boundary condition and using Eqs. (28) and (29), one may solve for all desired angles of parallel beam incidence simultaneously and avoid entirely the computation of a particular solution. (Details about the numerical implementation of this procedure are provided in Refs. [6] and [9].)

## IV. Summary

We have described an advanced discrete ordinate algorithm for radiative transfer calculations, including multiple scattering, in vertically inhomogeneous nonisothermal plane-parallel media. This algorithm is the product of over seven years of development by ourselves and our collaborators. It has been designed to be as versatile and general as possible and should find applications from the UV through to the radar region of the electromagnetic spectrum.

The algorithm has the following important features:

- It is unconditionally stable for an arbitrarily large number of quadrature angles (streams) and arbitrarily large optical depths.
- It can be forced by any combination of parallel beam or diffuse incidence and thermal sources in the medium and at the boundaries.
- It allows for an arbitrary bidirectional reflectivity at the lower boundary; directional emissivity is computed from this reflectivity.
- It offers rapid computation of bulk albedo and transmissivity when thermal sources are absent.
- Unlike the popular doubling method, computing time for individual layers is independent of optical thickness.
- Because the solution is analytic, and by using the iteration-of-the-source-function method, intensities can be returned at arbitrary angles and optical depths, unrelated to the computational meshes for these quantities.
- It has been thoroughly tested against a wide variety of published solutions.
- It is thoroughly documented both in Ref. [9] and in the code itself (with extensive references to published equation numbers).

Copies of the fortran program are available from the third author on microcomputer diskettes (IBM or MacIntosh).

This work was supported by the National Science Foundation through grant DPP 84-06093 to the University of Alaska. One of us (K.S.) would like to acknowledge partial support of this work from the Royal Norwegian Council on Science and Technology.

## Figure

## References

**1. **S. Chandrasekhar, *Radiative Transfer* (Dover, New York, 1960).

**2. **K.-N. Liou, “A Numerical Experiment on Chandrasekhar’s Discrete-Ordinate Method for Radiative Transfer: Applications to Cloudy and Hazy Atmosphere,” J. Atmos. Sci. **30**, 1303 (1973). [CrossRef]

**3. **K. Stamnes and R. A. Swanson, “A New Look at the Discrete Ordinate Method for Radiative Transfer Calculations in Anisotropically Scattering Atmospheres,” J. Atmos. Sci. **38**, 387 (1981). [CrossRef]

**4. **K. Stamnes and H. Dale, “A New Look at the Discrete Ordinate Method for Radiative Transfer Calculation in Anisotropically Scattering Atmospheres, II: Intensity Computations,” J. Atmos. Sci. **38**, 2696 (1981). [CrossRef]

**5. **K. Stamnes, “On the Computation of Angular Distributions of Radiation in Planetary Atmospheres,” J. Quant. Spectrosc. Radiat. Transfer **28**, 47 (1982a). [CrossRef]

**6. **K. Stamnes, “Reflection and Transmission by a Vertically Inhomogeneous Planetary Atmosphere,” Planet. Space Sci. **30**, 727 (1982b). [CrossRef]

**7. **K. Stamnes and P. Conklin, “A New Multi-Layer Discrete Ordinate Approach to Radiative Transfer in Vertically Inhomogeneous Atmospheres,” J. Quant. Spectrosc. Radiat. Transfer **31**, 273 (1984). [CrossRef]

**8. **K. Stamnes, S.-C. Tsay, and T. Nakajima, “Computation of Eigenvaluesand Eigenvectors for Discrete Ordinate and Matrix Operator Method Radiative Transfer,” J. Quant. Spectrosc. Radiat. Transfer in press (1988). [CrossRef]

**9. **K. Stamnes, S.-C. Tsay, W. J. Wiscombe, and K. Jayaweera, “An Improved, Numerically Stable Computer Code for Discrete-Ordinate-Method Radiative Transfer in Scattering and Emitting Layered Meda,” NASA Report, in press (1988).

**10. **K. N. Liou, *An Introduction to Atmospheric Radiation* (Academic, Orlando, FL, 1980).

**11. **J. E. Hansen and L. D. Travis, “Light Scattering in Planetary Atmospheres,” Space Sci. Rev. **16**, 527 (1974). [CrossRef]

**12. **W. M. Irvine, “Multiple Scattering in Planetary Atmospheres,” Icarus **25**, 175 (1975). [CrossRef]

**13. **J. Lenoble, Ed., *Radiative Transfer in Scattering and Absorbing Atmospheres: Standard Computational Procedures* (A. Deepak, Hampton, VA, 1985).

**14. **W. J. Wiscombe, “Atmospheric Radiation: 1975–1983,” Rev. Geophys. **21**, 957 (1983). [CrossRef]

**15. **K. Stamnes, “The Theory of Multiple Scattering of Radiation in Plane Parallel Atmospheres,” Rev. Geophys. **24**, 299 (1986). [CrossRef]

**16. **T. Nakajima and M. Tanaka, “Matrix Formulations for the Transfer of Solar Radiation in a Plane-Parallel Atmosphere,” J. Quant. Spectrosc. Radiat. Transfer **35**, 13 (1986). [CrossRef]

**17. **W. J. Wiscombe, “The Delta-M Method: Rapid Yet Accurate Radiative Flux Calculations for Strongly Asymmetric Phase Functions,” J. Atmos. Sci. **34**, 1408 (1977). [CrossRef]

**18. **J. B. Sykes, “Approximate Integration of the Equation of Transfer,” Mon. Not. Roy. Astron. Soc. **11**, 377 (1951).

**19. **W. J. Wiscombe, “Extension of the Doubling Method to Inhomogeneous Sources,” J. Quant. Spectrosc. Radiat. Transfer **16**, 477 (1976). [CrossRef]

**20. **W. R. Cowell, Ed., *Sources and Developments of Mathematical Software* (Prentice Hall, Englewood Cliffs, NJ, 1980).

**21. **J. J. Dongarra, C. B. Moler, J. R. Bunch, and G. W. Stewart, *linpack User’s Guide* (SIAM, Philadelphia, 1979).

**22. **V. Kourganoff, *Basic Methods in Transfer Problems* (Dover, New York, 1963).