45
46
47
48 USE eos_param_mod , ONLY : eos_param_
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72 IMPLICIT NONE
73#include "my_real.inc"
74
75
76
77 INTEGER,INTENT(IN) :: IFLAG, NEL
78 my_real,
INTENT(IN) :: off(nel) ,mu(nel) ,dvol(nel) ,df(nel) ,vnew(nel), rho0(nel), pmin
79 my_real,
INTENT(INOUT) :: pnew(nel), dpdm(nel), dpde(nel),psh(nel),eint(nel),espe(nel)
80 TYPE(EOS_PARAM_),INTENT(IN):: EOS_STRUCT
81
82
83
84 INTEGER I
86 my_real :: a1,a2,b0,b1,b2,c0,c1,d0, a2_
88
89
90
91 e0 = eos_struct%E0
92 psh(1:nel) = eos_struct%PSH
93
94 a1 = eos_struct%UPARAM(1)
95 a2 = eos_struct%UPARAM(2)
96 b0 = eos_struct%UPARAM(3)
97 b1 = eos_struct%UPARAM(4)
98 b2 = eos_struct%UPARAM(5)
99 c0 = eos_struct%UPARAM(6)
100 c1 = eos_struct%UPARAM(7)
101 d0 = eos_struct%UPARAM(8)
102 p0 = eos_struct%UPARAM(9)
103
104 IF(iflag == 0) THEN
105 DO i=1,nel
106 a2_=a2
107 IF(mu(i)<zero)a2_=-a2
108 denom = (espe(i)+d0)
109 pp = (a1*mu(i)+a2_*mu(i)*mu(i)+(b0+b1*mu(i)+b2*mu(i)*mu(i))*espe(i)+(c1*mu(i)+c0)*espe(i)*espe(i))/denom
110 pp =
max(pp,pmin) * off(i)
111 dpdmu = (a1+2*a2_*mu(i)+(two*b2*mu(i)+b1)*espe(i)+c1*espe(i)*espe(i))/denom
112 dpde(i) = (((b2*mu(i)+b1)*mu(i)+b0)+(two*(c1*mu(i)+c0))*espe(i) - pp/denom)/denom
113 dpdm(i) = dpdmu + dpde(i)*df(i)*df(i)*(pp)
114 pnew(i) = pp
115 pnew(i) = pnew(i) - psh(i)
116 ENDDO
117
118 ELSEIF(iflag == 1) THEN
119 DO i=1,nel
120 a2_=a2
121 IF(mu(i)<zero)a2_=-a2
122 dvv = dvol(i)*df(i) /
max(em15,vnew(i))
123 dvv = half*dvv
124 denom = (espe(i)+d0)
125 pp = (a1*mu(i)+a2_*mu(i)*mu(i)+(b0+b1*mu(i)+b2*mu(i)*mu(i))*espe(i)+(c1*mu(i)+c0)*espe(i)*espe(i))/denom
126 pp =
max(pp,pmin) * off(i)
127 espe(i) = espe(i) - (pp)*dvv
128 denom = (espe(i)+d0)
129 pp = (a1*mu(i)+a2_*mu(i)*mu(i)+(b0+b1*mu(i)+b2*mu(i)*mu(i))*espe(i)+(c1*mu(i)+c0)*espe(i)*espe(i))/denom
130 pp =
max(pp,pmin) * off(i)
131 espe(i) = espe(i) - (pp)*dvv
132 pnew(i) = pp
133 eint(i) = eint(i) - half*dvol(i)*(pnew(i))
134 pnew(i) = pnew(i) - psh(i)
135 dpde(i) = (((b2*mu(i)+b1)*mu(i)+b0)+(two*(c1*mu(i)+c0))*espe(i) - pp/denom)/denom
136 ENDDO
137
138 ELSEIF(iflag == 2) THEN
139 DO i=1, nel
140 IF (vnew(i) > zero) THEN
141 a2_=a2
142 IF(mu(i)<zero)a2_=-a2
143 denom = (espe(i)+d0)
144 pp = (a1*mu(i)+a2_*mu(i)*mu(i)+(b0+b1*mu(i)+b2*mu(i)*mu(i))*espe(i)+(c1*mu(i)+c0)*espe(i)*espe(i))/denom
145 pp =
max(pp,pmin)*off(i)
146 dpdmu = (a1+2*a2_*mu(i)+(two*b2*mu(i)+b1)*espe(i)+c1*espe(i)*espe(i))/denom
147 dpde(i) = (((b2*mu(i)+b1)*mu(i)+b0)+(two*(c1*mu(i)+c0))*espe(i) - pp/denom)/denom
148 dpdm(i) = dpdmu + dpde(i)*df(i)*df(i)*(pp)
149 pnew(i) = pp
150 pnew(i) = pnew(i) - psh(i)
151 ENDIF
152 ENDDO
153
154 ENDIF
155
156
157 RETURN