module_phot_mad.F
References to this file elsewhere.
1
2 MODULE module_phot_mad
3 ! preliminary fixed values for T, p, O3 and caer
4 ! use values in correct units!!
5 ! so3t(kl,i) - ozone cross sect. temperature dependence coefficients
6 ! wl(kl) - array of nominal center wavelengths of spectral intervals
7 ! f(kl) - extraterrestrial solar irradiance
8 ! xs(kl,nr) - cross sections for nr species.
9 ! xqy(kl,nr) - quantum yields. some are read, others are computed.
10 ! schumann-runge - part a
11 ! schumann-runge parameters - part b
12 ! ozone temperature coefficients for xsection
13 ! -----------------------------------------------------------------------
14 ! read stand profiles
15 ! aerosol
16 ! close (incvol)
17 ! wavelengths and extraterrestrial irradiance
18 ! -----------------------------------------------------
19 ! wave_length
20 ! WAVE LENGTHS USED BY PHOTOLYSIS PROGRAMS
21 ! FILE CREATED AUGUST 19, 1994
22 ! FROM MADRONICH 1989 DATA FILE
23 ! et_flux
24 ! EXTRA-TERRESTIAL IRRADIANCE
25 ! FILE CREATED AUGUST 19, 1994
26 ! FROM MADRONICH 1989 DATA FILE
27 ! -----------------------------------------------------------------------
28 ! current jindex assignments - if calculation order changes
29 ! these change - subroutine yield must be changed!
30 ! process jindex notes
31 ! wavelength none center values, modified wmo grid
32 ! e.t.irradian none photons per bin, not per nm
33 ! absorption cross sections:
34 ! o2 absorp 1 schumman-runge corrected in srband
35 ! o3 -> 1d 2 at 275 k. correct t-dep in subgrid
36 ! o3 -> 3p 3 at 275 k. correct t-dep in subgrid
37 ! no2 4
38 ! no3 -> no+o2 5
39 ! no3 -> no2+o 6
40 ! hno2 7
41 ! hno3 8
42 ! hno4 9
43 ! h2o2 10
44 ! ch2o -> rad 11
45 ! ch2o -> mol 12
46 ! ch3cho 13
47 ! ch3coch3 14
48 ! ch3coc2h5 15
49 ! hcocho proc a 16
50 ! ch3cocho 17
51 ! hcoch=chcho 18 estimate, no reliable measurement
52 ! ch3o2h 19
53 ! ch3coo2h 20 actually use 0.28*(h2o2 value)
54 ! ch3ono2 21
55 ! hcocho proc b 22
56 ! quantum yields:
57 ! no2 4
58 ! ch2o -> rad 11
59 ! ch2o -> mol 12
60 ! ch3cho 13
61 ! hcoch=chcho 18 (energetic threshold)
62 ! cross section and quantum yield data. this section assigns
63 ! yields for ntp air. yields are read from data file
64 ! some yields are temperature and/or pressure dependent:ch3cho, ch2o_b,
65 ! o3.
66 ! for ch2o_b and ch3cho, the data read above are stp values. these will
67 ! be
68 ! corrected in subgrid for t & ad dependence, after the altitude
69 ! dependent
70 ! values of t and ad are computed at each layer or level as appropriate.
71 ! the o3->o(1d) yield is also calculated later, in subgrid.
72 ! -----------------------------------------------------
73 ! o2 cross section
74 ! -----------------------------------------------------
75 ! o3 -> o1d cross section
76 ! -----------------------------------------------------
77 ! o3 -> o3p cross section
78 ! -----------------------------------------------------
79 ! no2 cross section
80 ! no2 quantum yield
81 ! -----------------------------------------------------
82 ! no3 -> no + o2 cross section
83 ! no3 -> no + o2 quantum yield
84 ! -----------------------------------------------------
85 ! no3 -> no2 + o cross section
86 ! no3 -> no2 + o quantum yield
87 ! -----------------------------------------------------
88 ! hono cross section
89 ! hono cross quantum yield
90 ! -----------------------------------------------------
91 ! hno3 cross section
92 ! hno3 cross quantum yield
93 ! HNO3 CROSS SECTION TEMPERATURE DEPENDENCE
94 ! -----------------------------------------------------
95 ! hno4 cross section
96 ! hno4 cross quantum yield
97 ! -----------------------------------------------------
98 ! h2o2 cross section
99 ! h2o2 cross quantum yield
100 ! -----------------------------------------------------
101 ! hcho -> ho2 cross section
102 ! hcho -> ho2 quantum yield
103 ! -----------------------------------------------------
104 ! hcho -> h2 cross section
105 ! hcho -> h2 quantum yield
106 ! -----------------------------------------------------
107 ! ch3cho cross section
108 ! ch3cho (ntp) quantum yield
109 ! -----------------------------------------------------
110 ! ch3coch3 cross section
111 ! ch3coch3 quantum yield
112 ! -----------------------------------------------------
113 ! ch3coc2h5 cross section
114 ! ch3coc2h5 quantum yield
115 ! -----------------------------------------------------
116 ! hcocho proc a cross section
117 ! hcocho a quantum yield
118 ! -----------------------------------------------------
119 ! ch3cocho cross section
120 ! ch3cocho quantum yield
121 ! -----------------------------------------------------
122 ! dcb cross section
123 ! dcb quantum yield
124 ! -----------------------------------------------------
125 ! ch3o2h cross section
126 ! ch3o2h cross quantum yield
127 ! -----------------------------------------------------
128 ! ch3coo2h cross section
129 ! ch3coo2h cross quantum yield
130 ! -----------------------------------------------------
131 ! ch3ono2 cross section
132 ! ch3ono2 cross quantum yield
133 ! -----------------------------------------------------
134 ! hcocho proc b cross section
135 ! hcocho b quantum yield
136 ! macr cross section
137 ! macr quantum yield
138 ! .. Parameters ..
139 INTEGER, PARAMETER :: kl0 = 30, kl1 = 130
140 INTEGER, PARAMETER :: kldif = (kl1-kl0+1)*3
141 INTEGER, PARAMETER :: nabv = 10, nj = 200, nreakj = 23
142 INTEGER, PARAMETER :: mj = 2*nj - 2
143 ! .. Local Scalars ..
144 INTEGER :: ip, kl
145 ! .. Local Arrays ..
146 REAL :: aerstd(51), airstd(51), albedoph(130), caabv(nabv), fext(130), &
147 o3abv(nabv), o3std(51), pabv(nabv), so3tx(70,3), sra(11,9), srb(11,5), &
148 tabv(nabv), tstd(51), txs(130,nreakj), wl(130), xqy(130,23), &
149 xs(130,nreakj), zabv(nabv), zstd(51)
150 ! .. Data Statements ..
151 DATA zabv/21., 22., 23., 24., 25., 30., 35., 40., 45., 50./
152 DATA tabv/215.19, 215.19, 215.19, 215.19, 215.19, 217.39, 227.80, &
153 243.19, 258.50, 265.70/
154 DATA pabv/1.57E18, 1.34E18, 1.14E18, 9.76E17, 8.33E17, 3.83E17, 1.76E17, &
155 8.31E16, 4.09E16, 2.14E16/
156 DATA o3abv/4.88E12, 4.86E12, 4.73E12, 4.54E12, 4.32E12, 2.52E12, &
157 1.40E12, 6.07E11, 2.03E11, 6.61E10/
158 DATA caabv/1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 1.90E-4, &
159 5.00E-5, 1.32E-5, 3.46E-6, 9.14E-7/
160 DATA zstd/0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., &
161 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., &
162 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., &
163 42., 43., 44., 45., 46., 47., 48., 49., 50./
164 DATA ((sra(kl,ip),ip=1,9),kl=1,11)/ -2.158311E+01, -4.164652E-01, &
165 5.266362E-02, 1.655877E-02, 0., 0., 0., 0., 0., -2.184813E+01, &
166 -4.753880E-01, 4.519945E-02, 3.228313E-02, 3.079373E-03, 0., 0., 0., &
167 0., -2.200507E+01, -4.628729E-01, -5.022541E-02, 2.545036E-02, &
168 5.791406E-02, 1.179966E-02, -8.296876E-03, -3.238368E-03, &
169 -3.069686E-04, -2.205527E+01, -4.400848E-01, -5.687308E-03, &
170 3.712279E-02, 6.025527E-03, 0., 0., 0., 0., -2.205261E+01, &
171 -5.707936E-01, -3.330207E-02, 5.959032E-02, 1.510540E-02, &
172 1.000376E-03, 0., 0., 0., -2.228000E+01, -3.960759E-01, -2.995798E-02, &
173 4.918104E-02, 9.269080E-03, -1.173411E-03, -2.599386E-04, 0., 0., &
174 -2.275796E+01, -2.054719E-01, -1.094205E-02, 2.079595E-02, &
175 3.769638E-03, 0., 0., 0., 0., -2.297610E+01, -5.823677E-02, &
176 -1.007612E-01, 2.404666E-02, 4.761876E-02, 4.169606E-03, &
177 -7.126663E-03, -2.263652E-03, -1.971653E-04, -2.506084E+01, &
178 3.442774E-02, -2.212047E-04, 6.186041E-07, -6.284394E-10, 0., 0., 0., &
179 0., -2.313436E+01, 1.177283E-04, 0., 0., 0., 0., 0., 0., 0., &
180 -2.312205E+01, 0., 0., 0., 0., 0., 0., 0., 0./
181 DATA ((srb(kl,ip),ip=1,5),kl=1,11)/ -2.431640E+03, 4.729722E+02, &
182 -3.452121E+01, 1.120677E+00, -1.365618E-02, -3.701955E+01, &
183 3.623290E+00, -8.929223E-02, 0., 0., -1.086239E+03, 1.981847E+02, &
184 -1.359057E+01, 4.155845E-01, -4.788462E-03, -1.213108E+03, &
185 2.277459E+02, -1.612207E+01, 5.101389E-01, -6.090518E-03, &
186 -8.334575E+01, 7.944254E+00, -1.898894E-01, 0., 0., -2.139117E+02, &
187 2.612729E+01, -1.036749E+00, 1.317695E-02, 0., -3.281301E+02, &
188 4.307004E+01, -1.870019E+00, 2.674331E-02, 0., 3.033416E+03, &
189 -5.978911E+02, 4.370384E+01, -1.406715E+00, 1.683967E-02, &
190 -2.535815E+00, 0., 0., 0., 0., -4.474937E+00, 0., 0., 0., 0., &
191 -2.996639E+00, 0., 0., 0., 0./
192 DATA ((so3tx(kl,ip),kl=33,61),ip=1,3)/9.630E+00, 8.320E+00, 6.880E+00, &
193 5.370E+00, 3.960E+00, 2.710E+00, 1.750E+00, 1.060E+00, 5.960E-01, &
194 3.330E-01, 2.400E-01, 2.100E-01, 1.800E-01, 1.600E-01, 1.400E-01, &
195 1.200E-01, 1.050E-01, 9.000E-02, 8.000E-02, 7.000E-02, 6.000E-02, &
196 5.500E-02, 4.000E-02, 2.190E-02, 1.010E-02, 5.080E-03, 2.120E-03, &
197 8.290E-04, 2.940E-04, 1.190E-03, 3.640E-04, 2.460E-04, 1.030E-03, &
198 1.690E-03, 1.450E-03, 8.940E-04, 7.830E-04, 4.940E-04, 3.550E-04, &
199 2.950E-04, 2.750E-04, 2.500E-04, 2.300E-04, 2.080E-04, 1.860E-04, &
200 1.640E-04, 1.450E-04, 1.280E-04, 1.121E-04, 1.000E-04, 9.200E-05, &
201 7.500E-05, 4.830E-05, 3.430E-05, 1.820E-05, 8.850E-06, 4.270E-06, &
202 5.300E-06, -1.740E-05, 2.470E-06, 1.170E-05, 1.260E-06, -6.860E-06, &
203 -2.890E-06, 3.590E-06, 2.000E-06, 3.660E-06, 2.600E-06, 2.170E-06, &
204 1.950E-06, 1.380E-06, 1.650E-06, 1.550E-06, 1.460E-06, 1.340E-06, &
205 1.210E-06, 1.130E-06, 1.060E-06, 9.400E-07, 8.700E-07, 7.500E-07, &
206 5.200E-07, 2.660E-07, 1.630E-07, 1.260E-07, 8.710E-08, 3.500E-08/
207 DATA aerstd/2.40E-1, 1.06E-1, 4.56E-2, 1.91E-2, 1.01E-2, 7.63E-3, &
208 5.38E-3, 5.00E-3, 5.15E-3, 4.94E-3, 4.82E-3, 4.51E-3, 4.74E-3, &
209 4.37E-3, 4.28E-3, 4.03E-3, 3.83E-3, 3.78E-3, 3.88E-3, 3.08E-3, &
210 2.26E-3, 1.64E-3, 1.23E-3, 9.45E-4, 7.49E-4, 6.30E-4, 5.50E-4, &
211 4.21E-4, 3.22E-4, 2.48E-4, 1.90E-4, 1.45E-4, 1.11E-4, 8.51E-5, &
212 6.52E-5, 5.00E-5, 3.83E-5, 2.93E-5, 2.25E-5, 1.72E-5, 1.32E-5, &
213 1.01E-5, 7.72E-6, 5.91E-6, 4.53E-6, 3.46E-6, 2.66E-6, 2.04E-6, &
214 1.56E-6, 1.19E-6, 9.14E-7/
215 DATA (wl(kl),kl=1,130)/1.861E+02, 1.878E+02, 1.896E+02, 1.914E+02, &
216 1.933E+02, 1.952E+02, 1.971E+02, 1.990E+02, 2.010E+02, 2.031E+02, &
217 2.052E+02, 2.073E+02, 2.094E+02, 2.117E+02, 2.139E+02, 2.162E+02, &
218 2.186E+02, 2.210E+02, 2.235E+02, 2.260E+02, 2.286E+02, 2.313E+02, &
219 2.340E+02, 2.367E+02, 2.396E+02, 2.425E+02, 2.454E+02, 2.485E+02, &
220 2.516E+02, 2.548E+02, 2.582E+02, 2.615E+02, 2.650E+02, 2.685E+02, &
221 2.722E+02, 2.759E+02, 2.798E+02, 2.837E+02, 2.878E+02, 2.920E+02, &
222 2.963E+02, 3.005E+02, 3.030E+02, 3.040E+02, 3.050E+02, 3.060E+02, &
223 3.070E+02, 3.080E+02, 3.090E+02, 3.100E+02, 3.110E+02, 3.120E+02, &
224 3.130E+02, 3.140E+02, 3.160E+02, 3.200E+02, 3.250E+02, 3.300E+02, &
225 3.350E+02, 3.400E+02, 3.450E+02, 3.500E+02, 3.550E+02, 3.600E+02, &
226 3.650E+02, 3.700E+02, 3.750E+02, 3.800E+02, 3.850E+02, 3.900E+02, &
227 3.950E+02, 4.000E+02, 4.050E+02, 4.100E+02, 4.150E+02, 4.200E+02, &
228 4.250E+02, 4.300E+02, 4.350E+02, 4.400E+02, 4.450E+02, 4.500E+02, &
229 4.550E+02, 4.600E+02, 4.650E+02, 4.700E+02, 4.750E+02, 4.800E+02, &
230 4.850E+02, 4.900E+02, 4.950E+02, 5.000E+02, 5.050E+02, 5.100E+02, &
231 5.150E+02, 5.200E+02, 5.250E+02, 5.300E+02, 5.350E+02, 5.400E+02, &
232 5.450E+02, 5.500E+02, 5.550E+02, 5.600E+02, 5.650E+02, 5.700E+02, &
233 5.750E+02, 5.800E+02, 5.850E+02, 5.900E+02, 5.950E+02, 6.000E+02, &
234 6.050E+02, 6.100E+02, 6.150E+02, 6.200E+02, 6.250E+02, 6.300E+02, &
235 6.350E+02, 6.400E+02, 6.448E+02, 6.510E+02, 6.600E+02, 6.700E+02, &
236 6.800E+02, 6.900E+02, 7.000E+02, 7.100E+02, 7.200E+02, 7.300E+02/
237 DATA (fext(kl),kl=1,130)/3.620E+11, 4.730E+11, 5.610E+11, 6.630E+11, &
238 6.900E+11, 9.560E+11, 1.150E+12, 1.270E+12, 1.520E+12, 1.780E+12, &
239 2.200E+12, 2.690E+12, 4.540E+12, 7.140E+12, 8.350E+12, 8.390E+12, &
240 1.080E+13, 1.180E+13, 1.600E+13, 1.340E+13, 1.410E+13, 1.570E+13, &
241 1.380E+13, 1.600E+13, 1.450E+13, 2.200E+13, 1.990E+13, 1.970E+13, &
242 1.940E+13, 2.910E+13, 4.950E+13, 4.530E+13, 1.070E+14, 1.200E+14, &
243 1.100E+14, 1.040E+14, 8.240E+13, 1.520E+14, 2.150E+14, 3.480E+14, &
244 3.396E+14, 2.730E+14, 9.109E+13, 8.745E+13, 9.577E+13, 8.507E+13, &
245 9.383E+13, 1.030E+14, 9.722E+13, 7.751E+13, 1.277E+14, 1.087E+14, &
246 1.102E+14, 1.184E+14, 3.153E+14, 5.930E+14, 6.950E+14, 8.150E+14, &
247 7.810E+14, 8.350E+14, 8.140E+14, 8.530E+14, 9.170E+14, 8.380E+14, &
248 1.040E+15, 1.100E+15, 9.790E+14, 1.130E+15, 8.890E+14, 1.140E+15, &
249 9.170E+14, 1.690E+15, 1.700E+15, 1.840E+15, 1.870E+15, 1.950E+15, &
250 1.810E+15, 1.670E+15, 1.980E+15, 2.020E+15, 2.180E+15, 2.360E+15, &
251 2.310E+15, 2.390E+15, 2.380E+15, 2.390E+15, 2.440E+15, 2.510E+15, &
252 2.300E+15, 2.390E+15, 2.480E+15, 2.400E+15, 2.460E+15, 2.490E+15, &
253 2.320E+15, 2.390E+15, 2.420E+15, 2.550E+15, 2.510E+15, 2.490E+15, &
254 2.550E+15, 2.530E+15, 2.540E+15, 2.500E+15, 2.570E+15, 2.580E+15, &
255 2.670E+15, 2.670E+15, 2.700E+15, 2.620E+15, 2.690E+15, 2.630E+15, &
256 2.680E+15, 2.660E+15, 2.590E+15, 2.690E+15, 2.610E+15, 2.620E+15, &
257 2.620E+15, 2.630E+15, 2.392E+15, 3.998E+15, 5.115E+15, 5.225E+15, &
258 5.215E+15, 5.105E+15, 5.140E+15, 5.010E+15, 4.930E+15, 4.895E+15/
259 DATA (xs(kl,1),kl=1,130)/7.040E-24, 7.360E-24, 7.640E-24, 7.870E-24, &
260 8.040E-24, 8.140E-24, 8.170E-24, 8.130E-24, 8.010E-24, 7.840E-24, &
261 7.630E-24, 7.330E-24, 6.990E-24, 6.450E-24, 5.810E-24, 5.230E-24, &
262 4.710E-24, 4.260E-24, 3.800E-24, 3.350E-24, 2.900E-24, 2.450E-24, &
263 2.050E-24, 1.690E-24, 1.300E-24, 9.300E-25, 0., 0., 0., 0., 0., 0., &
264 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
265 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
266 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
267 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
268 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., &
269 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0./
270 DATA (xs(kl,2),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
271 0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
272 0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
273 0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
274 0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
275 0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
276 0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
277 0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
278 0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
279 0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
280 0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
281 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
282 0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
283 0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
284 0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
285 0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
286 0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
287 0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
288 0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
289 0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
290 0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
291 0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
292 DATA (xs(kl,3),kl=1,130)/0.622E-18, 0.576E-18, 0.526E-18, 0.476E-18, &
293 0.428E-18, 0.383E-18, 0.347E-18, 0.323E-18, 0.314E-18, 0.326E-18, &
294 0.364E-18, 0.434E-18, 0.542E-18, 0.699E-18, 0.921E-18, 0.119E-17, &
295 0.155E-17, 0.199E-17, 0.256E-17, 0.323E-17, 0.400E-17, 0.483E-17, &
296 0.579E-17, 0.686E-17, 0.797E-17, 0.900E-17, 0.100E-16, 0.108E-16, &
297 0.113E-16, 0.115E-16, 0.112E-16, 0.106E-16, 0.963E-17, 0.836E-17, &
298 0.695E-17, 0.545E-17, 0.404E-17, 0.280E-17, 0.183E-17, 0.112E-17, &
299 0.647E-18, 0.369E-18, 0.270E-18, 0.238E-18, 0.203E-18, 0.183E-18, &
300 0.161E-18, 0.139E-18, 0.122E-18, 0.105E-18, 0.939E-19, 0.825E-19, &
301 0.711E-19, 0.653E-19, 0.486E-19, 0.276E-19, 0.137E-19, 0.707E-20, &
302 0.330E-20, 0.152E-20, 0.816E-21, 0.266E-21, 0.109E-21, 0.549E-22, &
303 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
304 0.000E+00, 0.000E+00, 0.000E+00, 0.291E-22, 0.314E-22, 0.399E-22, &
305 0.654E-22, 0.683E-22, 0.866E-22, 0.125E-21, 0.149E-21, 0.171E-21, &
306 0.212E-21, 0.357E-21, 0.368E-21, 0.406E-21, 0.489E-21, 0.711E-21, &
307 0.843E-21, 0.828E-21, 0.909E-21, 0.122E-20, 0.162E-20, 0.158E-20, &
308 0.160E-20, 0.178E-20, 0.207E-20, 0.255E-20, 0.274E-20, 0.288E-20, &
309 0.307E-20, 0.317E-20, 0.336E-20, 0.388E-20, 0.431E-20, 0.467E-20, &
310 0.475E-20, 0.455E-20, 0.435E-20, 0.442E-20, 0.461E-20, 0.489E-20, &
311 0.484E-20, 0.454E-20, 0.424E-20, 0.390E-20, 0.360E-20, 0.343E-20, &
312 0.317E-20, 0.274E-20, 0.261E-20, 0.240E-20, 0.207E-20, 0.172E-20, &
313 0.137E-20, 0.111E-20, 0.913E-21, 0.793E-21, 0.640E-21, 0.514E-21/
314 DATA (xs(kl,4),kl=1,130)/0.259E-18, 0.272E-18, 0.286E-18, 0.273E-18, &
315 0.251E-18, 0.244E-18, 0.246E-18, 0.246E-18, 0.282E-18, 0.415E-18, &
316 0.448E-18, 0.445E-18, 0.464E-18, 0.487E-18, 0.482E-18, 0.502E-18, &
317 0.444E-18, 0.471E-18, 0.377E-18, 0.393E-18, 0.274E-18, 0.278E-18, &
318 0.169E-18, 0.162E-18, 0.882E-19, 0.747E-19, 0.391E-19, 0.275E-19, &
319 0.201E-19, 0.197E-19, 0.211E-19, 0.236E-19, 0.270E-19, 0.325E-19, &
320 0.379E-19, 0.503E-19, 0.588E-19, 0.700E-19, 0.815E-19, 0.972E-19, &
321 0.115E-18, 0.128E-18, 0.154E-18, 0.159E-18, 0.158E-18, 0.156E-18, &
322 0.164E-18, 0.166E-18, 0.182E-18, 0.184E-18, 0.192E-18, 0.204E-18, &
323 0.204E-18, 0.202E-18, 0.224E-18, 0.248E-18, 0.281E-18, 0.313E-18, &
324 0.343E-18, 0.380E-18, 0.407E-18, 0.431E-18, 0.472E-18, 0.483E-18, &
325 0.517E-18, 0.532E-18, 0.551E-18, 0.564E-18, 0.576E-18, 0.593E-18, &
326 0.585E-18, 0.602E-18, 0.578E-18, 0.600E-18, 0.565E-18, 0.581E-18, &
327 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
328 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
329 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
330 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
331 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
332 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
333 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
334 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
335 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
336 DATA ((xqy(kl,ip),kl=kl0,kl1),ip=1,3)/kldif*1./
337 DATA (xqy(kl,4),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
338 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
339 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
340 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
341 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
342 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
343 0.999000, 0.998000, 0.997000, 0.996500, 0.996000, 0.996000, 0.996000, &
344 0.996000, 0.995000, 0.995000, 0.995000, 0.995000, 0.995000, 0.994000, &
345 0.994000, 0.994000, 0.993000, 0.992000, 0.991000, 0.990000, 0.989000, &
346 0.988000, 0.987000, 0.986000, 0.984000, 0.983000, 0.981000, 0.979000, &
347 0.975000, 0.969000, 0.960000, 0.927000, 0.694000, 0.355000, 0.134000, &
348 0.060000, 0.018000, 0.000900, 0.000000, 0.000000, 0.000000, 0.000000, &
349 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
350 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
351 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
352 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
353 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
354 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
355 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
356 DATA (xs(kl,5),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
357 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
358 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
359 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
360 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
361 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
362 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
363 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
364 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
365 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
366 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
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.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
369 0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
370 0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
371 0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
372 0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
373 0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
374 0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
375 0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
376 0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
377 0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
378 DATA (xqy(kl,5),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
379 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
380 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
381 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
382 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
383 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
384 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
385 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
386 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
387 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
388 0.000000, 0.000000, 0.000000, 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.038000, &
394 0.191000, 0.326000, 0.311000, 0.272000, 0.233000, 0.194000, 0.156000, &
395 0.117000, 0.078000, 0.039000, 0.000000, 0.000000, 0.000000, 0.000000, &
396 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
397 DATA (xs(kl,6),kl=1,130)/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.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
409 0.000E+00, 0.000E+00, 0.160E-19, 0.240E-19, 0.520E-19, 0.760E-19, &
410 0.110E-18, 0.136E-18, 0.170E-18, 0.198E-18, 0.220E-18, 0.288E-18, &
411 0.358E-18, 0.396E-18, 0.502E-18, 0.598E-18, 0.694E-18, 0.716E-18, &
412 0.828E-18, 0.984E-18, 0.110E-17, 0.114E-17, 0.125E-17, 0.153E-17, &
413 0.156E-17, 0.168E-17, 0.169E-17, 0.217E-17, 0.229E-17, 0.208E-17, &
414 0.213E-17, 0.261E-17, 0.299E-17, 0.329E-17, 0.278E-17, 0.281E-17, &
415 0.307E-17, 0.334E-17, 0.322E-17, 0.554E-17, 0.441E-17, 0.314E-17, &
416 0.365E-17, 0.188E-17, 0.233E-17, 0.473E-17, 0.100E-16, 0.584E-17, &
417 0.180E-17, 0.135E-17, 0.822E-18, 0.640E-18, 0.777E-17, 0.134E-17, &
418 0.337E-18, 0.175E-19, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
419 DATA (xqy(kl,6),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
420 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
421 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
422 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
423 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
424 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
425 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
426 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
427 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
428 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
429 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
430 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
431 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
432 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
433 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
434 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.962000, &
435 0.809000, 0.661000, 0.578000, 0.506000, 0.433000, 0.361000, 0.289000, &
436 0.217000, 0.144000, 0.072000, 0.000000, 0.000000, 0.000000, 0.000000, &
437 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
438 DATA (xs(kl,7),kl=1,130)/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.130E-19, 0.190E-19, 0.280E-19, &
447 0.220E-19, 0.360E-19, 0.340E-19, 0.536E-19, 0.534E-19, 0.111E-18, &
448 0.786E-19, 0.189E-18, 0.116E-18, 0.130E-18, 0.279E-18, 0.954E-19, &
449 0.179E-18, 0.260E-18, 0.590E-19, 0.101E-18, 0.176E-18, 0.304E-19, &
450 0.775E-20, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
451 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
452 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
453 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
454 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
455 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
456 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
457 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
458 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
459 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
460 DATA (xqy(kl,7),kl=1,130)/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, 1.000000, &
475 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
476 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
477 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
478 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
479 DATA (xs(kl,8),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.127E-16, &
480 0.114E-16, 0.100E-16, 0.847E-17, 0.679E-17, 0.518E-17, 0.382E-17, &
481 0.270E-17, 0.182E-17, 0.120E-17, 0.730E-18, 0.451E-18, 0.283E-18, &
482 0.195E-18, 0.134E-18, 0.102E-18, 0.802E-19, 0.650E-19, 0.518E-19, &
483 0.414E-19, 0.321E-19, 0.265E-19, 0.230E-19, 0.209E-19, 0.199E-19, &
484 0.196E-19, 0.195E-19, 0.193E-19, 0.188E-19, 0.180E-19, 0.170E-19, &
485 0.152E-19, 0.134E-19, 0.113E-19, 0.924E-20, 0.719E-20, 0.532E-20, &
486 0.371E-20, 0.249E-20, 0.188E-20, 0.167E-20, 0.150E-20, 0.133E-20, &
487 0.119E-20, 0.105E-20, 0.932E-21, 0.814E-21, 0.721E-21, 0.628E-21, &
488 0.547E-21, 0.465E-21, 0.362E-21, 0.197E-21, 0.975E-22, 0.452E-22, &
489 0.222E-22, 0.110E-22, 0.604E-23, 0.420E-23, 0.000E+00, 0.000E+00, &
490 0.000E+00, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
501 DATA (xqy(kl,8),kl=1,130)/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 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
520 DATA (txs(kl,8),kl=1,130)/0.000000, 0.000000, 0.000000, 1.700000, &
521 1.650000, 1.660000, 1.690000, 1.740000, 1.770000, 1.850000, 1.970000, &
522 2.080000, 2.170000, 2.170000, 2.210000, 2.150000, 2.060000, 1.960000, &
523 1.840000, 1.780000, 1.800000, 1.860000, 1.900000, 1.970000, 1.970000, &
524 1.970000, 1.880000, 1.750000, 1.610000, 1.440000, 1.340000, 1.230000, &
525 1.180000, 1.140000, 1.120000, 1.140000, 1.140000, 1.180000, 1.220000, &
526 1.250000, 1.450000, 1.490000, 1.560000, 1.640000, 1.690000, 1.780000, &
527 1.870000, 1.940000, 2.040000, 2.150000, 2.270000, 2.380000, 2.620000, &
528 2.700000, 2.920000, 3.100000, 3.240000, 3.520000, 3.770000, 3.910000, &
529 4.230000, 4.700000, 5.150000, 5.250000, 5.740000, 6.450000, 6.700000, &
530 7.160000, 7.550000, 8.160000, 9.750000, 9.930000, 9.600000, 10.50000, &
531 10.80000, 11.80000, 11.80000, 9.300000, 12.10000, 11.90000, 9.300000, &
532 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
533 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
534 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
535 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
536 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
537 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
538 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
539 DATA (xs(kl,9),kl=1,130)/0.000E+00, 0.000E+00, 0.100E-16, 0.980E-17, &
540 0.918E-17, 0.828E-17, 0.723E-17, 0.618E-17, 0.517E-17, 0.426E-17, &
541 0.352E-17, 0.292E-17, 0.245E-17, 0.205E-17, 0.175E-17, 0.151E-17, &
542 0.131E-17, 0.115E-17, 0.102E-17, 0.916E-18, 0.827E-18, 0.752E-18, &
543 0.687E-18, 0.631E-18, 0.578E-18, 0.529E-18, 0.484E-18, 0.439E-18, &
544 0.396E-18, 0.353E-18, 0.311E-18, 0.271E-18, 0.231E-18, 0.194E-18, &
545 0.158E-18, 0.125E-18, 0.946E-19, 0.694E-19, 0.485E-19, 0.325E-19, &
546 0.210E-19, 0.135E-19, 0.104E-19, 0.938E-20, 0.846E-20, 0.763E-20, &
547 0.689E-20, 0.623E-20, 0.564E-20, 0.511E-20, 0.463E-20, 0.420E-20, &
548 0.382E-20, 0.347E-20, 0.287E-20, 0.191E-20, 0.101E-20, 0.000E+00, &
549 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
550 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
551 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
552 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
553 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
554 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
555 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
556 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
557 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
558 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
559 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
560 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
561 DATA (xqy(kl,9),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
562 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
563 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
564 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
565 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
566 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
567 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
568 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
569 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
570 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
571 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
572 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
573 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
574 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
575 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
576 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
577 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
578 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
579 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
580 DATA (xs(kl,10),kl=1,130)/0.000E+00, 0.000E+00, 0.695E-18, 0.635E-18, &
581 0.586E-18, 0.547E-18, 0.515E-18, 0.488E-18, 0.462E-18, 0.436E-18, &
582 0.412E-18, 0.388E-18, 0.365E-18, 0.341E-18, 0.318E-18, 0.295E-18, &
583 0.272E-18, 0.250E-18, 0.229E-18, 0.209E-18, 0.190E-18, 0.172E-18, &
584 0.155E-18, 0.140E-18, 0.126E-18, 0.112E-18, 0.999E-19, 0.881E-19, &
585 0.774E-19, 0.675E-19, 0.582E-19, 0.500E-19, 0.425E-19, 0.358E-19, &
586 0.298E-19, 0.247E-19, 0.201E-19, 0.164E-19, 0.131E-19, 0.105E-19, &
587 0.828E-20, 0.658E-20, 0.573E-20, 0.543E-20, 0.514E-20, 0.486E-20, &
588 0.460E-20, 0.435E-20, 0.412E-20, 0.390E-20, 0.369E-20, 0.349E-20, &
589 0.330E-20, 0.312E-20, 0.279E-20, 0.220E-20, 0.160E-20, 0.130E-20, &
590 0.100E-20, 0.700E-21, 0.500E-21, 0.400E-21, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
602 DATA (xqy(kl,10),kl=1,130)/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 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
621 DATA (xs(kl,11),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
622 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
623 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
624 0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
625 0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
626 0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
627 0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
628 0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
629 0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
630 0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
631 0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
643 DATA (xqy(kl,11),kl=1,130)/0.000000, 0.000000, 0.000000, 0.000000, &
644 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
645 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
646 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.264091, &
647 0.288940, 0.297038, 0.295059, 0.289384, 0.285483, 0.288087, 0.301000, &
648 0.326819, 0.366764, 0.420506, 0.485961, 0.559106, 0.633887, 0.702103, &
649 0.733457, 0.740762, 0.747845, 0.753000, 0.754000, 0.754800, 0.754000, &
650 0.753000, 0.752000, 0.751000, 0.749500, 0.745000, 0.739600, 0.731700, &
651 0.723300, 0.690300, 0.593100, 0.458100, 0.305000, 0.122300, 0.003429, &
652 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
653 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
654 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
655 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
656 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
657 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
658 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
659 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
660 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
661 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
662 DATA (xs(kl,12),kl=1,130)/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.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
665 0.293E-21, 0.342E-21, 0.102E-20, 0.456E-21, 0.527E-21, 0.537E-21, &
666 0.347E-21, 0.759E-21, 0.628E-21, 0.974E-21, 0.104E-20, 0.219E-20, &
667 0.228E-20, 0.357E-20, 0.374E-20, 0.584E-20, 0.651E-20, 0.102E-19, &
668 0.114E-19, 0.176E-19, 0.180E-19, 0.259E-19, 0.227E-19, 0.275E-19, &
669 0.318E-19, 0.160E-19, 0.245E-19, 0.637E-19, 0.426E-19, 0.399E-19, &
670 0.186E-19, 0.131E-19, 0.310E-19, 0.182E-19, 0.596E-20, 0.111E-19, &
671 0.911E-20, 0.457E-19, 0.423E-19, 0.142E-19, 0.243E-19, 0.178E-19, &
672 0.129E-20, 0.213E-19, 0.661E-20, 0.139E-20, 0.827E-20, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
684 DATA (xqy(kl,12),kl=1,130)/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.000000, &
687 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.492085, &
688 0.483681, 0.483325, 0.487471, 0.492514, 0.495532, 0.493893, 0.485473, &
689 0.468839, 0.443373, 0.409405, 0.368400, 0.323132, 0.307820, 0.294564, &
690 0.280920, 0.266885, 0.253277, 0.249000, 0.247000, 0.245600, 0.248000, &
691 0.251000, 0.254000, 0.257000, 0.260200, 0.264500, 0.269000, 0.273500, &
692 0.278900, 0.310300, 0.394100, 0.508100, 0.676100, 0.759300, 0.636100, &
693 0.501500, 0.373400, 0.229000, 0.103600, 0.005906, 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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
703 DATA (xs(kl,13),kl=1,130)/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.593E-21, 0.548E-21, &
705 0.552E-21, 0.513E-21, 0.480E-21, 0.482E-21, 0.482E-21, 0.482E-21, &
706 0.536E-21, 0.593E-21, 0.734E-21, 0.948E-21, 0.125E-20, 0.171E-20, &
707 0.234E-20, 0.321E-20, 0.434E-20, 0.582E-20, 0.770E-20, 0.991E-20, &
708 0.127E-19, 0.159E-19, 0.200E-19, 0.237E-19, 0.287E-19, 0.326E-19, &
709 0.376E-19, 0.408E-19, 0.444E-19, 0.463E-19, 0.466E-19, 0.465E-19, &
710 0.432E-19, 0.406E-19, 0.372E-19, 0.348E-19, 0.342E-19, 0.342E-19, &
711 0.336E-19, 0.333E-19, 0.314E-19, 0.293E-19, 0.276E-19, 0.253E-19, &
712 0.247E-19, 0.243E-19, 0.210E-19, 0.169E-19, 0.108E-19, 0.651E-20, &
713 0.314E-20, 0.138E-20, 0.224E-21, 0.947E-22, 0.425E-22, 0.229E-22, &
714 0.275E-23, 0.192E-23, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
725 DATA (xqy(kl,13),kl=1,130)/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.000000, &
728 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
729 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.310000, &
730 0.349300, 0.377900, 0.430300, 0.501600, 0.566200, 0.561500, 0.541100, &
731 0.512100, 0.473200, 0.430000, 0.392000, 0.376000, 0.360000, 0.344000, &
732 0.328000, 0.312000, 0.296000, 0.279800, 0.262500, 0.245000, 0.227500, &
733 0.210000, 0.175000, 0.109300, 0.051880, 0.006266, 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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
744 DATA (xs(kl,14),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
745 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.003E-20, 1.873E-21, &
746 1.408E-21, 1.084E-21, 1.030E-21, 1.071E-21, 1.176E-21, 1.371E-21, &
747 1.690E-21, 2.153E-21, 2.779E-21, 3.561E-21, 4.575E-21, 5.906E-21, &
748 7.585E-21, 9.622E-21, 1.212E-20, 1.516E-20, 1.863E-20, 2.252E-20, &
749 2.676E-20, 3.143E-20, 3.616E-20, 4.058E-20, 4.475E-20, 4.774E-20, &
750 5.029E-20, 5.065E-20, 5.049E-20, 4.790E-20, 4.415E-20, 3.940E-20, &
751 3.308E-20, 2.717E-20, 2.371E-20, 2.244E-20, 2.106E-20, 1.953E-20, &
752 1.801E-20, 1.663E-20, 1.538E-20, 1.408E-20, 1.277E-20, 1.173E-20, &
753 1.081E-20, 9.675E-21, 7.783E-21, 4.712E-21, 2.078E-21, 7.065E-22, &
754 1.973E-22, 5.745E-23, 2.137E-23, 8.235E-24, 5.686E-24, 2.157E-24, &
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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
766 DATA (xqy(kl,14),kl=1,130)/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.000000, &
770 0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
771 0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
772 0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
773 0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
774 0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
785 DATA (xs(kl,15),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
786 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.234E-19, 0.685E-20, &
787 0.215E-20, 0.168E-20, 0.159E-20, 0.164E-20, 0.181E-20, 0.203E-20, &
788 0.230E-20, 0.268E-20, 0.320E-20, 0.387E-20, 0.471E-20, 0.582E-20, &
789 0.728E-20, 0.913E-20, 0.115E-19, 0.145E-19, 0.180E-19, 0.222E-19, &
790 0.267E-19, 0.321E-19, 0.374E-19, 0.430E-19, 0.479E-19, 0.525E-19, &
791 0.555E-19, 0.576E-19, 0.575E-19, 0.555E-19, 0.518E-19, 0.460E-19, &
792 0.388E-19, 0.319E-19, 0.269E-19, 0.251E-19, 0.233E-19, 0.217E-19, &
793 0.202E-19, 0.188E-19, 0.173E-19, 0.158E-19, 0.142E-19, 0.128E-19, &
794 0.114E-19, 0.101E-19, 0.796E-20, 0.463E-20, 0.196E-20, 0.705E-21, &
795 0.207E-21, 0.545E-22, 0.145E-22, 0.431E-23, 0.471E-23, 0.157E-23, &
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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
807 DATA (xqy(kl,15),kl=1,130)/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.000000, 0.000000, 0.000000, 0.000000, &
811 0.000000, 0.000000, 0.000000, 0.766400, 0.779200, 0.792800, 0.776000, &
812 0.720000, 0.672000, 0.620200, 0.586900, 0.551800, 0.457500, 0.355000, &
813 0.270000, 0.205500, 0.145000, 0.120000, 0.110000, 0.100000, 0.090000, &
814 0.080000, 0.070000, 0.060000, 0.050000, 0.047800, 0.045600, 0.043400, &
815 0.041200, 0.036800, 0.028000, 0.030500, 0.033000, 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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
826 DATA (xs(kl,16),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
827 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
828 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
829 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
830 0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
831 0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
832 0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
833 0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
834 0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
835 0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
836 0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
837 0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
838 0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
839 0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
840 0.201E-18, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
848 DATA (xqy(kl,16),kl=1,130)/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.000000, 0.000000, 0.000000, 0.000000, &
852 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
853 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
854 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
855 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
856 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
857 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
858 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
859 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
860 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
861 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
862 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
863 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
864 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
865 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, &
866 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000, 0.029000/
867 DATA (xs(kl,17),kl=1,130)/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.000E+00, 0.000E+00, 0.000E+00, 0.131E-19, 0.142E-19, 0.156E-19, &
871 0.174E-19, 0.189E-19, 0.205E-19, 0.219E-19, 0.233E-19, 0.252E-19, &
872 0.269E-19, 0.285E-19, 0.313E-19, 0.338E-19, 0.362E-19, 0.394E-19, &
873 0.427E-19, 0.450E-19, 0.486E-19, 0.476E-19, 0.479E-19, 0.465E-19, &
874 0.420E-19, 0.371E-19, 0.352E-19, 0.344E-19, 0.336E-19, 0.316E-19, &
875 0.296E-19, 0.276E-19, 0.256E-19, 0.237E-19, 0.227E-19, 0.218E-19, &
876 0.208E-19, 0.199E-19, 0.182E-19, 0.151E-19, 0.938E-20, 0.652E-20, &
877 0.482E-20, 0.323E-20, 0.300E-20, 0.394E-20, 0.560E-20, 0.695E-20, &
878 0.108E-19, 0.148E-19, 0.191E-19, 0.243E-19, 0.322E-19, 0.403E-19, &
879 0.473E-19, 0.566E-19, 0.692E-19, 0.846E-19, 0.968E-19, 0.103E-18, &
880 0.102E-18, 0.101E-18, 0.106E-18, 0.104E-18, 0.994E-19, 0.813E-19, &
881 0.395E-19, 0.109E-19, 0.327E-20, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
889 DATA (xqy(kl,17),kl=1,130)/1.000000, 1.000000, 1.000000, 1.000000, &
890 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
891 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
892 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
893 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
894 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
895 1.000000, 1.000000, 1.000000, 1.000000, 0.990000, 0.990000, 0.990000, &
896 0.980000, 0.980000, 0.980000, 0.980000, 0.970000, 0.970000, 0.970000, &
897 0.960000, 0.960000, 0.940000, 0.920000, 0.880000, 0.825000, 0.750000, &
898 0.660000, 0.560000, 0.480000, 0.400000, 0.320000, 0.250000, 0.200000, &
899 0.150000, 0.120000, 0.100000, 0.080000, 0.060000, 0.050000, 0.040000, &
900 0.030000, 0.020000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, &
901 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.005000, 0.000000, &
902 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
903 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
904 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
905 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
906 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
907 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
908 DATA (xs(kl,18),kl=1,130)/0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
909 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
910 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
911 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
912 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
913 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
914 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
915 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
916 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
917 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
918 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, 0.790E-19, &
919 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
920 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
921 0.000E+00, 0.000E+00, 0.000E+00, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
930 DATA (xqy(kl,18),kl=1,130)/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, 1.000000, 1.000000, 1.000000, &
936 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
937 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
938 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
939 1.000000, 1.000000, 0.500000, 0.000000, 0.000000, 0.000000, 0.000000, &
940 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
941 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
949 DATA (xs(kl,19),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
950 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
951 0.000E+00, 0.000E+00, 0.320E-18, 0.268E-18, 0.226E-18, 0.193E-18, &
952 0.167E-18, 0.147E-18, 0.129E-18, 0.115E-18, 0.102E-18, 0.899E-19, &
953 0.797E-19, 0.708E-19, 0.623E-19, 0.548E-19, 0.483E-19, 0.422E-19, &
954 0.369E-19, 0.321E-19, 0.278E-19, 0.242E-19, 0.209E-19, 0.180E-19, &
955 0.154E-19, 0.131E-19, 0.111E-19, 0.925E-20, 0.763E-20, 0.622E-20, &
956 0.501E-20, 0.402E-20, 0.352E-20, 0.333E-20, 0.316E-20, 0.299E-20, &
957 0.283E-20, 0.268E-20, 0.254E-20, 0.240E-20, 0.227E-20, 0.215E-20, &
958 0.204E-20, 0.193E-20, 0.172E-20, 0.138E-20, 0.105E-20, 0.801E-21, &
959 0.612E-21, 0.467E-21, 0.356E-21, 0.270E-21, 0.206E-21, 0.160E-21, &
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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
971 DATA (xqy(kl,19),kl=1,130)/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, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
980 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
981 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
982 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
983 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
984 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
985 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
986 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
987 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
988 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, &
989 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
990 DATA (xs(kl,20),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.898E-19, &
991 0.841E-19, 0.784E-19, 0.148E-18, 0.138E-18, 0.129E-18, 0.122E-18, &
992 0.114E-18, 0.107E-18, 0.998E-19, 0.932E-19, 0.868E-19, 0.807E-19, &
993 0.748E-19, 0.690E-19, 0.633E-19, 0.579E-19, 0.529E-19, 0.480E-19, &
994 0.433E-19, 0.390E-19, 0.349E-19, 0.313E-19, 0.278E-19, 0.247E-19, &
995 0.217E-19, 0.190E-19, 0.162E-19, 0.138E-19, 0.118E-19, 0.991E-20, &
996 0.829E-20, 0.690E-20, 0.566E-20, 0.452E-20, 0.361E-20, 0.288E-20, &
997 0.229E-20, 0.181E-20, 0.157E-20, 0.147E-20, 0.138E-20, 0.131E-20, &
998 0.125E-20, 0.118E-20, 0.112E-20, 0.106E-20, 0.100E-20, 0.947E-21, &
999 0.893E-21, 0.838E-21, 0.743E-21, 0.581E-21, 0.428E-21, 0.332E-21, &
1000 0.262E-21, 0.192E-21, 0.141E-21, 0.910E-22, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1012 DATA (xqy(kl,20),kl=1,130)/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 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1031 DATA (xs(kl,21),kl=1,130)/0.180E-16, 0.181E-16, 0.179E-16, 0.174E-16, &
1032 0.167E-16, 0.159E-16, 0.146E-16, 0.133E-16, 0.118E-16, 0.101E-16, &
1033 0.850E-17, 0.695E-17, 0.540E-17, 0.411E-17, 0.301E-17, 0.215E-17, &
1034 0.163E-17, 0.105E-17, 0.754E-18, 0.524E-18, 0.382E-18, 0.272E-18, &
1035 0.202E-18, 0.152E-18, 0.112E-18, 0.876E-19, 0.677E-19, 0.596E-19, &
1036 0.541E-19, 0.510E-19, 0.489E-19, 0.469E-19, 0.448E-19, 0.419E-19, &
1037 0.377E-19, 0.336E-19, 0.285E-19, 0.238E-19, 0.190E-19, 0.145E-19, &
1038 0.106E-19, 0.752E-20, 0.610E-20, 0.553E-20, 0.496E-20, 0.457E-20, &
1039 0.418E-20, 0.380E-20, 0.341E-20, 0.302E-20, 0.279E-20, 0.256E-20, &
1040 0.232E-20, 0.209E-20, 0.171E-20, 0.110E-20, 0.644E-21, 0.416E-21, &
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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1053 DATA (xqy(kl,21),kl=1,130)/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 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000/
1072 DATA (xs(kl,22),kl=1,130)/0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1073 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1074 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1075 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1076 0.340E-20, 0.401E-20, 0.486E-20, 0.573E-20, 0.662E-20, 0.754E-20, &
1077 0.906E-20, 0.112E-19, 0.133E-19, 0.162E-19, 0.202E-19, 0.224E-19, &
1078 0.248E-19, 0.276E-19, 0.289E-19, 0.318E-19, 0.322E-19, 0.321E-19, &
1079 0.338E-19, 0.343E-19, 0.307E-19, 0.290E-19, 0.275E-19, 0.272E-19, &
1080 0.272E-19, 0.272E-19, 0.272E-19, 0.273E-19, 0.280E-19, 0.283E-19, &
1081 0.268E-19, 0.249E-19, 0.212E-19, 0.151E-19, 0.127E-19, 0.142E-19, &
1082 0.229E-19, 0.358E-20, 0.000E+00, 0.000E+00, 0.286E-21, 0.208E-20, &
1083 0.344E-20, 0.764E-20, 0.107E-19, 0.159E-19, 0.166E-19, 0.303E-19, &
1084 0.263E-19, 0.336E-19, 0.366E-19, 0.456E-19, 0.643E-19, 0.546E-19, &
1085 0.922E-19, 0.677E-19, 0.599E-19, 0.117E-18, 0.715E-19, 0.730E-19, &
1086 0.201E-18, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1094 DATA (xqy(kl,22),kl=1,130)/0.400000, 0.400000, 0.400000, 0.400000, &
1095 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1096 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1097 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1098 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1099 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1100 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1101 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1102 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, 0.400000, &
1103 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1104 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1105 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1106 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1107 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1108 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1109 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1110 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1111 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1112 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1113 DATA (xs(kl,23),kl=1,130)/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.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1117 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.164E-20, &
1118 0.211E-20, 0.227E-20, 0.264E-20, 0.340E-20, 0.465E-20, 0.598E-20, &
1119 0.803E-20, 0.986E-20, 0.118E-19, 0.137E-19, 0.160E-19, 0.189E-19, &
1120 0.228E-19, 0.275E-19, 0.307E-19, 0.321E-19, 0.335E-19, 0.349E-19, &
1121 0.363E-19, 0.378E-19, 0.393E-19, 0.407E-19, 0.421E-19, 0.435E-19, &
1122 0.449E-19, 0.462E-19, 0.487E-19, 0.527E-19, 0.564E-19, 0.589E-19, &
1123 0.616E-19, 0.556E-19, 0.553E-19, 0.543E-19, 0.365E-19, 0.318E-19, &
1124 0.316E-19, 0.156E-19, 0.428E-20, 0.113E-20, 0.000E+00, 0.000E+00, &
1125 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
1126 0.000E+00, 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 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
1135 DATA (xqy(kl,23),kl=1,130)/0.019894, 0.019716, 0.019528, 0.019339, &
1136 0.019140, 0.018941, 0.018742, 0.018543, 0.018333, 0.018113, 0.017893, &
1137 0.017673, 0.017453, 0.017212, 0.016982, 0.016741, 0.016531, 0.016238, &
1138 0.015976, 0.015714, 0.015442, 0.015159, 0.014876, 0.014593, 0.014290, &
1139 0.013986, 0.013682, 0.013357, 0.013032, 0.012697, 0.012341, 0.011995, &
1140 0.011629, 0.011314, 0.010874, 0.010487, 0.010078, 0.009670, 0.009240, &
1141 0.008800, 0.008350, 0.007910, 0.007648, 0.007543, 0.007438, 0.007333, &
1142 0.007229, 0.007124, 0.007019, 0.006914, 0.006810, 0.006705, 0.006600, &
1143 0.006495, 0.006286, 0.005867, 0.005343, 0.004819, 0.004295, 0.003772, &
1144 0.003248, 0.002724, 0.002200, 0.001676, 0.001153, 0.000629, 0.000105, &
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 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/
1154 ! END MODULE module_data_photmad
1155 !
1156 CONTAINS
1157 subroutine madronich1_driver(id,ktau,dtstep,config_flags,haveaer, &
1158 gmt,julday,t_phy,moist,aerwrf,p8w,t8w,p_phy, &
1159 chem,rho_phy,dz8w, &
1160 xlat,xlong,z_at_w,gd_cloud,gd_cloud2, &
1161 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
1162 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
1163 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
1164 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
1165 pm2_5_dry,pm2_5_water,uvrad, &
1166 ids,ide, jds,jde, kds,kde, &
1167 ims,ime, jms,jme, kms,kme, &
1168 its,ite, jts,jte, kts,kte )
1169 USE module_configure
1170 USE module_state_description
1171 USE module_model_constants
1172 USE module_data_radm2
1173 implicit none
1174 INTEGER, INTENT(IN ) :: id,julday, &
1175 ids,ide, jds,jde, kds,kde, &
1176 ims,ime, jms,jme, kms,kme, &
1177 its,ite, jts,jte, kts,kte
1178 INTEGER, INTENT(IN ) :: &
1179 ktau
1180 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
1181 INTENT(IN ) :: moist
1182 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1183 INTENT(INOUT ) :: &
1184 pm2_5_dry,pm2_5_water,gd_cloud,gd_cloud2
1185 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1186 INTENT(INOUT ) :: &
1187 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,&
1188 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
1189 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
1190 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob
1191 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
1192 INTENT(IN ) :: chem
1193 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1194 INTENT(IN ) :: &
1195 t_phy, &
1196 p_phy, &
1197 dz8w, &
1198 t8w,p8w,z_at_w , &
1199 aerwrf , &
1200 rho_phy
1201 REAL, DIMENSION( ims:ime , jms:jme ) , &
1202 INTENT(IN ) :: &
1203 xlat, &
1204 xlong
1205 REAL, DIMENSION( ims:ime , jms:jme ) , &
1206 INTENT(INOUT ) :: uvrad
1207 REAL, INTENT(IN ) :: &
1208 dtstep,gmt
1209
1210 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1211
1212 LOGICAL, INTENT(IN) :: haveaer
1213 !
1214 !
1215 ! LOCAL VAR
1216
1217 INTEGER :: ki,i,j,k,ixhour,n,iprt
1218 ! photolysis input
1219 !
1220 real tt(kts:kte),o33(kts:kte),rhoa(kts:kte),aerext(kts:kte),qll(kts:kte), &
1221 phizz(kts:kte),phot1(nreakj-1,kts:kte)
1222 real :: xtime,xhour,xmin,gmtp,uvb_dd1,uvb_du1,uvb_dir1
1223 real :: zenith,zenita,azimuth,dobsi
1224 real :: bext340,bexth2o,ctr
1225 integer :: naerspec
1226 ! print *,'gmt,julday in madronich1= ',gmt,julday
1227 xtime=ktau*dtstep/60.
1228 ixhour=ifix(gmt+.01)+ifix(xtime/60.)
1229 xhour=float(ixhour)
1230 xmin=60.*gmt+(xtime-xhour*60.)
1231 gmtp=mod(xhour,24.)
1232 gmtp=gmtp+xmin/60.
1233 ! print *,'gmtp = ',gmtp,xhour,xmin
1234 do 100 j=jts,jte
1235 do 100 i=its,ite
1236 ! write(0,*)i,j
1237 do k=kts,kte
1238 do n=1,nreakj-1
1239 phot1(n,k)=0.
1240 END DO
1241 END DO
1242 iprt = 0
1243 zenith=0.
1244 zenita=0.
1245 azimuth=0.
1246 call calc_zenith(xlat(i,j),-xlong(i,j),julday,gmtp,azimuth,zenith)
1247 ! if nighttime, skip radiative transfer calculation
1248 if(zenith.eq.90.) zenith = 89.9
1249 if(zenith.ge.90.) go to 199
1250 zenita = cos(zenith*pi/180.)
1251 if(zenith.gt.75.) zenita=1./chap(zenith)
1252 ! photmad berechnet photolysefrequenzen nach Madronich in folgender Reihenfolge
1253
1254 ! o2 absorp 1 schumman-runge corrected in srband
1255 ! o3 -> 1d 2 at 275 k. correct t-dep in subgrid
1256 ! o3 -> 3p 3 at 275 k. correct t-dep in subgrid
1257 ! no2 4
1258 ! no3 -> no+o2 5
1259 ! no3 -> no2+o 6
1260 ! hno2 7
1261 ! hno3 8
1262 ! hno4 9
1263 ! h2o2 10
1264 ! ch2o -> rad 11
1265 ! ch2o -> mol 12
1266 ! ch3cho 13
1267 ! ch3coch3 14
1268 ! ch3coc2h5 15
1269 ! hcocho proc a 16
1270 ! ch3cocho 17
1271 ! hcoch=chcho 18 estimate, no reliable measurement
1272 ! ch3o2h 19
1273 ! ch3coo2h 20 actually use 0.28*(h2o2 value)
1274 ! ch3ono2 21
1275 ! hcocho proc b 22
1276 do k=kts,kte-1
1277 aerext(k+1)=aerwrf(i,k+1,j)
1278 !
1279 !--- if you have aerosols
1280 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1281 bext340=5.E-6
1282 bexth2o=5.E-6
1283 if(haveaer.and.ktau.gt.1)then
1284
1285 ! dry aerosol mass
1286 !rf check ki (or ki+1 or ki-1 ?)
1287 aerext(k)=pm2_5_dry(i,k,j)*bext340+ &
1288 pm2_5_water(i,k,j)*bexth2o
1289 aerext(k)=aerext(k)*1.E3
1290 ! if(i.eq.70.and.j.eq.70) write(06,*) 'aerext',k,aerext(k)
1291
1292 endif
1293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1294 qll(k+1)=0.
1295 tt(k+1) = t_phy(i,k,j)
1296 rhoa(k+1) = rho_phy(i,k,j)
1297 o33(k+1) = max(1.e-3,chem(i,k,j,p_o3))
1298 qll(k+1) = 1.e3*(moist(i,k,j,p_qc)+moist(i,k,j,p_qi)+ &
1299 gd_cloud(i,k,j)+gd_cloud2(i,k,j)) &
1300 *rho_phy(i,k,j)
1301 if(qll(k+1).lt.1.e-5)qll(k+1) = 0.
1302 phizz(k+1) = z_at_w(i,k+1,j)*.001-z_at_w(i,1,j)*.001
1303 ! if((i.eq.1.and.j.eq.17))then
1304 ! 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)
1305 ! 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)
1306 ! endif
1307 END DO
1308 tt(1)=t8w(i,kts,j)
1309 o33(1)=max(1.e-3,chem(i,kts,j,p_o3))
1310 qll(1)=0.
1311 phizz(1)=0.
1312 aerext(1)=aerwrf(i,1,j)
1313 k=0
1314 ! 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)
1315 !
1316 ! if you have aerosols....
1317 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1318 !
1319 if(haveaer)then
1320 aerext(1)=pm2_5_dry(i,1,j)*bext340+ &
1321 pm2_5_water(i,1,j)*bexth2o
1322 aerext(1)=aerext(1)*1.e3
1323 endif
1324 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1325 rhoa(1)=p8w(i,kts,j)/t8w(i,kts,j)/r_d
1326 dobsi=350.
1327 ! if((i.eq.87.and.j.eq.66).or.(i.eq.105.and.j.eq.70))then
1328 ! print *,'before photmad, i,j = ',i,j,pm2_5_dry(i,1,j),pm2_5_water(i,1,j)
1329 ! print *,k,rhoa(1),phizz(1),qll(1),aerext(1),o33(1),tt(1),zenita
1330 ! endif
1331 !write(0,*)'calling photolysis_mad --------',i,j
1332 call photolysis_mad(kte-1,zenita,phizz,tt,rhoa,o33,aerext,qll,dobsi,phot1, &
1333 nreacj,iprt,uvb_dd1,uvb_du1,uvb_dir1)
1334 !write(0,*)'back from photolysis_mad ---------------------------- '
1335 ! print *,'after photmad, i,j = ',i,j
1336 uvrad(i,j)=uvb_dd1+uvb_dir1-uvb_du1
1337 do k=kts,kte-1
1338 do n=1,nreakj-1
1339 phot1(n,k)=60.*phot1(n,k)
1340 END DO
1341 END DO
1342 199 continue
1343 do k=kts,kte-1
1344 !
1345 !
1346 ph_o31d(i,k,j) = phot1(1,k)
1347 ph_o33p(i,k,j) = phot1(2,k)
1348 ph_no2(i,k,j) = phot1(3,k)
1349 ph_no3o2(i,k,j) = phot1(4,k)
1350 ph_no3o(i,k,j) = phot1(5,k)
1351 ph_hno2(i,k,j) = phot1(6,k)
1352 ph_hno3(i,k,j) = phot1(7,k)
1353 ph_hno4(i,k,j) = phot1(8,k)
1354 ph_h2o2(i,k,j) = phot1(9,k)
1355 ph_ch2or(i,k,j) = phot1(10,k)
1356 ph_ch2om(i,k,j) = phot1(11,k)
1357 ph_ch3cho(i,k,j) = phot1(12,k)
1358 ph_ch3coch3(i,k,j) = phot1(13,k)
1359 ph_ch3coc2h5(i,k,j) = phot1(14,k)
1360 ph_hcocho(i,k,j) = phot1(15,k)
1361 ph_ch3cocho(i,k,j) = phot1(16,k)
1362 ph_hcochest(i,k,j) = phot1(17,k)
1363 ph_ch3o2h(i,k,j) = phot1(18,k)
1364 ph_ch3coo2h(i,k,j) = phot1(19,k)
1365 ph_ch3ono2(i,k,j) = phot1(20,k)
1366 ph_hcochob(i,k,j) = phot1(21,k)
1367 ph_macr(i,k,j) = phot1(22,k)
1368 ! if(i.eq.5.and.j.eq.5)print *,i,j,k,phot1(3,k),phot1(4,k),phot1(17,k)
1369 END DO
1370 100 continue
1371
1372 END SUBROUTINE madronich1_driver
1373
1374 SUBROUTINE photolysis_mad(mkxcc,a,zmm5,tmm5,pmm5,o3mm5,aerext,wlmm5, &
1375 dobsnew,phot1,nrtest,iprt,uvb_dd1,uvb_du1,uvb_dir1)
1376 ! Note that nreakj can differ from nreacj in csolve1!!!!!!
1377 ! ----- INPUT ------------------------------------------------
1378 ! Input from MM5
1379 ! ------ OUTPUT ----------------------------------------------
1380 ! -----------------------------------------------------------------------
1381 ! This program calculates J-values for selected atmospheric molecules
1382 ! the program and subroutine structure is as follows
1383 ! runph -main program and menu.
1384 ! calc_zenith -zenith angle calc.
1385 ! in addition, the chapman function ch(zenith) is used between
1386 ! readd -reads spectra and standard profiles
1387 ! photmat -main subroutine (calls the routines below)
1388 ! o3scal -rescales ozone profile to user-selected new dobson
1389 ! subgrid -regrids the altitude profiles (
1390 ! srband -computes effective ozone cross sections in the
1391 ! schumamkxcc+1+nabv-runge region (if necessary)
1392 ! trapez -interpolation subroutine
1393 ! optics -computes optical parameters and calls delta-eddington
1394 ! delted -delta-eddington code of wiscombe
1395 ! nlayde(ib) -multylayer calc
1396 ! leqt1b(a,n,ncl,nuc,ia,b,m,ib,ijob,kl) -solves
1397 ! pentadiag.system
1398 ! Original program written by S. Madronich, NCAR. The code was
1399 ! last modified by him on 18 Aug. 1987.
1400 ! The driver for the chemistry has been heavily modified by
1401 ! W. R. Stockwell at IFU Germany. (The modifications allow
1402 ! the simulation conditions to be directly modified
1403 ! without recompiling the program [this does not hold for the present
1404 ! version here any more]). The extremely large data base
1405 ! has bee replaced with separate files to allow cross sections and
1406 ! photolysis rates to be more easily updated. However, the order
1407 ! of the input files should be modified with extreme caution.
1408 ! The present version is made for the use of computed
1409 ! fields of T, P, O3, wl, and aerosol (if available)
1410 ! Modifications were made by Renate Forkel in Aug./Sept 1995
1411 ! The following modifications were made:
1412 ! - No climatological input any more
1413 ! - Arbitrary vertical grid
1414 ! - Several seperate cloud layers are possible. Number of additional
1415 ! layers within the clouds depends on the optical depth
1416 ! - Vertically inhogeneous clouds are possible now
1417 ! - No variable values in commons any more
1418 ! -----------------------------------------------------------------------
1419 ! xnkg: Molecules per kg air
1420 ! .. Scalar Arguments ..
1421 REAL :: a, dobsnew,uvb_dd1,uvb_du1,uvb_dir1
1422 INTEGER :: iprt, mkxcc, nrtest
1423 ! .. Local Scalars ..
1424 REAL :: aeru, airu, bextn, df, dj, dz, ff, gaer, gcld, gray, haer, &
1425 hair, ho3, o3u, omaer, omcld, omray, reff, tu, wmicron, xnkg, xx, &
1426 znorm, zu
1427 INTEGER :: i, j, k, kk, kl, lev, nlayer, nlevel, nn, nr, nsurf
1428 ! .. Array Arguments ..
1429 REAL :: aerext(mkxcc+1), o3mm5(mkxcc+1), phot1(nreakj-1,mkxcc+1), &
1430 pmm5(mkxcc+1), tmm5(mkxcc+1), wlmm5(mkxcc+1), zmm5(mkxcc+1)
1431 ! .. Local Arrays ..
1432 REAL :: aaer(130), aer(mkxcc+1+nabv), air(mkxcc+1+nabv), ao2(nj,130), &
1433 ao3(nj,130), arayl(130), cloud(mkxcc+1+nabv), cvo2(nj), &
1434 d(nreakj,nj), endir(nj,130), endn(nj,130), enup(nj,130), &
1435 hilf1(mkxcc+1), hilfd(nj), o3(mkxcc+1+nabv), qy(nj,130,nreakj), &
1436 s(nj,130,nreakj), t(mkxcc+1+nabv), vaer(nj), vair(nj), vcld(nj), &
1437 vo3(nj), vt(nj), z(nj), zkm(mkxcc+1), zmid(nj), zz(mkxcc+1+nabv)
1438 ! .. Data Statements ..
1439 DATA xnkg/2.143E25/
1440 ! EXTERNAL o3scal, optics, srband, subgrid, trapez
1441 ! wave length range !!! If photolysis is also desired for levels above
1442 ! 2
1443 ! kl0 should be set equal to 1 again!!!!!!!!!!!!!!
1444 i = 1
1445 j = 1
1446 ! albedoph ************************ specify ground albedoph
1447 ! use best estimate albedoph of demerjian et al.,
1448 ! adv.env.sci.tech.,v.10,p.369, (1980)
1449 DO kl = kl0, kl1
1450
1451 IF (wl(kl)<400.) albedoph(kl) = 0.05
1452
1453 IF ((wl(kl)>=400.) .AND. (wl(kl)<450.)) albedoph(kl) = 0.06
1454
1455 IF ((wl(kl)>=450.) .AND. (wl(kl)<500.)) albedoph(kl) = 0.08
1456
1457 IF ((wl(kl)>=500.) .AND. (wl(kl)<550.)) albedoph(kl) = 0.10
1458
1459 IF ((wl(kl)>=550.) .AND. (wl(kl)<600.)) albedoph(kl) = 0.11
1460
1461 IF ((wl(kl)>=600.) .AND. (wl(kl)<640.)) albedoph(kl) = 0.12
1462
1463 IF ((wl(kl)>=640.) .AND. (wl(kl)<660.)) albedoph(kl) = 0.135
1464
1465 IF (wl(kl)>=660.) albedoph(kl) = 0.15
1466 END DO
1467
1468
1469
1470 nn = mkxcc + 1 + nabv
1471
1472 ! omray = single scattering albedoph, rayleigh. use 1.00
1473 ! gray = asymetry factor for rayleigh scattering. use 0.0
1474 ! arayl(kl) = rayleigh scattering cross section, from
1475 ! frohlich and shaw, appl.opt. v.11, p.1773 (1980).
1476 ! overrides tabulation of jdata.base
1477 hair = 8.05
1478 omray = 1.0
1479 gray = 0.0
1480
1481 DO 10 kl = kl0, kl1
1482 wmicron = wl(kl)/1.E3
1483 xx = 3.916 + 0.074*wmicron + 0.050/wmicron
1484 arayl(kl) = 3.90E-28/(wmicron)**xx
1485 10 CONTINUE
1486
1487 ! aerosol*********************** specify aerosols
1488 ! aaer(kl) = aerosol total vertical optical depth variation with
1489 ! wavelength. estimated from elterman (1968)
1490 ! aer(i) = attenuation (per km) profile from elterman (1968).
1491 ! given in data statement in beginning of code, for 340 nm (kl=6
1492 ! same vertical shape at all wavelengths.
1493 ! normalized later (in subroutine subgrid) to total vertical dep
1494 ! this wavelength.
1495 ! omaer = aerosol single scattering albedoph. use 0.99 for now.
1496 ! gaer = aerosol asymetry factor. use 0.61 (hansen and travis 197
1497 ! (these are assuming particles of about 0.1 micron radius
1498 ! index of refraction of about 1.65 + 0.002i.
1499 ! haer = the aerosol scale height at top of atmosphere
1500 ! use equal to air (8.05 km)
1501
1502 DO 20 kl = kl0, kl1
1503 aaer(kl) = 0.379*(340./wl(kl))
1504 20 CONTINUE
1505 omaer = 0.990
1506 gaer = 0.610
1507 haer = 8.05
1508
1509
1510 omcld = 1.000
1511 gcld = 0.860
1512
1513 ho3 = 4.50
1514
1515 nsurf = 1
1516
1517 ! Vertikales Gitter ohne Wolken: Niveaus zz(i), i=1,mkxcc+1+nabv
1518
1519
1520 ! Transformation of MM5-values
1521 ! bextn: cloud extinction coeff per g/m**3 [1/km]
1522 DO k = 1, mkxcc + 1
1523 zz(k) = zmm5(k)
1524 !write(0,*)' here7a zz = zmm5 ',k,zmm5(k)
1525 t(k) = tmm5(k)
1526 air(k) = xnkg*pmm5(k)*1.E-6
1527
1528 ! falls o3mm5 in ppm
1529
1530 o3(k) = o3mm5(k)*1.E-6*air(k)
1531 aer(k) = aerext(k)
1532
1533 ! bextn: cloud extinction coeff per g/m**3 [1/km]
1534 ! **** warning: parameterization for bextn only good
1535 ! for continental clouds
1536
1537 IF (wlmm5(k)>0.) THEN
1538
1539 ! vereinfachte Version, falls kein
1540 ! Sulfat uebergeben wird.
1541
1542 reff = 9.6*wlmm5(k)**0.333
1543 bextn = (0.0275+1.3/reff)*1000.
1544 cloud(k) = wlmm5(k)*bextn
1545 ELSE
1546 cloud(k) = 0.
1547 END IF
1548 ! if(iprt.eq.1)write(6,'(i3,e12.4)')k,cloud(k)
1549 END DO
1550
1551 znorm = (50.-zz(mkxcc+1))/(50.-20.)
1552
1553 DO k = 1, nabv
1554 zz(mkxcc+1+k) = 50. - znorm*(50.-zabv(k))
1555 !write(0,*)' here8a ',k,' zz(',mkxcc+1+k,') ',zz(mkxcc+1+k),znorm,zabv(k)
1556 END DO
1557
1558 zu = zz(mkxcc+1)
1559 tu = t(mkxcc+1)
1560 o3u = o3(mkxcc+1)
1561 airu = air(mkxcc+1)
1562 aeru = aer(mkxcc+1)
1563 kk = 1
1564
1565 ! Zufuegen von Werten oberhalb von MM5
1566
1567 ! Die 'abv'-Werte sind bereits in den richtigen Einheiten
1568 DO k = mkxcc + 1 + 1, mkxcc + 1 + nabv
1569 ! write(0,'(2i3,5e12.3)') k,kk,zz(k),zabv(kk)
1570 30 IF (zz(k)<=zabv(kk)) THEN
1571 dz = zz(k) - zu
1572 ff = dz/(zabv(kk)-zu)
1573 t(k) = tu + ff*(tabv(kk)-tu)
1574 o3(k) = o3u + ff*(o3abv(kk)-o3u)
1575 air(k) = airu + ff*(pabv(kk)-airu)
1576 aer(k) = aeru + ff*(caabv(kk)-aeru)
1577 cloud(k) = 0.
1578 ! if(iprt.eq.1)then
1579 ! write(0,'(2i3,5e12.3)') k,kk,zz(k),zabv(kk),ff,tabv(kk),air(k)
1580 ! endif
1581 ELSE
1582 40 zu = zabv(kk)
1583 tu = tabv(kk)
1584 o3u = o3abv(kk)
1585 airu = pabv(kk)
1586 aeru = caabv(kk)
1587 kk = kk + 1
1588
1589 IF (zabv(kk)<zz(k)) GO TO 40
1590 GO TO 30
1591 END IF
1592
1593 END DO
1594
1595 IF (iprt==1) THEN
1596 !WRITE (06,'('k, zz(k), t(k), air(k), o3(k), aer(k) cloud(k)')')
1597
1598 !DO k = 1, mkxcc + 1 + nabv
1599 ! WRITE (06,'(i3,6e12.4)') k, zz(k), t(k), air(k), o3(k), aer(k), &
1600 ! cloud(k)
1601 !END DO
1602
1603 END IF
1604 ! if(iprt.eq.1)then
1605 ! do k=1,mkxcc+1
1606 ! write(06,'(i3,6e12.4)') k,zz(k),t(k),air(k),o3mm5(k), &
1607 ! aerext(k),cloud(k)
1608 ! enddo
1609 ! endif
1610 ! Scaling of O3-column with specified Dobson units
1611 ! -----------------------------------------------------------------------
1612 CALL o3scal(dobsnew,ho3,zz,o3,mkxcc+1+nabv)
1613 ! stop
1614 ! -----------------------------------------------------------------------
1615
1616 ! Einfuegen zusaetzlicher Schichten in Wolken,
1617 ! T-abh. Absorptionsquerschnitte, etc
1618
1619 ! -----------------------------------------------------------------------
1620 nlevel = 0
1621
1622 IF (iprt==1) PRINT *, 'nlevel = ', nlevel
1623
1624 CALL subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
1625 vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,mkxcc+1,mkxcc+1+nabv)
1626
1627 IF (iprt==1) PRINT *, 'nlevel = ', nlevel
1628 ! -----------------------------------------------------------------------
1629
1630 ! subroutine srband computes effective ozone cross sections in the
1631 ! schumann-runge region, using the parameterization of allen and
1632 ! freder
1633
1634 ! ----------------------------------------------------------------------
1635 CALL srband(a,zmid,cvo2,ao2,vt,nlevel)
1636 ! ----------------------------------------------------------------------
1637
1638
1639 ! subroutine optics is the driver for the flux calculation. output is
1640 ! endir(lev,kl) -irradiance of direct solar beam
1641 ! endn(lev,kl) -irradiance of down-welling diffuse light
1642 ! enup(lev,kl) -irradiance of up-welling diffuse light
1643 ! for each level lev and wavelength bin kl. conversion to actinic
1644 ! flux is in loop 90
1645
1646 ! ----------------------------------------------------------------------
1647 CALL optics(iprt,a,vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld,omcld, &
1648 vaer,aaer,gaer,omaer,nlayer,nlevel,nsurf,endir,endn,enup)
1649 ! ----------------------------------------------------------------------
1650
1651 ! Actinic flux is simply computed from irradiance, using a factor of 2
1652 ! for the diffuse part. This is only a very rough estimate according
1653 ! to
1654 ! Ruggaber et al. (1993), JGR 98D, 1151-1162
1655
1656 DO nr = 1, nreakj
1657
1658 DO k = 1, nlevel
1659 d(nr,k) = 0.
1660 END DO
1661
1662 END DO
1663
1664 DO 50 kl = kl0, kl1
1665
1666 DO 50 lev = 1, nlevel
1667 df = fext(kl)*(endir(lev,kl)/a+2.*(endn(lev,kl)+enup(lev,kl)))
1668
1669 DO 50 nr = 1, nreakj
1670 dj = df*s(lev,kl,nr)*qy(lev,kl,nr)
1671 d(nr,lev) = d(nr,lev) + dj
1672 50 CONTINUE
1673 uvb_dd1=0.
1674 uvb_du1=0.
1675 uvb_dir1=0.
1676 do kl = 38,54
1677 ! radfakt=fext(kl) ! photons/cm**2/s/bin
1678 ! ! UV in Photonen/cm**2/s
1679 radfakt=fext(kl)*1e4 * 6.63e-34 * 2.9979e8 / (wl(kl)*1e-9)
1680 ! UV in W/m**2
1681 uvb_dd1=uvb_dd1+radfakt*endn(1,kl)
1682 uvb_du1=uvb_du1+radfakt*enup(1,kl)
1683 uvb_dir1=uvb_dir1+radfakt*endir(1,kl)
1684 enddo
1685
1686 ! print *,'d(4,1) = ',d(4,1)
1687
1688 DO nr = 1, nreakj - 1
1689 ! write(06,'(' nr',i3)') nr
1690 DO lev = 1, mkxcc + 1
1691 zkm(lev) = zmm5(lev)
1692 END DO
1693
1694 DO lev = 1, nlevel
1695 hilfd(lev) = d(nr+1,lev)
1696 ! if(iprt.eq.1.and.nr.eq.3)
1697 ! 1 write(06,'(i3,2e12.4)') lev,z(lev),d(4,lev)
1698 END DO
1699 ! hilfd muss noch auf mm5-hoehen zurueckinterpoliert werden
1700 CALL trapez(z,hilfd,1,nlevel,zkm,hilf1,1,mkxcc+1,nj,nj,mkxcc+1, &
1701 mkxcc+1)
1702
1703 IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '
1704
1705 IF (iprt==1 .AND. nr==3) WRITE (06,'(5e12.6)') (hilf1(k),k=1, &
1706 mkxcc+1)
1707
1708 IF (iprt==1 .AND. nr==3) PRINT *, 'nach trapez '
1709
1710 DO lev = 1, mkxcc + 1
1711 phot1(nr,lev) = hilf1(lev)
1712 END DO
1713
1714 END DO
1715
1716 RETURN
1717
1718 END SUBROUTINE photolysis_mad
1719
1720 ! #################################################################
1721
1722 SUBROUTINE srband(a,zmid,cvo2,ao2,vt,nlevel)
1723 ! this subroutine calculates effective o2 cross sections in the
1724 ! schumann-runge band, using the formulae from allen and frederick,
1725 ! j.atmos.sci., v39, p2066 (1982). the coefficients are those modifi
1726 ! wuebbles (lll report, 1982).
1727 ! .. Scalar Arguments ..
1728 REAL :: a
1729 INTEGER :: nlevel
1730 ! .. Local Scalars ..
1731 REAL :: ao20, ao20lg, c, clog, e10, x1, x2, x3, x4, x5, x6, x7, x8, &
1732 y1, y2, y3, y4, zendep
1733 INTEGER :: lay, n20, nlayer
1734 ! .. Intrinsic Functions ..
1735 INTRINSIC alog
1736 ! .. Array Arguments ..
1737 REAL :: ao2(nj,130), cvo2(nj), vt(nj), zmid(nj)
1738
1739 nlayer = nlevel - 1
1740 ! initialize cross sections:
1741 DO 10 lay = 1, nlayer
1742
1743 DO 10 kl = kl0, kl1
1744 10 ao2(lay,kl) = xs(kl,1)
1745
1746
1747 ! correct as needed
1748 ! use formula for 20-50 km.
1749 ! below 20 km, use 20 km value
1750 ! find layer near 20 km
1751 DO 20 lay = 1, nlayer
1752
1753 IF (zmid(lay)>20.) THEN
1754 n20 = lay
1755 GO TO 30
1756 END IF
1757
1758 20 CONTINUE
1759 30 CONTINUE
1760
1761 e10 = alog(10.)
1762
1763 DO 60 kl = kl0, kl1
1764
1765 IF (wl(kl)>205.) RETURN
1766
1767 DO 40 lay = n20, nlayer
1768 x1 = alog(4.696E-23*cvo2(lay)/0.2095)/e10
1769
1770 IF (wl(kl)>=200.) x1 = vt(lay)
1771 x2 = x1*x1
1772 x3 = x2*x1
1773 x4 = x3*x1
1774 x5 = x4*x1
1775 x6 = x5*x1
1776 x7 = x6*x1
1777 x8 = x7*x1
1778 ao20lg = sra(kl,1) + sra(kl,2)*x1 + sra(kl,3)*x2 + sra(kl,4)*x3 + &
1779 sra(kl,5)*x4 + sra(kl,6)*x5 + sra(kl,7)*x6 + sra(kl,8)*x7 + &
1780 sra(kl,9)*x8
1781 ao20 = 10.**ao20lg
1782
1783 y1 = alog(cvo2(lay))/e10
1784 y2 = y1*y1
1785 y3 = y2*y1
1786 y4 = y3*y1
1787 clog = srb(kl,1) + srb(kl,2)*y1 + srb(kl,3)*y2 + srb(kl,4)*y3 + &
1788 srb(kl,5)*y4
1789 c = 10.**clog
1790 zendep = a**c
1791
1792 ao2(lay,kl) = ao20*zendep
1793 40 CONTINUE
1794 ! assign values below 20 km
1795 DO 50 lay = 1, n20 - 1
1796 ao2(lay,kl) = ao2(n20,kl)
1797 50 CONTINUE
1798 60 CONTINUE
1799
1800 RETURN
1801
1802 END SUBROUTINE srband
1803
1804 ! ######################################################################
1805
1806 SUBROUTINE o3scal(dobsnew,ho3,zz,o3,nn)
1807 ! adjustment of o3 profiles to a user-selected dobson value.
1808 ! select value of dobnew in main program
1809 ! if don't want to use, don't call this subroutine
1810 ! .. Scalar Arguments ..
1811 REAL :: dobsnew, ho3
1812 INTEGER :: nn
1813 ! .. Local Scalars ..
1814 REAL :: dobsref
1815 INTEGER :: i
1816 ! .. Intrinsic Functions ..
1817 INTRINSIC max, min
1818 ! .. Array Arguments ..
1819 REAL :: o3(nn), zz(nn)
1820 ! write(6,*) o3
1821 dobsref = o3(nn)*1.E5*ho3
1822 ! write(06,'('nn: dobsref,dobsnew',2e12.4)')
1823 ! & dobsref/2.687e16,dobsnew
1824 DO 10 i = 1, nn
1825 10 dobsref = dobsref + o3(i)*0.5*(zz(min(i+1,nn))-zz(max(i-1,1)))*1.E5
1826 dobsref = dobsref/2.687E16
1827 ! write(06,'('dobsref,dobsnew',2e12.4)') dobsref,dobsnew
1828 DO 20 i = 1, nn
1829 o3(i) = o3(i)*dobsnew/dobsref
1830 20 CONTINUE
1831 ! write(06,*) o3
1832
1833 RETURN
1834
1835 END SUBROUTINE o3scal
1836
1837 ! #######################################################################
1838
1839 SUBROUTINE leqt1b(a,n,nlc,nuc,ia,b,m,ib,ijob,xl)
1840 ! -leqt1b--------s-------library
1841 ! 3---------------------------------------
1842 ! function - matrix decomposition, linear equation
1843 ! solution - space economizer solution -
1844 ! band storage mode
1845 ! usage - call leqt1b (a,n,nlc,nuc,ia,b,m,ib,ijob,xl,
1846 ! ier)
1847 ! parameters a - input/output matrix of dimension n by
1848 ! (nuc+nlc+1). see parameter ijob.
1849 ! n - order of matrix a and the number of rows in
1850 ! b. (input)
1851 ! nlc - number of lower codiagonals in matrix a.
1852 ! (input)
1853 ! nuc - number of upper codiagonals in matrix a.
1854 ! (input)
1855 ! ia - row dimension of a as specified in the
1856 ! calling program. (input)
1857 ! b - input/output matrix of dimension n by m.
1858 ! on input, b contains the m right-hand sides
1859 ! of the equation ax = b. on output, the
1860 ! solution matrix x replaces b. if ijob = 1,
1861 ! b is not used.
1862 ! m - number of right hand sides (columns in b).
1863 ! (input)
1864 ! ib - row dimension of b as specified in the
1865 ! calling program. (input)
1866 ! ijob - input option parameter. ijob = i implies when
1867 ! i = 0, factor the matrix a and solve the
1868 ! equation ax = b. on input, a contains the
1869 ! coefficient matrix of the equation ax = b,
1870 ! where a is assumed to be an n by n band
1871 ! matrix. a is stored in band storage mode
1872 ! and therefore has dimension n by
1873 ! (nlc+nuc+1). on output, a is replaced
1874 ! by the u matrix of the l-u decomposition
1875 ! of a rowwise permutation of matrix a. u is
1876 ! stored in band storage mode.
1877 ! i = 1, factor the matrix a. a contains the
1878 ! same input/output information as if
1879 ! ijob = 0.
1880 ! i = 2, solve the equation ax = b. this
1881 ! option implies that leqt1b has already
1882 ! been called using ijob = 0 or 1 so that
1883 ! the matrix a has already been factored.
1884 ! in this case, output matrices a and xl
1885 ! must have been saved for reuse in the
1886 ! call to leqt1b.
1887 ! xl - work area of dimension n*(nlc+1). the first
1888 ! nlc*n locations of xl contain components of
1889 ! the l matrix of the l-u decomposition of a
1890 ! rowwise permutation of a. the last n
1891 ! locations contain the pivot indices.
1892 ! -----------------------------------------------------------------------
1893 ! latest revision - november 27,1973
1894 ! dimension a(ia,1),xl(n,1),b(ib,1) ! Urspr. Zustand, fun
1895 ! .. Scalar Arguments ..
1896 INTEGER :: ia, ib, ijob, m, n, nlc, nuc
1897 ! .. Array Arguments ..
1898 REAL :: a(ia,5), b(ib,5), xl(n,5)
1899 ! .. Local Scalars ..
1900 REAL :: one, p, q, rn, zero
1901 INTEGER :: i, ik, j, jbeg, jend, k, k1, kk, l, nc, nlc1, nn
1902 ! .. Intrinsic Functions ..
1903 INTRINSIC abs
1904 ! .. Data Statements ..
1905 DATA zero/0./, one/1.0/
1906
1907 p = 0
1908 jbeg = nlc + 1
1909 nlc1 = jbeg
1910
1911 IF (ijob==2) GO TO 170
1912 rn = n
1913 ! restructure the matrix
1914 ! find reciprocal of the largest
1915 ! absolute value in row i
1916 i = 1
1917 nc = jbeg + nuc
1918 nn = nc
1919 jend = nc
1920
1921 IF (n==1 .OR. nlc==0) GO TO 50
1922 10 k = 1
1923 p = zero
1924
1925 DO 20 j = jbeg, jend
1926 a(i,k) = a(i,j)
1927 q = abs(a(i,k))
1928
1929 IF (q>p) p = q
1930 k = k + 1
1931 20 CONTINUE
1932
1933 IF (p==zero) GO TO 280
1934 xl(i,nlc1) = one/p
1935
1936 IF (k>nc) GO TO 40
1937
1938 DO 30 j = k, nc
1939 a(i,j) = zero
1940 30 CONTINUE
1941 40 i = i + 1
1942 jbeg = jbeg - 1
1943
1944 IF (jend-jbeg==n) jend = jend - 1
1945
1946 IF (i<=nlc) GO TO 10
1947 jbeg = i
1948 nn = jend
1949 50 jend = n - nuc
1950
1951 DO 90 i = jbeg, n
1952 p = zero
1953
1954 DO 60 j = 1, nn
1955 q = abs(a(i,j))
1956
1957 IF (q>p) p = q
1958 60 CONTINUE
1959
1960 IF (p==zero) GO TO 280
1961 xl(i,nlc1) = one/p
1962
1963 IF (i==jend) GO TO 80
1964
1965 IF (i<jend) GO TO 90
1966 k = nn + 1
1967
1968 DO 70 j = k, nc
1969 a(i,j) = zero
1970 70 CONTINUE
1971 80 nn = nn - 1
1972 90 CONTINUE
1973 l = nlc
1974 ! l-u decomposition
1975 DO 160 k = 1, n
1976 p = abs(a(k,1))*xl(k,nlc1)
1977 i = k
1978
1979 IF (l<n) l = l + 1
1980 k1 = k + 1
1981
1982 IF (k1>l) GO TO 110
1983
1984 DO 100 j = k1, l
1985 q = abs(a(j,1))*xl(j,nlc1)
1986
1987 IF (q<=p) GO TO 100
1988 p = q
1989 i = j
1990 100 CONTINUE
1991 110 xl(i,nlc1) = xl(k,nlc1)
1992 xl(k,nlc1) = i
1993 ! singularity found
1994 IF ((rn+p)==rn) GO TO 280
1995 ! interchange rows i and k
1996 IF (k==i) GO TO 130
1997
1998 DO 120 j = 1, nc
1999 p = a(k,j)
2000 a(k,j) = a(i,j)
2001 a(i,j) = p
2002 120 CONTINUE
2003
2004 130 IF (k1>l) GO TO 160
2005
2006 DO 150 i = k1, l
2007 p = a(i,1)/a(k,1)
2008 ik = i - k
2009 xl(k1,ik) = p
2010
2011 DO 140 j = 2, nc
2012 a(i,j-1) = a(i,j) - p*a(k,j)
2013 140 CONTINUE
2014 a(i,nc) = zero
2015 150 CONTINUE
2016 160 CONTINUE
2017
2018 IF (ijob==1) GO TO 270
2019 ! forward substitution
2020 170 l = nlc
2021
2022 DO 220 k = 1, n
2023 i = xl(k,nlc1)
2024
2025 IF (i==k) GO TO 190
2026
2027 DO 180 j = 1, m
2028 p = b(k,j)
2029 b(k,j) = b(i,j)
2030 b(i,j) = p
2031 180 CONTINUE
2032
2033 190 IF (l<n) l = l + 1
2034 k1 = k + 1
2035
2036 IF (k1>l) GO TO 220
2037
2038 DO 210 i = k1, l
2039 ik = i - k
2040 p = xl(k1,ik)
2041
2042 DO 200 j = 1, m
2043 b(i,j) = b(i,j) - p*b(k,j)
2044 200 CONTINUE
2045 210 CONTINUE
2046 220 CONTINUE
2047 ! backward substitution
2048 jbeg = nuc + nlc
2049
2050 DO 260 j = 1, m
2051 l = 1
2052 k1 = n + 1
2053
2054 DO 250 i = 1, n
2055 k = k1 - i
2056 p = b(k,j)
2057
2058 IF (l==1) GO TO 240
2059
2060 DO 230 kk = 2, l
2061 ik = kk + k
2062 p = p - a(k,kk)*b(ik-1,j)
2063 230 CONTINUE
2064 240 b(k,j) = p/a(k,1)
2065
2066 IF (l<=jbeg) l = l + 1
2067 250 CONTINUE
2068 260 CONTINUE
2069
2070 270 RETURN
2071
2072 280 CONTINUE
2073 CALL wrf_error_fatal ( ' leqt1b error--matrix algorithmically singular')
2074 END SUBROUTINE leqt1b
2075
2076 ! #######################################################################
2077
2078 FUNCTION chap(zeta)
2079 ! chapman function is used when the solar zenith angle exceeds
2080 ! 75 deg.
2081 ! interpolates between values given in, e.g., mccartney (1976).
2082 ! .. Scalar Arguments ..
2083 REAL :: zeta
2084 ! .. Local Scalars ..
2085 REAL :: rm
2086 INTEGER :: i
2087 ! .. Local Arrays ..
2088 REAL :: y(22)
2089 ! .. Function Return Value ..
2090 REAL :: chap
2091 ! .. Data Statements ..
2092 DATA (y(i),i=1,22)/3.800, 4.055, 4.348, 4.687, 5.083, 5.551, 6.113, &
2093 6.799, 7.650, 8.732, 10.144, 12.051, 14.730, 18.686, 24.905, 35.466, &
2094 55.211, 96.753, 197., 485., 1476., 9999./
2095
2096 DO 10 i = 75, 96
2097 rm = i
2098
2099 IF (zeta<rm) GO TO 20
2100 10 CONTINUE
2101 20 chap = y(i-75) + (y(i-74)-y(i-75))*(zeta-(rm-1.))
2102 END FUNCTION chap
2103
2104
2105 ! #################################################################
2106 SUBROUTINE subgrid(zz,t,air,aer,o3,cloud,hair,haer,z,zmid,vair,vo3,vaer, &
2107 vcld,vt,ao3,nlevel,nlayer,s,qy,cvo2,omray,nz,nn)
2108 ! Output
2109 ! this subroutine computes all altitude dependent quantities over the
2110 ! sub-divided grid
2111 ! the following quantities are computed at each level (altitude)
2112 ! z(lev) = altitude (km) of each level
2113 ! zair(lev) = air concentration at each level
2114 ! zt(lev) = temperature at each level
2115 ! s(lev,kl,nr) = abs. cross sect at each altitude (t,p corrected)
2116 ! qy(lev,kl,nr) = quantum yield at each altitude (t,p corrected)
2117 ! the following quantities are computed at each layer (thickness)
2118 ! zmid(lay) = altitude of midpoint of layer
2119 ! vair(lay) = air column in layer (vertical)
2120 ! vo3(lay) = ozone column ' '
2121 ! vaer(lay) = aerosol ' ' '
2122 ! vcld(lay) = cloud ' ' '
2123 ! vt(lay) = average temperature of column
2124 ! ao3(lay,kl) = average o3 cross sect in layer (with ave. layer te
2125 ! .. Scalar Arguments ..
2126 REAL :: haer, hair, omray
2127 INTEGER :: nlayer, nlevel, nn, nz
2128 ! .. Local Scalars ..
2129 REAL :: a, ak300, akt, b, c, dz, dzo, dzt, dzu, fnum, hlocal, phi1, &
2130 phi2, phi20, sum, tau, tdiffx, x0, xl, xl0
2131 INTEGER :: i, idt, ii, il, lay, lev, llbt, nr
2132 ! .. Intrinsic Functions ..
2133 ! EXTERNAL trapez
2134 INTRINSIC alog, atan, exp, float, ifix, max
2135 ! .. Array Arguments ..
2136 REAL :: aer(nn), air(nn), ao3(nj,130), cloud(nn), cvo2(nj), o3(nn), &
2137 qy(nj,130,nreakj), s(nj,130,nreakj), t(nn), vaer(nj), vair(nj), &
2138 vcld(nj), vo3(nj), vt(nj), z(nj), zmid(nj), zz(nn)
2139 ! .. Local Arrays ..
2140 REAL :: zair(nj), zt(nj)
2141 ! ------------------------------------------------------ levels
2142 ! fnum: factor for calculation of additional levels in clouds
2143 fnum = 0.02
2144
2145 DO i = 1, nj
2146 vcld(i) = 0.
2147 END DO
2148
2149 lev = 0
2150 llbt = 0
2151
2152 DO 30 i = 1, nn
2153 ! write(06,'('cloud',i3,2e12.4)') i,zz(i),cloud(i)
2154 IF (cloud(i)<=0.) THEN
2155 lev = lev + 1
2156 z(lev) = zz(i)
2157 zair(lev) = air(i)
2158 zt(lev) = t(i)
2159 ! write(06,'('z,t,air ',i3,3e12.4)')
2160 ! lev,z(lev),zt(lev),zair(lev)
2161 ELSE
2162
2163 IF (i==1 .AND. lev==0) THEN
2164 llbt = llbt + 1
2165 ! lbase(llbt)=1
2166 lev = 1
2167 z(lev) = zz(i)
2168 zair(lev) = air(i)
2169 zt(lev) = t(i)
2170 END IF
2171
2172 IF (i==1) GO TO 10
2173
2174 IF (cloud(i-1)<=0.) THEN
2175 lev = lev + 1
2176 llbt = llbt + 1
2177 ! lbase(llbt)=lev
2178 z(lev) = 0.5*(zz(i)+zz(i-1))
2179 zt(lev) = 0.5*(t(i)+t(i-1))
2180 if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
2181 zair(lev) = air(i-1)
2182 else
2183 hlocal = 1./alog(air(i-1)/air(i))
2184 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
2185 zair(lev) = air(i-1)*exp(-x0/hlocal)
2186 endif
2187
2188 ! write(06,'('u:z,t,air ',i3,3e12.4)')
2189 ! lev,z(lev),zt(lev),zair(lev)
2190 END IF
2191
2192 dzt = zz(i) - z(lev)
2193 ! number of layers depents on optical depth
2194 idt = max(ifix(cloud(i)*dzt*fnum),1)
2195 dzu = dzt/float(idt)
2196 vcld(lev) = cloud(i)
2197
2198 DO il = 1, idt
2199 lev = lev + 1
2200
2201 IF (lev>nj) THEN
2202 CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
2203 END IF
2204
2205 z(lev) = z(lev-1) + dzu
2206 zt(lev) = t(i-1) + (z(lev)-zz(i-1))/(zz(i)-zz(i-1))*(t(i)-t(i-1) &
2207 )
2208
2209 if(abs(air(i)-air(i-1)).lt.air(i-1)/1.e5)then
2210 zair(lev) = air(i-1)
2211 else
2212 hlocal = 1./alog(air(i-1)/air(i))
2213 x0 = (z(lev)-zz(i-1))/(zz(i)+zz(i-1))
2214 zair(lev) = air(i-1)*exp(-x0/hlocal)
2215 endif
2216 vcld(lev) = cloud(i)
2217 ! write(06,'('u:z,t,air ',i3,3e12.4)')
2218 ! lev,z(lev),zt(lev),zair(lev)
2219 END DO
2220
2221 10 CONTINUE
2222
2223 IF (i==nn) GO TO 20
2224 ! if(lev.ne.1) vcld(lev)=cloud(i)
2225 dzt = (zz(i+1)-zz(i))*.5
2226 ! number of layers depents on optical depth
2227 idt = max(ifix(cloud(i)*dzt*fnum),1)
2228 dzo = dzt/float(idt)
2229
2230 DO il = 1, idt
2231 lev = lev + 1
2232
2233 IF (lev>nj) THEN
2234 CALL wrf_error_fatal ( 'LEV > NJ, NJ GROESSER WAEHLEN')
2235 END IF
2236
2237 z(lev) = z(lev-1) + dzo
2238 zt(lev) = t(i) + (z(lev)-zz(i))/(zz(i+1)-zz(i))*(t(i+1)-t(i))
2239 if(abs(air(i)-air(i+1)).lt.air(i)/1.e5)then
2240 zair(lev) = air(i)
2241 else
2242 hlocal = 1./alog(air(i)/air(i+1))
2243 x0 = (z(lev)-zz(i))/(zz(i+1)+zz(i))
2244 zair(lev) = air(i)*exp(-x0/hlocal)
2245 endif
2246 vcld(lev-1) = cloud(i)
2247 ! write(06,'('o:z,t,air ',i3,3e12.4)')
2248 ! lev,z(lev),zt(lev),zair(lev)
2249 END DO
2250
2251 20 CONTINUE
2252 END IF
2253 ! write(06,'('o:z,t,air ',i3,3e12.4)') lev,z(lev),zt(lev),zair(lev)
2254 30 CONTINUE
2255
2256 ! number of levels including additional cloud levels
2257
2258 nlevel = lev
2259
2260 IF (nlevel>nj) print *, ' NLEVEL > NJ, NJ GROESSER WAEHLEN ', &
2261 nlevel
2262
2263 ! write(06,'(' nlevel',i3)') nlevel
2264
2265
2266 ! assign default yields
2267 DO 40 nr = 1, nreakj
2268
2269 DO 40 kl = kl0, kl1
2270
2271 DO 40 lev = 1, nlevel
2272 40 qy(lev,kl,nr) = xqy(kl,nr)
2273 ! assign default absorption cross sections
2274 DO 50 kl = kl0, kl1
2275
2276 DO 50 lev = 1, nlevel
2277
2278 DO 50 nr = 1, nreakj
2279 50 s(lev,kl,nr) = xs(kl,nr)
2280 ! --------------------------------------------------------------------
2281 ! re-calculate altitude dependent quantum yields. this currently
2282 ! applies to
2283 ! 2=o3->o(1d)
2284 ! 12=ch2o->h2+co
2285 ! 13=ch3cho->ch3+cho
2286 ! 14=ch3coch3
2287 ! 15=ch3coch2ch3
2288 ! 16=hcocho -> 0.13 hcho + 1.87 co process a
2289 ! 17=ch3cocho
2290 ! 22=hcocho -> 0.45 hcho + 1.55 co + 0.80 ho2 process b
2291
2292 ! o3 and ketones yield is calculated from fit equations,
2293 ! while for ch3cho and the dicarbonyls yields are calculated
2294 ! from the ntp yield by linear adjustment to 1/yield. the ch2o yield
2295 ! recalculated only for wavelengths longer than 329 nm. the yields for
2296 ! 1=o3->o(3p) are calculated as (1.- singlet d yield).
2297
2298 DO 60 lev = 1, nlevel
2299 ! o3 ozone:
2300 tau = zt(lev) - 230.
2301 a = 0.9*(0.369+2.85E-4*tau+1.28E-5*tau*tau+2.57E-8*tau*tau*tau)
2302 b = -0.575 + 5.59E-3*tau - 1.439E-5*tau*tau - 3.27E-8*tau*tau*tau
2303 c = 0.9*(0.518+9.87E-4*tau-3.94E-5*tau*tau+3.91E-7*tau*tau*tau)
2304 xl0 = 308.20 + 4.4871E-2*tau + 6.9380E-5*tau*tau - &
2305 2.5452E-6*tau*tau*tau
2306
2307 DO 60 kl = kl0, kl1
2308 xl = wl(kl)
2309 qy(lev,kl,2) = a*atan(b*(xl-xl0)) + c
2310
2311 IF (qy(lev,kl,2)<0.) qy(lev,kl,2) = 0.0
2312
2313 IF (qy(lev,kl,2)>0.9) qy(lev,kl,2) = 0.9
2314 qy(lev,kl,3) = 1.0 - qy(lev,kl,2)
2315 ! ch2o formaldehyde:
2316 IF ((xl>=330.) .AND. qy(lev,kl,12)>0.) THEN
2317 phi1 = qy(lev,kl,11)
2318 phi2 = qy(lev,kl,12)
2319 phi20 = 1. - phi1
2320 ak300 = ((1./phi2)-(1./phi20))/2.54E+19
2321 akt = ak300*(1.+61.69*(1.-zt(lev)/300.)*(xl/329.-1.))
2322 qy(lev,kl,12) = 1./((1./phi20)+zair(lev)*akt)
2323 END IF
2324
2325 IF (qy(lev,kl,12)>1.) qy(lev,kl,12) = 1.0
2326
2327 IF (qy(lev,kl,12)<0.) qy(lev,kl,12) = 0.0
2328 ! ch3cho acetaldehyde:
2329 IF (xqy(kl,13)/=0.) THEN
2330 qy(lev,kl,13) = 1./(1.+(1./xqy(kl,13)-1.)*zair(lev)/2.465E19)
2331 END IF
2332 ! ch3coch3 acetone:
2333 qy(lev,kl,14) = 0.0766 + 0.09415*exp(-zair(lev)/3.222E18)
2334 ! ch3coch2ch3 methyl ethyl ketone:
2335 qy(lev,kl,15) = qy(lev,kl,14)
2336 ! hcocho glyoxal process a:
2337 IF (xqy(kl,16)/=0.) THEN
2338 qy(lev,kl,16) = 1./(1.+(1./xqy(kl,16)-1.)*zair(lev)/2.465E19)
2339 END IF
2340 ! ch3cocho methylglyoxal:
2341 IF (xqy(kl,17)/=0.) THEN
2342 qy(lev,kl,17) = 1./(1.+(1./xqy(kl,17)-1.)*zair(lev)/2.465E19)
2343 END IF
2344 ! hcocho glyoxal process b:
2345 IF (xqy(kl,22)/=0.) THEN
2346 qy(lev,kl,22) = 1./(1.+(1./xqy(kl,22)-1.)*zair(lev)/2.465E19)
2347 END IF
2348
2349 60 CONTINUE
2350 ! _______________________________________________________________________
2351 ! correct absorption cross sections for t and p dep. for now, do
2352 ! 2=ozone
2353 DO 90 kl = kl0, kl1
2354
2355 DO 80 lev = 1, nlevel
2356
2357 IF (kl<33 .OR. kl>61) GO TO 70
2358 tdiffx = zt(lev) - 230.
2359 s(lev,kl,2) = (so3tx(kl,1)+so3tx(kl,2)*tdiffx+so3tx(kl,3)*tdiffx* &
2360 tdiffx)*1.0E-18
2361 s(lev,kl,3) = s(lev,kl,2)
2362 70 CONTINUE
2363 80 CONTINUE
2364 90 CONTINUE
2365
2366 ! ----------------------------------------------* layers
2367 nlayer = nlevel - 1
2368 ! write(06,'(' Layers ',i3)') nlayer
2369 lay = 0
2370
2371 DO 100 i = 1, nlayer
2372 lay = lay + 1
2373 dz = z(i+1) - z(i)
2374 zmid(lay) = z(i) + 0.5*dz
2375 vt(lay) = (zt(lay+1)+zt(lay))/2.
2376 vair(lay) = dz*1.E5*(zair(i+1)+zair(i))/2.
2377 100 CONTINUE
2378
2379 ! vo3(lay) = dz*1.e5*(o3(i+1) + o3(i))/2. ! umr. dz in cm
2380 ! vcld(lay) = 0. *dz
2381 ! vaer(lay) = (aer(i+1)-aer(i))/alog(aer(i+1)/aer(i)) *dz ! bei
2382
2383 CALL trapez(zz,o3,1,nn,zmid,vo3,1,nlayer,nn,nn,nj,nj)
2384
2385
2386 CALL trapez(zz,aer,1,nn,zmid,vaer,1,nlayer,nn,nn,nj,nj)
2387
2388 DO 110 i = 1, nlayer
2389 lay = i
2390 dz = z(i+1) - z(i)
2391 ! write(06,'('layer ',i3,6e12.4)') lay,zmid(lay),vt(lay),
2392 ! & vair(lay)/dz,vo3(lay),vaer(lay),vcld(lay)
2393 dz = z(i+1) - z(i)
2394
2395 ! umr. dz in cm
2396
2397 vo3(i) = dz*1.E5*vo3(i)
2398 vcld(i) = vcld(i)*dz
2399 110 vaer(i) = vaer(i)*dz
2400
2401
2402 ! normalize aerosol optical depth to unity sum
2403 sum = 0.
2404
2405 DO 120 lay = 1, nlayer
2406 sum = sum + vaer(lay)
2407 120 CONTINUE
2408
2409 DO 130 lay = 1, nlayer
2410 vaer(lay) = vaer(lay)/sum
2411 130 CONTINUE
2412
2413 ! calculated vertical column of o2 above the midpoint of each layer:
2414 ! want to use this for computing the average schumann-runge cross
2415 ! secti
2416 ! in each layer.
2417 ! so use half of current layer and half of previous higher layer
2418
2419
2420 cvo2(nlayer) = 0.2095*vair(nlayer)/2.
2421
2422 DO 140 ii = 2, nlayer
2423 lay = nlayer - ii + 1
2424 cvo2(lay) = cvo2(lay+1) + 0.2095*(vair(lay)+vair(lay+1))/2.
2425 140 CONTINUE
2426
2427 ! correct attenuation coefficients for pressure and/or temperature
2428 ! dep.
2429 ! for now do only ozone absorption.
2430 DO 160 kl = kl0, kl1
2431
2432 DO 150 lay = 1, nlayer
2433 tdiffx = vt(lay) - 230.
2434 ao3(lay,kl) = xs(kl,2)
2435
2436 IF (kl>=33 .AND. kl<=61) ao3(lay,kl) = (so3tx(kl,1)+so3tx(kl,2)* &
2437 tdiffx+so3tx(kl,3)*tdiffx*tdiffx)*1.0E-18
2438
2439 150 CONTINUE
2440 160 CONTINUE
2441
2442 ! write(06,'(' z ')')
2443 ! write(06,'(5e12.4 )') z
2444 ! write(06,'(' zmid ')')
2445 ! write(06,'(5e12.4 )') zmid
2446 ! write(06,'(' zt ')')
2447 ! write(06,'(5e12.4 )') zt
2448 ! write(06,'(' vt ')')
2449 ! write(06,'(5e12.4 )') vt
2450 ! write(06,'(' vo3 ')')
2451 ! write(06,'(5e12.4 )') vo3
2452 ! write(06,'(' vair ')')
2453 ! write(06,'(5e12.4 )') vair
2454 ! write(06,'(' vaer ')')
2455 ! write(06,'(5e12.4 )') vaer
2456
2457
2458 RETURN
2459
2460 END SUBROUTINE subgrid
2461
2462 SUBROUTINE optics(iprt,a,vair,arayl,gray,omray,ao2,vo3,ao3,vcld,gcld, &
2463 omcld,vaer,aaer,gaer,omaer,nlayer,nlevel,nsurf,endir,endn,enup)
2464 ! sm this subroutine prepares the data needed for the flux
2465 ! calculation,
2466 ! t
2467 ! sm calls the scattering subroutine delted. it returns values of
2468 ! the
2469 ! sm flux flux(lev,kl) for altitude lev-1, wavelength kl.
2470 ! sm it calculates the optical depths (vertical)
2471 ! sm for each layer, from the vertical profiles of o2, o3,
2472 ! sm air, cloud, and aerosol, and from the associated 'cross
2473 ! sections
2474 ! implicit real*4 (a-h,o-z)
2475 ! .. Scalar Arguments ..
2476 REAL :: a, gaer, gcld, gray, omaer, omcld, omray
2477 INTEGER :: iprt, nlayer, nlevel, nsurf
2478 ! .. Local Scalars ..
2479 REAL :: dtabs, dtaer, dtair, dtcld, dto2, dto3, dtscat, fdinc, fuinc, &
2480 solflx, sumtau
2481 INTEGER :: ii, lay, lev, nlev, nz, nzm1
2482 LOGICAL :: mudept
2483 ! .. Intrinsic Functions ..
2484 ! EXTERNAL delted
2485 INTRINSIC amin1
2486 ! .. Parameters ..
2487 INTEGER, PARAMETER :: nsol = 1
2488 ! .. Array Arguments ..
2489 REAL :: aaer(130), ao2(nj,130), ao3(nj,130), arayl(130), &
2490 endir(nj,130), endn(nj,130), enup(nj,130), vaer(nj), vair(nj), &
2491 vcld(nj), vo3(nj)
2492 ! .. Local Arrays ..
2493 REAL :: alb(nsol), dir(nj,nsol), dtau(nj), flxd(nj,nsol), &
2494 flxu(nj,nsol), g(nj), musun(nsol), om(nj)
2495 ! if(iprt.eq.1)then
2496 ! write(06,*) a,kl0,kl1,gray,omray,gcld,omcld,gaer,omaer,
2497 ! &nlayer,nlevel,nsurf
2498 ! do kl=1,nlevel
2499 ! write(06,'(5e12.4)') vair(kl),vo3(kl),vcld(kl),vaer(kl)
2500 ! write(06,*) ao2
2501 ! enddo
2502 ! endif
2503 ! stop
2504 ! loop over wavelengths
2505 DO 30 kl = kl0, kl1
2506
2507 ! calculate optical depths for all layers (including cloud
2508 ! sublayers)
2509
2510 sumtau = 0.
2511
2512 DO 10 lay = nsurf, nlayer
2513 ii = nlayer + 1 - lay
2514 dtair = vair(lay)*arayl(kl)
2515 dto2 = 0.2095*vair(lay)*ao2(lay,kl)
2516 dto3 = vo3(lay)*ao3(lay,kl)
2517 dtcld = vcld(lay)
2518 dtaer = vaer(lay)*aaer(kl)
2519
2520 dtscat = dtair + dtcld + dtaer
2521 dtabs = dto2 + dto3
2522
2523 dtau(ii) = dtabs + dtscat
2524 om(ii) = (omray*dtair+omcld*dtcld+omaer*dtaer)/dtau(ii)
2525 g(ii) = (gray*dtair+gcld*dtcld+gaer*dtaer)/dtscat
2526 sumtau = sumtau + dtau(ii)
2527 ! if(kl.eq.103.and.iprt.eq.1) ! 55um
2528 ! & write(06,'(' Ope',i4,3e12.4)')
2529 ! 2 ii,arayl(kl),dtcld,dtau(ii)
2530 10 CONTINUE
2531
2532 ! if(kl.eq.103.and.iprt.eq.1) ! 55um
2533 ! &write(06,'(' Opt. Dicke',i4,2e12.4)') kl,sumtau
2534
2535 ! initialize fluxes and other delted parameters
2536 solflx = 1.
2537 fdinc = 0.
2538 fuinc = 0.
2539 alb(1) = albedoph(kl)
2540 musun(1) = amin1(a,0.999)
2541 mudept = .FALSE.
2542 nz = nlevel - nsurf + 1
2543 nzm1 = nz - 1
2544 nlev = nj
2545 ! if(iprt.eq.1)print *,'nz = ',nz
2546 ! ---------------------------------------------------------------------
2547 CALL delted(dtau,om,g,musun,alb,solflx,fdinc,fuinc,mudept,nz,dir, &
2548 flxu,flxd)
2549 ! if(iprt.eq.1)print *,'nz = ',nz
2550 ! ---------------------------------------------------------------------
2551 ! return to upright grid
2552 DO 20 ii = 1, nlevel - nsurf + 1
2553 lev = nlevel + 1 - ii
2554 endir(lev,kl) = dir(ii,nsol)
2555 endn(lev,kl) = flxd(ii,nsol)
2556 enup(lev,kl) = flxu(ii,nsol)
2557 20 CONTINUE
2558
2559 30 CONTINUE
2560
2561 RETURN
2562
2563 END SUBROUTINE optics
2564
2565 ! #######################################################################
2566
2567 SUBROUTINE nlayde(ib,musun,alb,fdinc,fuinc,mudept,dir,flxu,flxd,tau,lm, &
2568 pp,ex,tx,ty,tz,isav,flxsun,nl2)
2569 ! multi-layer delta-eddington
2570 ! ib = number of levels
2571 ! the top and bottom boundary conditions plus the flux continuity
2572 ! conditions at each interior level form a penta-diagonal system
2573 ! of 2*ib-2 equations for the unknown constants (2 for each layer).
2574 ! the columns of the -ss- array contain the diagonals of the coeffi
2575 ! cient matrix, the lowermost diagonal in column 1, etc. (this is t
2576 ! so-called band storage mode required by imsl routine leqt2b).
2577 ! ueberflussig, falls man nur eine Sonnenhoehe vorsieht.
2578 ! .. Scalar Arguments ..
2579 REAL :: fdinc, fuinc
2580 INTEGER :: ib, nl2
2581 LOGICAL :: mudept
2582 ! .. Local Scalars ..
2583 REAL :: albdo, rmu0, t1
2584 INTEGER :: i, ibm1, ic, ijob, ip1, ir, j, ktr, last, lastm2, np
2585 ! .. Intrinsic Functions ..
2586 ! EXTERNAL leqt1b
2587 INTRINSIC exp
2588 ! .. Parameters ..
2589 INTEGER, PARAMETER :: nsol = 1
2590 ! .. Array Arguments ..
2591 REAL :: alb(1), dir(nj,1), ex(nj), flxd(nj,1), flxsun(1), flxu(nj,1), &
2592 isav(nj), lm(nj), musun(1), pp(nj), tau(nj), tx(nj), ty(nj), tz(nj)
2593 ! .. Local Arrays ..
2594 REAL :: alph(nj), beta(nj), cc(mj,5), exsun(nj), ss(mj,5), work(mj,3), &
2595 x(mj)
2596
2597 ibm1 = ib - 1
2598 last = 2*ib - 2
2599 ss(1,5) = 0.
2600 ss(1,1) = ss(1,5)
2601 ss(1,2) = ss(1,5)
2602 ss(1,3) = (1.-pp(1))/ex(1)
2603 ss(1,4) = (1.+pp(1))*ex(1)
2604 lastm2 = last - 2
2605
2606 DO 10 j = 2, lastm2, 2
2607 i = j/2
2608 ip1 = i + 1
2609 ss(j,1) = 0.
2610 ss(j,2) = 1.0
2611 ss(j,3) = 1.0
2612 ss(j,4) = -1.0/ex(ip1)
2613 ss(j,5) = -ex(ip1)
2614 ss(j+1,1) = -pp(i)
2615 ss(j+1,2) = pp(i)
2616 ss(j+1,3) = pp(ip1)/ex(ip1)
2617 ss(j+1,4) = -pp(ip1)*ex(ip1)
2618 ss(j+1,5) = 0.
2619 10 CONTINUE
2620 ss(last,5) = 0.
2621 ss(last,1) = ss(last,5)
2622 ss(last,4) = ss(last,5)
2623
2624 IF (mudept) GO TO 30
2625 ss(last,2) = 1. + pp(ibm1) - alb(1)*(1.-pp(ibm1))
2626 ss(last,3) = 1. - pp(ibm1) - alb(1)*(1.+pp(ibm1))
2627
2628 ! calculate the l-u decomposition of penta-diagonal matrix -ss-
2629
2630 ! leqt2b call for testing purposes
2631 ! call leqt2b(ss,last,2,2,nl2,x,1,nl2, 1 ,work,nl2,work(1,8))
2632
2633 ! leqt1b destroys the input coeff matrix, so since we must
2634 ! preserve -ss-, we must let it destroy -cc- instead.
2635 DO 20 ic = 1, 5
2636
2637 DO 20 ir = 1, last
2638 20 cc(ir,ic) = ss(ir,ic)
2639 ! ---------------------------------------------------------------------
2640 CALL leqt1b(cc,last,2,2,nl2,x,1,nl2,1,work)
2641 ! ---------------------------------------------------------------------
2642
2643 ! for each sun angle, calculate the r.h.s. of the banded system,
2644 ! solve, and use the solution to construct the fluxes at each level
2645
2646 30 DO 90 np = 1, nsol
2647 rmu0 = 1./musun(np)
2648 t1 = rmu0**2 - lm(1)**2
2649
2650 IF (t1==0.) THEN
2651 t1 = 1.E-7
2652 PRINT *, 'ACHTUNG t1=0 fuer lm(1)'
2653 END IF
2654
2655 alph(1) = tx(1)/t1
2656 beta(1) = ty(1)*(musun(np)*tz(1)+rmu0)/t1
2657 x(1) = alph(1) + beta(1) + fdinc
2658
2659 DO 40 j = 2, lastm2, 2
2660 i = j/2
2661 ip1 = i + 1
2662 t1 = rmu0**2 - lm(ip1)**2
2663
2664 IF (t1==0.) THEN
2665 t1 = 1.E-7
2666 PRINT *, 'ACHTUNG t1=0 fuer lm(', ip1, ')'
2667 END IF
2668
2669 alph(ip1) = tx(ip1)/t1
2670 beta(ip1) = ty(ip1)*(musun(np)*tz(ip1)+rmu0)/t1
2671 exsun(ip1) = exp(-rmu0*tau(ip1))
2672 x(j) = (alph(i)-alph(ip1))*exsun(ip1)
2673 x(j+1) = (beta(i)-beta(ip1))*exsun(ip1)
2674 40 CONTINUE
2675 exsun(ib) = exp(-rmu0*tau(ib))
2676
2677 IF (mudept) GO TO 50
2678 albdo = alb(1)
2679 ijob = 2
2680 GO TO 70
2681 50 albdo = alb(np)
2682 ijob = 0
2683 ss(last,2) = 1. + pp(ibm1) - alb(np)*(1.-pp(ibm1))
2684 ss(last,3) = 1. - pp(ibm1) - alb(np)*(1.+pp(ibm1))
2685
2686 DO 60 ic = 1, 5
2687
2688 DO 60 ir = 1, last
2689 60 cc(ir,ic) = ss(ir,ic)
2690
2691 70 x(last) = (alph(ibm1)-beta(ibm1)+albdo*(flxsun(np)-alph(ibm1)-beta( &
2692 ibm1)))*exsun(ib) + fuinc
2693
2694 ! solve penta-diagonal system with r.h.s. -x-. soln goes into -x-
2695
2696 ! call leqt2b(ss,last,2,2,nl2,x,1,nl2,ijob,work,nl2,work(1,8))
2697
2698 ! ---------------------------------------------------------------------
2699 CALL leqt1b(cc,last,2,2,nl2,x,1,nl2,ijob,work)
2700 ! ---------------------------------------------------------------------
2701
2702 dir(1,np) = flxsun(np)
2703 flxd(1,np) = fdinc
2704 flxu(1,np) = (1.+pp(1))/ex(1)*x(1) + (1.-pp(1))*ex(1)*x(2) - &
2705 alph(1) + beta(1)
2706 ktr = 2
2707
2708 DO 80 i = 1, ibm1
2709
2710 IF (i+1/=isav(ktr)) GO TO 80
2711 dir(ktr,np) = flxsun(np)*exsun(i+1)
2712 flxd(ktr,np) = (1.-pp(i))*x(2*i-1) + (1.+pp(i))*x(2*i) - &
2713 (alph(i)+beta(i))*exsun(i+1)
2714 flxu(ktr,np) = (1.+pp(i))*x(2*i-1) + (1.-pp(i))*x(2*i) - &
2715 (alph(i)-beta(i))*exsun(i+1)
2716 ktr = ktr + 1
2717 80 CONTINUE
2718 90 CONTINUE
2719
2720 RETURN
2721
2722 END SUBROUTINE nlayde
2723
2724 ! #######################################################################
2725
2726 SUBROUTINE delted(dtau,om,g,musun,alb,solflx,fdinc,fuinc,mudept,nz,dir, &
2727 flxu,flxd)
2728 ! calculate up- and down-fluxes of radiation in a vertically inhomo-
2729 ! geneous atmosphere using the delta-eddington approximation
2730 ! author-- w.j. wiscombe
2731 ! national center for atmospheric research
2732 ! p.o. box 3000
2733 ! boulder, colorado 80303
2734 ! input variables
2735 ! nz = number of levels (level 1 is the top of the atmosphere,
2736 ! level nz is the surface)
2737 ! dtau(i), i=1,...,nz-1 = optical depth of layer between levels i an
2738 ! om(i), i=1,...,nz-1 = single-scattering albedoph of layer between
2739 ! levels i and i+1
2740 ! g(i), i=1,...,nz-1 = asymmetry factor for layer between levels i a
2741 ! nsol = number of incident-beam zenith angles
2742 ! musun(i),i=1,...,nsol = cosine(s) of incident-beam zenith angle(s)
2743 ! alb(i), i=1,...,nsol = surface albedoph
2744 ! mudept = true, alb(i) corresponds to musun(i). false, alb(1) is
2745 ! used for all values of musun(i).
2746 ! solflx = incident-beam flux (normal to beam) at level 1
2747 ! the beam) at the top of the atmosphere
2748 ! fdinc = incident diffuse down-flux at level 1
2749 ! fuinc = incident diffuse up-flux at level nz
2750 ! nlev = level dimension (of arrays dtau, etc.)
2751 ! output variables (in same units as solflx, fdinc, and fuinc)
2752 ! dir(i,np) direct flux at level -i- for sun angle -np-
2753 ! (note--in the delta-eddington approxn, because of t
2754 ! truncation of the forward scattering peak, this
2755 ! quantity includes scattered radiation travelling in
2756 ! very nearly the same direction as the actual direct
2757 ! flux. e.g., it includes the aureole around the sun
2758 ! flxd(i,np) diffuse down-flux at level -i- for sun angle -np-
2759 ! (note--this will be less than the actual diffuse
2760 ! down-flux by the same amount that the direct flux
2761 ! -dir- is augmented.)
2762 ! flxu(i,np) diffuse up-flux at level -i- for sun angle -np-
2763 ! internal
2764 ! code variable description (or name in write-up)
2765 ! gp g-prime (transformed asymmetry parameter)
2766 ! omp omega-prime (transformed single scattering albedoph)
2767 ! dtaup delta-tau-prime (transformed layer optical depth)
2768 ! tau(i) cumulative optical depth from top (i=1) to level i
2769 ! lm(i) lambda-sub-i
2770 ! pp(i) p-sub-i
2771 ! lmdtau lm(i)*dtaup
2772 ! ex(i) exp(lmdtau)
2773 ! exsun(i) exp(-tau(i)/musun(np))
2774 ! alph(i) alpha-sub-i
2775 ! beta(i) beta-sub-i
2776 ! tx(i) 0.75*solflx*omp*(1.+gp *(1.-omp))
2777 ! ty(i) 0.5*solflx*omp
2778 ! tz(i) 3.*gp*(1.-omp)
2779 ! (tx,ty,tz are merely intermediate quantities for computing alph,
2780 ! isav array of level indices. fluxes are calculated only
2781 ! at these levels.
2782 ! flxsun(np) incident flux musun(np)*solflx at level 1
2783 ! prec a number somewhat larger than the computer precisio
2784 ! subtracted from any single-scattering albedophs which
2785 ! are exactly equal to one.
2786 ! cutpt any layer for which lmdtau.gt.cutpt is subdivided
2787 ! into equal sublayers, all of which have lmdtau.lt.c
2788 ! nsub number of sublayers into which an offending layer
2789 ! is divided (the whole process being transparent
2790 ! to the user)
2791 ! nl2 2*nlev-2 (input to banded matrix subroutine)
2792 ! ss(nl2,5) the penta-diagonal matrix c, in band storage mode
2793 ! cc(nl2,5) same as -ss- array. used to submit -ss- to leqt1b.
2794 ! work(nl2,3) a temporary storage array used by subroutine leqt1
2795 ! x(nl2) the vector (x-hat). also temporarily stores
2796 ! r.h.s. d in linear system c*(x-hat) = d.
2797 ! --note-- this code is not perfectly optimized, either in terms of
2798 ! core storage or execution speed, but it should be noted t
2799 ! in general, the lions share of computing time is occupied
2800 ! the exponentials and the penta-diagonal solution routines
2801 ! so eliminating a few operations here or there has almost
2802 ! no impact on running time.
2803 ! ******************* 10 x computer precision **********************
2804 ! *********** cut-off point for lm(i)*dtaup ************
2805 ! .. Scalar Arguments ..
2806 REAL :: fdinc, fuinc, solflx
2807 INTEGER :: nz
2808 LOGICAL :: mudept
2809 ! .. Local Scalars ..
2810 REAL :: aux, aux2, c1, cutpt, dtaup, ff, gp, lmdtau, omp, prec, scale, &
2811 t1
2812 INTEGER :: i, ii, ip1, iup, iupm1, ktr, layers, nl2, nlev, np, nsub, &
2813 nzm1
2814 ! .. Intrinsic Functions ..
2815 ! EXTERNAL nlayde
2816 INTRINSIC exp, float, sqrt
2817 ! .. Parameters ..
2818 INTEGER, PARAMETER :: nsol = 1
2819 ! .. Array Arguments ..
2820 REAL :: alb(1), dir(nj,1), dtau(nj), flxd(nj,1), flxu(nj,1), g(nj), &
2821 musun(1), om(nj)
2822 ! .. Local Arrays ..
2823 REAL :: ex(nj), flxsun(1), isav(nj), lm(nj), pp(nj), tau(nj), tx(nj), &
2824 ty(nj), tz(nj)
2825 ! .. Data Statements ..
2826 DATA c1/0.66666666666667/
2827 DATA prec/1.E-7/
2828 DATA cutpt/7./
2829 ! set incident flux at top of atmosphere
2830 DO 10 np = 1, nsol
2831 10 flxsun(np) = musun(np)*solflx
2832
2833 nzm1 = nz - 1
2834 nlev = nj
2835 nl2 = 2*nlev - 2
2836
2837 ! scale optical depth, sing-scat albedoph, and asymmetry factor
2838 ! and calculate various fcns of these variables
2839
2840 nzm1 = nz - 1
2841 tau(1) = 0.
2842
2843 DO 20 i = 1, nzm1
2844 ff = g(i)**2
2845 gp = g(i)/(1.+g(i))
2846 scale = 1. - ff*om(i)
2847 omp = (1.-ff)*om(i)/scale
2848
2849 IF (om(i)==1.0) omp = 1. - prec
2850 t1 = 1. - omp*gp
2851 aux = 3.*(1.-omp)*t1
2852
2853 IF (om(i)==1.0) aux = 3.*prec*t1
2854 aux2 = 1. - om(i)
2855 lm(i) = sqrt(aux)
2856 pp(i) = c1*lm(i)/t1
2857 t1 = gp*(1.-omp)
2858 tx(i) = 0.75*solflx*omp*(1.+t1)
2859 ty(i) = 0.5*solflx*omp
2860 tz(i) = 3.*t1
2861 isav(i) = i
2862 dtaup = scale*dtau(i)
2863 lmdtau = lm(i)*dtaup
2864 ! test for a layer which is so highly absorbing that it would
2865 ! cause ill-conditioning in the penta-diagonal matrix. if one
2866 ! is found, subdivide it appropriately.
2867 IF (lmdtau>cutpt) GO TO 30
2868 ex(i) = exp(lmdtau)
2869 tau(i+1) = tau(i) + dtaup
2870 20 CONTINUE
2871 isav(nz) = nz
2872
2873
2874 ! normal calculation
2875
2876 ! ---------------------------------------------------------------------
2877 ! Nur noch nz > 2 moeglich
2878 CALL nlayde(nz,musun,alb,fdinc,fuinc,mudept,dir,flxu,flxd,tau,lm,pp, &
2879 ex,tx,ty,tz,isav,flxsun,nl2)
2880 ! ---------------------------------------------------------------------
2881 GO TO 80
2882
2883 ! sidestep potential ill-conditioning by subdividing offending
2884 ! layer, and any others like it.
2885
2886 30 layers = nzm1
2887 ktr = i
2888
2889 40 nsub = lmdtau/cutpt + 1.
2890 dtaup = dtaup/float(nsub)
2891 ex(i) = exp(lm(i)*dtaup)
2892 tau(i+1) = tau(i) + dtaup
2893 ip1 = i + 1
2894 iup = i + nsub
2895 iupm1 = iup - 1
2896
2897 DO 50 ii = ip1, iupm1
2898 lm(ii) = lm(i)
2899 pp(ii) = pp(i)
2900 tx(ii) = tx(i)
2901 ty(ii) = ty(i)
2902 tz(ii) = tz(i)
2903 ex(ii) = ex(i)
2904 50 tau(ii+1) = tau(ii) + dtaup
2905 ktr = ktr + 1
2906 isav(ktr) = isav(ktr-1) + nsub
2907 layers = layers + nsub - 1
2908
2909 IF (layers>nlev) then
2910 CALL wrf_error_fatal ( 'layers>nlev')
2911 endif
2912
2913 IF (iup>layers) GO TO 70
2914
2915 DO 60 i = iup, layers
2916 ff = g(ktr)**2
2917 gp = g(ktr)/(1.+g(ktr))
2918 scale = 1. - ff*om(ktr)
2919 omp = (1.-ff)*om(ktr)/scale
2920
2921 IF (om(ktr)==1.0) omp = 1. - prec
2922 t1 = 1. - omp*gp
2923 lm(i) = sqrt(3.*(1.-omp)*t1)
2924 pp(i) = c1*lm(i)/t1
2925 t1 = gp*(1.-omp)
2926 tx(i) = 0.75*solflx*omp*(1.+t1)
2927 ty(i) = 0.5*solflx*omp
2928 tz(i) = 3.*t1
2929 dtaup = scale*dtau(ktr)
2930 lmdtau = lm(i)*dtaup
2931 ! test for a layer which is so highly absorbing that it would
2932 ! cause ill-conditioning in the penta-diagonal matrix. if one
2933 ! is found, subdivide it appropriately.
2934 IF (lmdtau>cutpt) GO TO 40
2935 ex(i) = exp(lmdtau)
2936 tau(i+1) = tau(i) + dtaup
2937 ktr = ktr + 1
2938 isav(ktr) = isav(ktr-1) + 1
2939 60 CONTINUE
2940
2941 ! ---------------------------------------------------------------------
2942 70 CALL nlayde(layers+1,musun,alb,fdinc,fuinc,mudept,dir,flxu,flxd,tau, &
2943 lm,pp,ex,tx,ty,tz,isav,flxsun,nl2)
2944 ! ---------------------------------------------------------------------
2945
2946 80 CONTINUE
2947 RETURN
2948
2949 END SUBROUTINE delted
2950
2951 ! #######################################################################
2952
2953 SUBROUTINE calc_zenith(lat,long,ijd,gmt,azimuth,zenith)
2954 ! this subroutine calculates solar zenith and azimuth angles for a
2955 ! part
2956 ! time and location. must specify:
2957 ! input:
2958 ! lat - latitude in decimal degrees
2959 ! long - longitude in decimal degrees
2960 ! gmt - greenwich mean time - decimal military eg.
2961 ! 22.75 = 45 min after ten pm gmt
2962 ! output
2963 ! zenith
2964 ! azimuth
2965 ! .. Scalar Arguments ..
2966 REAL :: azimuth, gmt, lat, long, zenith
2967 INTEGER :: ijd
2968 ! .. Local Scalars ..
2969 REAL :: caz, csz, cw, d, decl, dr, ec, epsi, eqt, eyt, feqt, feqt1, &
2970 feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
2971 pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, &
2972 yt, zpt, zr
2973 INTEGER :: jd
2974 ! .. Intrinsic Functions ..
2975 INTRINSIC acos, atan, cos, min, sin, tan
2976 ! convert to radians
2977 pi = 3.1415926535590
2978 dr = pi/180.
2979 rlt = lat*dr
2980 rphi = long*dr
2981
2982 ! print julian days current 'ijd'
2983
2984 ! ???? + (yr - yref)
2985
2986 jd = ijd
2987
2988
2989
2990 d = jd + gmt/24.0
2991 ! calc geom mean longitude
2992 ml = 279.2801988 + .9856473354*d + 2.267E-13*d*d
2993 rml = ml*dr
2994
2995 ! calc equation of time in sec
2996 ! w = mean long of perigee
2997 ! e = eccentricity
2998 ! epsi = mean obliquity of ecliptic
2999 w = 282.4932328 + 4.70684E-5*d + 3.39E-13*d*d
3000 wr = w*dr
3001 ec = 1.6720041E-2 - 1.1444E-9*d - 9.4E-17*d*d
3002 epsi = 23.44266511 - 3.5626E-7*d - 1.23E-15*d*d
3003 pepsi = epsi*dr
3004 yt = (tan(pepsi/2.0))**2
3005 cw = cos(wr)
3006 sw = sin(wr)
3007 ssw = sin(2.0*wr)
3008 eyt = 2.*ec*yt
3009 feqt1 = sin(rml)*(-eyt*cw-2.*ec*cw)
3010 feqt2 = cos(rml)*(2.*ec*sw-eyt*sw)
3011 feqt3 = sin(2.*rml)*(yt-(5.*ec**2/4.)*(cw**2-sw**2))
3012 feqt4 = cos(2.*rml)*(5.*ec**2*ssw/4.)
3013 feqt5 = sin(3.*rml)*(eyt*cw)
3014 feqt6 = cos(3.*rml)*(-eyt*sw)
3015 feqt7 = -sin(4.*rml)*(.5*yt**2)
3016 feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
3017 eqt = feqt*13751.0
3018
3019 ! convert eq of time from sec to deg
3020 reqt = eqt/240.
3021 ! calc right ascension in rads
3022 ra = ml - reqt
3023 rra = ra*dr
3024 ! calc declination in rads, deg
3025 tab = 0.43360*sin(rra)
3026 rdecl = atan(tab)
3027 decl = rdecl/dr
3028 ! calc local hour angle
3029 lbgmt = 12.0 - eqt/3600. + long*24./360.
3030 lzgmt = 15.0*(gmt-lbgmt)
3031 zpt = lzgmt*dr
3032 csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
3033 if(csz.gt.1)print *,'calczen,csz ',csz
3034 csz = min(1.,csz)
3035 ! zr = acos(csz)
3036 ! zenith = zr/dr
3037 zr = acos(csz)
3038 zenith = zr/dr
3039 ! calc local solar azimuth
3040 caz = (sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr))
3041 if(caz.lt.-0.999999)then
3042 azimuth=180.
3043 elseif(caz.gt.0.999999)then
3044 azimuth=0.
3045 else
3046 raz = acos(caz)
3047 azimuth = raz/dr
3048 endif
3049 ! caz = min(1.,(sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr)))
3050 ! if(caz.lt.-1)print *,'calczen ',caz
3051 ! caz = max(-1.,caz)
3052 ! raz = acos(caz)
3053 ! azimuth = raz/dr
3054
3055 IF (lzgmt>0) azimuth = azimuth + (2*(180.-azimuth))
3056 ! 200 format(' ',f7.2,2(12x,f7.2))
3057 RETURN
3058
3059 END SUBROUTINE calc_zenith
3060
3061 SUBROUTINE trapez(x,y,ianfa,ienda,u,v,ianfn,iendn,ix,iy,iu,iv)
3062 ! * implemented 1992 by ansgar ruggaber, university of munich, frg
3063 ! * funded by the german minister of research and technology (bmft)
3064 ! * under contract no. 521-4007-07eu-738 8
3065 ! lineare interpolation of referencefield (x(i),y(i))
3066 ! to (u(i),v(i))
3067 ! save
3068 ! .. Scalar Arguments ..
3069 INTEGER :: ianfa, ianfn, ienda, iendn, iu, iv, ix, iy
3070 ! .. Array Arguments ..
3071 REAL :: u(iu), v(iv), x(ix), y(iy)
3072 ! .. Local Scalars ..
3073 REAL :: uumord, vumord, xumord, yumord
3074 INTEGER :: i, ianf, ianfnn, idrehu, idrehx, iendnn, iordu, iordx, j
3075
3076 idrehx = 0
3077
3078 IF (x(ianfa)>=x(ienda)) THEN
3079 ! das x-feld wird ansteigend geordnet, das y-feld entsprechend
3080 iordx = (ienda-ianfa+1)/2
3081
3082 DO i = ianfa, iordx
3083 xumord = x(i)
3084 x(i) = x(ienda+1-i)
3085 x(ienda+1-i) = xumord
3086 yumord = y(i)
3087 y(i) = y(ienda+1-i)
3088 y(ienda+1-i) = yumord
3089 END DO
3090
3091 idrehx = 1
3092 END IF
3093
3094 idrehu = 0
3095
3096 IF (u(ianfn)>=u(iendn)) THEN
3097 ! u-field increasing
3098 iordu = (iendn-ianfn+1)/2
3099
3100 DO i = ianfn, iordu
3101 uumord = u(i)
3102 u(i) = u(iendn+1-i)
3103 u(iendn+1-i) = uumord
3104 END DO
3105
3106 idrehu = 1
3107 END IF
3108
3109 ianfnn = ianfn
3110 10 CONTINUE
3111
3112 IF (u(ianfnn)<x(ianfa)) THEN
3113 ! no extrapolation at x(ianfa)
3114 v(ianfnn) = 1.0E-12
3115 ianfnn = ianfnn + 1
3116 GO TO 10
3117 END IF
3118
3119 iendnn = iendn
3120 20 CONTINUE
3121
3122 IF ((1.-1.0E-10)*u(iendnn)>x(ienda)) THEN
3123 ! no extrapolation at x(ienda)
3124 v(iendnn) = 1.0E-12
3125 iendnn = iendnn - 1
3126 GO TO 20
3127 END IF
3128
3129 ianf = ianfa
3130
3131 DO j = ianfnn, iendnn
3132
3133 DO i = ianf, ienda
3134
3135 IF (x(i)-u(j)) 30, 50, 40
3136 30 END DO
3137
3138 GO TO 70
3139 40 v(j) = y(i-1) + (y(i)-y(i-1))/(x(i)-x(i-1))*(u(j)-x(i-1))
3140 GO TO 60
3141 50 v(j) = y(i)
3142 60 ianf = i
3143 70 END DO
3144
3145 IF (idrehx/=0) THEN
3146 ! x- und y-field in starting position
3147 DO i = ianfa, iordx
3148 xumord = x(i)
3149 x(i) = x(ienda+1-i)
3150 x(ienda+1-i) = xumord
3151 yumord = y(i)
3152 y(i) = y(ienda+1-i)
3153 y(ienda+1-i) = yumord
3154 END DO
3155
3156 END IF
3157
3158 IF (idrehu/=0) THEN
3159
3160 DO i = ianfn, iordu
3161 uumord = u(i)
3162 u(i) = u(iendn+1-i)
3163 u(iendn+1-i) = uumord
3164 vumord = v(i)
3165 v(i) = v(iendn+1-i)
3166 v(iendn+1-i) = vumord
3167 END DO
3168
3169 END IF
3170
3171 END SUBROUTINE trapez
3172
3173 SUBROUTINE photmad_init(z_at_w,aerwrf,g,ids,ide,jds,jde,kds,kde,ims,ime, &
3174 jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
3175 ! local stuff
3176 ! .. Scalar Arguments ..
3177 REAL, INTENT (IN) :: g
3178 INTEGER, INTENT (IN) :: ide, ids, ime, ims, ite, its, jde, jds, jme, &
3179 jms, jte, jts, kde, kds, kme, kms, kte, kts
3180 ! .. Array Arguments ..
3181 REAL, INTENT (INOUT) :: aerwrf(ims:ime,kms:kme,jms:jme)
3182 REAL, INTENT (IN) :: z_at_w(ims:ime,kms:kme,jms:jme)
3183 ! .. Local Scalars ..
3184 REAL :: z1
3185 INTEGER :: i, j, k
3186 ! .. Local Arrays ..
3187 REAL :: aerext(kts:kte), phizz(kts:kte), z(kts:kte)
3188
3189 DO j = jts, jte
3190
3191 IF (j>jde-1) GO TO 20
3192
3193 DO i = its, ite
3194
3195 IF (i>ide-1) GO TO 10
3196
3197 ! z at w points
3198
3199 z1 = z_at_w(i,kts,j)
3200 z(kts)=0.
3201
3202 DO k = kts+1, kte
3203 z(k) = z_at_w(i,k,j) - z1
3204 END DO
3205
3206 DO k = kts, kte
3207 phizz(k) = .001*z(k)
3208 aerext(k) = 0.
3209 ! if(i.eq.its.and.j.eq.jts)print *,phizz(k),aerstd(k),kts,kte
3210 ! print *,phizz(k),kts,kte,ite,jte
3211 END DO
3212
3213 ! IF (phizz(kte-1)>20.) THEN
3214 ! CALL wrf_error_fatal ( 'phizz(kte-1)>20., set kl0 to 1')
3215 ! END IF
3216
3217 CALL trapez(zstd,aerstd,1,51,phizz,aerext,kts,kte,51,51,kte,kte)
3218
3219 DO k = kts, kte
3220 aerwrf(i,k,j) = aerext(k)
3221 ! if(i.eq.its)print *,k,i,j,aerext(k),phizz(k)
3222 ! print *,k,i,j,aerext(k),phizz(k)
3223 END DO
3224
3225 10 CONTINUE
3226 END DO
3227
3228 20 CONTINUE
3229 END DO
3230
3231 END SUBROUTINE photmad_init
3232
3233 END MODULE module_phot_mad