38
39
40
42
43
44
45#include "implicit_f.inc"
46#include "comlock.inc"
47
48
49
50#include "mvsiz_p.inc"
51
52
53
54#include "param_c.inc"
55#include "vect01_c.inc"
56#include "cong1_c.inc"
57#include "com04_c.inc"
58#include "com06_c.inc"
59#include "com08_c.inc"
60#include "scr14_c.inc"
61#include "scr16_c.inc"
62
63
64
65 my_real pm(npropm,nummat),geo(npropg,numgeo), off(*), rho(*),eani(*),
66 . y1(*),y2(*),y3(*),y4(*),z1(*),z2(*),z3(*),z4(*),
67 . vy1(*), vy2(*), vy3(*), vy4(*), vz1(*), vz2(*), vz3(*),
68 . py1(*), py2(*), pz1(*), pz2(*),ehou(*),
69 . t11(*), t12(*), t13(*), t14(*), t21(*), t22(*), t23(*), t24(*),
70 . vz4(*),
area(*), cxx(*),vd2(*),vis(*), partsav(npsav,*)
71 INTEGER MAT(*),PID(*),IPARTQ(NUMELQ), IPARG(63:63)
72
73
74
75 INTEGER I,MX, ISFLUID
77 . fcl(mvsiz) , fcq(mvsiz),
78 . g11(mvsiz) , g21(mvsiz) , g31(mvsiz) , g41(mvsiz),
79 . hgy(mvsiz), hgz(mvsiz),
80 . hy,hz,fac,px1h1,px2h1,ehourt, are
81
82
83
84 IF(mtn == 11 .OR. ((mtn == 17 .OR. mtn == 47) .AND.
ale%UPWIND%UPWM == 0))
THEN
85 DO i=lft,llt
86 t11(i) = zero
87 t12(i) = zero
88 t13(i) = zero
89 t14(i) = zero
90 t21(i) = zero
91 t22(i) = zero
92 t23(i) = zero
93 t24(i) = zero
94 ENDDO
95 RETURN
96 ENDIF
97
98
99
100 IF(invstr >= 35)THEN
101 DO i=lft,llt
102 caq(i)=geo(13,pid(i))
103 ENDDO
104 ELSE
105 DO i=lft,llt
106 caq(i)=pm(4,mat(i))
107 ENDDO
108 ENDIF
109
110 DO i=lft,llt
112 fcq(i)=rho(i)*sqrt(are)
113 fcl(i)=caq(i)*fcq(i)
114 ENDDO
115
116 isfluid=iparg(63)
117
118 IF(isfluid == 1 .AND.
ale%UPWIND%UPWM == 0)
THEN
119 DO i=lft,llt
120 fcl(i)=fcl(i)*cxx(i)
121 fcq(i)=zero
122 ENDDO
123 ELSEIF(isfluid == 1 .AND.
ale%UPWIND%UPWM == 1)
THEN
124 DO i=lft,llt
125 fcl(i)=
min(fcl(i)*cxx(i),
max(20.*caq(i)*vis(i),fcl(i)*sqrt(vd2(i))))
126 fcq(i)=zero
127 ENDDO
128 ELSEIF(isfluid == 1 .AND.
ale%UPWIND%UPWM > 0)
THEN
129 DO i=lft,llt
130 IF(vis(i) > zero)THEN
131 fcq(i)=zero
132 fcl(i)=twenty*caq(i)*vis(i)
133 ELSE
134 fcq(i)=fcl(i)*caq(i)*hundred
135 fcl(i)=fcl(i)*cxx(i)
136 ENDIF
137 ENDDO
138 ELSE
139 DO i=lft,llt
140 fcq(i)=fcl(i)*caq(i)*hundred
141 fcl(i)=fcl(i)*cxx(i)
142 ENDDO
143 ENDIF
144 IF(impl /= zero)THEN
145 DO i=lft,llt
146 fcq(i)=zero
147 ENDDO
148 ENDIF
149
150 DO i=lft,llt
151 IF(off(i) < one)THEN
152 fcl(i)=zero
153 fcq(i)=zero
154 ENDIF
155 ENDDO
156
157 IF(jhbe == 0)THEN
158
159
160
161 DO i=lft,llt
162 hgy(i)=half*(vy1(i)-vy2(i)+vy3(i)-vy4(i))
163 hgz(i)=half*(vz1(i)-vz2(i)+vz3(i)-vz4(i))
164 ENDDO
165 DO i=lft,llt
166 t11(i)=hgy(i)*(fcl(i)+abs(hgy(i))*fcq(i))
167 t12(i)=-t11(i)
168 t13(i)= t11(i)
169 t14(i)=-t11(i)
170 t21(i)=hgz(i)*(fcl(i)+abs(hgz(i))*fcq(i))
171 t22(i)=-t21(i)
172 t23(i)= t21(i)
173 t24(i)=-t21(i)
174 ehou(i)= two*dt1*(t11(i)*hgy(i) + t21(i)*hgz(i))
175 ENDDO
176 ELSE
177
178
179
180 DO i=lft,llt
181 hy=y1(i)-y2(i)+y3(i)-y4(i)
182 hz=z1(i)-z2(i)+z3(i)-z4(i)
184 px1h1=fac*(py1(i)*hy+pz1(i)*hz)
185 px2h1=fac*(py2(i)*hy+pz2(i)*hz)
186 g11(i)= one -px1h1
187 g21(i)=-one -px2h1
188 g31(i)= one +px1h1
189 g41(i)=-one +px2h1
190 ENDDO
191 DO i=lft,llt
192 hgy(i)=half*(g11(i)*vy1(i)+g21(i)*vy2(i)+g31(i)*vy3(i)+g41(i)*vy4(i))
193 hgz(i)=half*(g11(i)*vz1(i)+g21(i)*vz2(i)+g31(i)*vz3(i)+g41(i)*vz4(i))
194 ENDDO
195 DO i=lft,llt
196 hy=hgy(i)*(fcl(i)+abs(hgy(i))*fcq(i))
197 hz=hgz(i)*(fcl(i)+abs(hgz(i))*fcq(i))
198 t11(i) =g11(i)*hy
199 t12(i) =g21(i)*hy
200 t13(i) =g31(i)*hy
201 t14(i) =g41(i)*hy
202 t21(i) =g11(i)*hz
203 t22(i) =g21(i)*hz
204 t23(i) =g31(i)*hz
205 t24(i) =g41(i)*hz
206 ehou(i)= two*dt1*(hy*hgy(i) + hz*hgz(i))
207 ENDDO
208 ENDIF
209
210 IF(jlag == 1)THEN
211 ehourt = zero
212 DO i=lft,llt
213 ehourt= ehourt+ehou(i)
214 ENDDO
215 DO i=lft,llt
216 mx = ipartq(i)
217 partsav(8,mx)=partsav(8,mx) + ehou(i)
218 ENDDO
219
220 ehour = ehour + ehourt
221 ENDIF
222
223
224 DO i=lft,llt
225 eani(nft+i) = eani(nft+i)+ehou(i)/
max(em30,rho(i)*
area(i))
226 ENDDO
227
228
229 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)