41
42
43
44#include "implicit_f.inc"
45
46
47
48#include "mvsiz_p.inc"
49
50
51
52#include "param_c.inc"
53#include "com08_c.inc"
54
55
56
57 INTEGER JLT, , FCOND,ITAB(*) ,NSV(*),
58 . NPC(*),,IELESI(MVSIZ) ,IELEMI(MVSIZ),
59 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ)
60 my_real,
intent(in) :: theaccfact
62 . pm(npropm,*)
64 . kthe, xthe, frad, drad, fheats, fheatm,
65 . tempi(mvsiz), xi(mvsiz), yi(mvsiz),temp(*),
66 . zi(mvsiz), phi(mvsiz), areas(mvsiz), asi(mvsiz),
67 . bsi(mvsiz), gapv(mvsiz), condint(mvsiz),
68 . fni(mvsiz), tf(*), efrict(mvsiz),
69 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
70 . phi1(mvsiz),phi2(mvsiz),phi3(mvsiz),phi4(mvsiz),
71 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
72 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
73 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz)
75 . finter
76 EXTERNAL finter
77
78
79
80 INTEGER I, MATS , MATM
81
83 . sx1 , sy1 , sz1 , sx2 , sy2 , sz2,
84 . ts, rstif, tstifm, tstift, cond, p, dydx,
85 . tm, ddcond, dd, hcond, dist, penrad,
norm, areac,
86 . phim,
area,conds,condm
87
88 phi(1:jlt) = zero
89 IF(ifunctk==0)THEN
90 rstif = one/
max(em30,kthe)
91 DO i=1,jlt
92 ts = tempi(i)
93 IF(iform ==0) THEN
94 tm = tint
95 ELSE
96 tm = h1(i)*temp(ix1(i)) + h2(i)*temp(ix2(i))
97 . + h3(i)*temp(ix3(i)) + h4(i)*temp(ix4(i))
98 ENDIF
99 condint(i) = zero
100 ddcond =
max(dcond-gapv(i),em20)
101
102
103
104
105 sx1=(y1(i)-y3(i))*(z2(i)-z4(i)) - (z1(i)-z3(i))*(y2(i)-y4(i))
106 sy1=(z1(i)-z3(i))*(x2(i)-x4(i)) - (x1(i)-x3(i))*(z2(i)-z4(i))
107 sz1=(x1(i)-x3(i))*(y2(i)-y4(i)) - (y1(i)-y3(i))*(x2(i)-x4(i))
108
109 norm = sqrt(sx1**2 + sy1**2 + sz1**2)
110
111
112
113 IF(ix3(i)/=ix4(i))THEN
114 sx2 = fourth*(x1(i) + x2(i) + x3(i) + x4(i)) - xi(i)
115 sy2 = fourth*(y1(i) + y2(i) + y3(i) + y4(i)) - yi(i)
116 sz2 = fourth*(z1(i) + z2(i) + z3(i) + z4(i)) - zi(i)
117 ELSE
118 sx2 = third
119 sy2 = third*(y1(i) + y2(i) + y3(i)) - yi(i)
120 sz2 = third*(z1(i) + z2(i) + z3(i)) - zi(i)
121 END IF
122
123
124
125
126
127 dist = -(sx2*sx1+sy2*sy1+sz2*sz1) /
max(em15,
norm)
128
129 penrad = dist - gapv(i)
130
131 IF(areas(i) == zero )THEN
132
134
135 ELSE
136
137 areac = areas(i)
138
139 ENDIF
140 tstifm = zero
141 cond = zero
142
143
144
145
146
147
148
149
150 IF(penrad <= zero)THEN
151
152
153
154
155 mats = ielesi(i)
156 IF(mats > 0) THEN
157 conds=pm(75,mats)+pm(76,mats)*ts
158 ELSE
159 conds =zero
160 ENDIF
161 IF(iform == 1 ) THEN
162 matm = ielemi(i)
163 IF(matm > 0) THEN
164 condm=pm(75,matm)+pm(76,matm)*tm
165 ELSE
166 condm = zero
167 ENDIF
168 IF(condm == zero) THEN
169 cond = conds
170 ELSEIF(conds == zero) THEN
171 cond = condm
172 ELSE
173
174 ENDIF
175 ELSE
176 cond = conds
177 ENDIF
178 IF(cond /= zero) tstifm = abs(dist) / cond
179 tstift = tstifm + rstif
180 condint(i) = areac * theaccfact /tstift
181
182 phi(i) = areac*(tm - ts)*dt1*theaccfact
183
184
185
186 ELSEIF(penrad <= ddcond)THEN
187
188 dist = gapv(i)
189
190 mats = ielesi(i)
191 IF(mats > 0) THEN
192 conds=pm(75,mats)+pm(76,mats)*ts
193 ELSE
194 conds =zero
195 ENDIF
196
197 IF(iform == 1 ) THEN
198 matm = ielemi(i)
199 IF(matm > 0) THEN
200 condm=pm(75,matm)+pm(76,matm)*tm
201 ELSE
202 condm = zero
203 ENDIF
204 IF(condm == zero) THEN
205 cond = conds
206 ELSEIF(conds == zero) THEN
207 cond = condm
208 ELSE
209 cond = two*conds*condm/(condm + conds)
210 ENDIF
211 ELSE
212 cond = conds
213 ENDIF
214
215 IF(cond /= zero) tstifm =
max(dist,zero) / cond
216 tstift = tstifm + rstif
217 dd = penrad /ddcond
218 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
219 condint(i) = areac*hcond*theaccfact
220
221 phi(i) = areac * (tm - ts)*dt1* hcond*theaccfact
222
223 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
224 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
225
226
227
228 ELSEIF(penrad <= drad)THEN
229
230 phi(i) = frad * areac * (tm*tm+ts*ts)
231 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
232 END IF
233
234 IF(iform == 1 )THEN
235
236 phi1(i) = -phi(i) *h1(i)
237 phi2(i) = -phi(i) *h2(i)
238 phi3(i) = -phi(i) *h3(i)
239 phi4(i) = -phi(i) *h4(i)
240
241
242
243 phim = fheatm * efrict(i)
244 phi1(i) = phi1(i) + phim*h1(i)
245 phi2(i) = phi2(i) + phim*h2(i)
246 phi3(i) = phi3(i) + phim*h3(i)
247 phi4(i) = phi4(i) + phim*h4(i)
248
249 ENDIF
250
251
252
253 phi(i) = phi(i) + fheats * efrict(i) * theaccfact
254 ENDDO
255
256 ELSE
257 DO i=1,jlt
258 ts = tempi(i)
259 IF(iform ==0) THEN
260 tm = tint
261 ELSE
262 tm = h1(i)*temp(ix1(i)) + h2(i)*temp(ix2(i))
263 . + h3(i)*temp(ix3(i)) + h4(i)*temp(ix4(i))
264 ENDIF
265 condint(i) = zero
266 ddcond =
max(dcond-gapv(i),em20)
267
268
269
270
271 sx1=(y1(i)-y3(i))*(z2(i)-z4(i)) - (z1(i)-z3(i))*(y2(i)-y4(i))
272 sy1=(z1(i)-z3(i))*(x2(i)-x4(i)) - (x1(i)-x3(i))*(z2(i)-z4(i))
273 sz1=(x1(i)-x3(i))*(y2(i)-y4(i)) - (y1(i)-y3(i))*(x2(i)-x4(i))
274
275 norm = sqrt(sx1**2 + sy1**2 + sz1**2)
276
277
278
279 IF(ix3(i)/=ix4(i))THEN
280 sx2 = fourth*(x1(i) + x2(i) + x3(i) + x4(i)) - xi(i)
281 sy2 = fourth*(y1(i) + y2(i) + y3(i) + y4(i)) - yi(i)
282 sz2 = fourth*(z1(i) + z2
283 ELSE
284 sx2 = third*(x1(i) + x2(i) + x3(i)) - xi(i)
285 sy2 = third*(y1(i) + y2(i) + y3(i)) - yi(i)
286 sz2 = third*(z1(i) + z2(i) + z3(i)) - zi(i)
287 END IF
288
289
290
291
292
293 dist = -(sx2*sx1+sy2*sy1+sz2*sz1) /
max(em15,
norm)
294
295 penrad = dist - gapv(i)
296
297 IF(areas(i) == zero )THEN
298
300
301 ELSE
302
303 areac = areas(i)
304
305 ENDIF
306 tstifm = zero
307 cond = zero
308
309
310
311 IF(penrad <= zero)THEN
312
313
314
315 p = xthe * abs(fni(i)) / areas(i)
316 rstif = one /
max(em30,kthe * finter(ifunctk,p,npc
317 mats = ielesi(i)
318 IF(mats > 0) THEN
319 conds=pm(75,mats)+pm(76,mats)*ts
320 ELSE
321 conds =zero
322 ENDIF
323
324 IF(iform == 1 ) THEN
325 matm = ielemi(i)
326 IF(matm > 0) THEN
327 condm=pm(75,matm)+pm(76,matm)*tm
328 ELSE
329 condm = zero
330 ENDIF
331 IF(condm == zero) THEN
332 cond = conds
333 ELSEIF(conds == zero) THEN
334 cond = condm
335 ELSE
336 cond = two*conds*condm/(condm + conds)
337 ENDIF
338 ELSE
339 cond = conds
340 ENDIF
341
342 IF(cond /= zero) tstifm = abs(dist) / cond
343 tstift = tstifm + rstif
344 condint(i) = areac*theaccfact /tstift
345
346 phi(i) = areac * (tm- ts)*dt1 * theaccfact / tstift
347
348
349
350
351 ELSEIF(penrad <= ddcond)THEN
352
353 dist = gapv(i)
354 mats = ielesi(i)
355 IF(mats > 0) THEN
356 conds=pm(75,mats)+pm(76,mats)*ts
357 ELSE
358 conds =zero
359 ENDIF
360
361 IF(iform == 1 ) THEN
362 matm = ielemi(i)
363 IF(matm > 0) THEN
364 condm=pm(75,matm)+pm(76,matm)*tm
365 ELSE
366 condm = zero
367 ENDIF
368 IF(condm == zero) THEN
369 cond = conds
370 ELSEIF(conds == zero) THEN
371 cond = condm
372 ELSE
373 cond = two*conds*condm/(condm + conds)
374 ENDIF
375 ELSE
376 cond = conds
377 ENDIF
378 p = zero
379 rstif = one /
max(em30,kthe * finter(ifunctk,p,npc,tf,dydx))
380
381 IF(cond /= zero)tstifm =
max(dist,zero) / cond
382 tstift = tstifm + rstif
383
384 dd = penrad /ddcond
385 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
386 condint(i) = areac*hcond*theaccfact
387
388 phi(i) = areac * (tm - ts)*dt1 * hcond * theaccfact
389
390 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
391 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
392
393
394
395
396 ELSEIF(penrad <= drad)THEN
397
398 phi(i) = frad * areac * (tm*tm+ts*ts)
399 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
400 END IF
401
402 IF(iform == 1 )THEN
403
404 phi1(i) = -phi(i) *h1(i)
405 phi2(i) = -phi(i) *h2(i)
406 phi3(i) = -phi(i) *h3(i)
407 phi4(i) = -phi(i) *h4(i)
408
409
410
411 phim = fheatm * efrict(i)
412 phi1(i) = phi1(i) + phim*h1(i)
413 phi2(i) = phi2(i) + phim*h2(i)
414 phi3(i) = phi3(i) + phim*h3(i)
415 phi4(i) = phi4(i) + phim*h4(i)
416 ENDIF
417
418
419
420 phi(i) = phi(i) + fheats * efrict(i) * theaccfact
421 ENDDO
422 END IF
423
424 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
subroutine area(d1, x, x2, y, y2, eint, stif0)