45
46
47
48 USE constant_mod ,
ONLY : zero, em15, half, one, two, three, three100
49 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) :: pmin, off(nel) ,mu(nel), dvol(nel) ,vnew(nel)
79 my_real,
INTENT(INOUT) :: psh(nel),pnew(nel) ,dpdm(nel),dpde(nel),eint(nel)
80 TYPE(EOS_PARAM_),INTENT(IN) :: EOS_STRUCT
81
82
83
84 INTEGER I
86
87
88
89 k0 = eos_struct%UPARAM(1)
90 p0 = eos_struct%UPARAM(2)
91 psh(1:nel) = eos_struct%PSH
92 k1 = eos_struct%UPARAM(3)
93
94 IF(iflag == 0) THEN
95 DO i=1,nel
96 dpdm(i) = k0*exp( k1*log(1+mu(i)) )/(one+mu(i))
97 dpde(i) = zero
98 pnew(i) = k0/k1 * (exp(k1*log(one+mu(i))) - one ) + p0
99 pnew(i) =
max(pnew(i),pmin)*off(i)
100 pnew(i) = pnew(i)-psh(i)
101 ENDDO
102
103 ELSEIF(iflag == 1) THEN
104 DO i=1,nel
105 pnew(i) = k0/k1 * (exp(k1*log(one+mu(i))) - one ) + p0
106 pnew(i) =
max(pnew(i),pmin)*off(i)
107 eint(i) = eint(i) - half*dvol(i)*(pnew(i)+psh(i))
108 pnew(i) = pnew(i)-psh(i)
109 dpde(i) = zero
110 ENDDO
111
112 ELSEIF(iflag == 2) THEN
113 DO i=1,nel
114 IF (vnew(i) > zero) THEN
115 pnew(i) = k0/k1 * (exp(k1*log(1+mu(i))) - one ) + p0
116 dpdm(i) = k0*exp( k1*log(1+mu(i)) )/(one+mu(i))
117 dpde(i) = zero
118 pnew(i) =
max(pnew(i),pmin)*off(i)
119 pnew(i) = pnew(i)-psh(i)
120 ENDIF
121 ENDDO
122 ENDIF
123
124 RETURN