| blob | b5e12825c788a7ec3a000d50e77652066f12ca55 |
1 ! convert zsurf to km
2 ! fix params
3 ! SM NOTES, Aug 2008
4 ! deleted subroutine delted
5 ! deleted subroutine nlayde
6 ! deleted subroutine leqt1b
7 ! added subroutines ps2str, tridiag
8 ! added subroutines sphers, airmas, for future
9 ! use with better Schumman-Runge parameterization
10 ! should be able do delete: nsurf, use 1 instead
11 MODULE module_phot_mad
12 ! preliminary fixed values for T, p, O3 and caer
13 ! use values in correct units!!
14 ! so3t(kl,i) - ozone cross sect. temperature dependence coefficients
15 ! wl(kl) - array of nominal center wavelengths of spectral intervals
16 ! f(kl) - extraterrestrial solar irradiance
17 ! xs(kl,nr) - cross sections for nr species.
18 ! xqy(kl,nr) - quantum yields. some are read, others are computed.
19 ! schumann-runge - part a
20 ! schumann-runge parameters - part b
21 ! ozone temperature coefficients for xsection
22 ! -----------------------------------------------------------------------
23 ! read stand profiles
24 ! aerosol
25 ! close (incvol)
26 ! wavelengths and extraterrestrial irradiance
27 ! -----------------------------------------------------
28 ! wave_length
29 ! WAVE LENGTHS USED BY PHOTOLYSIS PROGRAMS
30 ! FILE CREATED AUGUST 19, 1994
31 ! FROM MADRONICH 1989 DATA FILE
32 ! et_flux
33 ! EXTRA-TERRESTIAL IRRADIANCE
34 ! FILE CREATED AUGUST 19, 1994
35 ! FROM MADRONICH 1989 DATA FILE
37 ! SAMcKeen 2/17/13 Update 296.309 to 305nm flux from SAO (K Chance et al., 2010)
38 ! 2.5% increase in JO1D at 17.34deg SZA, ~1% at 58.8deg SZA, both at 500m AGL
39 ! SAM 2/15/13 Dobson value set to 325., OMI - 5 days of eval. in May 2010 over LA
40 ! 14.2% dobson reduction leads to 25% increase in surf. JO1D at 17deg SZA, ~50% at 80deg SZA
41 !
42 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation
43 ! O3->O1D+O2 quantum yields taken from FTUV function fo3qy, based on
44 ! Matsumi, Y., F. J. Comes, G. Hancock, A. Hofzumahaus, A. J. Hynes,
45 ! M. Kawasaki, and A. R. Ravishankara, Quantum yields for production of O(1D)
46 ! in the ultraviolet photolysis of ozone: Recommendation based on evaluation
47 ! of laboratory data, J. Geophys. Res., 107, 10.1029/2001JD000510, 2002.
48 ! Revised numbers give 23.3% JO1D increase at 17.34deg SZA, ~50% at 58.8deg, both at 500m AGL.
49 ! -----------------------------------------------------------------------
50 ! current jindex assignments - if calculation order changes
51 ! these change - subroutine yield must be changed!
52 ! process jindex notes
53 ! wavelength none center values, modified wmo grid
54 ! e.t.irradian none photons per bin, not per nm
55 ! absorption cross sections:
56 ! o2 absorp 1 schumman-runge corrected in srband
57 ! o3 -> 1d 2 at 275 k. correct t-dep in subgrid
58 ! o3 -> 3p 3 at 275 k. correct t-dep in subgrid
59 ! no2 4
60 ! no3 -> no+o2 5
61 ! no3 -> no2+o 6
62 ! hno2 7
63 ! hno3 8
64 ! hno4 9
65 ! h2o2 10
66 ! ch2o -> rad 11
67 ! ch2o -> mol 12
68 ! ch3cho 13
69 ! ch3coch3 14
70 ! ch3coc2h5 15
71 ! hcocho proc a 16
72 ! ch3cocho 17
73 ! hcoch=chcho 18 estimate, no reliable measurement
74 ! ch3o2h 19
75 ! ch3coo2h 20 actually use 0.28*(h2o2 value)
76 ! ch3ono2 21
77 ! hcocho proc b 22
78 ! quantum yields:
79 ! no2 4
80 ! ch2o -> rad 11
81 ! ch2o -> mol 12
82 ! ch3cho 13
83 ! hcoch=chcho 18 (energetic threshold)
84 ! cross section and quantum yield data. this section assigns
85 ! yields for ntp air. yields are read from data file
86 ! some yields are temperature and/or pressure dependent:ch3cho, ch2o_b,
87 ! o3.
88 ! for ch2o_b and ch3cho, the data read above are stp values. these will
89 ! be
90 ! corrected in subgrid for t & ad dependence, after the altitude
91 ! dependent
92 ! values of t and ad are computed at each layer or level as appropriate.
93 ! the o3->o(1d) yield is also calculated later, in subgrid.
94 ! -----------------------------------------------------
95 ! o2 cross section
96 ! -----------------------------------------------------
97 ! o3 -> o1d cross section
98 ! -----------------------------------------------------
99 ! o3 -> o3p cross section
100 ! -----------------------------------------------------
101 ! no2 cross section
102 ! no2 quantum yield
103 ! -----------------------------------------------------
104 ! no3 -> no + o2 cross section
105 ! no3 -> no + o2 quantum yield
106 ! -----------------------------------------------------
107 ! no3 -> no2 + o cross section
108 ! no3 -> no2 + o quantum yield
109 ! -----------------------------------------------------
110 ! hono cross section
111 ! hono cross quantum yield
112 ! -----------------------------------------------------
113 ! hno3 cross section
114 ! hno3 cross quantum yield
115 ! HNO3 CROSS SECTION TEMPERATURE DEPENDENCE
116 ! -----------------------------------------------------
117 ! hno4 cross section
118 ! hno4 cross quantum yield
119 ! -----------------------------------------------------
120 ! h2o2 cross section
121 ! h2o2 cross quantum yield
122 ! -----------------------------------------------------
123 ! hcho -> ho2 cross section
124 ! hcho -> ho2 quantum yield
125 ! -----------------------------------------------------
126 ! hcho -> h2 cross section
127 ! hcho -> h2 quantum yield
128 ! -----------------------------------------------------
129 ! ch3cho cross section
130 ! ch3cho (ntp) quantum yield
131 ! -----------------------------------------------------
132 ! ch3coch3 cross section
133 ! ch3coch3 quantum yield
134 ! -----------------------------------------------------
135 ! ch3coc2h5 cross section
136 ! ch3coc2h5 quantum yield
137 ! -----------------------------------------------------
138 ! hcocho proc a cross section
139 ! hcocho a quantum yield
140 ! -----------------------------------------------------
141 ! ch3cocho cross section
142 ! ch3cocho quantum yield
143 ! -----------------------------------------------------
144 ! dcb cross section
145 ! dcb quantum yield
146 ! -----------------------------------------------------
147 ! ch3o2h cross section
148 ! ch3o2h cross quantum yield
149 ! -----------------------------------------------------
150 ! ch3coo2h cross section
151 ! ch3coo2h cross quantum yield
152 ! -----------------------------------------------------
153 ! ch3ono2 cross section
154 ! ch3ono2 cross quantum yield
155 ! -----------------------------------------------------
156 ! hcocho proc b cross section
157 ! hcocho b quantum yield
158 ! macr cross section
159 ! macr quantum yield
160 ! .. Parameters ..
161 INTEGER, PARAMETER :: kl0 = 30, kl1 = 130
162 INTEGER, PARAMETER :: kldif = (kl1-kl0+1)*3
163 !liqy. change nreakj from 23 to 27, adding 4 species.
164 INTEGER, PARAMETER :: nabv = 10, nj = 200, nreakj = 27
165 !liqy-20140908
166 INTEGER, PARAMETER :: mj = 2*nj - 2
167 ! .. Local Scalars ..
168 INTEGER :: ip, kl
169 ! .. Local Arrays ..
170 REAL :: aerstd(51), airstd(51), albedoph(130), caabv(nabv), fext(130), &
171 o3abv(nabv), o3std(51), pabv(nabv), jpl295(130), jpl218(130), sra(11,9), &
172 !liqy-change xqy from xqy(130,23) to xqy(130,nreakj)
173 srb(11,5), tabv(nabv), tstd(51), txs(130,nreakj),wl(130),xqy(130,nreakj), &
174 !liqy-20140908
175 xs(130,nreakj), zabv(nabv), zstd(51)
176 ! .. Data Statements ..
177 DATA zabv/21., 22., 23., 24., 25., 30., 35., 40., 45., 50./
178 DATA tabv/215.19, 215.19, 215.19, 215.19, 215.19, 217.39, 227.80, &
179 243.19, 258.50, 265.70/
180 DATA pabv/1.57E18, 1.34E18, 1.14E18, 9.76E17, 8.33E17, 3.83E17, 1.76E17, &
181 8.31E16, 4.09E16, 2.14E16/
182 DATA o3abv/4.88E12, 4.86E12, 4.73E12, 4.54E12, 4.32E12, 2.52E12, &
183 1.40E12, 6.07E11, 2.03E11, 6.61E10/
184 DATA caabv/1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 1.90E-4, &
185 5.00E-5, 1.32E-5, 3.46E-6, 9.14E-7/
186 DATA zstd/0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., &
187 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., &
188 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., &
189 42., 43., 44., 45., 46., 47., 48., 49., 50./
190 DATA ((sra(kl,ip),ip=1,9),kl=1,11)/ -2.158311E+01, -4.164652E-01, &
191 5.266362E-02, 1.655877E-02, 0., 0., 0., 0., 0., -2.184813E+01, &
192 -4.753880E-01, 4.519945E-02, 3.228313E-02, 3.079373E-03, 0., 0., 0., &
193 0., -2.200507E+01, -4.628729E-01, -5.022541E-02, 2.545036E-02, &
194 5.791406E-02, 1.179966E-02, -8.296876E-03, -3.238368E-03, &
195 -3.069686E-04, -2.205527E+01, -4.400848E-01, -5.687308E-03, &
196 3.712279E-02, 6.025527E-03, 0., 0., 0., 0., -2.205261E+01, &
197 -5.707936E-01, -3.330207E-02, 5.959032E-02, 1.510540E-02, &
198 1.000376E-03, 0., 0., 0., -2.228000E+01, -3.960759E-01, -2.995798E-02, &
199 4.918104E-02, 9.269080E-03, -1.173411E-03, -2.599386E-04, 0., 0., &
200 -2.275796E+01, -2.054719E-01, -1.094205E-02, 2.079595E-02, &
201 3.769638E-03, 0., 0., 0., 0., -2.297610E+01, -5.823677E-02, &
202 -1.007612E-01, 2.404666E-02, 4.761876E-02, 4.169606E-03, &
203 -7.126663E-03, -2.263652E-03, -1.971653E-04, -2.506084E+01, &
204 3.442774E-02, -2.212047E-04, 6.186041E-07, -6.284394E-10, 0., 0., 0., &
205 0., -2.313436E+01, 1.177283E-04, 0., 0., 0., 0., 0., 0., 0., &
206 -2.312205E+01, 0., 0., 0., 0., 0., 0., 0., 0./
207 DATA ((srb(kl,ip),ip=1,5),kl=1,11)/ -2.431640E+03, 4.729722E+02, &
208 -3.452121E+01, 1.120677E+00, -1.365618E-02, -3.701955E+01, &
209 3.623290E+00, -8.929223E-02, 0., 0., -1.086239E+03, 1.981847E+02, &
210 -1.359057E+01, 4.155845E-01, -4.788462E-03, -1.213108E+03, &
211 2.277459E+02, -1.612207E+01, 5.101389E-01, -6.090518E-03, &
212 -8.334575E+01, 7.944254E+00, -1.898894E-01, 0., 0., -2.139117E+02, &
213 2.612729E+01, -1.036749E+00, 1.317695E-02, 0., -3.281301E+02, &
214 4.307004E+01, -1.870019E+00, 2.674331E-02, 0., 3.033416E+03, &
215 -5.978911E+02, 4.370384E+01, -1.406715E+00, 1.683967E-02, &
216 -2.535815E+00, 0., 0., 0., 0., -4.474937E+00, 0., 0., 0., 0., &
217 -2.996639E+00, 0., 0., 0., 0./
218 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation
219 DATA (jpl295(kl),kl=1,130)/ 62.2, 57.5844, 52.5359, 47.7, 42.8638, 38.4591, &
220 34.9, 32.4255, 31.5, 32.6204, 36.3, 42.2149, 53.9, 69.3001, 90.3, 118., &
221 154., 199., 255., 322., 401., 490., 590., 693., 802., 908., 1001., 1080., &
222 1125., 1148., 1122., 1064., 968., 840., 698.001, 547., 406., 282., 184., &
223 113., 65.1, 37.2611, 26.2, 23.4, 20.1, 17.9, 15.5, 13.5, 12.2, 10.2, 9.24, &
224 7.95, 6.91, 6.25, 4.66, 2.82778, 1.378, 0.697, 0.32, 0.146, 0.0779, 0.0306, &
225 0.0136, 0.00694, 0.00305, 0.0013, 0.00085, 0.000572, 0.000542, 0.000668, &
226 0.000956, 0.00115, 0.00158, 0.00258, 0.00295, 0.00393, 0.00656, 0.00697, &
227 0.00882, 0.0137, 0.0165, 0.0185, 0.0218, 0.0366, 0.0367, 0.041, 0.0481, &
228 0.0754, 0.0813, 0.0816, 0.0908, 0.121, 0.16, 0.158, 0.166, 0.183, 0.219, &
229 0.267, 0.287, 0.295, 0.319, 0.337, 0.358, 0.398, 0.439, 0.467, 0.481, 0.464, &
230 0.446, 0.447, 0.476, 0.513, 0.514, 0.478, 0.438, 0.406, 0.382, 0.356, 0.327, &
231 0.297, 0.271, 0.245684, 0.21025, 0.17025, 0.13775, 0.112725, 0.087725, 0.07355, &
232 0.062075, 0.048525/
233 DATA (jpl218(kl),kl=1,130)/ 62.2, 57.5844, 52.5359, 47.7, 42.8638, 38.4535, &
234 34.4, 32.0245, 31.2, 32.4203, 36.2, 42.1149, 54.2, 69.6001, 90.6, 119., &
235 155., 201., 256., 323., 403., 492., 589., 692., 799., 905., 995., 1074., &
236 1116., 1136., 1105., 1047., 952., 823., 681.001, 531., 391., 271., 175., &
237 105., 59.4, 33.3104, 22.9, 20.6, 17.3, 15.6, 13.3, 11.5, 10.4, 8.5, 7.76, &
238 6.53, 5.62, 5.05, 3.67, 2.15, 0.9852, 0.483, 0.204, 0.0797, 0.0779, 0.0306, &
239 0.0136, 0.00694, 0.00305, 0.0013, 0.00085, 0.000572, 0.000542, 0.000668, &
240 0.000956, 0.00115, 0.00158, 0.00258, 0.00295, 0.00393, 0.00656, 0.00697, &
241 0.00882, 0.0137, 0.0165, 0.0185, 0.0218, 0.0366, 0.0367, 0.041, 0.0481, &
242 0.0754, 0.0813, 0.0816, 0.0908, 0.121, 0.16, 0.158, 0.166, 0.183, 0.219, &
243 0.267, 0.287, 0.295, 0.319, 0.337, 0.358, 0.398, 0.439, 0.467, 0.481, 0.464, &
244 0.446, 0.447, 0.476, 0.513, 0.514, 0.478, 0.438, 0.406, 0.382, 0.356, 0.327, &
245 0.297, 0.271, 0.245684, 0.21025, 0.17025, 0.13775, 0.112725, 0.087725, 0.07355, &
246 0.062075, 0.048525/
247 DATA aerstd/2.40E-1, 1.06E-1, 4.56E-2, 1.91E-2, 1.01E-2, 7.63E-3, &
248 5.38E-3, 5.00E-3, 5.15E-3, 4.94E-3, 4.82E-3, 4.51E-3, 4.74E-3, &
249 4.37E-3, 4.28E-3, 4.03E-3, 3.83E-3, 3.78E-3, 3.88E-3, 3.08E-3, &
250 2.26E-3, 1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 5.50E-4, &
251 4.21E-4, 3.22E-4, 2.48E-4, 1.90E-4, 1.45E-4, 1.11E-4, 8.51E-5, &
252 6.52E-5, 5.00E-5, 3.83E-5, 2.93E-5, 2.25E-5, 1.72E-5, 1.32E-5, &
253 1.01E-5, 7.72E-6, 5.91E-6, 4.53E-6, 3.46E-6, 2.66E-6, 2.04E-6, &
254 1.56E-6, 1.19E-6, 9.14E-7/
255 DATA (wl(kl),kl=1,130)/1.861E+02, 1.878E+02, 1.896E+02, 1.914E+02, &
256 1.933E+02, 1.952E+02, 1.971E+02, 1.990E+02, 2.010E+02, 2.031E+02, &
257 2.052E+02, 2.073E+02, 2.094E+02, 2.117E+02, 2.139E+02, 2.162E+02, &
258 2.186E+02, 2.210E+02, 2.235E+02, 2.260E+02, 2.286E+02, 2.313E+02, &
259 2.340E+02, 2.367E+02, 2.396E+02, 2.425E+02, 2.454E+02, 2.485E+02, &
260 2.516E+02, 2.548E+02, 2.582E+02, 2.615E+02, 2.650E+02, 2.685E+02, &
261 2.722E+02, 2.759E+02, 2.798E+02, 2.837E+02, 2.878E+02, 2.920E+02, &
262 2.963E+02, 3.005E+02, 3.030E+02, 3.040E+02, 3.050E+02, 3.060E+02, &
263 3.070E+02, 3.080E+02, 3.090E+02, 3.100E+02, 3.110E+02, 3.120E+02, &
264 3.130E+02, 3.140E+02, 3.160E+02, 3.200E+02, 3.250E+02, 3.300E+02, &
265 3.350E+02, 3.400E+02, 3.450E+02, 3.500E+02, 3.550E+02, 3.600E+02, &
266 3.650E+02, 3.700E+02, 3.750E+02, 3.800E+02, 3.850E+02, 3.900E+02, &
267 3.950E+02, 4.000E+02, 4.050E+02, 4.100E+02, 4.150E+02, 4.200E+02, &
268 4.250E+02, 4.300E+02, 4.350E+02, 4.400E+02, 4.450E+02, 4.500E+02, &
269 4.550E+02, 4.600E+02, 4.650E+02, 4.700E+02, 4.750E+02, 4.800E+02, &
270 4.850E+02, 4.900E+02, 4.950E+02, 5.000E+02, 5.050E+02, 5.100E+02, &
271 5.150E+02, 5.200E+02, 5.250E+02, 5.300E+02, 5.350E+02, 5.400E+02, &
272 5.450E+02, 5.500E+02, 5.550E+02, 5.600E+02, 5.650E+02, 5.700E+02, &
273 5.750E+02, 5.800E+02, 5.850E+02, 5.900E+02, 5.950E+02, 6.000E+02, &
274 6.050E+02, 6.100E+02, 6.150E+02, 6.200E+02, 6.250E+02, 6.300E+02, &
275 6.350E+02, 6.400E+02, 6.448E+02, 6.510E+02, 6.600E+02, 6.700E+02, &
276 6.800E+02, 6.900E+02, 7.000E+02, 7.100E+02, 7.200E+02, 7.300E+02/
277 DATA (fext(kl),kl=1,130)/3.620E+11, 4.730E+11, 5.610E+11, 6.630E+11, &
278 6.900E+11, 9.560E+11, 1.150E+12, 1.270E+12, 1.520E+12, 1.780E+12, &
279 2.200E+12, 2.690E+12, 4.540E+12, 7.140E+12, 8.350E+12, 8.390E+12, &
280 1.080E+13, 1.180E+13, 1.600E+13, 1.340E+13, 1.410E+13, 1.570E+13, &
281 1.380E+13, 1.600E+13, 1.450E+13, 2.200E+13, 1.990E+13, 1.970E+13, &
282 1.940E+13, 2.910E+13, 4.950E+13, 4.530E+13, 1.070E+14, 1.200E+14, &
283 1.100E+14, 1.040E+14, 8.240E+13, 1.520E+14, 2.150E+14, 3.480E+14, &
284 3.790E+14, 2.990E+14, 1.003E+14, 9.294E+13, 1.024E+14, 8.507E+13, &
285 9.383E+13, 1.030E+14, 9.722E+13, 7.751E+13, 1.277E+14, 1.087E+14, &
286 1.102E+14, 1.184E+14, 3.153E+14, 5.930E+14, 6.950E+14, 8.150E+14, &
287 7.810E+14, 8.350E+14, 8.140E+14, 8.530E+14, 9.170E+14, 8.380E+14, &
288 1.040E+15, 1.100E+15, 9.790E+14, 1.130E+15, 8.890E+14, 1.140E+15, &
289 9.170E+14, 1.690E+15, 1.700E+15, 1.840E+15, 1.870E+15, 1.950E+15, &
290 1.810E+15, 1.670E+15, 1.980E+15, 2.020E+15, 2.180E+15, 2.360E+15, &
291 2.310E+15, 2.390E+15, 2.380E+15, 2.390E+15, 2.440E+15, 2.510E+15, &
292 2.300E+15, 2.390E+15, 2.480E+15, 2.400E+15, 2.460E+15, 2.490E+15, &
293 2.320E+15, 2.390E+15, 2.420E+15, 2.550E+15, 2.510E+15, 2.490E+15, &
294 2.550E+15, 2.530E+15, 2.540E+15, 2.500E+15, 2.570E+15, 2.580E+15, &
295 2.670E+15, 2.670E+15, 2.700E+15, 2.620E+15, 2.690E+15, 2.630E+15, &
296 2.680E+15, 2.660E+15, 2.590E+15, 2.690E+15, 2.610E+15, 2.620E+15, &
297 2.620E+15, 2.630E+15, 2.392E+15, 3.998E+15, 5.115E+15, 5.225E+15, &
298 5.215E+15, 5.105E+15, 5.140E+15, 5.010E+15, 4.930E+15, 4.895E+15/
299 DATA (xs(kl,1),kl=1,130)/7.040E-24, 7.360E-24, 7.640E-24, 7.870E-24, &
300 8.040E-24, 8.140E-24, 8.170E-24, 8.130E-24, 8.010E-24, 7.840E-24, &
301 7.630E-24, 7.330E-24, 6.990E-24, 6.450E-24, 5.810E-24, 5.230E-24, &
302 4.710E-24, 4.260E-24, 3.800E-24, 3.350E-24, 2.900E-24, 2.450E-24, &
303 2.050E-24, 1.690E-24, 1.300E-24, 9.300E-25, 0., 0., 0., 0., 0., 0., &
304 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
305 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
306 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
307 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
308 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
309 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0./
310 DATA (xs(kl,2),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
311 0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
312 0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
313 0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
314 0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
315 0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
316 0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
317 0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
318 0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
319 0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
320 0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
321 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
322 0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
323 0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
324 0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
325 0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
326 0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
327 0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
328 0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
329 0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
330 0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
331 0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
332 DATA (xs(kl,3),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
333 0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
334 0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
335 0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
336 0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
337 0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
338 0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
339 0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
340 0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
341 0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
342 0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
343 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
344 0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
345 0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
346 0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
347 0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
348 0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
349 0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
350 0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
351 0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
352 0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
353 0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
354 DATA (xs(kl,4),kl=1,130)/0.259E-18, 0.272E-18, 0.286E-18, 0.273E-18, &
355 0.251E-18, 0.244E-18, 0.246E-18, 0.246E-18, 0.282E-18, 0.415E-18, &
356 0.448E-18, 0.445E-18, 0.464E-18, 0.487E-18, 0.482E-18, 0.502E-18, &
357 0.444E-18, 0.471E-18, 0.377E-18, 0.393E-18, 0.274E-18, 0.278E-18, &
358 0.169E-18, 0.162E-18, 0.882E-19, 0.747E-19, 0.391E-19, 0.275E-19, &
359 0.201E-19, 0.197E-19, 0.211E-19, 0.236E-19, 0.270E-19, 0.325E-19, &
360 0.379E-19, 0.503E-19, 0.588E-19, 0.700E-19, 0.815E-19, 0.972E-19, &
361 0.115E-18, 0.128E-18, 0.154E-18, 0.159E-18, 0.158E-18, 0.156E-18, &
362 0.164E-18, 0.166E-18, 0.182E-18, 0.184E-18, 0.192E-18, 0.204E-18, &
363 0.204E-18, 0.202E-18, 0.224E-18, 0.248E-18, 0.281E-18, 0.313E-18, &
364 0.343E-18, 0.380E-18, 0.407E-18, 0.431E-18, 0.472E-18, 0.483E-18, &
365 0.517E-18, 0.532E-18, 0.551E-18, 0.564E-18, 0.576E-18, 0.593E-18, &
366 0.585E-18, 0.602E-18, 0.578E-18, 0.600E-18, 0.565E-18, 0.581E-18, &
367 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
368 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
369 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
370 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
371 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
372 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
373 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
374 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
375 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
376 DATA ((xqy(kl,ip),kl=kl0,kl1),ip=1,3)/kldif*1./
377 DATA (xqy(kl,4),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
378 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
379 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
380 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
381 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
382 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
383 0.999000, 0.998000, 0.997000, 0.996500, 0.996000, 0.996000, 0.996000, &
384 0.996000, 0.995000, 0.995000, 0.995000, 0.995000, 0.995000, 0.994000, &
385 0.994000, 0.994000, 0.993000, 0.992000, 0.991000, 0.990000, 0.989000, &
386 0.988000, 0.987000, 0.986000, 0.984000, 0.983000, 0.981000, 0.979000, &
387 0.975000, 0.969000, 0.960000, 0.927000, 0.694000, 0.355000, 0.134000, &
388 0.060000, 0.018000, 0.000900, 0.000000, 0.000000, 0.000000, 0.000000, &
389 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
390 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
391 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
392 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
393 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
394 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
395 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
396 DATA (xs(kl,5),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
397 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
398 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
399 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
400 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
401 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
402 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
403 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
404 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
405 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
406 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
407 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
408 0.000E+00, 0.000E+00, 0.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
409 0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
410 0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
411 0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
412 0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
413 0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
414 0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
415 0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
416 0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
417 0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
418 DATA (xqy(kl,5),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
419 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
420 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
421 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
422 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
423 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
424 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
425 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
426 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
427 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
428 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
429 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
430 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
431 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
432 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
433 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.038000, &
434 0.191000, 0.326000, 0.311000, 0.272000, 0.233000, 0.194000, 0.156000, &
435 0.117000, 0.078000, 0.039000, 0.000000, 0.000000, 0.000000, 0.000000, &
436 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
437 DATA (xs(kl,6),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
438 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
439 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
440 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
441 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
442 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
443 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
444 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
445 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
446 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
447 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
448 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
449 0.000E+00, 0.000E+00, 0.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
450 0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
451 0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
452 0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
453 0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
454 0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
455 0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
456 0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
457 0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
458 0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
459 DATA (xqy(kl,6),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
460 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
461 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
462 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
463 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
464 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
465 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
466 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
467 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
468 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
469 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
470 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
471 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
472 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
473 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
474 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.962000, &
475 0.809000, 0.661000, 0.578000, 0.506000, 0.433000, 0.361000, 0.289000, &
476 0.217000, 0.144000, 0.072000, 0.000000, 0.000000, 0.000000, 0.000000, &
477 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
478 DATA (xs(kl,7),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
479 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
480 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
481 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
482 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
483 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
484 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
485 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
486 0.000E+00, 0.000E+00, 0.000E+00, 0.130E-19, 0.190E-19, 0.280E-19, &
487 0.220E-19, 0.360E-19, 0.340E-19, 0.536E-19, 0.534E-19, 0.111E-18, &
488 0.786E-19, 0.189E-18, 0.116E-18, 0.130E-18, 0.279E-18, 0.954E-19, &
489 0.179E-18, 0.260E-18, 0.590E-19, 0.101E-18, 0.176E-18, 0.304E-19, &
490 0.775E-20, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
491 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
492 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
493 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
494 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
495 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
496 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
497 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
498 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
499 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
500 DATA (xqy(kl,7),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
501 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
502 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
503 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
504 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
505 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
506 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
507 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
508 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
509 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
510 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
511 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
512 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
513 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
514 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
515 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
516 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
517 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
518 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
519 DATA (xs(kl,8),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.127E-16, &
520 0.114E-16, 0.100E-16, 0.847E-17, 0.679E-17, 0.518E-17, 0.382E-17, &
521 0.270E-17, 0.182E-17, 0.120E-17, 0.730E-18, 0.451E-18, 0.283E-18, &
522 0.195E-18, 0.134E-18, 0.102E-18, 0.802E-19, 0.650E-19, 0.518E-19, &
523 0.414E-19, 0.321E-19, 0.265E-19, 0.230E-19, 0.209E-19, 0.199E-19, &
524 0.196E-19, 0.195E-19, 0.193E-19, 0.188E-19, 0.180E-19, 0.170E-19, &
525 0.152E-19, 0.134E-19, 0.113E-19, 0.924E-20, 0.719E-20, 0.532E-20, &
526 0.371E-20, 0.249E-20, 0.188E-20, 0.167E-20, 0.150E-20, 0.133E-20, &
527 0.119E-20, 0.105E-20, 0.932E-21, 0.814E-21, 0.721E-21, 0.628E-21, &
528 0.547E-21, 0.465E-21, 0.362E-21, 0.197E-21, 0.975E-22, 0.452E-22, &
529 0.222E-22, 0.110E-22, 0.604E-23, 0.420E-23, 0.000E+00, 0.000E+00, &
530 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
531 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
532 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
533 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
534 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
535 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
536 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
537 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
538 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
539 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
540 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
541 DATA (xqy(kl,8),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
542 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
543 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
544 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
545 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
546 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
547 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
548 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
549 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
550 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
551 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
552 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
553 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
554 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
555 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
556 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
557 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
558 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
559 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
560 DATA (txs(kl,8),kl=1,130)/0.000000, 0.000000, 0.000000, 1.700000, &
561 1.650000, 1.660000, 1.690000, 1.740000, 1.770000, 1.850000, 1.970000, &
562 2.080000, 2.170000, 2.170000, 2.210000, 2.150000, 2.060000, 1.960000, &
563 1.840000, 1.780000, 1.800000, 1.860000, 1.900000, 1.970000, 1.970000, &
564 1.970000, 1.880000, 1.750000, 1.610000, 1.440000, 1.340000, 1.230000, &
565 1.180000, 1.140000, 1.120000, 1.140000, 1.140000, 1.180000, 1.220000, &
566 1.250000, 1.450000, 1.490000, 1.560000, 1.640000, 1.690000, 1.780000, &
567 1.870000, 1.940000, 2.040000, 2.150000, 2.270000, 2.380000, 2.620000, &
568 2.700000, 2.920000, 3.100000, 3.240000, 3.520000, 3.770000, 3.910000, &
569 4.230000, 4.700000, 5.150000, 5.250000, 5.740000, 6.450000, 6.700000, &
570 7.160000, 7.550000, 8.160000, 9.750000, 9.930000, 9.600000, 10.50000, &
571 10.80000, 11.80000, 11.80000, 9.300000, 12.10000, 11.90000, 9.300000, &
572 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
573 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
574 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
575 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
576 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
577 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
578 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
579 DATA (xs(kl,9),kl=1,130)/0.000E+00, 0.000E+00, 0.100E-16, 0.980E-17, &
580 0.918E-17, 0.828E-17, 0.723E-17, 0.618E-17, 0.517E-17, 0.426E-17, &
581 0.352E-17, 0.292E-17, 0.245E-17, 0.205E-17, 0.175E-17, 0.151E-17, &
582 0.131E-17, 0.115E-17, 0.102E-17, 0.916E-18, 0.827E-18, 0.752E-18, &
583 0.687E-18, 0.631E-18, 0.578E-18, 0.529E-18, 0.484E-18, 0.439E-18, &
584 0.396E-18, 0.353E-18, 0.311E-18, 0.271E-18, 0.231E-18, 0.194E-18, &
585 0.158E-18, 0.125E-18, 0.946E-19, 0.694E-19, 0.485E-19, 0.325E-19, &
586 0.210E-19, 0.135E-19, 0.104E-19, 0.938E-20, 0.846E-20, 0.763E-20, &
587 0.689E-20, 0.623E-20, 0.564E-20, 0.511E-20, 0.463E-20, 0.420E-20, &
588 0.382E-20, 0.347E-20, 0.287E-20, 0.191E-20, 0.101E-20, 0.000E+00, &
589 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
590 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
591 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
592 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
593 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
594 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
595 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
596 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
597 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
598 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
599 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
600 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
601 DATA (xqy(kl,9),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
602 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
603 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
604 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
605 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
606 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
607 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
608 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
609 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
610 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
611 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
612 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
613 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
614 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
615 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
616 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
617 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
618 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
619 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
620 DATA (xs(kl,10),kl=1,130)/0.000E+00, 0.000E+00, 0.695E-18, 0.635E-18, &
621 0.586E-18, 0.547E-18, 0.515E-18, 0.488E-18, 0.462E-18, 0.436E-18, &
622 0.412E-18, 0.388E-18, 0.365E-18, 0.341E-18, 0.318E-18, 0.295E-18, &
623 0.272E-18, 0.250E-18, 0.229E-18, 0.209E-18, 0.190E-18, 0.172E-18, &
624 0.155E-18, 0.140E-18, 0.126E-18, 0.112E-18, 0.999E-19, 0.881E-19, &
625 0.774E-19, 0.675E-19, 0.582E-19, 0.500E-19, 0.425E-19, 0.358E-19, &
626 0.298E-19, 0.247E-19, 0.201E-19, 0.164E-19, 0.131E-19, 0.105E-19, &
627 0.828E-20, 0.658E-20, 0.573E-20, 0.543E-20, 0.514E-20, 0.486E-20, &
628 0.460E-20, 0.435E-20, 0.412E-20, 0.390E-20, 0.369E-20, 0.349E-20, &
629 0.330E-20, 0.312E-20, 0.279E-20, 0.220E-20, 0.160E-20, 0.130E-20, &
630 0.100E-20, 0.700E-21, 0.500E-21, 0.400E-21, 0.000E+00, 0.000E+00, &
631 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
632 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
633 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
634 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
635 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
636 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
637 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
638 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
639 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
640 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
641 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
642 DATA (xqy(kl,10),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
643 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
644 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
645 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
646 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
647 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
648 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
649 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
650 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
651 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
652 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
653 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
654 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
655 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
656 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
657 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
658 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
659 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
660 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
661 DATA (xs(kl,11),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
662 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
663 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
664 0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
665 0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
666 0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
667 0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
668 0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
669 0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
670 0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
671 0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 0.000E+00, &
672 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
673 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
674 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
675 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
676 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
677 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
678 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
679 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
680 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
681 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
682 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
683 DATA (xqy(kl,11),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
684 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
685 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
686 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.264091, &
687 0.288940, 0.297038, 0.295059, 0.289384, 0.285483, 0.288087, 0.301000, &
688 0.326819, 0.366764, 0.420506, 0.485961, 0.559106, 0.633887, 0.702103, &
689 0.733457, 0.740762, 0.747845, 0.753000, 0.754000, 0.754800, 0.754000, &
690 0.753000, 0.752000, 0.751000, 0.749500, 0.745000, 0.739600, 0.731700, &
691 0.723300, 0.690300, 0.593100, 0.458100, 0.305000, 0.122300, 0.003429, &
692 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
693 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
694 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
695 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
696 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
697 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
698 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
699 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
700 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
701 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
702 DATA (xs(kl,12),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
703 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
704 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
705 0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
706 0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
707 0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
708 0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
709 0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
710 0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
711 0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
712 0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 0.000E+00, &
713 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
714 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
715 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
716 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
717 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
718 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
719 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
720 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
721 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
722 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
723 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
724 DATA (xqy(kl,12),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
725 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
726 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
727 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.492085, &
728 0.483681, 0.483325, 0.487471, 0.492514, 0.495532, 0.493893, 0.485473, &
729 0.468839, 0.443373, 0.409405, 0.368400, 0.323132, 0.307820, 0.294564, &
730 0.280920, 0.266885, 0.253277, 0.249000, 0.247000, 0.245600, 0.248000, &
731 0.251000, 0.254000, 0.257000, 0.260200, 0.264500, 0.269000, 0.273500, &
732 0.278900, 0.310300, 0.394100, 0.508100, 0.676100, 0.759300, 0.636100, &
733 0.501500, 0.373400, 0.229000, 0.103600, 0.005906, 0.000000, 0.000000, &
734 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
735 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
736 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
737 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
738 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
739 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
740 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
741 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
742 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
743 DATA (xs(kl,13),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
744 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.593E-21, 0.548E-21, &
745 0.552E-21, 0.513E-21, 0.480E-21, 0.482E-21, 0.482E-21, 0.482E-21, &
746 0.536E-21, 0.593E-21, 0.734E-21, 0.948E-21, 0.125E-20, 0.171E-20, &
747 0.234E-20, 0.321E-20, 0.434E-20, 0.582E-20, 0.770E-20, 0.991E-20, &
748 0.127E-19, 0.159E-19, 0.200E-19, 0.237E-19, 0.287E-19, 0.326E-19, &
749 0.376E-19, 0.408E-19, 0.444E-19, 0.463E-19, 0.466E-19, 0.465E-19, &
750 0.432E-19, 0.406E-19, 0.372E-19, 0.348E-19, 0.342E-19, 0.342E-19, &
751 0.336E-19, 0.333E-19, 0.314E-19, 0.293E-19, 0.276E-19, 0.253E-19, &
752 0.247E-19, 0.243E-19, 0.210E-19, 0.169E-19, 0.108E-19, 0.651E-20, &
753 0.314E-20, 0.138E-20, 0.224E-21, 0.947E-22, 0.425E-22, 0.229E-22, &
754 0.275E-23, 0.192E-23, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
755 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
756 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
757 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
758 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
759 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
760 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
761 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
762 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
763 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
764 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
765 DATA (xqy(kl,13),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
766 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
767 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
768 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
769 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.310000, &
770 0.349300, 0.377900, 0.430300, 0.501600, 0.566200, 0.561500, 0.541100, &
771 0.512100, 0.473200, 0.430000, 0.392000, 0.376000, 0.360000, 0.344000, &
772 0.328000, 0.312000, 0.296000, 0.279800, 0.262500, 0.245000, 0.227500, &
773 0.210000, 0.175000, 0.109300, 0.051880, 0.006266, 0.000000, 0.000000, &
774 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
775 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
776 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
777 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
778 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
779 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
780 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
781 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
782 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
783 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
784 DATA (xs(kl,14),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
785 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.003E-20, 1.873E-21, &
786 1.408E-21, 1.084E-21, 1.030E-21, 1.071E-21, 1.176E-21, 1.371E-21, &
787 1.690E-21, 2.153E-21, 2.779E-21, 3.561E-21, 4.575E-21, 5.906E-21, &
788 7.585E-21, 9.622E-21, 1.212E-20, 1.516E-20, 1.863E-20, 2.252E-20, &
789 2.676E-20, 3.143E-20, 3.616E-20, 4.058E-20, 4.475E-20, 4.774E-20, &
790 5.029E-20, 5.065E-20, 5.049E-20, 4.790E-20, 4.415E-20, 3.940E-20, &
791 3.308E-20, 2.717E-20, 2.371E-20, 2.244E-20, 2.106E-20, 1.953E-20, &
792 1.801E-20, 1.663E-20, 1.538E-20, 1.408E-20, 1.277E-20, 1.173E-20, &
793 1.081E-20, 9.675E-21, 7.783E-21, 4.712E-21, 2.078E-21, 7.065E-22, &
794 1.973E-22, 5.745E-23, 2.137E-23, 8.235E-24, 5.686E-24, 2.157E-24, &
795 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
796 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
797 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
798 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
799 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
800 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
801 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
802 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
803 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
804 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
805 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
806 DATA (xqy(kl,14),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
807 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
808 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
809 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
810 0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
811 0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
812 0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
813 0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
814 0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 0.000000, 0.000000, &
815 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
816 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
817 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
818 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
819 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
820 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
821 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
822 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
823 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
824 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
825 DATA (xs(kl,15),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
826 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.234E-19, 0.685E-20, &
827 0.215E-20, 0.168E-20, 0.159E-20, 0.164E-20, 0.181E-20, 0.203E-20, &
828 0.230E-20, 0.268E-20, 0.320E-20, 0.387E-20, 0.471E-20, 0.582E-20, &
829 0.728E-20, 0.913E-20, 0.115E-19, 0.145E-19, 0.180E-19, 0.222E-19, &
830 0.267E-19, 0.321E-19, 0.374E-19, 0.430E-19, 0.479E-19, 0.525E-19, &
831 0.555E-19, 0.576E-19, 0.575E-19, 0.555E-19, 0.518E-19, 0.460E-19, &
832 0.388E-19, 0.319E-19, 0.269E-19, 0.251E-19, 0.233E-19, 0.217E-19, &
833 0.202E-19, 0.188E-19, 0.173E-19, 0.158E-19, 0.142E-19, 0.128E-19, &
834 0.114E-19, 0.101E-19, 0.796E-20, 0.463E-20, 0.196E-20, 0.705E-21, &
835 0.207E-21, 0.545E-22, 0.145E-22, 0.431E-23, 0.471E-23, 0.157E-23, &
836 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
837 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
838 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
839 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
840 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
841 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
842 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
843 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
844 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
845 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
846 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
847 DATA (xqy(kl,15),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
848 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
849 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
850 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
851 0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
852 0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
853 0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
854 0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
855 0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 0.000000, 0.000000, &
856 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
857 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
858 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
859 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
860 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
861 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
862 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
863 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
864 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
865 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
866 DATA (xs(kl,16),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
867 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
868 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
869 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
870 0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
871 0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
872 0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
873 0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
874 0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
875 0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
876 0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
877 0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
878 0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
879 0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
880 0.201E-18, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
881 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
882 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
883 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
884 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
885 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
886 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
887 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
888 DATA (xqy(kl,16),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
889 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
890 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
891 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
892 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
893 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
894 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
895 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
896 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
897 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
898 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
899 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
900 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
901 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
902 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
903 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
904 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
905 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
906 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000/
907 DATA (xs(kl,17),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
908 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
909 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
910 0.000E+00, 0.000E+00, 0.000E+00, 0.131E-19, 0.142E-19, 0.156E-19, &
911 0.174E-19, 0.189E-19, 0.205E-19, 0.219E-19, 0.233E-19, 0.252E-19, &
912 0.269E-19, 0.285E-19, 0.313E-19, 0.338E-19, 0.362E-19, 0.394E-19, &
913 0.427E-19, 0.450E-19, 0.486E-19, 0.476E-19, 0.479E-19, 0.465E-19, &
914 0.420E-19, 0.371E-19, 0.352E-19, 0.344E-19, 0.336E-19, 0.316E-19, &
915 0.296E-19, 0.276E-19, 0.256E-19, 0.237E-19, 0.227E-19, 0.218E-19, &
916 0.208E-19, 0.199E-19, 0.182E-19, 0.151E-19, 0.938E-20, 0.652E-20, &
917 0.482E-20, 0.323E-20, 0.300E-20, 0.394E-20, 0.560E-20, 0.695E-20, &
918 0.108E-19, 0.148E-19, 0.191E-19, 0.243E-19, 0.322E-19, 0.403E-19, &
919 0.473E-19, 0.566E-19, 0.692E-19, 0.846E-19, 0.968E-19, 0.103E-18, &
920 0.102E-18, 0.101E-18, 0.106E-18, 0.104E-18, 0.994E-19, 0.813E-19, &
921 0.395E-19, 0.109E-19, 0.327E-20, 0.000E+00, 0.000E+00, 0.000E+00, &
922 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
923 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
924 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
925 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
926 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
927 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
928 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
929 DATA (xqy(kl,17),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
930 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
931 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
932 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
933 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
934 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
935 1.000000, 1.000000, 1.000000, 1.000000, 0.990000, 0.990000, 0.990000, &
936 0.980000, 0.980000, 0.980000, 0.980000, 0.970000, 0.970000, 0.970000, &
937 0.960000, 0.960000, 0.940000, 0.920000, 0.880000, 0.825000, 0.750000, &
938 0.660000, 0.560000, 0.480000, 0.400000, 0.320000, 0.250000, 0.200000, &
939 0.150000, 0.120000, 0.100000, 0.080000, 0.060000, 0.050000, 0.040000, &
940 0.030000, 0.020000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, &
941 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.000000, &
942 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
943 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
944 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
945 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
946 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
947 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
948 DATA (xs(kl,18),kl=1,130)/0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
949 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
950 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
951 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
952 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
953 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
954 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
955 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
956 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
957 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
958 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
959 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
960 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
961 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
962 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
963 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
964 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
965 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
966 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
967 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
968 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
969 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
970 DATA (xqy(kl,18),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
971 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
972 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
973 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
974 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
975 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
976 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
977 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
978 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
979 1.000000, 1.000000, 0.500000, 0.000000, 0.000000, 0.000000, 0.000000, &
980 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
981 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
982 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
983 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
984 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
985 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
986 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
987 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
988 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
989 DATA (xs(kl,19),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
990 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
991 0.000E+00, 0.000E+00, 0.320E-18, 0.268E-18, 0.226E-18, 0.193E-18, &
992 0.167E-18, 0.147E-18, 0.129E-18, 0.115E-18, 0.102E-18, 0.899E-19, &
993 0.797E-19, 0.708E-19, 0.623E-19, 0.548E-19, 0.483E-19, 0.422E-19, &
994 0.369E-19, 0.321E-19, 0.278E-19, 0.242E-19, 0.209E-19, 0.180E-19, &
995 0.154E-19, 0.131E-19, 0.111E-19, 0.925E-20, 0.763E-20, 0.622E-20, &
996 0.501E-20, 0.402E-20, 0.352E-20, 0.333E-20, 0.316E-20, 0.299E-20, &
997 0.283E-20, 0.268E-20, 0.254E-20, 0.240E-20, 0.227E-20, 0.215E-20, &
998 0.204E-20, 0.193E-20, 0.172E-20, 0.138E-20, 0.105E-20, 0.801E-21, &
999 0.612E-21, 0.467E-21, 0.356E-21, 0.270E-21, 0.206E-21, 0.160E-21, &
1000 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1001 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1002 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1003 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1004 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1005 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1006 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1007 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1008 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1009 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1010 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1011 DATA (xqy(kl,19),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
1012 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1013 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1014 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1015 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1016 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1017 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1018 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1019 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1020 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1021 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1022 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1023 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1024 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1025 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1026 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1027 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1028 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1029 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1030 DATA (xs(kl,20),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.898E-19, &
1031 0.841E-19, 0.784E-19, 0.148E-18, 0.138E-18, 0.129E-18, 0.122E-18, &
1032 0.114E-18, 0.107E-18, 0.998E-19, 0.932E-19, 0.868E-19, 0.807E-19, &
1033 0.748E-19, 0.690E-19, 0.633E-19, 0.579E-19, 0.529E-19, 0.480E-19, &
1034 0.433E-19, 0.390E-19, 0.349E-19, 0.313E-19, 0.278E-19, 0.247E-19, &
1035 0.217E-19, 0.190E-19, 0.162E-19, 0.138E-19, 0.118E-19, 0.991E-20, &
1036 0.829E-20, 0.690E-20, 0.566E-20, 0.452E-20, 0.361E-20, 0.288E-20, &
1037 0.229E-20, 0.181E-20, 0.157E-20, 0.147E-20, 0.138E-20, 0.131E-20, &
1038 0.125E-20, 0.118E-20, 0.112E-20, 0.106E-20, 0.100E-20, 0.947E-21, &
1039 0.893E-21, 0.838E-21, 0.743E-21, 0.581E-21, 0.428E-21, 0.332E-21, &
1040 0.262E-21, 0.192E-21, 0.141E-21, 0.910E-22, 0.000E+00, 0.000E+00, &
1041 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1042 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1043 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1044 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1045 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1046 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1047 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1048 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1049 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1050 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1051 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1052 DATA (xqy(kl,20),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
1053 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1054 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1055 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1056 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1057 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1058 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1059 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1060 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1061 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1062 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1063 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1064 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1065 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1066 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1067 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1068 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1069 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1070 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1071 DATA (xs(kl,21),kl=1,130)/0.180E-16, 0.181E-16, 0.179E-16, 0.174E-16, &
1072 0.167E-16, 0.159E-16, 0.146E-16, 0.133E-16, 0.118E-16, 0.101E-16, &
1073 0.850E-17, 0.695E-17, 0.540E-17, 0.411E-17, 0.301E-17, 0.215E-17, &
1074 0.163E-17, 0.105E-17, 0.754E-18, 0.524E-18, 0.382E-18, 0.272E-18, &
1075 0.202E-18, 0.152E-18, 0.112E-18, 0.876E-19, 0.677E-19, 0.596E-19, &
1076 0.541E-19, 0.510E-19, 0.489E-19, 0.469E-19, 0.448E-19, 0.419E-19, &
1077 0.377E-19, 0.336E-19, 0.285E-19, 0.238E-19, 0.190E-19, 0.145E-19, &
1078 0.106E-19, 0.752E-20, 0.610E-20, 0.553E-20, 0.496E-20, 0.457E-20, &
1079 0.418E-20, 0.380E-20, 0.341E-20, 0.302E-20, 0.279E-20, 0.256E-20, &
1080 0.232E-20, 0.209E-20, 0.171E-20, 0.110E-20, 0.644E-21, 0.416E-21, &
1081 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1082 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1083 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1084 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1085 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1086 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1087 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1088 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1089 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1090 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1091 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1092 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1093 DATA (xqy(kl,21),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
1094 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1095 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1096 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1097 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1098 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1099 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1100 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1101 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1102 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1103 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1104 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1105 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1106 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1107 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1108 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1109 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1110 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1111 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1112 DATA (xs(kl,22),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1113 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1114 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1115 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1116 0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
1117 0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
1118 0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
1119 0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
1120 0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
1121 0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
1122 0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
1123 0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
1124 0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
1125 0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
1126 0.201E-18, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1127 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1128 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1129 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1130 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1131 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1132 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1133 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1134 DATA (xqy(kl,22),kl=1,130)/0.400000, 0.400000, 0.400000, 0.400000, &
1135 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1136 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1137 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1138 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1139 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1140 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1141 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1142 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1143 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1144 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1145 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1146 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1147 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1148 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1149 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1150 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1151 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1152 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1153 DATA (xs(kl,23),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1154 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1155 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1156 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1157 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.164E-20, &
1158 0.211E-20, 0.227E-20, 0.264E-20, 0.340E-20, 0.465E-20, 0.598E-20, &
1159 0.803E-20, 0.986E-20, 0.118E-19, 0.137E-19, 0.160E-19, 0.189E-19, &
1160 0.228E-19, 0.275E-19, 0.307E-19, 0.321E-19, 0.335E-19, 0.349E-19, &
1161 0.363E-19, 0.378E-19, 0.393E-19, 0.407E-19, 0.421E-19, 0.435E-19, &
1162 0.449E-19, 0.462E-19, 0.487E-19, 0.527E-19, 0.564E-19, 0.589E-19, &
1163 0.616E-19, 0.556E-19, 0.553E-19, 0.543E-19, 0.365E-19, 0.318E-19, &
1164 0.316E-19, 0.156E-19, 0.428E-20, 0.113E-20, 0.000E+00, 0.000E+00, &
1165 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1166 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1167 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1168 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1169 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1170 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1171 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1172 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1173 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1174 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1175 DATA (xqy(kl,23),kl=1,130)/0.019894, 0.019716, 0.019528, 0.019339, &
1176 0.019140, 0.018941, 0.018742, 0.018543, 0.018333, 0.018113, 0.017893, &
1177 0.017673, 0.017453, 0.017212, 0.016982, 0.016741, 0.016531, 0.016238, &
1178 0.015976, 0.015714, 0.015442, 0.015159, 0.014876, 0.014593, 0.014290, &
1179 0.013986, 0.013682, 0.013357, 0.013032, 0.012697, 0.012341, 0.011995, &
1180 0.011629, 0.011314, 0.010874, 0.010487, 0.010078, 0.009670, 0.009240, &
1181 0.008800, 0.008350, 0.007910, 0.007648, 0.007543, 0.007438, 0.007333, &
1182 0.007229, 0.007124, 0.007019, 0.006914, 0.006810, 0.006705, 0.006600, &
1183 0.006495, 0.006286, 0.005867, 0.005343, 0.004819, 0.004295, 0.003772, &
1184 0.003248, 0.002724, 0.002200, 0.001676, 0.001153, 0.000629, 0.000105, &
1185 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1186 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1187 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1188 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1189 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1190 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1191 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1192 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1193 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1194 !liqy
1195 !!CL2
1196 DATA (xs(kl,24),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1197 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1198 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1199 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1200 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1201 0.000E+00, 0.000E+00, 0.000E+00, 2.930E-21, 5.100E-21,7.270E-21,&
1202 1.212E-20, 1.870E-20, 2.564E-20, 3.932E-20, 5.408E-20,7.340E-20,&
1203 9.791E-20, 1.223E-19, 1.388E-19, 1.454E-19, 1.520E-19,1.586E-19,&
1204 1.652E-19, 1.718E-19, 1.784E-19, 1.850E-19, 1.902E-19,1.954E-19,&
1205 2.006E-19, 2.058E-19, 2.162E-19, 2.370E-19, 2.460E-19,2.550E-19,&
1206 2.450E-19, 2.350E-19, 2.115E-19, 1.880E-19, 1.600E-19,1.320E-19,&
1207 1.080E-19, 8.400E-20, 6.700E-20, 5.000E-20, 3.950E-20,2.900E-20,&
1208 2.350E-20, 1.800E-20, 1.550E-20, 1.300E-20, 1.130E-20,9.600E-21,&
1209 8.450E-21, 7.300E-21, 6.350E-21, 5.400E-21, 4.600E-21,3.800E-21,&
1210 3.200E-21, 2.600E-21, 2.100E-21, 1.600E-21, 0.000E+00,0.000E+00,&
1211 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1212 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1213 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1214 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1215 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1216 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1217 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00/
1219 DATA (xqy(kl,24),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
1220 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1221 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1222 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1223 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1224 0.000000, 0.000000, 0.000000, 1.000000, 1.000000, 1.000000, &
1225 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1226 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1227 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1228 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1229 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1230 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1231 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1232 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1233 1.000000, 1.000000, 1.000000, 1.000000, 0.000000, 0.000000, &
1234 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1235 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1236 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1237 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1238 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1239 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1240 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1242 !!HOCL
1243 DATA (xs(kl,25),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,&
1244 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 6.785E-20,6.071E-20,&
1245 5.600E-20, 5.402E-20, 5.489E-20, 5.914E-20, 6.645E-20,7.748E-20,&
1246 9.227E-20, 1.090E-19, 1.280E-19, 1.470E-19, 1.659E-19,1.828E-19,&
1247 1.960E-19, 2.031E-19, 2.058E-19, 2.018E-19, 1.924E-19,1.783E-19,&
1248 1.604E-19, 1.408E-19, 1.198E-19, 1.002E-19, 8.215E-20,6.768E-20,&
1249 5.650E-20, 4.958E-20, 4.650E-20, 4.671E-20, 4.934E-20,5.330E-20,&
1250 5.733E-20, 6.013E-20, 6.100E-20, 6.120E-20, 6.120E-20,6.120E-20,&
1251 6.095E-20, 6.070E-20, 6.020E-20, 5.970E-20, 5.905E-20,5.840E-20,&
1252 5.750E-20, 5.660E-20, 5.450E-20, 4.950E-20, 4.235E-20,3.500E-20,&
1253 2.810E-20, 2.220E-20, 1.765E-20, 1.430E-20, 1.205E-20,1.060E-20,&
1254 9.680E-21, 8.880E-21, 8.040E-21, 7.080E-21, 6.020E-21,4.910E-21,&
1255 3.845E-21, 2.880E-21, 2.080E-21, 1.440E-21, 9.700E-22,6.300E-22,&
1256 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1257 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1258 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1259 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1260 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1261 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1262 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1263 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1264 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00/
1266 DATA (xqy(kl,25),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
1267 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 1.000000, &
1268 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1269 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1270 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1271 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1272 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1273 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1274 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1275 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1276 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1277 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1278 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1279 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1280 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1281 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1282 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1283 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1284 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1285 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1286 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1287 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1289 !!CLNO2
1290 DATA (xs(kl,26),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,&
1291 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1292 3.105E-18, 3.173E-18, 3.242E-18, 3.330E-18, 3.418E-18,3.413E-18,&
1293 3.314E-18, 3.171E-18, 2.956E-18, 2.706E-18, 2.392E-18,2.090E-18,&
1294 1.814E-18, 1.591E-18, 1.385E-18, 1.236E-18, 1.100E-18,9.845E-19,&
1295 8.750E-19, 7.679E-19, 6.597E-19, 5.606E-19, 4.626E-19,3.846E-19,&
1296 3.159E-19, 2.629E-19, 2.268E-19, 1.976E-19, 1.782E-19,1.657E-19,&
1297 1.556E-19, 1.462E-19, 1.399E-19, 1.373E-19, 1.348E-19,1.319E-19,&
1298 1.290E-19, 1.261E-19, 1.232E-19, 1.203E-19, 1.171E-19,1.139E-19,&
1299 1.106E-19, 1.074E-19, 1.010E-19, 8.809E-20, 7.246E-20,5.831E-20,&
1300 4.600E-20, 3.572E-20, 2.734E-20, 2.077E-20, 1.564E-20,1.170E-20,&
1301 8.726E-21, 6.501E-21, 4.836E-21, 3.593E-21, 2.687E-21,2.008E-21,&
1302 1.498E-21, 1.125E-21, 8.471E-22, 6.419E-22, 4.858E-22,3.627E-22,&
1303 2.777E-22, 2.168E-22, 1.668E-22, 1.234E-22, 9.395E-23,7.482E-23,&
1304 6.059E-23, 4.981E-23, 3.794E-23, 3.169E-23, 0.000E+00,0.000E+00,&
1305 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1306 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1307 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1308 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1309 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1310 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1311 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00/
1313 DATA (xqy(kl,26),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
1314 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1315 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1316 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1317 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1318 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1319 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1320 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1321 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1322 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1323 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1324 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1325 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1326 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1327 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1328 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1329 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1330 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1331 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1332 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1333 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1334 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1336 !!FMCL
1337 DATA (xs(kl,27),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,&
1338 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1339 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1340 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1341 0.000E+00, 3.922E-20, 4.513E-20, 5.021E-20, 5.371E-20,5.541E-20,&
1342 5.452E-20, 5.817E-20, 5.800E-20, 5.645E-20, 5.236E-20,4.550E-20,&
1343 4.033E-20, 3.512E-20, 2.308E-20, 2.046E-20, 8.900E-21,8.214E-21,&
1344 3.511E-21, 2.221E-21, 1.498E-21, 1.181E-21, 8.634E-22,6.538E-22,&
1345 4.710E-22, 2.883E-22, 2.250E-22, 2.061E-22, 2.006E-22,1.790E-22,&
1346 1.557E-22, 1.323E-22, 9.346E-23, 0.000E+00, 0.000E+00,0.000E+00,&
1347 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1348 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1349 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1350 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1351 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1352 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1353 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1354 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1355 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1356 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1357 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00,&
1358 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,0.000E+00/
1360 DATA (xqy(kl,27),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
1361 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1362 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1363 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1364 0.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1365 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1366 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1367 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1368 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
1369 1.000000, 1.000000, 1.000000, 0.000000, 0.000000, 0.000000, &
1370 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1371 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1372 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1373 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1374 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1375 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1376 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1377 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1378 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1379 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1380 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1381 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1384 !liqy-20140911
1386 ! END MODULE module_data_photmad
1388 CONTAINS
1389 subroutine madronich1_driver(id,curr_secs,ktau,config_flags,haveaer,&
1390 gmt,julday,t_phy,moist,aerwrf,p8w,t8w,p_phy, &
1391 chem,rho_phy,dz8w, &
1392 xlat,xlong,z_at_w,gd_cloud,gd_cloud2, &
1393 ph_macr,ph_o31d,ph_o33p,ph_no2, &
1394 !liqy
1395 ph_cl2,ph_hocl,ph_clno2,ph_fmcl, &
1396 !liqy-20140911
1397 ph_no3o2,ph_no3o,ph_hno2, &
1398 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
1399 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
1400 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
1401 pm2_5_dry,pm2_5_water,uvrad, &
1402 ids,ide, jds,jde, kds,kde, &
1403 ims,ime, jms,jme, kms,kme, &
1404 its,ite, jts,jte, kts,kte )
1405 USE module_configure
1406 USE module_state_description
1407 USE module_model_constants
1408 USE module_data_radm2
1409 implicit none
1410 INTEGER, INTENT(IN ) :: id,julday, &
1411 ids,ide, jds,jde, kds,kde, &
1412 ims,ime, jms,jme, kms,kme, &
1413 its,ite, jts,jte, kts,kte
1414 INTEGER, INTENT(IN ) :: ktau
1415 REAL(KIND=8), INTENT(IN ) :: curr_secs
1416 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
1417 INTENT(IN ) :: moist
1418 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1419 INTENT(INOUT ) :: &
1420 pm2_5_dry,pm2_5_water
1421 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1422 OPTIONAL, &
1423 INTENT(INOUT ) :: &
1424 gd_cloud,gd_cloud2
1425 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1426 INTENT(INOUT ) :: &
1427 ph_macr,ph_o31d,ph_o33p,ph_no2, &
1428 !liqy
1429 ph_cl2,ph_hocl,ph_clno2,ph_fmcl, &
1430 !liqy-20140911
1431 ph_no3o2,ph_no3o,ph_hno2,&
1432 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
1433 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
1434 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob
1435 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
1436 INTENT(IN ) :: chem
1437 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1438 INTENT(IN ) :: &
1439 t_phy, &
1440 p_phy, &
1441 dz8w, &
1442 t8w,p8w,z_at_w , &
1443 aerwrf , &
1444 rho_phy
1445 REAL, DIMENSION( ims:ime , jms:jme ) , &
1446 INTENT(IN ) :: &
1447 xlat, &
1448 xlong
1449 REAL, DIMENSION( ims:ime , jms:jme ) , &
1450 INTENT(INOUT ) :: uvrad
1451 REAL, INTENT(IN ) :: gmt
1453 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1455 LOGICAL, INTENT(IN) :: haveaer
1458 ! LOCAL VAR
1460 INTEGER(KIND=8) :: ixhour
1461 INTEGER :: ki,i,j,k,n,iprt
1462 ! photolysis input
1464 real tt(kts:kte+1),o33(kts:kte+1),rhoa(kts:kte+1),aerext(kts:kte+1),qll(kts:kte+1), &
1465 phizz(kts:kte+1),phot1(nreakj-1,kts:kte+1)
1466 real(kind=8) :: xtime,xhour
1467 real :: xmin,gmtp,uvb_dd1,uvb_du1,uvb_dir1
1468 real :: zenith,zenita,azimuth,dobsi
1469 real :: bext340,bexth2o,ctr
1470 integer :: naerspec
1471 real :: zsurf
1473 ! print *,'gmt,julday in madronich1= ',gmt,julday
1474 xtime=curr_secs/60._8 !min since simulation start
1475 ixhour=int(gmt+.01,8)+int(xtime/60._8,8)
1476 xhour=real(ixhour,8)
1477 xmin=60.*gmt+real(xtime-xhour*60._8,8)
1478 gmtp=mod(xhour,24._8)
1479 gmtp=gmtp+xmin/60.
1480 ! print *,'gmtp = ',gmtp,xhour,xmin
1481 do 100 j=jts,jte
1482 do 100 i=its,ite
1483 ! write(0,*)i,j
1484 do k=kts,kte+1
1485 do n=1,nreakj-1
1486 phot1(n,k)=0.
1487 END DO
1488 END DO
1489 iprt = 0
1490 zenith=0.
1491 zenita=0.
1492 azimuth=0.
1493 call calc_zenith(xlat(i,j),-xlong(i,j),julday,gmtp,azimuth,zenith)
1494 ! if nighttime, skip radiative transfer calculation
1495 if(zenith.eq.90.) zenith = 89.9
1496 if(zenith.ge.90.) go to 199
1497 zenita = cos(zenith*pi/180.)
1498 if(zenith.gt.75.) zenita=1./chap(zenith)
1499 zsurf = z_at_w(i,1,j)*0.001
1501 ! photmad berechnet photolysefrequenzen nach Madronich in folgender Reihenfolge
1503 ! o2 absorp 1 schumman-runge corrected in srband
1504 ! o3 -> 1d 2 at 275 k. correct t-dep in subgrid
1505 ! o3 -> 3p 3 at 275 k. correct t-dep in subgrid
1506 ! no2 4
1507 ! no3 -> no+o2 5
1508 ! no3 -> no2+o 6
1509 ! hno2 7
1510 ! hno3 8
1511 ! hno4 9
1512 ! h2o2 10
1513 ! ch2o -> rad 11
1514 ! ch2o -> mol 12
1515 ! ch3cho 13
1516 ! ch3coch3 14
1517 ! ch3coc2h5 15
1518 ! hcocho proc a 16
1519 ! ch3cocho 17
1520 ! hcoch=chcho 18 estimate, no reliable measurement
1521 ! ch3o2h 19
1522 ! ch3coo2h 20 actually use 0.28*(h2o2 value)
1523 ! ch3ono2 21
1524 ! hcocho proc b 22
1525 do k=kts,kte
1526 aerext(k+1)=aerwrf(i,k+1,j)
1527 !
1528 !--- if you have aerosols
1529 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1530 bext340=5.E-6
1531 bexth2o=5.E-6
1532 if(haveaer.and.ktau.gt.1)then
1534 ! dry aerosol mass
1535 !rf check ki (or ki+1 or ki-1 ?)
1536 aerext(k)=pm2_5_dry(i,k,j)*bext340+ &
1537 pm2_5_water(i,k,j)*bexth2o
1538 aerext(k)=aerext(k)*1.E3
1539 ! if(i.eq.70.and.j.eq.70) write(06,*) 'aerext',k,aerext(k)
1541 endif
1542 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1543 qll(k+1)=0.
1544 tt(k+1) = t_phy(i,k,j)
1545 rhoa(k+1) = rho_phy(i,k,j)
1546 o33(k+1) = max(1.e-3,chem(i,k,j,p_o3))
1548 IF (config_flags%cu_diag==1) THEN
1549 qll(k+1) = 1.e3*(moist(i,k,j,p_qc)+moist(i,k,j,p_qi)+ &
1550 gd_cloud(i,k,j)+gd_cloud2(i,k,j)) &
1551 *rho_phy(i,k,j)
1552 else
1553 qll(k+1) = 1.e3*(moist(i,k,j,p_qc)+moist(i,k,j,p_qi)) &
1554 *rho_phy(i,k,j)
1555 ENDIF
1557 if(qll(k+1).lt.1.e-5)qll(k+1) = 0.
1558 phizz(k+1) = z_at_w(i,k+1,j)*.001-z_at_w(i,1,j)*.001
1559 ! if((i.eq.1.and.j.eq.17))then
1560 ! write(0,*)k+1,phizz(k+1),qll(k+1),moist(i,k,j,p_qc),moist(i,k,j,p_qi),rqccuten(i,k,j),rqicuten(i,k,j),rho_phy(i,k,j)
1561 ! write(0,*)k+1,z_at_w(i,k+1,j),z_at_w(i,1,j),tt(k+1),o33(k+1),qll(k+1),rhoa(k+1)
1562 ! endif
1563 END DO
1564 tt(1)=t8w(i,kts,j)
1565 o33(1)=max(1.e-3,chem(i,kts,j,p_o3))
1566 qll(1)=0.
1567 phizz(1)=0.
1568 aerext(1)=aerwrf(i,1,j)
1569 k=0
1570 ! write(0,*)k+1,z_at_w(i,k+1,j),z_at_w(i,1,j),tt(k+1),o33(k+1),qll(k+1),rhoa(k+1)
1571 !
1572 ! if you have aerosols....
1573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1574 !
1575 if(haveaer.and.ktau.gt.1)then
1576 aerext(1)=pm2_5_dry(i,1,j)*bext340+ &
1577 pm2_5_water(i,1,j)*bexth2o
1578 aerext(1)=aerext(1)*1.e3
1579 endif
1580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1581 rhoa(1)=p8w(i,kts,j)/t8w(i,kts,j)/r_d
1582 ! dobsi=350.
1583 dobsi=325.
1584 ! if((i.eq.87.and.j.eq.66).or.(i.eq.105.and.j.eq.70))then
1585 ! print *,'before photmad, i,j = ',i,j,pm2_5_dry(i,1,j),pm2_5_water(i,1,j)
1586 ! print *,k,rhoa(1),phizz(1),qll(1),aerext(1),o33(1),tt(1),zenita
1587 ! endif
1588 !write(0,*)'calling photolysis_mad --------',i,j
1589 call photolysis_mad(kte,zenita,zenith,zsurf, &
1590 phizz,tt,rhoa,o33,aerext,qll,dobsi,phot1, &
1591 nreacj,iprt,uvb_dd1,uvb_du1,uvb_dir1)
1592 !write(0,*)'back from photolysis_mad ---------------------------- '
1593 ! print *,'after photmad, i,j = ',i,j
1594 uvrad(i,j)=uvb_dd1+uvb_dir1-uvb_du1
1595 do k=kts,kte
1596 do n=1,nreakj-1
1597 phot1(n,k)=60.*phot1(n,k)
1598 END DO
1599 END DO
1600 199 continue
1601 do k=kts,kte
1602 !
1603 !
1604 ph_o31d(i,k,j) = phot1(1,k)
1605 ph_o33p(i,k,j) = phot1(2,k)
1606 ph_no2(i,k,j) = phot1(3,k)
1607 ph_no3o2(i,k,j) = phot1(4,k)
1608 ph_no3o(i,k,j) = phot1(5,k)
1609 ph_hno2(i,k,j) = phot1(6,k)
1610 ph_hno3(i,k,j) = phot1(7,k)
1611 ph_hno4(i,k,j) = phot1(8,k)
1612 ph_h2o2(i,k,j) = phot1(9,k)
1613 ph_ch2or(i,k,j) = phot1(10,k)
1614 ph_ch2om(i,k,j) = phot1(11,k)
1615 ph_ch3cho(i,k,j) = phot1(12,k)
1616 ph_ch3coch3(i,k,j) = phot1(13,k)
1617 ph_ch3coc2h5(i,k,j) = phot1(14,k)
1618 ph_hcocho(i,k,j) = phot1(15,k)
1619 ph_ch3cocho(i,k,j) = phot1(16,k)
1620 ph_hcochest(i,k,j) = phot1(17,k)
1621 ph_ch3o2h(i,k,j) = phot1(18,k)
1622 ph_ch3coo2h(i,k,j) = phot1(19,k)
1623 ph_ch3ono2(i,k,j) = phot1(20,k)
1624 ph_hcochob(i,k,j) = phot1(21,k)
1625 ph_macr(i,k,j) = phot1(22,k)
1626 !liqy
1627 ph_cl2(i,k,j) = phot1(23,k)
1628 ph_hocl(i,k,j) = phot1(24,k)
1629 ph_clno2(i,k,j) = phot1(25,k)
1630 ph_fmcl(i,k,j) = phot1(26,k)
1632 !liqy-20140911
1633 ! if(i.eq.5.and.j.eq.5)print *,i,j,k,phot1(3,k),phot1(4,k),phot1(17,k)
1634 END DO
1636 100 continue
1638 END SUBROUTINE madronich1_driver
1640 SUBROUTINE photolysis_mad(mkxcc,a,zenith,zsurf, &
1641 zmm5,tmm5,pmm5,o3mm5,aerext,wlmm5, &
1642 dobsnew,phot1,nrtest,iprt,uvb_dd1,uvb_du1,uvb_dir1)
1643 ! Note that nreakj can differ from nreacj in csolve1!!!!!!
1644 ! ----- INPUT ------------------------------------------------
1645 ! Input from MM5
1646 ! ------ OUTPUT ----------------------------------------------
1647 ! -----------------------------------------------------------------------
1648 ! This program calculates J-values for selected atmospheric molecules
1649 ! the program and subroutine structure is as follows
1650 ! runph -main program and menu.
1651 ! calc_zenith -zenith angle calc.
1652 ! in addition, the chapman function ch(zenith) is used between
1653 ! readd -reads spectra and standard profiles
1654 ! photmat -main subroutine (calls the routines below)
1655 ! o3scal -rescales ozone profile to user-selected new dobson
1656 ! subgrid -regrids the altitude profiles (
1657 ! srband -computes effective ozone cross sections in the
1658 ! schumamkxcc+1+nabv-runge region (if necessary)
1659 ! trapez -interpolation subroutine
1660 ! optics -computes optical parameters and calls delta-eddington
1661 ! delted -delta-eddington code of wiscombe
1662 ! nlayde(ib) -multylayer calc
1663 ! leqt1b(a,n,ncl,nuc,ia,b,m,ib,ijob,kl) -solves
1664 ! pentadiag.system
1665 ! Original program written by S. Madronich, NCAR. The code was
1666 ! last modified by him on 18 Aug. 1987.
1667 ! The driver for the chemistry has been heavily modified by
1668 ! W. R. Stockwell at IFU Germany. (The modifications allow
1669 ! the simulation conditions to be directly modified
1670 ! without recompiling the program [this does not hold for the present
1671 ! version here any more]). The extremely large data base
1672 ! has bee replaced with separate files to allow cross sections and
1673 ! photolysis rates to be more easily updated. However, the order
1674 ! of the input files should be modified with extreme caution.
1675 ! The present version is made for the use of computed
1676 ! fields of T, P, O3, wl, and aerosol (if available)
1677 ! Modifications were made by Renate Forkel in Aug./Sept 1995
1678 ! The following modifications were made:
1679 ! - No climatological input any more
1680 ! - Arbitrary vertical grid
1681 ! - Several seperate cloud layers are possible. Number of additional
1682 ! layers within the clouds depends on the optical depth
1683 ! - Vertically inhogeneous clouds are possible now
1684 ! - No variable values in commons any more
1685 ! -----------------------------------------------------------------------
1686 ! xnkg: Molecules per kg air
1687 ! .. Scalar Arguments ..
1688 REAL :: a, dobsnew,uvb_dd1,uvb_du1,uvb_dir1
1689 INTEGER :: iprt, mkxcc, nrtest
1690 ! .. Local Scalars ..
1691 REAL :: aeru, airu, bextn, df, dj, dz, ff, gaer, gcld, gray, haer, &
1692 hair, ho3, o3u, omaer, omcld, omray, reff, tu, wmicron, xnkg, xx, &
1693 znorm, zu
1694 INTEGER :: i, j, k, kk, kl, lev, nlayer, nlevel, nn, nr, nsurf
1695 ! .. Array Arguments ..
1696 REAL :: aerext(mkxcc+1), o3mm5(mkxcc+1), phot1(nreakj-1,mkxcc+1), &
1697 pmm5(mkxcc+1), tmm5(mkxcc+1), wlmm5(mkxcc+1), zmm5(mkxcc+1)
1698 ! .. Local Arrays ..
1699 REAL :: aaer(130), aer(mkxcc+1+nabv), air(mkxcc+1+nabv), ao2(nj,130), &
1700 ao3(nj,130), arayl(130), cloud(mkxcc+1+nabv), cvo2(nj), &
1701 d(nreakj,nj), endir(nj,130), endn(nj,130), enup(nj,130), &
1702 hilf1(mkxcc+1), hilfd(nj), o3(mkxcc+1+nabv), qy(nj,130,nreakj), &
1703 s(nj,130,nreakj), t(mkxcc+1+nabv), vaer(nj), vair(nj), vcld(nj), &
1704 vo3(nj), vt(nj), z(nj), zkm(mkxcc+1), zmid(nj), zz(mkxcc+1+nabv)
1706 REAL :: fldir(nj,130), fldn(nj,130), flup(nj,130)
1708 ! .. Data Statements ..
1709 DATA xnkg/2.143E25/
1711 REAL :: zenith, zsurf
1712 REAL zasl(nj)
1714 ! EXTERNAL sphers, airmas
1715 INTEGER, PARAMETER :: kzmax = 200
1716 REAL :: dsdh(0:kzmax,kzmax)
1717 INTEGER :: nid(0:kzmax)
1718 INTEGER :: ii
1719 REAL vcol(nj), scol(nj)
1721 ! EXTERNAL o3scal, optics, srband, subgrid, trapez
1722 ! wave length range !!! If photolysis is also desired for levels above
1723 ! 2
1724 ! kl0 should be set equal to 1 again!!!!!!!!!!!!!!
1725 i = 1
1726 j = 1
1727 ! albedoph ************************ specify ground albedoph
1728 ! use best estimate albedoph of demerjian et al.,
1729 ! adv.env.sci.tech.,v.10,p.369, (1980)
1730 DO kl = kl0, kl1
1732 IF (wl(kl)<400.) albedoph(kl) = 0.05
1734 IF ((wl(kl)>=400.) .AND. (wl(kl)<450.)) albedoph(kl) = 0.06
1736 IF ((wl(kl)>=450.) .AND. (wl(kl)<500.)) albedoph(kl) = 0.08
1738 IF ((wl(kl)>=500.) .AND. (wl(kl)<550.)) albedoph(kl) = 0.10
1740 IF ((wl(kl)>=550.) .AND. (wl(kl)<600.)) albedoph(kl) = 0.11
1742 IF ((wl(kl)>=600.) .AND. (wl(kl)<640.)) albedoph(kl) = 0.12
1744 IF ((wl(kl)>=640.) .AND. (wl(kl)<660.)) albedoph(kl) = 0.135
1746 IF (wl(kl)>=660.) albedoph(kl) = 0.15
1747 END DO
1751 nn = mkxcc + 1 + nabv
1753 ! omray = single scattering albedoph, rayleigh. use 1.00
1754 ! gray = asymetry factor for rayleigh scattering. use 0.0
1755 ! arayl(kl) = rayleigh scattering cross section, from
1756 ! frohlich and shaw, appl.opt. v.11, p.1773 (1980).
1757 ! overrides tabulation of jdata.base
1758 hair = 8.05
1759 omray = 1.0
1760 gray = 0.0
1762 DO 10 kl = kl0, kl1
1763 wmicron = wl(kl)/1.E3
1764 xx = 3.916 + 0.074*wmicron + 0.050/wmicron
1765 arayl(kl) = 3.90E-28/(wmicron)**xx
1766 10 CONTINUE
1768 ! aerosol*********************** specify aerosols
1769 ! aaer(kl) = aerosol total vertical optical depth variation with
1770 ! wavelength. estimated from elterman (1968)
1771 ! aer(i) = attenuation (per km) profile from elterman (1968).
1772 ! given in data statement in beginning of code, for 340 nm (kl=6
1773 ! same vertical shape at all wavelengths.
1774 ! normalized later (in subroutine subgrid) to total vertical dep
1775 ! this wavelength.
1776 ! omaer = aerosol single scattering albedoph. use 0.99 for now.
1777 ! gaer = aerosol asymetry factor. use 0.61 (hansen and travis 197
1778 ! (these are assuming particles of about 0.1 micron radius
1779 ! index of refraction of about 1.65 + 0.002i.
1780 ! haer = the aerosol scale height at top of atmosphere
1781 ! use equal to air (8.05 km)
1783 DO 20 kl = kl0, kl1
1784 aaer(kl) = 0.379*(340./wl(kl))
1785 20 CONTINUE
1786 omaer = 0.990
1787 gaer = 0.610
1788 haer = 8.05
1791 omcld = 1.000
1792 gcld = 0.860
1794 ho3 = 4.50
1796 nsurf = 1
1798 ! Vertikales Gitter ohne Wolken: Niveaus zz(i), i=1,mkxcc+1+nabv
1801 ! Transformation of MM5-values
1802 ! bextn: cloud extinction coeff per g/m**3 [1/km]
1803 DO k = 1, mkxcc + 1
1804 zz(k) = zmm5(k)
1805 !write(0,*)' here7a zz = zmm5 ',k,zmm5(k)
1806 t(k) = tmm5(k)
1807 air(k) = xnkg*pmm5(k)*1.E-6
1809 ! falls o3mm5 in ppm
1811 o3(k) = o3mm5(k)*1.E-6*air(k)
1812 aer(k) = aerext(k)
1814 ! bextn: cloud extinction coeff per g/m**3 [1/km]
1815 ! **** warning: parameterization for bextn only good
1816 ! for continental clouds
1818 IF (wlmm5(k)>0.) THEN
1820 ! vereinfachte Version, falls kein
1821 ! Sulfat uebergeben wird.
1823 reff = 9.6*wlmm5(k)**0.333
1824 bextn = (0.0275+1.3/reff)*1000.
1825 cloud(k) = wlmm5(k)*bextn
1826 ELSE
1827 cloud(k) = 0.
1828 END IF
1829 ! if(iprt.eq.1)write(6,'(i3,e12.4)')k,cloud(k)
1830 END DO
1832 znorm = (50.-zz(mkxcc+1))/(50.-20.)
1834 DO k = 1, nabv
1835 zz(mkxcc+1+k) = 50. - znorm*(50.-zabv(k))
1836 !write(0,*)' here8a ',k,' zz(',mkxcc+1+k,') ',zz(mkxcc+1+k),znorm,zabv(k)
1837 END DO
1839 zu = zz(mkxcc+1)
1840 tu = t(mkxcc+1)
1841 o3u = o3(mkxcc+1)
1842 airu = air(mkxcc+1)
1843 aeru = aer(mkxcc+1)
1844 kk = 1
1846 ! Zufuegen von Werten oberhalb von MM5
1848 ! Die 'abv'-Werte sind bereits in den richtigen Einheiten
1849 DO k = mkxcc + 1 + 1, mkxcc + 1 + nabv
1850 ! write(0,'(2i3,5e12.3)') k,kk,zz(k),zabv(kk)
1851 30 IF (zz(k)<=zabv(kk)) THEN
1852 dz = zz(k) - zu
1853 ff = dz/(zabv(kk)-zu)
1854 t(k) = tu + ff*(tabv(kk)-tu)
1855 o3(k) = o3u + ff*(o3abv(kk)-o3u)
1856 air(k) = airu + ff*(pabv(kk)-airu)
1857 aer(k) = aeru + ff*(caabv(kk)-aeru)
1858 cloud(k) = 0.
1859 ! if(iprt.eq.1)then
1860 ! write(0,'(2i3,5e12.3)') k,kk,zz(k),zabv(kk),ff,tabv(kk),air(k)
1861 ! endif
1862 ELSE
1863 40 zu = zabv(kk)
1864 tu = tabv(kk)
1865 o3u = o3abv(kk)
1866 airu = pabv(kk)
1867 aeru = caabv(kk)
1868 kk = kk + 1
1870 IF (zabv(kk)<zz(k)) GO TO 40
1871 GO TO 30
1872 END IF
1874 END DO
1876 IF (iprt==1) THEN
1877 !WRITE (06,'('k, zz(k), t(k), air(k), o3(k), aer(k) cloud(k)')')
1879 !DO k = 1, mkxcc + 1 + nabv
1880 ! WRITE (06,'(i3,6e12.4)') k, zz(k), t(k), air(k), o3(k), aer(k), &
1881 ! cloud(k)
1882 !END DO
1884 END IF
1885 ! if(iprt.eq.1)then
1886 ! do k=1,mkxcc+1
1887 ! write(06,'(i3,6e12.4)') k,zz(k),t(k),air(k),o3mm5(k), &
1888 ! aerext(k),cloud(k)
1889 ! enddo
1890 ! endif
1891 ! Scaling of O3-column with specified Dobson units
1892 ! -----------------------------------------------------------------------
1893 CALL o3scal(dobsnew,ho3,zz,o3,mkxcc+1+nabv)
1894 ! stop
1895 ! -----------------------------------------------------------------------
1897 ! Einfuegen zusaetzlicher Schichten in Wolken,
1898 ! T-abh. Absorptionsquerschnitte, etc
1900 ! -----------------------------------------------------------------------
1901 nlevel = 0
1903 IF (iprt==1) PRINT *, 'nlevel = ', nlevel
1905 CALL subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
1906 vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,mkxcc+1,mkxcc+1+nabv)
1908 IF (iprt==1) PRINT *, 'nlevel = ', nlevel
1909 ! -----------------------------------------------------------------------
1911 ! compute altitudes relative to sea level
1912 DO ii = 1, nlevel
1913 zasl(ii) = z(ii) + zsurf
1914 ENDDO
1915 ! compute spherical increments
1916 CALL sphers(kzmax, nlevel, z, zenith, dsdh, nid)
1917 ! here can also call airmass to compute air mass factors for Lyman alpha and
1918 ! Schmumann-Runge parameterizations
1919 ! CALL airmas(nlevel, dsdh, nid, vair, vcol, scol)
1921 ! subroutine srband computes effective ozone cross sections in the
1922 ! schumann-runge region, using the parameterization of allen and
1923 ! freder
1925 ! ----------------------------------------------------------------------
1926 CALL srband(a,zmid,cvo2,ao2,vt,nlevel)
1927 ! ----------------------------------------------------------------------
1930 ! subroutine optics is the driver for the flux calculation. output is
1931 ! endir(lev,kl) -irradiance of direct solar beam
1932 ! endn(lev,kl) -irradiance of down-welling diffuse light
1933 ! enup(lev,kl) -irradiance of up-welling diffuse light
1934 ! for each level lev and wavelength bin kl. conversion to actinic
1935 ! flux is in loop 90
1937 ! ----------------------------------------------------------------------
1939 CALL optics(iprt, &
1940 zenith, a, dsdh, nid, &
1941 vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld,omcld, &
1942 vaer,aaer,gaer,omaer,nlevel, kzmax, &
1943 endir,endn,enup, fldir,fldn,flup)
1945 ! ----------------------------------------------------------------------
1947 DO nr = 1, nreakj
1949 DO k = 1, nlevel
1950 d(nr,k) = 0.
1951 END DO
1953 END DO
1955 DO 50 kl = kl0, kl1
1957 DO 50 lev = 1, nlevel
1959 ! calculate actinic flux, direct + diffuse-dn + diffuse-up
1960 df = fext(kl)*(fldir(lev,kl) + fldn(lev,kl) + flup(lev,kl))
1962 ! calculate contribution to j-vals, sum over wavelengths
1964 DO 50 nr = 1, nreakj
1965 dj = df*s(lev,kl,nr)*qy(lev,kl,nr)
1966 d(nr,lev) = d(nr,lev) + dj
1967 50 CONTINUE
1969 ! calculate spectral irradiance at surface:
1970 uvb_dd1=0.
1971 uvb_du1=0.
1972 uvb_dir1=0.
1973 do kl = 38,54
1974 ! radfakt=fext(kl) ! photons/cm**2/s/bin
1975 ! ! UV in Photonen/cm**2/s
1976 radfakt=fext(kl)*1e4 * 6.63e-34 * 2.9979e8 / (wl(kl)*1e-9)
1977 ! UV in W/m**2
1978 uvb_dd1=uvb_dd1+radfakt*endn(1,kl)
1979 uvb_du1=uvb_du1+radfakt*enup(1,kl)
1980 uvb_dir1=uvb_dir1+radfakt*endir(1,kl)
1981 enddo
1983 ! print *,'d(4,1) = ',d(4,1)
1985 ! interpolate j-values to new altitudes:
1987 DO nr = 1, nreakj - 1
1988 ! write(06,'(' nr',i3)') nr
1989 DO lev = 1, mkxcc + 1
1990 zkm(lev) = zmm5(lev)
1991 END DO
1993 DO lev = 1, nlevel
1994 hilfd(lev) = d(nr+1,lev)
1995 ! if(iprt.eq.1.and.nr.eq.3)
1996 ! 1 write(06,'(i3,2e12.4)') lev,z(lev),d(4,lev)
1997 END DO
1998 ! hilfd muss noch auf mm5-hoehen zurueckinterpoliert werden
1999 CALL trapez(z,hilfd,1,nlevel,zkm,hilf1,1,mkxcc+1,nj,nj,mkxcc+1, &
2000 mkxcc+1)
2002 IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '
2004 IF (iprt==1 .AND. nr==3) WRITE (06,'(5e12.6)') (hilf1(k),k=1, &
2005 mkxcc+1)
2007 IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '
2009 DO lev = 1, mkxcc + 1
2010 phot1(nr,lev) = hilf1(lev)
2011 END DO
2013 END DO
2015 RETURN
2017 END SUBROUTINE photolysis_mad
2019 ! #################################################################
2021 SUBROUTINE srband(a,zmid,cvo2,ao2,vt,nlevel)
2022 ! this subroutine calculates effective o2 cross sections in the
2023 ! schumann-runge band, using the formulae from allen and frederick,
2024 ! j.atmos.sci., v39, p2066 (1982). the coefficients are those modifi
2025 ! wuebbles (lll report, 1982).
2026 ! .. Scalar Arguments ..
2027 REAL :: a
2028 INTEGER :: nlevel
2029 ! .. Local Scalars ..
2030 REAL :: ao20, ao20lg, c, clog, e10, x1, x2, x3, x4, x5, x6, x7, x8, &
2031 y1, y2, y3, y4, zendep
2032 INTEGER :: lay, n20, nlayer
2033 ! .. Intrinsic Functions ..
2034 INTRINSIC alog
2035 ! .. Array Arguments ..
2036 REAL :: ao2(nj,130), cvo2(nj), vt(nj), zmid(nj)
2038 nlayer = nlevel - 1
2039 ! initialize cross sections:
2040 DO 10 lay = 1, nlayer
2042 DO 10 kl = kl0, kl1
2043 10 ao2(lay,kl) = xs(kl,1)
2046 ! correct as needed
2047 ! use formula for 20-50 km.
2048 ! below 20 km, use 20 km value
2049 ! find layer near 20 km
2050 DO 20 lay = 1, nlayer
2052 IF (zmid(lay)>20.) THEN
2053 n20 = lay
2054 GO TO 30
2055 END IF
2057 20 CONTINUE
2058 30 CONTINUE
2060 e10 = alog(10.)
2062 DO 60 kl = kl0, kl1
2064 IF (wl(kl)>205.) RETURN
2066 DO 40 lay = n20, nlayer
2067 x1 = alog(4.696E-23*cvo2(lay)/0.2095)/e10
2069 IF (wl(kl)>=200.) x1 = vt(lay)
2070 x2 = x1*x1
2071 x3 = x2*x1
2072 x4 = x3*x1
2073 x5 = x4*x1
2074 x6 = x5*x1
2075 x7 = x6*x1
2076 x8 = x7*x1
2077 ao20lg = sra(kl,1) + sra(kl,2)*x1 + sra(kl,3)*x2 + sra(kl,4)*x3 + &
2078 sra(kl,5)*x4 + sra(kl,6)*x5 + sra(kl,7)*x6 + sra(kl,8)*x7 + &
2079 sra(kl,9)*x8
2080 ao20 = 10.**ao20lg
2082 y1 = alog(cvo2(lay))/e10
2083 y2 = y1*y1
2084 y3 = y2*y1
2085 y4 = y3*y1
2086 clog = srb(kl,1) + srb(kl,2)*y1 + srb(kl,3)*y2 + srb(kl,4)*y3 + &
2087 srb(kl,5)*y4
2088 c = 10.**clog
2089 zendep = a**c
2091 ao2(lay,kl) = ao20*zendep
2092 40 CONTINUE
2093 ! assign values below 20 km
2094 DO 50 lay = 1, n20 - 1
2095 ao2(lay,kl) = ao2(n20,kl)
2096 50 CONTINUE
2097 60 CONTINUE
2099 RETURN
2101 END SUBROUTINE srband
2103 ! ######################################################################
2105 SUBROUTINE o3scal(dobsnew,ho3,zz,o3,nn)
2106 ! adjustment of o3 profiles to a user-selected dobson value.
2107 ! select value of dobnew in main program
2108 ! if don't want to use, don't call this subroutine
2109 ! .. Scalar Arguments ..
2110 REAL :: dobsnew, ho3
2111 INTEGER :: nn
2112 ! .. Local Scalars ..
2113 REAL :: dobsref
2114 INTEGER :: i
2115 ! .. Intrinsic Functions ..
2116 INTRINSIC max, min
2117 ! .. Array Arguments ..
2118 REAL :: o3(nn), zz(nn)
2119 ! write(6,*) o3
2120 dobsref = o3(nn)*1.E5*ho3
2121 ! write(06,'('nn: dobsref,dobsnew',2e12.4)')
2122 ! & dobsref/2.687e16,dobsnew
2123 DO 10 i = 1, nn
2124 10 dobsref = dobsref + o3(i)*0.5*(zz(min(i+1,nn))-zz(max(i-1,1)))*1.E5
2125 dobsref = dobsref/2.687E16
2126 ! write(06,'('dobsref,dobsnew',2e12.4)') dobsref,dobsnew
2127 DO 20 i = 1, nn
2128 o3(i) = o3(i)*dobsnew/dobsref
2129 20 CONTINUE
2130 ! write(06,*) o3
2132 RETURN
2134 END SUBROUTINE o3scal
2136 ! #######################################################################
2138 FUNCTION chap(zeta)
2139 ! chapman function is used when the solar zenith angle exceeds
2140 ! 75 deg.
2141 ! interpolates between values given in, e.g., mccartney (1976).
2142 ! .. Scalar Arguments ..
2143 REAL :: zeta
2144 ! .. Local Scalars ..
2145 REAL :: rm
2146 INTEGER :: i
2147 ! .. Local Arrays ..
2148 REAL :: y(22)
2149 ! .. Function Return Value ..
2150 REAL :: chap
2151 ! .. Data Statements ..
2152 DATA (y(i),i=1,22)/3.800, 4.055, 4.348, 4.687, 5.083, 5.551, 6.113, &
2153 6.799, 7.650, 8.732, 10.144, 12.051, 14.730, 18.686, 24.905, 35.466, &
2154 55.211, 96.753, 197., 485., 1476., 9999./
2156 DO 10 i = 75, 96
2157 rm = i
2159 IF (zeta<rm) GO TO 20
2160 10 CONTINUE
2161 20 chap = y(i-75) + (y(i-74)-y(i-75))*(zeta-(rm-1.))
2162 END FUNCTION chap
2165 ! #################################################################
2166 SUBROUTINE subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
2167 vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,nz,nn)
2168 ! Output
2169 ! this subroutine computes all altitude dependent quantities over the
2170 ! sub-divided grid
2171 ! the following quantities are computed at each level (altitude)
2172 ! z(lev) = altitude (km) of each level
2173 ! zair(lev) = air concentration at each level
2174 ! zt(lev) = temperature at each level
2175 ! s(lev,kl,nr) = abs. cross sect at each altitude (t,p corrected)
2176 ! qy(lev,kl,nr) = quantum yield at each altitude (t,p corrected)
2177 ! the following quantities are computed at each layer (thickness)
2178 ! zmid(lay) = altitude of midpoint of layer
2179 ! vair(lay) = air column in layer (vertical)
2180 ! vo3(lay) = ozone column ' '
2181 ! vaer(lay) = aerosol ' ' '
2182 ! vcld(lay) = cloud ' ' '
2183 ! vt(lay) = average temperature of column
2184 ! ao3(lay,kl) = average o3 cross sect in layer (with ave. layer te
2185 ! .. Scalar Arguments ..
2186 REAL :: haer, hair, omray
2187 INTEGER :: nlayer, nlevel, nn, nz
2188 ! .. Local Scalars ..
2189 REAL :: a, ak300, akt, b, c, dz, dzo, dzt, dzu, fnum, hlocal, phi1, &
2190 phi2, phi20, sum, tau, tdiffx, x0, xl, xl0
2191 REAL :: kt, q1, q2 ! SAMcKeen add for updated JPL recommended JO1D 2/14/13
2192 INTEGER :: i, idt, ii, il, lay, lev, llbt, nr
2193 ! .. Intrinsic Functions ..
2194 ! EXTERNAL trapez
2195 INTRINSIC alog, atan, exp, float, ifix, max
2196 ! .. Array Arguments ..
2197 REAL :: aer(nn), air(nn), ao3(nj,130), cloud(nn), cvo2(nj), o3(nn), &
2198 qy(nj,130,nreakj), s(nj,130,nreakj), t(nn), vaer(nj), vair(nj), &
2199 vcld(nj), vo3(nj), vt(nj), z(nj), zmid(nj), zz(nn)
2200 ! .. Local Arrays ..
2201 REAL :: zair(nj), zt(nj)
2202 ! ------------------------------------------------------ levels
2203 ! fnum: factor for calculation of additional levels in clouds
2204 fnum = 0.02
2206 DO i = 1, nj
2207 vcld(i) = 0.
2208 END DO
2210 lev = 0
2211 llbt = 0
2213 DO 30 i = 1, nn
2214 ! write(06,'('cloud',i3,2e12.4)') i,zz(i),cloud(i)
2215 IF (cloud(i)<=0.) THEN
2216 lev = lev + 1
2217 z(lev) = zz(i)
2218 zair(lev) = air(i)
2219 zt(lev) = t(i)
2220 ! write(06,'('z,t,air ',i3,3e12.4)')
2221 ! lev,z(lev),zt(lev),zair(lev)
2222 ELSE
2224 IF (i==1 .AND. lev==0) THEN
2225 llbt = llbt + 1
2226 ! lbase(llbt)=1
2227 lev = 1
2228 z(lev) = zz(i)
2229 zair(lev) = air(i)
2230 zt(lev) = t(i)
2231 END IF
2233 IF (i==1) GO TO 10
2235 IF (cloud(i-1)<=0.) THEN
2236 lev = lev + 1
2237 llbt = llbt + 1
2238 ! lbase(llbt)=lev
2239 z(lev) = 0.5*(zz(i)+zz(i-1))
2240 zt(lev) = 0.5*(t(i)+t(i-1))
2241 if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
2242 zair(lev) = air(i-1)
2243 else
2244 hlocal = 1./alog(air(i-1)/air(i))
2245 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
2246 zair(lev) = air(i-1)*exp(-x0/hlocal)
2247 endif
2249 ! write(06,'('u:z,t,air ',i3,3e12.4)')
2250 ! lev,z(lev),zt(lev),zair(lev)
2251 END IF
2253 dzt = zz(i) - z(lev)
2254 ! number of layers depents on optical depth
2255 idt = max(ifix(cloud(i)*dzt*fnum),1)
2256 dzu = dzt/float(idt)
2257 vcld(lev) = cloud(i)
2259 DO il = 1, idt
2260 lev = lev + 1
2262 IF (lev>nj) THEN
2263 CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
2264 END IF
2266 z(lev) = z(lev-1) + dzu
2267 zt(lev) = t(i-1) + (z(lev)-zz(i-1))/(zz(i)-zz(i-1))*(t(i)-t(i-1) &
2268 )
2270 if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
2271 zair(lev) = air(i-1)
2272 else
2273 hlocal = 1./alog(air(i-1)/air(i))
2274 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
2275 zair(lev) = air(i-1)*exp(-x0/hlocal)
2276 endif
2277 vcld(lev) = cloud(i)
2278 ! write(06,'('u:z,t,air ',i3,3e12.4)')
2279 ! lev,z(lev),zt(lev),zair(lev)
2280 END DO
2282 10 CONTINUE
2284 IF (i==nn) GO TO 20
2285 ! if(lev.ne.1) vcld(lev)=cloud(i)
2286 dzt = (zz(i+1)-zz(i))*.5
2287 ! number of layers depents on optical depth
2288 idt = max(ifix(cloud(i)*dzt*fnum),1)
2289 dzo = dzt/float(idt)
2291 DO il = 1, idt
2292 lev = lev + 1
2294 IF (lev>nj) THEN
2295 CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
2296 END IF
2298 z(lev) = z(lev-1) + dzo
2299 zt(lev) = t(i) + (z(lev)-zz(i))/(zz(i+1)-zz(i))*(t(i+1)-t(i))
2300 if(abs(air(i)-air(i+1)).lt.air(i)/1.e5)then
2301 zair(lev) = air(i)
2302 else
2303 hlocal = 1./alog(air(i)/air(i+1))
2304 x0 = (z(lev)-zz(i))/(zz(i+1)+zz(i))
2305 zair(lev) = air(i)*exp(-x0/hlocal)
2306 endif
2307 vcld(lev-1) = cloud(i)
2308 ! write(06,'('o:z,t,air ',i3,3e12.4)')
2309 ! lev,z(lev),zt(lev),zair(lev)
2310 END DO
2312 20 CONTINUE
2313 END IF
2314 ! write(06,'('o:z,t,air ',i3,3e12.4)') lev,z(lev),zt(lev),zair(lev)
2315 30 CONTINUE
2317 ! number of levels including additional cloud levels
2319 nlevel = lev
2321 IF (nlevel>nj) print *, ' NLEVEL > NJ, NJ GROESSER WAEHLEN ', &
2322 nlevel
2324 ! write(06,'(' nlevel',i3)') nlevel
2327 ! assign default yields
2328 DO 40 nr = 1, nreakj
2330 DO 40 kl = kl0, kl1
2332 DO 40 lev = 1, nlevel
2333 40 qy(lev,kl,nr) = xqy(kl,nr)
2334 ! assign default absorption cross sections
2335 DO 50 kl = kl0, kl1
2337 DO 50 lev = 1, nlevel
2339 DO 50 nr = 1, nreakj
2340 50 s(lev,kl,nr) = xs(kl,nr)
2341 ! --------------------------------------------------------------------
2342 ! re-calculate altitude dependent quantum yields. this currently
2343 ! applies to
2344 ! 2=o3->o(1d)
2345 ! 12=ch2o->h2+co
2346 ! 13=ch3cho->ch3+cho
2347 ! 14=ch3coch3
2348 ! 15=ch3coch2ch3
2349 ! 16=hcocho -> 0.13 hcho + 1.87 co process a
2350 ! 17=ch3cocho
2351 ! 22=hcocho -> 0.45 hcho + 1.55 co + 0.80 ho2 process b
2353 ! o3 and ketones yield is calculated from fit equations,
2354 ! while for ch3cho and the dicarbonyls yields are calculated
2355 ! from the ntp yield by linear adjustment to 1/yield. the ch2o yield
2356 ! recalculated only for wavelengths longer than 329 nm. the yields for
2357 ! 1=o3->o(3p) are calculated as (1.- singlet d yield).
2359 DO 60 lev = 1, nlevel
2360 ! o3 ozone:
2361 !orig tau = zt(lev) - 230.
2362 !orig a = 0.9*(0.369+2.85E-4*tau+1.28E-5*tau*tau+2.57E-8*tau*tau*tau)
2363 !orig b = -0.575 + 5.59E-3*tau - 1.439E-5*tau*tau - 3.27E-8*tau*tau*tau
2364 !orig c = 0.9*(0.518+9.87E-4*tau-3.94E-5*tau*tau+3.91E-7*tau*tau*tau)
2365 !orig xl0 = 308.20 + 4.4871E-2*tau + 6.9380E-5*tau*tau - &
2366 !orig 2.5452E-6*tau*tau*tau
2368 DO 60 kl = kl0, kl1
2369 xl = wl(kl)
2370 !orig qy(lev,kl,2) = a*atan(b*(xl-xl0)) + c
2371 !orig IF (qy(lev,kl,2)<0.) qy(lev,kl,2) = 0.0
2372 !orig IF (qy(lev,kl,2)>0.9) qy(lev,kl,2) = 0.9
2373 qy(lev,kl,2) = 0.
2374 kt = 0.695 * zt(lev)
2375 q1 = 1.
2376 q2 = exp( -825.518/kt )
2377 if( xl <= 305. ) then
2378 qy(lev,kl,2) = .90
2379 else if( xl > 305. .and. xl <= 328. ) then
2380 tau = zt(lev)/300.
2381 qy(lev,kl,2) = 0.0765 &
2382 + 0.8036* (q1/(q1+q2))*exp( -((304.225 - xl)/5.576)**4 ) &
2383 + 8.9061*(tau)**2 *(q2/(q1+q2))*exp( -((314.957 - xl)/6.601)**2 ) &
2384 + 0.1192*(tau)**1.5 *exp( -((310.737 - xl)/2.187)**2 )
2385 else if( xl > 328. .and. xl <= 340. ) then
2386 qy(lev,kl,2) = 0.08
2387 else if( xl > 340. ) then
2388 qy(lev,kl,2) = 0.
2389 end if
2390 qy(lev,kl,3) = 1.0 - qy(lev,kl,2)
2391 ! ch2o formaldehyde:
2392 IF ((xl>=330.) .AND. qy(lev,kl,12)>0.) THEN
2393 phi1 = qy(lev,kl,11)
2394 phi2 = qy(lev,kl,12)
2395 phi20 = 1. - phi1
2396 ak300 = ((1./phi2)-(1./phi20))/2.54E+19
2397 akt = ak300*(1.+61.69*(1.-zt(lev)/300.)*(xl/329.-1.))
2398 qy(lev,kl,12) = 1./((1./phi20)+zair(lev)*akt)
2399 END IF
2401 IF (qy(lev,kl,12)>1.) qy(lev,kl,12) = 1.0
2403 IF (qy(lev,kl,12)<0.) qy(lev,kl,12) = 0.0
2404 ! ch3cho acetaldehyde:
2405 IF (xqy(kl,13)/=0.) THEN
2406 qy(lev,kl,13) = 1./(1.+(1./xqy(kl,13)-1.)*zair(lev)/2.465E19)
2407 END IF
2408 ! ch3coch3 acetone:
2409 qy(lev,kl,14) = 0.0766 + 0.09415*exp(-zair(lev)/3.222E18)
2410 ! ch3coch2ch3 methyl ethyl ketone:
2411 qy(lev,kl,15) = qy(lev,kl,14)
2412 ! hcocho glyoxal process a:
2413 IF (xqy(kl,16)/=0.) THEN
2414 qy(lev,kl,16) = 1./(1.+(1./xqy(kl,16)-1.)*zair(lev)/2.465E19)
2415 END IF
2416 ! ch3cocho methylglyoxal:
2417 IF (xqy(kl,17)/=0.) THEN
2418 qy(lev,kl,17) = 1./(1.+(1./xqy(kl,17)-1.)*zair(lev)/2.465E19)
2419 END IF
2420 ! hcocho glyoxal process b:
2421 IF (xqy(kl,22)/=0.) THEN
2422 qy(lev,kl,22) = 1./(1.+(1./xqy(kl,22)-1.)*zair(lev)/2.465E19)
2423 END IF
2425 60 CONTINUE
2426 ! _______________________________________________________________________
2427 ! correct absorption cross sections for t and p dep. for now, do
2428 ! 2=ozone
2429 !orig DO 90 kl = kl0, kl1
2430 DO 90 kl = 1, 130
2431 DO 80 lev = 1, nlevel
2432 !orig IF (kl<33 .OR. kl>61) GO TO 70
2433 !orig tdiffx = zt(lev) - 230.
2434 !orig s(lev,kl,2) = (so3tx(kl,1)+so3tx(kl,2)*tdiffx+so3tx(kl,3)*tdiffx* &
2435 !orig tdiffx)*1.0E-18
2436 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation
2437 s(lev,kl,2) = jpl218(kl) + &
2438 (jpl295(kl) - jpl218(kl))*(zt(lev) - 218.)/(295. - 218.)
2439 IF (zt(lev) .LE. 218.) s(lev,kl,2) = jpl218(kl)
2440 IF (zt(lev) .GE. 295.) s(lev,kl,2) = jpl295(kl)
2441 s(lev,kl,2)=s(lev,kl,2)*1.E-20
2442 s(lev,kl,3) = s(lev,kl,2)
2443 !orig 70 CONTINUE
2444 80 CONTINUE
2445 90 CONTINUE
2447 ! ---------------------------------------------- layers
2448 nlayer = nlevel - 1
2449 ! write(06,'(' Layers ',i3)') nlayer
2450 lay = 0
2452 DO 100 i = 1, nlayer
2453 lay = lay + 1
2454 dz = z(i+1) - z(i)
2455 zmid(lay) = z(i) + 0.5*dz
2456 vt(lay) = (zt(lay+1)+zt(lay))/2.
2457 vair(lay) = dz*1.E5*(zair(i+1)+zair(i))/2.
2458 100 CONTINUE
2460 ! vo3(lay) = dz*1.e5*(o3(i+1) + o3(i))/2. ! umr. dz in cm
2461 ! vcld(lay) = 0. *dz
2462 ! vaer(lay) = (aer(i+1)-aer(i))/alog(aer(i+1)/aer(i)) *dz ! bei
2464 CALL trapez(zz,o3,1,nn,zmid,vo3,1,nlayer,nn,nn,nj,nj)
2467 CALL trapez(zz,aer,1,nn,zmid,vaer,1,nlayer,nn,nn,nj,nj)
2469 DO 110 i = 1, nlayer
2470 lay = i
2471 dz = z(i+1) - z(i)
2472 ! write(06,'('layer ',i3,6e12.4)') lay,zmid(lay),vt(lay),
2473 ! & vair(lay)/dz,vo3(lay),vaer(lay),vcld(lay)
2474 dz = z(i+1) - z(i)
2476 ! umr. dz in cm
2478 vo3(i) = dz*1.E5*vo3(i)
2479 vcld(i) = vcld(i)*dz
2480 110 vaer(i) = vaer(i)*dz
2483 ! normalize aerosol optical depth to unity sum
2484 sum = 0.
2486 DO 120 lay = 1, nlayer
2487 sum = sum + vaer(lay)
2488 120 CONTINUE
2490 DO 130 lay = 1, nlayer
2491 vaer(lay) = vaer(lay)/sum
2492 130 CONTINUE
2494 ! calculated vertical column of o2 above the midpoint of each layer:
2495 ! want to use this for computing the average schumann-runge cross
2496 ! secti
2497 ! in each layer.
2498 ! so use half of current layer and half of previous higher layer
2501 cvo2(nlayer) = 0.2095*vair(nlayer)/2.
2503 DO 140 ii = 2, nlayer
2504 lay = nlayer - ii + 1
2505 cvo2(lay) = cvo2(lay+1) + 0.2095*(vair(lay)+vair(lay+1))/2.
2506 140 CONTINUE
2508 ! correct attenuation coefficients for pressure and/or temperature
2509 ! dep.
2510 ! for now do only ozone absorption.
2511 !orig DO 160 kl = kl0, kl1
2512 DO 160 kl = 1, 130
2513 DO 150 lay = 1, nlayer
2514 !orig tdiffx = vt(lay) - 230.
2515 !orig ao3(lay,kl) = xs(kl,2)
2516 !orig IF (kl>=33 .AND. kl<=61) ao3(lay,kl) = (so3tx(kl,1)+so3tx(kl,2)* &
2517 !orig tdiffx+so3tx(kl,3)*tdiffx*tdiffx)*1.0E-18
2518 ! SAMcKeen, 2/14/13 update O3 cross sections and quantum yields, from NASA-JPL-2011 recommendation
2519 ao3(lay,kl) = jpl218(kl) + &
2520 (jpl295(kl) - jpl218(kl))*(vt(lay) - 218.)/(295. - 218.)
2521 IF (vt(lay) .LE. 218.) ao3(lay,kl) = jpl218(kl)
2522 IF (vt(lay) .GE. 295.) ao3(lay,kl) = jpl295(kl)
2523 ao3(lay,kl)=ao3(lay,kl)*1.E-20
2524 150 CONTINUE
2525 160 CONTINUE
2527 ! write(06,'(' z ')')
2528 ! write(06,'(5e12.4 )') z
2529 ! write(06,'(' zmid ')')
2530 ! write(06,'(5e12.4 )') zmid
2531 ! write(06,'(' zt ')')
2532 ! write(06,'(5e12.4 )') zt
2533 ! write(06,'(' vt ')')
2534 ! write(06,'(5e12.4 )') vt
2535 ! write(06,'(' vo3 ')')
2536 ! write(06,'(5e12.4 )') vo3
2537 ! write(06,'(' vair ')')
2538 ! write(06,'(5e12.4 )') vair
2539 ! write(06,'(' vaer ')')
2540 ! write(06,'(5e12.4 )') vaer
2543 RETURN
2545 END SUBROUTINE subgrid
2547 SUBROUTINE optics(iprt, &
2548 zenith, a, dsdh, nid, &
2549 vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld, &
2550 omcld,vaer,aaer,gaer,omaer,nlevel, kzmax, &
2551 endir,endn,enup, fldir,fldn,flup)
2553 ! sm this subroutine prepares the data needed for the flux
2554 ! calculation,
2555 ! sm calls the scattering subroutine ps2str. it returns values of
2556 ! sm flux(lev,kl) for altitude lev-1, wavelength kl.
2557 ! sm it calculates the optical depths (vertical)
2558 ! sm for each layer, from the vertical profiles of o2, o3,
2559 ! sm air, cloud, and aerosol, and from the associated 'cross
2560 ! sections
2562 ! .. Scalar Arguments ..
2564 INTEGER :: iprt, nlevel
2565 REAL :: zenith
2566 REAL :: gaer, gcld, gray, omaer, omcld, omray
2568 ! .. Local Scalars ..
2569 REAL :: dtabs, dtaer, dtair, dtcld, dto2, dto3, dtscat
2570 INTEGER :: ii, lay, lev
2572 ! .. Intrinsic Functions ..
2573 INTRINSIC amax1
2575 ! .. Array Arguments ..
2577 INTEGER :: kzmax
2578 REAL :: dsdh(0:kzmax,kzmax)
2579 INTEGER :: nid(0:kzmax)
2580 REAL :: aaer(130), ao2(nj,130), ao3(nj,130), arayl(130), &
2581 vaer(nj), vair(nj), vcld(nj), vo3(nj)
2582 REAL :: fldir(nj,130), fldn(nj,130), flup(nj,130)
2583 REAL :: endir(nj,130), endn(nj,130), enup(nj,130)
2584 LOGICAL, PARAMETER :: delta = .TRUE.
2586 ! .. Local Arrays ..
2587 REAL :: alb, dtau(nj), g(nj), om(nj)
2588 REAL :: fdr(nj), fup(nj), fdn(nj), edr(nj), eup(nj), edn(nj)
2589 REAL, PARAMETER :: largest = 1.E+36
2593 ! loop over wavelengths
2594 DO 30 kl = kl0, kl1
2596 ! calculate optical depths for all layers (including cloud
2597 ! sublayers)
2599 DO 10 lay = 1, nlevel - 1
2600 ii = nlevel - lay
2602 dtair = vair(lay)*arayl(kl)
2603 dto2 = 0.2095*vair(lay)*ao2(lay,kl)
2604 dto3 = vo3(lay)*ao3(lay,kl)
2605 dtcld = vcld(lay)
2606 dtaer = vaer(lay)*aaer(kl)
2608 dtscat = dtair + dtcld + dtaer
2609 dtabs = dto2 + dto3
2611 dtscat = AMAX1(dtscat, 1./largest)
2612 dtabs = AMAX1(dtabs, 1./largest)
2614 dtau(ii) = dtabs + dtscat
2615 om(ii) = dtscat/dtau(ii)
2616 om(ii) = AMAX1(om(ii), 1./largest)
2618 g(ii) = (gray*dtair+gcld*dtcld+gaer*dtaer)/dtscat
2620 10 CONTINUE
2622 alb = albedoph(kl)
2624 ! ---------------------------------------------------------------------
2626 ! call pseudo-spherical 2-stream, set to delta-Eddington in subroutine.
2628 CALL ps2str(kzmax,nlevel,zenith,a,alb,dtau,om,g, &
2629 dsdh, nid, delta, &
2630 fdr, fup, fdn, edr, eup, edn)
2632 ! return to upright grid
2633 DO 20 ii = 1, nlevel
2634 lev = nlevel - ii + 1
2636 endir(lev,kl) = edr(ii)
2637 endn(lev,kl) = edn(ii)
2638 enup(lev,kl) = eup(ii)
2639 fldir(lev,kl) = fdr(ii)
2640 fldn(lev,kl) = fdn(ii)
2641 flup(lev,kl) = fup(ii)
2643 20 CONTINUE
2645 30 CONTINUE
2647 RETURN
2649 END SUBROUTINE optics
2651 ! #######################################################################
2653 SUBROUTINE calc_zenith(lat,long,ijd,gmt,azimuth,zenith)
2654 ! this subroutine calculates solar zenith and azimuth angles for a
2655 ! part
2656 ! time and location. must specify:
2657 ! input:
2658 ! lat - latitude in decimal degrees
2659 ! long - longitude in decimal degrees
2660 ! gmt - greenwich mean time - decimal military eg.
2661 ! 22.75 = 45 min after ten pm gmt
2662 ! output
2663 ! zenith
2664 ! azimuth
2665 ! .. Scalar Arguments ..
2666 REAL :: azimuth, gmt, lat, long, zenith
2667 INTEGER :: ijd
2668 ! .. Local Scalars ..
2669 REAL :: caz, csz, cw, d, decl, dr, ec, epsi, eqt, eyt, feqt, feqt1, &
2670 feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
2671 pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, &
2672 yt, zpt, zr
2673 INTEGER :: jd
2674 ! .. Intrinsic Functions ..
2675 INTRINSIC acos, atan, cos, min, sin, tan
2676 ! convert to radians
2677 pi = 3.1415926535590
2678 dr = pi/180.
2679 rlt = lat*dr
2680 rphi = long*dr
2682 ! print julian days current 'ijd'
2684 ! ???? + (yr - yref)
2686 jd = ijd
2690 d = jd + gmt/24.0
2691 ! calc geom mean longitude
2692 ml = 279.2801988 + .9856473354*d + 2.267E-13*d*d
2693 rml = ml*dr
2695 ! calc equation of time in sec
2696 ! w = mean long of perigee
2697 ! e = eccentricity
2698 ! epsi = mean obliquity of ecliptic
2699 w = 282.4932328 + 4.70684E-5*d + 3.39E-13*d*d
2700 wr = w*dr
2701 ec = 1.6720041E-2 - 1.1444E-9*d - 9.4E-17*d*d
2702 epsi = 23.44266511 - 3.5626E-7*d - 1.23E-15*d*d
2703 pepsi = epsi*dr
2704 yt = (tan(pepsi/2.0))**2
2705 cw = cos(wr)
2706 sw = sin(wr)
2707 ssw = sin(2.0*wr)
2708 eyt = 2.*ec*yt
2709 feqt1 = sin(rml)*(-eyt*cw-2.*ec*cw)
2710 feqt2 = cos(rml)*(2.*ec*sw-eyt*sw)
2711 feqt3 = sin(2.*rml)*(yt-(5.*ec**2/4.)*(cw**2-sw**2))
2712 feqt4 = cos(2.*rml)*(5.*ec**2*ssw/4.)
2713 feqt5 = sin(3.*rml)*(eyt*cw)
2714 feqt6 = cos(3.*rml)*(-eyt*sw)
2715 feqt7 = -sin(4.*rml)*(.5*yt**2)
2716 feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
2717 eqt = feqt*13751.0
2719 ! convert eq of time from sec to deg
2720 reqt = eqt/240.
2721 ! calc right ascension in rads
2722 ra = ml - reqt
2723 rra = ra*dr
2724 ! calc declination in rads, deg
2725 tab = 0.43360*sin(rra)
2726 rdecl = atan(tab)
2727 decl = rdecl/dr
2728 ! calc local hour angle
2729 lbgmt = 12.0 - eqt/3600. + long*24./360.
2730 lzgmt = 15.0*(gmt-lbgmt)
2731 zpt = lzgmt*dr
2732 csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
2733 if(csz.gt.1)print *,'calczen,csz ',csz
2734 csz = min(1.,csz)
2735 ! zr = acos(csz)
2736 ! zenith = zr/dr
2737 zr = acos(csz)
2738 zenith = zr/dr
2739 ! calc local solar azimuth
2740 caz = (sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr))
2741 if(caz.lt.-0.999999)then
2742 azimuth=180.
2743 elseif(caz.gt.0.999999)then
2744 azimuth=0.
2745 else
2746 raz = acos(caz)
2747 azimuth = raz/dr
2748 endif
2749 ! caz = min(1.,(sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr)))
2750 ! if(caz.lt.-1)print *,'calczen ',caz
2751 ! caz = max(-1.,caz)
2752 ! raz = acos(caz)
2753 ! azimuth = raz/dr
2755 IF (lzgmt>0) azimuth = azimuth + (2*(180.-azimuth))
2756 ! 200 format(' ',f7.2,2(12x,f7.2))
2757 RETURN
2759 END SUBROUTINE calc_zenith
2761 SUBROUTINE trapez(x,y,ianfa,ienda,u,v,ianfn,iendn,ix,iy,iu,iv)
2762 ! * implemented 1992 by ansgar ruggaber, university of munich, frg
2763 ! * funded by the german minister of research and technology (bmft)
2764 ! * under contract no. 521-4007-07eu-738 8
2765 ! lineare interpolation of referencefield (x(i),y(i))
2766 ! to (u(i),v(i))
2767 ! save
2768 ! .. Scalar Arguments ..
2769 INTEGER :: ianfa, ianfn, ienda, iendn, iu, iv, ix, iy
2770 ! .. Array Arguments ..
2771 REAL :: u(iu), v(iv), x(ix), y(iy)
2772 ! .. Local Scalars ..
2773 REAL :: uumord, vumord, xumord, yumord
2774 INTEGER :: i, ianf, ianfnn, idrehu, idrehx, iendnn, iordu, iordx, j
2776 idrehx = 0
2778 IF (x(ianfa)>=x(ienda)) THEN
2779 ! das x-feld wird ansteigend geordnet, das y-feld entsprechend
2780 iordx = (ienda-ianfa+1)/2
2782 DO i = ianfa, iordx
2783 xumord = x(i)
2784 x(i) = x(ienda+1-i)
2785 x(ienda+1-i) = xumord
2786 yumord = y(i)
2787 y(i) = y(ienda+1-i)
2788 y(ienda+1-i) = yumord
2789 END DO
2791 idrehx = 1
2792 END IF
2794 idrehu = 0
2796 IF (u(ianfn)>=u(iendn)) THEN
2797 ! u-field increasing
2798 iordu = (iendn-ianfn+1)/2
2800 DO i = ianfn, iordu
2801 uumord = u(i)
2802 u(i) = u(iendn+1-i)
2803 u(iendn+1-i) = uumord
2804 END DO
2806 idrehu = 1
2807 END IF
2809 ianfnn = ianfn
2810 10 CONTINUE
2812 IF (u(ianfnn)<x(ianfa)) THEN
2813 ! no extrapolation at x(ianfa)
2814 v(ianfnn) = 1.0E-12
2815 ianfnn = ianfnn + 1
2816 GO TO 10
2817 END IF
2819 iendnn = iendn
2820 20 CONTINUE
2822 IF ((1.-1.0E-10)*u(iendnn)>x(ienda)) THEN
2823 ! no extrapolation at x(ienda)
2824 v(iendnn) = 1.0E-12
2825 iendnn = iendnn - 1
2826 GO TO 20
2827 END IF
2829 ianf = ianfa
2831 DO j = ianfnn, iendnn
2833 DO i = ianf, ienda
2835 IF (x(i)-u(j)) 30, 50, 40
2836 30 END DO
2838 GO TO 70
2839 40 v(j) = y(i-1) + (y(i)-y(i-1))/(x(i)-x(i-1))*(u(j)-x(i-1))
2840 GO TO 60
2841 50 v(j) = y(i)
2842 60 ianf = i
2843 70 END DO
2845 IF (idrehx/=0) THEN
2846 ! x- und y-field in starting position
2847 DO i = ianfa, iordx
2848 xumord = x(i)
2849 x(i) = x(ienda+1-i)
2850 x(ienda+1-i) = xumord
2851 yumord = y(i)
2852 y(i) = y(ienda+1-i)
2853 y(ienda+1-i) = yumord
2854 END DO
2856 END IF
2858 IF (idrehu/=0) THEN
2860 DO i = ianfn, iordu
2861 uumord = u(i)
2862 u(i) = u(iendn+1-i)
2863 u(iendn+1-i) = uumord
2864 vumord = v(i)
2865 v(i) = v(iendn+1-i)
2866 v(iendn+1-i) = vumord
2867 END DO
2869 END IF
2871 END SUBROUTINE trapez
2873 SUBROUTINE photmad_init(z_at_w,aerwrf,g,ids,ide,jds,jde,kds,kde,ims,ime, &
2874 jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
2875 ! local stuff
2876 ! .. Scalar Arguments ..
2877 REAL, INTENT (IN) :: g
2878 INTEGER, INTENT (IN) :: ide, ids, ime, ims, ite, its, jde, jds, jme, &
2879 jms, jte, jts, kde, kds, kme, kms, kte, kts
2880 ! .. Array Arguments ..
2881 REAL, INTENT (INOUT) :: aerwrf(ims:ime,kms:kme,jms:jme)
2882 REAL, INTENT (IN) :: z_at_w(ims:ime,kms:kme,jms:jme)
2883 ! .. Local Scalars ..
2884 REAL :: z1
2885 INTEGER :: i, j, k
2886 ! .. Local Arrays ..
2887 REAL :: aerext(kts:kte), phizz(kts:kte), z(kts:kte)
2889 DO j = jts, jte
2891 IF (j>jde-1) GO TO 20
2893 DO i = its, ite
2895 IF (i>ide-1) GO TO 10
2897 ! z at w points
2899 z1 = z_at_w(i,kts,j)
2900 z(kts)=0.
2902 DO k = kts+1, kte
2903 z(k) = z_at_w(i,k,j) - z1
2904 END DO
2906 DO k = kts, kte
2907 phizz(k) = .001*z(k)
2908 aerext(k) = 0.
2909 ! if(i.eq.its.and.j.eq.jts)print *,phizz(k),aerstd(k),kts,kte
2910 ! print *,phizz(k),kts,kte,ite,jte
2911 END DO
2913 ! IF (phizz(kte-1)>20.) THEN
2914 ! CALL wrf_error_fatal ( 'phizz(kte-1)>20., set kl0 to 1')
2915 ! END IF
2917 CALL trapez(zstd,aerstd,1,51,phizz,aerext,kts,kte,51,51,kte,kte)
2919 DO k = kts, kte
2920 aerwrf(i,k,j) = aerext(k)
2921 ! if(i.eq.its)print *,k,i,j,aerext(k),phizz(k)
2922 ! print *,k,i,j,aerext(k),phizz(k)
2923 END DO
2925 10 CONTINUE
2926 END DO
2928 20 CONTINUE
2929 END DO
2931 END SUBROUTINE photmad_init
2933 ! ===========================================================================
2935 SUBROUTINE ps2str(kzmax, kz,zen,mu,rsfc,tauu,omu,gu, &
2936 dsdh, nid, delta, &
2937 fdr, fup, fdn, edr, eup, edn)
2939 ! ----------------------------------------------------------------------------
2940 ! PURPOSE:
2941 ! Solve two-stream equations for multiple layers. The subroutine is based
2942 ! on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.
2943 ! It contains 9 two-stream methods to choose from. A pseudo-spherical
2944 ! correction has also been added.
2945 ! ----------------------------------------------------------------------------
2946 ! PARAMETERS:
2947 ! KZ - INTEGER, number of specified altitude levels in the working (I)
2948 ! grid
2949 ! ZEN - REAL, solar zenith angle (degrees) (I)
2950 ! RSFC - REAL, surface albedo at current wavelength (I)
2951 ! TAUU - REAL, unscaled optical depth of each layer (I)
2952 ! OMU - REAL, unscaled single scattering albedo of each layer (I)
2953 ! GU - REAL, unscaled asymmetry parameter of each layer (I)
2954 ! DSDH - REAL, slant path of direct beam through each layer crossed (I)
2955 ! when travelling from the top of the atmosphere to layer i;
2956 ! DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1
2957 ! NID - INTEGER, number of layers crossed by the direct beam when (I)
2958 ! travelling from the top of the atmosphere to layer i;
2959 ! NID(i), i = 0..NZ-1
2960 ! DELTA - LOGICAL, switch to use delta-scaling (I)
2961 ! .TRUE. -> apply delta-scaling
2962 ! .FALSE.-> do not apply delta-scaling
2963 ! FDR - REAL, contribution of the direct component to the total (O)
2964 ! actinic flux at each altitude level
2965 ! FUP - REAL, contribution of the diffuse upwelling component to (O)
2966 ! the total actinic flux at each altitude level
2967 ! FDN - REAL, contribution of the diffuse downwelling component to (O)
2968 ! the total actinic flux at each altitude level
2969 ! EDR - REAL, contribution of the direct component to the total (O)
2970 ! spectral irradiance at each altitude level
2971 ! EUP - REAL, contribution of the diffuse upwelling component to (O)
2972 ! the total spectral irradiance at each altitude level
2973 ! EDN - REAL, contribution of the diffuse downwelling component to (O)
2974 ! the total spectral irradiance at each altitude level
2975 ! ----------------------------------------------------------------------------
2977 IMPLICIT NONE
2979 !******
2980 ! input:
2981 !******
2982 INTEGER kz
2983 REAL zen, rsfc
2984 REAL tauu(kz), omu(kz), gu(kz)
2985 INTEGER kzmax
2986 REAL dsdh(0:kzmax,kzmax)
2987 INTEGER nid(0:kzmax)
2988 LOGICAL delta
2990 !******
2991 ! output:
2992 !******
2993 REAL fup(kz),fdn(kz),fdr(kz)
2994 REAL eup(kz),edn(kz),edr(kz)
2996 !******
2997 ! local:
2998 !******
2999 INTEGER nlevel
3000 REAL tausla(0:kz), tauc(0:kz)
3001 REAL mu2(0:kz), mu, sum
3003 ! internal coefficients and matrix
3005 REAL lam(kz),taun(kz),bgam(kz)
3006 REAL e1(kz),e2(kz),e3(kz),e4(kz)
3007 REAL cup(kz),cdn(kz),cuptn(kz),cdntn(kz)
3008 REAL mu1(kz)
3009 INTEGER row
3010 REAL a(2*kz),b(2*kz),d(2*kz),e(2*kz),y(2*kz)
3012 !******
3013 ! other:
3014 !******
3015 REAL pifs, fdn0
3016 REAL gi(kz), omi(kz), tempg
3017 REAL f, g, om
3018 REAL gam1, gam2, gam3, gam4
3020 ! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
3021 ! in delta-function, modified quadrature, hemispheric constant,
3022 ! Hybrid modified Eddington-delta function metods, p633,Table1.
3023 ! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
3024 ! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
3025 ! uncomment the following two lines and the appropriate statements further
3026 ! down.
3027 ! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
3028 ! > BETA1, BETAn, amu1, subd
3030 REAL expon, expon0, expon1, divisr, temp, up, dn
3031 REAL ssfc
3032 INTEGER nlayer, mrows, lev
3034 INTEGER i, j
3036 ! Some additional program constants:
3038 REAL pi
3039 PARAMETER (pi=3.1415926535898)
3040 REAL largest
3041 PARAMETER(largest=1.E+36)
3042 REAL precis
3043 PARAMETER(precis = 1.e-7)
3044 REAL eps
3045 PARAMETER (eps = 1.E-3)
3046 !_______________________________________________________________________
3048 ! MU = cosine of solar zenith angle
3049 ! RSFC = surface albedo
3050 ! TAUU = unscaled optical depth of each layer
3051 ! OMU = unscaled single scattering albedo
3052 ! GU = unscaled asymmetry factor
3053 ! KLEV = max dimension of number of layers in atmosphere
3054 ! NLAYER = number of layers in the atmosphere
3055 ! NLEVEL = nlayer + 1 = number of levels
3057 ! initial conditions: pi*solar flux = 1; diffuse incidence = 0
3059 nlevel = kz
3060 pifs = 1.
3061 fdn0 = 0.
3062 nlayer = nlevel - 1
3063 ! mu = COS(zen*pi/180.)
3065 !************* compute coefficients for each layer:
3066 ! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
3067 ! EXPON0 = calculation of e when TAU is zero
3068 ! EXPON1 = calculation of e when TAU is TAUN
3069 ! CUP and CDN = calculation when TAU is zero
3070 ! CUPTN and CDNTN = calc. when TAU is TAUN
3071 ! DIVISR = prevents division by zero
3073 DO j = 0, kz
3074 tauc(j) = 0.
3075 tausla(j) = 0.
3076 mu2(j) = 1./SQRT(largest)
3077 END DO
3079 IF( .NOT. delta ) THEN
3080 DO i = 1, nlayer
3081 gi(i) = gu(i)
3082 omi(i) = omu(i)
3083 taun(i) = tauu(i)
3084 ENDDO
3085 ELSE
3087 ! delta-scaling. Has to be done for delta-Eddington approximation,
3088 ! delta discrete ordinate, Practical Improved Flux Method, delta function,
3089 ! and Hybrid modified Eddington-delta function methods approximations
3091 DO i = 1, nlayer
3092 f = gu(i)*gu(i)
3093 gi(i) = (gu(i) - f)/(1 - f)
3094 omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)
3095 taun(i) = (1 - omu(i)*f)*tauu(i)
3096 ENDDO
3097 END IF
3100 ! calculate slant optical depth at the top of the atmosphere when zen>90.
3101 ! in this case, higher altitude of the top layer is recommended which can
3102 ! be easily changed in gridz.f.
3104 IF(zen .GT. 90.0) THEN
3105 IF(nid(0) .LT. 0) THEN
3106 tausla(0) = largest
3107 ELSE
3108 sum = 0.0
3109 DO j = 1, nid(0)
3110 sum = sum + 2.*taun(j)*dsdh(0,j)
3111 END DO
3112 tausla(0) = sum
3113 END IF
3114 END IF
3117 DO 11, i = 1, nlayer
3119 g = gi(i)
3120 om = omi(i)
3121 tauc(i) = tauc(i-1) + taun(i)
3123 ! stay away from 1 by precision. For g, also stay away from -1
3125 tempg = AMIN1(abs(g),1. - precis)
3126 g = SIGN(tempg,g)
3127 om = AMIN1(om,1.-precis)
3130 ! calculate slant optical depth
3131 !
3132 IF(nid(i) .LT. 0) THEN
3133 tausla(i) = largest
3134 ELSE
3135 sum = 0.0
3136 DO j = 1, MIN(nid(i),i)
3137 sum = sum + taun(j)*dsdh(i,j)
3138 ENDDO
3139 DO j = MIN(nid(i),i)+1,nid(i)
3140 sum = sum + 2.*taun(j)*dsdh(i,j)
3141 ENDDO
3142 tausla(i) = sum
3143 IF(tausla(i) .EQ. tausla(i-1)) THEN
3144 mu2(i) = SQRT(largest)
3145 ELSE
3146 mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
3147 mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)), &
3148 mu2(i) )
3149 END IF
3150 END IF
3152 !** the following gamma equations are from pg 16,289, Table 1
3153 !** save mu1 for each approx. for use in converting irradiance to actinic flux
3155 ! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
3157 gam1 = (7. - om*(4. + 3.*g))/4.
3158 gam2 = -(1. - om*(4. - 3.*g))/4.
3159 gam3 = (2. - 3.*g*mu)/4.
3160 gam4 = 1. - gam3
3161 mu1(i) = 0.5
3163 ! quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
3165 ! gam1 = 1.7320508*(2. - om*(1. + g))/2.
3166 ! gam2 = 1.7320508*om*(1. - g)/2.
3167 ! gam3 = (1. - 1.7320508*g*mu)/2.
3168 ! gam4 = 1. - gam3
3169 ! mu1(i) = 1./sqrt(3.)
3171 ! hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
3173 ! gam1 = 2. - om*(1. + g)
3174 ! gam2 = om*(1. - g)
3175 ! gam3 = (2. - g*mu)/4.
3176 ! gam4 = 1. - gam3
3177 ! mu1(i) = 0.5
3179 ! PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
3180 ! GAM1 = 0.25*(8. - OM*(5. + 3.*G))
3181 ! GAM2 = 0.75*OM*(1.-G)
3182 ! GAM3 = 0.25*(2.-3.*G*MU)
3183 ! GAM4 = 1. - GAM3
3184 ! mu1(i) = 0.5
3186 ! delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
3187 ! GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
3188 ! GAM2 = 0.5*1.7320508*OM*(1.-G)
3189 ! GAM3 = 0.5*(1.-1.7320508*G*MU)
3190 ! GAM4 = 1. - GAM3
3191 ! mu1(i) = 1./sqrt(3.)
3193 ! Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
3194 ! in delta-function, modified quadrature, hemispheric constant,
3195 ! Hybrid modified Eddington-delta function metods, p633,Table1.
3196 ! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
3197 ! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
3198 ! YLM0 = 2.
3199 ! YLM2 = -3.*G*MU
3200 ! YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
3201 ! YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
3202 ! YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
3203 ! *+429.*MU**6)
3204 ! YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
3205 ! *-25740.*MU**6+12155.*MU**8)
3206 ! YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
3207 ! *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
3208 ! YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
3209 ! YLMS=0.25*YLMS
3210 ! BETA0 = YLMS
3212 ! amu1=1./1.7320508
3213 ! YLM0 = 2.
3214 ! YLM2 = -3.*G*amu1
3215 ! YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
3216 ! YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
3217 ! YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
3218 ! *+429.*amu1**6)
3219 ! YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
3220 ! *-25740.*amu1**6+12155.*amu1**8)
3221 ! YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
3222 ! *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
3223 ! YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
3224 ! YLMS=0.25*YLMS
3225 ! BETA1 = YLMS
3227 ! BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
3228 ! *-0.045776*G**7)
3231 ! Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
3232 ! subd=4.*(1.-G*G*(1.-MU))
3233 ! GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
3234 ! GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
3235 ! GAM3 = BETA0
3236 ! GAM4 = 1. - GAM3
3237 ! mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)
3239 !****
3240 ! delta function (Meador, and Weaver, 1980, JAS, 37, 630):
3241 ! GAM1 = (1. - OM*(1. - beta0))/MU
3242 ! GAM2 = OM*BETA0/MU
3243 ! GAM3 = BETA0
3244 ! GAM4 = 1. - GAM3
3245 ! mu1(i) = mu
3246 !****
3247 ! modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
3248 ! GAM1 = 1.7320508*(1. - OM*(1. - beta1))
3249 ! GAM2 = 1.7320508*OM*beta1
3250 ! GAM3 = BETA0
3251 ! GAM4 = 1. - GAM3
3252 ! mu1(i) = 1./sqrt(3.)
3254 ! hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
3255 ! GAM1 = 2.*(1. - OM*(1. - betan))
3256 ! GAM2 = 2.*OM*BETAn
3257 ! GAM3 = BETA0
3258 ! GAM4 = 1. - GAM3
3259 ! mu1(i) = 0.5
3261 !****
3263 ! lambda = pg 16,290 equation 21
3264 ! big gamma = pg 16,290 equation 22
3265 ! if gam2 = 0., then bgam = 0.
3267 lam(i) = sqrt(gam1*gam1 - gam2*gam2)
3269 IF( gam2 .NE. 0.) THEN
3270 bgam(i) = (gam1 - lam(i))/gam2
3271 ELSE
3272 bgam(i) = 0.
3273 ENDIF
3275 expon = EXP(-lam(i)*taun(i))
3277 ! e1 - e4 = pg 16,292 equation 44
3279 e1(i) = 1. + bgam(i)*expon
3280 e2(i) = 1. - bgam(i)*expon
3281 e3(i) = bgam(i) + expon
3282 e4(i) = bgam(i) - expon
3284 ! the following sets up for the C equations 23, and 24
3285 ! found on page 16,290
3286 ! prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
3287 ! which is approx equiv to shifting MU by 0.5*EPS* (MU)**3
3289 expon0 = EXP(-tausla(i-1))
3290 expon1 = EXP(-tausla(i))
3292 divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
3293 temp = AMAX1(eps,abs(divisr))
3294 divisr = SIGN(temp,divisr)
3296 up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
3297 dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
3299 ! cup and cdn are when tau is equal to zero
3300 ! cuptn and cdntn are when tau is equal to taun
3302 cup(i) = up*expon0
3303 cdn(i) = dn*expon0
3304 cuptn(i) = up*expon1
3305 cdntn(i) = dn*expon1
3307 11 CONTINUE
3309 !**************** set up matrix ******
3310 ! ssfc = pg 16,292 equation 37 where pi Fs is one (unity).
3312 ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs
3314 ! MROWS = the number of rows in the matrix
3316 mrows = 2*nlayer
3318 ! the following are from pg 16,292 equations 39 - 43.
3319 ! set up first row of matrix:
3321 i = 1
3322 a(1) = 0.
3323 b(1) = e1(i)
3324 d(1) = -e2(i)
3325 e(1) = fdn0 - cdn(i)
3327 row=1
3329 ! set up odd rows 3 thru (MROWS - 1):
3330 i = 0
3331 DO 20, row = 3, mrows - 1, 2
3332 i = i + 1
3333 a(row) = e2(i)*e3(i) - e4(i)*e1(i)
3334 b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
3335 d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
3336 e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + &
3337 e1(i)*(cdntn(i) - cdn(i + 1))
3338 20 CONTINUE
3340 ! set up even rows 2 thru (MROWS - 2):
3342 i = 0
3343 DO 30, row = 2, mrows - 2, 2
3344 i = i + 1
3345 a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
3346 b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
3347 d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
3348 e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - &
3349 (cdn(i + 1) - cdntn(i))*e4(i + 1)
3350 30 CONTINUE
3352 ! set up last row of matrix at MROWS:
3354 row = mrows
3355 i = nlayer
3357 a(row) = e1(i) - rsfc*e3(i)
3358 b(row) = e2(i) - rsfc*e4(i)
3359 d(row) = 0.
3360 e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
3362 ! solve tri-diagonal matrix:
3364 CALL tridag(a, b, d, e, y, mrows)
3366 !*** unfold solution of matrix, compute output fluxes:
3368 row = 1
3369 lev = 1
3370 j = 1
3372 ! the following equations are from pg 16,291 equations 31 & 32
3374 fdr(lev) = EXP( -tausla(0) )
3375 edr(lev) = mu * fdr(lev)
3376 edn(lev) = fdn0
3377 eup(lev) = y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
3378 fdn(lev) = edn(lev)/mu1(lev)
3379 fup(lev) = eup(lev)/mu1(lev)
3381 DO 60, lev = 2, nlayer + 1
3382 fdr(lev) = EXP(-tausla(lev-1))
3383 edr(lev) = mu *fdr(lev)
3384 edn(lev) = y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
3385 eup(lev) = y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
3386 fdn(lev) = edn(lev)/mu1(j)
3387 fup(lev) = eup(lev)/mu1(j)
3389 row = row + 2
3390 j = j + 1
3391 60 CONTINUE
3392 !_______________________________________________________________________
3394 END SUBROUTINE ps2str
3396 ! ===========================================================================
3398 SUBROUTINE tridag(a,b,c,r,u,n)
3400 !_______________________________________________________________________
3401 ! solves tridiagonal system. From Numerical Recipies, p. 40
3402 !_______________________________________________________________________
3404 IMPLICIT NONE
3406 ! input:
3407 INTEGER n
3408 REAL a, b, c, r
3409 DIMENSION a(n),b(n),c(n),r(n)
3411 ! output:
3412 REAL u
3413 DIMENSION u(n)
3415 ! local:
3416 INTEGER j
3418 REAL bet, gam
3419 DIMENSION gam(2*n)
3420 !_______________________________________________________________________
3422 IF (b(1) .EQ. 0.) STOP 'subroutine tridag failed, 1001'
3423 bet = b(1)
3424 u(1) = r(1)/bet
3425 DO 11, j = 2, n
3426 gam(j) = c(j - 1)/bet
3427 bet = b(j) - a(j)*gam(j)
3428 IF (bet .EQ. 0.) STOP 'subroutine tridag failed, 2002'
3429 u(j) = (r(j) - a(j)*u(j - 1))/bet
3430 11 CONTINUE
3431 DO 12, j = n - 1, 1, -1
3432 u(j) = u(j) - gam(j + 1)*u(j + 1)
3433 12 CONTINUE
3434 !_______________________________________________________________________
3436 END SUBROUTINE tridag
3438 ! ===========================================================================
3440 SUBROUTINE sphers(kz, nz, z, zen, dsdh, nid)
3442 ! ----------------------------------------------------------------------------
3443 ! PURPOSE:
3444 ! Calculate slant path over vertical depth ds/dh in spherical geometry.
3445 ! Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model
3446 ! for computing the radiation field available for photolysis and heating
3447 ! at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)
3448 ! ----------------------------------------------------------------------------
3449 ! PARAMETERS:
3450 ! NZ - INTEGER, number of specified altitude levels in the working (I)
3451 ! grid
3452 ! Z - REAL, specified altitude working grid (km) (I)
3453 ! ZEN - REAL, solar zenith angle (degrees) (I)
3454 ! DSDH - REAL, slant path of direct beam through each layer crossed (O)
3455 ! when travelling from the top of the atmosphere to layer i;
3456 ! DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1
3457 ! NID - INTEGER, number of layers crossed by the direct beam when (O)
3458 ! travelling from the top of the atmosphere to layer i;
3459 ! NID(i), i = 0..NZ-1
3460 ! ----------------------------------------------------------------------------
3461 ! EDIT HISTORY:
3462 ! double precision fix for shallow layers - Julia Lee-Taylor Dec 2000
3463 ! ----------------------------------------------------------------------------
3464 ! This program is free software; you can redistribute it and/or modify
3465 ! it under the terms of the GNU General Public License as published by the
3466 ! Free Software Foundation; either version 2 of the license, or (at your
3467 ! option) any later version.
3468 ! The TUV package is distributed in the hope that it will be useful, but
3469 ! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTIBI-
3470 ! LITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
3471 ! License for more details.
3472 ! To obtain a copy of the GNU General Public License, write to:
3473 ! Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
3474 ! ----------------------------------------------------------------------------
3475 ! To contact the authors, please mail to:
3476 ! Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA or
3477 ! send email to: sasha@ucar.edu
3478 ! ----------------------------------------------------------------------------
3479 ! Copyright (C) 1994-2008 University Corporation for Atmospheric Research
3480 ! ----------------------------------------------------------------------------
3482 IMPLICIT NONE
3484 ! input
3485 INTEGER kz, nz
3486 REAL zen, z(kz)
3488 ! output
3489 INTEGER nid(0:kz)
3490 REAL dsdh(0:kz,kz)
3492 ! more program constants
3493 ! pi:
3494 REAL pi
3495 PARAMETER(pi=3.1415926535898)
3496 REAL dr
3497 PARAMETER ( dr = pi/180.)
3499 ! radius of the earth:
3500 REAL re
3501 REAL radius
3502 PARAMETER(radius=6.371E+3)
3504 REAL ze(kz)
3506 ! local
3508 DOUBLE PRECISION zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
3509 INTEGER i, j, k
3510 INTEGER id
3512 INTEGER nlayer
3513 REAL zd(0:kz-1)
3515 ! ----------------------------------------------------------------------------
3517 zenrad = zen*dr
3519 ! number of layers:
3520 nlayer = nz - 1
3522 ! include the elevation above sea level to the radius of the earth:
3523 re = radius + z(1)
3524 ! correspondingly z changed to the elevation above earth surface:
3525 DO k = 1, nz
3526 ze(k) = z(k) - z(1)
3527 END DO
3529 ! inverse coordinate of z
3530 zd(0) = ze(nz)
3531 DO k = 1, nlayer
3532 zd(k) = ze(nz - k)
3533 END DO
3535 ! initialize dsdh(i,j), nid(i)
3536 DO i = 0, kz
3537 nid(i) = 0
3538 DO j = 1, kz
3539 dsdh(i,j) = 0.
3540 END DO
3541 END DO
3543 ! calculate ds/dh of every layer
3544 DO 100 i = 0, nlayer
3546 rpsinz = (re + zd(i)) * SIN(zenrad)
3548 IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
3549 nid(i) = -1
3550 ELSE
3553 ! Find index of layer in which the screening height lies
3555 id = i
3556 IF( zen .GT. 90.0 ) THEN
3557 DO 10 j = 1, nlayer
3558 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND. &
3559 (rpsinz .GE. ( zd(j) + re )) ) id = j
3560 10 CONTINUE
3561 END IF
3563 DO 20 j = 1, id
3565 sm = 1.0
3566 IF(j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) &
3567 sm = -1.0
3569 rj = re + zd(j-1)
3570 rjp1 = re + zd(j)
3572 dhj = zd(j-1) - zd(j)
3574 ga = rj*rj - rpsinz*rpsinz
3575 gb = rjp1*rjp1 - rpsinz*rpsinz
3576 IF (ga .LT. 0.0) ga = 0.0
3577 IF (gb .LT. 0.0) gb = 0.0
3579 IF(id.GT.i .AND. j.EQ.id) THEN
3580 dsj = SQRT( ga )
3581 ELSE
3582 dsj = SQRT( ga ) - sm*SQRT( gb )
3583 END IF
3584 dsdh(i,j) = dsj / dhj
3585 20 CONTINUE
3587 nid(i) = id
3589 END IF
3591 100 CONTINUE
3593 ! ----------------------------------------------------------------------------
3595 END SUBROUTINE sphers
3597 ! ===========================================================================
3599 SUBROUTINE airmas(kz, nz, dsdh, nid, cz, &
3600 vcol, scol)
3602 ! ----------------------------------------------------------------------------
3603 ! PURPOSE:
3604 ! Calculate vertical and slant air columns, in spherical geometry, as a
3605 ! function of altitude.
3606 ! ----------------------------------------------------------------------------
3607 ! PARAMETERS:
3608 ! kZ - INTEGER, number of specified altitude levels in the working (I)
3609 ! grid
3610 ! DSDH - REAL, slant path of direct beam through each layer crossed (O)
3611 ! when travelling from the top of the atmosphere to layer i;
3612 ! DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1
3613 ! NID - INTEGER, number of layers crossed by the direct beam when (O)
3614 ! travelling from the top of the atmosphere to layer i;
3615 ! NID(i), i = 0..NZ-1
3616 ! VCOL - REAL, output, vertical air column, molec cm-2, above level iz
3617 ! SCOL - REAL, output, slant air column in direction of sun, above iz
3618 ! also in molec cm-2
3619 ! ----------------------------------------------------------------------------
3621 IMPLICIT NONE
3623 ! Input:
3625 INTEGER kz, nz
3626 INTEGER nid(0:kz)
3627 REAL dsdh(0:kz,kz)
3628 REAL cz(nz)
3630 ! output:
3632 REAL vcol(nz), scol(nz)
3634 ! internal:
3636 INTEGER id, j
3637 REAL sum, vsum
3638 REAL largest
3639 PARAMETER(largest=1.E+36)
3641 ! calculate vertical and slant column from each level:
3642 ! work downward
3644 vsum = 0.
3645 DO id = 0, nz - 1
3646 vsum = vsum + cz(nz-id)
3647 vcol(nz-id) = vsum
3648 sum = 0.
3649 IF(nid(id) .LT. 0) THEN
3650 sum = largest
3651 ELSE
3653 ! single pass layers:
3655 DO j = 1, MIN(nid(id), id)
3656 sum = sum + cz(nz-j)*dsdh(id,j)
3657 ENDDO
3659 ! double pass layers:
3661 DO j = MIN(nid(id),id)+1, nid(id)
3662 sum = sum + 2.*cz(nz-j)*dsdh(id,j)
3663 ENDDO
3665 ENDIF
3666 scol(nz - id) = sum
3668 ENDDO
3670 END SUBROUTINE airmas
3672 END MODULE module_phot_mad