41
42
43
44 USE eos_param_mod , ONLY : eos_param_
45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "com01_c.inc"
57#include "com08_c.inc"
58#include "param_c.inc"
59#include "scr14_c.inc"
60
61
62
63 INTEGER, INTENT(IN) :: JPOR
64 INTEGER MAT(*),NEL
66 . pm(npropm,*), off(*), sig(nel,6), eint(*), rho(*), rk(*), re(*),
67 . vorti(*),wxx(*),wyy(*),wzz(*),voln(mvsiz),vis(*),
68 . d1(*), d2(*), d3(*), d4(*), d5(*), d6(*), alogey(*), ssp(*),
69 . rho0(*), tmu(*), amu(*), amu2(*), psh(*), pc(*), espe(*),
70 . c1(*), c2(*), c3(*), c4(*), c5(*), c6(*), df(*), dpdm(*)
71 TYPE(EOS_PARAM_),INTENT(INOUT) :: EOS_STRUCT
72
73
74
75 INTEGER I, MX
77 . vis2(mvsiz),
78 . dav(mvsiz),
79 . yp0, cmu, ax, e, a, xmu,
80 . xm, xk, xe, yplus, rk2t, fac,
81 . rho0_1, vis_1, pc_1, c1_1, c2_1,
82 . c3_1, c4_1, c5_1, c6_1, psh_1
83
84
85
86 IF((anim_e(10)==1 .OR. anim_se(10)==1).AND. dt1/=0.)THEN
87
88 fac=four/dt1
89 IF(n2d==0)THEN
90 DO 5 i=1,nel
91 5 vorti(i)=fac*sqrt(wxx(i)**2+wyy(i)**2+wzz(i)**2)
92 ELSE
93 DO 6 i=1,nel
94 6 vorti(i)=fac*wzz(i)
95 ENDIF
96 ENDIF
97
98 mx =mat(1)
99
100 rho0_1=pm( 1,mx)
101 vis_1 =pm(24,mx)
102 pc_1 =pm(37,mx)
103 c1_1 = eos_struct%UPARAM(1)
104 c2_1 = eos_struct%UPARAM(2)
105 c3_1 = eos_struct%UPARAM(3)
106 c4_1 = eos_struct%UPARAM(4)
107 c5_1 = eos_struct%UPARAM(5)
108 c6_1 = eos_struct%UPARAM(6)
109 psh_1 =pm(88,mx)
110
111 DO 10 i=1,nel
112 rho0(i)=rho0_1
113 vis(i) =vis_1
114 pc(i) =pc_1
115 c1(i) =c1_1
116 c2(i) =c2_1
117 c3(i) =c3_1
118 c4(i) =c4_1
119 c5(i) =c5_1
120 c6(i) =c6_1
121 psh(i) =psh_1
122 10 CONTINUE
123
124 mx =mat(1)
125 IF(jpor/=2)THEN
126 DO 20 i=1,nel
127 yp0=pm(51,mx)
128 cmu=pm(81,mx)
129 ax =pm(47,mx)
130 e =pm(48,mx)
131 a =pm(49,mx)
132 xmu =rho(i)*vis(i)
133 xm =rho(i)*voln(i)
134 xk =rk(i)/xm
135 xe =
max(em15,re(i)/xm)
136 yplus =cmu*xk**2/(ax*xe*vis(i))
137 yplus =
max(yplus,yp0)
138 alogey(i)= a * log(e*yplus)
139 vis(i) =xmu*ax*yplus/alogey(i)
140 tmu(i) =vis(i)-xmu
141 20 CONTINUE
142 ELSE
143 DO 21 i=1,nel
144 vis(i)=zero
145 tmu(i)=zero
146 21 CONTINUE
147 ENDIF
148
149 DO 30 i=1,nel
150 vis2(i)=two*vis(i)
151 30 dav(i) =-(d1(i)+d2(i)+d3(i))*third
152
153
154
155 DO 40 i=1,nel
156 sig(i,1)=vis2(i)*(d1(i)+dav(i))
157 sig(i,2)=vis2(i)*(d2(i)+dav(i))
158 sig(i,3)=vis2(i)*(d3(i)+dav(i))
159 sig(i,4)=vis(i) *d4(i)
160 sig(i,5)=vis(i) *d5(i)
161 40 sig(i,6)=vis(i) *d6(i)
162
163 DO 50 i=1,nel
164 df(i) =rho0(i)/rho(i)
165 amu(i) =one/df(i)-one
166 amu2(i)=
max(zero,amu(i))**2
167 50 espe(i)=df(i)*eint(i)/voln(i)
168
169 DO 60 i=1,nel
170 rk2t=two*rk(i)/(three*voln(i))
171 dpdm(i) = dpdm(i)
172 . +(c5(i)+c6(i)*amu(i))*df(i)*df(i)*rk2t + rk2t*df(i)
173 60 CONTINUE
174
175 DO 70 i=1,nel
176 70 ssp(i)=sqrt(abs(dpdm(i))/rho0(i))
177
178 RETURN