31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57#include "implicit_f.inc"
58#include "comlock.inc"
59
60
61
62#include "param_c.inc"
63#include "com04_c.inc"
64#include "com06_c.inc"
65#include "com08_c.inc"
66#include "vect01_c.inc"
67#include "scr06_c.inc"
68
69
70
71 INTEGER MAT(NEL), IFLAG, NEL
73 . off(nel) ,eint(nel) ,mu(nel) ,
74 . dvol(nel) ,
75 . vnew(nel) ,pnew(nel) ,dpdm(nel),
76 . dpde(nel)
77 my_real,
INTENT(INOUT) :: psh(nel)
78
79
80
81 INTEGER I, MX
83
84
85
86 IF(iflag == 0) THEN
87 mx = mat(1)
88 k0 = pm(36,mx)
89 p0 = pm(32,mx
90 psh(1:nel
91 k1 = pm
92 pc = pm(37,mx)
93 DO i=1,nel
94 dpdm(i) = k0*exp( k1*log(1+mu(i)) )/(one+mu(i))
95 dpde(i) = zero
96 pnew(i) = k0/k1 * (exp(k1*log(one+mu(i))) - one ) + p0
97 pnew(i) =
max(pnew(i),pc)*off(i)
98 pnew(i) = pnew(i)-psh(i)
99 ENDDO
100
101 ELSEIF(iflag == 1) THEN
102 mx = mat(1)
103 k0 = pm(36,mx)
104 p0 = pm(32,mx)
105 psh(1:nel) = pm(88,mx)
106 k1 = pm(35,mx)
107 pc = pm(37,mx)
108 DO i=1,nel
109 pnew(i) = k0/k1 * (exp(k1*log(one+mu(i))) - one ) + p0
110 pnew(i) =
max(pnew(i),pc)*off(i)
111 eint(i) = eint(i) - half*dvol(i)*(pnew(i)+psh(i))
112 pnew(i) = pnew(i)-psh(i)
113 dpde(i) = zero
114 ENDDO
115
116 ELSEIF(iflag == 2) THEN
117 mx = mat(1)
118 k0 = pm(36,mx)
119 p0 = pm(32,mx)
120 psh(1:nel) = pm(88,mx)
121 k1 = pm(35,mx)
122 pc = pm(37,mx)
123 DO i=1,nel
124 IF (vnew(i) > zero) THEN
125 pnew(i) = k0/k1 * (exp(k1*log(1+mu(i))) - one ) + p0
126 dpdm(i) = k0*exp( k1*log(1+mu(i)) )/(one+mu(i))
127 dpde(i) = zero
128 pnew(i) =
max(pnew(i),pc)*off(i)
129 pnew(i) = pnew(i)-psh(i)
130 ENDIF
131 ENDDO
132 ENDIF
133
134 RETURN