44
45
46
47#include "implicit_f.inc"
48#include "comlock.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "param_c.inc"
59#include "vect01_c.inc"
60
61
62
63 INTEGER L_TEMP, L_PLAS, L_BFRAC,L_BULK,NEL
64 INTEGER MAT(*), NC1(*), NC2(*), NC3(*), NC4(*)
66 . pm(npropm,*), alph1(*), alph2(*), volt(*), fill(numnod,*),
67 . sig1(nel,6), eint1(*), volo1(*), rhon1(*), flux1(4,*), flu11(*),
68 . off1(*), sig2(nel,6), eint2(*), volo2(*), rhon2(*), flux2(4,*),
69 . flu12(*), off2(*), sigt(nel,6), eintt(*), rhot(*), tempt(*),
70 . plast(*), bfract(*),voln(mvsiz),bulkt(*),
71 . aire(*), aires(*),
72 . d1(*), d2(*), d3(*), d4(*), d5(*), d6(*),
73 . d1s(*), d2s(*), d3s(*), d4s(*), d5s(*), d6s(*),
74 . dalph1(*), dalph2(*)
75
76
77
78 INTEGER I, MX1, MX2
80 . alph1n(mvsiz), alph2n(mvsiz), alpn1,
81 . alpn2, alpn3, alpn4, alpd1, alpd2, alpd3, alpd4, alpd, alpn,
82 . eps, soff1, soff2, alphn, alphnn, c11, c12, alpht, da,
83 . ialph1, ialph2
84
85 DO i=1,nel
86 sigt(i,1)=zero
87 sigt(i,2)=zero
88 sigt(i,3)=zero
89 sigt(i,4)=zero
90 sigt(i,5)=zero
91 sigt(i,6)=zero
92 eintt(i) =zero
93 rhot(i) =zero
94 volt(i) =voln(i)
95 aires(i) =aire(i)
96 d1s(i) =d1(i)
97 d2s(i) =d2(i)
98 d3s(i) =d3(i)
99 d4s(i) =d4(i)
100 d5s(i) =d5(i)
101 d6s(i) =d6(i)
102 ENDDO
103
104 IF(l_temp>0)THEN
105 DO i=1,nel
106 tempt(i)=zero
107 ENDDO
108 ENDIF
109
110 IF(l_plas>0)THEN
111 DO i=1,nel
112 plast(i)=zero
113 ENDDO
114 ENDIF
115
116 IF(l_bfrac>0)THEN
117 DO i=1,nel
118 bfract(i)=zero
119 ENDDO
120 ENDIF
121
122 IF(l_bulk>0)THEN
123 DO i=1,nel
124 bulkt(i)=zero
125 ENDDO
126 ENDIF
127
128
129
130 DO i=1,nel
131 alph1n(i)=zero
132 alph2n(i)=zero
133 alpn1=
max(zero,fill(nc1(i),1))
134 alpn2=
max(zero,fill(nc2(i),1))
135 alpn3=
max(zero,fill(nc3(i),1))
136 alpn4=
max(zero,fill(nc4(i),1))
137 alpd1=abs(fill(nc1(i),1))
138 alpd2=abs(fill(nc2(i),1))
139 alpd3=abs(fill(nc3(i),1))
140 alpd4=abs(fill(nc4(i),1))
141 alpd=alpd1+alpd2+alpd3+alpd4
142 alpn=alpn1+alpn2+alpn3+alpn4
143 IF (alpd>em20) alph1n(i)=alpn/alpd
144 ENDDO
145
146 IF(jmult>1)THEN
147 DO i=1,nel
148 alpn1=
max(zero,fill(nc1(i),2))
149 alpn2=
max(zero,fill(nc2(i),2))
150 alpn3=
max(zero,fill(nc3(i),2))
151 alpn4=
max(zero,fill(nc4(i),2))
152 alpd1=abs(fill(nc1(i),2))
153 alpd2=abs(fill(nc2(i),2))
154 alpd3=abs(fill(nc3(i),2))
155 alpd4=abs(fill(nc4(i),2))
156 alpd=alpd1+alpd2+alpd3+alpd4
157 alpn=alpn1+alpn2+alpn3+alpn4
158 IF(alpd>em20)alph2n(i)=alpn/alpd
159 ENDDO
160 ENDIF
161
162 eps=em15
163 soff1=zero
164 soff2=zero
165
166 IF(jmult==1)THEN
167
168 DO i=1,nel
169 dalph1(i)=zero
170 IF(alph1n(i)/=zero.AND.alph1n(i)/=one)THEN
171 alphn=(volo1(i)-dt1
172 . +flux1(1,i)+flux1(2,i)+flux1
173 alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1
174 IF((sig1(i,1)+sig1(i,2)+sig1(i,3))>zero)THEN
175 alphn=
min(alphn,alphnn)
176 ELSE
177 alphn=
max(alphn,alphnn)
178 ENDIF
179 IF(alphn<=zep99.AND.alphn>zero.AND.rhon1(i)<=zero)THEN
180 IF(rhon1(i)/=zero)THEN
181#include "lockon.inc"
182 WRITE(6,*)' ***NEGATIVE RHO****ALPH1,RHON1'
183 WRITE(6,*)i+nft,alph1n(i),rhon1(i)
184#include "lockoff.inc"
185 rhon1(i)=zero
186 ENDIF
187 alphn=zero
188 ENDIF
189 dalph1(i)=alphn-alph1n(i)
190 alph1(i)=
min(one,alphn)
191 alph1(i)=
max(zero,alph1(i))
192 ELSEIF(alph1n(i)==one.AND.alph1(i)<zep999)THEN
193 alphn=(volo1(i)-dt1*half*(flu11(i)
194 . +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
195 alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
196 alphn=
max(alphn,alphnn,alph1(i))
197 IF(rhon1(i)<=zero)THEN
198 rhon1(i)=zero
199 alphn=zero
200 ENDIF
201 alph1(i)=
min(one,alphn)
202 ELSEIF(alph1n(i)==zero.AND.alph1(i)>em3)THEN
203 alph1(i)=alph1n(i)
204 ELSE
205 alph1(i)=alph1n(i)
206 ENDIF
207 ENDDO
208
209 DO i=1,nel
210 IF(alph1(i)<eps)THEN
211 alph1(i)=zero
212 off1(i) =zero
213 volo1(i)=zero
214 eint1(i)=zero
215 ELSE
216 off1(i)=one
217 ENDIF
218 soff1=soff1+off1(i)
219 IF(alph1(i)==zero)THEN
220 dalph1(i)=zero
221 ialph1=1
222 ELSEIF(alph1(i)==one)THEN
223 dalph1(i)=zero
224 ialph1=1
225 ENDIF
226 ENDDO
227
228 IF(ialph1==1) THEN
230 DO i=1,nel
231 IF(alph1(i)==zero)THEN
232 fill(nc1(i),1)=
min(zero,fill(nc1(i),1))
233 fill(nc2(i),1)=
min(zero,fill(nc2(i),1))
234 fill(nc3(i),1)=
min(zero,fill(nc3(i),1))
235 fill(nc4(i),1)=
min(zero,fill(nc4(i),1))
236 ELSEIF(alph1(i)==one)THEN
237 fill(nc1(i),1)=
max(zero,fill(nc1(i),1))
238 fill(nc2
239 fill(nc3(i),1)=
max(zero,fill(nc3(i),1))
240 fill(nc4(i),1)=
max(zero,fill(nc4(i),1))
241 ENDIF
242 ENDDO
244 ENDIF
245
246 ELSEIF(jmult==2)THEN
247 ialph1 = 0
248 ialph2 = 0
249 DO i=1,nel
250 dalph1(i)=zero
251 dalph2(i)=zero
252 mx1=nint(pm(21,mat(i)))
253 c11 =pm(32,mx1)
254 mx2=nint(pm(22,mat(i)))
255 c12 =pm(32,mx2)
256 IF(alph1n(i)/=zero.AND.alph1n(i)/=one)THEN
257 alphn=(volo1(i)-dt1*half*(flu11(i)
258 . +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
259 alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
260 IF((sig1(i,1)+sig1(i,2)+sig1(i,3))>zero)THEN
261 alphn=
min(alphn,alphnn)
262 ELSE
263 alphn=
max(alphn,alphnn)
264 ENDIF
265 IF(alphn<=zep99.AND.alphn>zero.AND.rhon1(i)<=zero)THEN
266 IF(rhon1(i)/=zero)THEN
267#include "lockon.inc"
268 WRITE(6,*)' ***NEGATIVE RHO****ALPH1,RHON1'
269 WRITE(6,*)i+nft,alph1n(i),rhon1(i)
270#include "lockoff.inc"
271 rhon1(i)=zero
272 ENDIF
273 alphn=zero
274 ENDIF
275 dalph1(i)=alphn-alph1n(i)
276 alph1(i)=
min(one,alphn)
277 alph1(i)=
max(zero,alph1(i))
278 ELSEIF(alph1n(i)==one.AND.alph1(i)<zep999)THEN
279 alphn=(volo1(i)-dt1*half*(flu11(i)
280 . +flux1(1,i)+flux1(2,i)+flux1(3,i)+flux1(4,i)))/voln(i)
281 alphnn=alphn*(one + (d1(i)+d2(i)+d3(i))*dt1)
282 alphn=
max(alphn,alphnn,alph1(i))
283 IF(rhon1(i)<=zero)THEN
284 rhon1(i)=zero
285 alphn=zero
286 ENDIF
287 alph1(i)=
min(one,alphn)
288 ELSEIF(alph1n(i)==zero.AND.alph1(i)>em3)THEN
289 alph1(i)=alph1n(i)
290 ELSE
291 alph1(i)=alph1n(i)
292 ENDIF
293 IF(alph2n(i)/=zero.AND.alph2n(i)/=one)THEN
294 alphn=(volo2(i)-dt1*half*(flu12(i)
295 . +flux2(1,i)+flux2(2,i)+flux2(3,i)+flux2(4,i)))/voln(i)
296 alphnn=alphn*(one+(d1(i)+d2(i)+d3(i))*dt1)
297 IF((sig2(i,1)+sig2(i,2)+sig2(i,3))>one)THEN
298 alphn=
min(alphn,alphnn)
299 ELSE
300 alphn=
max(alphn,alphnn)
301 ENDIF
302 IF(alphn<=zep99.AND.alphn>zero.AND.rhon2(i)<=zero)THEN
303 IF(rhon2(i)/=zero)THEN
304#include "lockon.inc"
305 WRITE(6,*)' ***NEGATIVE RHO****ALPH2,RHON2'
306 WRITE(6,*)i+nft,alph2n(i),rhon2(i)
307#include "lockoff.inc"
308 rhon2(i)=zero
309 ENDIF
310 alphn=zero
311 ENDIF
312 dalph2(i)=alphn-alph2n(i)
313 alph2(i)=
min(one,alphn)
314 alph2(i)=
max(zero,alph2(i))
315 ELSEIF(alph2n(i)==one.AND.alph2(i)<=zep999)THEN
316 alphn=(volo2(i)-dt1*half*(flu12(i)
317 . +flux2(1,i)+flux2(2,i)+flux2(3,i)+flux2(4,i)))/voln(i)
318 alphnn=alphn*(one + (d1(i)+d2(i)+d3(i))*dt1)
319 alphn=
max(alphn,alphnn,alph2(i))
320 IF(rhon2(i)<=zero)THEN
321 rhon2(i)=zero
322 alphn=zero
323 ENDIF
324 alph2(i)=
min(one,alphn)
325 ELSEIF(alph2n(i)==zero.AND.alph2(i)>em3)THEN
326 alph2(i)=alph2n(i)
327 ELSE
328 alph2(i)=alph2n(i)
329 ENDIF
330 alpht=alph1(i)+alph2(i)
331 IF(alpht>one)THEN
332 da=(alpht-one)/(c11*alph2(i)+c12*alph1(i))
333 alph1(i)=alph1(i)*(one -c12*da)
334 alph2(i)=alph2(i)*(one -c11*da)
335 dalph1(i)=dalph1(i)-c12*da
336 dalph2(i)=dalph2(i)-c11*da
337 ENDIF
338 ENDDO
339
340 DO i=1,nel
341 IF(alph1(i)<eps)THEN
342 alph1(i)=zero
343 off1(i) =zero
344 volo1(i)=zero
345 eint1(i)=zero
346 ELSE
347 off1(i) =one
348 ENDIF
349 soff1=soff1+off1(i)
350 IF(alph2(i)<eps)THEN
351 alph2(i)=zero
352 off2(i) =zero
353 volo2(i)=zero
354 eint2(i)=zero
355 ELSE
356 off2(i)=one
357 ENDIF
358 soff2=soff2+off2(i)
359 IF(alph1(i)==zero)THEN
360 dalph1(i)=zero
361 ialph1 = 1
362 ELSEIF(alph1(i)==one)THEN
363 dalph1(i)=zero
364 ialph1 = 1
365 ENDIF
366 IF(alph2(i)==zero)THEN
367 dalph2(i)=zero
368 ialph2 = 1
369 ELSEIF(alph2(i)==one)THEN
370 dalph2(i)=zero
371 ialph2 = 1
372 ENDIF
373 ENDDO
374
375 IF(ialph1==1) THEN
377 DO i=1,nel
378 IF(alph1(i)==zero)THEN
379 fill(nc1(i),1)=
min(zero,fill(nc1(i),1))
380 fill(nc2(i),1)=
min(zero,fill(nc2(i),1))
381 fill(nc3(i),1)=
min(zero,fill(nc3(i),1))
382 fill(nc4(i),1)=
min(zero,fill(nc4(i),1))
383 ELSEIF(alph1(i)==one)THEN
384 fill(nc1(i),1)=
max(zero,fill(nc1(i),1))
385 fill(nc2(i),1)=
max(zero,fill(nc2(i),1))
386 fill(nc3(i),1)=
max(zero,fill(nc3(i),1))
387 fill(nc4(i),1)=
max(zero,fill(nc4(i),1))
388 ENDIF
389 ENDDO
391 ENDIF
392
393 IF(ialph2==1) THEN
395 DO i=1,nel
396 IF(alph2(i)==zero)THEN
397 fill(nc1(i),2)=
min(zero,fill(nc1(i),2))
398 fill(nc2(i),2)=
min(zero,fill(nc2(i),2))
399 fill(nc3(i),2)=
min(zero,fill(nc3(i),2))
400 fill(nc4(i),2)=
min(zero,fill(nc4(i),2))
401 ELSEIF(alph2(i)==one)THEN
402 fill(nc1(i),2)=
max(zero,fill(nc1(i),2))
403 fill(nc2(i),2)=
max(zero,fill(nc2(i),2))
404 fill(nc3(i),2)=
max(zero,fill(nc3(i),2))
405 fill(nc4(i),2)=
max(zero,fill(nc4(i),2))
406 ENDIF
407 ENDDO
409 ENDIF
410 ENDIF
411
412 RETURN