-
Notifications
You must be signed in to change notification settings - Fork 12
/
mass_interior.F90
125 lines (103 loc) · 4.83 KB
/
mass_interior.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
!-*-f90-*-
subroutine mass_interior
! calculate baryonic mass interior to cell centers
! for innermost zone: use that zone inner radius and assume that its
! density is constant interior to the inner cell radius
use GR1D_module
implicit none
integer i
real*8 dmass,dx1,dx2,discrim
real*8 h,rhot,massint,phi_bound
real*8 include_nus
mass(:) = 0.0d0
if (GR) then
mass1 = X * rho * W * volume
else
mass1 = rho * volume
endif
if (do_M1) then
include_nus = 1.0d0
else
include_nus = 0.0d0
endif
if (GR) then
mass(ghosts1+1) = mass1(ghosts1+1) - 4.0d0/3.0d0*pi*rho(ghosts1+1) * X(ghosts1+1) * &
W(ghosts1+1)*( x1i(ghosts1+1+1)**3 - x1(ghosts1+1)**3 )
do i=ghosts1+2,n1-1
mass(i) = mass(i-1) + mass1(i) &
- 4.0d0/3.0d0*pi*rho(i)*X(i)*W(i) * ( x1i(i+1)**3 - x1(i)**3 ) &
+ 4.0d0/3.0d0*pi*rho(i-1)*X(i-1)*W(i-1) * (x1i(i)**3 - x1(i-1)**3)
enddo
mass(n1) = mass(n1-1) &
+ 4.0d0/3.0d0*pi*rho(n1-1)*X(n1-1)*W(n1-1)*(x1i(n1)**3 - x1(n1-1)**3) &
+ 4.0d0/3.0d0*pi*rho(n1)*X(n1)*W(n1)*( x1(n1)**3 - x1i(n1)**3 )
else
if (do_effectivepotential) then
mass(ghosts1+1) = 4.0d0*pi/3.0d0*x1(ghosts1+1)**3*(rho(ghosts1+1)* &
(1.0d0 + eps(ghosts1+1)) + (energy_nu(ghosts1+1)+mom_nu(ghosts1+1)*v1(ghosts1+1))*include_nus)
mass(ghosts1+1) = mass(ghosts1+1)*sqrt(1.0d0-2.0d0*mass(ghosts1+1)/x1(ghosts1+1))
dphidr(ghosts1+1) = (mass(ghosts1+1) + 4.0d0*pi*x1(ghosts1+1)**3* &
(press(ghosts1+1)+press_nu(ghosts1+1)))/ &
(x1(ghosts1+1)**2*(1.0d0+v1(ghosts1+1)**2-2.0d0*mass(ghosts1+1)/x1(ghosts1+1)))* &
(rho(ghosts1+1)+eps(ghosts1+1)*rho(ghosts1+1)+press(ghosts1+1))/rho(ghosts1+1)
do i=ghosts1+2,n1-1
mass(i) = mass(i-1) + &
4.0d0/3.0d0*pi*(rho(i-1)*(1.0d0+eps(i-1))+(energy_nu(i-1)+mom_nu(i-1)*v1(i-1))*include_nus) &
* ( x1i(i)**3 - x1(i-1)**3 )*sqrt(1.0d0-2.0d0*mass(i-1)/x1(i-1))
mass(i) = mass(i) + &
4.0d0/3.0d0*pi*(rho(i)*(1.0d0+eps(i))+(energy_nu(i)+mom_nu(i)*v1(i))*include_nus) * &
(x1(i)**3 - x1i(i)**3)*sqrt(1.0d0-2.0d0*mass(i)/x1i(i))
dphidr(i) = (mass(i) + 4.0d0*pi*x1(i)**3*(press(i)+press_nu(i)))/ &
(x1(i)**2*(1.0d0+v1(i)**2-2.0d0*mass(i)/x1(i)))* &
(rho(i)+eps(i)*rho(i)+press(i))/rho(i)
enddo
mass(n1) = mass(n1-1) + &
4.0d0/3.0d0*pi*(rho(n1-1)*(1.0d0+eps(n1-1))+(energy_nu(n1-1)+mom_nu(n1-1)*v1(n1-1))*include_nus)* &
(x1i(n1)**3 - x1(n1-1)**3)*sqrt(1.0d0-2.0d0*mass(n1-1)/x1(n1-1))
mass(n1) = mass(n1) + &
4.0d0/3.0d0*pi*(rho(n1)*(1.0d0+eps(n1))+(energy_nu(n1)+mom_nu(n1)*v1(n1))*include_nus)* &
(x1(n1)**3 - x1i(n1)**3)*sqrt(1.0d0-2.0d0*mass(n1)/x1i(n1))
dphidr(n1) = (mass(n1) + 4.0d0*pi*x1(n1)**3*(press(n1)+press_nu(n1)))/ &
(x1(n1)**2*(1.0d0+v1(n1)**2-2.0d0*mass(n1)/x1(n1)))* &
(rho(n1)+eps(n1)*rho(n1)+press(n1))/rho(n1)
phi(ghosts1+1) = 0.0d0 + x1(ghosts1+1)*(dphidr(ghosts1+1))
phii(ghosts1+1) = 0.0d0
do i=ghosts1+2,n1-ghosts1+1
phi(i) = phi(i-1) + (x1i(i)-x1(i-1))*dphidr(i-1) &
+ (x1(i)-x1i(i))*dphidr(i)
phii(i) = phi(i) - (x1(i)-x1i(i))*dphidr(i)
enddo
phii(n1-ghosts1+2) = phi(n1-ghosts1+1) + &
(x1i(n1-ghosts1+2)-x1(n1-ghosts1+1))*dphidr(n1-ghosts1+1)
phi_bound = 0.5d0 * &
log(1.0d0 - 2.0d0 * mass(n1-ghosts1+1)/x1(n1-ghosts1+1))
do i=ghosts1+1,n1-ghosts1+1
phi(i) = phi(i) + (phi_bound-phi(n1-ghosts1+1))
phii(i) = phii(i) + (phi_bound-phi(n1-ghosts1+1))
enddo
do i=ghosts1+1,n1-ghosts1+1
alp(i) = exp(phi(i))
alpm(i) = exp(phii(i))
alpp(i) = exp(phii(i+1))
enddo
dphidr(ghosts1+1) = (phi(ghosts1+2)-phi(ghosts1+1))/(x1(ghosts1+2)-x1(ghosts1+1))
do i=ghosts1+2,n1-ghosts1
dphidr(i) = (phi(i+1)-phi(i-1))/(x1(i+1)-x1(i-1))
enddo
else
mass(ghosts1+1) = mass1(ghosts1+1) - 4.0d0/3.0d0*pi*rho(ghosts1+1)* &
( x1i(ghosts1+1+1)**3 - x1(ghosts1+1)**3 )
dphidr(ghosts1+1) = mass(ghosts1+1)/x1(ghosts1+1)**2
do i=ghosts1+2,n1-1
mass(i) = mass(i-1) + mass1(i) &
- 4.0d0/3.0d0*pi*rho(i) * ( x1i(i+1)**3 - x1(i)**3 ) &
+ 4.0d0/3.0d0*pi*rho(i-1) * (x1i(i)**3 - x1(i-1)**3)
dphidr(i) = mass(i)/x1(i)**2
enddo
mass(n1) = mass(n1-1) &
+ 4.0d0/3.0d0*pi*rho(n1-1)*(x1i(n1)**3 - x1(n1-1)**3) &
+ 4.0d0/3.0d0*pi*rho(n1)*( x1(n1)**3 - x1i(n1)**3 )
dphidr(n1) = mass(n1)/x1(n1)**2
endif
endif
end subroutine mass_interior