**Colorado State University 3-D Monte Carlo Model**

Philip T. Partain and Graeme L. Stephens

Department of Atmospheric Science, Colorado State University

Fort Collins, Colorado

**1. The Origins of the Model**

The model we employed for this intercomparison is a modified version of
a forward Monte Carlo scheme used for the calculation of high-spectral
resolution, domain average shortwave fluxes in three-dimensional
scattering atmospheres. The complete model is separated into two parts:
the Monte Carlo calculation and an equivalence theorem calculation. The
Monte Carlo portion computes domain average fluxes and path length
information for a model atmosphere without gaseous absorption. The
equivalence theorem portion then uses the fluxes and path length
information to incorporate gaseous absorption and produce the final
fluxes. Further information about the equivalence theorem method can be
found in Irvine (1964), Van de Hulst (1980), and Partain *et* *al*.
(1999).

The original model was constructed to run on the Cray T3E supercomputers
located at the Lawrence Berkley National Laboratory (LBNL) and the
Geophysical Fluid Dynamics Laboratory (GFDL). The Cray T3E computers are
multi-processor machines that enable portions of programs to run in
parallel, reducing overall run time. In this application, 1 to 128
processors are employed which introduce photons into an associated grid
box on the top of the model atmosphere. Each processor then tracks the
photons, one at a time, through the atmosphere until they exit the
bottom (unless the surface albedo is non-zero) or top of the atmosphere
(boundaries are periodic in the x and y dimensions). Meanwhile, photons
are counted as they pass through levels at which fluxes are to be
calculated and their corresponding path lengths are stored.

Other details of the Monte Carlo portion of the model include the
following. The distribution of photons on the top of the atmosphere is
determined with a Halton sequence, whereas random numbers used for
scattering purposes are generated with the ran3 subroutine from
Numerical Recipes. The model can accommodate single and double
Henyey-Greenstein phase functions as well as the input of an external
phase function. The model can also accommodate Rayleigh extinction.

**2. Changes to the model for I3RC**

In order to participate in the Intercomparison of Three-Dimensional
Radiation Codes (I3RC), several modifications to the original model were
required. First of all, this phase of I3RC does not consider gaseous
absorption, so the equivalence theorem routine and all associated path
length variables were eliminated. Second, I3RC requires pixel fluxes in
addition to domain average fluxes. The ability to track photons crossing
flux measurement levels within individual pixels was added. In addition,
the capability to calculate nadir, zenith, and oblique radiances were
added. The method for doing so followed that in Marshak *et al*. (1995),
Appendix A.2:

$E_{j}(\Omega \ast ) = \sum_{k = 1}^{N}{\omega_{o}^{k}P(\Omega_{k},\Omega \ast )\text{exp}\lbrack - \sigma_{j}\tau} \ast \rbrack$

where E~j~(\*) is the contribution to the radiance at grid point j
exiting the cloud in the direction \*, ~o~ is the single-scattering
albedo of the cloud, P(~k~,\*) is the value of the phase function when
the photon scatters at an angle between ~k~ and \*, and \* is the
optical path traversed as the ray is traced out of the cloud in the
direction \*. This entailed the creation of code to retrieve the proper
phase function value and integrate optical depth from each scattering
point in the direction required for each radiance.

Rather than using the Cray T3E systems, the new model is run on a
cluster of 9 workstations at Colorado State University. Each workstation
possesses a Pentium II 450 MHz processor. The functional difference
between this cluster and the T3E systems is in the decreased speed of
communication between processors. However, this is not an issue for the
Monte Carlo model because processor \"cross-talk\" is kept at a minimum.
On the other hand, the equivalence theorem routine in the original model
requires much more information exchange between processors.

Perhaps the largest source of uncertainty is the use of a fixed number
of photons for each case rather than tracking error and terminating an
experiment when the error levels are below a specified amount. For
domain average fluxes, errors are generally smaller than those for pixel
fluxes, especially when large numbers of pixels are used. Therefore, the
original model did not have error calculation built in. Due to time
constraints, error calculation routines were not added to the new model.

**3. Experiments**

Three different cloud fields were used in this phase of I3RC including a
two-dimensional step cloud, a two-dimensional cloud retrieved from radar
observations (MMCR), and a three-dimensional cloud retrieved from
Landsat observations.

# A. step cloud

This cloud field consists of 32 cloud cells in the x-direction and 1
cell in the z-direction. The first 16 cells have an optical depth of 2
and the last 16 cells have an optical depth of 18. For this field
9,000,000 photons were simulated resulting in an average of 281,250
photons injected into each cell. Resulting domain average albedo,
transmittance and nadir reflectivity values along with associated
standard deviations agree well with the mean of the results from the
other participants in this intercomparison. Although this is not the
same as comparing absolute accuracy against a single benchmark, it does
provide confidence that this model is very nearly getting the same
answer as most of the other models in this study. The results for the
step cloud and other cases can be seen on the I3RC website
(http://i3rc.gsfc.nasa.gov).

# B. MMCR

This cloud field consists of 640 cloud cells in the x-direction and 54
cells in the z-direction. Here, 27,000,000 photons were simulated
resulting in 42,188 photons being introduced into each cell at the top
of the atmosphere. We would expect errors to increase for pixel-level
quantities over those seen for the step cloud case. As suspected, domain
averaged results for reflectance, transmittance and absorptance agree
fairly well with the other Monte Carlo models. The pixel-level
comparisons, however, show more noise than they did for the step cloud,
although agreement is still very good with the other models.

# C. Landsat

The Landsat cloud field consists of 128 cells in the x- and
y-directions. Cloud top heights were provided for the third dimension.
For use in this model, the cloud top heights were binned into a
100-meter vertical grid resulting in 24 cells in the z-direction.
90,000,000 photons were used, translating to 5,493 photons per cell at
the top of the atmosphere. Errors for pixel-level quantities using this
small number of photons per pixel were expected to be larger than for
the MMCR case. As suspected, the 2-D images of the radiative quantities
are fairly noisy compared to Monte Carlo codes that used more photons
and as compared to the analytic solution methods. However, the noise is
certainly not unacceptable. The level of agreement for the domain
averaged quantities is very good as it was for the previous cases.

**4. Future Considerations**

The modified Monte Carlo model used for this intercomparison is not
ideally suited for computing pixel fluxes. No attempts to accelerate the
scheme with analytical techniques have been made. However, this
intercomparison has been useful to test this model\'s ability to
calculate domain average quantities and show its capabilities concerning
pixel-level fluxes. The basic schemes involving scattering and photon
distance calculations are no different than those in the original model.
Up to this point, only plane-parallel atmosphere results had been
verified. We are very pleased with the overall performance of this code
in comparison with the other participants.

The degree of involvement for this model in phase 2 will be determined
at a future date, depending upon the experiments chosen. At the very
least, this model along with the equivalence theorem scheme can
contribute domain average broadband results. Broadband results for pixel
quantities, aside from the large error in pixel fluxes, may not be
possible as a result of the large amount of path length information that
is required. Options for experiments in phase 2 might include complex
layered clouds with different phase function properties (ice and water
clouds in the same atmosphere) and Rayleigh scattering in addition to
gaseous absorption.

**5. References**

Irvine, W. M., 1964: The formation of absorption bands and the
distribution of photon optical paths in a scattering atmosphere. *B. A.
N.*, **17**, 266-279.

Marshak, A., A. Davis, W. Wiscombe, and R. Cahalan, 1995: Radiative
smoothing in fractal clouds. *J. Geophys. Res.*, **100**, 26,247-26,261.

Partain, P. T., A. K. Heidinger, and G. L. Stephens, 2000: High
resolution atmospheric radiative transfer: Application of the
equivalence theorem. *J. Geophys. Res.*, **105**, pp. 2163-2177.

Van de Hulst, H. C., 1980: Chap. 17, Multiple Light Scattering; Tables,
Formulas, and Applications. Academic Press.