41 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC ,
42 2 NPF ,TF , TIME , TIMESTEP, UPARAM,
43 3 RHO0 , RHO ,VOLUME , EINT , NGL ,
44 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
45 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
46 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
47 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
48 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
49 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
50 A MFXX ,MFXY ,MFXZ,MFYX, MFYY , MFYZ ,
52 C SOUNDSP, VISCMAX, UVAR , OFF , ISMSTR, ET ,
53 D IHET ,OFFG , EPSTH3 , IEXPAN, NPARF ,
54 E UPARAMF,UVARF ,NVARF ,JCVT ,GAMA_R)
58#include "implicit_f.inc"
66 INTEGER ,
INTENT(IN) :: NPARF
67 INTEGER ,
INTENT(IN) :: NEL,JCVT,NVARF,NUPARAM,NUVAR,ISMSTR,IHET,IEXPAN
68 INTEGER ,
INTENT(IN) :: NGL(NEL)
69 my_real :: TIME,TIMESTEP
70 my_real ,
DIMENSION(NUPARAM) ,
INTENT(IN) :: UPARAM
71 my_real ,
DIMENSION(NPARF) ,
INTENT(IN) :: UPARAMF
72 my_real ,
DIMENSION(NEL) :: RHO,RHO0,VOLUME,EINT,EPSTH3,
73 . EPSPXX, EPSPYY, EPSPZZ, EPSPXY, EPSPYZ, EPSPZX,
74 . DEPSXX, DEPSYY, DEPSZZ, DEPSXY, DEPSYZ, DEPSZX,
75 . EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
76 . SIGOXX, SIGOYY, SIGOZZ, SIGOXY, SIGOYZ, SIGOZX,OFFG,
77 . MFXX,MFXY,MFXZ,MFYX,MFYY,MFYZ,MFZX,MFZY,MFZZ
79 my_real :: GAMA_R(NEL,6)
84 . signxx(nel), signyy(nel), signzz(nel),
85 . signxy(nel), signyz(nel), signzx(nel),
86 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
87 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
88 . soundsp(nel), viscmax(nel), et(nel)
93 . uvar(nel,nuvar), off(nel) ,uvarf(nel,nvarf)
97 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
98 my_real FINTER,FINTTE,TF(*),FINT2V
99 EXTERNAL FINTER,FINTTE
103 INTEGER I,J,K,KK,LL,,DIRECT,ITER,NITER,FLAG_MUL,IAVIS,
106 . et1,et2,et3,a1,a10,ksi,mu,d,lam,g,rbulk,aa,bb,cc,sb, temp,
107 . temp2,maxl,stiff0,dsig,deps,coef1,coef2,coef3,coef4,coef5,coef6,
108 . aa1,b1,c1,dd1,aa2,b2,c2,dd2,aa3,b3,c3,dd3,inv12,inv22,inv32,
109 . evb12,evb22,evb32,jdet2,jm1
110 . c10,c01,c20,c11,c02,c30,c21,c12
111 . bi1(nel),bi2(nel),jdet(nel),stiff(nel),
112 . ta1(nel), ta2(nel),ta3(nel),t1(nel), t2(nel),t3(nel),lpchain(nel),
113 . tb1(nel), tb2(nel),tb3(nel),trace(nel),traceb(nel),eta(nel),etb(nel),
114 . sb1(nel), sb2(nel),sb3(nel),sb4(nel), sb5(nel),sb6(nel),
115 . tbnorm(nel),dgamma(nel),
116 . devb1(nel),devb2(nel),devb3(nel),devb4(nel),devb5(nel),devb6(nel
117 . r1x(nel),r1y(nel),r1z(nel),r2x(nel),r2y(nel),r2z(nel),r3x(nel),r3y(nel),r3z(nel),
119 . f(nel,3,3),ft(nel,3,3),fe(nel,3,3),invfe(nel,3,3),fet(nel,3,3),fp(nel,3,3),
120 . fft(nel,3,3),invfpo(nel,3,3),matb(nel,3,3),fpo(nel,3,3),bpo(nel,3,3),invsn(nel,3,3),
121 . sig(nel,3,3),sigb(nel,3,3),siga(nel,3,3),sn(nel,3,3),fedp(nel,3,3),
122 . dfp(nel,3,3),lb(nel,3,3),dfp2(nel,3,3),fpdot(nel,3,3),
123 . caii(nel,3),cbii(nel,3),c_max(nel),cii(3)
153 iform = nint(uparam(23))
155 stiff0 = four_over_3*g + rbulk
158 IF (a10*sb==zero) iavis=0
160 flag_mul = uparam(21)
161 IF (flag_mul > zero)
THEN
171 fpo(i,1,1) = uvar(i,1)
172 fpo(i,2,2) = uvar(i,2)
173 fpo(i,3,3) = uvar(i,3)
174 fpo(i,1,2) = uvar(i,4)
175 fpo(i,2,3) = uvar(i,5)
176 fpo(i,3,1) = uvar(i,6)
177 fpo(i,2,1) = uvar(i,7)
178 fpo(i,3,2) = uvar(i,8)
179 fpo(i,1,3) = uvar(i,9)
183 r1x(1:nel) = gama_r(1:nel,1)
184 r1y(1:nel) = gama_r(1:nel,2)
185 r1z(1:nel) = gama_r(1:nel,3)
186 r2x(1:nel) = gama_r(1:nel,4)
187 r2y(1:nel) = gama_r(1:nel,5)
188 r2z(1:nel) = gama_r(1:nel,6)
190 r3x(1:nel) = r1y(1:nel)* r2z(1:nel) - r1z(1:nel) * r2y(1:nel)
191 r3y(1:nel) = r1z(1:nel)* r2x(1:nel) - r1x(1:nel) * r2z(1:nel)
192 r3z(1:nel) = r1x(1:nel)* r2y(1:nel) - r1y(1:nel) * r2x(1:nel)
194 r3r3 = sqrt(r3x(i)*r3x(i) + r3y(i)*r3y(i) + r3z(i)*r3z(i))
195 IF (r3r3 /= zero)
THEN
202 . r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
208 f(i,1,1) = one+mfxx(i)
209 f(i,2,2) = one+mfyy(i)
210 f(i,3,3) = one+mfzz(i)
223 !-------------------------
228 1 nel , fft , c10 , c01, c20,
229 2 c11 ,c02 , c30 , c21, c12,
230 3 c03 ,d1 , d2 , d3, siga ,
231 4 bi1 ,bi2 ,jdet ,flag_mul,
232 5 nvarf,coefr,betaf,coefm ,
233 6 uvarf,caii ,rbulk,iform)
252 CALL prodmat (f , invfpo, fe, nel)
260 1 nel , matb , c10, c01, c20,
261 2 c11 ,c02, c30, c21, c12,
262 3 c03 ,d1 ,d2 , d3, sigb ,
263 4 bi1,bi2,jdet ,flag_mul,
264 5 nvarf,coefr, betaf,coefm ,
269 sigb(i,1,1) = sb * sigb(i,1,1)
270 sigb(i,2,2) = sb * sigb(i,2,2)
271 sigb(i,3,3) = sb * sigb(i,3,3)
272 sigb(i,1,2) = sb * sigb(i,1,2)
273 sigb(i,2,3) = sb * sigb(i,2,3)
274 sigb(i,3,1) = sb * sigb(i,3,1)
279 traceb(i) = third*(sigb(i,1,1) +sigb(i,2,2) + sigb(i,3,3))
281 sb1(i) = sigb(i,1,1) - traceb(i)
282 sb2(i) = sigb(i,2,2) - traceb(i)
283 sb3(i) = sigb(i,3,3) - traceb(i)
288 tbnorm(i) = sqrt(
max(em20,sb1(i)**2+sb2(i)**2+sb3(i)**2
289 . + two*(sb4(i)**2+sb5(i)**2+sb6(i)**2 )) )
295 CALL viscbb ( nel , fpo, tbnorm , a1 , expc,
296 . expm, ksi, tauref , dgamma )
299 IF (time == zero) dgamma(1:nel) = uvar(1:nel,10)
301 uvar(i,10)= dgamma(i)
302 factor = dgamma(i)/tbnorm(i)
303 lb(i,1,1) = factor *sb1(i)
304 lb(i,2,2) = factor *sb2(i)
305 lb(i,3,3) = factor *sb3(i)
306 lb(i,1,2) = factor *sb4(i)
307 lb(i,2,3) = factor *sb5(i)
308 lb(i,3,1) = factor *sb6(i)
309 lb(i,2,1) = lb(i,1,2)
310 lb(i,3,2) = lb(i,2,3)
311 lb(i,1,3) = lb(i,3,1)
319 CALL prodmat (lb , fe, fedp, nel) ! fedp = lb * fe
320 CALL prodmat(invfe ,fedp, dfp , nel)
325 sn(i,1,1) = one + dfp(i,1,1)
326 sn(i,2,2) = one + dfp(i,2,2)
327 sn(i,3,3) = one + dfp(i,3,3)
328 sn(i,1,2) = dfp(i,1,2)
329 sn(i,2,3) = dfp(i,2,3)
330 sn(i,3,1) = dfp(i,3,1)
331 sn(i,2,1) = dfp(i,2,1)
332 sn(i,3,2) = dfp(i,3,2)
333 sn(i,1,3) = dfp(i,1,3)
342 1 nel , matb , c10, c01, c20,
343 2 c11 ,c02, c30, c21, c12,
344 3 c03 ,d1 ,d2 , d3, sigb ,
345 4 bi1,bi2,jdet ,flag_mul,
346 5 nvarf,coefr, betaf,coefm ,
347 6 uvarf,cbii,rbulk,iform)
350 sigb(i,1,1) = sb * sigb(i,1,1)
351 sigb(i,2,2) = sb * sigb(i,2,2)
352 sigb(i,3,3) = sb * sigb(i,3,3)
353 sigb(i,1,2) = sb * sigb(i,1,2)
354 sigb(i,2,3) = sb * sigb(i,2,3)
355 sigb(i,3,1) = sb * sigb(i,3,1)
359 . r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
363 uvar(i,1) = fp(i,1,1)
364 uvar(i,2) = fp(i,2,2)
365 uvar(i,3) = fp(i,3,3)
366 uvar(i,4) = fp(i,1,2)
367 uvar(i,5) = fp(i,2,3)
368 uvar(i,6) = fp(i,3,1)
369 uvar(i,7) = fp(i,2,1)
370 uvar(i,8) = fp(i,3,2)
371 uvar(i,9) = fp(i,1,3)
376 sigb(i,1:3,1:3) = sb * siga(i,1:3,1:3)
377 cbii(i,1:3) = caii(i,1:3)
385 signxx(i) = siga(i,1,1) + sigb(i,1,1)
386 signyy(i) = siga(i,2,2) + sigb(i,2,2)
387 signzz(i) = siga(i,3,3) + sigb(i,3,3)
388 signxy(i) = siga(i,1,2) + sigb(i,1,2)
389 signyz(i) = siga(i,2,3) + sigb(i,2,3)
390 signzx(i) = siga(i,3,1) + sigb(i,3,1)
395 cii(1:3)=caii(i,1:3)+sb*cbii(i,1:3)
396 c_max(i)=
max(stiff0,cii
414 soundsp(i)=sqrt(c_max(i)/rho(i))
417 IF (impl_s > 0 .OR. ihet > 1)
THEN