| 1 | c ************************************************************ |
|---|
| 2 | c |
|---|
| 3 | c chapman function for grazing incidence |
|---|
| 4 | c |
|---|
| 5 | c ************************************************************ |
|---|
| 6 | |
|---|
| 7 | function ch(x,z) |
|---|
| 8 | |
|---|
| 9 | integer ierr |
|---|
| 10 | real ch,x |
|---|
| 11 | |
|---|
| 12 | double precision ch1,sublim,z,error,aerr,rerr |
|---|
| 13 | double precision dcadre,fch |
|---|
| 14 | external fch,dcadre |
|---|
| 15 | |
|---|
| 16 | common/integ/ax,az |
|---|
| 17 | double precision ax,az |
|---|
| 18 | |
|---|
| 19 | c |
|---|
| 20 | |
|---|
| 21 | ax=x |
|---|
| 22 | az=z |
|---|
| 23 | |
|---|
| 24 | aerr= 0.0d0 |
|---|
| 25 | rerr= 0.0001d0 |
|---|
| 26 | sublim = 1.0d-4 |
|---|
| 27 | ch1 = dcadre( fch, sublim, az, aerr, rerr, error, ierr ) |
|---|
| 28 | ch = ax * sin(az) * ch1 |
|---|
| 29 | |
|---|
| 30 | return |
|---|
| 31 | end |
|---|
| 32 | |
|---|
| 33 | c ************************************************************ |
|---|
| 34 | c |
|---|
| 35 | c integrand for the chapman function |
|---|
| 36 | c |
|---|
| 37 | c ************************************************************ |
|---|
| 38 | |
|---|
| 39 | double precision function fch(gi) |
|---|
| 40 | |
|---|
| 41 | double precision gi |
|---|
| 42 | |
|---|
| 43 | common/integ/x,z |
|---|
| 44 | double precision x,z |
|---|
| 45 | |
|---|
| 46 | ! real x,z,a |
|---|
| 47 | |
|---|
| 48 | fch = (1./sin(gi))**2. * exp( x *(1.- sin(z)/sin(gi)) ) |
|---|
| 49 | |
|---|
| 50 | return |
|---|
| 51 | end |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | ****************************************************************************** |
|---|
| 57 | |
|---|
| 58 | double precision function dcadre (f,a,b,aerr,rerr,error,ier) |
|---|
| 59 | c specifications for arguments |
|---|
| 60 | integer ier |
|---|
| 61 | double precision f,a,b,aerr,rerr,error |
|---|
| 62 | c specifications for local variables |
|---|
| 63 | integer ibegs(30),maxts,maxtbl,mxstge,ibeg,ii,nnleft |
|---|
| 64 | integer i,n2,iii,istep2,iend,istep,l,lm1,it,istage,n |
|---|
| 65 | double precision t(10,10),r(10),ait(10),dif(10),rn(4),ts(2049) |
|---|
| 66 | double precision begin(30),finis(30),est(30) |
|---|
| 67 | double precision h2tol,aittol,length,jumptl,zero,p1,half,one |
|---|
| 68 | double precision two,four,fourp5,ten,hun,cadre,aitlow |
|---|
| 69 | double precision stepmn,stepnm,stage,curest,fnsize,hrerr |
|---|
| 70 | double precision prever,beg,fbeg,edn,fend,step,astep,tabs,hovn |
|---|
| 71 | double precision fn,sum,sumabs,vint,tabtlm,ergl,ergoal |
|---|
| 72 | double precision erra,errr,fextrp,errer,diff,sing,fextm1 |
|---|
| 73 | double precision h2next,singnx,slope,fbeg2,erret,h2tfex,fi |
|---|
| 74 | logical h2conv,aitken,right,reglar,reglsv(30) |
|---|
| 75 | data aitlow,h2tol,aittol,jumptl,maxts,maxtbl,mxstge |
|---|
| 76 | 1 /1.1d0,.15d0,.1d0,.01d0,2049,10,30/ |
|---|
| 77 | data rn(1),rn(2),rn(3),rn(4)/ |
|---|
| 78 | 1 .7142005d0,.3466282d0,.843751d0,.1263305d0/ |
|---|
| 79 | data zero,p1,half,one,two,four,fourp5,ten,hun |
|---|
| 80 | 1 /0.0d0,0.1d0,0.5d0,1.0d0,2.0d0,4.0d0, |
|---|
| 81 | 2 4.5d0,10.0d0,100.0d0/ |
|---|
| 82 | c first executable statement |
|---|
| 83 | ier = 0 |
|---|
| 84 | cadre = zero |
|---|
| 85 | error = zero |
|---|
| 86 | curest = zero |
|---|
| 87 | vint = zero |
|---|
| 88 | length = dabs(b-a) |
|---|
| 89 | |
|---|
| 90 | if (length .eq. zero) go to 215 |
|---|
| 91 | if (rerr .gt. p1 .or. rerr .lt. zero) go to 210 |
|---|
| 92 | hrerr = rerr+hun |
|---|
| 93 | if (aerr .eq. zero .and. hrerr .le. hun) go to 210 |
|---|
| 94 | errr = rerr |
|---|
| 95 | erra = dabs(aerr) |
|---|
| 96 | stepmn = length/(two**mxstge) |
|---|
| 97 | stepnm = dmax1(length,dabs(a),dabs(b))*ten |
|---|
| 98 | stage = half |
|---|
| 99 | istage = 1 |
|---|
| 100 | fnsize = zero |
|---|
| 101 | prever = zero |
|---|
| 102 | reglar = .false. |
|---|
| 103 | c the given interval of integration |
|---|
| 104 | c is the first interval considered. |
|---|
| 105 | beg = a |
|---|
| 106 | fbeg = f(beg)*half |
|---|
| 107 | ts(1) = fbeg |
|---|
| 108 | ibeg = 1 |
|---|
| 109 | edn = b |
|---|
| 110 | fend = f(edn)*half |
|---|
| 111 | ts(2) = fend |
|---|
| 112 | iend = 2 |
|---|
| 113 | 5 right = .false. |
|---|
| 114 | c investigation of a particular |
|---|
| 115 | c subinterval begins at this point. |
|---|
| 116 | 10 step = edn - beg |
|---|
| 117 | astep = dabs(step) |
|---|
| 118 | if (astep .lt. stepmn) go to 205 |
|---|
| 119 | hrerr = stepnm+astep |
|---|
| 120 | if (hrerr .eq. stepnm) go to 205 |
|---|
| 121 | t(1,1) = fbeg + fend |
|---|
| 122 | tabs = dabs(fbeg) + dabs(fend) |
|---|
| 123 | l = 1 |
|---|
| 124 | n = 1 |
|---|
| 125 | h2conv = .false. |
|---|
| 126 | aitken = .false. |
|---|
| 127 | 15 lm1 = l |
|---|
| 128 | l = l + 1 |
|---|
| 129 | c calculate the next trapezoid sum, |
|---|
| 130 | c t(l,1), which is based on *n2* + 1 |
|---|
| 131 | c equispaced points. here, |
|---|
| 132 | c n2 = n*2 = 2**(l-1). |
|---|
| 133 | n2 = n+n |
|---|
| 134 | fn = n2 |
|---|
| 135 | istep = (iend - ibeg)/n |
|---|
| 136 | if (istep .gt. 1) go to 25 |
|---|
| 137 | ii = iend |
|---|
| 138 | iend = iend + n |
|---|
| 139 | if (iend .gt. maxts) go to 200 |
|---|
| 140 | hovn = step/fn |
|---|
| 141 | iii = iend |
|---|
| 142 | fi = one |
|---|
| 143 | do 20 i=1,n2,2 |
|---|
| 144 | ts(iii) = ts(ii) |
|---|
| 145 | ts(iii-1) = f(edn - fi * hovn) |
|---|
| 146 | fi = fi+two |
|---|
| 147 | iii = iii-2 |
|---|
| 148 | ii = ii-1 |
|---|
| 149 | 20 continue |
|---|
| 150 | istep = 2 |
|---|
| 151 | 25 istep2 = ibeg + istep/2 |
|---|
| 152 | sum = zero |
|---|
| 153 | sumabs = zero |
|---|
| 154 | do 30 i=istep2,iend,istep |
|---|
| 155 | sum = sum + ts(i) |
|---|
| 156 | sumabs = sumabs + dabs(ts(i)) |
|---|
| 157 | 30 continue |
|---|
| 158 | t(l,1) = t(l-1,1)*half+sum/fn |
|---|
| 159 | tabs = tabs*half+sumabs/fn |
|---|
| 160 | n = n2 |
|---|
| 161 | c get preliminary value for *vint* |
|---|
| 162 | c from last trapezoid sum and update |
|---|
| 163 | c the error requirement *ergoal* |
|---|
| 164 | c for this subinterval. |
|---|
| 165 | it = 1 |
|---|
| 166 | vint = step*t(l,1) |
|---|
| 167 | tabtlm = tabs*ten |
|---|
| 168 | fnsize = dmax1(fnsize,dabs(t(l,1))) |
|---|
| 169 | ergl = astep*fnsize*ten |
|---|
| 170 | ergoal = stage*dmax1(erra,errr*dabs(curest+vint)) |
|---|
| 171 | c complete row l and column l of *t* |
|---|
| 172 | c array. |
|---|
| 173 | fextrp = one |
|---|
| 174 | do 35 i=1,lm1 |
|---|
| 175 | fextrp = fextrp*four |
|---|
| 176 | t(i,l) = t(l,i) - t(l-1,i) |
|---|
| 177 | t(l,i+1) = t(l,i) + t(i,l)/(fextrp-one) |
|---|
| 178 | 35 continue |
|---|
| 179 | errer = astep*dabs(t(1,l)) |
|---|
| 180 | c preliminary decision procedure |
|---|
| 181 | c if l = 2 and t(2,1) = t(1,1), |
|---|
| 182 | c go to 135 to follow up the |
|---|
| 183 | c impression that intergrand is |
|---|
| 184 | c straight line. |
|---|
| 185 | if (l .gt. 2) go to 40 |
|---|
| 186 | hrerr = tabs+p1*dabs(t(1,2)) |
|---|
| 187 | if (hrerr .eq. tabs) go to 135 |
|---|
| 188 | go to 15 |
|---|
| 189 | c caculate next ratios for |
|---|
| 190 | c columns 1,...,l-2 of t-table |
|---|
| 191 | c ratio is set to zero if difference |
|---|
| 192 | c in last two entries of column is |
|---|
| 193 | c about zero |
|---|
| 194 | 40 do 45 i=2,lm1 |
|---|
| 195 | diff = zero |
|---|
| 196 | hrerr = tabtlm+dabs(t(i-1,l)) |
|---|
| 197 | if (hrerr .ne. tabtlm) diff = t(i-1,lm1)/t(i-1,l) |
|---|
| 198 | t(i-1,lm1) = diff |
|---|
| 199 | 45 continue |
|---|
| 200 | if (dabs(four-t(1,lm1)) .le. h2tol) go to 60 |
|---|
| 201 | if (t(1,lm1) .eq. zero) go to 55 |
|---|
| 202 | if (dabs(two-dabs(t(1,lm1))) .lt. jumptl) go to 130 |
|---|
| 203 | if (l .eq. 3) go to 15 |
|---|
| 204 | h2conv = .false. |
|---|
| 205 | if (dabs((t(1,lm1)-t(1,l-2))/t(1,lm1)) .le. aittol) go to 75 |
|---|
| 206 | 50 if (reglar) go to 55 |
|---|
| 207 | if (l .eq. 4) go to 15 |
|---|
| 208 | hrerr = ergl+errer |
|---|
| 209 | 55 if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 175 |
|---|
| 210 | go to 145 |
|---|
| 211 | c cautious romberg extrapolation |
|---|
| 212 | 60 if (h2conv) go to 65 |
|---|
| 213 | aitken = .false. |
|---|
| 214 | h2conv = .true. |
|---|
| 215 | 65 fextrp = four |
|---|
| 216 | 70 it = it + 1 |
|---|
| 217 | vint = step*t(l,it) |
|---|
| 218 | errer = dabs(step/(fextrp-one)*t(it-1,l)) |
|---|
| 219 | if (errer .le. ergoal) go to 160 |
|---|
| 220 | hrerr = ergl+errer |
|---|
| 221 | if (hrerr .eq. ergl) go to 160 |
|---|
| 222 | if (it .eq. lm1) go to 125 |
|---|
| 223 | if (t(it,lm1) .eq. zero) go to 70 |
|---|
| 224 | if (t(it,lm1) .le. fextrp) go to 125 |
|---|
| 225 | if (dabs(t(it,lm1)/four-fextrp)/fextrp .lt. aittol) |
|---|
| 226 | 1 fextrp = fextrp*four |
|---|
| 227 | go to 70 |
|---|
| 228 | c integrand may have x**alpha type |
|---|
| 229 | c singularity |
|---|
| 230 | c resulting in a ratio of *sing* = |
|---|
| 231 | c 2**(alpha + 1) |
|---|
| 232 | 75 if (t(1,lm1) .lt. aitlow) go to 175 |
|---|
| 233 | if (aitken) go to 80 |
|---|
| 234 | h2conv = .false. |
|---|
| 235 | aitken = .true. |
|---|
| 236 | 80 fextrp = t(l-2,lm1) |
|---|
| 237 | if (fextrp .gt. fourp5) go to 65 |
|---|
| 238 | if (fextrp .lt. aitlow) go to 175 |
|---|
| 239 | if (dabs(fextrp-t(l-3,lm1))/t(1,lm1) .gt. h2tol) go to 175 |
|---|
| 240 | sing = fextrp |
|---|
| 241 | fextm1 = one/(fextrp - one) |
|---|
| 242 | ait(1) = zero |
|---|
| 243 | do 85 i=2,l |
|---|
| 244 | ait(i) = t(i,1) + (t(i,1)-t(i-1,1))*fextm1 |
|---|
| 245 | r(i) = t(1,i-1) |
|---|
| 246 | dif(i) = ait(i) - ait(i-1) |
|---|
| 247 | 85 continue |
|---|
| 248 | it = 2 |
|---|
| 249 | 90 vint = step*ait(l) |
|---|
| 250 | errer = errer*fextm1 |
|---|
| 251 | hrerr = ergl+errer |
|---|
| 252 | if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 95 |
|---|
| 253 | ier = max0(ier,65) |
|---|
| 254 | go to 160 |
|---|
| 255 | 95 it = it + 1 |
|---|
| 256 | if (it .eq. lm1) go to 125 |
|---|
| 257 | if (it .gt. 3) go to 100 |
|---|
| 258 | h2next = four |
|---|
| 259 | singnx = sing+sing |
|---|
| 260 | 100 if (h2next .lt. singnx) go to 105 |
|---|
| 261 | fextrp = singnx |
|---|
| 262 | singnx = singnx+singnx |
|---|
| 263 | go to 110 |
|---|
| 264 | 105 fextrp = h2next |
|---|
| 265 | h2next = four*h2next |
|---|
| 266 | 110 do 115 i=it,lm1 |
|---|
| 267 | r(i+1) = zero |
|---|
| 268 | hrerr = tabtlm+dabs(dif(i+1)) |
|---|
| 269 | if (hrerr .ne. tabtlm) r(i+1) = dif(i)/dif(i+1) |
|---|
| 270 | 115 continue |
|---|
| 271 | h2tfex = -h2tol*fextrp |
|---|
| 272 | if (r(l) - fextrp .lt. h2tfex) go to 125 |
|---|
| 273 | if (r(l-1)-fextrp .lt. h2tfex) go to 125 |
|---|
| 274 | errer = astep*dabs(dif(l)) |
|---|
| 275 | fextm1 = one/(fextrp - one) |
|---|
| 276 | do 120 i=it,l |
|---|
| 277 | ait(i) = ait(i) + dif(i)*fextm1 |
|---|
| 278 | dif(i) = ait(i) - ait(i-1) |
|---|
| 279 | 120 continue |
|---|
| 280 | go to 90 |
|---|
| 281 | c current trapezoid sum and resulting |
|---|
| 282 | c extrapolated values did not give |
|---|
| 283 | c a small enough *errer*. |
|---|
| 284 | c note -- having prever .lt. errer |
|---|
| 285 | c is an almost certain sign of |
|---|
| 286 | c beginning trouble with in the func- |
|---|
| 287 | c tion values. hence, a watch for, |
|---|
| 288 | c and control of, noise should |
|---|
| 289 | c begin here. |
|---|
| 290 | 125 fextrp = dmax1(prever/errer,aitlow) |
|---|
| 291 | prever = errer |
|---|
| 292 | if (l .lt. 5) go to 15 |
|---|
| 293 | if (l-it .gt. 2 .and. istage .lt. mxstge) go to 170 |
|---|
| 294 | erret = errer/(fextrp**(maxtbl-l)) |
|---|
| 295 | hrerr = ergl+erret |
|---|
| 296 | if (erret .gt. ergoal .and. hrerr .ne. ergl) go to 170 |
|---|
| 297 | go to 15 |
|---|
| 298 | c integrand has jump (see notes) |
|---|
| 299 | 130 hrerr = ergl+errer |
|---|
| 300 | if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 170 |
|---|
| 301 | c note that 2*fn = 2**l |
|---|
| 302 | diff = dabs(t(1,l))*(fn+fn) |
|---|
| 303 | go to 160 |
|---|
| 304 | c integrand is straight line |
|---|
| 305 | c test this assumption by comparing |
|---|
| 306 | c the value of the integrand at |
|---|
| 307 | c four *randomly chosen* points with |
|---|
| 308 | c the value of the straight line |
|---|
| 309 | c interpolating the integrand at the |
|---|
| 310 | c two end points of the sub-interval. |
|---|
| 311 | c if test is passed, accept *vint* |
|---|
| 312 | 135 slope = (fend-fbeg)*two |
|---|
| 313 | fbeg2 = fbeg+fbeg |
|---|
| 314 | do 140 i=1,4 |
|---|
| 315 | diff = dabs(f(beg+rn(i)*step) - fbeg2-rn(i)*slope) |
|---|
| 316 | hrerr = tabtlm+diff |
|---|
| 317 | if(hrerr .ne. tabtlm) go to 155 |
|---|
| 318 | 140 continue |
|---|
| 319 | go to 160 |
|---|
| 320 | c noise may be dominant feature |
|---|
| 321 | c estimate noise level by comparing |
|---|
| 322 | c the value of the integrand at |
|---|
| 323 | c four *randomly chosen* points with |
|---|
| 324 | c the value of the straight line |
|---|
| 325 | c interpolating the integrand at the |
|---|
| 326 | c two endpoints. if small enough, |
|---|
| 327 | c accept *vint* |
|---|
| 328 | 145 slope = (fend-fbeg)*two |
|---|
| 329 | fbeg2 = fbeg+fbeg |
|---|
| 330 | i = 1 |
|---|
| 331 | 150 diff = dabs(f(beg+rn(i)*step) - fbeg2-rn(i)*slope) |
|---|
| 332 | 155 errer = dmax1(errer,astep*diff) |
|---|
| 333 | hrerr = ergl+errer |
|---|
| 334 | if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 175 |
|---|
| 335 | i = i+1 |
|---|
| 336 | if (i .le. 4) go to 150 |
|---|
| 337 | ier = 66 |
|---|
| 338 | c intergration over current sub- |
|---|
| 339 | c interval successful |
|---|
| 340 | c add *vint* to *dcadre* and *errer* |
|---|
| 341 | c to *error*, then set up next sub- |
|---|
| 342 | c interval, if any. |
|---|
| 343 | 160 cadre = cadre + vint |
|---|
| 344 | error = error + errer |
|---|
| 345 | if (right) go to 165 |
|---|
| 346 | istage = istage - 1 |
|---|
| 347 | if (istage .eq. 0) go to 220 |
|---|
| 348 | reglar = reglsv(istage) |
|---|
| 349 | beg = begin(istage) |
|---|
| 350 | edn = finis(istage) |
|---|
| 351 | curest = curest - est(istage+1) + vint |
|---|
| 352 | iend = ibeg - 1 |
|---|
| 353 | fend = ts(iend) |
|---|
| 354 | ibeg = ibegs(istage) |
|---|
| 355 | go to 180 |
|---|
| 356 | 165 curest = curest + vint |
|---|
| 357 | stage = stage+stage |
|---|
| 358 | iend = ibeg |
|---|
| 359 | ibeg = ibegs(istage) |
|---|
| 360 | edn = beg |
|---|
| 361 | beg = begin(istage) |
|---|
| 362 | fend = fbeg |
|---|
| 363 | fbeg = ts(ibeg) |
|---|
| 364 | go to 5 |
|---|
| 365 | c integration over current subinterval |
|---|
| 366 | c is unsuccessful. mark subinterval |
|---|
| 367 | c for further subdivision. set up |
|---|
| 368 | c next subinterval. |
|---|
| 369 | 170 reglar = .true. |
|---|
| 370 | 175 if (istage .eq. mxstge) go to 205 |
|---|
| 371 | if (right) go to 185 |
|---|
| 372 | reglsv(istage+1) = reglar |
|---|
| 373 | begin(istage) = beg |
|---|
| 374 | ibegs(istage) = ibeg |
|---|
| 375 | stage = stage*half |
|---|
| 376 | 180 right = .true. |
|---|
| 377 | beg = (beg+edn)*half |
|---|
| 378 | ibeg = (ibeg+iend)/2 |
|---|
| 379 | ts(ibeg) = ts(ibeg)*half |
|---|
| 380 | fbeg = ts(ibeg) |
|---|
| 381 | go to 10 |
|---|
| 382 | 185 nnleft = ibeg - ibegs(istage) |
|---|
| 383 | if (iend+nnleft .ge. maxts) go to 200 |
|---|
| 384 | iii = ibegs(istage) |
|---|
| 385 | ii = iend |
|---|
| 386 | do 190 i=iii,ibeg |
|---|
| 387 | ii = ii + 1 |
|---|
| 388 | ts(ii) = ts(i) |
|---|
| 389 | 190 continue |
|---|
| 390 | do 195 i=ibeg,ii |
|---|
| 391 | ts(iii) = ts(i) |
|---|
| 392 | iii = iii + 1 |
|---|
| 393 | 195 continue |
|---|
| 394 | iend = iend + 1 |
|---|
| 395 | ibeg = iend - nnleft |
|---|
| 396 | fend = fbeg |
|---|
| 397 | fbeg = ts(ibeg) |
|---|
| 398 | finis(istage) = edn |
|---|
| 399 | edn = beg |
|---|
| 400 | beg = begin(istage) |
|---|
| 401 | begin(istage) = edn |
|---|
| 402 | reglsv(istage) = reglar |
|---|
| 403 | istage = istage + 1 |
|---|
| 404 | reglar = reglsv(istage) |
|---|
| 405 | est(istage) = vint |
|---|
| 406 | curest = curest + est(istage) |
|---|
| 407 | go to 5 |
|---|
| 408 | c failure to handle given integra- |
|---|
| 409 | c tion problem |
|---|
| 410 | 200 ier = 131 |
|---|
| 411 | go to 215 |
|---|
| 412 | 205 ier = 132 |
|---|
| 413 | go to 215 |
|---|
| 414 | 210 ier = 133 |
|---|
| 415 | 215 cadre = curest + vint |
|---|
| 416 | 220 dcadre = cadre |
|---|
| 417 | 9000 continue |
|---|
| 418 | if (ier .ne. 0) call uertst (ier,6hdcadre) |
|---|
| 419 | 9005 return |
|---|
| 420 | end |
|---|
| 421 | |
|---|
| 422 | |
|---|
| 423 | ****************************************************************************** |
|---|
| 424 | |
|---|
| 425 | subroutine uertst (ier,name) |
|---|
| 426 | c specifications for arguments |
|---|
| 427 | integer ier |
|---|
| 428 | integer name(2) |
|---|
| 429 | c specifications for local variables |
|---|
| 430 | integer i,ieq,ieqdf,iounit,level,levold,nameq(6), |
|---|
| 431 | * namset(6),namupk(6),nin,nmtb |
|---|
| 432 | data namset/1hu,1he,1hr,1hs,1he,1ht/ |
|---|
| 433 | data nameq/6*1h / |
|---|
| 434 | data level/4/,ieqdf/0/,ieq/1h=/ |
|---|
| 435 | c unpack name into namupk |
|---|
| 436 | c first executable statement |
|---|
| 437 | call uspkd (name,6,namupk,nmtb) |
|---|
| 438 | c get output unit number |
|---|
| 439 | call ugetio(1,nin,iounit) |
|---|
| 440 | c check ier |
|---|
| 441 | if (ier.gt.999) go to 25 |
|---|
| 442 | if (ier.lt.-32) go to 55 |
|---|
| 443 | if (ier.le.128) go to 5 |
|---|
| 444 | if (level.lt.1) go to 30 |
|---|
| 445 | c print terminal message |
|---|
| 446 | if (ieqdf.eq.1) write(iounit,35) ier,nameq,ieq,namupk |
|---|
| 447 | if (ieqdf.eq.0) write(iounit,35) ier,namupk |
|---|
| 448 | go to 30 |
|---|
| 449 | 5 if (ier.le.64) go to 10 |
|---|
| 450 | if (level.lt.2) go to 30 |
|---|
| 451 | c print warning with fix message |
|---|
| 452 | c if (ieqdf.eq.1) write(iounit,40) ier,nameq,ieq,namupk |
|---|
| 453 | c if (ieqdf.eq.0) write(iounit,40) ier,namupk |
|---|
| 454 | if (ieqdf.eq.1) continue |
|---|
| 455 | if (ieqdf.eq.0) continue |
|---|
| 456 | go to 30 |
|---|
| 457 | 10 if (ier.le.32) go to 15 |
|---|
| 458 | c print warning message |
|---|
| 459 | if (level.lt.3) go to 30 |
|---|
| 460 | if (ieqdf.eq.1) write(iounit,45) ier,nameq,ieq,namupk |
|---|
| 461 | if (ieqdf.eq.0) write(iounit,45) ier,namupk |
|---|
| 462 | go to 30 |
|---|
| 463 | 15 continue |
|---|
| 464 | c check for uerset call |
|---|
| 465 | do 20 i=1,6 |
|---|
| 466 | if (namupk(i).ne.namset(i)) go to 25 |
|---|
| 467 | 20 continue |
|---|
| 468 | levold = level |
|---|
| 469 | level = ier |
|---|
| 470 | ier = levold |
|---|
| 471 | if (level.lt.0) level = 4 |
|---|
| 472 | if (level.gt.4) level = 4 |
|---|
| 473 | go to 30 |
|---|
| 474 | 25 continue |
|---|
| 475 | if (level.lt.4) go to 30 |
|---|
| 476 | c print non-defined message |
|---|
| 477 | if (ieqdf.eq.1) write(iounit,50) ier,nameq,ieq,namupk |
|---|
| 478 | if (ieqdf.eq.0) write(iounit,50) ier,namupk |
|---|
| 479 | 30 ieqdf = 0 |
|---|
| 480 | return |
|---|
| 481 | 35 format(19h *** terminal error,10x,7h(ier = ,i3, |
|---|
| 482 | 1 20h) from imsl routine ,6a1,a1,6a1) |
|---|
| 483 | 40 format(27h *** warning with fix error,2x,7h(ier = ,i3, |
|---|
| 484 | 1 20h) from imsl routine ,6a1,a1,6a1) |
|---|
| 485 | 45 format(18h *** warning error,11x,7h(ier = ,i3, |
|---|
| 486 | 1 20h) from imsl routine ,6a1,a1,6a1) |
|---|
| 487 | 50 format(20h *** undefined error,9x,7h(ier = ,i5, |
|---|
| 488 | 1 20h) from imsl routine ,6a1,a1,6a1) |
|---|
| 489 | c |
|---|
| 490 | c save p for p = r case |
|---|
| 491 | c p is the page namupk |
|---|
| 492 | c r is the routine namupk |
|---|
| 493 | 55 ieqdf = 1 |
|---|
| 494 | do 60 i=1,6 |
|---|
| 495 | 60 nameq(i) = namupk(i) |
|---|
| 496 | 65 return |
|---|
| 497 | end |
|---|
| 498 | |
|---|
| 499 | |
|---|
| 500 | ************************************************************************ |
|---|
| 501 | |
|---|
| 502 | subroutine uspkd (packed,nchars,unpakd,nchmtb) |
|---|
| 503 | c specifications for arguments |
|---|
| 504 | integer nc,nchars,nchmtb |
|---|
| 505 | c |
|---|
| 506 | c integer unpakd(1),iblank |
|---|
| 507 | integer unpakd(nchars),iblank |
|---|
| 508 | character packed(1) ! Modificado el 29-Ago-2001 porque |
|---|
| 509 | ! integer packed(1) esto deberia ser CHARACTER, no?? |
|---|
| 510 | data iblank /1h / |
|---|
| 511 | c initialize nchmtb |
|---|
| 512 | nchmtb = 0 |
|---|
| 513 | c return if nchars is le zero |
|---|
| 514 | if(nchars.le.0) return |
|---|
| 515 | c set nc=number of chars to be decoded |
|---|
| 516 | nc = min0 (129,nchars) |
|---|
| 517 | ! decode (nc,150,packed) (unpakd(i),i=1,nc) |
|---|
| 518 | read (packed,150) (unpakd(i),i=1,nc) |
|---|
| 519 | 150 format (129a1) |
|---|
| 520 | c check unpakd array and set nchmtb |
|---|
| 521 | c based on trailing blanks found |
|---|
| 522 | do 200 n = 1,nc |
|---|
| 523 | nn = nc - n + 1 |
|---|
| 524 | if(unpakd(nn) .ne. 536870912) go to 210 |
|---|
| 525 | 200 continue |
|---|
| 526 | 210 nchmtb = nn |
|---|
| 527 | return |
|---|
| 528 | end |
|---|
| 529 | |
|---|
| 530 | |
|---|
| 531 | **************************************************************************** |
|---|
| 532 | |
|---|
| 533 | subroutine ugetio(iopt,nin,nout) |
|---|
| 534 | c specifications for arguments |
|---|
| 535 | integer iopt,nin,nout |
|---|
| 536 | c specifications for local variables |
|---|
| 537 | integer nind,noutd |
|---|
| 538 | data nind/5/,noutd/6/ |
|---|
| 539 | c first executable statement |
|---|
| 540 | if (iopt.eq.3) go to 10 |
|---|
| 541 | if (iopt.eq.2) go to 5 |
|---|
| 542 | if (iopt.ne.1) go to 9005 |
|---|
| 543 | nin = nind |
|---|
| 544 | nout = noutd |
|---|
| 545 | go to 9005 |
|---|
| 546 | 5 nind = nin |
|---|
| 547 | go to 9005 |
|---|
| 548 | 10 noutd = nout |
|---|
| 549 | 9005 return |
|---|
| 550 | end |
|---|
