34
35
36
37#include "implicit_f.inc"
38
39
40
41#include "com08_c.inc"
42#include "param_c.inc"
43
44
45
46 INTEGER NC, ISKIP, NCF_S,
47 . IMPCNC(*),IMPCNN(*),IMPCDL(*),IMPCSK(*),LLL(*),JLL(*),
48 . SLL(*),IADLL(*),COMNTAG(*),ICFTAG(*),JCFTAG(*)
49
51 . xll(*),skew(lskew,*),rbmpc(*),ms(*),in(*),v(3,*),vr(3,*),
52 . a(3,*),ar(3,*)
53
54
55
56 INTEGER I, J, JJ, LL, IK, NK, KF, ISK, NN, NDL, IDL, NUMC, IC
58
59
60
61
62
63
64
65
66
67
68
69
70 s1 = 0
71 s2 = 0
72 kf = 0
73 nn = 0
74 s3 = -huge(s3)
75 DO i=1,nummpc
76 hh = zero
77 r = zero
78 nc = nc + 1
79 ik = iadll(nc)-1
80 numc = impcnc(i)
81 nk = 0
82 DO j=1,numc
83 kf = kf+1
84 nn = impcnn(kf)
85 ndl = impcdl(kf)
86 isk = impcsk(kf)
87 coef= rbmpc(kf)
88
89 IF(isk==1)THEN
90 nk = nk + 1
91 ik = ik + 1
92 lll(ik) = nn
93 jll(ik) = ndl
94 sll(ik) = 0
95 xll(ik) = coef
96 IF (ndl>3) THEN
97 ndl = ndl - 3
98 hh = hh + coef*coef / in(nn)
99 r = r + coef*(vr(ndl,nn) / dt12 + ar(ndl,nn))
100 ELSE
101 hh = hh + coef*coef / ms(nn)
102 r = r + coef*(v(ndl,nn) / dt12 + a(ndl,nn))
103 ENDIF
104 ELSE
105 idl = ndl
106 IF (ndl>3) idl = ndl - 3
107 IF (idl==1) THEN
108 s1 = coef*skew(1,isk)
109 s2 = coef*skew(2,isk)
110 s3 = coef*skew(3,isk)
111 ELSEIF (idl==2) THEN
112 s1 = coef*skew(4,isk)
113 s2 = coef*skew(5,isk)
114 s3 = coef*skew(6,isk)
115 ELSEIF (idl==3) THEN
116 s1 = coef*skew(7,isk)
117 s2 = coef*skew(8,isk)
118 s3 = coef*skew(9,isk)
119 ENDIF
120 nk = nk + 3
121 IF (ndl>3) THEN
122 ik = ik + 1
123 lll(ik) = nn
124 jll(ik) = 4
125 sll(ik) = 0
126 xll(ik) = s1
127 ik = ik + 1
128 lll(ik) = nn
129 jll(ik) = 5
130 sll(ik) = 0
131 xll(ik) = s2
132 ik = ik + 1
133 lll(ik) = nn
134 jll(ik) = 6
135 sll(ik) = 0
136 xll(ik) = s3
137 hh = hh + (s1*s1 + s2*s2 + s3*s3) / in(nn)
138 r = r + s1*(vr(1,nn) / dt12 + ar(1,nn))
139 . + s2*(vr(2,nn) / dt12 + ar(2,nn))
140 . + s3*(vr(3,nn) / dt12 + ar(3,nn))
141 ELSE
142 ik = ik + 1
143 lll(ik) = nn
144 jll(ik) = 1
145 sll(ik) = 0
146 xll(ik) = s1
147 ik = ik + 1
148 lll(ik) = nn
149 jll(ik) = 2
150 sll(ik) = 0
151 xll(ik) = s2
152 ik = ik + 1
153 lll(ik) = nn
154 jll(ik) = 3
155 sll(ik) = 0
156 xll(ik) = s3
157 hh = hh + (s1*s1 + s2*s2 + s3*s3) / ms(nn)
158 r = r + s1*(v(1,nn) / dt12 + a(1,nn))
159 . + s2*(v(2,nn) / dt12 + a(2,nn))
160 . + s3*(v(3,nn) / dt12 + a(3,nn))
161 ENDIF
162 ENDIF
163 ENDDO
164 iadll(nc+1) = iadll(nc) + nk
165
166
167 IF (hh/=zero) r = r / hh
168 DO ik=iadll(nc),iadll(nc+1)-1
169 jj = jll(ik)
170 ll = lll(ik)
171 IF (jj>3) THEN
172 jj = jj - 3
173 ar(jj,ll) = ar(jj,ll) - xll(ik)*r / in(ll)
174 ELSE
175 a(jj,ll) = a(jj,ll) - xll(ik)*r / ms(ll)
176 ENDIF
177 ENDDO
178 IF (comntag(nn)==1) THEN
179 iskip = iskip + 1
180 nc = nc - 1
181 ELSE
182 ic = nc - ncf_s
183 icftag(ic) = ic + iskip
184 jcftag(ic+iskip) = nc
185 ENDIF
186 ENDDO
187
188 RETURN