This appendix contains a FORTRAN program which uses information on location, time, and weather conditions to estimate parameters used in race simulation. It estimates the position of the sun in the sky, scattering and absorption of solar radiation by a cloudless atmosphere, and transmittance of solar radiation by clouds. Transmittance and backscattering by the cloudless atmosphere is considered to occur before the influence of clouds. Although being somewhat unrealistic, this provides reasonably accurate estimates.
The solar position algorithm is from the Astronomical Almanac. It calculates solar position accurate to within .1° in the sky hemisphere, which is more accurate than is necessary for my purposes. The clear sky scattering and absorption algorithm is from Hoyt. It models the effects of "dry air", aerosols, water vapor, and ozone on "broadband" insolation, separating it into direct and diffuse components, and includes multiple scattering between the earth and sky. The algorithm estimates the CSSF accurate to within 3% in most cases, and within 5% in almost all cases. I have modified this somewhat to separate the aerosol-scattered diffuse solar radiation, which is anisotropic, into equivalent beam and isotropic components. This is useful for solar array modeling, which requires some specification of the directional dependence of diffuse solar radiation. Molecular, cloud-scattered, and sky backscattered solar radiation are all assumed to be isotropic. This is probably a worse approximation for cloud-scattering than for the other two. The effects of clouds are considered by estimating an effective cloud shadowing fraction and mean transmittance. When working on this modeling, I considered "cloud layer models", such as those by Davies and McKay. However, in practical applications, I have found it most reasonable to just estimate a single total cloud shadowing fraction and mean transmittance. The effect of multiple layers can be estimated as shown in figure A1. The cloud shadowing fraction should be thought of as the fraction of direct solar radiation that is extinct from the solar beam by clouds at a given time or over a given time interval. This preserves the total amount of direct solar radiation in a given situation, which is important when considering collection by surfaces of specific orientations, such as a solar array. For example, consider a Cs cloud cover of .8 with an average transmittance of .8 (figure A2). In this case, much of the solar beam will penetrate the clouds intact, resulting in a total direct solar radiation amount of considerably greater than 20% of that which would be seen under a cloudless sky. A cloud shadowing fraction of .4 with an average transmittance of .6 is more appropriate, which increases the direct solar radiation amount to 60% of its cloudless sky value and preserves the global solar radiation amount. In some cases, the cloud transmittance can be estimated as the average cloud-shadowed GSR divided by its cloudless sky value, and then decreased slightly, since there will typically be a higher diffuse solar flux from the surrounding clear sky areas than from cloudy areas. When possible, these can be reasonably estimated from measurements.
Below is the program with numbered statements, and brief identifications and explanations of some program variables. The results of a sample run are also shown. All solar flux variables are in W/m².
1 OPEN(1,FILE="WEA-B.DAT")
2 OPEN(2,FILE="FLUX.DAT")
3 PI=3.141593
4 DTR=PI/180.
5 READ(1,*)ID,ZH,ZL
6 1 READ(1,*,END=2)RLAT,RLON,B,W,O3,A0,TIME,CS,TR,P0,IWIND,T,DP,
&RHO,HEAD,ELEV
7 RLAT=RLAT*DTR
8 HEAD=HEAD*DTR
9 TK=T+3.*ELEV/500.+273.15
10 DP=DP+3.*ELEV/500.
11 GA=9.80616*(1.-.00259*COS(RLAT+RLAT))*(1.-.314E-6*ELEV)
12 H$=287.05*TK/GA
13 P$=P0*EXP(-ELEV/H$)
14 E=6.112*EXP(17.67*DP/(DP+243.5))
15 R=.62197*E/(P$-E)
16 TV=TK*(1.+.608*R)
17 H=287.05*TV/GA
18 P=P0*EXP(-ELEV/H)
19 IWDIR=IWIND/100
20 IWSPD=IWIND-IWDIR*100
21 WN=IWSPD*.5144*COS(IWDIR*10.*DTR)
22 WW=-IWSPD*.5144*SIN(IWDIR*10.*DTR)
23 BD=-2557.5+ID+(TIME+ZH)/24.
24 SLONG=(280.46+.9856474*BD)*DTR
25 L=-SLONG/2./PI+1
26 SLONG=SLONG+L*2.*PI
27 G=(357.528+.9856003*BD)*DTR
28 L=-G/2./PI+1
29 G=G+L*2.*PI
30 SE=1.00014-.01671*COS(G)-.0014*COS(2.*G)
31 ELO=SLONG+(1.915*SIN(G)+.02*SIN(2.*G))*DTR
32 OE=(23.439-.0000004*BD)*DTR
33 RNUM=COS(OE)*SIN(ELO)
34 DENOM=COS(ELO)
35 RA=ATAN2(RNUM,DENOM)
36 IF(RA.GT.0..AND.SIN(ELO).LE.0.)RA=RA+PI
37 IF(RA.LT.0..AND.SIN(ELO).GE.0.)RA=RA+PI
38 IF(RA.LT.0..AND.SIN(ELO).LE.0.)RA=RA+2.*PI
39 IF(RA.EQ.0..AND.SIN(ELO).LE.0.)RA=RA+PI
40 EQTM=(SLONG-RA)/DTR/15.
41 DECL=ASIN(SIN(OE)*SIN(ELO))
42 STIME=TIME-(RLON-ZL)/15.+EQTM-12.
43 HAS=-STIME*15.*DTR
44 Q=SIN(RLAT)*SIN(DECL)+COS(RLAT)*COS(DECL)*COS(HAS)
45 ZED=ACOS(Q)
46 IF(ZED.EQ.0.)THEN
47 ASUN=0.
48 ELSE
49 ASUN=ASIN(COS(DECL)*SIN(HAS)/SIN(ZED))
50 ENDIF
51 IF(RLAT.EQ.0.)GOTO3
52 IF(Q.GE.SIN(DECL)/SIN(RLAT))GOTO4
53 3 IF(ASUN.LT.0.)ASUN=ASUN+2.*PI
54 ASUN=PI-ASUN
55 4 COR=1./(.955+20.267*Q)-.047121
56 IF(COR.GE.2.)COR=0.
57 COSZ=Q+.0083*COR
58 IF(COSZ.LT.0.)THEN
59 COSZ=Q
60 ASUN=0.
61 ENDIF
62 ZSUN=ACOS(COSZ)
63 ZSUND=ZSUN/DTR
64 IF(COSZ.LT.0.)THEN
65 RBEAM=0.
66 RISO=0.
67 DIR=0.
68 DIF=0.
69 GLB=0.
70 GOTO5
71 ELSE
72 AM=1./(COSZ+.15*(93.885-ZSUND)**(-1.253))*(P/1013.25)
73 WAM=1./(COSZ+.0548*(92.65-ZSUND)**(-1.452))
74 OAM=1.00314/SQRT(COSZ*COSZ+.0063)
75 EW=W*WAM
76 EO3=O3*OAM
77 AW=.1*(.75*EW+.000631)**.3-.0121
78 AO3=.045*(EO3+.000834)**.38-.0031
79 AG=.00235*(126.*AM+.0129)**.26+.0075*AM**.875-.00075
80 TM=.616+.3756*EXP(-.2212*AM)
81 TAS=(1.909*(EXP(-.667*B)-1.)+1.)**AM
82 AA=A0*(1.-TAS)
83 TMA=1.-AW-AO3-AG-AA
84 TDIR=TMA*TM*TAS
85 AER=.71*TMA*(1.-TAS)
86 FB=1./SQRT(1.+AM)
87 FI=1.-FB
88 TDIFB=FB*AER
89 TDIFI=.46*TMA*(1.-TM)+FI*AER
90 AMA=1.67*(P/1013.25)
91 AGZ=.00235*(126.*AMA+.0129)**.26+.0075*AMA**.875-.00075
92 TMZ=.616+.3756*EXP(-.2212*AMA)
93 TASZ=(1.909*(EXP(-.667*B)-1.)+1.)**AMA
94 AAZ=A0*TASZ
95 RMR=RHO*(1.-AW-AO3-AGZ-AAZ)*(.54*(1.-TMZ)+.29*(1.-TASZ))
96 ETRAD=1372.*COSZ/SE/SE
97 DIRC=ETRAD*(TDIR+TDIFB)
98 DFFI=ETRAD*TDIFI
99 RN=DIRC/COSZ
100 DOWN=DIRC+DFFI
101 DFFB=DOWN*RMR
102 CL=1.-CS
103 RBEAM=RN*CL
104 DFFC=DOWN*CS*TR
105 RISO=(DFFI+DFFB)*CL+DFFC
106 DIR=ETRAD*TDIR*CL
107 GLB=RBEAM*COSZ+RISO
108 DIF=GLB-DIR
109 ENDIF
110 5 WRITE(2,6)TIME,AM,ZSUN,ASUN,RBEAM,RISO,GLB,WN,WW,T,HEAD,RHO
111 GOTO1
112 2 CLOSE(1)
113 CLOSE(2)
114 6 FORMAT(F6.3,F8.4,F7.4,F8.4,3F7.1,2F6.1,F6.1,F7.3,F5.2)
115 END
Lines 1-2 open the input file, WEA-B.DAT (WEA stands for weather, and B for broadband), and the
output file FLUX.DAT
Lines 3-4 are numbers used for pi and conversions from degrees to radians.
Lines 5-6 read input data:
ID Day of year, 1993A
ZH Difference in hours from UTC (5 for EST)
ZL Time meridian in °W (75 for EST)
RLAT Location latitude in °N (- for deg S)
RLON Location longitude in °W
B Aerosol scattering parameter; this is about .18 for typical conditions
W Atmospheric water vapor thickness in cm
O3 Atmospheric ozone thickness in cm
A0 Aerosol fractional absorption parameter; this is typically about .07
TIME Local time as defined by ZH and ZL
CS Cloud shadowing fraction
TR Average cloud transmittance
P0 Sea level atmospheric pressure in mb
IWIND Wind (on solar car) in airways code format (e.g., 2618 is 260° at 18 kt (9.3 m/s))
T Temperature (around solar car) in °C
DP Dew point (around solar car) in °C
RHO Surrounding albedo
HEAD Heading of solar car (+ clockwise from N)
ELEV Location elevation in m
Lines 7-8 convert latitude and heading from degrees to radians.
Lines 9-18 estimate surface pressure. The variables there refer to values used in the pressure
correction:
TK Temperature in °KB
DP Dew point in °KB
GA Acceleration due to gravity in m/s²
E Water vapor pressure in mb
R Water vapor mixing ratio
TV Virtual temperature in °K
H Scale height in m
P Surface atmospheric pressure in mb
Lines 19-22 convert wind in airways code format into north-to-south and west-to-east components
in m/s:
WN North-to-south wind component in m/s
WW West-to-east wind component in m/s
Lines 23-63 estimate solar position:
BD Number of days before J.2000 (1200 UTC 31 December 1999)
SE Sun to earth distance in astronomical units
OE Obliquity of the ecliptic in radiansC
EQTM Correction factor for standard time in hours ("equation of time")
DECL Solar declination in radians
STIME Solar time in hours; defined here using a standard second time interval with 0 being
solar noon
HAS Hour angle of sun in radians (actually, same as solar time defined here)
Q Cosine of actual solar zenith angle
ZED Actual solar zenith angle in radians
ASUN Solar azimuth angle (0 is S, + to the E)
COSZ Cosine of apparent solar zenith angle (due to atmospheric refraction); Refers to the
position at which we actually see the sun in the sky
ZSUN Apparent solar zenith angle in radians
ZSUND Apparent solar zenith angle in °
Lines 64-70 set solar flux amounts equal to 0 if the sun (actually the center of the
sunD) is below the horizon.
Lines 71-95 estimate the effects of cloudless sky scattering and absorption on solar radiation:
AM Optical air mass
WAM Water vapor equivalent optical air mass
OAM Ozone equivalent optical air mass
AW Fractional absorption of solar radiation by water vapor
AO3 Fractional absorption of solar radiation by ozone
AG Fractional absorption of solar radiation by major atmospheric gasses
TM Transmittance of solar radiation due to molecular scattering
TAS Transmittance of solar radiation due to aerosol scattering
AA Fractional absorption of solar radiation due to aerosols
TMA Transmittance of solar radiation by cloudless sky absorbers
TDIR Transmittance of direct solar radiation
AER Fractional amount of forward scattering by aerosols
FB Equivalent beam component factor for aerosol scattering
FI Equivalent isotropic component factor for aerosol scattering
TDIFB Transmittance of diffuse solar radiation - equivalent beam component
TDIFI Transmittance of diffuse solar radiation - equivalent isotropic component
RMR Backscattering fraction of horizontal downward incident solar radiation
Lines 97-108 estimate equivalent beam and isotropic solar radiation components, and calculate
direct, diffuse, and global horizontal solar radiation.
ETRAD Extraterrestrial horizontal solar radiation; A "solar constant" of 1372 W/m² is used
DIRC Cloudless sky horizontal beam component of solar radiation
DFFI Cloudless sky horizontal downward isotropic component of solar radiation
RN Cloudless sky normal beam solar radiation
DOWN Total horizontal downward solar radiation
DFFB Backscattered horizontal solar radiation
CL Equivalent cloudless sky fraction
RBEAM Normal beam component of solar radiation
DFFC Horizontal solar radiation transmitted by clouds
RISO Horizontal isotropic component of solar radiation
DIR Direct solar radiation
DIF Diffuse solar radiation
GLB Global solar radiation
Line 110 writes values of parameters used to estimate solar array power and in race simulation.
Lines 112-113 close the input and output files.
Line 114 formats the program output.
Line 115 ends execution of the program.
Below are results of a sample run of the program shown above. This can be thought of as a forecast for a solar car outing from the Detroit area to Big Bay (in the Upper Peninsula) in Michigan. Cumulus clouds form in the afternoon, with a couple of very brief rainshowers in the northern Lower Peninsula, where sea breezes from Lakes Michigan and Huron augment the convection.
Input file: WEA-B.DAT (letters not in actual file)
ID ZH ZL
172 5. 75.
TIME RLON W A0 TR IWIND DP HEAD
RLAT B O3 CS P T RHO ELEV
10. 42.3 83.3 .21 2.55 .28 .075 .00 .00 1007.3 2004 28.2 17.2 .25 323. 173.
11. 42.7 83.9 .21 2.62 .29 .073 .00 .00 1006.7 2005 29.6 17.6 .26 315. 191.
12. 43.2 84.4 .20 2.68 .28 .072 .08 1.13 1005.9 2106 30.8 17.1 .26 320. 157.
13. 43.6 85.0 .20 2.73 .27 .070 .20 .53 1005.2 2107 31.1 16.7 .26 344. 88.
14. 44.2 84.8 .19 2.71 .28 .069 .33 .35 1004.8 2207 30.7 16.2 .26 5. 144.
15. 45.0 84.8 .19 2.89 .29 .068 .89 .23 1003.9 2309 30.0 16.8 .27 11. 161.
16. 45.5 84.5 .18 2.92 .29 .067 .78 .20 1002.6 2310 27.7 17.0 .28 0. 137.
17. 45.9 84.8 .18 2.72 .29 .066 .50 .29 1002.2 2412 26.5 16.5 .28 312. 54.
18. 46.1 85.6 .15 2.58 .29 .066 .22 .32 1003.0 2409 24.7 16.1 .27 281. 97.
19. 46.2 86.5 .14 2.42 .29 .065 .12 .35 1003.7 2506 23.1 15.6 .28 277. 143.
20. 46.3 87.3 .13 2.23 .29 .065 .07 .38 1004.3 3404 21.7 15.0 .29 300. 116.
21. 46.7 87.7 .11 2.17 .29 .065 .03 .46 1004.8 3202 20.4 14.6 .30 339. 10.
Output file: FLUX.DAT
TIME AM ZSUN ASUN RBEAM RISO GLB WN WW T HEAD RHO
10.000 1.2218 .6483 1.2553 843.5 125.0 797.4 -1.9 .7 28.2 5.637 .25
11.000 1.0991 .4862 .9429 873.0 128.7 900.5 -2.4 .9 29.6 5.498 .26
12.000 1.0467 .3729 .4452 821.3 201.0 965.8 -2.7 1.5 30.8 5.585 .26
13.000 1.0480 .3578 -.2093 714.3 201.7 870.8 -3.1 1.8 31.1 6.004 .26
14.000 1.0882 .4591 -.7745 595.5 185.4 719.2 -2.8 2.3 30.7 .087 .26
15.000 1.1935 .6190 -1.1220 94.6 175.7 252.8 -3.0 3.5 30.0 .192 .27
16.000 1.3959 .7998 -1.3723 180.7 128.7 254.6 -3.3 3.9 27.7 .000 .28
17.000 1.7556 .9784 -1.5656 376.3 129.4 339.5 -3.1 5.3 26.5 5.445 .28
18.000 2.3792 1.1492 -1.7348 530.1 98.0 314.9 -2.3 4.0 24.7 4.904 .27
19.000 3.7732 1.3135 -1.8974 464.1 77.5 195.6 -1.1 2.9 23.1 4.835 .28
20.000 8.8259 1.4685 -2.0643 243.6 44.5 69.3 1.9 .7 21.7 5.236 .29
21.000 8.8259 1.6132 .0000 .0 .0 .0 .8 .7 20.4 5.917 .30
Text and embedded images are copyright of Joseph Bartlo, though may be used with proper crediting.