DOWNWARD SURFACE FLUX COMPUTATIONS IN A VERTICALLY INHOMOGENEOUS GREY PLANETARY ATMOSPHERE MaRcos

We describe an efficient computational scheme for downward surface flux computations in a vertically inhomogeneous grey planetary atmosphere for different values of solar zenith angle. We start with the basic equations of a recently developed discrete ordinates spectral nodal method, and we derive suitable bidirectional functions whose diffuse components do not depend on the solar zenith angle. We then make use of these bidirectional functions to construct an efficient scheme for computing the downward surface fluxes in a given model atmosphere for a number of solar zenith angles. We illustrate the merit of the computational scheme described here with downward surface flux computations in a three-layer grey model atmosphere for four values of solar zenith angle, and we conclude this article with general remarks and directions for future work.


INTRODUCTION
We have recently developed a discrete ordinates �doR� methodology for efficiently and accurately solving radiative transfer problems on a stratified plane-parallel (multislab) grey planetary atmosphere �de aBReu, 2004a. in a continuous effort to further develop our methodology, we have recently derived discrete ordinates bidirectional functions for the boundary layers of a multislab model atmosphere �de aBReu, 2005b,c�.these bidirectional functions are response matrix-like relations that can be used to replace the boundary layers in multislab radiative transfer computations without degrading the numerical results.Moreover, they have shown to increase the efficiency of our DOR methodology in solving a set of multislab problems in which the optical properties of the boundary layers do not change from one problem to another in the set.
More recently, we have derived discrete ordinates bidirectional functions for an internal layer �de aBReu, 2006).As with our boundary layer functions, our internal layer functions are response matrix-like relations that can be used to replace without loss of accuracy an internal layer in multislab radiative transfer computations, and they have also shown to increase the efficiency of our DOR methodology in solving a set of multislab problems with an optically stationary internal layer.
In this article, we briefly discuss our discrete bidirectional functions and, in connection to our doR methodology, we describe a computational scheme for efficiently computing discrete ordinates approximations to the downward surface fluxes (CHANDRASEKHAR, 1950) in a vertically inhomogeneous grey model atmosphere for a number of solar zenith angles.Efficient and accurate downward surface flux computations are central to a number of relevant solar-atmospheric applications such as environmental remote sensing �tHoMas and staMnes, 1999; Liou, 2002�, solar radiation dosimetry (KRZYSCIN ET AL., 2003;BOLDEMANN ET AL., 2006), and solar power plant design (SEGAL AND EPSTEIN, 2003;FRICKER, 2004).
An outline of the remainder of this article follows.In Section 2, we introduce the governing equations of the atmospheric radiation problems considered in this article.in Section 3, we briefly discuss our discrete bidirectional functions with the help of adequate referencing, and we show how they can be used to construct an efficient scheme for downward surface flux computations in a grey model atmosphere for a number of solar zenith angles.In Section 4, we present numerical results for downward surface fluxes for a three-layer model atmosphere, and in Section 5 we close this article with concluding remarks and directions for further research.

THE GOVERNING EQUATIONS
Let us consider a vertically inhomogeneous grey planetary atmosphere illuminated with a monodirectional beam of solar radiation (CHANDRASEKHAR, 1950).We assume that the atmosphere overlies a dark surface of a terrestrial planet, and that it can be represented by a multislab domain consisting of R contiguous, disjoint, and uniform layers each so that layer 1 �one� is the uppermost layer and layer R is the lowermost one.We assume further that the frequency-integrated diffuse component of the radiation field in the shortwave can be accurately described by the frequency-integrated discrete ordinates scalar equations (DE ABREU, 2004a) where the integer N is the order of the angular quadrature, m =1 : N/2 holds for the downward directions µm > 0, m = N/2 +1 : N holds for the upward directions µm < 0, and the remaining notation in equations ( 1) is standard.Details about the uncollided component of the radiation field can be found elsewhere (THOMAS AND STAMNES, 1999;LIOU, 2002).We assume that changes in the refractive index are negligible between two adjacent layers, and so continuity conditions for the diffuse intensity at layer interfaces can be considered.
In addition, we use the discrete homogeneous boundary conditions �2� Equations ( 1) and ( 2) can be numerically solved with our new DOR methodology without optical truncation error.a detailed discussion on the advantages and limitations of this methodology can be found in earlier work (DE ABREU, 2004a,b, 2005a�.one of the advantages is that it provides us with the means of deriving nonstandard discrete bidirectional functions, which replace without loss of accuracy an atmospheric layer in radiative transfer computations, as we shall see in the next section.6) and ( 7) can be found in earlier work (DE ABREU, 2004a,b, 2005a�.

DISCRETE BIDIRECTIONAL FUNCTIONS AND A COMPUTATIONAL SCHEME
Equations ( 3) and ( 6) are the starting point for the derivation of our discrete bidirectional functions.upon substitution of equations ( 6) into equations (3) for an arbitrary layer r, and after half-range splitting and suitable reformulation, we get the discrete bidirectional functions: where definitions for H r , m , p and 0 m , r Q can be found in de aBReu �2006�.the discrete bidirectional functions �8� and �9� give us the emerging diffuse intensities at the bottom and top edges of layer r, respectively, in terms of the incident diffuse intensities at top and bottom and of corresponding first collided source terms.For r varying from two to R-1, these bidirectional functions play for each r the role of discontinuous conditions, and can advantageously be used to replace an internal layer in some multislab radiative transfer computations (DE ABREU, 2006).For r = R and in view of the discrete boundary conditions �2� at t r , the second summation terms in the bidirectional functions ( 8) and ( 9) vanish, and so we may write: The coefficients H R , m , p in the discrete bidirectional functions �10� and �11� play the role of discrete bidirectional transmittance and reflectance functions, respectively, for the lowermost layer of a multislab model atmosphere.As in the case of an internal layer, the discrete bidirectional functions �10� and (11) can advantageously be used to replace the lowermost layer in some multislab radiative transfer computations �de aBReu, 2005b).For r = 1 and in view of the discrete boundary conditions �2� at t 0 , the first summation terms in the bidirectional functions (8) and ( 9) vanish, and so we may write: The coefficients H 1 , m , p in the discrete bidirectional functions �12� and �13� play the role of discrete bidirectional reflectance and transmittance functions, respectively, for the uppermost layer of a multislab model atmosphere.as in the lowermost layer case, the discrete bidirectional functions ( 12) and �13� can advantageously be used to replace the uppermost layer in some multislab radiative transfer computations �de aBReu, 2005c�.
Let us now consider an application basic to radiative transfer in planetary atmospheres.Suppose that we wish to compute a discrete ordinates approximation to the downward surface flux in a vertically inhomogeneous grey atmosphere for a number of solar zenith angles.Since the coefficients H r,m,p depend only on the optical properties of the corresponding layers and on the quadrature points and weights used in the discrete ordinates scalar equations ( 1) and ( 2), we can do the following: we compute and store the coefficients H r,m,p for all layers of the multislab model atmosphere.We then replace the layers with the bidirectional functions (8-9), (10-11), or (12-13), where appropriate, and we compute the downward surface flux with our DOR methodology, as described in DE aBReu �2004a�, for each solar zenith angle.We do not have to recompute the coefficients H r,m,p each time the downward surface flux is computed with our DOR methodology for a new value of solar zenith angle.this integrated strategy prevents unnecessary (re)computation of available quantities (H r,m,p �, and thereby increases the efficiency of our DOR methodology for downward surface flux computations, as we shall illustrate next.

NUMERICAL RESULTS
to illustrate the computational scheme outlined in the previous section for efficient downward surface flux computations, we consider a three-layer model (R = 3) for a vertically inhomogeneous grey atmosphere: a conservative tropospheric air layer �r = 2, Dt 2 = 1, v 2 = 1� is surrounded by an overlying background stratospheric aerosol layer �r = 1, Dt 1 = 0.2, v 1 = 0.95� and by an underlying land haze layer �r = 3, Dt 3 = 0.3, v 3 = 0.9�.a greytone representation of this grey model atmosphere is shown in Figure 1.We assume that the aerosol layer consists of spherical (Lorentz-Mie) particles with an asymmetry factor of 0.9 (LIOU, 2002�.the Legendre scattering phase function data for the tropospheric layer �L 2 = 8� and for the haze layer �L 3 = 82) were extracted from a work of GARCIA AND SIEWERT (1985).The Legendre scattering data for the stratospheric aerosol layer (we set L 1 = 100 for simplicity� can be computed from the Legendre expansion of the Henyey-Greenstein phase function model for the aerosol particles �Liou, 2002�. In Table 1, we present N = 200 results for the total (uncollided plus diffuse) downward surface flux (F + � for four different values of m 0 , the cosine of the solar zenith angle. in Table 2, we show computer memory location and execution time for the four runs with (integrated computation) and without (independent computation) the scheme outlined in the last paragraph of section 3. We should note that the numerical results presented in this section were obtained executing our DOR computer program written in standard FORTRAN on an iBM-compatible Pc �1.4 GHz-clock intel Pentium 4 processor, 256 Mbytes of RAM) running on a GNU/Linux platform.
In order to ensure the quality of our results, we note that the numerical values in Table 1 all agree to within ±1 in the last figures given with corresponding results generated on a very fine-mesh grid by a (less efficient) discrete ordinates computer program based on an unconditionally stable discretization method -diamond difference -that is (still) widely used in radiation transport computations �LeWis and MiLLeR, 1993�.it is apparent that our computational scheme increased the efficiency of our DOR methodology for downward surface flux computations, since significant reductions in both computer memory (a 31.6%reduction) and execution time (a 22.1% reduction) were achieved.

CONCLUDING REMARKS
We have described in connection to our doR methodology an efficient scheme for computing the downward surface flux in a vertically inhomogeneous grey planetary atmosphere for different values of solar zenith angle.We have assumed that the planetary surface is dark, that the atmosphere can be represented by a multi-layered medium, and that the optical properties of the layers are known quantities.Then, we used the basic equations of our DOR methodology to derive discrete bidirectional functions for each layer.these functions play the role of response matrices to incident diffuse radiation and to first-collided sources.Since their diffuse components �H r,m,p ) do not depend on the solar zenith angle, we were able to use these bidirectional functions to construct an efficient scheme for computing the downward surface fluxes for an arbitrary number of solar zenith angles.We plan to extend this computational scheme to the more realistic case of a reflecting (rather than a dark) planetary surface, and we intend to report our findings on these matters in a forthcoming article.
We conclude by noting that the discrete bidirectional functions �8-13� can be incorporated into discrete ordinates methodologies other than our doR methodology, e.g. the methods described by staMnes et aL.�1988�, sieWeRt (2000, and BARICHELLO AND SIEWERT (2001) as follows: the discrete bidirectional functions (10-11) and/or (12-13), where appropriate, can be used to replace corresponding boundary layers at bottom and/or top, while the discrete bidirectional functions �8-9� can be used to replace corresponding internal layers, playing the role of fictitious discontinuous conditions �de aBReu, 2006� for the diffuse intensity of the radiation field.
We begin this section by setting forth the equations of our DOR methodology.These equations are two sets of optically discretized equations.The first set consists of the radiative heat balance equations�3�where Dtr L tr -tr-1 is the optical thickness of layer r; and i d r-1,m i d r,m are the diffuse intensities in direction m at top and bottom of layer r, respectively, �4� is the rth layer-average diffuse intensity in direction m, and: �5� is the rth layer-average first collided source for direction m.The quantity (tr-1-t0) in expression (5) is the relative optical depth of layer r.The second set consists of the spectral auxiliary equations �6� where �7� a detailed description of the numerical schemes employed for determining the coefficients q r,m,p

Figure 1
Figure 1 -a three-layer model atmosphere

Table 2 .
Computer memory and execution (CPU) time.