Abstract

The spherical harmonics discrete ordinate method for plane-parallel data assimilation (SHDOMPPDA) model is an unpolarized plane-parallel radiative transfer forward model, with corresponding tangent linear and adjoint models, suitable for use in assimilating cloudy sky visible and infrared radiances. It is derived from the spherical harmonics discrete ordinate method plane-parallel (SHDOMPP, also described in this article) version of the spherical harmonics discrete ordinate method (SHDOM) model for three-dimensional atmospheric radiative transfer. The inputs to the SHDOMPPDA forward model are profiles of pressure, temperature, water vapor, and mass mixing ratio and number concentration for a number of hydrometeor species. Hydrometeor optical properties, including detailed phase functions, are determined from lookup tables as a function of mass mean radius. The SHDOMPP and SHDOMPPDA algorithms and construction of the tangent-linear and adjoint models are described. The SHDOMPPDA forward model is validated against the Discrete Ordinate Radiative Transfer Model (DISORT) by comparing upwelling radiances in multiple directions from 100 cloud model columns at visible and midinfrared wavelengths. For this test in optically thick clouds the computational time for SHDOMPPDA is comparable to DISORT for visible reflection, and roughly 5 times faster for thermal emission. The tangent linear and adjoint models are validated by comparison to finite differencing of the forward model.

1. Introduction

Variational assimilation of visible and infrared radiances by numerical models in cloudy skies requires forward and adjoint radiative transfer models capable of handling scattering. When cloud properties are the target of the assimilation, visible and near-infrared satellite radiances should be considered because reflected solar radiation provides important information about cloud water path and particle size (e.g., Twomey and Cocks 1982). Due to heavy computational costs and the crude representation of clouds in numerical weather prediction models, cloudy visible and infrared satellite radiances have generally not been assimilated in the past. Recently, experiments assimilating Geostationary Operational Environmental Satellite (GOES) imager channels with a numerical cloud model have been performed (Vukicevic et al. 2004, 2006).

A few radiative transfer models that include scattering and are specifically designed for assimilation have been developed recently. Mostly these have been for microwave radiance assimilation (e.g., Bauer et al. 2005) and seek improvement upon the standard Eddington second approximation method. Weng and Liu (2003) modified their fully polarized vector discrete-ordinate radiative transfer code to calculate the radiance Jacobian or adjoint for solar and thermal sources. Greenwald et al. (2005) developed a very fast multistream successive order of scattering method designed for microwave radiative transfer and compared it with the Eddington method. Liou et al. (2005) introduce a four-stream polarized radiative transfer model intended for assimilation, though they did not present the adjoint model in that paper.

Forward and adjoint models for thermal infrared and solar radiative transfer have been developed for the Regional Atmospheric Modeling and Data Assimilation System (RAMDAS) GOES radiance assimilation system (Greenwald et al. 2002, 2004). The solar radiative transfer part of this observational operator was based on the spherical harmonics discrete-ordinate method (SHDOM; Evans 1998), which is a three-dimensional (3D) radiative transfer model. SHDOM was chosen for this 1D application because it can be faster than the standard discrete ordinate method. Owing to the complexities of the 3D model and its adaptive grid, a true adjoint was not created; instead, the adjoint was performed by application of the tangent linear model for each element of the input vector.

This use of SHDOM for 1D radiative transfer was the motivation to develop a plane-parallel version (SHDOMPP) and then to develop a data assimilation version (SHDOMPPDA) with corresponding tangent linear and adjoint models. Like its parent, SHDOMPP performs both solar and thermal radiative transfer using an iteration between spherical harmonics and discrete ordinate representations of the radiance field on an adaptive grid. The next section describes the SHDOMPP forward radiative model and development of the SHDOMPPDA forward, tangent linear, and adjoint models. Section 3 describes testing of the forward, tangent linear, and adjoint models, and shows examples of the adjoint sensitivities for cloud model fields.

2. The SHDOMPP and SHDOMPPDA algorithms

a. SHDOMPP

SHDOMPP calculates unpolarized radiative transfer in a plane-parallel medium for either collimated solar and/or thermal emission sources of radiation. The optical properties of the medium input to SHDOMPP are assumed to be uniform in each layer, which is different from the gridpoint medium specification in SHDOM. The specified optical properties for each layer are the optical depth, single-scattering albedo, and Legendre series of the phase function (delta-M scaling is always applied). The surface reflection options are Lambertian or Rahman–Pinty–Verstraete (RPV; Rahman et al. 1993) land surface bidirectional reflectance. The SHDOMPP calculation may be either monochromatic or spectrally integrated over molecular absorption lines with a k distribution. The output radiative quantities are upwelling and downwelling hemispheric fluxes and mean radiances at the layer boundaries and radiances at specified levels and directions.

SHDOMPP represents the radiation field using the source function from which the radiance field can be obtained by integrating the radiative transfer equation. The angular aspects of the source function are represented with spherical harmonics series, while the spatial aspects are represented with a discrete vertical grid. Only cosine terms are needed for the azimuthal part of the spherical harmonic functions because of the symmetry about the solar beam in the plane-parallel geometry. Like SHDOM, the vertical grid is adaptive; however, the spherical harmonic series truncation is fixed for all grid levels. The “base” grid (before new levels are adaptively added) has levels at the boundaries of the medium property layers and evenly spaced levels within the property layers so that the optical depth of each sublayer is less than a specified value (with a default setting of 1). Property layers with optical depths greater than 0.001 have at least two sublayers.

As with SHDOM the solution is achieved by iteration of the source function in four steps:

  1. the spherical harmonics source function is transformed to discrete ordinates;

  2. the source function is used in the integral form of the radiative transfer equation to obtain the radiance field in discrete ordinates;

  3. the radiance field is transformed to spherical harmonics;

  4. the source function is computed from the radiance field in spherical harmonics.

The same series acceleration technique as in SHDOM is used to improve the convergence of these iterations, which is important for optically thick media with little absorption. The iterations continue until the normalized rms change in the source function in one iteration (the solution criterion) is below a specified value (the solution accuracy).

The discrete ordinates are defined by double Gaussian quadrature in zenith angle and equal spacing in azimuth, but with a reduced number of azimuth angles near the nadir and zenith directions (as in SHDOM). The number of discrete angles is specified in zenith (Nμ) and azimuth (Nϕ for 0 to 2π), and the spherical harmonic truncation is set to L = Nμ − 1 and M = Nϕ/2 − 1. The transforms between spherical harmonics and discrete ordinates representations are performed as in SHDOM; however fast Fourier transforms in azimuth are not used in the SHDOMPPDA version to simplify the adjoint calculation. A substantial departure from SHDOM is that the source function is interpolated within each property layer by a cubic polynomial (quadratic at the layer boundaries) in the integration of the radiative transfer equation. This gives higher accuracy in the source function integration compared to the linear interpolation in SHDOM and allows fewer levels to be used for the same accuracy.

The adaptive grid is implemented using a different criterion from SHDOM. Sublayers are split in half if their splitting criterion exceeds the current splitting accuracy (which decreases as the solution iterations proceed). The splitting criterion for a sublayer is the rms difference in radiance over outgoing discrete ordinates between the original cubic source function interpolation integration and a linear interpolation integration. Unlike SHDOM, the splitting criterion is normalized by the rms radiance over all ordinates and levels (i.e., it does not depend on the incident solar flux).

The hemispheric fluxes and mean radiances are calculated with quadrature integration of the discrete ordinate radiances during the solution iterations. The radiances are calculated by integration of the source function using the Nakajima and Tanaka (1988) TMS method. SHDOMPP is distributed from the SHDOM Web site along with a program for preparing the optical properties file using Mie theory.

b. SHDOMPPDA

1) Forward model

SHDOMPPDA is a version of SHDOMPP designed to be useful for assimilating cloudy satellite radiances. The forward model of this code version includes the ability to calculate optical properties from the mixing ratio and number concentration of any number of hydrometeor species. The scattering properties of the particle species are stored in tables as a function of mass mean radius for a mass content of 1 g m−2. Mass mean radius is used instead of the more common effective radius because it can be calculated from the mixing ratio and number concentration without any assumptions about the size distribution or particle shape. The assumptions about particle size distribution and shape are still made in generating the scattering tables, but do not need to be input to SHDOMPPDA. The advantage of using scattering tables instead of single-scattering parameterizations (as in Greenwald et al. 2002) is that the high accuracy of Mie theory (for liquid clouds) is preserved and the computational cost is negligible. The SHDOMPPDA distribution package contains two programs for generating the scattering tables. One is a Mie program for gamma or lognormal size distributions of spherical particles. The other generates scattering tables for gamma size distributions of mixtures of six ice crystal shapes at wavelengths from 0.2 to 100 μm using a precomputed database (Yang et al. 2005) of optical properties for individual particle lengths from 2 to 9500 μm.

The input hydrometeor mixing ratios and concentrations, temperatures, and pressures are assumed to be defined at the boundaries of the radiative transfer layers. Thus the water path of each layer is derived from the average mixing ratio at the layer boundaries and the pressure change across the layer. Height is input, but is only used to select which hydrometeors are active according to a specified height range for each hydrometeor species. A subroutine interface is provided for the calculation of the channel-averaged molecular absorption profile from temperature, pressure, water vapor, and ozone profiles. Since the SHDOMPPDA package does not include routines for calculating molecular absorption, the molecular absorption is set to zero in the examples in section 3. For wavelengths shorter than 1.0 μm, molecular Rayleigh scattering is calculated using the pressure change across each layer and the wavelength. For efficiency, the channel radiances are currently calculated assuming monochromatic radiative transfer, which is a reasonable approximation for window channels appropriate for cloud assimilation. The forward radiative transfer routine calculates the top-of-atmosphere radiances (in reflectance units or brightness temperature) for a number of specified directions for a single input column of properties.

2) Tangent linear and adjoint models

The tangent linear is the gradient of the forward model in the direction of an input perturbation vector δxi:

 
formula

where x is the state vector of atmospheric and surface properties, yj(x) is the jth radiance computed by the forward model, and ∂yj/∂xi is the Jacobian of the forward model. The adjoint is the transpose of the Jacobian dotted into the radiance perturbation “forcing” δyj:

 
formula

The adjoint is a linearized expression of how much sensitivity each element of the input state vector has to giving a particular radiance pattern.

The tangent linear and adjoint models of SHDOMPPDA were developed using the Transformation of Algorithms in FORTRAN (TAF) software operated by FastOpt (Giering and Kaminski 1998). Adjoint compilers like TAF are not perfect, however, and one has to carefully monitor the produced code for efficiency and correctness. The philosophy behind building the tangent linear and adjoint models was to minimize hand coding by modifying the forward model so that TAF would generate an efficient adjoint.

First the phase function inputs to the SHDOMPP radiative transfer model were simplified. While the scattering tables represent the phase function with a Legendre series potentially thousands of terms long, only the first L + 1 terms are used in the SHDOMPP solution iterations. When calculating radiances from a solar source function, the whole Legendre series is summed, but only for one scattering angle per outgoing direction. Therefore the phase function inputs to forward radiative transfer were modified to be the Legendre series coefficients for l = 1, . . . , L + 1 and the phase function evaluated at each of the required scattering angles. This approach allowed the adjoint to be computed for a reasonable number of optical property inputs without resorting to the inaccurate assumption of a Henyey–Greenstein phase function (Greenwald et al. 2002). These phase function quantities along with the extinction and single-scattering albedo are computed for each layer from the water path and mean mass radius of each hydrometeor by interpolation in the scattering tables input from files.

To avoid operating TAF on the very complex dependences of the full SHDOMPP solution algorithm, the algorithm was divided into three parts. The first routine performs the adaptive grid portion of the SHDOMPP solution iterations, that is, until the layer splitting criterion is below the specified splitting accuracy. The output of this first routine is the adaptive grid structure and the spherical harmonic series of the source function at every adaptive grid level. The second routine performs the rest of the solution iterations on the now fixed grid, starting with the source function from the first routine. The output of the second routine are the desired radiances and the series acceleration step sizes for each iteration taken. The forward radiative transfer model of SHDOMPPDA consists of generating the optical properties and calling the first and second solution routines in sequence. The third routine performs the solution iterations on the fixed grid starting with the stored source function and using the precomputed series acceleration parameters. The tangent linear and adjoint models were developed with almost no hand coding by operating TAF on only the third routine. Thus, the forward model must be run (storing the adaptive grid structure, the intermediate source function, and the series acceleration parameters) before the tangent linear and adjoint models are run. This process of dividing the solution into three parts avoids complexity in the adjoint model and also improves the efficiency of the adjoint because only the ending iterations are performed in the adjoint. Good accuracy in the tangent linear and adjoint models requires more than just one iteration, however, so the solution accuracy should be smaller than needed for the forward model alone (e.g., 10−6 instead of 10−5).

The radiance adjoint should show sensitivity to cloud at clear levels as well as cloudy levels because perturbing a clear level with a small amount of cloud will indeed affect the outgoing radiance. If the forward model excludes levels with zero hydrometeor mixing ratio from the radiative transfer, then the adjoint will only show sensitivity where cloud already exists. To avoid this problem, very small values of mixing ratio and concentration are added to all layers in the process of converting to layer water path and mean mass radius. A priori information about where various hydrometeor species should exist may be provided by specifying the height range over which each species is active (the adjoint sensitivity will be zero outside this range).

3. Testing SHDOMPPDA

This section describes testing the forward, tangent linear, and adjoint model of SHDOMPPDA for correctness. In addition, the trade-off of radiance accuracy versus computational cost of the SHDOMPPDA forward model is compared with that of DISORT (Stamnes et al. 1988). These tests are performed using columns from an XZ slice of a 3D cloud model simulation. The cloud fields are from a Regional Atmospheric Modeling System (RAMS) simulation of the atmosphere centered on the Atmospheric Radiation Measurement program central facility in Oklahoma on 16 March 2001. The RAMS grid is 100 × 100 × 84 with 6-km horizontal spacing and a stretched vertical grid extending to above 16 km. The cloud fields used here are from the end of a 6-h run. This version of RAMS uses the new two-moment liquid cloud scheme (Saleeby and Cotton 2004), in which the mixing ratio and concentration of small (<40 μm) and large (40 to 80 μm) cloud droplet modes are prognostic. The ice hydrometeors are also predicted with a two-moment ice scheme (Meyers et al. 1997). The small cloud mode and the pristine ice hydrometeors are used in the radiative transfer calculations below. Figure 1 shows the mixing ratio and number concentration for the cloud water and pristine ice fields used in the radiative transfer simulations.

Fig. 1.

(top) Mixing ratio and (bottom) concentration of the liquid cloud and pristine ice hydrometeors from the RAMS slice used in the simulation. The liquid cloud is generally below 4.0-km altitude and the ice cloud mostly above 2.5 km, with the liquid cloud having the high concentration values.

Fig. 1.

(top) Mixing ratio and (bottom) concentration of the liquid cloud and pristine ice hydrometeors from the RAMS slice used in the simulation. The liquid cloud is generally below 4.0-km altitude and the ice cloud mostly above 2.5 km, with the liquid cloud having the high concentration values.

a. Forward model comparisons

Radiances are simulated at two wavelengths: solar reflectances at 0.63 μm and infrared brightness temperatures at 10.73 μm. The single-scattering properties are obtained from Mie theory for liquid cloud droplets and geometrical-optics and finite-difference time-domain methods (Yang et al. 2000, 2005) for ice crystals. Gamma particle size distributions are assumed, with ice crystal shapes following a size-dependent mixture of droxtals, solid and hollow columns, plates, bullet rosettes, and aggregates (mixture function from Baum et al. 2005, p. 1893). There is no molecular absorption, though Rayleigh scattering at 0.63 μm is included. The solar reflection simulation has a solar zenith angle of 45° and six viewing angles with zenith angles of 45°, 30°, and 15° and solar-viewing relative azimuth angles of 150° and 30°. The Lambertian surface albedo is 0.20. The infrared simulation assumes a fixed surface temperature of 290 K and surface emissivity of 0.98. Brightness temperatures are output for five viewing angles (60°, 45°, 30°, 15°, 0°). Both simulations allowed liquid cloud droplets from 0 to 4.5 km and ice crystals from 3.0 to 8.0 km. Figure 2 plots the solar reflectance and infrared brightness temperature computed with DISORT for the 100 columns in the RAMS slice. The liquid cloud is optically thick with reflectance above 0.75 for nearly three-quarters of the columns.

Fig. 2.

Solar reflectances at 0.63 μm and infrared brightness temperatures at 10.73 μm computed with DISORT. The reflectances are for a viewing zenith of 15° and relative azimuth of 150°, while the brightness temperatures are for a viewing angle of 45°.

Fig. 2.

Solar reflectances at 0.63 μm and infrared brightness temperatures at 10.73 μm computed with DISORT. The reflectances are for a viewing zenith of 15° and relative azimuth of 150°, while the brightness temperatures are for a viewing angle of 45°.

Radiances computed by SHDOMPPDA are compared to those from DISORT to validate the SHDOMPPDA forward model and to determine how the computational time and accuracy depend on the angular resolution of SHDOMPP and DISORT. The same procedure is used to convert the RAMS hydrometeor input to the optical properties input to SHDOMPP and DISORT. For the solar cases SHDOMPPDA is run with Nϕ = 2Nμ and a layer-splitting accuracy of 0.003. For the infrared cases SHDOMPPDA is run with Nϕ = 1 and a layer-splitting accuracy of 0.002. The solution accuracy is 10−5 in both cases. The DISORT reference runs have Nμ = 64 for solar and Nμ = 32 for infrared. The DISORT ACCUR parameter is set to zero so that all terms in the Fourier series in azimuth angle are used. Note that DISORT does not allow Nμ = 2.

Table 1 lists the computational times and the rms differences of the radiances from the reference runs as a function of Nμ for SHDOMPPDA and DISORT. The first conclusion from the comparisons is that SHDOMPPDA radiances agree with DISORT radiances, with rms differences as low as 0.0008 in reflectance and 0.007 K in brightness temperature. There is a steady convergence of SHDOMPPDA rms differences from the DISORT reference with increasing angular resolution, except perhaps at the highest Nμ where SHDOMPPDA accuracy is limited by the layer-splitting or solution accuracy.

Table 1.

Radiance accuracy and computing time as a function of the number of computational angles (Nμ) for DISORT and SHDOMPPDA. For the 0.63-μm case the absolute rms solar reflectance difference over six viewing angles and 100 columns (col) is listed. For the 10.73-μm case the rms brightness temperature difference (K) over five viewing angles and 100 columns is listed. DISORT is the reference for both DISORT and SHDOMPPDA. The computational times are on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler v.5.2.

Radiance accuracy and computing time as a function of the number of computational angles (Nμ) for DISORT and SHDOMPPDA. For the 0.63-μm case the absolute rms solar reflectance difference over six viewing angles and 100 columns (col) is listed. For the 10.73-μm case the rms brightness temperature difference (K) over five viewing angles and 100 columns is listed. DISORT is the reference for both DISORT and SHDOMPPDA. The computational times are on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler v.5.2.
Radiance accuracy and computing time as a function of the number of computational angles (Nμ) for DISORT and SHDOMPPDA. For the 0.63-μm case the absolute rms solar reflectance difference over six viewing angles and 100 columns (col) is listed. For the 10.73-μm case the rms brightness temperature difference (K) over five viewing angles and 100 columns is listed. DISORT is the reference for both DISORT and SHDOMPPDA. The computational times are on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler v.5.2.

Table 1 also shows that the computational time for SHDOMPPDA is comparable to or less than DISORT for the solar case and much less for the thermal infrared case. The solar example is close to the worst case in terms of computational time for SHDOMPPDA because it has a large amount of optically very thick cloud and is conservative scattering. The computational time for DISORT is independent of the optical thickness, but SHDOMPPDA requires many more computational layers and takes more iterations for high optical depth clouds. The quick iteration convergence for the highly absorbing clouds explains the timing advantage of SHDOMPPDA for the infrared case. What is most important for practical use of a model is the trade-off in computational time versus radiance accuracy. Figure 3 shows this trade-off graphically for SHDOMPPDA and DISORT. For reasonable reflectance accuracies (>0.001) SHDOMPPDA offers comparable performance to DISORT for solar reflectance from optically thick clouds. In the thermal infrared SHDOMPPDA is five times faster than DISORT for 0.1-K rms accuracy. Further optimization of the SHDOMPPDA performance might be made by adjusting the layer splitting accuracy, which was not done here.

Fig. 3.

Computational time as a function of radiance accuracy for (top) solar and (bottom) infrared radiative transfer comparing SHDOMPPDA and DISORT.

Fig. 3.

Computational time as a function of radiance accuracy for (top) solar and (bottom) infrared radiative transfer comparing SHDOMPPDA and DISORT.

b. Tangent linear and adjoint testing

The SHDOMPPDA tangent linear and adjoint models are validated by comparison to finite differencing of the forward model. This testing has been performed for various columns, but is illustrated here with column number 51 or x = 300 km in the RAMS slice. The radiative transfer parameters are the same as before, but the SHDOMPPDA parameters are Nμ = 6, Nϕ = 12 for the 0.63-μm simulation and Nμ = 4, Nϕ = 1 for the 10.73-μm simulation. The layer splitting accuracy is 0.001, and the solution accuracy is 10−7. The SHDOMPDDA adjoint is calculated with a unit “forcing” in reflectance or brightness temperature for a single radiance direction (θ = 30°, ϕ = 120°). The adjoint is also approximately calculated with one-sided finite differencing by perturbing one element at a time (e.g., the pristine ice concentration at one height level) by a small amount δx, calling the forward model, and subtracting the original reflectance or brightness temperature from the perturbed one to obtain δy. The finite difference adjoint is δy/δxi for each element i perturbed in turn. The tangent linear is compared with the adjoint by successively setting only one element of the input perturbation vector to unity (δxi = 1) and calling the SHDOMPPDA tangent linear model to calculate the exact radiance gradient for that element of the state vector (δy/δxi).

Figures 4 –6 show the comparison of the SHDOMPPDA adjoint with finite differencing for this one column. The SHDOMPPDA adjoint agrees with the finite different adjoint usually to well within 1% except where the adjoint is near zero. For example the finite difference is a very poor approximation for 10.73 μm below the liquid cloud top because there is virtually no sensitivity of the outgoing midinfrared radiance to changes in the lower portion of the cloud. The comparison of the tangent linear to the adjoint is not shown because the differences are almost always less than one part in 105.

Fig. 4.

Comparison of the SHDOMPPDA 0.63-μm reflectance adjoint with finite differencing for (top) mixing ratio and (bottom) number concentration. The left side of each panel shows the SHDOMPPDA adjoint, while the right side plots the ratio of the adjoint to the finite difference.

Fig. 4.

Comparison of the SHDOMPPDA 0.63-μm reflectance adjoint with finite differencing for (top) mixing ratio and (bottom) number concentration. The left side of each panel shows the SHDOMPPDA adjoint, while the right side plots the ratio of the adjoint to the finite difference.

Fig. 6.

Comparison of the SHDOMPPDA 10.73-μm brightness temperature adjoint with finite differencing for the temperature profile.

Fig. 6.

Comparison of the SHDOMPPDA 10.73-μm brightness temperature adjoint with finite differencing for the temperature profile.

For practical use of the SHDOMPPDA adjoint in variational assimilation it is important that the computing time not be excessive. Table 2 lists the computational times to run the SHDOMPPDA forward and adjoint models per column averaged over the 100 columns in the RAMS slice. As noted in section 2b, the forward model and adjoint are run with a smaller solution accuracy (10−6) for good accuracy of the adjoint. The smaller solution accuracy results in more iterations on the fixed grid and is part of the reason that the solar reflection computational times are four to six times longer than those in Table 1. The thermal emission calculation has few iterations and the forward plus adjoint times are only about twice as long as the forward times listed in Table 1. The results in Table 2 illustrate that assimilation of solar reflectances will be computationally challenging. Depending on the particular input parameters it is likely that the computational times could be reduced somewhat by using fewer azimuthal discrete ordinates and higher layer splitting accuracy. While the accuracy results definitely depend on the solar-viewing geometry, wavelength, and cloud situation, it appears from Table 1 that Nμ = 6 may be adequate for solar reflection calculations and that Nμ = 4 should be sufficient for midinfrared radiative transfer.

Table 2.

Computational time to run the forward and adjoint SHDOMPPDA models on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler.

Computational time to run the forward and adjoint SHDOMPPDA models on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler.
Computational time to run the forward and adjoint SHDOMPPDA models on a 3.2-GHz Intel Xeon processor with Portland Group FORTRAN compiler.

Figures 7 and 8 show the 0.63-μm reflectance cloud water adjoints and the 10.73-μm brightness temperature pristine ice adjoints, respectively. The reflectance sensitivity to cloud mixing ratio is positive, indicating that more cloud mass causes more reflectance. The reflectance sensitivity is highest where the background is dark, such as the clear region on the right side of the domain, and lowest in the liquid cloud layer. The increase in reflectance sensitivity to mixing ratio with height in the clear region is due to the larger spacing in pressure between levels, which results in a given mixing ratio resulting in more liquid water path in a layer. The reflectance sensitivity to cloud concentration is also positive because more cloud droplets for a fixed mixing ratio results in smaller droplets and thus higher optical depth. The reflectance sensitivity to concentration is highest where the cloud droplet concentration is small. The 10.73-μm brightness temperature adjoints to ice cloud mixing ratio and concentration are negative because more optical depth decreases the infrared brightness temperature. The larger brightness temperature sensitivity to mixing ratio at high altitudes is probably due to lower temperatures giving more contrast with the background brightness temperature. The brightness temperature sensitivity is near zero below where the cloud optical depth is high.

Fig. 7.

The 0.63-μm reflectance adjoint field for liquid cloud mixing ratio and number concentration computed by SHDOMPPDA.

Fig. 7.

The 0.63-μm reflectance adjoint field for liquid cloud mixing ratio and number concentration computed by SHDOMPPDA.

Fig. 8.

The 10.73-μm brightness temperature adjoint field for pristine ice cloud mixing ratio and number concentration computed by SHDOMPPDA.

Fig. 8.

The 10.73-μm brightness temperature adjoint field for pristine ice cloud mixing ratio and number concentration computed by SHDOMPPDA.

4. Summary

This article describes two closely related unpolarized plane-parallel radiative transfer models for atmospheres with particle scattering: SHDOMPP and SHDOMPPDA. SHDOMPP is a plane-parallel version of the 3D radiative transfer model SHDOM. SHDOMPP reads optical properties from a file and calculates fluxes and radiances for solar and/or thermal radiation sources. SHDOMPPDA is based on SHDOMPP, but includes forward, tangent linear, and adjoint models used for data assimilation. SHDOMPPDA inputs mixing ratio and number concentration profiles for any number of hydrometeor species and calculates upwelling radiances.

SHDOMPP represents the radiation field with the source function expressed as spherical harmonics series (with a fixed truncation) at discrete grid levels. The optical properties input to SHDOMPP are assumed uniform in each layer, like DISORT, not linearly interpolated between grid points as in SHDOM. The solution is achieved by iterations in which the source function is transformed from spherical harmonics to discrete ordinates, the source function is integrated in the radiative transfer equation to obtain the new radiance field, the radiance field is transformed from discrete ordinates to spherical harmonics, and the radiance field is then used to compute the source function in spherical harmonics. The integration along discrete ordinates is more accurate than in SHDOM because the source function is interpolated with a cubic (instead of linear) polynomial between the grid levels. The vertical grid is adaptive so that, as the solution iterations proceed, new levels are added where needed to improve the accuracy of the integration.

The SHDOMPPDA forward model includes optical property calculations using scattering tables as a function of mass mean radius. Mass mean radius is chosen instead of effective radius because it can be directly calculated from the input hydrometeor mass mixing ratio and number concentration without further assumptions about the particle size distribution and shape. An interface to user-provided molecular absorption routines is included. The SHDOMPPDA tangent linear and adjoint models were built using the Transformation of Algorithms in FORTRAN (TAF) software operating on the forward model. To facilitate an efficient adjoint, a number of changes were made in the SHDOMPP forward model. These include dividing the SHDOMPP solution method into three parts and having TAF operate only on the third routine that finishes the solution iterations on a fixed grid using precomputed iteration acceleration parameters.

The SHDOMPPDA forward model was tested by comparison with DISORT for optical properties calculated from a 100 × 84 slice of a 3D field of cloud droplets and ice crystal simulated with RAMS. Radiative transfer simulations were performed for solar reflection at 0.63 μm and thermal emission at 10.73 μm to calculate the upwelling radiance in a number of directions. Radiances from SHDOMPPDA and DISORT were compared over the five or six directions and 100 columns as a function of the angular resolution of the models, Nμ. The agreement in rms differences is as low as 0.0008 in 0.63-μm reflectance and 0.007 K in 10.73-μm brightness temperature. SHDOMPPDA is comparable to DISORT in computational time for the solar reflection calculation in this mostly cloudy, optically thick slice (which favors DISORT). SHDOMPPDA is considerably faster than DISORT (five times for 0.1-K accuracy) for the thermal emission calculation.

The SHDOMPPDA tangent linear and adjoint were tested by comparison to finite differencing of the forward model and to each other for all relevant inputs. The adjoint and the finite differencing agree to within a few tenths of a percent except where the adjoint sensitivity is very small. The adjoint and the tangent linear agree to within one part in 105 for almost all elements of the input temperature, hydrometeor mixing ratio, and concentration profiles. The computational time for the forward-plus-adjoint models is two to six times that of the forward model alone.

One advantage of SHDOMPPDA is that the user can decide the trade-off point between the desired computational speed and the appropriate accuracy by adjusting the number of discrete ordinates, the layer splitting accuracy, and the solution accuracy parameters. While SHDOMPPDA may not be the fastest possible radiative transfer model in the midinfrared (though it could be competitive), it has the ability to handle solar reflection as well as thermal emission wavelengths. For assimilation of both solar and midinfrared wavelength, such as GOES imager radiances, the solar reflection wavelengths will dominate the computational costs of the observational operator. Thus there is no penalty to using SHDOMPPDA for all wavelengths.

SHDOMPP and SHDOMPPDA are freely distributed from a Web page linked to the SHDOM Web site. Both models are written in FORTRAN 90. The SHDOMPPDA distribution package contains example scripts and datafiles for using the model and includes programs for computing the particle-scattering tables.

Fig. 5.

Comparison of the SHDOMPPDA 10.73-μm brightness temperature adjoint with finite differencing for (top) mixing ratio and (bottom) number concentration. The left side of each panel shows the SHDOMPPDA adjoint, while the right side plots the ratio of the adjoint to the finite difference.

Fig. 5.

Comparison of the SHDOMPPDA 10.73-μm brightness temperature adjoint with finite differencing for (top) mixing ratio and (bottom) number concentration. The left side of each panel shows the SHDOMPPDA adjoint, while the right side plots the ratio of the adjoint to the finite difference.

Acknowledgments

The author is grateful to Tomi Vukicevic for involving him in the RAMDAS project and providing encouragement and financial support through a Joint Center for Satellite Data Assimilation project (NOAA Award NA17RJ1228 to Colorado State University). Laura Fowler kindly providing the the RAMS fields used in the radiative transfer simulations. Istvan Laszlo helped choose parameters for high-accuracy DISORT computations. Financial support was also provided by the Environmental Sciences Division of the U.S. Department of Energy (under Grant DE-A1005-90ER61069 to Warren Wiscombe at NASA Goddard Space Flight Center) as part of the ARM program.

REFERENCES

REFERENCES
Bauer
,
P.
,
E.
Moreau
, and
S.
Di Michele
,
2005
:
Hydrometeor retrieval accuracy using microwave window and sounding channel observations.
J. Appl. Meteor.
,
44
,
1016
1032
.
Baum
,
B. A.
,
A. J.
Heymsfield
,
P.
Yang
, and
S. T.
Bedka
,
2005
:
Bulk scattering properties for the remote sensing of ice clouds. Part I: Microphysical data and models.
J. Appl. Meteor.
,
44
,
1885
1895
.
Evans
,
K. F.
,
1998
:
The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer.
J. Atmos. Sci.
,
55
,
429
446
.
Giering
,
R.
, and
T.
Kaminski
,
1998
:
Recipes for adjoint code construction.
ACM Trans. Math. Software
,
24
,
437
474
.
Greenwald
,
T. J.
,
R.
Hertenstein
, and
T.
Vukicevic
,
2002
:
An all-weather observational operator for radiance data assimilation with mesoscale forecast models.
Mon. Wea. Rev.
,
130
,
1882
1897
.
Greenwald
,
T. J.
,
T.
Vukicevic
,
L. D.
Grasso
, and
T. H.
Vonder Haar
,
2004
:
Adjoint sensitivity analysis of an observational operator for visible and infrared cloudy-sky radiance assimilation.
Quart. J. Roy. Meteor. Soc.
,
130
,
685
705
.
Greenwald
,
T.
,
R.
Bennartz
,
C.
O’Dell
, and
A.
Heidinger
,
2005
:
Fast computation of microwave radiances for data assimilation using the “successive order of scattering” method.
J. Appl. Meteor.
,
44
,
960
966
.
Liou
,
K. N.
,
S. C.
Ou
,
Y.
Takano
, and
Q.
Liu
,
2005
:
A polarized delta-four-stream approximation for infrared and microwave radiative transfer: Part I.
J. Atmos. Sci.
,
62
,
2542
2554
.
Meyers
,
M. P.
,
R. L.
Walko
,
J. Y.
Harrington
, and
W. R.
Cotton
,
1997
:
New RAMS cloud microphysics parameterization. Part II: The two-moment scheme.
Atmos. Res.
,
45
,
3
39
.
Nakajima
,
T.
, and
M.
Tanaka
,
1988
:
Algorithms for radiative intensity calculations in moderately thick atmospheres using a truncation approximation.
J. Quant. Spectrosc. Radiat. Transfer
,
40
,
51
69
.
Rahman
,
H.
,
B.
Pinty
, and
M. M.
Verstraete
,
1993
:
Coupled surface atmosphere reflectance (CSAR) model. 2. Semiempirical surface model usable with NOAA Advanced Very High Resolution Radiometer data.
J. Geophys. Res.
,
98
,
20791
20802
.
Saleeby
,
S. M.
, and
W. R.
Cotton
,
2004
:
A large-droplet mode and prognostic number concentration of cloud droplets in the Colorado State University Regional Atmospheric Modeling System (RAMS). Part I: Module descriptions and supercell test simulations.
J. Appl. Meteor.
,
43
,
182
195
.
Stamnes
,
K.
,
S.
Tsay
,
W.
Wiscombe
, and
K.
Jayaweera
,
1988
:
Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media.
Appl. Opt.
,
27
,
2502
2509
.
Twomey
,
S.
, and
T.
Cocks
,
1982
:
Spectral reflectance of clouds in the near-infrared—Comparison of measurements and calculations.
Meteor. Soc. Japan
,
60
,
583
592
.
Vukicevic
,
T.
,
T. J.
Greenwald
,
M.
Zupanski
,
D.
Zupanski
,
T.
Vonder Haar
, and
A. S.
Jones
,
2004
:
Mesoscale cloud state estimation from visible and infrared satellite radiances.
Mon. Wea. Rev.
,
132
,
3066
3077
.
Vukicevic
,
T.
,
M.
Sengupta
,
A. S.
Jones
, and
T.
Vonder Haar
,
2006
:
Cloud-resolving satellite data assimilation: Information content of IR window observations and uncertainties in estimation.
J. Atmos. Sci.
,
63
,
901
919
.
Weng
,
F.
, and
Q.
Liu
,
2003
:
Satellite data assimilation in numerical weather prediction models. Part I: Forward radiative transfer and Jacobian modeling in cloudy atmospheres.
J. Atmos. Sci.
,
60
,
2633
2646
.
Yang
,
P.
,
K. N.
Liou
,
K.
Wyser
, and
D.
Mitchell
,
2000
:
Parameterization of the scattering and absorption properties of individual ice crystals.
J. Geophys. Res.
,
105
,
4699
4718
.
Yang
,
P.
,
H.
Wei
,
H-L.
Huang
,
B. A.
Baum
,
Y. X.
Hu
,
G. W.
Kattawar
,
M. I.
Mischenko
, and
Q.
Fu
,
2005
:
Scattering and absorption property database for nonspherical ice particles in the near- through far-infrared spectral region.
Appl. Opt.
,
44
,
5512
5523
.

Footnotes

Corresponding author address: Dr. K. Franklin Evans, Department of Atmospheric and Oceanic Sciences, University of Colorado, 311 UCB, Boulder, CO 80309-0311. Email: evans@nit.colorado.edu

This article included in the Assimilation of Satellite Cloud and Precipitation Observations special collection.