Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 2 : program sample_eos
21 2 : use eos_def
22 : use eos_lib
23 : use chem_def
24 : use chem_lib
25 : use const_def, only: dp, ln10
26 : use const_lib, only: const_init
27 : use math_lib
28 :
29 : implicit none
30 :
31 4 : real(dp) :: X, Z, Y, abar, zbar, z2bar, z53bar, ye
32 : integer, parameter :: species = 7
33 : integer, parameter :: h1 = 1, he4 = 2, c12 = 3, n14 = 4, o16 = 5, ne20 = 6, mg24 = 7
34 2 : integer, pointer, dimension(:) :: net_iso, chem_id
35 16 : real(dp) :: xa(species)
36 : character(len=256) :: my_mesa_dir
37 :
38 2 : call Sample
39 :
40 : contains
41 :
42 4 : subroutine Sample
43 :
44 : integer :: handle
45 4 : real(dp) :: Rho, T, Pgas, log10Rho, log10T
46 2 : real(dp) :: dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
47 160 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
48 44 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa
49 : integer :: ierr
50 : character(len=*), parameter :: fmt1 = "(a20, 3x, e20.12)"
51 :
52 2 : ierr = 0
53 :
54 2 : my_mesa_dir = '..' ! if empty string, uses environment variable MESA_DIR
55 2 : call const_init(my_mesa_dir, ierr)
56 2 : if (ierr /= 0) then
57 0 : write (*, *) 'const_init failed'
58 0 : call mesa_error(__FILE__, __LINE__)
59 : end if
60 :
61 2 : call math_init()
62 :
63 2 : call chem_init('isotopes.data', ierr)
64 2 : if (ierr /= 0) then
65 0 : write (*, *) 'failed in chem_init'
66 0 : call mesa_error(__FILE__, __LINE__)
67 : end if
68 :
69 : ! allocate and initialize the eos tables
70 2 : call Setup_eos(handle)
71 :
72 2 : allocate (net_iso(num_chem_isos), chem_id(species), stat=ierr)
73 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'allocate failed')
74 2 : X = 0.70_dp
75 2 : Z = 0.02_dp
76 2 : call Init_Composition
77 :
78 2 : Rho = 1.3519d2
79 2 : log10T = 8.15_dp
80 2 : T = exp10(log10T)
81 :
82 : ! get a set of results for given temperature and density
83 : call eosDT_get( &
84 : handle, &
85 : species, chem_id, net_iso, xa, &
86 : Rho, log10(Rho), T, log10T, &
87 : res, d_dlnd, d_dlnT, &
88 2 : d_dxa, ierr)
89 :
90 2 : Pgas = exp(res(i_lnPgas))
91 :
92 : ! the indices for the results are defined in eos_def.f90
93 2 : write (*, '(A)')
94 2 : write (*, fmt1) 'temperature', T
95 2 : write (*, fmt1) 'density', Rho
96 2 : write (*, fmt1) 'logT', log10T
97 2 : write (*, fmt1) 'logRho', log10(Rho)
98 2 : write (*, '(A)')
99 2 : write (*, fmt1) 'Z', Z
100 2 : write (*, fmt1) 'X', X
101 2 : write (*, fmt1) 'abar', abar
102 2 : write (*, fmt1) 'zbar', zbar
103 2 : write (*, '(A)')
104 2 : write (*, fmt1) 'Pgas', Pgas
105 2 : write (*, fmt1) 'logPgas', res(i_lnPgas)/ln10
106 2 : write (*, fmt1) 'grad_ad', res(i_grad_ad)
107 2 : write (*, fmt1) 'c_P', res(i_Cp)
108 2 : write (*, '(A)')
109 :
110 : call eosPT_get( &
111 : handle, &
112 : species, chem_id, net_iso, xa, &
113 : Pgas, log10(Pgas), T, log10(T), &
114 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
115 2 : res, d_dlnd, d_dlnT, d_dxa, ierr)
116 :
117 : ! the indices for the results are defined in eos_def.f90
118 2 : write (*, '(A)')
119 2 : write (*, fmt1) 'temperature', T
120 2 : write (*, fmt1) 'Pgas', Pgas
121 2 : write (*, fmt1) 'logT', log10(T)
122 2 : write (*, fmt1) 'logPgas', res(i_lnPgas)/ln10
123 2 : write (*, '(A)')
124 2 : write (*, fmt1) 'Z', Z
125 2 : write (*, fmt1) 'X', X
126 2 : write (*, fmt1) 'abar', abar
127 2 : write (*, fmt1) 'zbar', zbar
128 2 : write (*, '(A)')
129 2 : write (*, fmt1) 'density', Rho
130 2 : write (*, fmt1) 'logRho', log10Rho
131 2 : write (*, fmt1) 'grad_ad', res(i_grad_ad)
132 2 : write (*, fmt1) 'c_P', res(i_Cp)
133 2 : write (*, '(A)')
134 :
135 : ! deallocate the eos tables
136 2 : call Shutdown_eos(handle)
137 :
138 2 : deallocate (net_iso, chem_id)
139 :
140 2 : if (ierr /= 0) then
141 0 : write (*, *) 'bad result from eos_get'
142 0 : call mesa_error(__FILE__, __LINE__)
143 : end if
144 :
145 8 : end subroutine Sample
146 :
147 2 : subroutine Setup_eos(handle)
148 : ! allocate and load the eos tables
149 : use eos_def
150 : use eos_lib
151 : integer, intent(out) :: handle
152 :
153 : integer :: ierr
154 : logical, parameter :: use_cache = .true.
155 :
156 2 : call eos_init('', use_cache, ierr)
157 2 : if (ierr /= 0) then
158 0 : write (*, *) 'eos_init failed in Setup_eos'
159 0 : call mesa_error(__FILE__, __LINE__)
160 : end if
161 :
162 2 : write (*, *) 'loading eos tables'
163 :
164 2 : handle = alloc_eos_handle(ierr)
165 2 : if (ierr /= 0) then
166 0 : write (*, *) 'failed trying to allocate eos handle'
167 0 : call mesa_error(__FILE__, __LINE__)
168 : end if
169 :
170 2 : end subroutine Setup_eos
171 :
172 2 : subroutine Shutdown_eos(handle)
173 2 : use eos_def
174 : use eos_lib
175 : integer, intent(in) :: handle
176 2 : call free_eos_handle(handle)
177 2 : call eos_shutdown
178 2 : end subroutine Shutdown_eos
179 :
180 2 : subroutine Init_Composition
181 2 : use chem_lib
182 :
183 : real(dp), parameter :: Zfrac_C = 0.173312d0
184 : real(dp), parameter :: Zfrac_N = 0.053177d0
185 : real(dp), parameter :: Zfrac_O = 0.482398d0
186 : real(dp), parameter :: Zfrac_Ne = 0.098675d0
187 :
188 32 : real(dp) :: dabar_dx(species), dzbar_dx(species), &
189 18 : sumx, xh, xhe, xz, mass_correction, dmc_dx(species)
190 :
191 15714 : net_iso(:) = 0
192 :
193 2 : chem_id(h1) = ih1; net_iso(ih1) = h1
194 2 : chem_id(he4) = ihe4; net_iso(ihe4) = he4
195 2 : chem_id(c12) = ic12; net_iso(ic12) = c12
196 2 : chem_id(n14) = in14; net_iso(in14) = n14
197 2 : chem_id(o16) = io16; net_iso(io16) = o16
198 2 : chem_id(ne20) = ine20; net_iso(ine20) = ne20
199 2 : chem_id(mg24) = img24; net_iso(img24) = mg24
200 :
201 2 : Y = 1 - (X + Z)
202 :
203 2 : xa(h1) = X
204 2 : xa(he4) = Y
205 2 : xa(c12) = Z*Zfrac_C
206 2 : xa(n14) = Z*Zfrac_N
207 2 : xa(o16) = Z*Zfrac_O
208 2 : xa(ne20) = Z*Zfrac_Ne
209 14 : xa(species) = 1 - sum(xa(1:species - 1))
210 :
211 : call composition_info( &
212 : species, chem_id, xa, xh, xhe, xz, abar, zbar, z2bar, z53bar, ye, &
213 2 : mass_correction, sumx, dabar_dx, dzbar_dx, dmc_dx)
214 :
215 2 : end subroutine Init_Composition
216 :
217 : end program sample_eos
|
