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
58
59
60#include "implicit_f.inc"
61#include "comlock.inc"
62
63
64
65#include "param_c.inc"
66#include "tabsiz_c.inc"
67#include "com04_c.inc"
68#include "com06_c.inc"
69#include "com08_c.inc"
70#include "vect01_c.inc"
71#include "scr06_c.inc"
72
73
74
75 INTEGER, INTENT(IN) :: MAT(NEL), IFLAG, NEL,NPF(SNPC)
76 my_real,
INTENT(INOUT) :: pm(npropm,nummat),
77 . off(nel) ,eint(nel) ,mu(nel) ,
78 . espe(nel) ,dvol(nel) ,df(nel) ,
79 . vnew(nel) ,pnew(nel) ,dpdm(nel),
80 . dpde(nel) ,
81 . psh(nel), tf(stf)
82
83
84
85 INTEGER I, MX
87 my_real :: xscale_a,xscale_b,fscale_a,fscale_b
88 INTEGER :: A_fun_id, B_fun_id
89 my_real :: res_a(nel),res_b(nel),deri_a(nel),deri_b(nel),pc
91
92
93
94 IF(iflag == 0) THEN
95 mx = mat(1)
96 e0 = pm(23,mx)
97 psh(1:nel) = pm(88,mx)
98 xscale_a = pm(33,mx)
99 xscale_b = pm(34,mx)
100 a_fun_id = pm(35,mx)
101 b_fun_id = pm(36,mx)
102 pc = pm( 37,mx)
103 fscale_a = pm(160,mx)
104 fscale_b = pm(161,mx)
105
106
107 IF(a_fun_id == 0)THEN
108 DO i=1,nel
109 res_a(i) = zero
110 deri_a(i) = zero
111 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
112 ENDDO
113 ELSEIF(b_fun_id == 0)THEN
114 DO i=1,nel
115 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
116 res_b(i) = zero
117 deri_b(i) = zero
118 ENDDO
119 ELSE
120 DO i=1,nel
121 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
122 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
123 ENDDO
124 ENDIF
125
126 DO i=1,nel
127 pp = res_a(i) + res_b(i) * espe(i) - psh(i)
128 dpdm(i) = deri_a(i)+deri_b(i)*espe(i) + res_b(i)*(pp+psh(i))/( (one+mu(i))*(one+mu(i)) )
129 dpde(i) = res_b(i)
130 pnew(i) =
max(pp,pc)*off(i)
131 ENDDO
132
133
134 ELSEIF(iflag == 1) THEN
135 mx = mat(1)
136 e0 = pm(23,mx)
137 psh(1:nel) = pm(88,mx)
138 xscale_a = pm(33,mx)
139 xscale_b = pm(34,mx)
140 a_fun_id = pm(35,mx)
141 b_fun_id = pm(36,mx)
142 pc = pm( 37,mx)
143 fscale_a = pm(160,mx)
144 fscale_b = pm(161,mx)
145
146 IF(a_fun_id == 0)THEN
147 DO i=1,nel
148 res_a(i) = zero
149 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
150 ENDDO
151 ELSEIF(b_fun_id == 0)THEN
152 DO i=1,nel
153 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
154 res_b(i) = zero
155 ENDDO
156 ELSE
157 DO i=1,nel
158 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
159 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
160 ENDDO
161 ENDIF
162 DO i=1,nel
163 aa = res_a(i)
164 bb = res_b(i)
165 dpde(i) = bb
166 dvv = half*dvol(i)*df(i) /
max(em15,vnew(i))
167 pp = aa + bb * espe(i)
168 pnew(i) = (aa+bb*(espe(i)-psh(i)*dvv))/(one+bb*dvv)
169 pnew(i) =
max(pnew(i),pc )*off(i)
170 eint(i) = eint(i) - half*dvol(i)*(pnew(i)+psh(i) )
171 ENDDO
172
173
174 ELSEIF (iflag == 2) THEN
175 mx = mat(1)
176 e0 = pm(23,mx)
177 psh(1:nel) = pm(88,mx)
178 xscale_a = pm(33,mx)
179 xscale_b = pm(34,mx)
180 a_fun_id = pm(35,mx)
181 b_fun_id = pm(36,mx)
182 pc = pm( 37,mx)
183 fscale_a = pm(160,mx)
184 fscale_b = pm(161,mx)
185
186 IF(a_fun_id == 0)THEN
187 DO i=1,nel
188 res_a(i) = zero
189 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
190 deri_a(i) = zero
191 ENDDO
192 ELSEIF(b_fun_id == 0)THEN
193 DO i=1,nel
194 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
195 res_b(i) = zero
196 deri_b(i) = zero
197 ENDDO
198 ELSE
199 DO i=1,nel
200 res_a(i) = fscale_a*finter(a_fun_id,mu(i),npf,tf,deri_a(i))
201 res_b(i) = fscale_b*finter(b_fun_id,mu(i),npf,tf,deri_b(i))
202 ENDDO
203 ENDIF
204 DO i=1, nel
205 IF (vnew(i) > zero) THEN
206 pp = res_a(i) + res_b(i)*espe(i) - psh(i)
207 dpdm(i) = deri_a(i)+deri_b(i)*espe(i) + res_b(i)*(pp+psh(i))/( (one+mu(i))*(one+mu(i)) )
208 dpde(i) = res_b(i)
209 pnew(i) = pp
210 ENDIF
211 ENDDO
212
213 ENDIF
214
215 RETURN