Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Bill Paxton & 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 : module mod_ionization
21 :
22 : use const_def, only: one_third, two_thirds, avo, dp, mesa_data_dir
23 : use math_lib
24 : use ionization_def
25 : use utils_lib, only: mesa_error
26 :
27 : implicit none
28 :
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 :
33 : contains
34 :
35 :
36 3 : subroutine do_init_ionization(ionization_cache_dir_in, use_cache, ierr)
37 : character (len=*), intent(in) :: ionization_cache_dir_in
38 : logical, intent(in) :: use_cache
39 : integer, intent(out) :: ierr
40 3 : ierr = 0
41 3 : table_is_initialized = .false.
42 3 : end subroutine do_init_ionization
43 :
44 :
45 2 : subroutine do_load(ierr)
46 : use ionization_def
47 : integer, intent(out) :: ierr
48 :
49 : integer :: io_log_ne, io_logT, io_z
50 : integer, pointer :: ibound(:,:), tmp_version(:)
51 : integer, parameter :: num_log_ne_fe56_he4 = 105, num_logT_fe56_he4 = 30
52 :
53 : ierr = 0
54 :
55 2 : fe_he_ptr => fe_he_info
56 :
57 : call load_table_summary( &
58 : 'log_ne_fe56_he4.data', 'logT_fe56_he4.data', 'z_fe56_he4.data', &
59 2 : num_log_ne_fe56_he4, num_logT_fe56_he4, fe_he_ptr, ierr)
60 2 : if (ierr /= 0) return
61 :
62 2 : call create_interpolants(fe_he_ptr,num_log_ne_fe56_he4,num_logT_fe56_he4,ierr)
63 2 : if (ierr /= 0) return
64 :
65 2 : table_is_initialized = .true.
66 :
67 : contains
68 :
69 6 : subroutine openfile(filename, iounit, ierr)
70 : character(len=*) :: filename
71 : integer, intent(inout) :: iounit
72 : integer, intent(out) :: ierr
73 : if (dbg) write(*,*) 'read ' // trim(filename)
74 6 : ierr = 0
75 6 : open(newunit=iounit,file=trim(filename),action='read',status='old',iostat=ierr)
76 6 : if (ierr/= 0) then
77 0 : write(*,*) 'table_ionization_init: missing ionization data'
78 0 : write(*,*) filename
79 0 : write(*,*)
80 0 : write(*,*)
81 0 : write(*,*)
82 0 : write(*,*)
83 0 : write(*,*)
84 0 : write(*,*) 'FATAL ERROR: missing or bad ionization data.'
85 0 : write(*,*) 'Please update by removing the directory mesa/data/ionization_data,'
86 0 : write(*,*) 'and rerunning the mesa ./install script.'
87 0 : write(*,*)
88 0 : call mesa_error(__FILE__,__LINE__)
89 : end if
90 6 : end subroutine openfile
91 :
92 :
93 2 : subroutine load_table_summary( &
94 : log_ne_fname, logT_fname, z_fname, num_log_ne, num_logT, p, ierr)
95 : character(len=*), intent(in) :: log_ne_fname, logT_fname, z_fname
96 : integer, intent(in) :: num_log_ne, num_logT
97 : type (Ionization_Info), pointer :: p
98 : integer, intent(out) :: ierr
99 :
100 : character(len=256) :: filename
101 2 : real(dp), pointer :: f(:,:,:)
102 : integer :: i, j
103 :
104 2 : ierr = 0
105 2 : p% have_interpolation_info = .false.
106 2 : p% num_log_ne = num_log_ne
107 2 : p% num_logT = num_logT
108 : allocate( &
109 : p% log_ne(num_log_ne), p% logT(num_logT), &
110 2 : p% f1(4*num_log_ne*num_logT), stat=ierr)
111 2 : if (ierr /= 0) then
112 0 : write(*,*) 'failed in allocate for ionization tables'
113 0 : call mesa_error(__FILE__,__LINE__)
114 : end if
115 2 : f(1:4,1:num_log_ne,1:num_logT) => p% f1(1:4*num_log_ne*num_logT)
116 :
117 2 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(z_fname)
118 2 : call openfile(filename, io_z, ierr)
119 2 : if (ierr /= 0) return
120 62 : do i=1,num_logT
121 60 : read(io_z,fmt=*,iostat=ierr) p% log_ne(1:num_log_ne)
122 60 : if (ierr /= 0) then
123 0 : write(*,*) 'failed in reading ionization z ' // trim(filename)
124 0 : call mesa_error(__FILE__,__LINE__)
125 : end if
126 : !p% f(1,1:num_log_ne,i) = p% log_ne(1:num_log_ne) << segfault on UBUNTU
127 6362 : do j=1,num_log_ne
128 6360 : f(1,j,i) = p% log_ne(j) ! sets p% f1
129 : end do
130 : end do
131 2 : close(io_z)
132 :
133 2 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(log_ne_fname)
134 2 : call openfile(filename, io_log_ne, ierr)
135 2 : if (ierr /= 0) return
136 212 : do i=1,num_log_ne
137 210 : read(io_log_ne,fmt=*,iostat=ierr) p% log_ne(i)
138 212 : if (ierr /= 0) then
139 0 : write(*,*) 'failed in reading ionization log_ne ' // trim(filename)
140 0 : call mesa_error(__FILE__,__LINE__)
141 : end if
142 : end do
143 2 : close(io_log_ne)
144 :
145 2 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(logT_fname)
146 2 : call openfile(filename, io_logT, ierr)
147 2 : if (ierr /= 0) return
148 62 : do i=1,num_logT
149 60 : read(io_logT,fmt=*,iostat=ierr) p% logT(i)
150 62 : if (ierr /= 0) then
151 0 : write(*,*) 'failed in reading ionization logT ' // trim(filename)
152 0 : call mesa_error(__FILE__,__LINE__)
153 : end if
154 : end do
155 2 : close(io_logT)
156 :
157 8 : end subroutine load_table_summary
158 :
159 :
160 : end subroutine do_load
161 :
162 :
163 2 : subroutine create_interpolants(p,nx,ny,ierr)
164 : use interp_2d_lib_db
165 : type (Ionization_Info), pointer :: p
166 : integer, intent(in) :: nx, ny
167 : integer, intent(out) :: ierr
168 : integer :: ibcxmin, ibcxmax, ibcymin, ibcymax
169 542 : real(dp) :: bcxmin(ny), bcxmax(ny), bcymin(nx), bcymax(nx)
170 : ! use "not a knot" bc's
171 62 : ibcxmin = 0; bcxmin(:) = 0d0
172 62 : ibcxmax = 0; bcxmax(:) = 0d0
173 212 : ibcymin = 0; bcymin(:) = 0d0
174 212 : ibcymax = 0; bcymax(:) = 0d0
175 : call interp_mkbicub_db( &
176 : p% log_ne, p% num_log_ne, p% logT, p% num_logT, &
177 : p% f1, p% num_log_ne, &
178 : ibcxmin, bcxmin, ibcxmax, bcxmax, &
179 : ibcymin, bcymin, ibcymax, bcymax, &
180 2 : p% ilinx, p% iliny, ierr )
181 2 : if (ierr /= 0) then
182 : if (dbg) write(*,*) 'interp_mkbicub_db failed for ionization interpolant'
183 : return
184 : end if
185 2 : p% have_interpolation_info = .true.
186 : end subroutine create_interpolants
187 :
188 6 : real(dp) function charge_of_Fe56_in_He4(log_ne, logT, ierr)
189 : use interp_2d_lib_db
190 : real(dp), intent(in) :: log_ne ! ne=avo*rho*free_e
191 : real(dp), intent(in) :: logT
192 : integer, intent(out) :: ierr
193 :
194 : integer :: ict(6) ! code specifying output desired
195 42 : real(dp) :: fval(6) ! output data
196 : type (Ionization_Info), pointer :: p
197 :
198 6 : ierr = 0
199 6 : charge_of_Fe56_in_He4 = 0
200 :
201 6 : if (.not. table_is_initialized) then
202 4 : !$omp critical (ionization_table)
203 2 : if (.not. table_is_initialized) call do_load(ierr)
204 : !$omp end critical (ionization_table)
205 2 : if (ierr /= 0) return
206 : end if
207 :
208 6 : ict = 0; ict(1) = 1 ! just the result; no partials
209 6 : p => fe_he_ptr
210 : call interp_evbicub_db( &
211 : log_ne, logT, p% log_ne, p% num_log_ne, p% logT, p% num_logT, &
212 6 : p% ilinx, p% iliny, p% f1, p% num_log_ne, ict, fval, ierr)
213 :
214 6 : charge_of_Fe56_in_He4 = fval(1)
215 :
216 6 : end function charge_of_Fe56_in_He4
217 :
218 0 : subroutine chi_info(a1, z1, T, log_T, rho, log_rho, chi, c0, c1, c2)
219 : real(dp), intent(in) :: a1, z1, T, log_T, rho, log_rho
220 : real(dp), intent(out) :: chi, c0, c1, c2
221 0 : chi = 1.987d-4*T*(-8.392d0 - log_rho + 1.5d0*log_T - log10(z1/a1)) ! eqn 20
222 : ! coef's used in eqn 21
223 0 : c0 = 1.085d-4*rho*T/a1
224 0 : c1 = 1.617d4*sqrt(rho*(z1*z1 + z1)/(T*a1))
225 0 : c2 = 29.38d0*z1*pow(rho/a1,one_third)
226 : ! c2 had a typo in eqn 21, now corrected to match Dupuis et al. (1992) eqn 3
227 0 : end subroutine chi_info
228 :
229 0 : real(dp) function chi_effective(chi, c0, c1, c2, z1, z2)
230 : real(dp), intent(in) :: chi, c0, c1, c2, z1, z2
231 : chi_effective = chi + c0/(z2*z2*z2) + &
232 0 : min(c1*z2, c2*(pow(z2/z1,two_thirds) + 0.6d0))
233 0 : end function chi_effective
234 :
235 :
236 : end module mod_ionization
237 :
|
