45
46
47
48 USE eos_param_mod , ONLY : eos_param_
49 USE constant_mod ,
ONLY : zero, em15, half, one, three_half, two, three, three100
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73 IMPLICIT NONE
74#include "my_real.inc"
75
76
77
78 INTEGER,INTENT(IN) :: IFLAG, NEL
79 my_real,
INTENT(IN) :: pmin,off(nel),mu(nel),espe(nel),dvol(nel),df(nel),vnew(nel)
80 my_real,
INTENT(INOUT) :: pnew(nel) ,dpdm(nel),dpde(nel),eint(nel),psh(nel)
81 TYPE(EOS_PARAM_),INTENT(IN) :: EOS_STRUCT
82
83
84
85 INTEGER I
86 my_real :: p0,gamma,e0,aa,bb,dvv,pp,ar0b,b
87
88
89
90 e0 = eos_struct%E0
91 psh(1:nel) = eos_struct%PSH
92
93 gamma = eos_struct%UPARAM(1)
94 p0 = eos_struct%UPARAM(2)
95 ar0b = eos_struct%UPARAM(3)
96 b = eos_struct%UPARAM(4)
97
98 IF(iflag == 0) THEN
99 DO i=1,nel
100 pp = (gamma-one)*(one+mu(i))*espe(i) + ar0b*exp(b*log(one+mu(i)))
101 dpdm(i) = (gamma-one)*espe(i)+ar0b*b/(one+mu(i))*exp((b-one)*log(one+mu(i)))+(gamma-one)*df(i)*pp
102 dpde(i) = (gamma-one)*(one+mu(i))
103 pnew(i) =
max(pp,pmin)*off(i)
104 pnew(i) = pnew(i) - psh(i)
105 ENDDO
106
107 ELSEIF(iflag == 1) THEN
108 DO i=1,nel
109 bb = (gamma-one)*(one+mu(i))
110 dpde(i) = bb
111 aa = ar0b*exp(b*log(one+mu(i)))
112 dvv = half*dvol(i)*df(i) /
max(em15,vnew(i))
113 pnew(i) = (aa + bb*espe(i) ) / (one+bb*dvv)
114 pnew(i) =
max(pnew(i),pmin)*off(i)
115 eint(i) = eint(i) - half*dvol(i)*(pnew(i))
116 pnew(i) = pnew(i) - psh(i)
117 ENDDO
118
119 ELSEIF (iflag == 2) THEN
120 DO i=1, nel
121 IF (vnew(i) > zero) THEN
122 pnew(i) = (gamma-one)*(one+mu(i))*espe(i) + ar0b*exp(b*log(one+mu(i)))
123 dpdm(i) = (gamma-one)*espe(i)+ar0b*b/(one+mu(i))*exp((b-one)*log(one+mu(i)))+(gamma-one)*df(i)*pnew(i)
124 dpde(i) = (gamma-one)*(one+mu(i))
125 pnew(i) = pnew(i)-psh(i)
126 ENDIF
127 ENDDO
128
129 ENDIF
130
131
132 RETURN