13 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
17 dimension alm(3), qms(4:6)
18 DATA alm / 0.177d0, 0.239d0, 0.247d0 /
19 DATA qmn / 1.60d0 /, (qms(
i),
i=4,6) / 1.60d0, 5.00d0, 180.0d0 /
24 IF(
q .LE. qms(ip))
THEN
35 sbl =
log(
q/alam) /
log(qi/alam)
43 a0=exp( 0.1907d+00+0.4205d-01*sb +0.2752d+00*sb2-
45 a1= 0.4611d+00+0.2331d-01*sb -0.3403d-01*sb2+0.3174d-01*sb3
46 a2= 0.3504d+01+0.5739d+00*sb +0.2676d+00*sb2-0.1553d+00*sb3
47 a3= 0.7452d+01-0.6742d+01*sb +0.2849d+01*sb2-0.1964d+00*sb3
48 a4= 0.1116d+01-0.3435d+00*sb +0.2865d+00*sb2-0.1288d+00*sb3
49 a5= 0.6659d-01+0.2714d+00*sb -0.2688d+00*sb2+0.2763d+00*sb3
50 ELSEIF(iprt .EQ. 1)
THEN
51 a0=exp( 0.1141d+00+0.4764d+00*sb -0.1745d+01*sb2+
53 a1= 0.4275d+00-0.1290d+00*sb +0.3609d+00*sb2-0.1689d+00*sb3
54 a2= 0.3000d+01+0.2946d+01*sb -0.4117d+01*sb2+0.1989d+01*sb3
55 a3=-0.1302d+01+0.2322d+01*sb -0.4258d+01*sb2+0.2109d+01*sb3
56 a4= 0.2586d+01-0.1920d+00*sb -0.3754d+00*sb2+0.2731d+00*sb3
57 a5=-0.2251d+00-0.5374d+00*sb +0.2245d+01*sb2-0.1034d+01*sb3
58 ELSEIF(iprt .EQ. 0)
THEN
59 a0=exp(-0.7631d+00-0.7241d+00*sb -0.1170d+01*sb2+
61 a1=-0.3573d+00+0.3469d+00*sb -0.3396d+00*sb2+0.9188d-01*sb3
62 a2= 0.5604d+01+0.7458d+00*sb -0.5082d+00*sb2+0.1844d+00*sb3
63 a3= 0.1549d+02-0.1809d+02*sb +0.1162d+02*sb2-0.3483d+01*sb3
64 a4= 0.9881d+00+0.1364d+00*sb -0.4421d+00*sb2+0.2051d+00*sb3
65 a5=-0.9505d-01+0.3259d+01*sb -0.1547d+01*sb2+0.2918d+00*sb3
66 ELSEIF(iprt .EQ. -1)
THEN
67 a0=exp(-0.2449d+01-0.3513d+01*sb +0.4529d+01*sb2-
69 a1=-0.4050d+00+0.3411d+00*sb -0.3669d+00*sb2+0.1109d+00*sb3
70 a2= 0.7470d+01-0.2982d+01*sb +0.5503d+01*sb2-0.2419d+01*sb3
71 a3= 0.1503d+02+0.1638d+01*sb -0.8772d+01*sb2+0.3852d+01*sb3
72 a4= 0.1137d+01-0.1006d+01*sb +0.1485d+01*sb2-0.6389d+00*sb3
73 a5=-0.5299d+00+0.3160d+01*sb -0.3104d+01*sb2+0.1219d+01*sb3
74 ELSEIF(iprt .EQ. -2)
THEN
75 a0=exp(-0.2740d+01-0.7987d-01*sb -0.9015d+00*sb2-
77 a1=-0.3909d+00+0.1244d+00*sb -0.4487d-01*sb2+0.1277d-01*sb3
78 a2= 0.9163d+01+0.2823d+00*sb -0.7720d+00*sb2-0.9360d-02*sb3
79 a3= 0.1080d+02-0.3915d+01*sb -0.1153d+01*sb2+0.2649d+01*sb3
80 a4= 0.9894d+00-0.1647d+00*sb -0.9426d-02*sb2+0.2945d-02*sb3
81 a5=-0.3395d+00+0.6998d+00*sb +0.7000d+00*sb2-0.6730d-01*sb3
82 ELSEIF(iprt .EQ. -3)
THEN
83 a0=exp(-0.3640d+01+0.1250d+01*sb -0.2914d+01*sb2+
85 a1=-0.3595d+00-0.5259d-01*sb +0.3122d+00*sb2-0.1642d+00*sb3
86 a2= 0.7305d+01+0.9727d+00*sb -0.9788d+00*sb2-0.5193d-01*sb3
87 a3= 0.1198d+02-0.1799d+02*sb +0.2614d+02*sb2-0.1091d+02*sb3
88 a4= 0.9882d+00-0.6101d+00*sb +0.9737d+00*sb2-0.4935d+00*sb3
89 a5=-0.1186d+00-0.3231d+00*sb +0.3074d+01*sb2-0.1274d+01*sb3
90 ELSEIF(iprt .EQ. -4)
THEN
91 a0=sb** 0.1122d+01*exp(-0.3718d+01-0.1335d+01*sb +
93 a1=-0.4719d+00+0.7509d+00*sb -0.8420d+00*sb2+0.2901d+00*sb3
94 a2= 0.6194d+01-0.1641d+01*sb +0.4907d+01*sb2-0.2523d+01*sb3
95 a3= 0.4426d+01-0.4270d+01*sb +0.6581d+01*sb2-0.3474d+01*sb3
96 a4= 0.2683d+00+0.9876d+00*sb -0.7612d+00*sb2+0.1780d+00*sb3
97 a5=-0.4547d+00+0.4410d+01*sb -0.3712d+01*sb2+0.1245d+01*sb3
98 ELSEIF(iprt .EQ. -5)
THEN
99 a0=sb** 0.9838d+00*exp(-0.2548d+01-0.7660d+01*sb +
101 a1=-0.3122d+00-0.2120d+00*sb +0.5716d+00*sb2-0.3773d+00*sb3
102 a2= 0.6257d+01-0.8214d-01*sb -0.2537d+01*sb2+0.2981d+01*sb3
103 a3=-0.6723d+00+0.2131d+01*sb +0.9599d+01*sb2-0.7910d+01*sb3
104 a4= 0.9169d-01+0.4295d-01*sb -0.5017d+00*sb2+0.3811d+00*sb3
105 a5= 0.2402d+00+0.2656d+01*sb -0.1586d+01*sb2+0.2880d+00*sb3
106 ELSEIF(iprt .EQ. -6)
THEN
107 a0=sb** 0.1001d+01*exp(-0.6934d+01+0.3050d+01*sb -
109 a1=-0.1713d+00-0.5167d+00*sb +0.1241d+01*sb2-0.1703d+01*sb3
110 a2= 0.6169d+01+0.3023d+01*sb -0.1972d+02*sb2+0.1069d+02*sb3
111 a3= 0.4439d+01-0.1746d+02*sb +0.1225d+02*sb2+0.8350d+00*sb3
112 a4= 0.5458d+00-0.4586d+00*sb +0.9089d+00*sb2-0.4049d+00*sb3
113 a5= 0.3207d+01-0.3362d+01*sb +0.5877d+01*sb2-0.7659d+01*sb3
117 ELSEIF(iset .EQ. 2)
THEN
119 a0=exp( 0.2259d+00+0.1237d+00*sb +0.3035d+00*sb2-
121 a1= 0.5085d+00+0.1651d-01*sb -0.3592d-01*sb2+0.2782d-01*sb3
122 a2= 0.3732d+01+0.4901d+00*sb +0.2218d+00*sb2-0.1116d+00*sb3
123 a3= 0.7011d+01-0.6620d+01*sb +0.2557d+01*sb2-0.1360d+00*sb3
124 a4= 0.8969d+00-0.2429d+00*sb +0.1811d+00*sb2-0.6888d-01*sb3
125 a5= 0.8636d-01+0.2558d+00*sb -0.3082d+00*sb2+0.2535d+00*sb3
126 ELSEIF(iprt .EQ. 1)
THEN
127 a0=exp(-0.7266d+00-0.1584d+01*sb +0.1259d+01*sb2-
129 a1= 0.5285d+00-0.3721d+00*sb +0.5150d+00*sb2-0.1697d+00*sb3
130 a2= 0.4075d+01+0.8282d+00*sb -0.4496d+00*sb2+0.2107d+00*sb3
131 a3= 0.3279d+01+0.5066d+01*sb -0.9134d+01*sb2+0.2897d+01*sb3
132 a4= 0.4399d+00-0.5888d+00*sb +0.4802d+00*sb2-0.1664d+00*sb3
133 a5= 0.3678d+00-0.8929d+00*sb +0.1592d+01*sb2-0.5713d+00*sb3
134 ELSEIF(iprt .EQ. 0)
THEN
135 a0=exp(-0.2318d+00-0.9779d+00*sb -0.3783d+00*sb2+
137 a1=-0.2916d+00+0.1754d+00*sb -0.1884d+00*sb2+0.6116d-01*sb3
138 a2= 0.5349d+01+0.7460d+00*sb +0.2319d+00*sb2-0.2622d+00*sb3
139 a3= 0.6920d+01-0.3454d+01*sb +0.2027d+01*sb2-0.7626d+00*sb3
140 a4= 0.1013d+01+0.1423d+00*sb -0.1798d+00*sb2+0.1872d-01*sb3
141 a5=-0.5465d-01+0.2303d+01*sb -0.9584d+00*sb2+0.3098d+00*sb3
142 ELSEIF(iprt .EQ. -1)
THEN
143 a0=exp(-0.2328d+01-0.3061d+01*sb +0.3620d+01*sb2-
145 a1=-0.3358d+00+0.3198d+00*sb -0.4210d+00*sb2+0.1571d+00*sb3
146 a2= 0.8478d+01-0.3112d+01*sb +0.5243d+01*sb2-0.2255d+01*sb3
147 a3= 0.1971d+02+0.3389d+00*sb -0.5268d+01*sb2+0.2099d+01*sb3
148 a4= 0.1128d+01-0.4701d+00*sb +0.7779d+00*sb2-0.3506d+00*sb3
149 a5=-0.4708d+00+0.3341d+01*sb -0.3375d+01*sb2+0.1353d+01*sb3
150 ELSEIF(iprt .EQ. -2)
THEN
151 a0=exp(-0.2906d+01-0.1069d+00*sb -0.1055d+01*sb2+
153 a1=-0.2875d+00+0.6571d-01*sb -0.1987d-01*sb2-0.1800d-02*sb3
154 a2= 0.9854d+01-0.2715d+00*sb -0.7407d+00*sb2+0.2888d+00*sb3
155 a3= 0.1583d+02-0.7687d+01*sb +0.3428d+01*sb2-0.3327d+00*sb3
156 a4= 0.9763d+00+0.7599d-01*sb -0.2128d+00*sb2+0.6852d-01*sb3
157 a5=-0.8444d-02+0.9434d+00*sb +0.4152d+00*sb2-0.1481d+00*sb3
158 ELSEIF(iprt .EQ. -3)
THEN
159 a0=exp(-0.3780d+01+0.2499d+01*sb -0.4962d+01*sb2+
161 a1=-0.2639d+00-0.1575d+00*sb +0.3584d+00*sb2-0.1646d+00*sb3
162 a2= 0.8082d+01+0.2794d+01*sb -0.5438d+01*sb2+0.2321d+01*sb3
163 a3= 0.1811d+02-0.2000d+02*sb +0.1951d+02*sb2-0.6904d+01*sb3
164 a4= 0.9822d+00+0.4972d+00*sb -0.8690d+00*sb2+0.3415d+00*sb3
165 a5= 0.1772d+00-0.6078d+00*sb +0.3341d+01*sb2-0.1473d+01*sb3
166 ELSEIF(iprt .EQ. -4)
THEN
167 a0=sb** 0.1122d+01*exp(-0.4232d+01-0.1808d+01*sb +
169 a1=-0.2824d+00+0.5846d+00*sb -0.7230d+00*sb2+0.2419d+00*sb3
170 a2= 0.5683d+01-0.2948d+01*sb +0.5916d+01*sb2-0.2560d+01*sb3
171 a3= 0.2051d+01+0.4795d+01*sb -0.4271d+01*sb2+0.4174d+00*sb3
172 a4= 0.1737d+00+0.1717d+01*sb -0.1978d+01*sb2+0.6643d+00*sb3
173 a5= 0.8689d+00+0.3500d+01*sb -0.3283d+01*sb2+0.1026d+01*sb3
174 ELSEIF(iprt .EQ. -5)
THEN
175 a0=sb** 0.9906d+00*exp(-0.1496d+01-0.6576d+01*sb +
177 a1=-0.2140d+00-0.6419d-01*sb -0.2741d-02*sb2+0.3185d-02*sb3
178 a2= 0.5781d+01+0.1049d+00*sb -0.3930d+00*sb2+0.5174d+00*sb3
179 a3=-0.9420d+00+0.5511d+00*sb +0.8817d+00*sb2+0.1903d+01*sb3
180 a4= 0.2418d-01+0.4232d-01*sb -0.1244d-01*sb2-0.2365d-01*sb3
181 a5= 0.7664d+00+0.1794d+01*sb -0.4917d+00*sb2-0.1284d+00*sb3
182 ELSEIF(iprt .EQ. -6)
THEN
183 a0=sb** 0.1000d+01*exp(-0.8460d+01+0.1154d+01*sb +
185 a1=-0.4316d-01-0.2976d+00*sb +0.3174d+00*sb2-0.1429d+01*sb3
186 a2= 0.4910d+01+0.2273d+01*sb +0.5631d+01*sb2-0.1994d+02*sb3
187 a3= 0.1190d+02-0.2000d+02*sb -0.2000d+02*sb2+0.1292d+02*sb3
188 a4= 0.5771d+00-0.2552d+00*sb +0.7510d+00*sb2+0.6923d+00*sb3
189 a5= 0.4402d+01-0.1627d+01*sb -0.2085d+01*sb2-0.6737d+01*sb3
193 ELSEIF(iset .EQ. 3)
THEN
195 a0=exp( 0.2148d+00+0.5814d-01*sb +0.2734d+00*sb2-
197 a1= 0.4810d+00+0.1657d-01*sb -0.3800d-01*sb2+0.3125d-01*sb3
198 a2= 0.3509d+01+0.3923d+00*sb +0.4010d+00*sb2-0.1932d+00*sb3
199 a3= 0.7055d+01-0.6552d+01*sb +0.3466d+01*sb2-0.5657d+00*sb3
200 a4= 0.1061d+01-0.3453d+00*sb +0.4089d+00*sb2-0.1817d+00*sb3
201 a5= 0.8687d-01+0.2548d+00*sb -0.2967d+00*sb2+0.2647d+00*sb3
202 ELSEIF(iprt .EQ. 1)
THEN
203 a0=exp( 0.3961d+00+0.4914d+00*sb -0.1728d+01*sb2+
205 a1= 0.4162d+00-0.1419d+00*sb +0.3680d+00*sb2-0.1618d+00*sb3
206 a2= 0.3248d+01+0.3028d+01*sb -0.4307d+01*sb2+0.1920d+01*sb3
207 a3=-0.1100d+01+0.2184d+01*sb -0.3820d+01*sb2+0.1717d+01*sb3
208 a4= 0.2082d+01-0.2756d+00*sb +0.3043d+00*sb2-0.1260d+00*sb3
209 a5=-0.4822d+00-0.5706d+00*sb +0.2243d+01*sb2-0.9760d+00*sb3
210 ELSEIF(iprt .EQ. 0)
THEN
211 a0=exp(-0.4665d+00-0.7554d+00*sb -0.3323d+00*sb2-
213 a1=-0.3359d+00+0.2395d+00*sb -0.2377d+00*sb2+0.7059d-01*sb3
214 a2= 0.5451d+01+0.6086d+00*sb +0.8606d-01*sb2-0.1425d+00*sb3
215 a3= 0.1026d+02-0.9352d+01*sb +0.4879d+01*sb2-0.1150d+01*sb3
216 a4= 0.9935d+00-0.5017d-01*sb -0.1707d-01*sb2-0.1464d-02*sb3
217 a5=-0.4160d-01+0.2305d+01*sb -0.1063d+01*sb2+0.3211d+00*sb3
218 ELSEIF(iprt .EQ. -1)
THEN
219 a0=exp(-0.2714d+01-0.2868d+01*sb +0.3700d+01*sb2-
221 a1=-0.3893d+00+0.3341d+00*sb -0.3897d+00*sb2+0.1420d+00*sb3
222 a2= 0.8359d+01-0.3267d+01*sb +0.5327d+01*sb2-0.2245d+01*sb3
223 a3= 0.2359d+02-0.5669d+01*sb -0.4602d+01*sb2+0.3153d+01*sb3
224 a4= 0.1106d+01-0.4745d+00*sb +0.7739d+00*sb2-0.3417d+00*sb3
225 a5=-0.5557d+00+0.3433d+01*sb -0.3390d+01*sb2+0.1354d+01*sb3
226 ELSEIF(iprt .EQ. -2)
THEN
227 a0=exp(-0.3323d+01+0.2296d+00*sb -0.1109d+01*sb2+
229 a1=-0.3410d+00+0.8847d-01*sb -0.1111d-01*sb2-0.5927d-02*sb3
230 a2= 0.9753d+01-0.5182d+00*sb -0.4670d+00*sb2+0.1921d+00*sb3
231 a3= 0.1977d+02-0.1600d+02*sb +0.9481d+01*sb2-0.1864d+01*sb3
232 a4= 0.9818d+00+0.2839d-02*sb -0.1188d+00*sb2+0.3584d-01*sb3
233 a5=-0.7934d-01+0.1004d+01*sb +0.3704d+00*sb2-0.1220d+00*sb3
234 ELSEIF(iprt .EQ. -3)
THEN
235 a0=exp(-0.3985d+01+0.2855d+01*sb -0.5208d+01*sb2+
237 a1=-0.3337d+00-0.1150d+00*sb +0.3691d+00*sb2-0.1709d+00*sb3
238 a2= 0.7968d+01+0.3641d+01*sb -0.6599d+01*sb2+0.2642d+01*sb3
239 a3= 0.1873d+02-0.1999d+02*sb +0.1734d+02*sb2-0.5813d+01*sb3
240 a4= 0.9731d+00+0.5082d+00*sb -0.8780d+00*sb2+0.3231d+00*sb3
241 a5=-0.5542d-01-0.4189d+00*sb +0.3309d+01*sb2-0.1439d+01*sb3
242 ELSEIF(iprt .EQ. -4)
THEN
243 a0=sb** 0.1105d+01*exp(-0.3952d+01-0.1901d+01*sb +
245 a1=-0.3543d+00+0.6055d+00*sb -0.6941d+00*sb2+0.2278d+00*sb3
246 a2= 0.5955d+01-0.2629d+01*sb +0.5337d+01*sb2-0.2300d+01*sb3
247 a3= 0.1933d+01+0.4882d+01*sb -0.3810d+01*sb2+0.2290d+00*sb3
248 a4= 0.1806d+00+0.1655d+01*sb -0.1893d+01*sb2+0.6395d+00*sb3
249 a5= 0.4790d+00+0.3612d+01*sb -0.3152d+01*sb2+0.9684d+00*sb3
250 ELSEIF(iprt .EQ. -5)
THEN
251 a0=sb** 0.9818d+00*exp(-0.1825d+01-0.7464d+01*sb +
253 a1=-0.2604d+00-0.1400d+00*sb +0.1702d+00*sb2-0.8476d-01*sb3
254 a2= 0.6005d+01+0.6275d+00*sb -0.2535d+01*sb2+0.2219d+01*sb3
255 a3=-0.9067d+00+0.1149d+01*sb +0.1974d+01*sb2+0.4716d+01*sb3
256 a4= 0.3915d-01+0.5945d-01*sb -0.9844d-01*sb2+0.2783d-01*sb3
257 a5= 0.5500d+00+0.1994d+01*sb -0.6727d+00*sb2-0.1510d+00*sb3
258 ELSEIF(iprt .EQ. -6)
THEN
259 a0=sb** 0.1002d+01*exp(-0.8553d+01+0.3793d+00*sb +
261 a1=-0.5870d-01-0.2792d+00*sb +0.6526d+00*sb2-0.1984d+01*sb3
262 a2= 0.4716d+01+0.4473d+00*sb +0.1128d+02*sb2-0.1937d+02*sb3
263 a3= 0.1289d+02-0.1742d+02*sb -0.1983d+02*sb2-0.9274d+00*sb3
264 a4= 0.5647d+00-0.2732d+00*sb +0.1074d+01*sb2+0.5981d+00*sb3
265 a5= 0.4390d+01-0.1262d+01*sb -0.9026d+00*sb2-0.9394d+01*sb3
270 pycteq =
max(0d0, a0 *(
x**a1) *((1d0-
x)**a2) *(1d0+a3*(
x**a4))
271 & *(
log(1d0+1d0/
x))**a5 )