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