## Numerical solution of the twostream equations

Eq. 5.27 is a coupled, linear system of ordinary differential equations requiring two boundary conditions. It takes the form of a two-point boundary value problem, because one specifies a boundary condition on I+ at t* form, the system can be written

In vector

where

The matrix M varies with t* because the single scattering albedo and also the asymmetry parameter in general will be functions of altitude. The forcing F depends on t* on account of the variation of T with altitude and the exponential factor in the direct-beam term. Because of these variations, the equations must generally be solved numerically. In most applications, the profiles of T and scatterer properties are given as functions of p, and often as tabulated values on a fixed grid rather than as functions. In order to carry out the needed integrations, the profiles must be re-expressed as functions of T* . The practicalities of how one goes about doing this will be discussed somewhat later. In this section we will show how to get the solution for a single frequency v, for which the absorption is characterized by a single absorption coefficient at each pressure level. The extension to the band-averaged case where we have to deal with a distribution of absorption coefficients will be dealt with in Section 5.9.

Because the system described by Eq. 5.51 is linear, the solution satisfying the boundary conditions can, in principle, be built up from a suitable superposition of solutions obtained by numerically integrating the system starting from the lower boundary. Using a numerical differential equation integrator, one first constructs the following three solutions:

• A particular solution Vpart which satisfies Eq. 5.51 subject to values of 1+ and 1— at t* =0 specified in whatever way proves convenient. We will assume 1+ = 1— =0 for the particular solution at the lower boundary, but almost any other choice would do as well.

• A pair of homogeneous solutions Vhom,o and Vhom,i which satisfy Eq. 5.51 with the forcing term F set equal to zero. We will construct these solutions by integrating the homogeneous form of the equation subject to lower boundary conditions 1+ = 1,1— = 0 for Vhomo and 1+ =0,1- = 1 for Vhom,1

The general solution to the original inhomogeneous system is then

where ao and a1 are any real numbers. At the upper boundary we require that the downward diffuse component vanish, i.e. 1—(t*) = 0. At the lower boundary, we require that Eq. 5.31 be satisfied. The superposition of basic solutions that satisfies both boundary conditions is determined by solving the system

0 = I— ,part (t* ) + 1-,hom,o((T* ))a0 + /-,hom,l((T* ))ai egnB(v, Tg) + a.gL© cos Zexp(—t*/ cos Z) = (1+,part(0) - agI-,part(0))

for the coefficients ao and a1. The conditions define a 2 x 2 linear system of equations which in general has a unique solution. Once the coefficients ao and a1 are known, the solution to the problem is complete. The procedure can equally well be carried out by starting at the top and integrating downward instead.

This approach works well and is very efficient as long as the atmosphere is not too optically thick. However, when t* becomes large, the method breaks down for the following reasons. The homogeneous version of Eq. 5.51 has solutions which grow or decay exponentially, with a local exponential growth or decay rate of ±\Jy2 — y2 . In the pure scattering limit y2 = y1 and the growth becomes merely algebraic, as evidenced by the analytic solutions considered previously.

However, in the general case where y2 < Yi, the exponentially growing solution causes numerical difficulties in the optically thick limit because the exponential growth will cause an overflow error when one attempts to integrate from the lower boundary to the upper boundary (or vice-versa). A related problem is that it becomes impossible to find two linearly independent homogeneous solutions because the the exponentially growing solution comes to dominate both homogenous solutions as one integrates sufficiently far from the boundary. As a simple example of this, consider that exp(T*) + aexp(—t*) eventually converges on exp(T*) to within computer roundoff error for sufficiently large t*, regardless of the value of a. Exactly how large t* needs to be before this problem becomes serious depends on the vertical profile of \pj\—YJ, but except when gaseous absorption is weak and scattering is dominant this quantity is order unity, and so the solution method tends to break down when t* « 10. Because of the exponential growth which is at the root of the problem, increasing the precision of the computer arithmetic only modestly extends this value.

Fortunately, there are a number of simple resolutions to the problem. First of all, if one only wants the OLR, then it is not necessary to integrate the equations all the way from the top of the atmosphere down to the ground. At a frequency where the atmosphere is optically thick by virtue of strong absorption, radiation from lower levels of the atmosphere is exponentially attenuated and does not significantly affect the OLR. In these circumstances, one starts the integration at a height t* — t* for some suitably chosen t* < taw*. At t* — t* one can impose any order unity boundary condition on I+ that proves convenient, since the boundary condition there won't affect the OLR as long as the fluxes are not made exponentially large. The lower boundary condition I+ =0 would suffice, though one could squeeze out a little more accuracy by using the optically thick limit for the boundary condition (i.e. I+ = nB(v, T(t* — t*) ). One carries out the rest of the procedure exactly as before, except for applying the lower boundary condition at the artificial lower boundary instead of t* =0. The choice of t* depends on the profile of \Jy2 — Y2. If the typical value of this quantity is A, then we would take t* « 10/A, based on the notion that e 10 + e_1 0 is accurately distinguishable from e1 0 with computer arithmetic having precision of at least twelve or thirteen decimal places. When \Jy2 — y2 varies greatly in the vertical it can be a bit tricky determining an appropriate "typical value," and so it is generally better to start the integration from the top, in which case one can simply integrate the two homogeneous solutions downward until one or the other grows by a factor of at least e10, whereupon one stops and applies the artificial lower boundary condition.

A similar idea can be applied if one wants to compute the fluxes deep in the interior of the atmosphere. The typical application of such a calculation would be to determine the radiative heating profile for use in a radiative-convective equilibrium calculation. Even in an all-troposphere model, for which the climate can be solved for knowing the OLR alone, one might wish to compute the interior infrared cooling rate so as to determine the magnitude of the convective heat transport required to balance radiative cooling.

Suppose one wishes to compute the fluxes in the vicinity of some level t*. The solution method exploits the fact that the fluxes near t* will be insensitive to conditions many optical thicknesses above or below t*, since the fluxes from distant layers will be exponentially attenuated by the time they reach t* . One then builds the solution from solutions constructed by integrating both upwards and downwards from taw* as follows:

• Integrating upward, one constructs a particular solution V>part and a homogeneous solution V>,h°m which satisfy the upper boundary condition I_ = 0. The upper boundary condition is applied at t* if the true upper boundary is reached before the exponentially growing mode dominates, but otherwise the integration is halted and the upper boundary condition is ap-

plied when the exponential dominance criterion is met. The particular solution satisfying the upper boundary condition is constructed from a superposition of an arbitrary particular solution and two independent homogeneous solutions, all obtained by upward integration starting from t* . Similarly, the homogeneous solution satisfying the upper boundary condition is obtained as a suitable superposition of two independent homogeneous solutions.

• Integrating downward , one constructs a particular solution V<,part and a homogeneous solution V<,hom which satisfy the appropriate lower boundary condition. If the true lower boundary t* = 0 is reached before the exponentially growing mode dominates, then the physical boundary condition in Eq. 5.31 is applied. Otherwise, an artificial boundary condition I+ = 0 is applied. As for the upward integration, the solutions are constructed through suitable superpositions of independent particular and homogeneous solutions obtained by downward integration.

The general solution satisfying the upper and lower boundary condition is then V>,part + aoV>,hom for t* > t* and V<,part + a1V<,hom for t* < t* . By imposing the requirement that the fluxes I+ and I- be continuous across t*, one obtains a 2 x 2 linear system for the coefficients ao and a1 much as we had before. This uniquely determines the solution. When the physical boundary isn't reached, then on the side of t* where this happens the solution will only be valid within a few optical thicknesses of t*, because the artificial boundary condition will begin to affect the solution as the artificial boundary is approached. To map out the entire flux profile, one only needs to carry out the procedure for a list of t* dense enough that the regions where the solutions are valid overlap. We'll refer to this solution method as the piecewise ODE method (ODE being shorthand for Ordinary Differential Equation). Note that in the optically thick limit, the fluxes typically vary very little over a unit optical thickness, because the fluxes are mostly determined by the local values of temperature and optical constants, which vary slowly when written in optical thickness coordinates. This means that it is usually possible to compute the flux at a rather coarse grid of t* and just use interpolation to get intermediate values if they are needed.

Once the fluxes have been determined, it is easy to get the heating rate, which is proportional to d(I+ — I-)/dT. Rather than numerically differentiating the fluxes, one can obtain the necessary derivatives by simply evaluating the right hand side of Eq. 5.51.

The procedure may sound complicated, but it is actually rather simple to implement on the computer, since it is built from the same basic integration operation carried out several times over, supplemented by a numerical routine that solves a 2 x 2 linear system. An example showing the results of this procedure is given in Fig. 5.8 . The solution is carried out for a single infrared wavenumber (600 cm-1) with gaseous absorption as described in the figure caption. In addition, a purely scattering cloud has been placed by the top of the atmosphere. For this profile, carrying out the procedure for four suitably chosen t* is sufficient to map out the flux profile over the entire domain. Note that the portion of the subintegrations which are contaminated by artificial boundary conditions are different for the upward and downward fluxes; for I+, it is the vicinity of the artificial lower boundary that is most affected, whereas for I- it is the vicinity of the upper boundary. Note also that the subintegration carried out nearest the top of the atmosphere covers a wider range of optical thickness than the others. This is because the optical thickness near the top is dominated by the purely scattering cloud, and pure scattering does not yield solutions with exponentially growing character.

In the calculations shown in Fig. 5.8 we used I+ = 0 or I- =0 for the artifical boundary conditions, so as to show more clearly the regions where the artificial conditions affect the solution. One can achieve considerably better accuracy close to the artificial boundaries by instead using the approximate optically thick solution I+ = I- = nB(v, T) as the boundary condition, where T

Upward Flux

Upward Flux Downward Flux Figure 5.8: An example of the results of computing the upward and downward flux using the piecewise ODE method in the optically thick case. The arrows show the values of to* for which the upward and downward integrations were performed. The shaded ellipses show the regions where the artificial boundary condition has significantly affected the solution. This calculation was carried out for a wavenumber of 600 cm—1 assuming the temperature profile to lie on the dry air adiabat with a ground temperature and surface air temperature of 270K. The surface pressure is 2 bar, the absorption coefficient is 1.5 • 10-4m2/kg at 100mb and linearly pressure-broadened elsewhere. An optically thick purely scattering cloud is included, centered at an optical thickness value of 35 . The calculation assumes a gravitational acceleration of 10 m/s2.

Figure 5.8: An example of the results of computing the upward and downward flux using the piecewise ODE method in the optically thick case. The arrows show the values of to* for which the upward and downward integrations were performed. The shaded ellipses show the regions where the artificial boundary condition has significantly affected the solution. This calculation was carried out for a wavenumber of 600 cm—1 assuming the temperature profile to lie on the dry air adiabat with a ground temperature and surface air temperature of 270K. The surface pressure is 2 bar, the absorption coefficient is 1.5 • 10-4m2/kg at 100mb and linearly pressure-broadened elsewhere. An optically thick purely scattering cloud is included, centered at an optical thickness value of 35 . The calculation assumes a gravitational acceleration of 10 m/s2.

is evaluated at the boundary.

The piecewise ODE method is flexible, accurate and easy to implement, and it deserves to be more widely used. A more customary approach to dealing with the numerical issues in the optically thick case is to turn the problem into a discrete set of linear equations and solve it using numerical linear algebra algorithms. To do this we set up a grid of discrete values Tj = jAT for j = 0,1, ...N, and approximate the derivative in Eq. 5.51 as (dV/dT*)j « (Vj+1 — V'-1)/At, where the integer subscripts stand for the value of the corresponding quantity at t* = jAT. With this approximation, Eq. 5.51 can be written as

where j = 1, 2, ...N — 1. The j = 0 and j = N lines must be left out because evaluation of the approximate derivative would require data off the end of the array. These two missing lines are replaced by the statement of the boundary condition at the bottom and top of the atmosphere respectively. This results in a 2N x 2N linear system if the upward and downward fluxes are merged into a single solution vector of the form [(I+,o, 1—,o, I+,1,1—,1,1+,2,1—,2,..-]. The corresponding matrix defining the coefficients of the system is band-diagonal, and has nonzero entries only within two spaces to the left or right of the diagonal. Such systems can be solved efficiently using standard methods of linear algebra available in virtually all numerical libraries 3. We'll refer to this solution technique a the matrix method.

Exercise 5.7.1 Write out the first three lines and the last three lines of the matrix problem outlined

3 See, for example the routine bandec in Numerical Recipes.

above, assuming to be constant.

Resolving the exponential growth or decay of the homogeneous solutions requires At < 1 at least, so it might be thought that a disadvantage of the matrix method is that it requires dealing with very large matrices when the t* is large - and with strong gaseous absorption, optical thicknesses of several thousand are not at all uncommon. However, the required resolution is not really determined by the exponential behavior, since, as noted previously, the fluxes are largely determined by local properties in the optically thick limit, but local properties such as temperature usually vary slowly when written as a function of t* in the optically thick case. Thus, the fluxes vary little over a unit optical depth, and the derivative of the flux can be accurately computed even if At is quite large. Another way of seeing this is to note that when At >> 1 the term Vj+i — Vj-i in Eq. 5.55 is negligible compared to the remainder of the terms, and the solution of what is left gives the correct lowest-order solution in the optically-thick limit. (Recall that the optically thick local solution is determined by treating the temperature and single-scatter albedo as if they were constant). Corrections to this approximate solution are of order 1/At, and will be properly computed in the linear solution algorithm. The only caveat needed is to note that besides the correct optically thick solution the linear system also has a spurious homogeneous two-gridpoint oscillating solution of the form I+(j) = [1, —1,1, —1,...] and similarly for I-. This spurious solution replaces the unresolved exponential homogenous solution, and care must be taken to use a linear system solver that correctly suppresses contamination by the spurious mode.

There are two other small technical details that need to be taken care of when using the matrix method in the optically thick limit. The first is that the conversion layer at the top of the atmosphere is not resolved, so that the direct-beam term is not correctly transformed into diffuse radiation. This is easily dealt with by eliminating the direct beam term and dumping the corresponding energy directly into I- (t* as an upper boundary condition, which explicitly captures the conversion of direct beam energy at the top of the atmosphere. Finally, if there is a significant temperature discontinuity between the ground temperature and the overlying air temperature, the exponential homogeneous solution comes into play in a kind of radiative boundary layer, and something must be done to explicitly resolve this layer, either analytically or by increasing the resolution near the boundary. When there is a temperature discontinuity, there is a strong, shallow radiative heating or cooling within a unit optical depth of the ground, and this will not be resolved by the matrix method using large At . With the ODE method, the radiative boundary layer is automatically resolved, but when using the matrix method one must generally use an analytical exponential solution to get the radiative fluxes near the ground.

Whichever method one chooses, it is often most convenient to work in optical thickness coordinates rather than transforming the equations to pressure or log-pressure space before performing the solution. Since the mapping between pressure and optical thickness is usually frequency dependent, it is necessary to express the fluxes as a function of pressure before summing up the fluxes and heating rates across frequencies. The easiest way to re-express the results is to use numerical integration to compute a suitable list of t* (p) from the differential equation defining t*, and then to use this list to define a numerical interpolation function giving p(t*) 4. Using this function a list of triples (I+ ,I-,t*) can be transformed to the list of triples (I+ ,I-,p) at the corresponding pressure levels. This list can then be used to interpolate the fluxes to a standard grid of pressure values; heating profiles are treated similarly. The interpolation function p(t*) fulfills an additional role when carrying out the integration for the fluxes. Namely, the temperature profile T and the scattering parameters are generally given as functions of p, or as tabulated values for a list of

4 By "interpolation function" we mean a computer-implemented function that takes a list of N values (xj ,yj) and an arbitrary argument x, and returns an interpolated estimate of the value of y corresponding to x. By putting in a list of (yj,xj) instead, the same function can be used to obtain x(y)

pressure levels. Using the interpolation function, however, one simply evaluates T(p(t*)) and so forth to get the required expression. If T(p) is specified as a table of values rather than a function, then one transforms the pressure entry to optical depth using the interpolation function, and then writes an interpolation function to get T(t*) from this table. When using the matrix method, it may not be necessary to write functions for T and the scattering parameters since tabulated values will generally be sufficient. When using the differential equation, method, however, it is usually necessary to provide functions, since the typical high order numerical integrator needs to be able to evaluate the right hand side of Eq. 5.51 at an arbitrary value of the vertical coordinate.

Clearly, in order to carry out the above procedure, it is necessary to have at hand an efficient, accurate and easy-to-use interpolation function. Some useful advice on interpolation methods has already been given in connection with the numerical analysis tutorial problems in the Workbook section of Chapter 1.