Line data Source code
1 : MODULE m_outint
2 : CONTAINS
3 362489242 : SUBROUTINE outint(
4 520234 : > msh,e,fkap,cs,cis,s,vr,z,rn,rnot,h,d,
5 520234 : < a0,b0,a,b,
6 : < ki,nodes)
7 : c
8 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
9 : c-----x perform the outward integration. this routine is a x----
10 : c-----x derivative of start1/diff1. much unnecessary indexing x----
11 : c-----x has been removed among other things. dale koelling x----
12 :
13 : c with adams' predictor and corrector. more stable than milne's
14 : c for GGA.
15 : cc Handbook of math.func. p.896. T.Asada Feb.20,1998.
16 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
17 :
18 : IMPLICIT NONE
19 :
20 : c .. Arguments ..
21 : INTEGER,INTENT(IN) :: msh
22 : INTEGER,INTENT(OUT) :: ki,nodes
23 : REAL ,INTENT(IN) :: e,fkap,cs,cis,s,z,rn,rnot,h,d
24 : REAL ,INTENT(IN) :: vr(msh)
25 : REAL ,INTENT(OUT) :: a0(5),b0(5),a(msh),b(msh)
26 : c ..
27 : c .. local scalars ..
28 : REAL astr,atk,aw,bstr,btk,bw,da1,da2,da3,da4,db1,db2,db3,db4,g,gm,
29 : + gp,h3,q11,q22,r,rp12,rp21,test,vt
30 : INTEGER i,k,kim,n
31 : REAL :: da(5),db(5)
32 : c ..
33 : c .. intrinsic functions ..
34 : INTRINSIC abs,sign,sqrt
35 : c ..
36 : c.....------------------------------------------------------------------
37 : c .. equivalences ..
38 : EQUIVALENCE (da1,da(1)),(da2,da(2)),(da3,da(3)),(da4,da(4)),
39 : & (da5,da(5))
40 : EQUIVALENCE (db1,db(1)),(db2,db(2)),(db3,db(3)),(db4,db(4)),
41 : & (db5,db(5))
42 :
43 : REAL t18,t58,t98,t198,t378,t558,t598, da5,db5
44 : c.....------------------------------------------------------------------
45 :
46 520234 : h3 = h/3.0e0
47 520234 : n = msh
48 : cta+
49 520234 : t558=55.e0/8.e0
50 520234 : t598=59.e0/8.e0
51 520234 : t378=37.e0/8.e0
52 520234 : t98=9.e0/8.e0
53 :
54 520234 : t198=19.e0/8.e0
55 520234 : t58=5.e0/8.e0
56 520234 : t18=1.e0/8.e0
57 : cta-
58 520234 : g = sqrt(fkap**2- (cis*z)**2)
59 520234 : gp = g + s*fkap
60 520234 : gm = g - s*fkap
61 520234 : q11 = -gm
62 520234 : q22 = -gp
63 : c 25 april,1978
64 : c simple self-consistent starting procedure:
65 520234 : astr = sqrt(abs(gp))
66 520234 : bstr = sqrt(abs(gm))
67 520234 : r = rnot
68 3121404 : do 10 i = 1,5
69 2601170 : a(i) = astr
70 2601170 : b(i) = bstr
71 2601170 : rp21 = cis* (e*r-vr(i))
72 2601170 : rp12 = -rp21 - 2*cs*r
73 2601170 : a0(i) = h* (q11*astr+rp12*bstr)
74 2601170 : b0(i) = h* (q22*bstr+rp21*astr)
75 2601170 : r = r*d
76 520234 : 10 continue
77 6190826 : test = 1.e-6
78 6711060 : 20 vt = 0
79 6711060 : r = rnot
80 6711060 : r = rnot*d
81 33555300 : do 30 i = 2,5
82 26844240 : aw = a(i-1) + .5e0* (a0(i-1)+a0(i))
83 26844240 : bw = b(i-1) + .5e0* (b0(i-1)+b0(i))
84 26844240 : vt = vt + abs(aw-a(i))/astr + abs(bw-b(i))/bstr
85 26844240 : a(i) = 0.5e0* (a(i)+aw)
86 26844240 : b(i) = 0.5e0* (b(i)+bw)
87 26844240 : rp21 = cis* (e*r-vr(i))
88 26844240 : rp12 = -rp21 - 2*cs*r
89 26844240 : a0(i) = h* (q11*a(i)+rp12*b(i))
90 26844240 : b0(i) = h* (q22*b(i)+rp21*a(i))
91 26844240 : r = r*d
92 6711060 : 30 continue
93 6711060 : if(vt.gt.test) go to 20
94 : r = rnot
95 3121404 : do 40 i = 1,5
96 2601170 : rp21 = cis* (e*r-vr(i))
97 2601170 : rp12 = -rp21 - 2*cs*r
98 2601170 : da(i) = h3* (q11*a(i)+rp12*b(i))
99 2601170 : db(i) = h3* (q22*b(i)+rp21*a(i))
100 2601170 : r = r*d
101 520234 : 40 continue
102 : c end of starting procedure
103 :
104 520234 : nodes = 0
105 520234 : r = rn
106 :
107 : c..... find the classical turning point........
108 :
109 520234 : ki = n
110 103074228 : 50 ki = ki - 1
111 103074228 : if(ki.le.10) go to 60
112 103074228 : r = r/d
113 103074228 : if(e*r.lt.vr(ki)) go to 50
114 520234 : 60 ki = ki + 1
115 520234 : if(ki.ge.n) ki = ki - 1
116 520234 : kim = ki + 1
117 : if(kim.gt.n) kim = n
118 520234 : r = rnot* (d**3)
119 520234 : k = 4
120 361969008 : 70 k = k + 1
121 361969008 : r = r*d
122 361969008 : rp21 = cis* (e*r-vr(k))
123 361969008 : rp12 = -rp21 - 2.0e0*cs*r
124 :
125 : c adams' extrapolation formula for predictor.
126 361969008 : atk=a(k-1)+ t558*da4-t598*da3+t378*da2-t98*da1
127 361969008 : btk=b(k-1)+ t558*db4-t598*db3+t378*db2-t98*db1
128 :
129 361969008 : da1 = da2
130 361969008 : da2 = da3
131 361969008 : da3 = da4
132 :
133 361969008 : db1 = db2
134 361969008 : db2 = db3
135 361969008 : db3 = db4
136 :
137 361969008 : da4 = h3* (q11*atk+rp12*btk)
138 361969008 : db4 = h3* (q22*btk+rp21*atk)
139 :
140 : c adams interpolation formula for corrector.
141 361969008 : a(k)=a(k-1)+ t98*da4+t198*da3-t58*da2+t18*da1
142 361969008 : b(k)=b(k-1)+ t98*db4+t198*db3-t58*db2+t18*db1
143 :
144 361969008 : da4=h3*(q11*a(k)+rp12*b(k))
145 361969008 : db4=h3*(q22*b(k)+rp21*a(k))
146 :
147 361969008 : if(k.ge.kim) go to 80
148 :
149 361448774 : nodes = nodes + 0.501e0*abs(sign(1.0e0,a(k))-sign(1.0e0,a(k-1)))
150 361969008 : go to 70
151 :
152 : 80 continue
153 :
154 520234 : RETURN
155 : END SUBROUTINE
156 : END
|
