chgres_cube  1.6.0
 All Data Structures Files Functions Variables
wam_climo_data.f90
Go to the documentation of this file.
1 
12 
13 !-----------------------------------------------------------------------
29  module wam_gtd7bk_mod
30 ! msise-00 01-feb-02
31 !
32  real :: pt1(50)
33  real :: pt2(50)
34  real :: pt3(50)
35  real :: pa1(50)
36  real :: pa2(50)
37  real :: pa3(50)
38  real :: pb1(50)
39  real :: pb2(50)
40  real :: pb3(50)
41  real :: pc1(50)
42  real :: pc2(50)
43  real :: pc3(50)
44  real :: pd1(50)
45  real :: pd2(50)
46  real :: pd3(50)
47  real :: pe1(50)
48  real :: pe2(50)
49  real :: pe3(50)
50  real :: pf1(50)
51  real :: pf2(50)
52  real :: pf3(50)
53  real :: pg1(50)
54  real :: pg2(50)
55  real :: pg3(50)
56  real :: ph1(50)
57  real :: ph2(50)
58  real :: ph3(50)
59  real :: pi1(50)
60  real :: pi2(50)
61  real :: pi3(50)
62  real :: pj1(50)
63  real :: pj2(50)
64  real :: pj3(50)
65  real :: pk1(50)
66  real :: pl1(50)
67  real :: pl2(50)
68  real :: pm1(50)
69  real :: pm2(50)
70  real :: pn1(50)
71  real :: pn2(50)
72  real :: po1(50)
73  real :: po2(50)
74  real :: pp1(50)
75  real :: pp2(50)
76  real :: pq1(50)
77  real :: pq2(50)
78  real :: pr1(50)
79  real :: pr2(50)
80  real :: ps1(50)
81  real :: ps2(50)
82  real :: pu1(50)
83  real :: pu2(50)
84  real :: pv1(50)
85  real :: pv2(50)
86  real :: pw1(50)
87  real :: pw2(50)
88  real :: px1(50)
89  real :: px2(50)
90  real :: py1(50)
91  real :: py2(50)
92  real :: pz1(50)
93  real :: pz2(50)
94  real :: paa1(50)
95  real :: paa2(50)
96 !
97  real :: ptm(10)
98  real :: pdm(10,8)
99 !
100  real :: pavgm(10)
101 !
102  character*4:: isdate(3)
103  character*4:: istime(2)
104  character*4:: name(2)
105 !
106  integer :: imr
107 !
108  real :: pr65(2,65)
109  real :: pr151(2,151)
110 
111  data imr/0/
112  data isdate/'01-f','eb-0','2 '/,istime/'15:4','9:27'/
113  data name/'msis','e-00'/
114 ! temperature
115  data pt1/ &
116  9.86573e-01, 1.62228e-02, 1.55270e-02,-1.04323e-01,-3.75801e-03,&
117  -1.18538e-03,-1.24043e-01, 4.56820e-03, 8.76018e-03,-1.36235e-01,&
118  -3.52427e-02, 8.84181e-03,-5.92127e-03,-8.61650e+00, 0.00000e+00,&
119  1.28492e-02, 0.00000e+00, 1.30096e+02, 1.04567e-02, 1.65686e-03,&
120  -5.53887e-06, 2.97810e-03, 0.00000e+00, 5.13122e-03, 8.66784e-02,&
121  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,-7.27026e-06,&
122  0.00000e+00, 6.74494e+00, 4.93933e-03, 2.21656e-03, 2.50802e-03,&
123  0.00000e+00, 0.00000e+00,-2.08841e-02,-1.79873e+00, 1.45103e-03,&
124  2.81769e-04,-1.44703e-03,-5.16394e-05, 8.47001e-02, 1.70147e-01,&
125  5.72562e-03, 5.07493e-05, 4.36148e-03, 1.17863e-04, 4.74364e-03/
126  data pt2/ &
127  6.61278e-03, 4.34292e-05, 1.44373e-03, 2.41470e-05, 2.84426e-03,&
128  8.56560e-04, 2.04028e-03, 0.00000e+00,-3.15994e+03,-2.46423e-03,&
129  1.13843e-03, 4.20512e-04, 0.00000e+00,-9.77214e+01, 6.77794e-03,&
130  5.27499e-03, 1.14936e-03, 0.00000e+00,-6.61311e-03,-1.84255e-02,&
131  -1.96259e-02, 2.98618e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
132  6.44574e+02, 8.84668e-04, 5.05066e-04, 0.00000e+00, 4.02881e+03,&
133  -1.89503e-03, 0.00000e+00, 0.00000e+00, 8.21407e-04, 2.06780e-03,&
134  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
135  -1.20410e-02,-3.63963e-03, 9.92070e-05,-1.15284e-04,-6.33059e-05,&
136  -6.05545e-01, 8.34218e-03,-9.13036e+01, 3.71042e-04, 0.00000e+00/
137  data pt3/ &
138  4.19000e-04, 2.70928e-03, 3.31507e-03,-4.44508e-03,-4.96334e-03,&
139  -1.60449e-03, 3.95119e-03, 2.48924e-03, 5.09815e-04, 4.05302e-03,&
140  2.24076e-03, 0.00000e+00, 6.84256e-03, 4.66354e-04, 0.00000e+00,&
141  -3.68328e-04, 0.00000e+00, 0.00000e+00,-1.46870e+02, 0.00000e+00,&
142  0.00000e+00, 1.09501e-03, 4.65156e-04, 5.62583e-04, 3.21596e+00,&
143  6.43168e-04, 3.14860e-03, 3.40738e-03, 1.78481e-03, 9.62532e-04,&
144  5.58171e-04, 3.43731e+00,-2.33195e-01, 5.10289e-04, 0.00000e+00,&
145  0.00000e+00,-9.25347e+04, 0.00000e+00,-1.99639e-03, 0.00000e+00,&
146  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
147  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
148 ! he density
149  data pa1/ &
150  1.09979e+00,-4.88060e-02,-1.97501e-01,-9.10280e-02,-6.96558e-03,&
151  2.42136e-02, 3.91333e-01,-7.20068e-03,-3.22718e-02, 1.41508e+00,&
152  1.68194e-01, 1.85282e-02, 1.09384e-01,-7.24282e+00, 0.00000e+00,&
153  2.96377e-01,-4.97210e-02, 1.04114e+02,-8.61108e-02,-7.29177e-04,&
154  1.48998e-06, 1.08629e-03, 0.00000e+00, 0.00000e+00, 8.31090e-02,&
155  1.12818e-01,-5.75005e-02,-1.29919e-02,-1.78849e-02,-2.86343e-06,&
156  0.00000e+00,-1.51187e+02,-6.65902e-03, 0.00000e+00,-2.02069e-03,&
157  0.00000e+00, 0.00000e+00, 4.32264e-02,-2.80444e+01,-3.26789e-03,&
158  2.47461e-03, 0.00000e+00, 0.00000e+00, 9.82100e-02, 1.22714e-01,&
159  -3.96450e-02, 0.00000e+00,-2.76489e-03, 0.00000e+00, 1.87723e-03/
160  data pa2/ &
161  -8.09813e-03, 4.34428e-05,-7.70932e-03, 0.00000e+00,-2.28894e-03,&
162  -5.69070e-03,-5.22193e-03, 6.00692e-03,-7.80434e+03,-3.48336e-03,&
163  -6.38362e-03,-1.82190e-03, 0.00000e+00,-7.58976e+01,-2.17875e-02,&
164  -1.72524e-02,-9.06287e-03, 0.00000e+00, 2.44725e-02, 8.66040e-02,&
165  1.05712e-01, 3.02543e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
166  -6.01364e+03,-5.64668e-03,-2.54157e-03, 0.00000e+00, 3.15611e+02,&
167  -5.69158e-03, 0.00000e+00, 0.00000e+00,-4.47216e-03,-4.49523e-03,&
168  4.64428e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
169  4.51236e-02, 2.46520e-02, 6.17794e-03, 0.00000e+00, 0.00000e+00,&
170  -3.62944e-01,-4.80022e-02,-7.57230e+01,-1.99656e-03, 0.00000e+00/
171  data pa3/ &
172  -5.18780e-03,-1.73990e-02,-9.03485e-03, 7.48465e-03, 1.53267e-02,&
173  1.06296e-02, 1.18655e-02, 2.55569e-03, 1.69020e-03, 3.51936e-02,&
174  -1.81242e-02, 0.00000e+00,-1.00529e-01,-5.10574e-03, 0.00000e+00,&
175  2.10228e-03, 0.00000e+00, 0.00000e+00,-1.73255e+02, 5.07833e-01,&
176  -2.41408e-01, 8.75414e-03, 2.77527e-03,-8.90353e-05,-5.25148e+00,&
177  -5.83899e-03,-2.09122e-02,-9.63530e-03, 9.77164e-03, 4.07051e-03,&
178  2.53555e-04,-5.52875e+00,-3.55993e-01,-2.49231e-03, 0.00000e+00,&
179  0.00000e+00, 2.86026e+01, 0.00000e+00, 3.42722e-04, 0.00000e+00,&
180  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
181  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
182 ! o density
183  data pb1/ &
184  1.02315e+00,-1.59710e-01,-1.06630e-01,-1.77074e-02,-4.42726e-03,&
185  3.44803e-02, 4.45613e-02,-3.33751e-02,-5.73598e-02, 3.50360e-01,&
186  6.33053e-02, 2.16221e-02, 5.42577e-02,-5.74193e+00, 0.00000e+00,&
187  1.90891e-01,-1.39194e-02, 1.01102e+02, 8.16363e-02, 1.33717e-04,&
188  6.54403e-06, 3.10295e-03, 0.00000e+00, 0.00000e+00, 5.38205e-02,&
189  1.23910e-01,-1.39831e-02, 0.00000e+00, 0.00000e+00,-3.95915e-06,&
190  0.00000e+00,-7.14651e-01,-5.01027e-03, 0.00000e+00,-3.24756e-03,&
191  0.00000e+00, 0.00000e+00, 4.42173e-02,-1.31598e+01,-3.15626e-03,&
192  1.24574e-03,-1.47626e-03,-1.55461e-03, 6.40682e-02, 1.34898e-01,&
193  -2.42415e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 6.13666e-04/
194  data pb2/ &
195  -5.40373e-03, 2.61635e-05,-3.33012e-03, 0.00000e+00,-3.08101e-03,&
196  -2.42679e-03,-3.36086e-03, 0.00000e+00,-1.18979e+03,-5.04738e-02,&
197  -2.61547e-03,-1.03132e-03, 1.91583e-04,-8.38132e+01,-1.40517e-02,&
198  -1.14167e-02,-4.08012e-03, 1.73522e-04,-1.39644e-02,-6.64128e-02,&
199  -6.85152e-02,-1.34414e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
200  6.07916e+02,-4.12220e-03,-2.20996e-03, 0.00000e+00, 1.70277e+03,&
201  -4.63015e-03, 0.00000e+00, 0.00000e+00,-2.25360e-03,-2.96204e-03,&
202  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
203  3.92786e-02, 1.31186e-02,-1.78086e-03, 0.00000e+00, 0.00000e+00,&
204  -3.90083e-01,-2.84741e-02,-7.78400e+01,-1.02601e-03, 0.00000e+00/
205  data pb3/ &
206  -7.26485e-04,-5.42181e-03,-5.59305e-03, 1.22825e-02, 1.23868e-02,&
207  6.68835e-03,-1.03303e-02,-9.51903e-03, 2.70021e-04,-2.57084e-02,&
208  -1.32430e-02, 0.00000e+00,-3.81000e-02,-3.16810e-03, 0.00000e+00,&
209  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
210  0.00000e+00,-9.05762e-04,-2.14590e-03,-1.17824e-03, 3.66732e+00,&
211  -3.79729e-04,-6.13966e-03,-5.09082e-03,-1.96332e-03,-3.08280e-03,&
212  -9.75222e-04, 4.03315e+00,-2.52710e-01, 0.00000e+00, 0.00000e+00,&
213  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
214  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
215  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
216 ! n2 density
217  data pc1/ &
218  1.16112e+00, 0.00000e+00, 0.00000e+00, 3.33725e-02, 0.00000e+00,&
219  3.48637e-02,-5.44368e-03, 0.00000e+00,-6.73940e-02, 1.74754e-01,&
220  0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,&
221  1.26733e-01, 0.00000e+00, 1.03154e+02, 5.52075e-02, 0.00000e+00,&
222  0.00000e+00, 8.13525e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
223  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
224  0.00000e+00,-2.50482e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
225  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.48894e-03,&
226  6.16053e-04,-5.79716e-04, 2.95482e-03, 8.47001e-02, 1.70147e-01,&
227  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
228  data pc2/ &
229  0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
230  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
231  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
232  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
233  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
234  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
235  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
236  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
237  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
238  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
239  data pc3/ &
240  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
241  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
242  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
243  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
244  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
245  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
246  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
247  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
248  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
249  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
250 ! tlb
251  data pd1/ &
252  9.44846e-01, 0.00000e+00, 0.00000e+00,-3.08617e-02, 0.00000e+00,&
253  -2.44019e-02, 6.48607e-03, 0.00000e+00, 3.08181e-02, 4.59392e-02,&
254  0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,&
255  2.13260e-02, 0.00000e+00,-3.56958e+02, 0.00000e+00, 1.82278e-04,&
256  0.00000e+00, 3.07472e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
257  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
258  0.00000e+00, 0.00000e+00, 3.83054e-03, 0.00000e+00, 0.00000e+00,&
259  -1.93065e-03,-1.45090e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
260  0.00000e+00,-1.23493e-03, 1.36736e-03, 8.47001e-02, 1.70147e-01,&
261  3.71469e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
262  data pd2/ &
263  5.10250e-03, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
264  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
265  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
266  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
267  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
268  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
269  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
270  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
271  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
272  0.00000e+00, 3.68756e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00/
273  data pd3/ &
274  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
275  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
276  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
277  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
278  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
279  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
280  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
281  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
282  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
283  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
284 ! o2 density
285  data pe1/ &
286  1.35580e+00, 1.44816e-01, 0.00000e+00, 6.07767e-02, 0.00000e+00,&
287  2.94777e-02, 7.46900e-02, 0.00000e+00,-9.23822e-02, 8.57342e-02,&
288  0.00000e+00, 0.00000e+00, 0.00000e+00, 2.38636e+01, 0.00000e+00,&
289  7.71653e-02, 0.00000e+00, 8.18751e+01, 1.87736e-02, 0.00000e+00,&
290  0.00000e+00, 1.49667e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
291  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
292  0.00000e+00,-3.67874e+02, 5.48158e-03, 0.00000e+00, 0.00000e+00,&
293  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
294  0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,&
295  1.22631e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
296  data pe2/ &
297  8.17187e-03, 3.71617e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
298  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
299  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.10826e-03,&
300  -3.13640e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
301  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
302  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
303  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
304  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
305  -7.35742e-02,-5.00266e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
306  0.00000e+00, 1.94965e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00/
307  data pe3/ &
308  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
309  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
310  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
311  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
312  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
313  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
314  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
315  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
316  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
317  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
318 ! ar density
319  data pf1/ &
320  1.04761e+00, 2.00165e-01, 2.37697e-01, 3.68552e-02, 0.00000e+00,&
321  3.57202e-02,-2.14075e-01, 0.00000e+00,-1.08018e-01,-3.73981e-01,&
322  0.00000e+00, 3.10022e-02,-1.16305e-03,-2.07596e+01, 0.00000e+00,&
323  8.64502e-02, 0.00000e+00, 9.74908e+01, 5.16707e-02, 0.00000e+00,&
324  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
325  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
326  0.00000e+00, 3.46193e+02, 1.34297e-02, 0.00000e+00, 0.00000e+00,&
327  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-3.48509e-03,&
328  -1.54689e-04, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,&
329  1.47753e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
330  data pf2/ &
331  1.89320e-02, 3.68181e-05, 1.32570e-02, 0.00000e+00, 0.00000e+00,&
332  3.59719e-03, 7.44328e-03,-1.00023e-03,-6.50528e+03, 0.00000e+00,&
333  1.03485e-02,-1.00983e-03,-4.06916e-03,-6.60864e+01,-1.71533e-02,&
334  1.10605e-02, 1.20300e-02,-5.20034e-03, 0.00000e+00, 0.00000e+00,&
335  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
336  -2.62769e+03, 7.13755e-03, 4.17999e-03, 0.00000e+00, 1.25910e+04,&
337  0.00000e+00, 0.00000e+00, 0.00000e+00,-2.23595e-03, 4.60217e-03,&
338  5.71794e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
339  -3.18353e-02,-2.35526e-02,-1.36189e-02, 0.00000e+00, 0.00000e+00,&
340  0.00000e+00, 2.03522e-02,-6.67837e+01,-1.09724e-03, 0.00000e+00/
341  data pf3/ &
342  -1.38821e-02, 1.60468e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
343  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.51574e-02,&
344  -5.44470e-04, 0.00000e+00, 7.28224e-02, 6.59413e-02, 0.00000e+00,&
345  -5.15692e-03, 0.00000e+00, 0.00000e+00,-3.70367e+03, 0.00000e+00,&
346  0.00000e+00, 1.36131e-02, 5.38153e-03, 0.00000e+00, 4.76285e+00,&
347  -1.75677e-02, 2.26301e-02, 0.00000e+00, 1.76631e-02, 4.77162e-03,&
348  0.00000e+00, 5.39354e+00, 0.00000e+00,-7.51710e-03, 0.00000e+00,&
349  0.00000e+00,-8.82736e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
350  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
351  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
352 ! h density
353  data pg1/ &
354  1.26376e+00,-2.14304e-01,-1.49984e-01, 2.30404e-01, 2.98237e-02,&
355  2.68673e-02, 2.96228e-01, 2.21900e-02,-2.07655e-02, 4.52506e-01,&
356  1.20105e-01, 3.24420e-02, 4.24816e-02,-9.14313e+00, 0.00000e+00,&
357  2.47178e-02,-2.88229e-02, 8.12805e+01, 5.10380e-02,-5.80611e-03,&
358  2.51236e-05,-1.24083e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
359  1.58727e-01,-3.48190e-02, 0.00000e+00, 0.00000e+00, 2.89885e-05,&
360  0.00000e+00, 1.53595e+02,-1.68604e-02, 0.00000e+00, 1.01015e-02,&
361  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.84552e-04,&
362  -1.22181e-03, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,&
363  -1.04927e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,-5.91313e-03/
364  data pg2/ &
365  -2.30501e-02, 3.14758e-05, 0.00000e+00, 0.00000e+00, 1.26956e-02,&
366  8.35489e-03, 3.10513e-04, 0.00000e+00, 3.42119e+03,-2.45017e-03,&
367  -4.27154e-04, 5.45152e-04, 1.89896e-03, 2.89121e+01,-6.49973e-03,&
368  -1.93855e-02,-1.48492e-02, 0.00000e+00,-5.10576e-02, 7.87306e-02,&
369  9.51981e-02,-1.49422e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
370  2.65503e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
371  0.00000e+00, 0.00000e+00, 0.00000e+00, 6.37110e-03, 3.24789e-04,&
372  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
373  6.14274e-02, 1.00376e-02,-8.41083e-04, 0.00000e+00, 0.00000e+00,&
374  0.00000e+00,-1.27099e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00/
375  data pg3/ &
376  -3.94077e-03,-1.28601e-02,-7.97616e-03, 0.00000e+00, 0.00000e+00,&
377  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
378  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
379  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
380  0.00000e+00,-6.71465e-03,-1.69799e-03, 1.93772e-03, 3.81140e+00,&
381  -7.79290e-03,-1.82589e-02,-1.25860e-02,-1.04311e-02,-3.02465e-03,&
382  2.43063e-03, 3.63237e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
383  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
384  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
385  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
386 ! n density
387  data ph1/ &
388  7.09557e+01,-3.26740e-01, 0.00000e+00,-5.16829e-01,-1.71664e-03,&
389  9.09310e-02,-6.71500e-01,-1.47771e-01,-9.27471e-02,-2.30862e-01,&
390  -1.56410e-01, 1.34455e-02,-1.19717e-01, 2.52151e+00, 0.00000e+00,&
391  -2.41582e-01, 5.92939e-02, 4.39756e+00, 9.15280e-02, 4.41292e-03,&
392  0.00000e+00, 8.66807e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
393  1.58727e-01, 9.74701e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
394  0.00000e+00, 6.70217e+01,-1.31660e-03, 0.00000e+00,-1.65317e-02,&
395  0.00000e+00, 0.00000e+00, 8.50247e-02, 2.77428e+01, 4.98658e-03,&
396  6.15115e-03, 9.50156e-03,-2.12723e-02, 8.47001e-02, 1.70147e-01,&
397  -2.38645e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.37380e-03/
398  data ph2/ &
399  -8.41918e-03, 2.80145e-05, 7.12383e-03, 0.00000e+00,-1.66209e-02,&
400  1.03533e-04,-1.68898e-02, 0.00000e+00, 3.64526e+03, 0.00000e+00,&
401  6.54077e-03, 3.69130e-04, 9.94419e-04, 8.42803e+01,-1.16124e-02,&
402  -7.74414e-03,-1.68844e-03, 1.42809e-03,-1.92955e-03, 1.17225e-01,&
403  -2.41512e-02, 1.50521e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
404  1.60261e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
405  0.00000e+00, 0.00000e+00, 0.00000e+00,-3.54403e-04,-1.87270e-02,&
406  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
407  2.76439e-02, 6.43207e-03,-3.54300e-02, 0.00000e+00, 0.00000e+00,&
408  0.00000e+00,-2.80221e-02, 8.11228e+01,-6.75255e-04, 0.00000e+00/
409  data ph3/ &
410  -1.05162e-02,-3.48292e-03,-6.97321e-03, 0.00000e+00, 0.00000e+00,&
411  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
412  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
413  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
414  0.00000e+00,-1.45546e-03,-1.31970e-02,-3.57751e-03,-1.09021e+00,&
415  -1.50181e-02,-7.12841e-03,-6.64590e-03,-3.52610e-03,-1.87773e-02,&
416  -2.22432e-03,-3.93895e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
417  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
418  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
419  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
420 ! hot o density
421  data pi1/ &
422  6.04050e-02, 1.57034e+00, 2.99387e-02, 0.00000e+00, 0.00000e+00,&
423  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-1.51018e+00,&
424  0.00000e+00, 0.00000e+00, 0.00000e+00,-8.61650e+00, 1.26454e-02,&
425  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
426  0.00000e+00, 5.50878e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,&
427  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
428  0.00000e+00, 0.00000e+00, 6.23881e-02, 0.00000e+00, 0.00000e+00,&
429  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
430  0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,&
431  -9.45934e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
432  data pi2/ &
433  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
434  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
435  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
436  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
437  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
438  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
439  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
440  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
441  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
442  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
443  data pi3/ &
444  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
445  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
446  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
447  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
448  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
449  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
450  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
451  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
452  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
453  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
454 ! s param
455  data pj1/ &
456  9.56827e-01, 6.20637e-02, 3.18433e-02, 0.00000e+00, 0.00000e+00,&
457  3.94900e-02, 0.00000e+00, 0.00000e+00,-9.24882e-03,-7.94023e-03,&
458  0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,&
459  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
460  0.00000e+00, 2.74677e-03, 0.00000e+00, 1.54951e-02, 8.66784e-02,&
461  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
462  0.00000e+00, 0.00000e+00, 0.00000e+00,-6.99007e-04, 0.00000e+00,&
463  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
464  0.00000e+00, 1.24362e-02,-5.28756e-03, 8.47001e-02, 1.70147e-01,&
465  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
466  data pj2/ &
467  0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
468  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
469  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
470  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
471  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
472  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
473  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
474  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
475  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
476  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
477  data pj3/ &
478  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
479  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
480  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
481  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
482  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
483  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
484  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
485  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
486  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
487  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
488 ! turbo
489  data pk1/ &
490  1.09930e+00, 3.90631e+00, 3.07165e+00, 9.86161e-01, 1.63536e+01,&
491  4.63830e+00, 1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
492  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
493  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
494  0.00000e+00, 0.00000e+00, 1.28840e+00, 3.10302e-02, 1.18339e-01,&
495  1.00000e+00, 7.00000e-01, 1.15020e+00, 3.44689e+00, 1.28840e+00,&
496  1.00000e+00, 1.08738e+00, 1.22947e+00, 1.10016e+00, 7.34129e-01,&
497  1.15241e+00, 2.22784e+00, 7.95046e-01, 4.01612e+00, 4.47749e+00,&
498  1.23435e+02,-7.60535e-02, 1.68986e-06, 7.44294e-01, 1.03604e+00,&
499  1.72783e+02, 1.15020e+00, 3.44689e+00,-7.46230e-01, 9.49154e-01/
500 ! lower boundary
501  data ptm/ &
502  1.04130e+03, 3.86000e+02, 1.95000e+02, 1.66728e+01, 2.13000e+02,&
503  1.20000e+02, 2.40000e+02, 1.87000e+02,-2.00000e+00, 0.00000e+00/
504  data pdm/ &
505  2.45600e+07, 6.71072e-06, 1.00000e+02, 0.00000e+00, 1.10000e+02,&
506  1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
507  8.59400e+10, 1.00000e+00, 1.05000e+02,-8.00000e+00, 1.10000e+02,&
508  1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,&
509  2.81000e+11, 0.00000e+00, 1.05000e+02, 2.80000e+01, 2.89500e+01,&
510  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
511  3.30000e+10, 2.68270e-01, 1.05000e+02, 1.00000e+00, 1.10000e+02,&
512  1.00000e+01, 1.10000e+02,-1.00000e+01, 0.00000e+00, 0.00000e+00,&
513  1.33000e+09, 1.19615e-02, 1.05000e+02, 0.00000e+00, 1.10000e+02,&
514  1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
515  1.76100e+05, 1.00000e+00, 9.50000e+01,-8.00000e+00, 1.10000e+02,&
516  1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,&
517  1.00000e+07, 1.00000e+00, 1.05000e+02,-8.00000e+00, 1.10000e+02,&
518  1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,&
519  1.00000e+06, 1.00000e+00, 1.05000e+02,-8.00000e+00, 5.50000e+02,&
520  7.60000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 4.00000e+03/
521 ! tn1(2)
522  data pl1/ &
523  1.00858e+00, 4.56011e-02,-2.22972e-02,-5.44388e-02, 5.23136e-04,&
524  -1.88849e-02, 5.23707e-02,-9.43646e-03, 6.31707e-03,-7.80460e-02,&
525  -4.88430e-02, 0.00000e+00, 0.00000e+00,-7.60250e+00, 0.00000e+00,&
526  -1.44635e-02,-1.76843e-02,-1.21517e+02, 2.85647e-02, 0.00000e+00,&
527  0.00000e+00, 6.31792e-04, 0.00000e+00, 5.77197e-03, 8.66784e-02,&
528  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
529  0.00000e+00,-8.90272e+03, 3.30611e-03, 3.02172e-03, 0.00000e+00,&
530  -2.13673e-03,-3.20910e-04, 0.00000e+00, 0.00000e+00, 2.76034e-03,&
531  2.82487e-03,-2.97592e-04,-4.21534e-03, 8.47001e-02, 1.70147e-01,&
532  8.96456e-03, 0.00000e+00,-1.08596e-02, 0.00000e+00, 0.00000e+00/
533  data pl2/ &
534  5.57917e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
535  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
536  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
537  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
538  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
539  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
540  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
541  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
542  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
543  0.00000e+00, 9.65405e-03, 0.00000e+00, 0.00000e+00, 2.00000e+00/
544 ! tn1(3)
545  data pm1/ &
546  9.39664e-01, 8.56514e-02,-6.79989e-03, 2.65929e-02,-4.74283e-03,&
547  1.21855e-02,-2.14905e-02, 6.49651e-03,-2.05477e-02,-4.24952e-02,&
548  0.00000e+00, 0.00000e+00, 0.00000e+00, 1.19148e+01, 0.00000e+00,&
549  1.18777e-02,-7.28230e-02,-8.15965e+01, 1.73887e-02, 0.00000e+00,&
550  0.00000e+00, 0.00000e+00,-1.44691e-02, 2.80259e-04, 8.66784e-02,&
551  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
552  0.00000e+00, 2.16584e+02, 3.18713e-03, 7.37479e-03, 0.00000e+00,&
553  -2.55018e-03,-3.92806e-03, 0.00000e+00, 0.00000e+00,-2.89757e-03,&
554  -1.33549e-03, 1.02661e-03, 3.53775e-04, 8.47001e-02, 1.70147e-01,&
555  -9.17497e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
556  data pm2/ &
557  3.56082e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
558  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
559  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
560  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
561  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
562  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
563  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
564  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
565  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
566  0.00000e+00,-1.00902e-02, 0.00000e+00, 0.00000e+00, 2.00000e+00/
567 ! tn1(4)
568  data pn1/ &
569  9.85982e-01,-4.55435e-02, 1.21106e-02, 2.04127e-02,-2.40836e-03,&
570  1.11383e-02,-4.51926e-02, 1.35074e-02,-6.54139e-03, 1.15275e-01,&
571  1.28247e-01, 0.00000e+00, 0.00000e+00,-5.30705e+00, 0.00000e+00,&
572  -3.79332e-02,-6.24741e-02, 7.71062e-01, 2.96315e-02, 0.00000e+00,&
573  0.00000e+00, 0.00000e+00, 6.81051e-03,-4.34767e-03, 8.66784e-02,&
574  1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
575  0.00000e+00, 1.07003e+01,-2.76907e-03, 4.32474e-04, 0.00000e+00,&
576  1.31497e-03,-6.47517e-04, 0.00000e+00,-2.20621e+01,-1.10804e-03,&
577  -8.09338e-04, 4.18184e-04, 4.29650e-03, 8.47001e-02, 1.70147e-01,&
578  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
579  data pn2/ &
580  -4.04337e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
581  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
582  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-9.52550e-04,&
583  8.56253e-04, 4.33114e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
584  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.21223e-03,&
585  2.38694e-04, 9.15245e-04, 1.28385e-03, 8.67668e-04,-5.61425e-06,&
586  1.04445e+00, 3.41112e+01, 0.00000e+00,-8.40704e-01,-2.39639e+02,&
587  7.06668e-01,-2.05873e+01,-3.63696e-01, 2.39245e+01, 0.00000e+00,&
588  -1.06657e-03,-7.67292e-04, 1.54534e-04, 0.00000e+00, 0.00000e+00,&
589  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
590 ! tn1(5) tn2(1)
591  data po1/ &
592  1.00320e+00, 3.83501e-02,-2.38983e-03, 2.83950e-03, 4.20956e-03,&
593  5.86619e-04, 2.19054e-02,-1.00946e-02,-3.50259e-03, 4.17392e-02,&
594  -8.44404e-03, 0.00000e+00, 0.00000e+00, 4.96949e+00, 0.00000e+00,&
595  -7.06478e-03,-1.46494e-02, 3.13258e+01,-1.86493e-03, 0.00000e+00,&
596  -1.67499e-02, 0.00000e+00, 0.00000e+00, 5.12686e-04, 8.66784e-02,&
597  1.58727e-01,-4.64167e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
598  4.37353e-03,-1.99069e+02, 0.00000e+00,-5.34884e-03, 0.00000e+00,&
599  1.62458e-03, 2.93016e-03, 2.67926e-03, 5.90449e+02, 0.00000e+00,&
600  0.00000e+00,-1.17266e-03,-3.58890e-04, 8.47001e-02, 1.70147e-01,&
601  0.00000e+00, 0.00000e+00, 1.38673e-02, 0.00000e+00, 0.00000e+00/
602  data po2/ &
603  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
604  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
605  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.60571e-03,&
606  6.28078e-04, 5.05469e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
607  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-1.57829e-03,&
608  -4.00855e-04, 5.04077e-05,-1.39001e-03,-2.33406e-03,-4.81197e-04,&
609  1.46758e+00, 6.20332e+00, 0.00000e+00, 3.66476e-01,-6.19760e+01,&
610  3.09198e-01,-1.98999e+01, 0.00000e+00,-3.29933e+02, 0.00000e+00,&
611  -1.10080e-03,-9.39310e-05, 1.39638e-04, 0.00000e+00, 0.00000e+00,&
612  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
613 ! tn2(2)
614  data pp1/ &
615  9.81637e-01,-1.41317e-03, 3.87323e-02, 0.00000e+00, 0.00000e+00,&
616  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-3.58707e-02,&
617  -8.63658e-03, 0.00000e+00, 0.00000e+00,-2.02226e+00, 0.00000e+00,&
618  -8.69424e-03,-1.91397e-02, 8.76779e+01, 4.52188e-03, 0.00000e+00,&
619  2.23760e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
620  0.00000e+00,-7.07572e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
621  -4.11210e-03, 3.50060e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
622  0.00000e+00, 0.00000e+00,-8.36657e-03, 1.61347e+01, 0.00000e+00,&
623  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
624  0.00000e+00, 0.00000e+00,-1.45130e-02, 0.00000e+00, 0.00000e+00/
625  data pp2/ &
626  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
627  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
628  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.24152e-03,&
629  6.43365e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
630  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.33255e-03,&
631  2.42657e-03, 1.60666e-03,-1.85728e-03,-1.46874e-03,-4.79163e-06,&
632  1.22464e+00, 3.53510e+01, 0.00000e+00, 4.49223e-01,-4.77466e+01,&
633  4.70681e-01, 8.41861e+00,-2.88198e-01, 1.67854e+02, 0.00000e+00,&
634  7.11493e-04, 6.05601e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
635  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
636 ! tn2(3)
637  data pq1/ &
638  1.00422e+00,-7.11212e-03, 5.24480e-03, 0.00000e+00, 0.00000e+00,&
639  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-5.28914e-02,&
640  -2.41301e-02, 0.00000e+00, 0.00000e+00,-2.12219e+01,-1.03830e-02,&
641  -3.28077e-03, 1.65727e-02, 1.68564e+00,-6.68154e-03, 0.00000e+00,&
642  1.45155e-02, 0.00000e+00, 8.42365e-03, 0.00000e+00, 0.00000e+00,&
643  0.00000e+00,-4.34645e-03, 0.00000e+00, 0.00000e+00, 2.16780e-02,&
644  0.00000e+00,-1.38459e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
645  0.00000e+00, 0.00000e+00, 7.04573e-03,-4.73204e+01, 0.00000e+00,&
646  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
647  0.00000e+00, 0.00000e+00, 1.08767e-02, 0.00000e+00, 0.00000e+00/
648  data pq2/ &
649  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
650  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-8.08279e-03,&
651  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.21769e-04,&
652  -2.27387e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
653  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.26769e-03,&
654  3.16901e-03, 4.60316e-04,-1.01431e-04, 1.02131e-03, 9.96601e-04,&
655  1.25707e+00, 2.50114e+01, 0.00000e+00, 4.24472e-01,-2.77655e+01,&
656  3.44625e-01, 2.75412e+01, 0.00000e+00, 7.94251e+02, 0.00000e+00,&
657  2.45835e-03, 1.38871e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
658  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
659 ! tn2(4) tn3(1)
660  data pr1/ &
661  1.01890e+00,-2.46603e-02, 1.00078e-02, 0.00000e+00, 0.00000e+00,&
662  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-6.70977e-02,&
663  -4.02286e-02, 0.00000e+00, 0.00000e+00,-2.29466e+01,-7.47019e-03,&
664  2.26580e-03, 2.63931e-02, 3.72625e+01,-6.39041e-03, 0.00000e+00,&
665  9.58383e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
666  0.00000e+00,-1.85291e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
667  0.00000e+00, 1.39717e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
668  0.00000e+00, 0.00000e+00, 9.19771e-03,-3.69121e+02, 0.00000e+00,&
669  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
670  0.00000e+00, 0.00000e+00,-1.57067e-02, 0.00000e+00, 0.00000e+00/
671  data pr2/ &
672  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
673  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-7.07265e-03,&
674  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.92953e-03,&
675  -2.77739e-03,-4.40092e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
676  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.47280e-03,&
677  2.95035e-04,-1.81246e-03, 2.81945e-03, 4.27296e-03, 9.78863e-04,&
678  1.40545e+00,-6.19173e+00, 0.00000e+00, 0.00000e+00,-7.93632e+01,&
679  4.44643e-01,-4.03085e+02, 0.00000e+00, 1.15603e+01, 0.00000e+00,&
680  2.25068e-03, 8.48557e-04,-2.98493e-04, 0.00000e+00, 0.00000e+00,&
681  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
682 ! tn3(2)
683  data ps1/ &
684  9.75801e-01, 3.80680e-02,-3.05198e-02, 0.00000e+00, 0.00000e+00,&
685  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.85575e-02,&
686  5.04057e-02, 0.00000e+00, 0.00000e+00,-1.76046e+02, 1.44594e-02,&
687  -1.48297e-03,-3.68560e-03, 3.02185e+01,-3.23338e-03, 0.00000e+00,&
688  1.53569e-02, 0.00000e+00,-1.15558e-02, 0.00000e+00, 0.00000e+00,&
689  0.00000e+00, 4.89620e-03, 0.00000e+00, 0.00000e+00,-1.00616e-02,&
690  -8.21324e-03,-1.57757e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
691  0.00000e+00, 0.00000e+00, 6.63564e-03, 4.58410e+01, 0.00000e+00,&
692  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
693  0.00000e+00, 0.00000e+00,-2.51280e-02, 0.00000e+00, 0.00000e+00/
694  data ps2/ &
695  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
696  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 9.91215e-03,&
697  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-8.73148e-04,&
698  -1.29648e-03,-7.32026e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
699  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-4.68110e-03,&
700  -4.66003e-03,-1.31567e-03,-7.39390e-04, 6.32499e-04,-4.65588e-04,&
701  -1.29785e+00,-1.57139e+02, 0.00000e+00, 2.58350e-01,-3.69453e+01,&
702  4.10672e-01, 9.78196e+00,-1.52064e-01,-3.85084e+03, 0.00000e+00,&
703  -8.52706e-04,-1.40945e-03,-7.26786e-04, 0.00000e+00, 0.00000e+00,&
704  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
705 ! tn3(3)
706  data pu1/ &
707  9.60722e-01, 7.03757e-02,-3.00266e-02, 0.00000e+00, 0.00000e+00,&
708  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.22671e-02,&
709  4.10423e-02, 0.00000e+00, 0.00000e+00,-1.63070e+02, 1.06073e-02,&
710  5.40747e-04, 7.79481e-03, 1.44908e+02, 1.51484e-04, 0.00000e+00,&
711  1.97547e-02, 0.00000e+00,-1.41844e-02, 0.00000e+00, 0.00000e+00,&
712  0.00000e+00, 5.77884e-03, 0.00000e+00, 0.00000e+00, 9.74319e-03,&
713  0.00000e+00,-2.88015e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
714  0.00000e+00, 0.00000e+00,-4.44902e-03,-2.92760e+01, 0.00000e+00,&
715  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
716  0.00000e+00, 0.00000e+00, 2.34419e-02, 0.00000e+00, 0.00000e+00/
717  data pu2/ &
718  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
719  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.36685e-03,&
720  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-4.65325e-04,&
721  -5.50628e-04, 3.31465e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
722  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.06179e-03,&
723  -3.08575e-03,-7.93589e-04,-1.08629e-04, 5.95511e-04,-9.05050e-04,&
724  1.18997e+00, 4.15924e+01, 0.00000e+00,-4.72064e-01,-9.47150e+02,&
725  3.98723e-01, 1.98304e+01, 0.00000e+00, 3.73219e+03, 0.00000e+00,&
726  -1.50040e-03,-1.14933e-03,-1.56769e-04, 0.00000e+00, 0.00000e+00,&
727  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
728 ! tn3(4)
729  data pv1/ &
730  1.03123e+00,-7.05124e-02, 8.71615e-03, 0.00000e+00, 0.00000e+00,&
731  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-3.82621e-02,&
732  -9.80975e-03, 0.00000e+00, 0.00000e+00, 2.89286e+01, 9.57341e-03,&
733  0.00000e+00, 0.00000e+00, 8.66153e+01, 7.91938e-04, 0.00000e+00,&
734  0.00000e+00, 0.00000e+00, 4.68917e-03, 0.00000e+00, 0.00000e+00,&
735  0.00000e+00, 7.86638e-03, 0.00000e+00, 0.00000e+00, 9.90827e-03,&
736  0.00000e+00, 6.55573e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
737  0.00000e+00, 0.00000e+00, 0.00000e+00,-4.00200e+01, 0.00000e+00,&
738  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
739  0.00000e+00, 0.00000e+00, 7.07457e-03, 0.00000e+00, 0.00000e+00/
740  data pv2/ &
741  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
742  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.72268e-03,&
743  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.04970e-04,&
744  1.21560e-03,-8.05579e-06, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
745  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.49941e-03,&
746  -4.57256e-04,-1.59311e-04, 2.96481e-04,-1.77318e-03,-6.37918e-04,&
747  1.02395e+00, 1.28172e+01, 0.00000e+00, 1.49903e-01,-2.63818e+01,&
748  0.00000e+00, 4.70628e+01,-2.22139e-01, 4.82292e-02, 0.00000e+00,&
749  -8.67075e-04,-5.86479e-04, 5.32462e-04, 0.00000e+00, 0.00000e+00,&
750  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
751 ! tn3(5) surface temp tsl
752  data pw1/ &
753  1.00828e+00,-9.10404e-02,-2.26549e-02, 0.00000e+00, 0.00000e+00,&
754  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-2.32420e-02,&
755  -9.08925e-03, 0.00000e+00, 0.00000e+00, 3.36105e+01, 0.00000e+00,&
756  0.00000e+00, 0.00000e+00,-1.24957e+01,-5.87939e-03, 0.00000e+00,&
757  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
758  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
759  0.00000e+00, 2.79765e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
760  0.00000e+00, 0.00000e+00, 0.00000e+00, 2.01237e+03, 0.00000e+00,&
761  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
762  0.00000e+00, 0.00000e+00,-1.75553e-02, 0.00000e+00, 0.00000e+00/
763  data pw2/ &
764  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
765  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
766  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.29699e-03,&
767  1.26659e-03, 2.68402e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
768  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.17894e-03,&
769  1.48746e-03, 1.06478e-04, 1.34743e-04,-2.20939e-03,-6.23523e-04,&
770  6.36539e-01, 1.13621e+01, 0.00000e+00,-3.93777e-01, 2.38687e+03,&
771  0.00000e+00, 6.61865e+02,-1.21434e-01, 9.27608e+00, 0.00000e+00,&
772  1.68478e-04, 1.24892e-03, 1.71345e-03, 0.00000e+00, 0.00000e+00,&
773  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
774 ! tgn3(2) surface grad tslg
775  data px1/ &
776  1.57293e+00,-6.78400e-01, 6.47500e-01, 0.00000e+00, 0.00000e+00,&
777  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-7.62974e-02,&
778  -3.60423e-01, 0.00000e+00, 0.00000e+00, 1.28358e+02, 0.00000e+00,&
779  0.00000e+00, 0.00000e+00, 4.68038e+01, 0.00000e+00, 0.00000e+00,&
780  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
781  0.00000e+00,-1.67898e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
782  0.00000e+00, 2.90994e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
783  0.00000e+00, 0.00000e+00, 0.00000e+00, 3.15706e+01, 0.00000e+00,&
784  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
785  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
786  data px2/ &
787  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
788  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
789  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
790  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
791  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
792  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
793  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
794  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
795  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
796  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
797 ! tgn2(1) tgn1(2)
798  data py1/ &
799  8.60028e-01, 3.77052e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
800  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-1.17570e+00,&
801  0.00000e+00, 0.00000e+00, 0.00000e+00, 7.77757e-03, 0.00000e+00,&
802  0.00000e+00, 0.00000e+00, 1.01024e+02, 0.00000e+00, 0.00000e+00,&
803  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
804  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
805  0.00000e+00, 6.54251e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
806  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
807  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
808  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
809  data py2/ &
810  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
811  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
812  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
813  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
814  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,-1.56959e-02,&
815  1.91001e-02, 3.15971e-02, 1.00982e-02,-6.71565e-03, 2.57693e-03,&
816  1.38692e+00, 2.82132e-01, 0.00000e+00, 0.00000e+00, 3.81511e+02,&
817  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
818  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
819  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
820 ! tgn3(1) tgn2(2)
821  data pz1/ &
822  1.06029e+00,-5.25231e-02, 3.73034e-01, 0.00000e+00, 0.00000e+00,&
823  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.31072e-02,&
824  -3.88409e-01, 0.00000e+00, 0.00000e+00,-1.65295e+02,-2.13801e-01,&
825  -4.38916e-02,-3.22716e-01,-8.82393e+01, 1.18458e-01, 0.00000e+00,&
826  -4.35863e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
827  0.00000e+00,-1.19782e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
828  0.00000e+00, 2.62229e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
829  0.00000e+00, 0.00000e+00, 0.00000e+00,-5.37443e+01, 0.00000e+00,&
830  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
831  0.00000e+00, 0.00000e+00,-4.55788e-01, 0.00000e+00, 0.00000e+00/
832  data pz2/ &
833  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
834  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
835  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.84009e-02,&
836  3.96733e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
837  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.05494e-02,&
838  7.39617e-02, 1.92200e-02,-8.46151e-03,-1.34244e-02, 1.96338e-02,&
839  1.50421e+00, 1.88368e+01, 0.00000e+00, 0.00000e+00,-5.13114e+01,&
840  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
841  5.11923e-02, 3.61225e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
842  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00/
843 ! semiannual mult sam
844  data paa1/ &
845  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
846  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
847  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
848  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
849  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
850  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
851  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
852  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
853  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,&
854  1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00/
855  data paa2/ &
856  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
857  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
858  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
859  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
860  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
861  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
862  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
863  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
864  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,&
865  0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00/
866 ! middle atmosphere averages
867  data pavgm/ &
868  2.61000e+02, 2.64000e+02, 2.29000e+02, 2.17000e+02, 2.17000e+02,&
869  2.23000e+02, 2.86760e+02,-2.93940e+00, 2.50000e+00, 0.00000e+00/
870  end module wam_gtd7bk_mod
871 
872 !-----------------------------------------------------------------------
879  module gettemp_mod
880 !
881  real :: tlb
882  real :: s
883  real :: db04
884  real :: db16
885  real :: db28
886  real :: db32
887  real :: db40
888  real :: db48
889  real :: db01
890  real :: za
891  real :: t0
892  real :: z0
893  real :: g0
894  real :: rl
895  real :: dd
896  real :: db14
897  real :: tr12
898 !
899  real :: tn1(5)
900  real :: tn2(4)
901  real :: tn3(5)
902  real :: tgn1(2)
903  real :: tgn2(2)
904  real :: tgn3(2)
905 !
906  real :: pt(150)
907  real :: pd(150,9)
908  real :: ps(150)
909  real :: pdl(25,2)
910  real :: ptl(100,4)
911  real :: pma(100,10)
912  real :: sam(100)
913 !
914  real :: sw(25)
915  real :: swc(25)
916 !
917  real :: dm04
918  real :: dm16
919  real :: dm28
920  real :: dm32
921  real :: dm40
922  real :: dm01
923  real :: dm14
924 !
925  real :: gsurf
926  real :: re
927 !
928  real :: tinfg
929  real :: tt(15)
930 !
931  real :: plg(9,4)
932  real :: ctloc
933  real :: stloc
934  real :: c2tloc
935  real :: s2tloc
936  real :: c3tloc
937  real :: s3tloc
938  real :: day
939  real :: df
940  real :: dfa
941  real :: apd
942  real :: apdf
943  real :: apt(4)
944  real :: xlong
945 !
946  integer :: isw
947  integer :: iyr
948 !
949  end module gettemp_mod
950 
951 ! ----------------------------------------------------------------------
967  subroutine gettemp(iday,nday,xlat,nlat,pr,np,temp,n_o,n_o2,n_n2)
968  implicit none
969  integer, intent(in) :: nday ! number of days
970  integer, intent(in) :: nlat ! number of latitudes
971  integer, intent(in) :: np ! number of pressure layer
972  real, intent(in) :: pr(np) ! pressure in mb
973  real, intent(in) :: xlat(nlat) ! latitude in degree
974  integer, intent(in) :: iday(nday) ! calender day
975  real, intent(out) :: temp(np,nlat,nday) ! temperature
976  real, intent(out) :: n_o(np,nlat,nday) ! number density of o
977  real, intent(out) :: n_o2(np,nlat,nday) ! number density of o2
978  real, intent(out) :: n_n2(np,nlat,nday) ! number density of n2
979  real :: alt(np,nlat,nday) ! altitude in km
980  real :: d(9),t(2),sw(25),ap(7),ut,xlong,xlst,f107, &
981  f107a
982  integer :: k,il,ip
983 ! set magnetic index average value
984  data ap/7*9./
985 ! set swich 7,8,10,14 zero to avoid diurnal changes in output temperatu
986 ! swich #7 is for diurnal,#8 is for semidiurnal,# 10 is for all ut/longi
987 ! effect,#14 is for terdiurnal
988  data sw/1.,1.,1.,1.,1.,1.,0.,0.,1.,0.,1.,1.,1.,0.,1.,1.,1.,1.,1., &
989  1.,1.,1.,1.,1.,1./
990 ! set 10.7cm flux be average value
991  f107=150.
992  f107a=150.
993 ! turn on swich
994  call tselec(sw)
995 ! set longitude, ut, local time , it should not make difference to outpu
996  ut=0.
997  xlong=0.
998  xlst=ut/3600.+xlong/15.
999 ! calculate temperature for each lat,pres level,day
1000  do k=1,nday
1001  do il=1,nlat
1002  do ip=1,np
1003  call ghp7(iday(k),ut,alt(ip,il,k),xlat(il),xlong,xlst,f107a,f107, &
1004  ap,d,t,pr(ip))
1005  temp(ip,il,k)=t(2)
1006  n_o(ip,il,k)=d(2)
1007  n_o2(ip,il,k)=d(4)
1008  n_n2(ip,il,k)=d(3)
1009  enddo
1010  enddo
1011  enddo
1012  end subroutine gettemp
1013 
1014 !-----------------------------------------------------------------------
1143  subroutine gtd7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
1144  use wam_gtd7bk_mod
1145  use gettemp_mod, only: dd, tn1,tn2,tn3,tgn1,tgn2,tgn3, &
1146  pt,pd,ps,pdl,ptl,pma,sam,sw,isw,&
1147  dm28,gsurf,re
1148 
1149  dimension d(9),t(2),ap(7),ds(9),ts(2)
1150  dimension zn3(5),zn2(4),sv(25)
1151 !
1152  save
1153  data mn3/5/,zn3/32.5,20.,15.,10.,0./
1154  data mn2/4/,zn2/72.5,55.,45.,32.5/
1155  data zmix/62.5/,alast/99999./,mssl/-999/
1156  data sv/25*1./
1157 ! ==== assign common/parm7/
1158  pt(1:50) =pt1(1:50); pt(51:100) =pt2(1:50); pt(101:150) =pt3(1:50)
1159  pd(1:50,1)=pa1(1:50); pd(51:100,1)=pa2(1:50); pd(101:150,1)=pa3(1:50)
1160  pd(1:50,2)=pb1(1:50); pd(51:100,2)=pb2(1:50); pd(101:150,2)=pb3(1:50)
1161  pd(1:50,3)=pc1(1:50); pd(51:100,3)=pc2(1:50); pd(101:150,3)=pc3(1:50)
1162  pd(1:50,4)=pd1(1:50); pd(51:100,4)=pd2(1:50); pd(101:150,4)=pd3(1:50)
1163  pd(1:50,5)=pe1(1:50); pd(51:100,5)=pe2(1:50); pd(101:150,5)=pe3(1:50)
1164  pd(1:50,6)=pf1(1:50); pd(51:100,6)=pf2(1:50); pd(101:150,6)=pf3(1:50)
1165  pd(1:50,7)=pg1(1:50); pd(51:100,7)=pg2(1:50); pd(101:150,7)=pg3(1:50)
1166  pd(1:50,8)=ph1(1:50); pd(51:100,8)=ph2(1:50); pd(101:150,8)=ph3(1:50)
1167  pd(1:50,9)=pi1(1:50); pd(51:100,9)=pi2(1:50); pd(101:150,9)=pi3(1:50)
1168  ps(1:50) =pj1(1:50); ps(51:100) =pj2(1:50); ps(101:150) =pj3(1:50)
1169  pdl(1:25,1)=pk1(1:25); pdl(1:25,2)=pk1(26:50)
1170  ptl(1:50,1)=pl1(1:50); ptl(51:100,1)=pl2(1:50)
1171  ptl(1:50,2)=pm1(1:50); ptl(51:100,2)=pm2(1:50)
1172  ptl(1:50,3)=pn1(1:50); ptl(51:100,3)=pn2(1:50)
1173  ptl(1:50,4)=po1(1:50); ptl(51:100,4)=po2(1:50)
1174  pma(1:50,1)=pp1(1:50); pma(51:100,1)=pp2(1:50)
1175  pma(1:50,2)=pq1(1:50); pma(51:100,2)=pq2(1:50)
1176  pma(1:50,3)=pr1(1:50); pma(51:100,3)=pr2(1:50)
1177  pma(1:50,4)=ps1(1:50); pma(51:100,4)=ps2(1:50)
1178  pma(1:50,5)=pu1(1:50); pma(51:100,5)=pu2(1:50)
1179  pma(1:50,6)=pv1(1:50); pma(51:100,6)=pv2(1:50)
1180  pma(1:50,7)=pw1(1:50); pma(51:100,7)=pw2(1:50)
1181  pma(1:50,8)=px1(1:50); pma(51:100,8)=px2(1:50)
1182  pma(1:50,9)=py1(1:50); pma(51:100,9)=py2(1:50)
1183  pma(1:50,10)=pz1(1:50); pma(51:100,10)=pz2(1:50)
1184  sam(1:50)=paa1(1:50); sam(51:100)=paa2(1:50)
1185 !
1186  if(isw.ne.64999) call tselec(sv)
1187 !
1188 ! test for changed input
1189  v1=vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,1)
1190 ! latitude variation of gravity (none for sw(2)=0)
1191  xlat=glat
1192  if(sw(2).eq.0) xlat=45.
1193  call glatf(xlat,gsurf,re)
1194 !
1195  xmm=pdm(5,3)
1196 !
1197 ! thermosphere/mesosphere (above zn2(1))
1198  altt=amax1(alt,zn2(1))
1199  mss=mass
1200 ! only calculate n2 in thermosphere if alt in mixed region
1201  if(alt.lt.zmix.and.mass.gt.0) mss=28
1202 ! only calculate thermosphere if input parameters changed
1203 ! or altitude above zn2(1) in mesosphere
1204  if(v1.eq.1..or.alt.gt.zn2(1).or.alast.gt.zn2(1).or.mss.ne.mssl) then
1205  call gts7(iyd,sec,altt,glat,glong,stl,f107a,f107,ap,mss,ds,ts)
1206  dm28m=dm28
1207 ! metric adjustment
1208  if(imr.eq.1) dm28m=dm28*1.e6
1209  mssl=mss
1210  endif
1211  t(1)=ts(1)
1212  t(2)=ts(2)
1213  if(alt.ge.zn2(1)) then
1214  do 5 j=1,9
1215  d(j)=ds(j)
1216  5 continue
1217  goto 10
1218  endif
1219 !
1220 ! lower mesosphere/upper stratosphere [between zn3(1) and zn2(1)]
1221 ! temperature at nodes and gradients at end nodes
1222 ! inverse temperature a linear function of spherical harmonics
1223 ! only calculate nodes if input changed
1224  if(v1.eq.1..or.alast.ge.zn2(1)) then
1225  tgn2(1)=tgn1(2)
1226  tn2(1)=tn1(5)
1227  tn2(2)=pma(1,1)*pavgm(1)/(1.-sw(20)*glob7s(pma(1,1)))
1228  tn2(3)=pma(1,2)*pavgm(2)/(1.-sw(20)*glob7s(pma(1,2)))
1229  tn2(4)=pma(1,3)*pavgm(3)/(1.-sw(20)*sw(22)*glob7s(pma(1,3)))
1230  tgn2(2)=pavgm(9)*pma(1,10)*(1.+sw(20)*sw(22)*glob7s(pma(1,10))) &
1231  *tn2(4)*tn2(4)/(pma(1,3)*pavgm(3))**2
1232  tn3(1)=tn2(4)
1233  endif
1234  if(alt.ge.zn3(1)) goto 6
1235 !
1236 ! lower stratosphere and troposphere [below zn3(1)]
1237 ! temperature at nodes and gradients at end nodes
1238 ! inverse temperature a linear function of spherical harmonics
1239 ! only calculate nodes if input changed
1240  if(v1.eq.1..or.alast.ge.zn3(1)) then
1241  tgn3(1)=tgn2(2)
1242  tn3(2)=pma(1,4)*pavgm(4)/(1.-sw(22)*glob7s(pma(1,4)))
1243  tn3(3)=pma(1,5)*pavgm(5)/(1.-sw(22)*glob7s(pma(1,5)))
1244  tn3(4)=pma(1,6)*pavgm(6)/(1.-sw(22)*glob7s(pma(1,6)))
1245  tn3(5)=pma(1,7)*pavgm(7)/(1.-sw(22)*glob7s(pma(1,7)))
1246  tgn3(2)=pma(1,8)*pavgm(8)*(1.+sw(22)*glob7s(pma(1,8))) &
1247  *tn3(5)*tn3(5)/(pma(1,7)*pavgm(7))**2
1248  endif
1249  6 continue
1250  if(mass.eq.0) goto 50
1251 ! linear transition to full mixing below zn2(1)
1252  dmc=0
1253  if(alt.gt.zmix) dmc=1.-(zn2(1)-alt)/(zn2(1)-zmix)
1254  dz28=ds(3)
1255 ! ***** n2 density ****
1256  dmr=ds(3)/dm28m-1.
1257  d(3)=densm(alt,dm28m,xmm,tz,mn3,zn3,tn3,tgn3,mn2,zn2,tn2,tgn2)
1258  d(3)=d(3)*(1.+dmr*dmc)
1259 ! ***** he density ****
1260  d(1)=0
1261  if(mass.ne.4.and.mass.ne.48) goto 204
1262  dmr=ds(1)/(dz28*pdm(2,1))-1.
1263  d(1)=d(3)*pdm(2,1)*(1.+dmr*dmc)
1264  204 continue
1265 ! **** o density ****
1266  d(2)=0
1267  d(9)=0
1268  216 continue
1269 ! ***** o2 density ****
1270  d(4)=0
1271  if(mass.ne.32.and.mass.ne.48) goto 232
1272  dmr=ds(4)/(dz28*pdm(2,4))-1.
1273  d(4)=d(3)*pdm(2,4)*(1.+dmr*dmc)
1274  232 continue
1275 ! ***** ar density ****
1276  d(5)=0
1277  if(mass.ne.40.and.mass.ne.48) goto 240
1278  dmr=ds(5)/(dz28*pdm(2,5))-1.
1279  d(5)=d(3)*pdm(2,5)*(1.+dmr*dmc)
1280  240 continue
1281 ! ***** hydrogen density ****
1282  d(7)=0
1283 ! ***** atomic nitrogen density ****
1284  d(8)=0
1285 !
1286 ! total mass density
1287 !
1288  if(mass.eq.48) then
1289  d(6) = 1.66e-24*(4.*d(1)+16.*d(2)+28.*d(3)+32.*d(4)+40.*d(5)+ &
1290  d(7)+14.*d(8))
1291  if(imr.eq.1) d(6)=d(6)/1000.
1292  endif
1293  t(2)=tz
1294  10 continue
1295  goto 90
1296  50 continue
1297  dd=densm(alt,1.,0.,tz,mn3,zn3,tn3,tgn3,mn2,zn2,tn2,tgn2)
1298  t(2)=tz
1299  90 continue
1300  alast=alt
1301  return
1302  end subroutine gtd7
1303 
1304 !-----------------------------------------------------------------------
1357  subroutine gtd7d(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
1358  use wam_gtd7bk_mod, only: imr
1359 !
1360  dimension d(9),t(2),ap(7)
1361  call gtd7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
1362 ! total mass density
1363 !
1364  if(mass.eq.48) then
1365  d(6) = 1.66e-24*(4.*d(1)+16.*d(2)+28.*d(3)+32.*d(4)+40.*d(5)+ &
1366  d(7)+14.*d(8)+16.*d(9))
1367  if(imr.eq.1) d(6)=d(6)/1000.
1368  endif
1369  return
1370  end subroutine gtd7d
1371 
1372 !-----------------------------------------------------------------------
1412  subroutine ghp7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,d,t,press)
1413  use gettemp_mod, only: gsurf,re
1414  use wam_gtd7bk_mod, only: imr
1415 !
1416  dimension d(9),t(2),ap(7)
1417  save
1418  data bm/1.3806e-19/,rgas/831.4/
1419  data test/.00043/,ltest/12/
1420 !! bm=1.3806e-19; rgas=831.4
1421 !! test=.00043; ltest=12
1422  pl=alog10(press)
1423 ! initial altitude estimate
1424  if(pl.ge.-5.) then
1425  if(pl.gt.2.5) zi=18.06*(3.00-pl)
1426  if(pl.gt..75.and.pl.le.2.5) zi=14.98*(3.08-pl)
1427  if(pl.gt.-1..and.pl.le..75) zi=17.8*(2.72-pl)
1428  if(pl.gt.-2..and.pl.le.-1.) zi=14.28*(3.64-pl)
1429  if(pl.gt.-4..and.pl.le.-2.) zi=12.72*(4.32-pl)
1430  if(pl.le.-4.) zi=25.3*(.11-pl)
1431  iday=mod(iyd,1000)
1432  cl=glat/90.
1433  cl2=cl*cl
1434  if(iday.lt.182) cd=1.-iday/91.25
1435  if(iday.ge.182) cd=iday/91.25-3.
1436  ca=0
1437  if(pl.gt.-1.11.and.pl.le.-.23) ca=1.0
1438  if(pl.gt.-.23) ca=(2.79-pl)/(2.79+.23)
1439  if(pl.le.-1.11.and.pl.gt.-3.) ca=(-2.93-pl)/(-2.93+1.11)
1440  z=zi-4.87*cl*cd*ca-1.64*cl2*ca+.31*ca*cl
1441  endif
1442  if(pl.lt.-5.) z=22.*(pl+4.)**2+110
1443 ! iteration loop
1444  l=0
1445  10 continue
1446  l=l+1
1447  call gtd7(iyd,sec,z,glat,glong,stl,f107a,f107,ap,48,d,t)
1448  xn=d(1)+d(2)+d(3)+d(4)+d(5)+d(7)+d(8)
1449  p=bm*xn*t(2)
1450  if(imr.eq.1) p=p*1.e-6
1451  diff=pl-alog10(p)
1452  if(abs(diff).lt.test .or. l.eq.ltest) goto 20
1453  xm=d(6)/xn/1.66e-24
1454  if(imr.eq.1) xm = xm*1.e3
1455  g=gsurf/(1.+z/re)**2
1456  sh=rgas*t(2)/(xm*g)
1457 ! new altitude estimate using scale height
1458  if(l.lt.6) then
1459  z=z-sh*diff*2.302
1460  else
1461  z=z-sh*diff
1462  endif
1463  goto 10
1464  20 continue
1465  if(l.eq.ltest) write(6,100) press,diff
1466  100 format(1x,29hghp7 not converging for press, 1pe12.2,e12.2)
1467  alt=z
1468  return
1469  end subroutine ghp7
1470 
1471 !-----------------------------------------------------------------------
1479  subroutine glatf(lat,gv,reff)
1480  real lat
1481  save
1482  data dgtr/1.74533e-2/
1483 !! dgtr=1.74533e-2
1484  c2 = cos(2.*dgtr*lat)
1485  gv = 980.616*(1.-.0026373*c2)
1486  reff = 2.*gv/(3.085462e-6 + 2.27e-9*c2)*1.e-5
1487  return
1488  end subroutine glatf
1489 
1490 !-----------------------------------------------------------------------
1518  function vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,ic)
1519  use gettemp_mod, only: sw,swc
1520 !
1521  dimension ap(7),iydl(2),secl(2),glatl(2),gll(2),stll(2)
1522  dimension fal(2),fl(2),apl(7,2),swl(25,2),swcl(25,2)
1523  save
1524  data iydl/2*-999/,secl/2*-999./,glatl/2*-999./,gll/2*-999./
1525  data stll/2*-999./,fal/2*-999./,fl/2*-999./,apl/14*-999./
1526  data swl/50*-999./,swcl/50*-999./
1527  vtst7=0
1528  if(iyd.ne.iydl(ic)) goto 10
1529  if(sec.ne.secl(ic)) goto 10
1530  if(glat.ne.glatl(ic)) goto 10
1531  if(glong.ne.gll(ic)) goto 10
1532  if(stl.ne.stll(ic)) goto 10
1533  if(f107a.ne.fal(ic)) goto 10
1534  if(f107.ne.fl(ic)) goto 10
1535  do 5 i=1,7
1536  if(ap(i).ne.apl(i,ic)) goto 10
1537  5 end do
1538  do 7 i=1,25
1539  if(sw(i).ne.swl(i,ic)) goto 10
1540  if(swc(i).ne.swcl(i,ic)) goto 10
1541  7 end do
1542  goto 20
1543  10 continue
1544  vtst7=1
1545  iydl(ic)=iyd
1546  secl(ic)=sec
1547  glatl(ic)=glat
1548  gll(ic)=glong
1549  stll(ic)=stl
1550  fal(ic)=f107a
1551  fl(ic)=f107
1552  do 15 i=1,7
1553  apl(i,ic)=ap(i)
1554  15 end do
1555  do 16 i=1,25
1556  swl(i,ic)=sw(i)
1557  swcl(i,ic)=swc(i)
1558  16 end do
1559  20 continue
1560  return
1561  end function vtst7
1562 
1563 !-----------------------------------------------------------------------
1606  subroutine gts7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
1607  use wam_gtd7bk_mod, only: &
1608  ptm,pdm, &
1609  imr
1610 
1611  use gettemp_mod, only: tlb,s,db04,db16,db28,db32,db40,db48,db01,za,t0, &
1612  z0,g0,rl,dd,db14,tr12,tn1,tgn1, &
1613  pt,pd,ps,pdl,ptl,pma,sw, &
1614  dm04,dm16,dm28,dm32,dm40,dm01,dm14
1615 
1616  dimension zn1(5),alpha(9)
1617  dimension d(9),t(2),mt(11),ap(*),altl(8)
1618  save
1619  data mt/48,0,4,16,28,32,40,1,49,14,17/
1620  data altl/200.,300.,160.,250.,240.,450.,320.,450./
1621  data mn1/5/,zn1/120.,110.,100.,90.,72.5/
1622  data dgtr/1.74533e-2/,dr/1.72142e-2/,alast/-999./
1623  data alpha/-0.38,0.,0.,0.,0.17,0.,-0.38,0.,0./
1624 ! test for changed input
1625  v2=vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,2)
1626 !
1627  yrd=iyd
1628  za=pdl(16,2)
1629  zn1(1)=za
1630  do 2 j=1,9
1631  d(j)=0.
1632  2 end do
1633 ! tinf variations not important below za or zn1(1)
1634  if(alt.gt.zn1(1)) then
1635  if(v2.eq.1..or.alast.le.zn1(1)) tinf=ptm(1)*pt(1) &
1636  *(1.+sw(16)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pt))
1637  else
1638  tinf=ptm(1)*pt(1)
1639  endif
1640  t(1)=tinf
1641 ! gradient variations not important below zn1(5)
1642  if(alt.gt.zn1(5)) then
1643  if(v2.eq.1.or.alast.le.zn1(5)) g0=ptm(4)*ps(1) &
1644  *(1.+sw(19)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,ps))
1645  else
1646  g0=ptm(4)*ps(1)
1647  endif
1648 ! calculate these temperatures only if input changed
1649  if(v2.eq.1. .or. alt.lt.300.) &
1650  tlb=ptm(2)*(1.+sw(17)*globe7(yrd,sec,glat,glong,stl, &
1651  f107a,f107,ap,pd(1,4)))*pd(1,4)
1652  s=g0/(tinf-tlb)
1653 ! lower thermosphere temp variations not significant for
1654 ! density above 300 km
1655  if(alt.lt.300.) then
1656  if(v2.eq.1..or.alast.ge.300.) then
1657  tn1(2)=ptm(7)*ptl(1,1)/(1.-sw(18)*glob7s(ptl(1,1)))
1658  tn1(3)=ptm(3)*ptl(1,2)/(1.-sw(18)*glob7s(ptl(1,2)))
1659  tn1(4)=ptm(8)*ptl(1,3)/(1.-sw(18)*glob7s(ptl(1,3)))
1660  tn1(5)=ptm(5)*ptl(1,4)/(1.-sw(18)*sw(20)*glob7s(ptl(1,4)))
1661  tgn1(2)=ptm(9)*pma(1,9)*(1.+sw(18)*sw(20)*glob7s(pma(1,9))) &
1662  *tn1(5)*tn1(5)/(ptm(5)*ptl(1,4))**2
1663  endif
1664  else
1665  tn1(2)=ptm(7)*ptl(1,1)
1666  tn1(3)=ptm(3)*ptl(1,2)
1667  tn1(4)=ptm(8)*ptl(1,3)
1668  tn1(5)=ptm(5)*ptl(1,4)
1669  tgn1(2)=ptm(9)*pma(1,9) &
1670  *tn1(5)*tn1(5)/(ptm(5)*ptl(1,4))**2
1671  endif
1672 !
1673  z0=zn1(4)
1674  t0=tn1(4)
1675  tr12=1.
1676 !
1677  if(mass.eq.0) go to 50
1678 ! n2 variation factor at zlb
1679  g28=sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107, &
1680  ap,pd(1,3))
1681  day=amod(yrd,1000.)
1682 ! variation of turbopause height
1683  zhf=pdl(25,2) &
1684  *(1.+sw(5)*pdl(25,1)*sin(dgtr*glat)*cos(dr*(day-pt(14))))
1685  yrd=iyd
1686  t(1)=tinf
1687  xmm=pdm(5,3)
1688  z=alt
1689 !
1690  do 10 j = 1,11
1691  if(mass.eq.mt(j)) go to 15
1692  10 end do
1693  write(6,100) mass
1694  go to 90
1695  15 if(z.gt.altl(6).and.mass.ne.28.and.mass.ne.48) go to 17
1696 !
1697 ! **** n2 density ****
1698 !
1699 ! diffusive density at zlb
1700  db28 = pdm(1,3)*exp(g28)*pd(1,3)
1701 ! diffusive density at alt
1702  d(3)=densu(z,db28,tinf,tlb, 28.,alpha(3),t(2),ptm(6),s,mn1,zn1, &
1703  tn1,tgn1)
1704  dd=d(3)
1705 ! turbopause
1706  zh28=pdm(3,3)*zhf
1707  zhm28=pdm(4,3)*pdl(6,2)
1708  xmd=28.-xmm
1709 ! mixed density at zlb
1710  b28=densu(zh28,db28,tinf,tlb,xmd,alpha(3)-1.,tz,ptm(6),s,mn1, &
1711  zn1,tn1,tgn1)
1712  if(z.gt.altl(3).or.sw(15).eq.0.) go to 17
1713 ! mixed density at alt
1714  dm28=densu(z,b28,tinf,tlb,xmm,alpha(3),tz,ptm(6),s,mn1, &
1715  zn1,tn1,tgn1)
1716 ! net density at alt
1717  d(3)=dnet(d(3),dm28,zhm28,xmm,28.)
1718  17 continue
1719  go to(20,50,20,25,90,35,40,45,25,48,46), j
1720  20 continue
1721 !
1722 ! **** he density ****
1723 !
1724 ! density variation factor at zlb
1725  g4 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,1))
1726 ! diffusive density at zlb
1727  db04 = pdm(1,1)*exp(g4)*pd(1,1)
1728 ! diffusive density at alt
1729  d(1)=densu(z,db04,tinf,tlb, 4.,alpha(1),t(2),ptm(6),s,mn1,zn1, &
1730  tn1,tgn1)
1731  dd=d(1)
1732  if(z.gt.altl(1).or.sw(15).eq.0.) go to 24
1733 ! turbopause
1734  zh04=pdm(3,1)
1735 ! mixed density at zlb
1736  b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha(1)-1., &
1737  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1738 ! mixed density at alt
1739  dm04=densu(z,b04,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1740  zhm04=zhm28
1741 ! net density at alt
1742  d(1)=dnet(d(1),dm04,zhm04,xmm,4.)
1743 ! correction to specified mixing ratio at ground
1744  rl=alog(b28*pdm(2,1)/b04)
1745  zc04=pdm(5,1)*pdl(1,2)
1746  hc04=pdm(6,1)*pdl(2,2)
1747 ! net density corrected at alt
1748  d(1)=d(1)*ccor(z,rl,hc04,zc04)
1749  24 continue
1750  if(mass.ne.48) go to 90
1751  25 continue
1752 !
1753 ! **** o density ****
1754 !
1755 ! density variation factor at zlb
1756  g16= sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,2))
1757 ! diffusive density at zlb
1758  db16 = pdm(1,2)*exp(g16)*pd(1,2)
1759 ! diffusive density at alt
1760  d(2)=densu(z,db16,tinf,tlb, 16.,alpha(2),t(2),ptm(6),s,mn1, &
1761  zn1,tn1,tgn1)
1762  dd=d(2)
1763  if(z.gt.altl(2).or.sw(15).eq.0.) go to 34
1764 ! corrected from pdm(3,1) to pdm(3,2) 12/2/85
1765 ! turbopause
1766  zh16=pdm(3,2)
1767 ! mixed density at zlb
1768  b16=densu(zh16,db16,tinf,tlb,16-xmm,alpha(2)-1., &
1769  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1770 ! mixed density at alt
1771  dm16=densu(z,b16,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1772  zhm16=zhm28
1773 ! net density at alt
1774  d(2)=dnet(d(2),dm16,zhm16,xmm,16.)
1775 ! 3/16/99 change form to match o2 departure from diff equil near 150
1776 ! km and add dependence on f10.7
1777 ! rl=alog(b28*pdm(2,2)*abs(pdl(17,2))/b16)
1778  rl=pdm(2,2)*pdl(17,2)*(1.+sw(1)*pdl(24,1)*(f107a-150.))
1779  hc16=pdm(6,2)*pdl(4,2)
1780  zc16=pdm(5,2)*pdl(3,2)
1781  hc216=pdm(6,2)*pdl(5,2)
1782  d(2)=d(2)*ccor2(z,rl,hc16,zc16,hc216)
1783 ! chemistry correction
1784  hcc16=pdm(8,2)*pdl(14,2)
1785  zcc16=pdm(7,2)*pdl(13,2)
1786  rc16=pdm(4,2)*pdl(15,2)
1787 ! net density corrected at alt
1788  d(2)=d(2)*ccor(z,rc16,hcc16,zcc16)
1789  34 continue
1790  if(mass.ne.48.and.mass.ne.49) go to 90
1791  35 continue
1792 !
1793 ! **** o2 density ****
1794 !
1795 ! density variation factor at zlb
1796  g32= sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,5))
1797 ! diffusive density at zlb
1798  db32 = pdm(1,4)*exp(g32)*pd(1,5)
1799 ! diffusive density at alt
1800  d(4)=densu(z,db32,tinf,tlb, 32.,alpha(4),t(2),ptm(6),s,mn1, &
1801  zn1,tn1,tgn1)
1802  if(mass.eq.49) then
1803  dd=dd+2.*d(4)
1804  else
1805  dd=d(4)
1806  endif
1807  if(sw(15).eq.0.) go to 39
1808  if(z.gt.altl(4)) go to 38
1809 ! turbopause
1810  zh32=pdm(3,4)
1811 ! mixed density at zlb
1812  b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha(4)-1., &
1813  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1814 ! mixed density at alt
1815  dm32=densu(z,b32,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1816  zhm32=zhm28
1817 ! net density at alt
1818  d(4)=dnet(d(4),dm32,zhm32,xmm,32.)
1819 ! correction to specified mixing ratio at ground
1820  rl=alog(b28*pdm(2,4)/b32)
1821  hc32=pdm(6,4)*pdl(8,2)
1822  zc32=pdm(5,4)*pdl(7,2)
1823  d(4)=d(4)*ccor(z,rl,hc32,zc32)
1824  38 continue
1825 ! correction for general departure from diffusive equilibrium above
1826  hcc32=pdm(8,4)*pdl(23,2)
1827  hcc232=pdm(8,4)*pdl(23,1)
1828  zcc32=pdm(7,4)*pdl(22,2)
1829  rc32=pdm(4,4)*pdl(24,2)*(1.+sw(1)*pdl(24,1)*(f107a-150.))
1830 ! net density corrected at alt
1831  d(4)=d(4)*ccor2(z,rc32,hcc32,zcc32,hcc232)
1832  39 continue
1833  if(mass.ne.48) go to 90
1834  40 continue
1835 !
1836 ! **** ar density ****
1837 !
1838 ! density variation factor at zlb
1839  g40= sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,6))
1840 ! diffusive density at zlb
1841  db40 = pdm(1,5)*exp(g40)*pd(1,6)
1842 ! diffusive density at alt
1843  d(5)=densu(z,db40,tinf,tlb, 40.,alpha(5),t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1844  dd=d(5)
1845  if(z.gt.altl(5).or.sw(15).eq.0.) go to 44
1846 ! turbopause
1847  zh40=pdm(3,5)
1848 ! mixed density at zlb
1849  b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha(5)-1., &
1850  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1851 ! mixed density at alt
1852  dm40=densu(z,b40,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1853  zhm40=zhm28
1854 ! net density at alt
1855  d(5)=dnet(d(5),dm40,zhm40,xmm,40.)
1856 ! correction to specified mixing ratio at ground
1857  rl=alog(b28*pdm(2,5)/b40)
1858  hc40=pdm(6,5)*pdl(10,2)
1859  zc40=pdm(5,5)*pdl(9,2)
1860 ! net density corrected at alt
1861  d(5)=d(5)*ccor(z,rl,hc40,zc40)
1862  44 continue
1863  if(mass.ne.48) go to 90
1864  45 continue
1865 !
1866 ! **** hydrogen density ****
1867 !
1868 ! density variation factor at zlb
1869  g1 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,7))
1870 ! diffusive density at zlb
1871  db01 = pdm(1,6)*exp(g1)*pd(1,7)
1872 ! diffusive density at alt
1873  d(7)=densu(z,db01,tinf,tlb,1.,alpha(7),t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1874  dd=d(7)
1875  if(z.gt.altl(7).or.sw(15).eq.0.) go to 47
1876 ! turbopause
1877  zh01=pdm(3,6)
1878 ! mixed density at zlb
1879  b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha(7)-1., &
1880  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1881 ! mixed density at alt
1882  dm01=densu(z,b01,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1883  zhm01=zhm28
1884 ! net density at alt
1885  d(7)=dnet(d(7),dm01,zhm01,xmm,1.)
1886 ! correction to specified mixing ratio at ground
1887  rl=alog(b28*pdm(2,6)*abs(pdl(18,2))/b01)
1888  hc01=pdm(6,6)*pdl(12,2)
1889  zc01=pdm(5,6)*pdl(11,2)
1890  d(7)=d(7)*ccor(z,rl,hc01,zc01)
1891 ! chemistry correction
1892  hcc01=pdm(8,6)*pdl(20,2)
1893  zcc01=pdm(7,6)*pdl(19,2)
1894  rc01=pdm(4,6)*pdl(21,2)
1895 ! net density corrected at alt
1896  d(7)=d(7)*ccor(z,rc01,hcc01,zcc01)
1897  47 continue
1898  if(mass.ne.48) go to 90
1899  48 continue
1900 !
1901 ! **** atomic nitrogen density ****
1902 !
1903 ! density variation factor at zlb
1904  g14 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,8))
1905 ! diffusive density at zlb
1906  db14 = pdm(1,7)*exp(g14)*pd(1,8)
1907 ! diffusive density at alt
1908  d(8)=densu(z,db14,tinf,tlb,14.,alpha(8),t(2),ptm(6),s,mn1, &
1909  zn1,tn1,tgn1)
1910  dd=d(8)
1911  if(z.gt.altl(8).or.sw(15).eq.0.) go to 49
1912 ! turbopause
1913  zh14=pdm(3,7)
1914 ! mixed density at zlb
1915  b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha(8)-1., &
1916  t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1917 ! mixed density at alt
1918  dm14=densu(z,b14,tinf,tlb,xmm,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1919  zhm14=zhm28
1920 ! net density at alt
1921  d(8)=dnet(d(8),dm14,zhm14,xmm,14.)
1922 ! correction to specified mixing ratio at ground
1923  rl=alog(b28*pdm(2,7)*abs(pdl(3,1))/b14)
1924  hc14=pdm(6,7)*pdl(2,1)
1925  zc14=pdm(5,7)*pdl(1,1)
1926  d(8)=d(8)*ccor(z,rl,hc14,zc14)
1927 ! chemistry correction
1928  hcc14=pdm(8,7)*pdl(5,1)
1929  zcc14=pdm(7,7)*pdl(4,1)
1930  rc14=pdm(4,7)*pdl(6,1)
1931 ! net density corrected at alt
1932  d(8)=d(8)*ccor(z,rc14,hcc14,zcc14)
1933  49 continue
1934  if(mass.ne.48) go to 90
1935  46 continue
1936 !
1937 ! **** anomalous oxygen density ****
1938 !
1939  g16h = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(1,9))
1940  db16h = pdm(1,8)*exp(g16h)*pd(1,9)
1941  tho=pdm(10,8)*pdl(7,1)
1942  dd=densu(z,db16h,tho,tho,16.,alpha(9),t2,ptm(6),s,mn1, &
1943  zn1,tn1,tgn1)
1944  zsht=pdm(6,8)
1945  zmho=pdm(5,8)
1946  zsho=scalh(zmho,16.,tho)
1947  d(9)=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.))
1948  if(mass.ne.48) go to 90
1949 !
1950 ! total mass density
1951 !
1952  d(6) = 1.66e-24*(4.*d(1)+16.*d(2)+28.*d(3)+32.*d(4)+40.*d(5)+ &
1953  d(7)+14.*d(8))
1954  db48=1.66e-24*(4.*db04+16.*db16+28.*db28+32.*db32+40.*db40+db01+ &
1955  14.*db14)
1956  go to 90
1957 ! temperature at altitude
1958  50 continue
1959  z=abs(alt)
1960  ddum = densu(z,1., tinf,tlb,0.,0.,t(2),ptm(6),s,mn1,zn1,tn1,tgn1)
1961  90 continue
1962 ! adjust densities from cgs to kgm
1963  if(imr.eq.1) then
1964  do 95 i=1,9
1965  d(i)=d(i)*1.e6
1966  95 continue
1967  d(6)=d(6)/1000.
1968  endif
1969  alast=alt
1970  return
1971  100 format(1x,'mass', i5, ' not valid')
1972  end subroutine gts7
1973 
1974 !-----------------------------------------------------------------------
1978  subroutine meters(meter)
1979  use wam_gtd7bk_mod, only: imr
1980 !
1981  logical meter
1982  save
1983  imr=0
1984  if(meter) imr=1
1985  end subroutine meters
1986 
1987 !-----------------------------------------------------------------------
1995  function scalh(alt,xm,temp)
1996  use gettemp_mod, only: gsurf,re
1997 !
1998  save
1999  data rgas/831.4/
2000  g=gsurf/(1.+alt/re)**2
2001  scalh=rgas*temp/(g*xm)
2002  return
2003  end function scalh
2004 
2005 !-----------------------------------------------------------------------
2029  function globe7(yrd,sec,lat,long,tloc,f107a,f107,ap,p)
2030  use gettemp_mod, only: tinf=>tinfg,t=>tt, &
2031  sw,swc,isw, &
2032  plg,ctloc,stloc,c2tloc,s2tloc,c3tloc,s3tloc, &
2033  day,df,dfa,apd,apdf,apt,xlong,iyr
2034  real lat, long
2035  dimension p(*),sv(25),ap(*)
2036 !---- functions ------
2037 ! 3hr magnetic activity functions
2038 ! eq. a24d
2039  g0(a)=(a-4.+(p(26)-1.)*(a-4.+(exp(-abs(p(25))*(a-4.))-1.)/ &
2040  abs(p(25))))
2041 ! eq. a24c
2042  sumex(ex)=1.+(1.-ex**19)/(1.-ex)*ex**(.5)
2043 ! eq. a24a
2044  sg0(ex)=(g0(ap(2))+(g0(ap(3))*ex+g0(ap(4))*ex*ex+g0(ap(5))*ex**3 &
2045  +(g0(ap(6))*ex**4+g0(ap(7))*ex**12)*(1.-ex**8)/(1.-ex)) &
2046  )/sumex(ex)
2047 !---------------------
2048  save
2049  data dgtr/1.74533e-2/,dr/1.72142e-2/, xl/1000./,tll/1000./
2050  data sw9/1./,dayl/-1./,p14/-1000./,p18/-1000./,p32/-1000./
2051  data hr/.2618/,sr/7.2722e-5/,sv/25*1./,nsw/14/,p39/-1000./
2052  if(isw.ne.64999) call tselec(sv)
2053  do 10 j=1,14
2054  t(j)=0
2055  10 end do
2056  if(sw(9).gt.0) sw9=1.
2057  if(sw(9).lt.0) sw9=-1.
2058  iyr = yrd/1000.
2059  day = yrd - iyr*1000.
2060  xlong=long
2061 ! eq. a22 (remainder of code)
2062  if(xl.eq.lat) go to 15
2063 ! calculate legendre polynomials
2064  c = sin(lat*dgtr)
2065  s = cos(lat*dgtr)
2066  c2 = c*c
2067  c4 = c2*c2
2068  s2 = s*s
2069  plg(2,1) = c
2070  plg(3,1) = 0.5*(3.*c2 -1.)
2071  plg(4,1) = 0.5*(5.*c*c2-3.*c)
2072  plg(5,1) = (35.*c4 - 30.*c2 + 3.)/8.
2073  plg(6,1) = (63.*c2*c2*c - 70.*c2*c + 15.*c)/8.
2074  plg(7,1) = (11.*c*plg(6,1) - 5.*plg(5,1))/6.
2075 ! plg(8,1) = (13.*c*plg(7,1) - 6.*plg(6,1))/7.
2076  plg(2,2) = s
2077  plg(3,2) = 3.*c*s
2078  plg(4,2) = 1.5*(5.*c2-1.)*s
2079  plg(5,2) = 2.5*(7.*c2*c-3.*c)*s
2080  plg(6,2) = 1.875*(21.*c4 - 14.*c2 +1.)*s
2081  plg(7,2) = (11.*c*plg(6,2)-6.*plg(5,2))/5.
2082 ! plg(8,2) = (13.*c*plg(7,2)-7.*plg(6,2))/6.
2083 ! plg(9,2) = (15.*c*plg(8,2)-8.*plg(7,2))/7.
2084  plg(3,3) = 3.*s2
2085  plg(4,3) = 15.*s2*c
2086  plg(5,3) = 7.5*(7.*c2 -1.)*s2
2087  plg(6,3) = 3.*c*plg(5,3)-2.*plg(4,3)
2088  plg(7,3)=(11.*c*plg(6,3)-7.*plg(5,3))/4.
2089  plg(8,3)=(13.*c*plg(7,3)-8.*plg(6,3))/5.
2090  plg(4,4) = 15.*s2*s
2091  plg(5,4) = 105.*s2*s*c
2092  plg(6,4)=(9.*c*plg(5,4)-7.*plg(4,4))/2.
2093  plg(7,4)=(11.*c*plg(6,4)-8.*plg(5,4))/3.
2094  xl=lat
2095  15 continue
2096  if(tll.eq.tloc) go to 16
2097  if(sw(7).eq.0.and.sw(8).eq.0.and.sw(14).eq.0) goto 16
2098  stloc = sin(hr*tloc)
2099  ctloc = cos(hr*tloc)
2100  s2tloc = sin(2.*hr*tloc)
2101  c2tloc = cos(2.*hr*tloc)
2102  s3tloc = sin(3.*hr*tloc)
2103  c3tloc = cos(3.*hr*tloc)
2104  tll = tloc
2105  16 continue
2106  if(day.ne.dayl.or.p(14).ne.p14) cd14=cos(dr*(day-p(14)))
2107  if(day.ne.dayl.or.p(18).ne.p18) cd18=cos(2.*dr*(day-p(18)))
2108  if(day.ne.dayl.or.p(32).ne.p32) cd32=cos(dr*(day-p(32)))
2109  if(day.ne.dayl.or.p(39).ne.p39) cd39=cos(2.*dr*(day-p(39)))
2110  dayl = day
2111  p14 = p(14)
2112  p18 = p(18)
2113  p32 = p(32)
2114  p39 = p(39)
2115 ! f10.7 effect
2116  df = f107 - f107a
2117  dfa=f107a-150.
2118  t(1) = p(20)*df*(1.+p(60)*dfa) + p(21)*df*df + p(22)*dfa &
2119  + p(30)*dfa**2
2120  f1 = 1. + (p(48)*dfa +p(20)*df+p(21)*df*df)*swc(1)
2121  f2 = 1. + (p(50)*dfa+p(20)*df+p(21)*df*df)*swc(1)
2122 ! time independent
2123  t(2) = &
2124  (p(2)*plg(3,1) + p(3)*plg(5,1)+p(23)*plg(7,1)) &
2125  +(p(15)*plg(3,1))*dfa*swc(1) &
2126  +p(27)*plg(2,1)
2127 ! symmetrical annual
2128  t(3) = &
2129  (p(19) )*cd32
2130 ! symmetrical semiannual
2131  t(4) = &
2132  (p(16)+p(17)*plg(3,1))*cd18
2133 ! asymmetrical annual
2134  t(5) = f1* &
2135  (p(10)*plg(2,1)+p(11)*plg(4,1))*cd14
2136 ! asymmetrical semiannual
2137  t(6) = p(38)*plg(2,1)*cd39
2138 ! diurnal
2139  if(sw(7).eq.0) goto 200
2140  t71 = (p(12)*plg(3,2))*cd14*swc(5)
2141  t72 = (p(13)*plg(3,2))*cd14*swc(5)
2142  t(7) = f2* &
2143  ((p(4)*plg(2,2) + p(5)*plg(4,2) + p(28)*plg(6,2) &
2144  + t71)*ctloc &
2145  + (p(7)*plg(2,2) + p(8)*plg(4,2) +p(29)*plg(6,2) &
2146  + t72)*stloc)
2147  200 continue
2148 ! semidiurnal
2149  if(sw(8).eq.0) goto 210
2150  t81 = (p(24)*plg(4,3)+p(36)*plg(6,3))*cd14*swc(5)
2151  t82 = (p(34)*plg(4,3)+p(37)*plg(6,3))*cd14*swc(5)
2152  t(8) = f2* &
2153  ((p(6)*plg(3,3) + p(42)*plg(5,3) + t81)*c2tloc &
2154  +(p(9)*plg(3,3) + p(43)*plg(5,3) + t82)*s2tloc)
2155  210 continue
2156 ! terdiurnal
2157  if(sw(14).eq.0) goto 220
2158  t(14) = f2* &
2159  ((p(40)*plg(4,4)+(p(94)*plg(5,4)+p(47)*plg(7,4))*cd14*swc(5))* &
2160  s3tloc &
2161  +(p(41)*plg(4,4)+(p(95)*plg(5,4)+p(49)*plg(7,4))*cd14*swc(5))* &
2162  c3tloc)
2163  220 continue
2164 ! magnetic activity based on daily ap
2165 
2166  if(sw9.eq.-1.) go to 30
2167  apd=(ap(1)-4.)
2168  p44=p(44)
2169  p45=p(45)
2170  if(p44.lt.0) p44=1.e-5
2171  apdf = apd+(p45-1.)*(apd+(exp(-p44 *apd)-1.)/p44)
2172  if(sw(9).eq.0) goto 40
2173  t(9)=apdf*(p(33)+p(46)*plg(3,1)+p(35)*plg(5,1)+ &
2174  (p(101)*plg(2,1)+p(102)*plg(4,1)+p(103)*plg(6,1))*cd14*swc(5)+ &
2175  (p(122)*plg(2,2)+p(123)*plg(4,2)+p(124)*plg(6,2))*swc(7)* &
2176  cos(hr*(tloc-p(125))))
2177  go to 40
2178  30 continue
2179  if(p(52).eq.0) go to 40
2180  exp1 = exp(-10800.*abs(p(52))/(1.+p(139)*(45.-abs(lat))))
2181  if(exp1.gt..99999) exp1=.99999
2182  if(p(25).lt.1.e-4) p(25)=1.e-4
2183  apt(1)=sg0(exp1)
2184 ! apt(2)=sg2(exp1)
2185 ! apt(3)=sg0(exp2)
2186 ! apt(4)=sg2(exp2)
2187  if(sw(9).eq.0) goto 40
2188  t(9) = apt(1)*(p(51)+p(97)*plg(3,1)+p(55)*plg(5,1)+ &
2189  (p(126)*plg(2,1)+p(127)*plg(4,1)+p(128)*plg(6,1))*cd14*swc(5)+ &
2190  (p(129)*plg(2,2)+p(130)*plg(4,2)+p(131)*plg(6,2))*swc(7)* &
2191  cos(hr*(tloc-p(132))))
2192  40 continue
2193  if(sw(10).eq.0.or.long.le.-1000.) go to 49
2194 ! longitudinal
2195  if(sw(11).eq.0) goto 230
2196  t(11)= (1.+p(81)*dfa*swc(1))* &
2197  ((p(65)*plg(3,2)+p(66)*plg(5,2)+p(67)*plg(7,2) &
2198  +p(104)*plg(2,2)+p(105)*plg(4,2)+p(106)*plg(6,2) &
2199  +swc(5)*(p(110)*plg(2,2)+p(111)*plg(4,2)+p(112)*plg(6,2))*cd14)* &
2200  cos(dgtr*long) &
2201  +(p(91)*plg(3,2)+p(92)*plg(5,2)+p(93)*plg(7,2) &
2202  +p(107)*plg(2,2)+p(108)*plg(4,2)+p(109)*plg(6,2) &
2203  +swc(5)*(p(113)*plg(2,2)+p(114)*plg(4,2)+p(115)*plg(6,2))*cd14)* &
2204  sin(dgtr*long))
2205  230 continue
2206 ! ut and mixed ut,longitude
2207  if(sw(12).eq.0) goto 240
2208  t(12)=(1.+p(96)*plg(2,1))*(1.+p(82)*dfa*swc(1))* &
2209  (1.+p(120)*plg(2,1)*swc(5)*cd14)* &
2210  ((p(69)*plg(2,1)+p(70)*plg(4,1)+p(71)*plg(6,1))* &
2211  cos(sr*(sec-p(72))))
2212  t(12)=t(12)+swc(11)* &
2213  (p(77)*plg(4,3)+p(78)*plg(6,3)+p(79)*plg(8,3))* &
2214  cos(sr*(sec-p(80))+2.*dgtr*long)*(1.+p(138)*dfa*swc(1))
2215  240 continue
2216 ! ut,longitude magnetic activity
2217  if(sw(13).eq.0) goto 48
2218  if(sw9.eq.-1.) go to 45
2219  t(13)= apdf*swc(11)*(1.+p(121)*plg(2,1))* &
2220  ((p( 61)*plg(3,2)+p( 62)*plg(5,2)+p( 63)*plg(7,2))* &
2221  cos(dgtr*(long-p( 64)))) &
2222  +apdf*swc(11)*swc(5)* &
2223  (p(116)*plg(2,2)+p(117)*plg(4,2)+p(118)*plg(6,2))* &
2224  cd14*cos(dgtr*(long-p(119))) &
2225  + apdf*swc(12)* &
2226  (p( 84)*plg(2,1)+p( 85)*plg(4,1)+p( 86)*plg(6,1))* &
2227  cos(sr*(sec-p( 76)))
2228  goto 48
2229  45 continue
2230  if(p(52).eq.0) goto 48
2231  t(13)=apt(1)*swc(11)*(1.+p(133)*plg(2,1))* &
2232  ((p(53)*plg(3,2)+p(99)*plg(5,2)+p(68)*plg(7,2))* &
2233  cos(dgtr*(long-p(98)))) &
2234  +apt(1)*swc(11)*swc(5)* &
2235  (p(134)*plg(2,2)+p(135)*plg(4,2)+p(136)*plg(6,2))* &
2236  cd14*cos(dgtr*(long-p(137))) &
2237  +apt(1)*swc(12)* &
2238  (p(56)*plg(2,1)+p(57)*plg(4,1)+p(58)*plg(6,1))* &
2239  cos(sr*(sec-p(59)))
2240  48 continue
2241 ! parms not used: 83, 90,100,140-150
2242  49 continue
2243  tinf=p(31)
2244  do i = 1,nsw
2245  tinf = tinf + abs(sw(i))*t(i)
2246  enddo
2247  globe7 = tinf
2248  return
2249  end function globe7
2250 
2251 !-----------------------------------------------------------------------
2263  subroutine tselec(sv)
2264  use gettemp_mod, only: sw,swc,isw
2265 !
2266  dimension sv(*),sav(25),svv(*)
2267  save
2268  do 100 i = 1,25
2269  sav(i)=sv(i)
2270  sw(i)=amod(sv(i),2.)
2271  if(abs(sv(i)).eq.1.or.abs(sv(i)).eq.2.) then
2272  swc(i)=1.
2273  else
2274  swc(i)=0.
2275  endif
2276  100 end do
2277  isw=64999
2278  return
2279  entry tretrv(svv)
2280  do 200 i=1,25
2281  svv(i)=sav(i)
2282  200 end do
2283  end subroutine tselec
2284 
2285 !-----------------------------------------------------------------------
2292  function glob7s(p)
2293  use mpi
2294  use gettemp_mod, only:plg,ctloc,stloc,c2tloc,s2tloc,c3tloc,s3tloc, &
2295  day,dfa,apdf,apt,long=>xlong,sw,swc
2296  dimension p(*),t(14)
2297 !
2298  integer :: rc
2299  save
2300  data dr/1.72142e-2/,dgtr/1.74533e-2/,pset/2./
2301  data dayl/-1./,p32,p18,p14,p39/4*-1000./
2302  if(p(100).eq.0) p(100)=pset
2303  if(p(100).ne.pset) then
2304  write(6,900) pset,p(100)
2305  900 format(1x,'FATAL ERROR: Wrong parameter set for glob7s',3f10.1)
2306  call mpi_abort(mpi_comm_world, 999, rc)
2307  endif
2308  do 10 j=1,14
2309  t(j)=0.
2310  10 end do
2311  if(day.ne.dayl.or.p32.ne.p(32)) cd32=cos(dr*(day-p(32)))
2312  if(day.ne.dayl.or.p18.ne.p(18)) cd18=cos(2.*dr*(day-p(18)))
2313  if(day.ne.dayl.or.p14.ne.p(14)) cd14=cos(dr*(day-p(14)))
2314  if(day.ne.dayl.or.p39.ne.p(39)) cd39=cos(2.*dr*(day-p(39)))
2315  dayl=day
2316  p32=p(32)
2317  p18=p(18)
2318  p14=p(14)
2319  p39=p(39)
2320 !
2321 ! f10.7
2322  t(1)=p(22)*dfa
2323 ! time independent
2324  t(2)=p(2)*plg(3,1)+p(3)*plg(5,1)+p(23)*plg(7,1) &
2325  +p(27)*plg(2,1)+p(15)*plg(4,1)+p(60)*plg(6,1)
2326 ! symmetrical annual
2327  t(3)=(p(19)+p(48)*plg(3,1)+p(30)*plg(5,1))*cd32
2328 ! symmetrical semiannual
2329  t(4)=(p(16)+p(17)*plg(3,1)+p(31)*plg(5,1))*cd18
2330 ! asymmetrical annual
2331  t(5)=(p(10)*plg(2,1)+p(11)*plg(4,1)+p(21)*plg(6,1))*cd14
2332 ! asymmetrical semiannual
2333  t(6)=(p(38)*plg(2,1))*cd39
2334 ! diurnal
2335  if(sw(7).eq.0) goto 200
2336  t71 = p(12)*plg(3,2)*cd14*swc(5)
2337  t72 = p(13)*plg(3,2)*cd14*swc(5)
2338  t(7) = &
2339  ((p(4)*plg(2,2) + p(5)*plg(4,2) &
2340  + t71)*ctloc &
2341  + (p(7)*plg(2,2) + p(8)*plg(4,2) &
2342  + t72)*stloc)
2343  200 continue
2344 ! semidiurnal
2345  if(sw(8).eq.0) goto 210
2346  t81 = (p(24)*plg(4,3)+p(36)*plg(6,3))*cd14*swc(5)
2347  t82 = (p(34)*plg(4,3)+p(37)*plg(6,3))*cd14*swc(5)
2348  t(8) = &
2349  ((p(6)*plg(3,3) + p(42)*plg(5,3) + t81)*c2tloc &
2350  +(p(9)*plg(3,3) + p(43)*plg(5,3) + t82)*s2tloc)
2351  210 continue
2352 ! terdiurnal
2353  if(sw(14).eq.0) goto 220
2354  t(14) = p(40)*plg(4,4)*s3tloc +p(41)*plg(4,4)*c3tloc
2355  220 continue
2356 ! magnetic activity
2357  if(sw(9).eq.0) goto 40
2358  if(sw(9).eq.1) &
2359  t(9)=apdf*(p(33)+p(46)*plg(3,1)*swc(2))
2360  if(sw(9).eq.-1) &
2361  t(9)=(p(51)*apt(1)+p(97)*plg(3,1)*apt(1)*swc(2))
2362  40 continue
2363  if(sw(10).eq.0.or.sw(11).eq.0.or.long.le.-1000.) go to 49
2364 ! longitudinal
2365  t(11)= (1.+plg(2,1)*(p(81)*swc(5)*cos(dr*(day-p(82))) &
2366  +p(86)*swc(6)*cos(2.*dr*(day-p(87)))) &
2367  +p(84)*swc(3)*cos(dr*(day-p(85))) &
2368  +p(88)*swc(4)*cos(2.*dr*(day-p(89)))) &
2369  *((p(65)*plg(3,2)+p(66)*plg(5,2)+p(67)*plg(7,2) &
2370  +p(75)*plg(2,2)+p(76)*plg(4,2)+p(77)*plg(6,2) &
2371  )*cos(dgtr*long) &
2372  +(p(91)*plg(3,2)+p(92)*plg(5,2)+p(93)*plg(7,2) &
2373  +p(78)*plg(2,2)+p(79)*plg(4,2)+p(80)*plg(6,2) &
2374  )*sin(dgtr*long))
2375  49 continue
2376  tt=0.
2377  do i=1,14
2378  tt=tt+abs(sw(i))*t(i)
2379  enddo
2380  glob7s=tt
2381  return
2382  end function glob7s
2383 
2384 !--------------------------------------------------------------------
2404  function densu(alt,dlb,tinf,tlb,xm,alpha,tz,zlb,s2, &
2405  mn1,zn1,tn1,tgn1)
2406  use gettemp_mod, only: gsurf,re
2407 !
2408  dimension zn1(mn1),tn1(mn1),tgn1(2),xs(5),ys(5),y2out(5)
2409 !function
2410  zeta(zz,zl)=(zz-zl)*(re+zl)/(re+zz)
2411  save
2412  data rgas/831.4/
2413 !! rgas=831.4
2414 !cccccwrite(6,*) 'db',alt,dlb,tinf,tlb,xm,alpha,zlb,s2,mn1,zn1,tn1
2415  densu=1.
2416 ! joining altitude of bates and spline
2417  za=zn1(1)
2418  z=amax1(alt,za)
2419 ! geopotential altitude difference from zlb
2420  zg2=zeta(z,zlb)
2421 ! bates temperature
2422  tt=tinf-(tinf-tlb)*exp(-s2*zg2)
2423  ta=tt
2424  tz=tt
2425  densu=tz
2426  if(alt.ge.za) go to 10
2427 !
2428 ! calculate temperature below za
2429 ! temperature gradient at za from bates profile
2430  dta=(tinf-ta)*s2*((re+zlb)/(re+za))**2
2431  tgn1(1)=dta
2432  tn1(1)=ta
2433  z=amax1(alt,zn1(mn1))
2434  mn=mn1
2435  z1=zn1(1)
2436  z2=zn1(mn)
2437  t1=tn1(1)
2438  t2=tn1(mn)
2439 ! geopotental difference from z1
2440  zg=zeta(z,z1)
2441  zgdif=zeta(z2,z1)
2442 ! set up spline nodes
2443  do 20 k=1,mn
2444  xs(k)=zeta(zn1(k),z1)/zgdif
2445  ys(k)=1./tn1(k)
2446  20 end do
2447 ! end node derivatives
2448  yd1=-tgn1(1)/(t1*t1)*zgdif
2449  yd2=-tgn1(2)/(t2*t2)*zgdif*((re+z2)/(re+z1))**2
2450 ! calculate spline coefficients
2451  call spline(xs,ys,mn,yd1,yd2,y2out)
2452  x=zg/zgdif
2453  call splint(xs,ys,y2out,mn,x,y)
2454 ! temperature at altitude
2455  tz=1./y
2456  densu=tz
2457  10 if(xm.eq.0.) go to 50
2458 !
2459 ! calculate density above za
2460  glb=gsurf/(1.+zlb/re)**2
2461  gamma=xm*glb/(s2*rgas*tinf)
2462  expl=exp(-s2*gamma*zg2)
2463  if(expl.gt.50.or.tt.le.0.) then
2464  expl=50.
2465  endif
2466 ! density at altitude
2467  densa=dlb*(tlb/tt)**(1.+alpha+gamma)*expl
2468  densu=densa
2469  if(alt.ge.za) go to 50
2470 !
2471 ! calculate density below za
2472  glb=gsurf/(1.+z1/re)**2
2473  gamm=xm*glb*zgdif/rgas
2474 ! integrate spline temperatures
2475  call splini(xs,ys,y2out,mn,x,yi)
2476  expl=gamm*yi
2477  if(expl.gt.50..or.tz.le.0.) then
2478  expl=50.
2479  endif
2480 ! density at altitude
2481  densu=densu*(t1/tz)**(1.+alpha)*exp(-expl)
2482  50 continue
2483  return
2484  end function densu
2485 
2486 !--------------------------------------------------------------------
2504  function densm(alt,d0,xm,tz,mn3,zn3,tn3,tgn3,mn2,zn2,tn2,tgn2)
2505  use gettemp_mod, only: gsurf,re
2506 !
2507  dimension zn3(mn3),tn3(mn3),tgn3(2),xs(10),ys(10),y2out(10)
2508  dimension zn2(mn2),tn2(mn2),tgn2(2)
2509 ! function
2510  zeta(zz,zl)=(zz-zl)*(re+zl)/(re+zz)
2511  save
2512  data rgas/831.4/
2513 !! rgas=831.4
2514  densm=d0
2515  if(alt.gt.zn2(1)) goto 50
2516 ! stratosphere/mesosphere temperature
2517  z=amax1(alt,zn2(mn2))
2518  mn=mn2
2519  z1=zn2(1)
2520  z2=zn2(mn)
2521  t1=tn2(1)
2522  t2=tn2(mn)
2523  zg=zeta(z,z1)
2524  zgdif=zeta(z2,z1)
2525 ! set up spline nodes
2526  do 210 k=1,mn
2527  xs(k)=zeta(zn2(k),z1)/zgdif
2528  ys(k)=1./tn2(k)
2529  210 end do
2530  yd1=-tgn2(1)/(t1*t1)*zgdif
2531  yd2=-tgn2(2)/(t2*t2)*zgdif*((re+z2)/(re+z1))**2
2532 ! calculate spline coefficients
2533  call spline(xs,ys,mn,yd1,yd2,y2out)
2534  x=zg/zgdif
2535  call splint(xs,ys,y2out,mn,x,y)
2536 ! temperature at altitude
2537  tz=1./y
2538  if(xm.eq.0.) go to 20
2539 !
2540 ! calculate stratosphere/mesosphere density
2541  glb=gsurf/(1.+z1/re)**2
2542  gamm=xm*glb*zgdif/rgas
2543 ! integrate temperature profile
2544  call splini(xs,ys,y2out,mn,x,yi)
2545  expl=gamm*yi
2546  if(expl.gt.50.) expl=50.
2547 ! density at altitude
2548  densm=densm*(t1/tz)*exp(-expl)
2549  20 continue
2550  if(alt.gt.zn3(1)) goto 50
2551 !
2552 ! troposphere/stratosphere temperature
2553  z=alt
2554  mn=mn3
2555  z1=zn3(1)
2556  z2=zn3(mn)
2557  t1=tn3(1)
2558  t2=tn3(mn)
2559  zg=zeta(z,z1)
2560  zgdif=zeta(z2,z1)
2561 ! set up spline nodes
2562  do 220 k=1,mn
2563  xs(k)=zeta(zn3(k),z1)/zgdif
2564  ys(k)=1./tn3(k)
2565  220 end do
2566  yd1=-tgn3(1)/(t1*t1)*zgdif
2567  yd2=-tgn3(2)/(t2*t2)*zgdif*((re+z2)/(re+z1))**2
2568 ! calculate spline coefficients
2569  call spline(xs,ys,mn,yd1,yd2,y2out)
2570  x=zg/zgdif
2571  call splint(xs,ys,y2out,mn,x,y)
2572 ! temperature at altitude
2573  tz=1./y
2574  if(xm.eq.0.) go to 30
2575 !
2576 ! calculate tropospheric/stratosphere density
2577 !
2578  glb=gsurf/(1.+z1/re)**2
2579  gamm=xm*glb*zgdif/rgas
2580 ! integrate temperature profile
2581  call splini(xs,ys,y2out,mn,x,yi)
2582  expl=gamm*yi
2583  if(expl.gt.50.) expl=50.
2584 ! density at altitude
2585  densm=densm*(t1/tz)*exp(-expl)
2586  30 continue
2587  50 continue
2588  if(xm.eq.0) densm=tz
2589  return
2590  end function densm
2591 
2592 !-----------------------------------------------------------------------
2605  subroutine spline(x,y,n,yp1,ypn,y2)
2606  parameter(nmax=100)
2607  dimension x(n),y(n),y2(n),u(nmax)
2608  save
2609  if(yp1.gt..99e30) then
2610  y2(1)=0
2611  u(1)=0
2612  else
2613  y2(1)=-.5
2614  u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
2615  endif
2616  do 11 i=2,n-1
2617  sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
2618  p=sig*y2(i-1)+2.
2619  y2(i)=(sig-1.)/p
2620  u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
2621  /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
2622  11 end do
2623  if(ypn.gt..99e30) then
2624  qn=0
2625  un=0
2626  else
2627  qn=.5
2628  un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
2629  endif
2630  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
2631  do 12 k=n-1,1,-1
2632  y2(k)=y2(k)*y2(k+1)+u(k)
2633  12 end do
2634  return
2635  end subroutine spline
2636 
2637 !-----------------------------------------------------------------------
2649  subroutine splint(xa,ya,y2a,n,x,y)
2650  dimension xa(n),ya(n),y2a(n)
2651  save
2652  klo=1
2653  khi=n
2654  1 continue
2655  if(khi-klo.gt.1) then
2656  k=(khi+klo)/2
2657  if(xa(k).gt.x) then
2658  khi=k
2659  else
2660  klo=k
2661  endif
2662  goto 1
2663  endif
2664  h=xa(khi)-xa(klo)
2665  if(h.eq.0) write(6,*) 'bad xa input to splint'
2666  a=(xa(khi)-x)/h
2667  b=(x-xa(klo))/h
2668  y=a*ya(klo)+b*ya(khi)+ &
2669  ((a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi))*h*h/6.
2670  return
2671  end subroutine splint
2672 
2673 !-----------------------------------------------------------------------
2684  subroutine splini(xa,ya,y2a,n,x,yi)
2685  dimension xa(n),ya(n),y2a(n)
2686  save
2687  yi=0
2688  klo=1
2689  khi=2
2690  1 continue
2691  if(x.gt.xa(klo).and.khi.le.n) then
2692  xx=x
2693  if(khi.lt.n) xx=amin1(x,xa(khi))
2694  h=xa(khi)-xa(klo)
2695  a=(xa(khi)-xx)/h
2696  b=(xx-xa(klo))/h
2697  a2=a*a
2698  b2=b*b
2699  yi=yi+((1.-a2)*ya(klo)/2.+b2*ya(khi)/2.+ &
2700  ((-(1.+a2*a2)/4.+a2/2.)*y2a(klo)+ &
2701  (b2*b2/4.-b2/2.)*y2a(khi))*h*h/6.)*h
2702  klo=klo+1
2703  khi=khi+1
2704  goto 1
2705  endif
2706  return
2707  end subroutine splini
2708 
2709 !-----------------------------------------------------------------------
2720  function dnet(dd,dm,zhm,xmm,xm)
2721  save
2722  a=zhm/(xmm-xm)
2723  if(dm.gt.0.and.dd.gt.0) goto 5
2724  write(6,*) 'dnet log error',dm,dd,xm
2725  if(dd.eq.0.and.dm.eq.0) dd=1.
2726  if(dm.eq.0) goto 10
2727  if(dd.eq.0) goto 20
2728  5 continue
2729  ylog=a*alog(dm/dd)
2730  if(ylog.lt.-10.) go to 10
2731  if(ylog.gt.10.) go to 20
2732  dnet=dd*(1.+exp(ylog))**(1/a)
2733  go to 50
2734  10 continue
2735  dnet=dd
2736  go to 50
2737  20 continue
2738  dnet=dm
2739  go to 50
2740  50 continue
2741  return
2742  end function dnet
2743 
2744 !-----------------------------------------------------------------------
2754  function ccor(alt, r,h1,zh)
2755  save
2756  e=(alt-zh)/h1
2757  if(e.gt.70.) go to 20
2758  if(e.lt.-70.) go to 10
2759  ex=exp(e)
2760  ccor=r/(1.+ex)
2761  go to 50
2762  10 ccor=r
2763  go to 50
2764  20 ccor=0.
2765  go to 50
2766  50 continue
2767  ccor=exp(ccor)
2768  return
2769  end function ccor
2770 
2771 !-----------------------------------------------------------------------
2782  function ccor2(alt, r,h1,zh,h2)
2783  e1=(alt-zh)/h1
2784  e2=(alt-zh)/h2
2785  if(e1.gt.70. .or. e2.gt.70.) go to 20
2786  if(e1.lt.-70. .and. e2.lt.-70) go to 10
2787  ex1=exp(e1)
2788  ex2=exp(e2)
2789  ccor2=r/(1.+.5*(ex1+ex2))
2790  go to 50
2791  10 ccor2=r
2792  go to 50
2793  20 ccor2=0.
2794  go to 50
2795  50 continue
2796  ccor2=exp(ccor2)
2797  return
2798  end function ccor2
densm
function densm(alt, d0, xm, tz, mn3, zn3, tn3, tgn3, mn2, zn2, tn2, tgn2)
Calculate temperature and density profiles for lower atmos.
Definition: wam_climo_data.f90:2504
gts7
subroutine gts7(iyd, sec, alt, glat, glong, stl, f107a, f107, ap, mass, d, t)
Thermospheric portion of nrlmsise-00.
Definition: wam_climo_data.f90:1606
globe7
function globe7(yrd, sec, lat, long, tloc, f107a, f107, ap, p)
Calculate g(l) function for upper thermosphere parameters.
Definition: wam_climo_data.f90:2029
ccor
function ccor(alt, r, h1, zh)
Chemistry/dissociation correction.
Definition: wam_climo_data.f90:2754
wam_gtd7bk_mod
Use moduke for blockdata gtd7bk.
Definition: wam_climo_data.f90:29
meters
subroutine meters(meter)
Convert outputs to kg & meters if meter true.
Definition: wam_climo_data.f90:1978
dnet
function dnet(dd, dm, zhm, xmm, xm)
Turbopause correction.
Definition: wam_climo_data.f90:2720
tselec
subroutine tselec(sv)
Set switches.
Definition: wam_climo_data.f90:2263
ghp7
subroutine ghp7(iyd, sec, alt, glat, glong, stl, f107a, f107, ap, d, t, press)
Find altitude of pressure surface (press) from gtd7.
Definition: wam_climo_data.f90:1412
glob7s
function glob7s(p)
Version of globe for lower atmosphere.
Definition: wam_climo_data.f90:2292
splint
subroutine splint(xa, ya, y2a, n, x, y)
Calculate cubic spline interp value.
Definition: wam_climo_data.f90:2649
vtst7
function vtst7(iyd, sec, glat, glong, stl, f107a, f107, ap, ic)
Test variable condition.
Definition: wam_climo_data.f90:1518
glatf
subroutine glatf(lat, gv, reff)
Calculate latitude variable.
Definition: wam_climo_data.f90:1479
spline
subroutine spline(x, y, n, yp1, ypn, y2)
Calculate 2nd derivatives of cubic spline interp function.
Definition: wam_climo_data.f90:2605
gettemp_mod
Use moduke for common blocks.
Definition: wam_climo_data.f90:879
densu
function densu(alt, dlb, tinf, tlb, xm, alpha, tz, zlb, s2, mn1, zn1, tn1, tgn1)
Calculate temperature and density profiles.
Definition: wam_climo_data.f90:2404
gettemp
subroutine gettemp(iday, nday, xlat, nlat, pr, np, temp, n_o, n_o2, n_n2)
Entry routine to get WAM needed temperature and composition profiles.
Definition: wam_climo_data.f90:967
gtd7d
subroutine gtd7d(iyd, sec, alt, glat, glong, stl, f107a, f107, ap, mass, d, t)
The nrlmsise-00 subroutine gtd7d.
Definition: wam_climo_data.f90:1357
ccor2
function ccor2(alt, r, h1, zh, h2)
O and O2 chemistry/dissociation correction.
Definition: wam_climo_data.f90:2782
scalh
function scalh(alt, xm, temp)
Calculate scale height (km)
Definition: wam_climo_data.f90:1995
gtd7
subroutine gtd7(iyd, sec, alt, glat, glong, stl, f107a, f107, ap, mass, d, t)
The nrlmsise-00 subroutine gtd7.
Definition: wam_climo_data.f90:1143
splini
subroutine splini(xa, ya, y2a, n, x, yi)
Integrate cubic spline function.
Definition: wam_climo_data.f90:2684