62
63
64
65 use glob_therm_mod
66
67
68
69#include "implicit_f.inc"
70
71
72
73#include "mvsiz_p.inc"
74
75
76
77#include "com08_c.inc"
78#include "param_c.inc"
79#include "sms_c.inc"
80#include "scr05_c.inc"
81
82
83
84 INTEGER, INTENT(IN) :: ITY
85 INTEGER, INTENT(IN) :: ISMSTR
86 INTEGER, INTENT(IN) :: JTUR
87 INTEGER, INTENT(IN) :: JTHE
88 INTEGER, INTENT(IN) :: JCVT
89 INTEGER, INTENT(IN) :: JSPH
90 INTEGER, INTENT(IN) :: JSMS,NPG
91
92 INTEGER NELTST,ITYPTST,PID(*),G_DT,NEL,ISELECT
93 INTEGER MAT(*),NGL(*),IPM(NPROPMI,*)
95 . dt2t
96
98 . pm(npropm,*), off(*), sig(nel,6), eint(*), rho(*), qold(*),
99 . vol(*),stifn(*), offg(*),geo(npropg,*),mumax(*),sigl(nel,6)
101 . vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*),
102 . psh(*), pnew(*),qnew(*) ,ssp_eq(*), dvol(*),
103 . sold1(*), sold2(*), sold3(*), sold4(*), sold5(*), sold6(*),
104 . d1(*), d2(*), d3(*), d4(*), d5(*), d6(*),
105 . mssa(*), dmels(*),conde(*),
106 . mfxx(*) ,mfxy(*) ,mfxz(*) ,mfyx(*) ,mfyy(*) ,
107 . mfyz(*) ,mfzx(*) ,mfzy(*) ,mfzz(*) ,offg0(*) ,vol_avg(*),
108 . epsth(*),
dtel(*),etotsh(nel,6), rhoref(*), rhosp(*),amu(*)
109 type (glob_therm_) ,intent(inout) :: glob_therm
110
111
112
113 INTEGER I, MX, J,IBID,NLAR,INDEX(MVSIZ),NSM,IPRES,ISTAB
115 . rho0(mvsiz),
116 . g(mvsiz), g1(mvsiz), g2(mvsiz),
117 . c1(mvsiz),
118 . es1(mvsiz), es2(mvsiz), es3(mvsiz),
119 . es4(mvsiz), es5(mvsiz), es6(mvsiz),
120 . df,dav,ekk, dpdm, p,
121 . e1, e2, e3, e4, e5, e6, einc,
122 . bid1, bid2, bid3, dta, ym, dpdmp,facq0,
123 . rho0_1,c1_1,g2_11,
124 . sv(3),p3,p2,nu_1,
125 . epsth11(mvsiz),sigtmp(mvsiz,6),facg,ff,frho
126
127 facq0 = one
128 mx = mat(1)
129 rho0_1 =pm( 1,mx)
130 c1_1 =pm(32,mx)
131 nu_1 =pm(21,mx)
132 g2_11 = two*pm(22,mx)
133 ipres=0
134 IF (nu_1>0.49)ipres=1
135 istab = 1
136 IF (ipres==1.OR.ismstr == 11) istab = 0
137 DO i=1,nel
138 rho0(i) =rho0_1
139 g(i) =pm(22,mx)*off(i)
140 c1(i) =c1_1
141 ENDDO
142
143
144 IF (ismstr == 10) THEN
145 IF (iselect>0) THEN
146 DO i=1,nel
147 es1(i)=etotsh(i,1)*off(i)
148 es2(i)=etotsh(i,2)*off(i)
149 es3(i)=etotsh(i,3)*off(i)
150 es4(i)=etotsh(i,4)*off(i)
151 es6(i)=etotsh(i,6)*off(i)
152 es5(i)=etotsh(i,5)*off(i)
153 ENDDO
154 ELSE
155 DO i=1,nel
156 es1(i)=(mfxx(i)*(two+mfxx(i))+mfxy(i)*mfxy(i)+mfxz(i)*mfxz(i))*off(i)
157 es2(i)=(mfyy(i)*(two+mfyy(i))+mfyx(i)*mfyx(i)+mfyz(i)*mfyz(i))*off(i)
158 es3(i)=(mfzz(i)*(two+mfzz(i))+mfzx(i)*mfzx(i)+mfzy(i)*mfzy(i))*off(i)
159 es4(i)=(mfxy(i)+mfyx(i)+mfxx(i)*mfyx(i)+mfxy(i)*mfyy(i)+mfxz(i)*mfyz(i))*off(i)
160 es6(i)=(mfxz(i)+mfzx(i)+mfxx(i)*mfzx(i)+mfxy(i)*mfzy(i)+mfxz(i)*mfzz(i))*off(i)
161 es5(i)=(mfzy(i)+mfyz(i)+mfzx(i)*mfyx(i)+mfzy(i)*mfyy(i)+mfzz(i)*mfyz(i))*off(i)
162 ENDDO
163 END IF
165 1 es5 ,es6 ,epsth, c1 ,g ,sigtmp)
166 IF(ipres==1.AND.iresp==1)THEN
167 DO i=1,nel
168 ff = -
min(sig(i,1),sig(i,2),sig(i,3))
169 facg =
max(one,ff/g2_11)
170 IF (facg>one) facg = onep2*facg
171 sigtmp(i,4)= facg*sigtmp(i,4)
172 sigtmp(i,5)= facg*sigtmp(i,5)
173 sigtmp(i,6)= facg*sigtmp(i,6)
174 g(i) = facg*g(i)
175 END DO
176 END IF
177 sig(1:nel,1:6)=sigtmp(1:nel,1:6)
178 ELSEIF (ismstr == 12) THEN
179 nlar = 0
180 IF (iselect>0) THEN
181 DO i=1,nel
182 IF (offg0(i)<=one.AND.off(i)/=zero) THEN
183 nlar =nlar+1
184 index(nlar)=i
185 es1(nlar)=etotsh(i,1)
186 es2(nlar)=etotsh(i,2)
187 es3(nlar)=etotsh(i,3)
188 es4(nlar)=etotsh(i,4)
189 es6(nlar)=etotsh(i,6)
190 es5(nlar)=etotsh(i,5)
191 epsth11(nlar)=epsth(i)
192 END IF
193 ENDDO
194 ELSE
195 DO i=1,nel
196 IF (offg0(i)<=one.AND.off(i)/=zero) THEN
197 nlar =nlar+1
198 index(nlar)=i
199 es1(nlar)=mfxx(i)*(two+mfxx(i))+mfxy(i)*mfxy(i)+mfxz(i)*mfxz(i)
200 es2(nlar)=mfyy(i)*(two+mfyy(i))+mfyx(i)*mfyx(i)+mfyz(i)*mfyz(i)
201 es3(nlar)=mfzz(i)*(two+mfzz(i))+mfzx(i)*mfzx(i)+mfzy(i)*mfzy(i)
202 es4(nlar)=mfxy(i)+mfyx(i)+mfxx(i)*mfyx(i)+mfxy(i)*mfyy(i)+mfxz(i)*mfyz(i)
203 es6(nlar)=mfxz(i)+mfzx(i)+mfxx(i)*mfzx(i)+mfxy(i)*mfzy(i)+mfxz(i)*mfzz(i)
204 es5(nlar)=mfzy(i)+mfyz(i)+mfzx(i)*mfyx(i)+mfzy(i)*mfyy(i)+mfzz(i)*mfyz(i)
205 epsth11(nlar)=epsth(i)
206 END IF
207 ENDDO
208 END IF
209
211 1 es5 ,es6 ,epsth11, c1 ,g ,sigtmp)
212#include "vectorize.inc"
213 DO j=1,nlar
214 i = index(j)
215 sig(i,1) =sigtmp(j,1)
216 sig(i,2) =sigtmp(j,2)
217 sig(i,3) =sigtmp(j,3)
218 sig(i,4) =sigtmp(j,4)
219 sig(i,5) =sigtmp(j,5)
220 sig(i,6) =sigtmp(j,6)
221 ENDDO
222
223 IF (nlar < nel) THEN
224
225 nsm = 0
226 IF (jcvt>0) THEN
227 DO i=1,nel
228 IF (offg0(i)>one) THEN
229 nsm = nsm +1
230 index(nsm) = i
231 es1(nsm)=mfxx(i)
232 es2(nsm)=mfyy(i)
233 es3(nsm)=mfzz(i)
234 es4(nsm)=mfxy(i)+mfyx(i)
235 es6(nsm)=mfxz(i)+mfzx(i)
236 es5(nsm)=mfzy(i)+mfyz(i)
237 epsth11(nsm)=epsth(i)
238 END IF
239 ENDDO
240 DO i=1,nsm
241 ekk=es1(i)+es2(i)+es3(i)-epsth11(i)
242 dav=-third*(es1(i)+es2(i)+es3(i))
243 sigtmp(i,1)=c1_1*ekk+g2_11*(es1(i)+dav)
244 sigtmp(i,2)=c1_1*ekk+g2_11*(es2(i)+dav)
245 sigtmp(i,3)=c1_1*ekk+g2_11*(es3(i)+dav)
246 sigtmp(i,4)=g(i)*es4(i)
247 sigtmp(i,5)=g(i)*es5(i)
248 sigtmp(i,6)=g(i)*es6(i)
249 ENDDO
250 IF(ipres==1)THEN
251#include "vectorize.inc"
252 DO j=1,nsm
253 i = index(j)
254 ff = -
min(sig(i,1),sig(i,2),sig(i,3))
255 facg =
max(one,ff/g2_11)
256 IF (facg>one) facg = onep2*facg
257 frho = rhoref(i)/rho0_1
258 IF (frho<onep01) facg = one
259 sig(i,1)=sigl(i,1) + sigtmp(j,1)
260 sig(i,2)=sigl(i,2) + sigtmp(j,2)
261 sig(i,3)=sigl(i,3) + sigtmp(j,3)
262 sig(i,3)=sigl(i,4) + facg*sigtmp(j,4)
263 sig(i,5)=sigl(i,5) + facg*sigtmp(j,5)
264 sig(i,6)=sigl(i,6) + facg*sigtmp(j,6)
265 g(i) = facg*g(i)
266 END DO
267 ELSE
268#include "vectorize.inc"
269 DO j=1,nsm
270 i = index(j)
271 sig(i,1)=sigl(i,1) + sigtmp(j,1)
272 sig(i,2)=sigl(i,2) + sigtmp(j,2)
273 sig(i,3)=sigl(i,3) + sigtmp(j,3)
274 sig(i,4)=sigl(i,4) + sigtmp(j,4)
275 sig(i,5)=sigl(i,5) + sigtmp(j,5)
276 sig(i,6)=sigl(i,6) + sigtmp(j,6)
277 END DO
279 ELSE
280 IF (iselect>0) THEN
281 DO i=1,nel
282 IF (offg0(i)>one) THEN
283 nsm = nsm +1
284 index(nsm) = i
285 es1(nsm)=etotsh(i,1)
286 es2(nsm)=etotsh(i,2)
287 es3(nsm)=etotsh(i,3)
288 es4(nsm)=etotsh(i,4)
289 es6(nsm)=etotsh(i,6)
290 es5(nsm)=etotsh(i,5)
291 epsth11(nsm)=epsth(i)
292 END IF
293 ENDDO
294 ELSE
295 DO i=1,nel
296 IF (offg0(i)>one) THEN
297 nsm = nsm +1
298 index(nsm) = i
299 es1(nsm)=mfxx(i)*(two+mfxx(i))+mfxy(i)*mfxy(i)+mfxz(i)*mfxz(i)
300 es2(nsm)=mfyy(i)*(two+mfyy(i))+mfyx(i)*mfyx(i)+mfyz(i)*mfyz(i)
301 es3(nsm)=mfzz(i)*(two+mfzz(i))+mfzx(i)*mfzx(i)+mfzy(i)*mfzy(i)
302 es4(nsm)=mfxy(i)+mfyx(i)+mfxx(i)*mfyx(i)+mfxy(i)*mfyy(i)+mfxz(i)*mfyz(i)
303 es6(nsm)=mfxz(i)+mfzx(i)+mfxx(i)*mfzx(i)+mfxy(i)*mfzy(i)+mfxz(i)*mfzz(i)
304 es5(nsm)=mfzy(i)+mfyz(i)+mfzx(i)*mfyx(i)+mfzy(i)*mfyy(i)+mfzz(i)*mfyz(i)
305 epsth11(nsm)=epsth(i)
306 END IF
307 END DO
308 END IF
310 1 es1 , es2 ,es3 ,es4 ,es5 ,
311 2 es6 ,epsth11, c1_1 ,g2_11 ,sigtmp )
312#include "vectorize.inc"
313 DO j=1,nsm
314 i = index(j)
315 sig(i,1)=sigl(i,1) + sigtmp(j,1)
316 sig(i,2)=sigl(i,2) + sigtmp(j,2)
317 sig(i,3)=sigl(i,3) + sigtmp(j,3)
318 sig(i,4)=sigl(i,4) + sigtmp(j,4)
319 sig(i,5)=sigl(i,5) + sigtmp(j,5)
320 sig(i,6)=sigl(i,6) + sigtmp(j,6)
321 pnew(i)= c1_1*amu(i)
322 p =-third*(sigtmp(j,1)+sigtmp(j,2)+sigtmp(j,3))
323 sig(i,1)= sig(i,1)-pnew(i) + p
324 sig(i,2)= sig(i,2)-pnew(i) + p
325 sig(i,3)= sig(i,3)-pnew(i) + p
326 END DO
327 END IF
328 END IF
329
330 ELSEIF (ismstr==11) THEN
331 DO i=1,nel
332 es1(i)=mfxx(i)
333 es2(i)=mfyy(i)
334 es3(i)=mfzz(i)
335 es4(i)=mfxy(i)+mfyx(i)
336 es6(i)=mfxz(i)+mfzx(i)
337 es5(i)=mfzy(i)+mfyz(i)
338 ENDDO
339 DO i=1,nel
340 g1(i)=g(i)
341 g2(i)=two*g1(i)
342 ekk=es1(i)+es2(i)+es3(i)-epsth(i)
343 dav=-third*(es1(i)+es2(i)+es3(i))
344 sig(i,1)=c1(i)*ekk+g2(i)*(es1(i)+dav)
345 sig(i,2)=c1(i)*ekk+g2(i)*(es2(i)+dav)
346 sig(i,3)=c1(i)*ekk+g2(i)*(es3(i)+dav)
347 sig(i,4)=g1(i)*es4(i)
348 sig(i,5)=g1(i)*es5(i)
349 sig(i,6)=g1(i)*es6(i)
350 ENDDO
353 . d2,d3 ,d4 ,d5 ,d6 ,
354 . rho0,npg ,nel)
355 ENDIF
356 IF (istab>0) THEN
357 IF (jcvt==0) THEN
358 IF (npg>1) THEN
361 . d2 ,d3 ,d4 ,d5 ,d6 ,
362 . rhoref ,g ,npg ,nel)
363 END IF
364 ELSE
365
366 IF (npg>1) THEN
369 . d2 ,d3 ,d4 ,d5 ,d6 ,
370 . rhoref,g ,npg ,nel )
371 ELSE
374 . d2 ,d3 ,d4 ,d5 ,d6 ,
375 . rhoref,g ,npg ,nel )
376 END IF
377 END IF
378 ENDIF
379 IF(ismstr==12)THEN
380 DO i=1,nel
381 IF(abs(offg0(i)) <= one)THEN
382 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho(i))
383 rhosp(i)=rho(i)
384 ELSE
385 ssp(i)=sqrt((onep333*g(i)+c1(i))/rhoref(i))
386 rhosp(i)=rhoref(i)
387 END IF
388 ENDDO
389 ELSEIF(idtmins/=2.OR.jsms==0)THEN
390
391 DO i=1,nel
392
393
394
395
396 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
397 rhosp(i)=rho0(i)
398 ENDDO
399 ELSEIF(ismstr==11)THEN
400 DO i=1,nel
401
402
403
404
405 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
406 rhosp(i)=rho0(i)
407 ENDDO
408 ELSE
409 DO i=1,nel
410 IF(abs(offg0(i)) <= one)THEN
411
412
413
414
415 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho(i))
416 rhosp(i)=rho(i)
417 ELSE
418
419
420
421
422 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
423 rhosp(i)=rho0(i)
424 END IF
425 ENDDO
426 END IF
427
428 IF (jsph==0)THEN
430 1 pm, off, rho, bid1,
431 2 bid2, ssp, bid3, stifn,
432 3 dt2t, neltst, ityptst, aire,
433 4 offg, geo, pid, vnew,
434 5 vd2, deltax, vis, d1,
435 6 d2, d3, pnew, psh,
436 7 mat, ngl, qnew, ssp_eq,
437 8 vol, mssa, dmels, ibid,
438 9 facq0, conde,
dtel, g_dt,
439 a ipm, rhoref, rhosp, nel,
440 b ity, ismstr, jtur, jthe,
441 c jsms, npg , glob_therm)
442 ELSE
444 1 pm, off, rho, bid1,
445 2 bid2, bid3, stifn, dt2t,
446 3 neltst, ityptst, offg, geo,
447 4 pid, mumax, ssp, vnew,
448 5 vd2, deltax, vis, d1,
449 6 d2, d3, pnew, psh,
450 7 mat, ngl, qnew, ssp_eq,
451 8 g_dt,
dtel, nel, ity,
452 9 jtur, jthe)
453 ENDIF
454
455 dta =half*dt1
456
457 DO i=1,nel
458 sig(i,1)=sig(i,1)*off(i)
459 sig(i,2)=sig(i,2)*off(i)
460 sig(i,3)=sig(i,3)*off(i)
461 sig(i,4)=sig(i,4)*off(i)
462 sig(i,5)=sig(i,5)*off(i)
463 sig(i,6)=sig(i,6)*off(i)
464 p2 = -(sold1(i)+sig(i,1)+sold2(i)+sig(i,2)+sold3(i)+sig(i,3))* third
465 e1=d1(i)*(sold1(i)+sig(i,1)+p2)
466 e2=d2(i)*(sold2(i)+sig(i,2)+p2)
467 e3=d3(i)*(sold3(i)+sig(i,3)+p2)
468 e4=d4(i)*(sold4(i)+sig(i,4))
469 e5=d5(i)*(sold5(i)+sig(i,5))
470 e6=d6(i)*(sold6(i)+sig(i,6))
471 einc= vol_avg(i)*(e1+e2+e3+e4+e5+e6)*dta - half*dvol(i)*(qold(i)+qnew(i)+p2)
472 eint(i)=(eint(i)+einc*off(i)) /
max(em15,vol(i))
473 qold(i)=qnew(i)
474 ENDDO
475
476 RETURN
if(complex_arithmetic) id
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine m1ismstr10_pon(nel, es1, es2, es3, es4, es5, es6, epsth, c1, g, sig)
subroutine m1ismstr11(nel, es1, es2, es3, es4, es5, es6, epsth, c1, g2, sig11)
subroutine m1tot_stab11(sig, g2, g, nel)
subroutine m1tot_stab18(sig, sigl, g2, g, offg, ismstr, nel)
subroutine m1tot_stab24(sig, sigl, g2, g, offg, ismstr, nel)
subroutine m1tot_stab_p(sig, g2, c1, nel)
subroutine mdtsph(pm, off, rho, rk, t, re, sti, dt2t, neltst, ityptst, offg, geo, pid, mumax, ssp, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, g_dt, dtsph, nel, ity, jtur, jthe)
subroutine mqviscb(pm, off, rho, rk, temp, ssp, re, sti, dt2t, neltst, ityptst, aire, offg, geo, pid, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, vol0, mssa, dmels, igeo, facq0, conde, dtel, g_dt, ipm, rhoref, rhosp, nel, ity, ismstr, jtur, jthe, jsms, npg, glob_therm)
subroutine nsvis_stab11(sig, c1, g2, vol, d1, d2, d3, d4, d5, d6, rhoref, npg, nel)
subroutine nsvis_stab18(sig, c1, g2, vol, d1, d2, d3, d4, d5, d6, rhoref, g, npg, nel)
subroutine nsvis_stab(sig, c1, g2, vol, d1, d2, d3, d4, d5, d6, rhoref, g, npg, nel)