OSA's Digital Library

Applied Optics

Applied Optics

APPLICATIONS-CENTERED RESEARCH IN OPTICS

  • Vol. 27, Iss. 12 — Jun. 15, 1988
  • pp: 2502–2509
« Show journal navigation

Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media

Knut Stamnes, S-Chee Tsay, Warren Wiscombe, and Kolf Jayaweera  »View Author Affiliations


Applied Optics, Vol. 27, Issue 12, pp. 2502-2509 (1988)
http://dx.doi.org/10.1364/AO.27.002502


View Full Text Article

Acrobat PDF (936 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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]

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

Computer implementations of that method were, however, plagued by numerical difficulties [2]

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]

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]

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]

[9]

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).

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]

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).

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]

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

, [10]

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

and review articles. [11]

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

[15]

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

)

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

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

μduν(τν,μ,ϕ)dτ=uν(τν,μ,ϕ)-Sν(τν,μ,ϕ),
(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

Sν(τν,μ,ϕ)=ων(τν)4π02πdϕ-11dμPν(τν,μ,ϕ;μ,ϕ)×uν(τν,μ,ϕ)+Qν(τν,μ,ϕ),
(2)

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

Qν(thermal)(τν)=[1-ων(τν)]Bν[T(τν)],
(3)

where Bν(T) is the Planck function at frequency ν and temperature T. If the usual diffuse–direct distinction is made (Ref. [1]

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

, 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,

Qν(beam)(τν,μ,ϕ)=ων(τν)I04πPν(τν,μ,ϕ;-μ0,ϕ0)exp(-τν/μ0),
(4)

where μ0I0 is the incident flux. We take

Qν(τν,μ,ϕ)=Qν(thermal)(τν)+Qν(beam)(τν,μ,ϕ).

(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, [1]

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

, [9]

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).

u(τ,μ,ϕ)=m=02N-1um(τ,μ)cosm(ϕ0-ϕ),
(5)

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

μdum(τ,μ)dτ=um(τ,μ)--11Dm(τ,μ,μ)um(τ,μ)dμ-Qm(τ,μ)         (m=0,1,2,,2N-1),
(6a)

where

Dm(τ,μ,μ)=ω(τ)2l=m2N-1(2l+1)glm(τ)Plm(μ)Plm(μ)
(6b)
Qm(τ,μ)=X0m(τ,μ)exp(-τ/μ0)+δm0Q(thermal)(τ),
(6c)
X0m(τ,μ)=ω(τ)I04π(2-δm0)l=02N-1(-1)l+m(2l+1)×glm(τ)Plm(μ)Plm(μ0),
(6d)
δm0=1ifm=0(0otherwise),glm(τ)=gl(τ)(l-m)!(l+m)!,gl(τ)=12-11Pl(cosθ)P(τ,cosθ)dcosθ.

Pl(cosθ) is the Legendre polynomial, Plm(μ) 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]

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]

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]

μidum(τ,μi)dτ=um(τ,μi)-j=-Nj0NwjDm(τ,μi,μj)×um(τ,μj)-Qm(τ,μi)         (i=±1,±N),
(7a)

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.

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

|du+dτdu-dτ|=|-α-ββα||u+u-|+|Q˜+Q˜-|,
(7b)

where [3]

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]

, [16]

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]

u±=[um(τ,±μi)]         i=1,,N,Q˜±=M-1Q±,Q±=Qm(τ,±μi)         i=1,,N,M={μiδij}         i,j=1,,N,α=M-1(D+W-I),β=M-1D-W,W=(wiδij)         i,j=1,,N,D+=[Dm(μi,μj)]=[Dm(-μi,-μj)]         i,j=1,,N,D-=[Dm(-μi,μj)]=[Dm(μi,-μj)]         i,j=1,,N.

C. Quadrature Rule

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

j=-Nj0NwjD0(τ,μi,μj)=i=-Ni0NwiD0(τ,μi,μj)=ω(τ),
(7c)

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

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]

The quadrature points and weights of the double-Gauss scheme adopted here satisfy μj = −μj and wj = wj. Double-Gauss simply refers to a quadrature rule suggested by Sykes [18]

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

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 ( Q˜ = O) of Eq. (7b) of the form

u±=G±exp(-kτ),
(8a)

we find

|αβ-β-α||G+G-|=k|G+G-|,
(8b)

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. [3]

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]

Breaking Eq. (8b) into

αG++βG-=kG+,βG++αG-=-kG-,

and then adding and subtracing these two equations, we find

(α-β)(G+-G-)=k(G++G-),
(8c)
(α+β)(G++G-)=k(G+-G-).
(8d)

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

(α-β)(α+β)(G++G-)=k2(G++G-),
(8e)

which completes the reduction of the order. Following Stamnes and Swanson [3]

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]

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)(τ,μ) = X0(μ) exp(−τ/μ0) [Eq. (6d)], it is easily verified that a particular solution of Eq. (7) is (omitting m superscript)

u(τ,μi)=Z0(μi)exp(-τ/μ0).
(9a)

For thermal sources the emitted radiation is isotropic, so that Qm = 0, m > 0, Q(τ) ≡ Q0(τ) = (1 − ω)B(τ). By approximating the Planck function for each layer by a polynomial in τ, [3]

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]

, [19]

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

i.e.,

B[T(τ)]=l=0Kblτl,
(9a)

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

u(τ,μi)=l=0KYl(μi)τl.
(9b)

The coefficients Z0(μi) and Yl(μi) are readily determined by solving systems of linear algebraic equations (see Refs. [3]

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]

and [9]

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).

for details).

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):

u(τ,μi)=j=-NNCjGj(μi)exp(-kjτ)+δm0l=0NYl(μi)τl.
(10a)

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

C0G0(μi)Z0(μi),
(10b)
k01/μ0.
(10c)

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

u(τ=0,-μ,ϕ)=u(μ,ϕ),
(11a)
u(τ=τL,+μ,ϕ)=ug(μ,ϕ),
(11b)

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

u(τL,+μ,ϕ)=(μ)B(Tg)+1π02πdϕ01ρd(μ,ϕ;-μ,ϕ)u(τL,-μ,ϕ)μdμ+μ0πI0exp(-τL/μ0)ρd(μ,ϕ;-μ0,ϕ0),
(12)

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

(μ)+1π02πdϕ01ρd(μ,ϕ;-μ,ϕ)μdμ=1

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:

ρd(μ,ϕ;-μ,ϕ)=l=02N-1(2l+1)hlPl(cosθ),
(13a)

where

cosθ=-μμ+(1-μ2)1/2(1-μ2)1/2cos(ϕ-ϕ).
(13b)

Using the addition theorem for spherical harmonics, we find

ρd(μ,ϕ;-μ,ϕ)=m=02N-1ρdm(μ,-μ)cosm(ϕ-ϕ),
(13c)
ρdm(μ,-μ)=(2-δ0m)l=m2N-1(2l+1)hlmPlm(-μ)Plm(μ),
(14)

where

hlm=hl(l-m)!(l+m)!,hl=12-11Pl(cosθ)ρd(cosθ)dcosθ.
(15)

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

um(τL,+μ)=ugm(μ),
(16a)

where

ugm(μ)=δ0m(μ)B(Tg)+(1+δ0m)01μρdm(μ,-μ)um(τL,-μ)dμ+μ0πI0exp(-τL/μ0)ρdm(μ,-μ0).
(16b)

In a multilayered medium we must also require the intensity to be continuous across layer interfaces. [7]

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]

Thus Eq. (7) must satisfy boundary and continuity conditions as follows:

u1m(0,-μi)=um(-μi),         i=1,,N,
(17a)
upm(τp,μi)=up+1m(τp,μi),         i=±1,,±N;p=1,,L-1;
(17b)
uLm(τL,+μi)=ugm(μi),         i=1,,N,
(17c)

where

ugm(μi)=δ0m(μi)B(Tg)+(1+δ0m)j=1Nwjμjρdm(μi,-μj)um(τL,-μj)+μ0πI0exp(-τL/μ0)ρdm(μi,-μ0).
(17d)

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 kjp = −kjp):

up(τ,μi)=j=1N[CjpGjp(μi)exp(-kjpτ)+C-jpG-jp(μi)exp(+kjpτ)]+Rp(τ,μi),
(18)

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)]

Rp(τ,μi)=Z0(μi)exp(-τ/μ0)+δm0[Y0(μi)+Y1τ].
(19)

Insertion of Eq. (18) into Eqs. (17a)(17c) yields (omitting the m superscript)

j=1N[Cj1Gj1(-μi)+C-j1G-j1(-μi)]=u(-μi)-Rp(0,-μi),         i=1,,N
(20a)
j=1N{CjpGjp(μi)exp(-kjpτp)+C-jpG-jp(μi)exp(kjpτp)-[Cj,p+1Gj,p+1(μi)exp(-kj,p+1τp)+C-j,p+1G-j,p+1(μi)exp(kj,p+1τp)]}=[Rp+1(τp,μi)-Rp(τp,μi)],         i=±1,±N;p=1,,L-1
(20b)
j=1N[CjLrjL(μi)exp(-kjLτL)+C-jLr-jL(μi)exp(kjLτL)]=γ(τL,μi),         i=1,,N,
(20c)

where

γ(τL,μi)=δ0m(μi)B(Tg)-RL(τL,+μi)+μ0πI0exp(-τL/μ0)ρd(μi,μ0)+(1+δ0m)j=1Nρd(μi,-μj)wjμjRL(τL,-μj),
(21a)
rjL(μi)=GjL(μi)-(1+δ0m)×n=1Nρd(μi,-μn)wnμnGjL(-μn).
(21b)

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. [7]

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]

III. Numerical Implementation

A. Computation of Eigenvalues and Eigenvectors

Stamnes and Swanson [3]

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]

obtained the eigenvalues and eigenvectors of Eq. (8e) using subroutines in eispack [20]

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

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]

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]

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]

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]

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]

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]

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]

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]

, [16]

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]

introduce matrix operations in which the effect of rounding error can seriously degrade the results. Stamnes et al. [8]

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]

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:

  1. 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.
  2. 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]

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]

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 [7]

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]

C+jp=C˜+jpexp(kjpτp-1),
(22a)
C-jp=C˜-jpexp(-kjpτp).
(22b)

Inserting Eqs. (22) in Eqs. (20) and solving for the C˜jp 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 C˜jp is unconditionally stable for arbitrary total optical depths and arbitrary individual layer thicknesses. [To solve Eqs. (20) for the C˜jp we use the linpack [21]

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

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.,

up(τ,μi)=j=1N{C˜jpGjp(μi)exp[-kjp(τ-τp-1)]+C˜-jpG-jp(μi)exp[-kjp(τp-τ)]}.
(23)

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)

u(τ,+μ)=u(τL,+μ)exp[-(τL-τ)μ]+ττLS(t,+μ)exp[-(t-τ)/μ]dtμ,
(24a)
u(τ,-μ)=u(0,-μ)exp(-τ/μ)+0τS(t,-μ)exp[-(τ-t)/μ]dtμ,
(24b)

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]

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]

[5]

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

, [9]

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).

This procedure is sometimes referred to as the iteration of the source-function technique, although as pointed out by Kourganoff [22]

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

and discussed in some detail by Stamnes, [5]

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

this procedure essentially amounts to an interpolation.

In a multilayered medium we may now integrate Eq. (24a) layer by layer [5]

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

to obtain the following expression for the upward intensity at arbitrary τ and μ:

up(τ,+μ)=u(τL,+μ)exp[-(τL-τ)/μ]+n=pL{j=1N[C˜-jnG-jn(+μ)1+kjnμE-jn(τ,+μ)+C˜+jnG+jn(+μ)1+kjnμE+jn(τ,+μ)]]+G0n1+μ/μ0E0n(τ,+μ)+δm0[V0n(μ)F0n(τ,μ)+b1nF1n(τ,μ)]},
(25a)

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

E-jn(τ,+μ)=exp[-(kjnΔτn-δτ/μ)]-exp[-(τn-τ)/μ],
(25b)
F0n(τ,+μ)=exp(-δτ/μ)-exp[-(τn-τ)/μ],
(25c)
F1n(τ,+μ)=(τn-1+μ)exp(-δτ/μ)-(τn+μ)exp[-(τn-τ)/μ],
(25d)

with

δτ=τn-1-τforn>pandδτ=0forn=p;Δτn=τn-τn-1forn>pandΔτp=τp-τforn=p,E+jn(τ,+μ)=exp[-(τn-1-τ)/μ]-exp{-[kjn(τn-τn-1)-(τ-τn)/μ]},
(25e)

for n > p and

E+jp(τ,+μ)=exp[-kjp(τ-τp-1)]-exp{-[kjp(τp-τp-1)+(τp-τ)/μ]}.
(25f)

Similarly, we find that the downward intensity becomes

up(τ,-μ)=u(0,-μ)exp(-τ/μ)+n=1p{j=1N[C˜-jnG-jn(-μ)1+kjnμE-jn(τ,-μ)+C˜+jnG+jn(-μ)1-kjnμE+jn(τ,-μ)]+G0n(-μs)1-μ/μ0E0n(τ,-μ)+δm0[V0n(μ)F0n(τ,-μ)+b1nF1n(τ,-μ)]},
(26a)

where

E+jn(τ,-μ)=exp[-(kjnΔτn+δτ/μ)]-exp[-(τ-τn-1)/μ]
(26b)
F0n(τ,-μ)=exp(-δτ/μ)-exp[-(τ-τn-1)/μ]
(26c)
F1n(τ,-μ)=(τn-μ)exp(-δτ/μ)-(τn-1-μ)exp[-(τ-τn-1)/μ],
(26d)

with

δτ=τ-τnforn<pandδτ=0forn=p,Δτn=τn-τn-1forn<pandΔτ=τp-1forn=p,E-jn(τ,-μ)=exp[-(τ-τn)/μ]-exp{-[kjn(τn-τn-1)+(τ-τn-1)/μ]},
(26e)

for n < p and

E-jp(τ,-μ)=exp[-kjp(τp-τ)]-exp{-[kjp(τp-τn-1)+(τ-τp-1)/μ]}.
(26f)

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 Gjn terms and V¯0n(μ) in Eqs. (25a) and (26a) are given by the following explicit expressions (omitting the n subscript):

Gj(μ)=i=-Ni0NwiD(μ,μi)Gj(μi)forj0;
(27a)
G0(μ)Z0(μ)=i=-Ni0NwiD(μ,μi)Z0(μi)+X0(μ);
(27b)
V0(μ)=i=-Ni0NwiD0(μ,μi)Y0(μi)+(1-ω)b0;
(27c)

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

  1. 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;
  2. 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 [6]

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

:

a(μ)=u0(0,μ);a*(μ)=u0*(τL,-μ);
(28)
t(μ)=u0*(0,μ);t*(μ)=u0(τL,-μ);
(29)

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]

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

and [9]

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).

.)

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:

  1. It is unconditionally stable for an arbitrarily large number of quadrature angles (streams) and arbitrarily large optical depths.
  2. It can be forced by any combination of parallel beam or diffuse incidence and thermal sources in the medium and at the boundaries.
  3. It allows for an arbitrary bidirectional reflectivity at the lower boundary; directional emissivity is computed from this reflectivity.
  4. It offers rapid computation of bulk albedo and transmissivity when thermal sources are absent.
  5. Unlike the popular doubling method, computing time for individual layers is independent of optical thickness.
  6. 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.
  7. It has been thoroughly tested against a wide variety of published solutions.
  8. It is thoroughly documented both in Ref. [9]

    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).

    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

Fig. 1 Schematic illustration of a multilayered medium.

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).

History
Original Manuscript: August 7, 1987
Published: June 15, 1988

Citation
Knut Stamnes, S-Chee Tsay, Warren Wiscombe, and Kolf Jayaweera, "Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media," Appl. Opt. 27, 2502-2509 (1988)
http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-27-12-2502


Sort:  Author  |  Year  |  Journal  |  Reset  

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, 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, 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, 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, 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, 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, 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, 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, G. W. Stewart, linpack User’s Guide (SIAM, Philadelphia, 1979).
  22. V. Kourganoff, Basic Methods in Transfer Problems (Dover, New York, 1963).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited