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
The discrete ordinate method for radiative transfer is commonly ascribed to Chandrasekhar. Computer implementations of that method were, however, plagued by numerical difficulties 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– 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  provides a more detailed account as well as documentation, test problems, and a listing of the fortran code.
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, and review articles.–)
The equation describing the transfer of monochromatic radiation at frequency ν through a plane-parallel medium is given by
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. , 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 μ0I0 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 2N Legendre polynomials and the intensity in a Fourier cosine series,,
Pl(cosθ) is the Legendre polynomial, 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
where μi and wi are quadrature points and weights.
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 2N 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.
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.
The quadrature points and weights of the double-Gauss scheme adopted here satisfy μ−j = −μj and w−j = wj. Double-Gauss simply refers to a quadrature rule suggested by Sykes in which the Gaussian formula is applied separately to the half-ranges −1 < μ < 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 2N-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 ( = O) of Eq. (7b) of the form
which is a standard algebraic eigenvalue problem of order 2N × 2N 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. Breaking Eq. (8b) into
and then adding and subtracing these two equations, we find
we find that the particular solution may also be expressed as a polynomial in τ:
The general solution to Eq. (7a) consists of a linear combination, with coefficients Cj, of all the homogeneous solutions plus the particular solutions for beam and thermal emission sources (omitting the m superscript):
The kj and the Gj(μi) are the eigenvalues and eigenvectors obtained as described in Sec. III.A, the μi the quadrature angles, and the Cj the constants of integration. We have for convenience included the beam solution in the sum by defining for 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 ug are the intensities incident at the top and bottom boundaries, respectively, and τL is the total optical depth.
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 Tg is the temperature of the bottom boundary; ρd(μ,ϕ;−μ′,ϕ′) is the bidirectional reflectivity, and I0 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 ρd is independent of angle it reduces to an albedo.
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 2N Legendre polynomials:
Using the addition theorem for spherical harmonics, we find
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 (kjp > 0 and k−jp = −kjp):
where the sum contains the homogeneous solution involving the unknown coefficients (the Cjp) to be determined, and Rp(τ,μi) is the particular solution given by [Eq. (10a)]
Equations (20a)–(20c) constitute a (2N × L) × (2N × L) system of linear algebraic equations from which the 2N × L unknown coefficients, the Cjp (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 (6N − 1) diagonal band matrix.
III. Numerical Implementation
A. Computation of Eigenvalues and Eigenvectors
Stamnes and Swanson obtained the eigenvalues and eigenvectors of Eq. (8e) using subroutines in eispack 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., 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 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. 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, introduce matrix operations in which the effect of rounding error can seriously degrade the results. Stamnes et al. 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, to avoid numerical ill-conditioning, it is necessary to remove the positive exponentials in Eq. (20) (remember kjp > 0 by convention). This is achieved by the scaling transformation
Inserting Eqs. (22) in Eqs. (20) and solving for the instead of the Cjp, we find that all the exponential terms in the coefficient matrix have the form exp[−kjp(τp − τp−1)];i.e., all the exponentials in the coefficient matrix have negative arguments (kjp > 0, τp > τp−1). Consequently, numerical ill-conditioning is avoided implying that our numerical scheme for finding the is unconditionally stable for arbitrary total optical depths and arbitrary individual layer thicknesses. [To solve Eqs. (20) for the we use the linpack routine sgbsl designed to solve a banded system of linear algebraic equations.]
C. Scaled Solutions—Angular Distributions
Since kjp > 0 and τp−1 ≤ τ ≤ τp, all exponentials in Eq. (23) have negative arguments as they should to avoid fatal overflows.
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.–, This procedure is sometimes referred to as the iteration of the source-function technique, although as pointed out by Kourganoff and discussed in some detail by Stamnes, this procedure essentially amounts to an interpolation.
where the scaling given by Eqs. (22) has been incorporated and
for n > p and
Similarly, we find that the downward intensity becomes
for n < p and
where b0 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 D0 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 u0(0,μ) for isotropic illumination of unit intensity incident at the top boundary;
- the transmissivity equals the azimuthally averaged transmitted intensity u0*(0,μ) for isotropic illumination of unit intensity incident at the bottom boundary.
Consequently, for an inhomogeneous medium the following duality relations apply:
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.  and .)
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.  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.
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).