45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57#include "com04_c.inc"
58#include "com06_c.inc"
59#include "param_c.inc"
60#include "scr02_c.inc"
61#include "scr06_c.inc"
62#include "scr14_c.inc"
63#include "scr16_c.inc"
64#include "scr17_c.inc"
65#include "sms_c.inc"
66
67
68
69 INTEGER IXC(NIXC,*),IPARTC(*), JFT, JLT,PID(*),
70 . IHBE ,NFT ,ISMSTR,IGEO(NPROPGI, *),NEL,MTN
71 INTEGER MAT(MVSIZ)
72
74 . pm(npropm,*), geo(npropg,*), thk(*), hour(nel,5), off
75 . px1(*), px2(*), py1(*), py2(*),dt1c(*),eani(*),
76 . ssp(mvsiz), rho(mvsiz),sti(mvsiz),stir(*),
77 . h1(mvsiz), h2(mvsiz),
78 . thk0(mvsiz),viscmx(mvsiz), alpe(mvsiz),partsav(npsav,*)
79
81 . b11(mvsiz), b12(mvsiz), b13(mvsiz), b14(mvsiz), b21(mvsiz),
82 . b22(mvsiz), b23(mvsiz), b24(mvsiz), h11(mvsiz), h12(mvsiz),
83 . h13(mvsiz), h14(mvsiz), h21(mvsiz), h22(mvsiz), h23(mvsiz),
84 . h24(mvsiz), h31(mvsiz), h32(mvsiz), h33(mvsiz), h34
85 . rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz), ry1(mvsiz),
86 . ry2(mvsiz), ry3(mvsiz), ry4(mvsiz), vhx(mvsiz), vhy(mvsiz),
87 . vx1(mvsiz), vx2(mvsiz), vx3(mvsiz),
88 . vx4(mvsiz), vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz),
89 . vz1(mvsiz), vz2(mvsiz), vz3(mvsiz), vz4(mvsiz),
area(mvsiz),
90 . a11r(mvsiz)
91 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: a1
92
93
94
95 INTEGER I, MX,IPID,IGTYP,IGMAT,IPGMAT
96 my_real h1l(mvsiz), h2l(mvsiz), h1q(mvsiz), h2q(mvsiz), hg1(mvsiz),
97 . hg2(mvsiz), fac(mvsiz), ym(mvsiz), pr(mvsiz),
98 . gama1(mvsiz), gama2(mvsiz), gama3(mvsiz), gama4(mvsiz),
99 . h4(mvsiz), h4l(mvsiz), h4q(mvsiz),thk02(mvsiz),
100 . g(mvsiz) , b1(mvsiz), b2(mvsiz),a11(mvsiz),ehou(mvsiz),
101 . px1v, px2v, py1v, py2v, ehourt,vv ,scale(mvsiz),fac1
102
103 IF(ismstr /=3 .AND. ihbe >= 1)THEN
104 DO i=jft,jlt
105 px1v = px1(i)*vhx(i)
106 px2v = px2(i)*vhx(i)
107 py1v = py1(i)*vhy(i)
108 py2v = py2(i)*vhy(i)
109 gama1(i)= off(i)*( one- px1v-py1v)
110 gama3(i)= off(i)*( one+ px1v+py1v)
111 gama2(i)= off(i)*(-one- px2v-py2v)
112 gama4(i)= off(i)*(-one+ px2v+py2v)
113 ENDDO
114 ELSE
115 DO i=jft,jlt
116 gama1(i)= off(i)
117 gama3(i)= off(i)
118 gama2(i)= -off(i)
119 gama4(i)= -off(i)
120 ENDDO
121 ENDIF
122
123 IF(invstr >= 35)THEN
124 mx=mat(jft)
125 DO i=jft,jlt
126 ipid=pid(i)
127 h4(i) =geo(17,ipid)
128 ENDDO
129
130 ELSE
131 mx=mat(jft)
132 DO i=jft,jlt
133 h4(i) =pm(91,mx)
134 ENDDO
135
136 ENDIF
137
138 DO i=jft,jlt
139 fac(i)=fourth*rho(i)*thk0(i)
140 h1l(i)= zero
141 h1q(i)= zero
142 h2l(i)= zero
143 h2q(i)= zero
144 h4l(i)=fac(i)*sqrt(hvisc*h4(i)*
area(i))
145 h4q(i)=sqrt(hvisc*h4(i))*h4l(i)*hundred
146 h4l(i)=h4l(i)*ssp(i)
147 ENDDO
148
149 DO i=jft,jlt
150 thk02(i)= thk0(i)*thk0(i)
151 b1(i) = px1(i)*px1(i)+py1(i)*py1(i)
152 b2(i) = px2(i)*px2(i)+py2(i)*py2(i)
153 fac(i)=fourth*ym(i)*thk0(i)*dt1c(i)*helas
154 h1(i)=h1(i)*fac(i)
155 h2(i)=h2(i)*fac(i)
156 ENDDO
157
158
159
160 DO i=jft,jlt
161 IF(ixc(4,i)/=ixc(5,i))cycle
162 h1q(i)=zero
163 h1l(i)=zero
164 h2q(i)=zero
165 h2l(i)=zero
166 h4l(i)=zero
167 h4q(i)=zero
168 h1(i)=zero
169 h2(i)=zero
170 ENDDO
171
172
173
174 DO i=jft,jlt
175 scale(i)= zero
176 ENDDO
177 ipid=pid(1)
178 igtyp = igeo(11,ipid)
179 igmat = igeo(98,ipid)
180 ipgmat = 700
181 IF(nodadt /= 0 .OR. idt1sh == 1.OR. idtmins == 2)THEN
182
183 DO i=jft,jlt
184 scale(i)=
max(gama1(i)*gama1(i),gama2(i)*gama2(i),gama3(i)*gama3(i),gama4(i)*gama4(i)) *
185 . dt1c(i)*
max(h1(i)+h1l(i),h2(i)+h2l(i),h4l(i)) /
max(dt1c(i)*dt1c(i),em20)
186 sti(i)=sti(i) + scale(i)
187 ENDDO
188 IF(igtyp == 11 .AND. igmat > 0) THEN
189 DO i=jft,jlt
190 a11(i) = geo(ipgmat +5 ,pid(i))
191 g(i) = geo(ipgmat+4,pid(i))
192 a11r(i) = geo(ipgmat+7,pid(i))
193 IF (off(i)==zero) THEN
194 sti(i) = zero
195 stir(i) = zero
196 ELSE
197 vv = viscmx(i) * viscmx(i) * alpe(i)
198 fac1 =
max(b1(i),b2(i)) / (
area(i) * vv)
199 sti(i) = sti(i) + fac1* thk0(i) * a11(i)
200 stir(i) = fac1 * a11r(i)*one_over_12*thk0(i)**3 +
201 . fac1 * a11(i)*thk0(i)*
area(i)*one_over_9 +
202 . fac1*scale(i)*(one_over_12*thk0(i)**2 +
area(i)*one_over_9)
203 ENDIF
204 ENDDO
205 ELSE
206 mx = mat(jft)
207 DO i=jft,jlt
208 a11(i) =pm(24,mx)
209 g(i) =pm(22,mx)
210 ENDDO
211 IF (mtn==58) a11(jft:jlt)=a1(jft:jlt)
212 DO i=jft,jlt
213 IF (off(i)==zero) THEN
214 sti(i) = zero
215 stir(i) = zero
216 ELSE
217 vv = viscmx(i) * viscmx(i) * alpe(i)
218 sti(i) = sti(i) +
max(b1(i),b2(i)) * thk0(i) * a11(i) / (
area(i) * vv)
219 stir(i) = sti(i)*(thk02(i) * one_over_12 +
area(i) * one_over_9)
220 ENDIF
221 ENDDO
222 ENDIF
223 ENDIF
224
225
226
227 IF(ismstr == 3 .OR. ihbe < 1)THEN
228 DO i=jft,jlt
229 hg1(i)=vx1(i)-vx2(i)+vx3(i)-vx4(i)
230 hg2(i)=vy1(i)-vy2(i)+vy3(i)-vy4(i)
231 ENDDO
232
233 DO i=jft,jlt
234 hour(i,1)=hour(i,1)+hg1(i)*h1(i)
235 hour(i,2)=hour(i,2)+hg2(i)*h1(i)
236 hg1(i)=hg1(i)*(h1l(i)+h1q(i)*abs(hg1(i)))
237 hg2(i)=hg2(i)*(h1l(i)+h1q(i)*abs(hg2(i)))
238 h11(i)= hour(i,1)+hg1(i)
239 h12(i)=-hour(i,1)-hg1(i)
240 h13(i)= hour(i,1)+hg1(i)
241 h14(i)=-hour(i,1)-hg1(i)
242 h21(i)= hour(i,2)+hg2(i)
243 h22(i)=-hour(i,2)-hg2(i)
244 h23(i)= hour(i,2)+hg2(i)
245 h24(i)=-hour(i,2)-hg2(i)
246 ENDDO
247 ELSE
248 DO i=jft,jlt
249 hg1(i)=vx1(i)*gama1(i)+vx2(i)*gama2(i)+vx3(i)*gama3(i)+vx4(i)*gama4(i)
250 hg2(i)=vy1(i)*gama1(i)+vy2(i)*gama2(i)+vy3(i)*gama3(i)+vy4(i)*gama4(i)
251 ENDDO
252 DO i=jft,jlt
253 hour(i,1)=hour(i,1)+hg1(i)*h1(i)
254 hour(i,2)=hour(i,2)+hg2(i)*h1(i)
255 hg1(i)=hg1(i)*(h1l(i)+h1q(i)*abs(hg1(i)))
256 hg2(i)=hg2(i)*(h1l(i)+h1q(i)*abs(hg2(i)))
257 h11(i)=(hour(i,1)+hg1(i))*gama1(i)
258 h12(i)=(hour(i,1)+hg1(i))*gama2(i)
259 h13(i)=(hour(i,1)+hg1(i))*gama3(i)
260 h14(i)=(hour(i,1)+hg1(i))*gama4(i)
261 h21(i)=(hour(i,2)+hg2(i))*gama1(i)
262 h22(i)=(hour(i,2)+hg2(i))*gama2(i)
263 h23(i)=(hour(i,2)+hg2(i))*gama3(i)
264 h24(i)=(hour(i,2)+hg2(i))*gama4(i)
265 ENDDO
266 ENDIF
267 ehourt = 0.
268 DO i=jft,jlt
269 ehou(i) = vx1(i)*h11(i) + vx2(i)*h12(i) + vx3(i)*h13(i) + vx4(i)*h14(i)
270 . + vy1(i)*h21(i) + vy2(i)*h22(i) + vy3(i)*h23(i) + vy4(i)*h24(i)
271 ENDDO
272
273
274
275 IF(ismstr==3.OR.ihbe<1)THEN
276 DO i=jft,jlt
277 hg1(i)=vz1(i)-vz2(i)+vz3(i)-vz4(i)
278 ENDDO
279
280 DO i=jft,jlt
281 hour(i,3)=hour(i,3)+hg1(i)*h2(i)
282 hg1(i)=hg1(i)*(h2l(i)+h2q(i)*abs(hg1(i)))
283 h31(i)= hour(i,3)+hg1(i)
284 h32(i)=-hour(i,3)-hg1(i)
285 h33(i)= hour(i,3)+hg1(i)
286 h34(i)=-hour(i,3)-hg1(i)
287 ENDDO
288 ELSE
289 DO i=jft,jlt
290 hg1(i)=vz1(i)*gama1(i)+vz2(i)*gama2(i)+vz3(i)*gama3(i)+vz4(i)*gama4(i)
291 ENDDO
292 DO i=jft,jlt
293 hour(i,3)=hour(i,3)+hg1(i)*h2(i)
294 hg1(i)=hg1(i)*(h2l(i)+h2q(i)*abs(hg1(i)))
295 h31(i)=(hour(i,3)+hg1(i))*gama1(i)
296 h32(i)=(hour(i,3)+hg1(i))*gama2(i)
297 h33(i)=(hour(i,3)+hg1(i))*gama3(i)
298 h34(i)=(hour(i,3)+hg1(i))*gama4(i)
299 ENDDO
300 ENDIF
301
302
303
304 DO i=jft,jlt
305 hg1(i)=+vz1(i)+vz2(i)-vz3(i)-vz4(i)
306 ENDDO
307
308 DO i=jft,jlt
309 hg1(i)=hg1(i)*(h4l(i)+h4q(i)*abs(hg1(i)))
310 h31(i)=h31(i) +hg1(i)
311 h32(i)=h32(i) +hg1(i)
312 h33(i)=h33(i) -hg1(i)
313 h34(i)=h34(i) -hg1(i)
314 ENDDO
315
316 DO i=jft,jlt
317 hg1(i)=vz1(i)-vz2(i)-vz3(i)+vz4(i)
318 ENDDO
319
320 DO i=jft,jlt
321 hg1(i)=hg1(i)*(h4l(i)+h4q(i)*abs(hg1(i)))
322 h31(i)=h31(i) +hg1(i)
323 h32(i)=h32(i) -hg1(i)
324 h33(i)=h33(i) -hg1(i)
325 h34(i)=h34(i) +hg1(i)
326 ENDDO
327
328 DO i=jft,jlt
329 b11(i)= zero
330 b12(i)= zero
331 b13(i)= zero
332 b14(i)= zero
333 b21(i)= zero
334 b22(i)= zero
335 b23(i)= zero
336 b24(i)= zero
337 ENDDO
338
339 DO i=jft,jlt
340 ehou(i) = ehou(i) + vz1(i)*h31(i) + vz2(i)*h32(i) + vz3(i)*h33(i) + vz4(i)*h34(i)
341 ehou(i) = dt1c(i) * ehou(i)
342 ehourt = ehourt + ehou(i)
343 ENDDO
344
345 DO i=jft,jlt
346 mx = ipartc(i)
347 partsav(8,mx)=partsav(8,mx) + ehou(i)
348 ENDDO
349
350
351 ehour = ehour + ehourt
352
353 DO i=jft,jlt
354 eani(nft+numels+i) = eani(nft+numels+i)+ehou(i)
355 ENDDO
356 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)