42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "mvsiz_p.inc"
50
51
52
53#include "com08_c.inc"
54
55
56
57 INTEGER JLT, NOINT, IFUNCTK, FCOND,
58 . NPC(*),IFORM,ITAB(*),
59 . NSV(*),MSR(*),IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ)
60 my_real,
intent(in) :: theaccfact
63 . tempi(mvsiz), xi(mvsiz), yi(mvsiz),temp(*),
64 . zi(mvsiz), phi(mvsiz), areas(mvsiz),
65 . asi(mvsiz), bsi(mvsiz), gapv(mvsiz), pene(mvsiz),
66 . kthe, xthe, fni(mvsiz), tf(*), frad, drad,
67 . penrad(mvsiz), tempm(mvsiz),fheat,efrict(mvsiz),
68 . condint(mvsiz),
69 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
70 . phi1(mvsiz),phi2(mvsiz),phi3(mvsiz),phi4(mvsiz),
71 . ax1,ay1,az1,ax2,ay2,az2,ax,ay,az,areac,phim,
area,
72 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
73 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
74 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz)
76 . finter
77 EXTERNAL finter
78
79
80
81 INTEGER I
82
84 . ts, rstif, tstifm, tstift, dist, cond, p, dydx,
85 . tm, ddcond, dd, hcond
86
87 IF(ifunctk==0)THEN
88 rstif = one/
max(em30,kthe)
89 DO i=1,jlt
90 ts = tempi(i)
91 tm = tempm(i)
92 condint(i) = zero
93 ddcond =
max(dcond-gapv(i),em20)
94
95 IF(areas(i) == zero )THEN
96
97 ax1 = x3(i) - x1(i)
98 ay1 = y3(i) - y1(i)
99 az1 = z3(i) - z1(i)
100 ax2 = x4(i) - x2(i)
101 ay2 = y4(i) - y2(i)
102 az2 = z4(i) - z2(i)
103
104 ax = ay1*az2 - az1*ay2
105 ay = az1*ax2 - ax1*az2
106 az = ax1*ay2 - ay1*ax2
107
108 area = half*sqrt(ax*ax+ay*ay+az*az)
109
111
112 ELSE
113
114 areac = areas(i)
115
116 ENDIF
117
118
119
120
121
122
123
124 IF(penrad(i) <= zero)THEN
125
126 dist = penrad(i)+ gapv(i)
127
128
129
130 cond = asi(i)+bsi(i)*ts
131 tstifm =
max(dist,zero) / cond
132 tstift = tstifm + rstif
133 condint(i) = areac * theaccfact /tstift
134
135 phi(i) = areac*(tm - ts)*dt1*theaccfact / tstift
136
137
138
139 ELSEIF(penrad(i) <= ddcond)THEN
140
141 dist = gapv(i)
142 cond = asi(i)+bsi(i)*ts
143 tstifm =
max(dist,zero) / cond
144 tstift = tstifm + rstif
145 dd = penrad(i) /ddcond
146 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
147 condint(i) = areac*hcond * theaccfact
148
149 phi(i) = areac * (tm - ts)*dt1* hcond * theaccfact
150
151 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
152 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
153
154
155
156 ELSEIF(penrad(i) <= drad)THEN
157
158 phi(i) = frad * areac * (tm*tm+ts*ts)
159 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
160
161 END IF
162
163 IF(iform == 1 )THEN
164
165 phi1(i) = -phi(i) *h1(i)
166 phi2(i) = -phi(i) *h2(i)
167 phi3(i) = -phi(i) *h3(i)
168 phi4(i) = -phi(i) *h4(i)
169
170
171
172 phim = fheat * efrict(i)
173 phi1(i) = phi1(i) + phim*h1(i) !
main heating
174 phi2(i) = phi2(i) + phim*h2(i)
175 phi3(i) = phi3(i) + phim*h3(i)
176 phi4(i) = phi4(i) + phim*h4(i)
177 ENDIF
178
179
180
181 phi(i) = phi(i) + fheat * efrict(i) * theaccfact
182 ENDDO
183
184 ELSE
185 DO i=1,jlt
186 ts = tempi(i)
187 tm = tempm(i)
188 condint(i) = zero
189 ddcond =
max(dcond-gapv(i),em20)
190
191 IF(areas(i) == zero )THEN
192
193 ax1 = x3(i) - x1(i)
194 ay1 = y3(i) - y1(i)
195 az1 = z3
196 ax2 = x4(i) - x2
197 ay2 = y4
198 az2 = z4
199
200 ax = ay1*az2 - az1*ay2
201 ay = az1*ax2 - ax1*az2
202 az = ax1*ay2 - ay1*ax2
203
204 area = half*sqrt(ax*ax+ay*ay+az*az)
205
207
208 ELSE
209
210 areac = areas(i)
211
212 ENDIF
213
214
215
216
217
218
219
220 IF(penrad(i) <= zero)THEN
221
222 dist = penrad(i) +gapv(i)
223
224
225
226 p = xthe * abs(fni(i)) / areas(i)
227 rstif = one /
max(em30,kthe * finter(ifunctk,p,npc,tf,dydx))
228 cond = asi(i)+bsi(i)*ts
229 tstifm =
max(dist,zero) / cond
230 tstift = tstifm + rstif
231 condint(i) = areac * theaccfact/tstift
232
233 phi(i) = areac * (tm- ts)*dt1*theaccfact / tstift
234
235
236
237
238 ELSEIF(penrad(i) <= ddcond)THEN
239
240 dist = gapv(i)
241 cond = asi(i)+bsi(i)*ts
242 p = zero
243 rstif = one /
max(em30,kthe * finter(ifunctk,p,npc,tf,dydx))
244 tstifm =
max(dist,zero) / cond
245 tstift = tstifm + rstif
246 dd = penrad(i) /ddcond
247 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
248 condint(i) = areac*hcond*theaccfact
249
250 phi(i) = areac * (tm - ts)*dt1 * hcond *theaccfact
251
252 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
253 . * (tm + ts) * (tm - ts) * dt1 *theaccfact
254
255
256
257
258 ELSEIF(penrad(i) <= drad)THEN
259
260 phi(i) = frad * areac * (tm*tm+ts*ts)
261 . * (tm + ts) * (tm - ts) * dt1 *theaccfact
262 END IF
263
264 IF(iform == 1 )THEN
265
266 phi1(i) = -phi(i) *h1(i)
267 phi2(i) = -phi(i) *h2(i)
268 phi3(i) = -phi(i) *h3(i)
269 phi4(i) = -phi(i) *h4(i)
270
271
272
273 phim = fheat * efrict(i)
274 phi1(i) = phi1(i) + phim*h1(i)
275 phi2(i) = phi2(i) + phim*h2(i)
276 phi3(i) = phi3(i) + phim*h3(i)
277 phi4(i) = phi4(i) + phim*h4(i)
278 ENDIF
279
280
281
282 phi(i) = phi(i) + fheat * efrict(i)*theaccfact
283 ENDDO
284 END IF
285
286 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
int main(int argc, char *argv[])