1C Output from Public domain Ratfor, version 1.01
2      subroutine qpsedg8xf(tgiyxdw1, dufozmt7, wy1vqfzu)
3      implicit logical (a-z)
4      integer wy1vqfzu, tgiyxdw1(*), dufozmt7(*)
5      integer urohxe6t, bpvaqm5z, ayfnwr1v
6      ayfnwr1v = 1
7      urohxe6t = wy1vqfzu
823000 if(.not.(urohxe6t .ge. 1))goto 23002
9      do23003 bpvaqm5z=1,urohxe6t
10      tgiyxdw1(ayfnwr1v) = bpvaqm5z
11      ayfnwr1v = ayfnwr1v+1
1223003 continue
1323004 continue
1423001 urohxe6t=urohxe6t-1
15      goto 23000
1623002 continue
17      ayfnwr1v = 1
18      do23005 urohxe6t=1,wy1vqfzu
19      do23007 bpvaqm5z=urohxe6t,wy1vqfzu
20      dufozmt7(ayfnwr1v) = bpvaqm5z
21      ayfnwr1v = ayfnwr1v+1
2223007 continue
2323008 continue
2423005 continue
2523006 continue
26      return
27      end
28      integer function viamf(cz8qdfyj, rvy1fpli, wy1vqfzu, tgiyxdw1, duf
29     *ozmt7)
30      integer cz8qdfyj, rvy1fpli, wy1vqfzu, tgiyxdw1(*), dufozmt7(*)
31      integer urohxe6t, imk5wjxg
32      imk5wjxg = wy1vqfzu*(wy1vqfzu+1)/2
33      do23009 urohxe6t=1,imk5wjxg
34      if((tgiyxdw1(urohxe6t).eq.cz8qdfyj .and. dufozmt7(urohxe6t).eq.rvy
35     *1fpli) .or. (tgiyxdw1(urohxe6t).eq.rvy1fpli .and. dufozmt7(urohxe6
36     *t).eq.cz8qdfyj))then
37      viamf = urohxe6t
38      return
39      endif
4023009 continue
4123010 continue
42      viamf = 0
43      return
44      end
45      subroutine vm2af(mat, a, dimm, tgiyxdw1, dufozmt7, kuzxj1lo, wy1vq
46     *fzu, rb1onzwu)
47      implicit logical (a-z)
48      integer dimm, tgiyxdw1(dimm), dufozmt7(dimm), kuzxj1lo, wy1vqfzu,
49     *rb1onzwu
50      double precision mat(dimm,kuzxj1lo), a(wy1vqfzu,wy1vqfzu,kuzxj1lo)
51      integer ayfnwr1v, yq6lorbx, gp1jxzuh, imk5wjxg
52      imk5wjxg = wy1vqfzu * (wy1vqfzu + 1) / 2
53      if(rb1onzwu .eq. 1 .or. dimm .ne. imk5wjxg)then
54      ayfnwr1v = 1
5523015 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23017
56      yq6lorbx = 1
5723018 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23020
58      gp1jxzuh = 1
5923021 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23023
60      a(gp1jxzuh,yq6lorbx,ayfnwr1v) = 0.0d0
6123022 gp1jxzuh=gp1jxzuh+1
62      goto 23021
6323023 continue
6423019 yq6lorbx=yq6lorbx+1
65      goto 23018
6623020 continue
6723016 ayfnwr1v=ayfnwr1v+1
68      goto 23015
6923017 continue
70      endif
71      do23024 ayfnwr1v=1,kuzxj1lo
72      do23026 yq6lorbx=1,dimm
73      a(tgiyxdw1(yq6lorbx),dufozmt7(yq6lorbx),ayfnwr1v) = mat(yq6lorbx,a
74     *yfnwr1v)
75      if(rb1onzwu .eq. 0)then
76      a(dufozmt7(yq6lorbx),tgiyxdw1(yq6lorbx),ayfnwr1v) = mat(yq6lorbx,a
77     *yfnwr1v)
78      endif
7923026 continue
8023027 continue
8123024 continue
8223025 continue
83      return
84      end
85      subroutine mux22f(wpuarq2m, tlgduey8, lfu2qhid, dimu, tgiyxdw1, du
86     *fozmt7, kuzxj1lo, wy1vqfzu, wk1200)
87      implicit logical (a-z)
88      integer dimu, tgiyxdw1(*), dufozmt7(*), kuzxj1lo, wy1vqfzu
89      double precision wpuarq2m(dimu,kuzxj1lo), tlgduey8(kuzxj1lo,wy1vqf
90     *zu), lfu2qhid(wy1vqfzu,kuzxj1lo), wk1200(wy1vqfzu,wy1vqfzu)
91      double precision q6zdcwxk
92      integer ayfnwr1v, yq6lorbx, bpvaqm5z, one, rb1onzwu
93      one = 1
94      rb1onzwu = 1
95      ayfnwr1v = 1
9623030 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23032
97      call vm2af(wpuarq2m(1,ayfnwr1v), wk1200, dimu, tgiyxdw1, dufozmt7,
98     * one, wy1vqfzu, rb1onzwu)
99      yq6lorbx = 1
10023033 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23035
101      q6zdcwxk = 0.0d0
102      bpvaqm5z = yq6lorbx
10323036 if(.not.(bpvaqm5z .le. wy1vqfzu))goto 23038
104      q6zdcwxk = q6zdcwxk + wk1200(yq6lorbx,bpvaqm5z) * tlgduey8(ayfnwr1
105     *v,bpvaqm5z)
10623037 bpvaqm5z=bpvaqm5z+1
107      goto 23036
10823038 continue
109      lfu2qhid(yq6lorbx,ayfnwr1v) = q6zdcwxk
11023034 yq6lorbx=yq6lorbx+1
111      goto 23033
11223035 continue
11323031 ayfnwr1v=ayfnwr1v+1
114      goto 23030
11523032 continue
116      return
117      end
118      subroutine vbksf(wpuarq2m, bvecto, wy1vqfzu, kuzxj1lo, wk1200, tgi
119     *yxdw1, dufozmt7, dimu)
120      implicit logical (a-z)
121      integer wy1vqfzu, kuzxj1lo, tgiyxdw1(*), dufozmt7(*), dimu
122      double precision wpuarq2m(dimu,kuzxj1lo), bvecto(wy1vqfzu,kuzxj1lo
123     *), wk1200(wy1vqfzu,wy1vqfzu)
124      double precision q6zdcwxk
125      integer ayfnwr1v, yq6lorbx, gp1jxzuh, rb1onzwu, one
126      rb1onzwu = 1
127      one = 1
128      ayfnwr1v = 1
12923039 if(.not.(ayfnwr1v .le. kuzxj1lo))goto 23041
130      call vm2af(wpuarq2m(1,ayfnwr1v), wk1200, dimu, tgiyxdw1, dufozmt7,
131     * one, wy1vqfzu, rb1onzwu)
132      yq6lorbx = wy1vqfzu
13323042 if(.not.(yq6lorbx .ge. 1))goto 23044
134      q6zdcwxk = bvecto(yq6lorbx,ayfnwr1v)
135      gp1jxzuh = yq6lorbx+1
13623045 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23047
137      q6zdcwxk = q6zdcwxk - wk1200(yq6lorbx,gp1jxzuh) * bvecto(gp1jxzuh,
138     *ayfnwr1v)
13923046 gp1jxzuh=gp1jxzuh+1
140      goto 23045
14123047 continue
142      bvecto(yq6lorbx,ayfnwr1v) = q6zdcwxk / wk1200(yq6lorbx,yq6lorbx)
14323043 yq6lorbx=yq6lorbx-1
144      goto 23042
14523044 continue
14623040 ayfnwr1v=ayfnwr1v+1
147      goto 23039
14823041 continue
149      return
150      end
151      subroutine vcholf(wmat, bvecto, wy1vqfzu, dvhw1ulq, isolve)
152      implicit logical (a-z)
153      integer isolve
154      integer wy1vqfzu, dvhw1ulq
155      double precision wmat(wy1vqfzu,wy1vqfzu), bvecto(wy1vqfzu)
156      double precision q6zdcwxk, dsqrt
157      integer ayfnwr1v, yq6lorbx, gp1jxzuh
158      dvhw1ulq=1
159      do23048 ayfnwr1v=1,wy1vqfzu
160      q6zdcwxk = 0d0
161      do23050 gp1jxzuh=1,ayfnwr1v-1
162      q6zdcwxk = q6zdcwxk + wmat(gp1jxzuh,ayfnwr1v) * wmat(gp1jxzuh,ayfn
163     *wr1v)
16423050 continue
16523051 continue
166      wmat(ayfnwr1v,ayfnwr1v) = wmat(ayfnwr1v,ayfnwr1v) - q6zdcwxk
167      if(wmat(ayfnwr1v,ayfnwr1v) .le. 0d0)then
168      dvhw1ulq = 0
169      return
170      endif
171      wmat(ayfnwr1v,ayfnwr1v) = dsqrt(wmat(ayfnwr1v,ayfnwr1v))
172      do23054 yq6lorbx=ayfnwr1v+1,wy1vqfzu
173      q6zdcwxk = 0d0
174      do23056 gp1jxzuh=1,ayfnwr1v-1
175      q6zdcwxk = q6zdcwxk + wmat(gp1jxzuh,ayfnwr1v) * wmat(gp1jxzuh,yq6l
176     *orbx)
17723056 continue
17823057 continue
179      wmat(ayfnwr1v,yq6lorbx) = (wmat(ayfnwr1v,yq6lorbx) - q6zdcwxk) / w
180     *mat(ayfnwr1v,ayfnwr1v)
18123054 continue
18223055 continue
18323048 continue
18423049 continue
185      if(isolve .eq. 0)then
186      do23060 ayfnwr1v=2,wy1vqfzu
187      do23062 yq6lorbx=1,ayfnwr1v-1
188      wmat(ayfnwr1v,yq6lorbx) = 0.0d0
18923062 continue
19023063 continue
191      return
19223060 continue
19323061 continue
194      endif
195      do23064 yq6lorbx=1,wy1vqfzu
196      q6zdcwxk = bvecto(yq6lorbx)
197      do23066 gp1jxzuh=1,yq6lorbx-1
198      q6zdcwxk = q6zdcwxk - wmat(gp1jxzuh,yq6lorbx) * bvecto(gp1jxzuh)
19923066 continue
20023067 continue
201      bvecto(yq6lorbx) = q6zdcwxk / wmat(yq6lorbx,yq6lorbx)
20223064 continue
20323065 continue
204      yq6lorbx = wy1vqfzu
20523068 if(.not.(yq6lorbx .ge. 1))goto 23070
206      q6zdcwxk = bvecto(yq6lorbx)
207      gp1jxzuh = yq6lorbx+1
20823071 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23073
209      q6zdcwxk = q6zdcwxk - wmat(yq6lorbx,gp1jxzuh) * bvecto(gp1jxzuh)
21023072 gp1jxzuh=gp1jxzuh+1
211      goto 23071
21223073 continue
213      bvecto(yq6lorbx) = q6zdcwxk / wmat(yq6lorbx,yq6lorbx)
21423069 yq6lorbx=yq6lorbx-1
215      goto 23068
21623070 continue
217      return
218      end
219      subroutine mux17f(wpuarq2m, he7mqnvy, wy1vqfzu, xjc4ywlh, kuzxj1lo
220     *, wk1200, wk3400, tgiyxdw1, dufozmt7, dimu, rutyk8mg)
221      implicit logical (a-z)
222      integer dimu, wy1vqfzu, xjc4ywlh, kuzxj1lo, tgiyxdw1(*), dufozmt7(
223     **), rutyk8mg
224      double precision wpuarq2m(dimu,kuzxj1lo), he7mqnvy(rutyk8mg,xjc4yw
225     *lh), wk1200(wy1vqfzu,wy1vqfzu), wk3400(wy1vqfzu,xjc4ywlh)
226      double precision q6zdcwxk
227      integer ayfnwr1v, yq6lorbx, gp1jxzuh, bpvaqm5z
228      do23074 yq6lorbx=1,wy1vqfzu
229      do23076 ayfnwr1v=1,wy1vqfzu
230      wk1200(ayfnwr1v,yq6lorbx) = 0.0d0
23123076 continue
23223077 continue
23323074 continue
23423075 continue
235      do23078 ayfnwr1v=1,kuzxj1lo
236      do23080 bpvaqm5z=1,dimu
237      wk1200(tgiyxdw1(bpvaqm5z), dufozmt7(bpvaqm5z)) = wpuarq2m(bpvaqm5z
238     *,ayfnwr1v)
23923080 continue
24023081 continue
241      do23082 gp1jxzuh=1,xjc4ywlh
242      do23084 yq6lorbx=1,wy1vqfzu
243      wk3400(yq6lorbx,gp1jxzuh) = he7mqnvy((ayfnwr1v-1)*wy1vqfzu+yq6lorb
244     *x,gp1jxzuh)
24523084 continue
24623085 continue
24723082 continue
24823083 continue
249      do23086 gp1jxzuh=1,xjc4ywlh
250      do23088 yq6lorbx=1,wy1vqfzu
251      q6zdcwxk = 0d0
252      do23090 bpvaqm5z=yq6lorbx,wy1vqfzu
253      q6zdcwxk = q6zdcwxk + wk1200(yq6lorbx,bpvaqm5z) * wk3400(bpvaqm5z,
254     *gp1jxzuh)
25523090 continue
25623091 continue
257      he7mqnvy((ayfnwr1v-1)*wy1vqfzu+yq6lorbx,gp1jxzuh) = q6zdcwxk
25823088 continue
25923089 continue
26023086 continue
26123087 continue
26223078 continue
26323079 continue
264      return
265      end
266      subroutine vrinvf9(wpuarq2m, ldr, wy1vqfzu, dvhw1ulq, ks3wejcv, wo
267     *rk)
268      implicit logical (a-z)
269      integer ldr, wy1vqfzu, dvhw1ulq
270      double precision wpuarq2m(ldr,wy1vqfzu), ks3wejcv(wy1vqfzu,wy1vqfz
271     *u), work(wy1vqfzu,wy1vqfzu)
272      double precision q6zdcwxk
273      integer yq6lorbx, gp1jxzuh, col, uaoynef0
274      dvhw1ulq = 1
275      yq6lorbx = 1
27623092 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23094
277      col = 1
27823095 if(.not.(col .le. wy1vqfzu))goto 23097
279      work(yq6lorbx,col) = 0.0d0
28023096 col=col+1
281      goto 23095
28223097 continue
28323093 yq6lorbx=yq6lorbx+1
284      goto 23092
28523094 continue
286      col = 1
28723098 if(.not.(col .le. wy1vqfzu))goto 23100
288      yq6lorbx = col
28923101 if(.not.(yq6lorbx .ge. 1))goto 23103
290      if(yq6lorbx .eq. col)then
291      q6zdcwxk = 1.0d0
292      else
293      q6zdcwxk = 0.0d0
294      endif
295      gp1jxzuh = yq6lorbx+1
29623106 if(.not.(gp1jxzuh .le. col))goto 23108
297      q6zdcwxk = q6zdcwxk - wpuarq2m(yq6lorbx,gp1jxzuh) * work(gp1jxzuh,
298     *col)
29923107 gp1jxzuh=gp1jxzuh+1
300      goto 23106
30123108 continue
302      if(wpuarq2m(yq6lorbx,yq6lorbx) .eq. 0.0d0)then
303      dvhw1ulq = 0
304      else
305      work(yq6lorbx,col) = q6zdcwxk / wpuarq2m(yq6lorbx,yq6lorbx)
306      endif
30723102 yq6lorbx=yq6lorbx-1
308      goto 23101
30923103 continue
31023099 col=col+1
311      goto 23098
31223100 continue
313      yq6lorbx = 1
31423111 if(.not.(yq6lorbx .le. wy1vqfzu))goto 23113
315      col = yq6lorbx
31623114 if(.not.(col .le. wy1vqfzu))goto 23116
317      if(yq6lorbx .lt. col)then
318      uaoynef0 = col
319      else
320      uaoynef0 = yq6lorbx
321      endif
322      q6zdcwxk = 0.0d0
323      gp1jxzuh = uaoynef0
32423119 if(.not.(gp1jxzuh .le. wy1vqfzu))goto 23121
325      q6zdcwxk = q6zdcwxk + work(yq6lorbx,gp1jxzuh) * work(col,gp1jxzuh)
32623120 gp1jxzuh=gp1jxzuh+1
327      goto 23119
32823121 continue
329      ks3wejcv(yq6lorbx,col) = q6zdcwxk
330      ks3wejcv(col,yq6lorbx) = q6zdcwxk
33123115 col=col+1
332      goto 23114
33323116 continue
33423112 yq6lorbx=yq6lorbx+1
335      goto 23111
33623113 continue
337      return
338      end
339      subroutine tldz5ion(xx, lfu2qhid)
340      implicit logical (a-z)
341      double precision xx, lfu2qhid
342      double precision x, y, hofjnx2e, q6zdcwxk, xd4mybgj(6)
343      integer yq6lorbx
344      xd4mybgj(1)= 76.18009172947146d0
345      xd4mybgj(2)= -86.50532032941677d0
346      xd4mybgj(3)= 24.01409824083091d0
347      xd4mybgj(4)= -1.231739572450155d0
348      xd4mybgj(5)= 0.1208650973866179d-2
349      xd4mybgj(6)= -0.5395239384953d-5
350      x = xx
351      y = xx
352      hofjnx2e = x+5.50d0
353      hofjnx2e = hofjnx2e - (x+0.50d0) * dlog(hofjnx2e)
354      q6zdcwxk=1.000000000190015d0
355      yq6lorbx=1
35623122 if(.not.(yq6lorbx .le. 6))goto 23124
357      y = y + 1.0d0
358      q6zdcwxk = q6zdcwxk + xd4mybgj(yq6lorbx)/y
35923123 yq6lorbx=yq6lorbx+1
360      goto 23122
36123124 continue
362      lfu2qhid = -hofjnx2e + dlog(2.5066282746310005d0 * q6zdcwxk / x)
363      return
364      end
365      subroutine enbin9(bzmd6ftv, hdqsx7bk, nm0eljqk, n2kersmx, n, dvhw1
366     *ulq, zy1mchbf, ux3nadiw, rsynp1go, sguwj9ty)
367      implicit logical (a-z)
368      integer n, dvhw1ulq, zy1mchbf, sguwj9ty
369      double precision bzmd6ftv(n, zy1mchbf), hdqsx7bk(n, zy1mchbf), nm0
370     *eljqk(n, zy1mchbf), n2kersmx, ux3nadiw, rsynp1go
371      integer ayfnwr1v, kij0gwer
372      double precision oxjgzv0e, btiehdm2, ydb, vjz5sxty, esql7umk, pvcj
373     *l2na, mwuvskg1, ft3ijqmy, hmayv1xt, q6zdcwxk, plo6hkdr
374      real csi9ydge
375      if(n2kersmx .le. 0.80d0 .or. n2kersmx .ge. 1.0d0)then
376      dvhw1ulq = 0
377      return
378      endif
379      btiehdm2 = 100.0d0 * rsynp1go
380      oxjgzv0e = 0.001d0
381      dvhw1ulq = 1
382      kij0gwer=1
38323127 if(.not.(kij0gwer.le.zy1mchbf))goto 23129
384      ayfnwr1v=1
38523130 if(.not.(ayfnwr1v.le.n))goto 23132
386      vjz5sxty = nm0eljqk(ayfnwr1v,kij0gwer) / hdqsx7bk(ayfnwr1v,kij0gwe
387     *r)
388      if((vjz5sxty .lt. oxjgzv0e) .or. (nm0eljqk(ayfnwr1v,kij0gwer) .gt.
389     * 1.0d5))then
390      bzmd6ftv(ayfnwr1v,kij0gwer) = -nm0eljqk(ayfnwr1v,kij0gwer) * (1.0d
391     *0 + hdqsx7bk(ayfnwr1v,kij0gwer)/(hdqsx7bk(ayfnwr1v,kij0gwer) + nm0
392     *eljqk(ayfnwr1v,kij0gwer))) / hdqsx7bk(ayfnwr1v,kij0gwer)**2
393      if(bzmd6ftv(ayfnwr1v,kij0gwer) .gt. -btiehdm2)then
394      bzmd6ftv(ayfnwr1v,kij0gwer) = -btiehdm2
395      endif
396      goto 20
397      endif
398      q6zdcwxk = 0.0d0
399      pvcjl2na = hdqsx7bk(ayfnwr1v,kij0gwer) / (hdqsx7bk(ayfnwr1v,kij0gw
400     *er) + nm0eljqk(ayfnwr1v,kij0gwer))
401      mwuvskg1 = 1.0d0 - pvcjl2na
402      csi9ydge = hdqsx7bk(ayfnwr1v,kij0gwer)
403      if(pvcjl2na .lt. btiehdm2)then
404      pvcjl2na = btiehdm2
405      endif
406      if(mwuvskg1 .lt. btiehdm2)then
407      mwuvskg1 = btiehdm2
408      endif
409      esql7umk = 100.0d0 + 15.0d0 * nm0eljqk(ayfnwr1v,kij0gwer)
410      if(esql7umk .lt. sguwj9ty)then
411      esql7umk = sguwj9ty
412      endif
413      ft3ijqmy = pvcjl2na ** csi9ydge
414      ux3nadiw = ft3ijqmy
415      plo6hkdr = (1.0d0 - ux3nadiw) / hdqsx7bk(ayfnwr1v,kij0gwer)**2
416      q6zdcwxk = q6zdcwxk + plo6hkdr
417      ydb = 1.0d0
418      ft3ijqmy = hdqsx7bk(ayfnwr1v,kij0gwer) * mwuvskg1 * ft3ijqmy
419      ux3nadiw = ux3nadiw + ft3ijqmy
420      plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb
421     *)**2
422      q6zdcwxk = q6zdcwxk + plo6hkdr
423      ydb = 2.0d0
42423143 if(((ux3nadiw .le. n2kersmx) .or. (plo6hkdr .gt. 1.0d-4)) .and. (y
425     *db .lt. esql7umk))then
426      ft3ijqmy = (hdqsx7bk(ayfnwr1v,kij0gwer) - 1.0d0 + ydb) * mwuvskg1
427     ** ft3ijqmy / ydb
428      ux3nadiw = ux3nadiw + ft3ijqmy
429      plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb
430     *)**2
431      q6zdcwxk = q6zdcwxk + plo6hkdr
432      ydb = ydb + 1.0d0
433      goto 23143
434      endif
43523144 continue
436      bzmd6ftv(ayfnwr1v,kij0gwer) = -q6zdcwxk
43720    hmayv1xt = 0.0d0
43823131 ayfnwr1v=ayfnwr1v+1
439      goto 23130
44023132 continue
44123128 kij0gwer=kij0gwer+1
442      goto 23127
44323129 continue
444      return
445      end
446      subroutine enbin8(bzmd6ftv, hdqsx7bk, hsj9bzaq, n2kersmx, kuzxj1lo
447     *, dvhw1ulq, zy1mchbf, ux3nadiw, rsynp1go)
448      implicit logical (a-z)
449      integer kuzxj1lo, dvhw1ulq, zy1mchbf
450      double precision bzmd6ftv(kuzxj1lo, zy1mchbf), hdqsx7bk(kuzxj1lo,
451     *zy1mchbf), hsj9bzaq(kuzxj1lo, zy1mchbf), n2kersmx, ux3nadiw, rsynp
452     *1go
453      integer ayfnwr1v, kij0gwer, esql7umk
454      double precision ft3ijqmy, tad5vhsu, o3jyipdf, pq0hfucn, q6zdcwxk,
455     * d1, d2, plo6hkdr, hnu1vjyw
456      logical pok1, pok2, pok12
457      double precision oxjgzv0e, onemse, nm0eljqk, btiehdm2, ydb, kbig
458      d1 = 0.0d0
459      d2 = 0.0d0
460      btiehdm2 = -100.0d0 * rsynp1go
461      esql7umk = 3000
462      if(n2kersmx .le. 0.80d0 .or. n2kersmx .ge. 1.0d0)then
463      dvhw1ulq = 0
464      return
465      endif
466      kbig = 1.0d4
467      oxjgzv0e = 0.001d0
468      hnu1vjyw = 1.0d0 - rsynp1go
469      onemse = 1.0d0 / (1.0d0 + oxjgzv0e)
470      dvhw1ulq = 1
471      kij0gwer=1
47223147 if(.not.(kij0gwer.le.zy1mchbf))goto 23149
473      ayfnwr1v=1
47423150 if(.not.(ayfnwr1v.le.kuzxj1lo))goto 23152
475      if(hdqsx7bk(ayfnwr1v,kij0gwer) .gt. kbig)then
476      hdqsx7bk(ayfnwr1v,kij0gwer) = kbig
477      endif
478      if(hsj9bzaq(ayfnwr1v,kij0gwer) .lt. oxjgzv0e)then
479      hsj9bzaq(ayfnwr1v,kij0gwer) = oxjgzv0e
480      endif
481      if((hsj9bzaq(ayfnwr1v,kij0gwer) .gt. onemse))then
482      nm0eljqk = hdqsx7bk(ayfnwr1v,kij0gwer) * (1.0d0/hsj9bzaq(ayfnwr1v,
483     *kij0gwer) - 1.0d0)
484      bzmd6ftv(ayfnwr1v,kij0gwer) = -nm0eljqk * (1.0d0 + hdqsx7bk(ayfnwr
485     *1v,kij0gwer)/(hdqsx7bk(ayfnwr1v,kij0gwer) + nm0eljqk)) / hdqsx7bk(
486     *ayfnwr1v,kij0gwer)**2
487      if(bzmd6ftv(ayfnwr1v,kij0gwer) .gt. btiehdm2)then
488      bzmd6ftv(ayfnwr1v,kij0gwer) = btiehdm2
489      endif
490      goto 20
491      endif
492      q6zdcwxk = 0.0d0
493      pok1 = .true.
494      pok2 = hsj9bzaq(ayfnwr1v,kij0gwer) .lt. (1.0d0-rsynp1go)
495      pok12 = pok1 .and. pok2
496      if(pok12)then
497      d2 = hdqsx7bk(ayfnwr1v,kij0gwer) * dlog(hsj9bzaq(ayfnwr1v,kij0gwer
498     *))
499      ux3nadiw = dexp(d2)
500      else
501      ux3nadiw = 0.0d0
502      endif
503      plo6hkdr = (1.0d0 - ux3nadiw) / hdqsx7bk(ayfnwr1v,kij0gwer)**2
504      q6zdcwxk = q6zdcwxk + plo6hkdr
505      call tldz5ion(hdqsx7bk(ayfnwr1v,kij0gwer), o3jyipdf)
506      ydb = 1.0d0
507      call tldz5ion(ydb + hdqsx7bk(ayfnwr1v,kij0gwer), tad5vhsu)
508      pq0hfucn = 0.0d0
509      if(pok12)then
510      d1 = dlog(1.0d0 - hsj9bzaq(ayfnwr1v,kij0gwer))
511      ft3ijqmy = dexp(ydb * d1 + d2 + tad5vhsu - o3jyipdf - pq0hfucn)
512      else
513      ft3ijqmy = 0.0d0
514      endif
515      ux3nadiw = ux3nadiw + ft3ijqmy
516      plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb
517     *)**2
518      q6zdcwxk = q6zdcwxk + plo6hkdr
519      ydb = 2.0d0
52023165 if((ux3nadiw .le. n2kersmx) .or. (plo6hkdr .gt. 1.0d-4))then
521      tad5vhsu = tad5vhsu + dlog(ydb + hdqsx7bk(ayfnwr1v,kij0gwer) - 1.0
522     *d0)
523      pq0hfucn = pq0hfucn + dlog(ydb)
524      if(pok12)then
525      ft3ijqmy = dexp(ydb * d1 + d2 + tad5vhsu - o3jyipdf - pq0hfucn)
526      else
527      ft3ijqmy = 0.0d0
528      endif
529      ux3nadiw = ux3nadiw + ft3ijqmy
530      plo6hkdr = (1.0d0 - ux3nadiw) / (hdqsx7bk(ayfnwr1v,kij0gwer) + ydb
531     *)**2
532      q6zdcwxk = q6zdcwxk + plo6hkdr
533      ydb = ydb + 1.0d0
534      if(ydb .gt. 1.0d3)then
535      goto 21
536      endif
537      goto 23165
538      endif
53923166 continue
54021    bzmd6ftv(ayfnwr1v,kij0gwer) = -q6zdcwxk
54120    tad5vhsu = 0.0d0
54223151 ayfnwr1v=ayfnwr1v+1
543      goto 23150
54423152 continue
54523148 kij0gwer=kij0gwer+1
546      goto 23147
54723149 continue
548      return
549      end
550      subroutine mbessi0(bvecto, kuzxj1lo, kpzavbj3, d0, d1, d2, zjkrtol
551     *8, qaltf0nz)
552      implicit logical (a-z)
553      integer kuzxj1lo, kpzavbj3, zjkrtol8, c5aesxkus
554      double precision bvecto(kuzxj1lo), d0(kuzxj1lo), d1(kuzxj1lo), d2(
555     *kuzxj1lo), qaltf0nz
556      integer ayfnwr1v, gp1jxzuh
557      double precision f0, t0, m0, f1, t1, m1, f2, t2, m2
558      double precision toobig
559      toobig = 20.0d0
560      zjkrtol8 = 0
561      if(.not.(kpzavbj3 .eq. 0 .or. kpzavbj3 .eq. 1 .or. kpzavbj3 .eq. 2
562     *))then
563      zjkrtol8 = 1
564      return
565      endif
566      do23173 gp1jxzuh=1,kuzxj1lo
567      if(dabs(bvecto(gp1jxzuh)) .gt. toobig)then
568      zjkrtol8 = 1
569      return
570      endif
571      t1 = bvecto(gp1jxzuh) / 2.0d0
572      f1 = t1
573      t0 = t1 * t1
574      f0 = 1.0d0 + t0
575      t2 = 0.50d0
576      f2 = t2
577      c5aesxkus = 15
578      if(dabs(bvecto(gp1jxzuh)) .gt. 10)then
579      c5aesxkus = 25
580      endif
581      if(dabs(bvecto(gp1jxzuh)) .gt. 15)then
582      c5aesxkus = 35
583      endif
584      if(dabs(bvecto(gp1jxzuh)) .gt. 20)then
585      c5aesxkus = 40
586      endif
587      if(dabs(bvecto(gp1jxzuh)) .gt. 30)then
588      c5aesxkus = 55
589      endif
590      do23185 ayfnwr1v=1,c5aesxkus
591      m0 = (bvecto(gp1jxzuh) / (2.0d0*(ayfnwr1v+1.0d0))) ** 2.0
592      m1 = m0 * (1.0d0 + 1.0d0/ayfnwr1v)
593      m2 = m1 * (2.0d0*ayfnwr1v + 1.0d0) / (2.0d0*ayfnwr1v - 1.0d0)
594      t0 = t0 * m0
595      t1 = t1 * m1
596      t2 = t2 * m2
597      f0 = f0 + t0
598      f1 = f1 + t1
599      f2 = f2 + t2
600      if((dabs(t0) .lt. qaltf0nz) .and. (dabs(t1) .lt. qaltf0nz) .and. (
601     *dabs(t2) .lt. qaltf0nz))then
602      goto 23186
603      endif
60423185 continue
60523186 continue
606      if(0 .le. kpzavbj3)then
607      d0(gp1jxzuh) = f0
608      endif
609      if(1 .le. kpzavbj3)then
610      d1(gp1jxzuh) = f1
611      endif
612      if(2 .le. kpzavbj3)then
613      d2(gp1jxzuh) = f2
614      endif
61523173 continue
61623174 continue
617      return
618      end
619