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