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 2 : 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 2 : ierr = 0
41 2 : table_is_initialized = .false.
42 2 : end subroutine do_init_ionization
43 :
44 :
45 1 : 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 1 : 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 1 : num_log_ne_fe56_he4, num_logT_fe56_he4, fe_he_ptr, ierr)
60 1 : if (ierr /= 0) return
61 :
62 1 : call create_interpolants(fe_he_ptr,num_log_ne_fe56_he4,num_logT_fe56_he4,ierr)
63 1 : if (ierr /= 0) return
64 :
65 1 : table_is_initialized = .true.
66 :
67 : contains
68 :
69 3 : 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 3 : ierr = 0
75 3 : open(newunit=iounit,file=trim(filename),action='read',status='old',iostat=ierr)
76 3 : 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 3 : end subroutine openfile
91 :
92 :
93 1 : 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 1 : real(dp), pointer :: f(:,:,:)
102 : integer :: i, j
103 :
104 1 : ierr = 0
105 1 : p% have_interpolation_info = .false.
106 1 : p% num_log_ne = num_log_ne
107 1 : p% num_logT = num_logT
108 : allocate( &
109 : p% log_ne(num_log_ne), p% logT(num_logT), &
110 1 : p% f1(4*num_log_ne*num_logT), stat=ierr)
111 1 : if (ierr /= 0) then
112 0 : write(*,*) 'failed in allocate for ionization tables'
113 0 : call mesa_error(__FILE__,__LINE__)
114 : end if
115 1 : f(1:4,1:num_log_ne,1:num_logT) => p% f1(1:4*num_log_ne*num_logT)
116 :
117 1 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(z_fname)
118 1 : call openfile(filename, io_z, ierr)
119 1 : if (ierr /= 0) return
120 31 : do i=1,num_logT
121 30 : read(io_z,fmt=*,iostat=ierr) p% log_ne(1:num_log_ne)
122 30 : 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 3181 : do j=1,num_log_ne
128 3180 : f(1,j,i) = p% log_ne(j) ! sets p% f1
129 : end do
130 : end do
131 1 : close(io_z)
132 :
133 1 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(log_ne_fname)
134 1 : call openfile(filename, io_log_ne, ierr)
135 1 : if (ierr /= 0) return
136 106 : do i=1,num_log_ne
137 105 : read(io_log_ne,fmt=*,iostat=ierr) p% log_ne(i)
138 106 : 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 1 : close(io_log_ne)
144 :
145 1 : filename = trim(mesa_data_dir) // '/ionization_data/' // trim(logT_fname)
146 1 : call openfile(filename, io_logT, ierr)
147 1 : if (ierr /= 0) return
148 31 : do i=1,num_logT
149 30 : read(io_logT,fmt=*,iostat=ierr) p% logT(i)
150 31 : 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 1 : close(io_logT)
156 :
157 4 : end subroutine load_table_summary
158 :
159 :
160 : end subroutine do_load
161 :
162 :
163 1 : 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 0 : real(dp) :: bcxmin(ny), bcxmax(ny), bcymin(nx), bcymax(nx)
170 : ! use "not a knot" bc's
171 31 : ibcxmin = 0; bcxmin(:) = 0d0
172 31 : ibcxmax = 0; bcxmax(:) = 0d0
173 106 : ibcymin = 0; bcymin(:) = 0d0
174 106 : 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 1 : p% ilinx, p% iliny, ierr )
181 1 : if (ierr /= 0) then
182 : if (dbg) write(*,*) 'interp_mkbicub_db failed for ionization interpolant'
183 : return
184 : end if
185 1 : 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 : real(dp) :: fval(6) ! output data
196 : type (Ionization_Info), pointer :: p
197 :
198 3 : ierr = 0
199 3 : charge_of_Fe56_in_He4 = 0
200 :
201 3 : if (.not. table_is_initialized) then
202 2 : !$omp critical (ionization_table)
203 1 : if (.not. table_is_initialized) call do_load(ierr)
204 : !$omp end critical (ionization_table)
205 1 : if (ierr /= 0) return
206 : end if
207 :
208 3 : ict = 0; ict(1) = 1 ! just the result; no partials
209 3 : 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 3 : p% ilinx, p% iliny, p% f1, p% num_log_ne, ict, fval, ierr)
213 :
214 3 : charge_of_Fe56_in_He4 = fval(1)
215 :
216 3 : 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 :
|
