41
42
43
45
46
47
48#include "implicit_f.inc"
49#include "param_c.inc"
50
51
52
53 INTEGER, INTENT(IN) :: NMN,NSN
54 INTEGER IRECT(4,*), LMSR(4,*), MSR(*), NSV(NSN), ILOC(*), IRTL(*),LCODE(*), ISKEW(*), ITAB(*)
56 . x(3,*), skew(lskew,*), a(*), fn(*), ft(*), msmn(*), msmt(*),
57 . crst(2,*), v(*), nor(3,*), ms(*), efric(*), fsav(nthvki),
58 . frigap(*)
59
60
61
62#include "com01_c.inc"
63#include "com08_c.inc"
64#include "scr08_a_c.inc"
65
66
67
68 INTEGER JBC(3), NIR, I, J, I3, J3, I2, J2, I1, J1, ISK, LCOD, II,
69 . L, JJ, NN, LK, IBC
71 . xn(3), yn(3), zn(3), h(4), n1, n2, n3, ax, ay, az,
72 . vx, vy, vz, amn, vmn, amod, vmod, bvz, baz, bvx, bvy, bax, bay,
73 . a11, a12, a13, a21, a22, a23, a31, a32, a33, axt, ayt, azt,
74 . vxt, vyt, vzt, axn, ayn, azn, vt, at, ftt, fac, vxn, vyn,
75 . vzn, fnn, det, fric, fheat
76
77 fric = frigap(1)
78 fheat= frigap(4)
79 nir=2
80 IF(n2d == 0)nir=4
81
82
83
84 DO i=1,nmn
85 j=msr(i)
86 i3=3*i
87 j3=3*j
88 i2=i3-1
89 j2=j3-1
90 i1=i2-1
91 j1=j2-1
92 fsav(1)=fsav(1)+fn(i1)*dt12
93 fsav(2)=fsav(2)+fn(i2)*dt12
94 fsav(3)=fsav(3)+fn(i3)*dt12
95 fsav(4)=fsav(4)+ft(i1)*dt12
96 fsav(5)=fsav(5)+ft(i2)*dt12
97 fsav(6)=fsav(6)+ft(i3)*dt12
98 IF(msmn(i1) > zero)THEN
99 a(j1)=a(j1) + fn(i1)/msmn(i1)
100 a(j2)=a(j2) + fn(i2)/msmn(i1)
101 a(j3)=a(j3) + fn(i3)/msmn(i1)
102 ENDIF
103 IF(msmt(i1) > zero)THEN
104 a(j1)=a(j1) + ft(i1)/msmt(i1)
105 a(j2)=a(j2) + ft(i2)/msmt(i1)
106 a(j3)=a(j3) + ft(i3)/msmt(i1)
107 ENDIF
108 isk=iskew(j)
109 lcod=lcode(j)
110 CALL bcs2(a(j1),skew(1,isk),isk,lcod)
111 ENDDO
112
113
114
115 DO ii=1,nsn
116 IF(iloc(ii) >= 1) THEN
117 i=nsv(ii)
118 l=irtl(ii)
119 DO jj=1,nir
120 nn=irect(jj,l)
121 ix(jj)=msr(nn)
122 ENDDO
123 n1=nor(1,ii)
124 n2=nor(2,ii)
125 n3=nor(3,ii)
126 IF(n2d == 0)THEN
127 CALL shapeh(h,crst(1,ii),crst(2,ii))
128 ELSE
129 h(1) = half*(one - crst(1,ii))
130 h(2) = half*(one + crst(1,ii))
131 ENDIF
132 i3=3*i
133 i2=i3-1
134 i1=i2-1
135 ax=zero
136 ay=zero
137 az=zero
138 vx=zero
139 vy=zero
140 vz=zero
141 DO jj=1,nir
142 j3=3*ix(jj)
143 j2=j3-1
144 j1=j2-1
145 ax=ax+a(j1)*h(jj)
146 ay=ay+a(j2)*h(jj)
147 az=az+a(j3)*h(jj)
148 vx=vx+v(j1)*h(jj)
149 vy=vy+v(j2)*h(jj)
150 vz=vz+v(j3)*h(jj)
151 ENDDO
152
153 amn=n1*ax+n2*ay+n3*az
154 vmn=n1*vx+n2*vy+n3*vz
155 amod=amn-n1*a(i1)-n2*a(i2)-n3*a(i3)
156 vmod=vmn-n1*v(i1)-n2*v(i2)-n3*v(i3)
157
158 axn = amod*n1
159 ayn = amod*n2
160 azn = amod*n3
161 vxn = vmod*n1
162 vyn = vmod*n2
163 vzn = vmod*n3
164 fnn = (vmod/dt12 + amod) * ms(i)
165
166 axt = ax - a(i1) - axn
167 ayt = ay - a(i2) - ayn
168 azt = az - a(i3) - azn
169 vxt = vx - v(i1) - vxn
170 vyt = vy - v(i2) - vyn
171 vzt = vz - v(i3) - vzn
172 vt = sqrt(vxt**2+vyt**2+vzt**2)
173 at = sqrt(axt**2+ayt**2+azt**2)
174 ftt = (vt/dt12 + at) * ms(i)
175
176 fac =
min(one,fric*fnn/
max(em30,ftt))
177 ftt = ftt*fac
178 axt = axt*fac
179 ayt = ayt*fac
180 azt = azt*fac
181 vxt = vxt*fac
182 vyt = vyt*fac
183 vzt = vzt*fac
184
185
186
187 a(i1)=a(i1)+axn+axt
188 a(i2)=a(i2)+ayn+ayt
189 a(i3)=a(i3)+azn+azt
190 v(i1)=v(i1)+vxn+vxt
191 v(i2)=v(i2)+vyn+vyt
192 v(i3)=v(i3)+vzn+vzt
193
194 lcod=lcode(i)
195 IF(lcod > 0)THEN
196
197
198
199 xn(1)=n1
200 yn(1)=n2
201 zn(1)=n3
202
203 jbc(1:3) = 0
204 IF(lcode(i) > 0) THEN
205 jbc(3) = iand(lcode(i), 1)
206 jbc(2) = iand(lcode(i), 2)
207 jbc(1) = iand(lcode(i), 4)
208 ENDIF
209
210 lk=iskew(i)
211 ibc=2
212
213 IF(jbc(1) /= 0)THEN
214 xn(ibc)=skew(1,lk)
215 yn(ibc)=skew(2,lk)
216 zn(ibc)=skew(3,lk)
217 ibc=ibc+1
218 ENDIF
219 IF(jbc(2) /= 0)THEN
220 xn(ibc)=skew(4,lk)
221 yn(ibc)=skew(5,lk)
222 zn(ibc)=skew(6,lk)
223 ibc=ibc+1
224 ENDIF
225 IF(jbc(3) /= 0)THEN
226 IF(ibc==4)THEN
227
228 CALL ancmsg(msgid=11,anmode=aninfo,i1=itab(i))
230 ENDIF
231 xn(ibc)=skew(7,lk)
232 yn(ibc)=skew(8,lk)
233 zn(ibc)=skew(9,lk)
234 ibc=ibc+1
235 ENDIF
236 IF(ibc==3)THEN
237
238
239
240 xn(3)=yn(1)*zn(2)-zn(1)*yn(2)
241 yn(3)=zn(1)*xn(2)-xn(1)*zn(2)
242 zn(3)=xn(1)*yn(2)-yn(1)*xn(2)
243 bvz=v(i1)*xn(3)+v(i2)*yn(3)+v(i3)*zn(3)
244 baz=a(i1)*xn(3)+a(i2)*yn(3)+a(i3)*zn(3)
245 ELSE
246
247
248
249 bvz=zero
250 baz=zero
251 ENDIF
252
253 bvx=vmn
254 bvy=zero
255 bax=amn
256 bay=zero
257
258 a11=yn(2)*zn(3)-zn(2)*yn(3)
259 a12=zn(2)*xn(3)-xn(2)*zn(3)
260 a13=xn(2)*yn(3)-yn(2)*xn(3)
261 a21=yn(3)*zn(1)-zn(3)*yn(1)
262 a22=zn(3)*xn(1)-xn(3)*zn(1)
263 a23=xn(3)*yn(1)-yn(3)*xn(1)
264 a31=yn(1)*zn(2)-zn(1)*yn(2)
265 a32=zn(1)*xn(2)-xn(1)*zn(2)
266 a33=xn(1)*yn(2)-yn(1)*xn(2)
267 det=xn(1)*a11+yn(1)*a12+zn(1)*a13
268
269
270
271 IF(det /= zero) THEN
272 v(i1)=(a11*bvx+a21*bvy+a31*bvz)/det
273 v(i2)=(a12*bvx+a22*bvy+a32*bvz)/det
274 v(i3)=(a13*bvx+a23*bvy+a33*bvz)/det
275
276 a(i1)=(a11*bax+a21*bay+a31*baz)/det
277 a(i2)=(a12*bax+a22*bay+a32*baz)/det
278 a(i3)=(a13*bax+a23*bay+a33*baz)/det
279 ENDIF
280
281 ENDIF
282
283
284
285
286 vt = sqrt((v(i1)-vx)**2 + (v(i2)-vy)**2 + (v(i3)-vz)**2 )
287 efric(ii) = fheat * ftt * vt * dt1
288
289 ENDIF
290
291 ENDDO
292
293 RETURN
subroutine bcs2(a, b, j, k)
subroutine shapeh(h, s, t)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)