-
Notifications
You must be signed in to change notification settings - Fork 1
/
disort.txt
executable file
·1267 lines (965 loc) · 49 KB
/
disort.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
***********************************************************************
Documentation File for 'DISORT', A Discrete Ordinates Radiative
Transfer Fortran-90 Program for a Multi-Layered Plane-Parallel 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.
FORTRAN-90: This version and v1.3 of DISORT require a Fortran-90
compiler. So far, the only Fortran-90 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 Fortran-90, 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 e-mail 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 Si-chee
Tsay in his thesis under Knut Stamnes, largely during an intense four-week
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 NON-UNIT
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 Fortran-90 version;
R1MACH, D1MACH use Fortran-90 intrinsic functions;
(Many thanks to Robert Pincus for the conversion to LAPACK)
version 2.0, 2000 (Mar): Introduced Nakajima-Tanaka 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) azimuthally-averaged intensity is no longer returned as an output;
(d) the mean intensities are not intensity-corrected; they are already
very accurate for small N using the usual delta-M approximation
and benefit little from the corrections; however, it is a technical
inconsistency since they are derived from the uncorrected
azimuthally-averaged 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 Legendre-polynomial 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 quadrature-angles 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 auto-doubling option now
available on many compilers (if you auto-double, however, you need to
change all instances of 'R1MACH' to 'D1MACH'). This is also the
version to use on Cray and other high-accuracy machines.
BDREF.f
A default FUNCTION subprogram to supply surface bidirectional
reflectivity. The user must replace this with their own version
when if a non-Lambertian 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 single-double precision
on an IEEE-arithmetic Unix workstation. Filename may have a
version number included.
RDI1MACH.f
Fortran-90 Routines for returning machine constants, from netlib author
Eric Grosse; also see file RDI1MACH.readme.
LINPAK.f : Routines used by DISORT from the linear-equation-solving
public-domain packages LINPACK and BLAS (available
at most computer sites). Included only to make package
self-contained. 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
super-vector-speed LINPACK in 'libsci'.)
The entire LINPACK can be obtained from NetLib,
http://www.netlib.org/.
ErrPack.f
Error-handling routines. Standard across almost all programs on
the wiscombe anonymous ftp site.
rundis
Example C-shell 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 stub-file 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 user-controllable
(through the number of streams NSTR) "truncation error". DISORT
is free of the *catastrophic* growth of roundoff error that
plagued pre-1980 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 32-bit-single-precision computers
(VAXen and IBMs) which were too inaccurate for those computations.
Running DISORT.f on a typical 32-bit-single-precision computer
usually gives results precise to at least 2-3 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 ([email protected])
Stevens Institute of Technology
Hoboken, New Jersey
Si-Chee Tsay ([email protected])
NASA Goddard Space Flight Center
Code 913
Greenbelt, MD 20771
Warren Wiscombe ([email protected])
NASA Goddard Space Flight Center
Code 913
Greenbelt, MD 20771
Stuart Freidenreich ([email protected])
NOAA Geophysical Fluid Dynamics Laboratory
Princeton, NJ 08540
Istvan Laszlo ([email protected])
University of Maryland
Dept. of Meteorology
College Park, MD 20742
Robert Pincus ([email protected])
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 Multi-Region Problems in Radiative Transfer,
Astrophys. Space Sci. 62, 225-233
GS: Garcia, R.D.M. and C.E. Siewert, 1985: Benchmark
Results in Radiative Transfer, Transport Theory
and Statistical Physics 14, 437-483
KS: Kylling, A. and K. Stamnes, 1992: Efficient yet accurate
solution of the linear transport equation in the
presence of internal sources: The exponential-linear-
in-depth approximation, J. Comp. Phys., 102, 265-276.
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, 51-69.
OS: Ozisik, M. and S. Shouman, 1980: Source Function
Expansion Method for Radiative Transfer in a Two-Layer
Slab, J.Q.S.R.T. 24, 441-449
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, 387-399
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, 2696-2706
S1: Stamnes, K., 1982: On the Computation of Angular
Distributions of Radiation in Planetary
Atmospheres, J.Q.S.R.T. 28, 47-51
S2: Stamnes, K., 1982: Reflection and Transmission by
a Vertically Inhomogeneous Planetary Atmosphere,
Planet. Space Sci. 30, 727-732
SC: Stamnes, K. and P. Conklin, 1984: A New Multi-Layer
Discrete Ordinate Approach to Radiative Transfer
in Vertically Inhomogeneous Atmospheres,
J.Q.S.R.T. 31, 273-282
SW: Sweigart, A., 1970: Radiative Transfer in Atmospheres
Scattering According to the Rayleigh Phase Function
with Absorption, The Astrophysical Journal
Supplement Series 22, 1-80
STWJ: Stamnes, K., S.-C. Tsay, W. Wiscombe and K. Jayaweera, 1988:
A Numerically Stable Algorithm for
Discrete-Ordinate-Method Radiative Transfer in
Multiple Scattering and Emitting Layered Media,
Appl. Opt. 27, 2502-2509.
STWL: Stamnes, K., S.C. Tsay, W. Wiscombe and I. Laszlo:
A General-Purpose Numerically Stable Computer
Code for Discrete-Ordinate-Method 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 Delta-M Method: Rapid Yet
Accurate Radiative Flux Calculations, J. Atmos. Sci.
34, 1408-1422
+---------------------------------------------------------------------+
PREFACE
DISORT was designed to be the most general and versatile
plane-parallel 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 5-10 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 delta-M transformation of Wiscombe (1977) is used in DISORT to
achieve optimum computational efficiency and accuracy for strongly
forward-peaked phase functions. The "delta-M-scaled" 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 6-month 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 full-featured to boot.
Of course, we would be foolhardy to claim that a package as
large and complex as this one is entirely error-free. 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
several-year 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
single-scatter 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 non-dimensionalized 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 long-wavelength limit;
in this case, all radiant quantities are in degrees Kelvin. If you
rewrite PLKAVG, however, you must also put in new self-test
'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 single-scatter-albedo=1 causes removable 0/0-type
singularities in our general-case formulae. These can be
eliminated by applying L'Hospital's Rule. However, this
creates large amounts of special-case code whose IF-statements
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 single-scatter albedo. This works
surprisingly well, but also causes some loss of accuracy,
which can be seen in the test problems for which
single-scatter-albedo=1. A better method would be to
work out some kind of Laurent expansion near s.s.-albedo=1
that merged smoothly with the general-case formulae.
(2) Another removable 0/0-type 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 flip-flopping 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,
Single-scatter 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 multi-layer 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 = K-th 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 user-specified 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 two-stream 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 user-specified 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 azimuthally-averaged 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
isotropically-incident 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
delta-M 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,
special-case formulae must be invoked to prevent
a 0/0-type 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 top-boundary 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 : bottom-boundary 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 non-convergence. 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 Not-a-Numbers out of unset parts of the
output arrays. With the modern emphasis on array-processing,
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) Direct-beam flux (without delta-M scaling)
RFLDN(LU) Diffuse down-flux (total minus direct-beam)
(without delta-M scaling)
FLUP(LU) Diffuse up-flux
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 delta-M-scaling 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 azimuthally-averaged intensity (equal to
albedo and/or transmissivity) at user angles