***********************************************************************
Documentation File for 'DISORT', A Discrete Ordinates Radiative
Transfer Fortran90 Program for a MultiLayered PlaneParallel Medium
(VERSION 2.0)
***********************************************************************
NOTE: If you have not received DISORT directly from its main
anonymous ftp site
ftp://climate.gsfc.nasa.gov/pub/wiscombe/Multiple_Scatt/
you should check there for the latest version, and any updates.
FORTRAN90: This version and v1.3 of DISORT require a Fortran90
compiler. So far, the only Fortran90 actually used is in RDI1MACH.f,
since these routines in their f77 form were by far the most unpopular
among DISORT users; their f90 implementation, by contrast, is clean and
simple. The rest of DISORT will gradually migrate to Fortran90, for
reasons given in RDI1MACH.readme and in the DISORT Report, Appendix B.
Note that any f77 code can be compiled by an f90 compiler, as long as
the filename has a .f extension indicating fixed source form.
Those without access to an f90 compiler can either
(a) use version 1.2; or
(b) use this version but substitute R1MACH.f and D1MACH.f from version 1.2.
Send us email if you find that, after a reasonable effort, you cannot
gain access to an f90 compiler; this will help us understand the extent
to which f90 compilers are available (or not) among our user base.

CODE HISTORY:
version 1.0, 1988 (Fall): built upon code foundation developed by Sichee
Tsay in his thesis under Knut Stamnes, largely during an intense fourweek
visit by Wiscombe to the University of Alaska, Fairbanks;
version 1.1 1993 (Jan): small bug fixes since 1988, mostly in ASYMTX
and IBCND=1 special case; NCUT option no longer applied if thermal
emission present; zero optical depth layers handled better; other
mostly cosmetic changes (printing formats, etc.)
version 1.2, 1997 (Feb): put under RCS control; many small cosmetic
changes; ALBTRN and SOLVE1 reorganized; fewer significant digits
printed in test problems, in order to reduce trivial differences
when comparing two outputs; in test problems, SIGNIFICANTLY NONUNIT
RATIOS only counted when flux or intensity above noise level;
SQT calculation moved up to DISORT;
version 1.3, 2000 (Mar): LAPACK used instead of LINPACK to do linear
algebra (but not for eigenvalue problem yet); first Fortran90 version;
R1MACH, D1MACH use Fortran90 intrinsic functions;
(Many thanks to Robert Pincus for the conversion to LAPACK)
version 2.0, 2000 (Mar): Introduced NakajimaTanaka TMS/IMS intensity
correction algorithm; a real surface BRDF is implemented.
The new features introduced in version 2 resulted in the following
major changes to the input/output arguments of DISORT:
(a) the full phase function Legendre expansion (all the significant
Legendre coefficients) must be provided, not just the Legendre
coefficients from 0 to 2N as in v1.x;
(b) the DELTAM option is no longer visible to the user and is always
turned on. It can still be turned off internally for testing;
however, inexperienced users are advised against turning it
off, because that will almost always result in worse results;
DELTAM is not needed in cases of weakly varying phase function,
e.g. Rayleigh scattering, but it does no harm in such cases either;
(c) azimuthallyaveraged intensity is no longer returned as an output;
(d) the mean intensities are not intensitycorrected; they are already
very accurate for small N using the usual deltaM approximation
and benefit little from the corrections; however, it is a technical
inconsistency since they are derived from the uncorrected
azimuthallyaveraged intensity not the corrected intensities;
(e) bidirectional reflectance as a function of polar (zenith) and
relative azimuth angles must be supplied in function BDREF
instead of coefficients in Legendrepolynomial expansion as in v1.x.
NOTE: Even though there have been many small changes to DISORT,
there have been no significant bug reports for many years. The
bugs found have tended to be obscure and unlikely to affect the vast
majority of users in any significant way.

This program package consists of the following files (in addition
to the file you are reading, DISORT.doc):
DISORT.f
The user entry point (subroutine DISORT), and most of the routines
it calls. All routines are in single precision except for the
eigenvalue routine ASYMTX and the quadratureangles routine QGAUSN.
DISORTsp.f
Like DISORT.f, but EVERYTHING is in single precision. Single
precision may be adequate for many applications, but this variant
is also necessary if you want to use the autodoubling option now
available on many compilers (if you autodouble, however, you need to
change all instances of 'R1MACH' to 'D1MACH'). This is also the
version to use on Cray and other highaccuracy machines.
BDREF.f
A default FUNCTION subprogram to supply surface bidirectional
reflectivity. The user must replace this with their own version
when if a nonLambertian lower boundary reflectivity is desired.
DISOTEST.f
A program for checking DISORT on a comprehensive set of test problems.
Be sure to run this before using DISORT.
DISOTEST.doc
Documentation for DISOTEST.f
DISOTEST.out
Output from running DISOTEST.f in normal mixed singledouble precision
on an IEEEarithmetic Unix workstation. Filename may have a
version number included.
RDI1MACH.f
Fortran90 Routines for returning machine constants, from netlib author
Eric Grosse; also see file RDI1MACH.readme.
LINPAK.f : Routines used by DISORT from the linearequationsolving
publicdomain packages LINPACK and BLAS (available
at most computer sites). Included only to make package
selfcontained. Slightly modified from originals
by upgrading to FORTRAN77. (The user is encouraged
to employ the local LINPACK in place of LINPAK.f,
because many computers have optimized their LINPACKs,
as for example CRAY computers with their
supervectorspeed LINPACK in 'libsci'.)
The entire LINPACK can be obtained from NetLib,
http://www.netlib.org/.
ErrPack.f
Errorhandling routines. Standard across almost all programs on
the wiscombe anonymous ftp site.
rundis
Example Cshell script to run DISORT. Shows how to put the various
pieces together, and how to handle single, double, and mixed single
double precision runs in an automated way.

It is HIGHLY recommended that you use the code in DISOTEST.f as a
template for creating your own CALLs to DISORT, rather than starting
from scratch. Then to compile DISORT, replace DISOTEST.f with your
calling program and add BDREF.f in the statements above. For problems
with Lambertian lower boundary you can use the stubfile BDREF.f
supplied with the code; in this case BDREF.f is not used, but it is
needed for compiling. For problems with bidirectionally reflecting
surfaces you must supply your own BDREF.f instead of the one
distributed with the code.

All the DISORT code has passed tests by 'flint', 'ftnchek', and
'nag_pfort'  tools that test the semantics of Fortran programs at a
higher level than compilers do. These tools were helpful in exposing
subtle bugs in DISORT in the early days, and remain useful for ensuring
that new bugs are not introduced. 'ftnchek' is free and available from
http://www.netlib.org; the others are commercial products (see the Fortran
Market, http://www.fortran.com/) but well worth the investment.
Remarks on Computer Precision

DISORT has certain intrinsic limitations because of computer
precision. These limitations are related to ordinary computer
"roundoff error" and have nothing to do with usercontrollable
(through the number of streams NSTR) "truncation error". DISORT
is free of the *catastrophic* growth of roundoff error that
plagued pre1980 discrete ordinate programs, but certain parts
of the calculation (the eigenvalue/vector and Gauss quadrature
rule computations, and to a much lesser extent the linear
equation solving computations) are just inherently more sensitive
to roundoff than the rest of DISORT, because they involve so many
arithmetic operations.
The reason DISORT.f does the eigenvalue/vector and Gauss
quadrature rule computations in double precision is that DISORT
was originally developed on 32bitsingleprecision computers
(VAXen and IBMs) which were too inaccurate for those computations.
Running DISORT.f on a typical 32bitsingleprecision computer
usually gives results precise to at least 23 significant digits,
although for certain special situations the precision can fall to
one significant digit. Results can even have NO significant digits
when they fall below about 10**(8) times the driving radiation,
as one can see in the fluxes in some of the test problems.
IEEE Underflow: On some Unix workstations (notably an HP9000),
turning on certain compiler options can also cause IEEE underflows
to trigger an error condition and abort the test run. On the
HP9000, this option was +T, which ostensibly produces a traceback
but as a side effect bombs the run. DISORT will never pass underflow
tests because it has potential underflows all over the place, esp.
in the computations of exponentials.
LINPACK

There are two versions of LINPACK, one in single precision
(routines begin with 'S') and one in double precision
(routines begin with 'D'). DISORT calls the 'S' versions.
These can still be used with an autodoubling compiler. But
if a local LINPACK rather than our file LINPAK.f is used in
running DISORT in double precision, change the first
letters of SGBCO, SGBFA, SGBSL, SGEFA, SGECO, SGESL to 'D'
(e.g. using a 'sed' command like that above).
Memory Usage

DISORT defines a considerable number of local arrays, whose default
sizes are rather large to accommodate the test problems. The size
of these arrays can be reduced by the user simply by altering the
PARAMETER statements (for MX...) just following the internal
variable documentation in subroutine DISORT. This can give a
dramatic reduction in memory usage for typical applications.
Elimination of this annoyance through use of dynamically dimensioned
arrays is one of the first goals of the f90 conversion of DISORT.

AUTHORS :
Knut Stamnes and Collaborators (kstamnes@stevenstech.edu)
Stevens Institute of Technology
Hoboken, New Jersey
SiChee Tsay (Tsay@gsfc.nasa.gov)
NASA Goddard Space Flight Center
Code 913
Greenbelt, MD 20771
Warren Wiscombe (Warren.Wiscombe@gsfc.nasa.gov)
NASA Goddard Space Flight Center
Code 913
Greenbelt, MD 20771
Stuart Freidenreich (vr@gfdl.gov)
NOAA Geophysical Fluid Dynamics Laboratory
Princeton, NJ 08540
Istvan Laszlo (laszlo@atmos.umd.edu)
University of Maryland
Dept. of Meteorology
College Park, MD 20742
Robert Pincus (Robert.Pincus@gsfc.nasa.gov)
Ran Song
NASA Goddard Space Flight Center
Code 913
Greenbelt, MD 20771
Collaborators: cf. REFERENCES

REFERENCES (cited in the programs using the acronyms shown) :
DGIS: Devaux, C., Grandjean, P., Ishiguro, Y. and C.E. Siewert,
1979: On MultiRegion Problems in Radiative Transfer,
Astrophys. Space Sci. 62, 225233
GS: Garcia, R.D.M. and C.E. Siewert, 1985: Benchmark
Results in Radiative Transfer, Transport Theory
and Statistical Physics 14, 437483
KS: Kylling, A. and K. Stamnes, 1992: Efficient yet accurate
solution of the linear transport equation in the
presence of internal sources: The exponentiallinear
indepth approximation, J. Comp. Phys., 102, 265276.
L: Lenoble, J., ed., 1985: Radiative Transfer in Absorbing
and Scattering Atmospheres: Standard Computational
Procedures, Deepak Publishing, Hampton, Virginia
NT: Nakajima, T. and M. Tanaka, 1988: Algorithms for Radiative
Intensity Calculations in Moderately Thick Atmospheres
Using a Truncation Approximation, J.Q.S.R.T. 40, 5169.
OS: Ozisik, M. and S. Shouman, 1980: Source Function
Expansion Method for Radiative Transfer in a TwoLayer
Slab, J.Q.S.R.T. 24, 441449
SS: Stamnes, K. and R. Swanson, 1981: A New Look at
the Discrete Ordinate Method for Radiative Transfer
Calculations in Anisotropically Scattering
Atmospheres, J. Atmos. Sci. 38, 387399
SD: Stamnes, K. and H. Dale, 1981: A New Look at the
Discrete Ordinate Method for Radiative Transfer
Calculations in Anisotropically Scattering
Atmospheres. II: Intensity Computations,
J. Atmos. Sci. 38, 26962706
S1: Stamnes, K., 1982: On the Computation of Angular
Distributions of Radiation in Planetary
Atmospheres, J.Q.S.R.T. 28, 4751
S2: Stamnes, K., 1982: Reflection and Transmission by
a Vertically Inhomogeneous Planetary Atmosphere,
Planet. Space Sci. 30, 727732
SC: Stamnes, K. and P. Conklin, 1984: A New MultiLayer
Discrete Ordinate Approach to Radiative Transfer
in Vertically Inhomogeneous Atmospheres,
J.Q.S.R.T. 31, 273282
SW: Sweigart, A., 1970: Radiative Transfer in Atmospheres
Scattering According to the Rayleigh Phase Function
with Absorption, The Astrophysical Journal
Supplement Series 22, 180
STWJ: Stamnes, K., S.C. Tsay, W. Wiscombe and K. Jayaweera, 1988:
A Numerically Stable Algorithm for
DiscreteOrdinateMethod Radiative Transfer in
Multiple Scattering and Emitting Layered Media,
Appl. Opt. 27, 25022509.
STWL: Stamnes, K., S.C. Tsay, W. Wiscombe and I. Laszlo:
A GeneralPurpose Numerically Stable Computer
Code for DiscreteOrdinateMethod Radiative
Transfer in Scattering and Emitting Layered Media,
DISORT Report v1.1 (2000)
VH1,VH2: Van de Hulst, H.C., 1980: Multiple Light Scattering,
Tables, Formulas and Applications, Volumes 1 and 2,
Academic Press, New York.
W: Wiscombe, W., 1977: The DeltaM Method: Rapid Yet
Accurate Radiative Flux Calculations, J. Atmos. Sci.
34, 14081422
++
PREFACE
DISORT was designed to be the most general and versatile
planeparallel radiative transfer program available, applicable
to problems from the ultraviolet to the radar regions of the
electromagnetic spectrum. As such, it has a rather large list
of input variables. This list is more easily comprehended if
several simple facts are borne in mind :
* there is one vertical coordinate, measured in optical depth
units, and two angular coordinates, one polar and one azimuthal;
* the layers and polar angles necessary for computational
purposes are *entirely decoupled* from the levels
and polar angles at which the user desires results.
The computational layering is usually constrained by the problem,
in the sense that each computational layer must be reasonably
homogeneous and not have a temperature variation of more than
about 510 K across it (if thermal sources are important).
For example, a dusty boundary layer topped by a cloud topped by
clear sky would suggest three computational layers, in the
absence of thermal sources.
Computational polar angles ('streams') are constrained by the
need for accuracy; for example, 4 streams may be enough for
accurate fluxes, while 16 streams or more may be necessary for
accurate intensities.
But the radiant quantities can be returned to the user at ANY
level and ANY angle. For example, the user may have picked 3
computational layers and 16 streams, but he can then request
intensities from only the middle of the 2nd layer, and only
in the nadir direction.
The deltaM transformation of Wiscombe (1977) is used in DISORT to
achieve optimum computational efficiency and accuracy for strongly
forwardpeaked phase functions. The "deltaMscaled" intensities are
"unscaled" by applying the TMS and IMS intensity corrections of
Nakajima and Tanaka (1988).
++
Developing a package such as this is a humbling experience.
Every time we thought it was finally free of errors, further
testing would reveal another. What seemed only a 6month project
thus stretched into 3 years; however, we think the result is worth
it. We believe this package to be freer of errors than any other
similar package available today, and more fullfeatured to boot.
Of course, we would be foolhardy to claim that a package as
large and complex as this one is entirely errorfree. We have
followed two cardinal principles of software development in an
effort to minimize errors:
(1) offloading hard but standard computational tasks onto excellent
software written by experts (in our case, the linear equation and
eigenvalue/vector computations)
(2) "unit testing" many subroutines outside the program, using
specially developed test drivers (for example, the Gauss and
Planck routines)
We are confident that the remaining errors are subtle and unlikely
to be encountered by the average user. If you do find any errors,
please report them to the authors, and we will do our best,
time permitting, to find a solution.
B E W A R E :
It is very easy to introduce errors into this package. We did
it many times ourselves in the course of developing it. The most
seemingly innocent, casual changes are fraught with danger. After a
severalyear debugging process, we are not prepared to find bugs
that YOU introduce. If you change the code, you are on your own.
++
INDEX CONVENTIONS ( for all variables described below ) :
IU : for user polar angles (where intensities are computed)
IQ : for computational polar angles ('quadrature angles')
J : for user azimuthal angles
K : for Legendre expansion coefficients
LU : for user levels (where fluxes and intensities
are computed)
LC : for computational layers (each having a different
singlescatter albedo and/or phase function)
LEV : for computational levels
LAYERING CONVENTION:
Layers are numbered from the top boundary down.
ANGLE CONVENTION:
Polar (zenith) angles are measured from the upward direction:
straight up is zero degrees and straight down is 180 degrees.
There is a small inconsistency in that, for historical reasons,
the cosine of the incident beam angle (UMU0) is positive,
whereas according to this convention it should be negative.
Azimuth angles are measured in an absolute frame of reference,
rather than from the plane of the incident beam; hence the
azimuth angle of the incident beam is an input variable.
There is nothing in this version of DISORT which can
introduce any further azimuth dependence, however, although
in Nature such things as plowed fields and oriented ice
crystals can introduce further absolute origins of azimuth
angle.
UNITS CONVENTION:
The radiant output units are determined by the sources of
radiation driving the problem. Lacking thermal emission, the
radiant output units are the same as the units of the beam
source FBEAM and the isotropic source FISOT. The whole problem
could then be nondimensionalized by setting these to unity.
If thermal emission of any kind is included, subprogram PLKAVG
determines the units. The default PLKAVG has MKS units (W/sq m).
Several users have rewritten PLKAVG to return just the temperature
(i.e. their only executable statement in PLKAVG is PLKAVG = temp.),
an approximation which is widely used in the longwavelength limit;
in this case, all radiant quantities are in degrees Kelvin. If you
rewrite PLKAVG, however, you must also put in new selftest
'correct answers' in subroutine SLFTST (or bypass it).
NOTE: make sure FBEAM and FISOT have the same units as PLKAVG
when thermal emission is present.
CAVEATS:
(1) The case singlescatteralbedo=1 causes removable 0/0type
singularities in our generalcase formulae. These can be
eliminated by applying L'Hospital's Rule. However, this
creates large amounts of specialcase code whose IFstatements
in some cases ruined the possibility of loop vectorization.
It also led to quantum jumps in results as one crossed from
s.s.albedo near unity to s.s.albedo exactly unity.
Thus, we chose instead to use a "dithering method" in which
a very small quantity (about 10x machine precision) is
subtracted from the singlescatter albedo. This works
surprisingly well, but also causes some loss of accuracy,
which can be seen in the test problems for which
singlescatteralbedo=1. A better method would be to
work out some kind of Laurent expansion near s.s.albedo=1
that merged smoothly with the generalcase formulae.
(2) Another removable 0/0type singularity condition arises
when the Sun (beam) angle coincides with one of the angles
at which output intensities are desired. In this
case, the user would be advised to slightly change their
sun angle or their output angle so that they no longer
coincide. The program handles this case using L'Hospital's
Rule to get the correct limit, but it is not sophisticated
and may amplify the error by not expanding to a higher
level of approximation.
This singularity also occurs when the beam angle coincides
with one of the quadrature angles, but the latter are not
under user control, and they take such unusual values that
the odds of such a coincidence are practically zero. The
problem is most likely to occur when NSTR/2 is odd and
UMU0 = 0.5, and it can easily be corrected by changing NSTR.
In general, it may be better to avoid requiring intensities
exactly at the beam angle. In the direct backscattering
region, real phase functions are most poorly known, especially
in the low order of Legendre approximation in which they are
represented in DISORT, and if one is looking for the
heiligenschein or opposition effect such as seen
when flying over vegetation, forget it  DISORT
doesn't calculate that. The region of direct forward
scattering is also difficult for DISORT, because in order to
do as well as it does at other angles it has to fiddle with
the photons scattered within a few degrees of the forward
direction; thus its exact forward intensity may be less
accurate than at other angles.
(3) If you flip between ONLYFL = TRUE and ONLYFL = FALSE in the
same run, your input UMU values in USRANG = TRUE cases will
be destroyed. Since such flipflopping is an extremely
unlikely usage scenario, no guards against this disaster
have been implemented. If you reset UMU before each call
to DISORT, this problem cannot occur.
++
I N P U T V A R I A B L E S
++
******** COMPUTATIONAL LAYER STRUCTURE ********
NLYR Number of computational layers
DTAUC(LC) LC = 1 to NLYR,
Optical depths of computational layers
SSALB(LC) LC = 1 to NLYR,
Singlescatter albedos of computational layers
NMOM Number of phase function moments (all the significant
Legendre coefficients) not including the zeroth moment.
Should be greater than or equal to NSTR in problems
with scattering. In problems without scattering,
NMOM is not used.
In a multilayer problem, NMOM number of moments
should be supplied for every layer, even when the
individual layers are characterized by different phase
functions.
PMOM(K,LC) K = 0 to NMOM, LC = 1 to NLYR,
Coefficients in Legendre polynomial expansions of
phase functions for computational layers :
P(mu) = sum,K=0 to NMOM( (2K+1) PMOM(K) PK(mu) )
WHERE P = phase function
mu = cos(scattering angle)
PK = Kth Legendre polynomial
The K = 0 coefficient should be unity (it will be
reset to unity in any case). Subroutine GETMOM,
supplied in the test problem file, may be used to
set coefficients in special cases.
TEMPER(LEV) LEV = 0 to NLYR, Temperatures (K) of levels.
(Note that temperature is specified at LEVELS
rather than for layers.) Be sure to put top level
temperature in TEMPER(0), not TEMPER(1). Top and
bottom level values do not need to agree with top and
bottom boundary temperatures (i.e. temperature
discontinuities are allowed).
Needed only if PLANK is TRUE.
******** USER LEVEL STRUCTURE ********
USRTAU = FALSE, Radiant quantities are to be returned
at boundary of every computational layer.
= TRUE, Radiant quantities are to be returned
at userspecified optical depths, as follows:
NTAU Number of optical depths
UTAU(LU) LU = 1 to NTAU, user optical depths,
in increasing order. UTAU(NTAU) must
not exceed the total optical depth
of the medium, as deduced from DTAUC.
******** COMPUTATIONAL POLAR ANGLE STRUCTURE ********
NSTR Number of computational polar angles to be used
(= number of 'streams') ( should be even and .GE. 2 ).
Note that these are Gaussian angles and hence the user
has no control over what values are used.
In general, the more streams used, the more accurate
the calculated fluxes and intensities will be. However,
there is no rigorous proof that increasing NSTR
produces a monotonic decrease in error; hence it is
possible that small increases in NSTR may make
the error slightly worse. Large increases in NSTR
(like doubling it), on the other hand, are almost
certain to reduce the error.
For NSTR = 2 a twostream program should be used
instead, since DISORT is not optimized for this case
except in the eigenvalue/vector routine. Also,
intensities will be totally unreliable for NSTR = 2,
since they are just extrapolations from a single point.
We only allow this case for our own consistency tests.
******** USER POLAR ANGLE STRUCTURE ********
USRANG = FALSE, Radiant quantities are to be returned
at computational polar angles. Also, UMU will
return the cosines of the computational polar
angles and NUMU will return their number
( = NSTR). UMU must be large enough to
contain NSTR elements (cf. MAXUMU).
= TRUE, Radiant quantities are to be returned
at userspecified polar angles, as follows:
NUMU No. of polar angles ( zero is a legal
value only when ONLYFL = TRUE )
UMU(IU) IU=1 to NUMU, cosines of output polar
angles in increasing order  starting
with negative (downward) values (if any)
and on through positive (upward) values;
*** MUST NOT HAVE ANY ZERO VALUES ***
** NOTE ** If only fluxes are desired (ONLYFL = TRUE), then
UMU will return the computational polar angles if it is
big enough to contain them (and NUMU will return the number of
such angles). This is so the user will know the angles that the
returned azimuthallyaveraged intensities refer to. But a bad
byproduct is that if the user flips between ONLYFL = TRUE and
ONLYFL = FALSE in the same run, his input UMU in USRANG = TRUE
cases will be destroyed. Thus, he should reset his input UMU
values prior to every DISORT call. (For USRANG = FALSE cases
there is no difficulty because UMU always returns computational
angles.)
********* AZIMUTHAL ANGLE STRUCTURE ***********
NPHI : Number of azimuthal angles at which to return
intensities ( zero is a legal value only when
ONLYFL = TRUE )
PHI(J) : J = 1 to NPHI, Azimuthal output angles (in degrees)
( not used when ONLYFL = TRUE )
********* TOP AND BOTTOM BOUNDARY CONDITIONS **********
IBCND = 0 : General case: boundary conditions any combination of:
* beam illumination from the top ( see FBEAM )
* isotropic illumination from the top ( see FISOT )
* thermal emission from the top ( see TEMIS, TTEMP )
* internal thermal emission sources ( see TEMPER )
* reflection at the bottom ( see LAMBER, ALBEDO, BDREF )
* thermal emission from the bottom ( see BTEMP )
= 1 : Return only albedo and transmissivity of the entire
medium vs. incident beam angle; see S2 for details.
(There can be no Planck sources in this case.)
Technically, this is accomplished by assuming an
isotropicallyincident source of radiation at the
top boundary, but this is of no real concern to the
user.
Many users overlook this option even though it
turns out to be exactly what they need.
The only input variables considered in this case
are NLYR, DTAUC, SSALB, PMOM, NSTR, USRANG, NUMU, UMU,
ALBEDO, PRNT, HEADER and the array dimensions (see below).
PLANK is assumed FALSE, LAMBER is assumed TRUE, and the
bottom boundary can have any ALBEDO. The sole output is
ALBMED, TRNMED; since these are just ratios, this option
does not use source strength information in FBEAM or
FISOT.
UMU is interpreted as the array of beam
angles in this case. If USRANG = TRUE they must be
positive and in increasing order, and will be returned
this way; internally, however, the negatives of the
UMU's are added, so MAXUMU must be at least 2*NUMU.
If USRANG = FALSE, UMU is returned as the NSTR/2
positive quadrature angle cosines, in increasing
order.
NOTE: The intensities involved here are not corrected for
deltaM scaling effects.
FBEAM : Intensity of incident parallel beam at top boundary.
[same units as PLKAVG (default W/sq m) if thermal
sources active, otherwise arbitrary units].
Corresponding incident flux is UMU0 times FBEAM.
Note that this is an infinitely wide beam, not a
searchlight beam.
UMU0 : Polar angle cosine of incident beam (positive).
If this equals the negative of one of the UMU values,
specialcase formulae must be invoked to prevent
a 0/0type removable singularity.
** WARNING ** If this equals one of the
computational polar angle cosines, a singularity
occurs; hence this is treated as a fatal
error. The problem is most likely to
occur when NSTR/2 is odd and UMU0 = 0.5;
otherwise, it is almost impossible to hit a
computational angle by chance. The problem can
easily be corrected by changing NSTR.
PHI0 : Azimuth angle of incident beam (0 to 360 degrees)
FISOT : Intensity of topboundary isotropic illumination.
[same units as PLKAVG (default W/sq m) if thermal
sources active, otherwise arbitrary units].
Corresponding incident flux is pi (3.14159...)
times FISOT.
LAMBER : TRUE, isotropically reflecting bottom boundary.
In this case must also specify :
ALBEDO : bottomboundary albedo
FALSE, bidirectionally reflecting bottom boundary.
In this case, the bottom bidirectional reflectivity
regarded as a function of the cosine of angle of
reflection MU, the cosine of angle of incidence
MUP, and the difference of azimuth angles of
incidence and reflection DPHI (in radians), must
be supplied in the function subprogram:
REAL FUNCTION BDREF( WVNMLO, WVNMHI, MU, MUP, DPHI )
In DISORT, the bidirectional reflectivity is defined
as (see STWL Eq. 39):
UU( MU, PHI ) = 1 / PI *
Integral over PHI from 0 to 2PI [
Integral over MUP from 0 to 1 [
BDREF( WVNMLO, WVNMHI, MU, MUP, DPHI ) *
UU( MUP, PHIP ) * MUP d MUP ] d PHIP ],
where all the variables are as defined above, UU and
PHIP are the incident intensity and the azimuth angle
of incidence, respectively.
Note that in BDREF the cosine of the incident angle
(MUP) is positive. (Another inconsistency: according
to the convention used in DISORT it should be
negative.)
The values of MU and MUP are determined by UMU and
the Gaussian angles corresponding to NSTR and to
those used for calculating the flux albedo. In BDREF,
the bidirectional reflectivities should be specified
for these angles. Because the user has no control
over the Gaussian angles, the simplest way is to
supply an analytical function of the bidirectional
reflectance. In case the reflectivities are only
available at discrete angles the user should include
an interpolation routine in BDREF which would
calculate the reflectivities at arbitrary angles.
In its default form, DISORT only includes a "stub"
version of BDREF in file BDREF.f which only prints a
message telling the user to include his/her own
function BDREF. The function BDREF supplied in file
DISOTEST.f, that is set to give bidirectional
reflectivity in a special case, can be used as an
example.
** NOTE 1 ** Flux albedos calculated from function
BDREF will be checked to be sure they lie between
zero and one for all possible incidence angles.
** NOTE 2 ** The lower and upper boundaries of the
spectral interval WVNMLO and WVNMHI may be used to
specify a spectral dependent lower boundary in
multiple DISORT runs.
In principle, the function BDREF could also be used to
specify an isotropically reflecting lower boundary.
However, this is NOT recommended since calculation of
the Fourier expansion coefficients increases execution
time; which can be avoided using the LAMBER=TRUE option.
** NOTE ** For DISORT to successfully compile, the
function BDREF must always be present, even though it is
not used when LAMBER is TRUE.
BTEMP : Temperature of bottom boundary (K) (bottom emissivity
is calculated from ALBEDO or function BDREF, so it need
not be specified).
Needed only if PLANK is TRUE.
TTEMP : Temperature of top boundary (K).
Needed only if PLANK is TRUE.
TEMIS : Emissivity of top boundary.
Needed only if PLANK is TRUE.
********** CONTROL FLAGS **************
PLANK TRUE, include thermal emission
FALSE, ignore all thermal emission (saves computer time)
( if PLANK = FALSE, it is not necessary to set any of
the variables having to do with thermal emission )
WVNMLO, Wavenumbers (inv cm) of spectral interval of interest
WVNMHI ( used only for calculating Planck function ).
Needed only if PLANK is TRUE, or in multiple runs, if
LAMBER is FALSE and BDREF depends on spectral interval.
ONLYFL TRUE, return fluxes, flux divergences, and mean
intensities.
FALSE, return fluxes, flux divergences, mean
intensities, AND intensities.
ACCUR Convergence criterion for azimuthal (Fourier cosine)
series. Will stop when the following occurs twice:
largest term being added is less than ACCUR times
total series sum. (Twice because there are cases where
terms are anomalously small but azimuthal series has
not converged.) Should be between 0 and 0.01 to avoid
risk of serious nonconvergence. Has no effect on
problems lacking a beam source, since azimuthal series
has only one term in that case.
PRNT(L) Array of LOGICAL print flags causing the following prints:
L quantities printed
 
1 input variables (except PMOM)
2 fluxes
3 intensities at user levels and angles
4 planar transmissivity and planar albedo
as a function solar zenith angle ( IBCND = 1 )
5 phase function moments PMOM for each layer
( only if PRNT(1) = TRUE, and only for layers
with scattering )
HEADER A 127 (or less) character header for prints, embedded in
a DISORT banner; setting HEADER = '' (the null string)
will eliminate both the banner and the header, and this
is the only way to do so (HEADER is not controlled by any
of the PRNT flags); HEADER can be used to mark the
progress of a calculation in which DISORT is called
many times, while leaving all other printing turned off.
****** ARRAY DIMENSIONS *******
MAXCLY : Dimension of DTAUC, SSALB, TEMPER. 2nd dimension of
PMOM. Should be .GE. NLYR.
Max. number of computational layers
MAXMOM : First dimension of PMOM. Should be .GE. NMOM.
Max. number of Legendre coefficients.
MAXULV : Dimension of UTAU, RFLDIR, RFLDN, FLUP, DFDT.
2nd dimension of U0U, UU.
Should be .GE. NTAU if USRTAU=TRUE,
.GE. NLYR+1 if USRTAU=FALSE.
Max. number of user levels
MAXUMU : Dimension of UMU. First dimension of UU, U0U.
Should be .GE. NUMU if USRANG=TRUE,
.GE. NSTR if USRANG=FALSE.
Max. number of user polar angles
MAXPHI : Dimension of PHI. 3rd dimension of UU.
Should be .GE. NPHI.
Max. number of user azimuth angles.
++
O U T P U T V A R I A B L E S
++
== NOTE ON UNITS ==
If thermal sources are specified, fluxes come out in the units
of PLKAVG, intensities in those units per steradian. Otherwise,
the flux and intensity units are determined by the units of
FBEAM and FISOT.
== NOTE ON ZEROING ==
All output arrays are completely zeroed on each call to DISORT
before being partially refilled with results. This keeps garbage
from accumulating in unused parts of the output arrays, and
keeps indefinites and NotaNumbers out of unset parts of the
output arrays. With the modern emphasis on arrayprocessing,
it is wise to keep entire arrays 'clean' even if only parts
contain useful results. Symbolic dumps also look cleaner.
If USRTAU = FALSE :
NTAU Number of optical depths at which radiant
quantities are evaluated ( = NLYR+1 )
UTAU(LU) LU = 1 to NTAU; optical depths, in increasing
order, corresponding to boundaries of
computational layers (see DTAUC)
If USRANG = FALSE or (ONLYFL=TRUE and MAXUMU.GE.NSTR):
NUMU No. of computational polar angles at which radiant
quantities are evaluated ( = NSTR )
UMU(IU) IU = 1 to NUMU; cosines of computational polar
angles, in increasing order. All positive if
IBCND = 1, otherwise half negative (downward)
and half positive (upward).
RFLDIR(LU) Directbeam flux (without deltaM scaling)
RFLDN(LU) Diffuse downflux (total minus directbeam)
(without deltaM scaling)
FLUP(LU) Diffuse upflux
DFDT(LU) Flux divergence d(net flux)/d(optical depth),
where 'net flux' includes the direct beam
(an exact result; not from differencing fluxes)
UAVG(LU) Mean intensity (including the direct beam)
(Not corrected for deltaMscaling effects)
UU(IU,LU,J) Intensity ( if ONLYFL = FALSE; zero otherwise )
ALBMED(IU) Albedo of the medium as a function of incident
beam angle cosine UMU(IU) (IBCND = 1 case only)
TRNMED(IU) Transmissivity of the medium as a function of incident
beam angle cosine UMU(IU) (IBCND = 1 case only)
++
>>>>>>>>>> ROUTINES AND PURPOSES <<<<<<<<<<<<<
file DISORT.f :

ALBTRN: manages the IBCND = 1 SPECIAL CASE
ALTRIN: calculates azimuthallyaveraged intensity (equal to
albedo and/or transmissivity) at user angles
for the IBCND = 1 special case
ASYMTX: solves eigenfunction problem for real asymmetric matrix
whose eigenvalues are known to be real
CHEKIN: checks input variables for errors
CMPINT: computes intensities at user levels and computational
angles
DREF : flux albedo as a function of incident angle when a
bottomboundary bidirectional reflectivity is
specified (LAMBER = FALSE)
FLUXES: computes upward and downward fluxes
INTCOR: corrects intensity field by using NakajimaTanaka
algorithm
LEPOLY: evaluates normalized associated Legendre polynomials
PLKAVG: computes integral of Planck function over a
wavenumber interval
PRALTR: prints flux albedo and transmissivity of medium
(IBCND = 1 special case)
PRAVIN: prints azimuthally averaged intensities
PRTINP: prints input variables
PRTINT: prints intensities at user angles
QGAUSN: computes Gaussian quadrature points and weights
RATIO : computes ratio of two numbers in failsafe way
SECSCA: calculates secondary scattered intensity for intensity
correction
SETDIS: performs miscellaneous settingup operations
SETMTX: calculates coefficient matrix for linear equations
embodying boundary and layer interface conditions
SINSCA: calculates singlescattered intensity for intensity
correction
SLFTST: sets input for selftest and checks for failure
SOLEIG: solves eigenfunction problem for a single layer
SOLVE0: solves linear equations embodying boundary and layer
interface conditions (general boundary conditions)
SOLVE1: solves linear equations embodying boundary and layer
interface conditions (IBCND = 1 case)
SPALTR: calculates spherical albedo, transmissivity
(IBCND = 1 special case)
SURFAC: calculates surface bidirectional reflectivity/emissivity
TERPEV: interpolates eigenvectors to user angles
TERPSO: interpolates particular solutions to user angles
UPBEAM: finds particular solution for beam source
UPISOT: finds particular solution for thermal emission source
USRINT: computes intensities at user levels and user angles
XIFUNC: calculates Xi function used in SECSCA
ZEROAL: zeros a group of matrices
ZEROIT: zeros a given matrix
file ErrPack.f :

ERRMSG: prints out error message and kills run if fatal.
also prevents errormessage runaway.
TSTBAD: prints message when selftest fails
WRTBAD: writes out names of erroneous input variables
WRTDIM: writes out names of toosmall symbolic dimensions
file LINPAK.f :

SGBCO,: LU decompose a banded matrix
SGBFA
SGBSL : solve linear equations with banded matrix of
coefficients after LU decomposition of matrix
SGECO,: LU decompose a general matrix
SGEFA
SGESL : solve linear equations with general matrix of
coefficients after LU decomposition of matrix
SASUM : sum of elements of vector
SAXPY : scalar times one vector plus another vector
SDOT : dot product of two vectors
SSCAL : scalar times a vector
ISAMAX: index of maximum element of a vector
R1MACH.f : returns machinespecific constants (singleprecision)
D1MACH.f : returns machinespecific constants (doubleprecision)
In addition to the routines above, the user also must supply the
function:
BDREF : sets bidirectional bottom boundary; only used when
LAMBER = FALSE, but must be supplied at all times
for DISORT to compile
Full Call Tree (nx means called n times):
DISORT+(R1MACH)

+SLFTST (1)+(TSTBAD) (4x)
 
 +(ERRMSG)

+ZEROIT

+CHEKIN+(WRTBAD)
 
 +(ERRMSG)
 
 +(WRTBAD) (31x)
 
 +DREF+(ERRMSG) (2x)
 
 +(WRTBAD) (7x)
 
 +(WRTDIM) (9x)
 
 +(ERRMSG) (3x)

+ZEROAL

+SETDIS+QGAUSN (2)+(D1MACH)
  
  +(ERRMSG) (2x)
 
 +(ERRMSG)

+PRTINP

+ALBTRN+LEPOLY (3)(ERRMSG)
 
 +LEPOLY see 3
 
 +ZEROIT
 
 +SOLEIG (4)+ASYMTX+(D1MACH)
   
   +(ERRMSG) (2x)
  
  +(ERRMSG)
 
 +TERPEV
 
 +SETMTX (5)ZEROIT
 
 +SOLVE1+ZEROIT
  
  +(SGBCO)
  
  +(ERRMSG)
  
  +(SGBSL)
 
 +ALTRIN
 
 +SPALTR (2x)
 
 +PRALTR

+ZEROIT

+PLKAVG (6)+(R1MACH) (2x)
 
 +(ERRMSG) (3x)

+PLKAVG see 6 (2x)

+LEPOLY see 3 (3x)

+SURFAC+QGAUSN see 2
 
 +ZEROIT (2x)
 
 +BDREF (3x)
 
 +ZEROIT (2x)
 
 +BDREF (3x)

+SOLEIG see 4

+UPBEAM+(SGECO)
 
 +(ERRMSG)
 
 +(SGESL)

+UPISOT+(SGECO)
 
 +(ERRMSG)
 
 +(SGESL)

+TERPEV

+TERPSO

+SETMTX see 5

+SOLVE0+ZEROIT
 
 +(SGBCO)
 
 +(ERRMSG)
 
 +(SGBSL)
 
 +ZEROIT

+FLUXES+ZEROIT (3x)

+ZEROIT

+USRINT

+CMPINT

+PRAVIN

+ZEROIT

+RATIO+(R1MACH) (2x)

+INTCOR+SINSCA (2x)
 
 +SECSCA+XIFUNC

+PRTINT

+SLFTST see 1