Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013 Josiah Schwab & 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 :
21 : module load_ecapture
22 :
23 : use rates_def
24 : use const_def, only: dp
25 : use chem_def, only: iso_name_length
26 :
27 : implicit none
28 :
29 : contains
30 :
31 1 : subroutine load_ecapture_data(ierr)
32 : integer, intent(out) :: ierr
33 1 : ierr = 0
34 1 : if (do_ecapture) then
35 1 : call load_ecapture_states_list(ierr)
36 1 : if (ierr /= 0) return
37 1 : call load_ecapture_transitions_list(ierr)
38 : end if
39 : end subroutine load_ecapture_data
40 :
41 :
42 1 : subroutine load_ecapture_states_list(ierr)
43 : use utils_lib
44 : use chem_lib, only: chem_get_iso_id
45 : use chem_def, only: iso_name_length
46 :
47 : integer, intent(out) :: ierr
48 : integer :: iounit, i, j, id
49 : character (len=256) :: filename, string
50 : integer :: nstates
51 : character(len=iso_name_length) :: iso
52 :
53 : real(dp) :: Ei, Ji
54 :
55 : logical, parameter :: dbg = .false.
56 :
57 : include 'formats'
58 :
59 1 : ierr = 0
60 :
61 1 : filename = trim(ecapture_states_file)
62 :
63 : if (dbg) then
64 : write(*,'(A)')
65 : write(*,*) 'ecapture states filename <' // trim(filename) // '>'
66 : write(*,'(A)')
67 : end if
68 :
69 : ierr = 0
70 1 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
71 1 : if (ierr /= 0) then
72 0 : write(*,*) 'failed to open special_weak_states_file ' // trim(filename)
73 0 : return
74 : end if
75 :
76 : allocate(ecapture_nuclide_name(max_num_ecapture_nuclei), &
77 1 : ecapture_nuclide_id(max_num_ecapture_nuclei))
78 :
79 : allocate(ecapture_states_data(max_num_ecapture_nuclei * &
80 1 : max_num_ecapture_states, num_states_data))
81 :
82 1 : nullify(ecapture_states_number_dict)
83 1 : nullify(ecapture_states_offset_dict)
84 :
85 1 : num_ecapture_nuclei = 0
86 1 : num_ecapture_states = 0
87 3 : do i = 1, max_num_ecapture_nuclei ! keep reading until end of file or max # nuclei
88 :
89 : do ! skip any blank or comment lines
90 9 : read(iounit,'(a)',iostat=ierr) string
91 9 : if (ierr /= 0) exit
92 9 : if ((index(string,"!") /= 0) .or. (len_trim(string) == 0)) then
93 : cycle ! comment or blank line
94 : else
95 6 : exit ! good line
96 : end if
97 : end do
98 3 : if (ierr /= 0) then
99 1 : ierr = 0; exit
100 : end if
101 :
102 : ! string now holds the first good line
103 2 : read(string,fmt=*,iostat=ierr) iso, nstates
104 2 : if (ierr /= 0) then
105 0 : ierr = 0; exit
106 : end if
107 2 : num_ecapture_nuclei = i
108 :
109 2 : if (nstates > max_num_ecapture_states) stop "ecapture: too many states"
110 2 : call integer_dict_define(ecapture_states_number_dict, iso, nstates, ierr)
111 :
112 2 : id = chem_get_iso_id(iso)
113 2 : if (id <= 0) then
114 0 : write(*,*) 'ecapture FATAL ERROR: unknown nuclide ' // iso
115 0 : call mesa_error(__FILE__,__LINE__)
116 : end if
117 :
118 2 : ecapture_nuclide_id(i) = id
119 2 : ecapture_nuclide_name(i) = iso
120 :
121 : if (dbg) write(*,'(a)') 'ecapture nucleus ' // trim(iso)
122 :
123 : ! store where this list of states starts
124 2 : call integer_dict_define(ecapture_states_offset_dict, iso, num_ecapture_states, ierr)
125 2 : if (failed('integer_dict_define')) return
126 :
127 14 : do j = 1, nstates
128 8 : read(iounit,fmt=*,iostat=ierr) Ei, Ji
129 8 : if (ierr /= 0) then
130 0 : ierr = 0; exit
131 : end if
132 8 : num_ecapture_states = num_ecapture_states + 1
133 :
134 : ! pack the data into the array
135 8 : ecapture_states_data(num_ecapture_states, i_E) = Ei
136 10 : ecapture_states_data(num_ecapture_states, i_J) = Ji
137 :
138 : end do
139 :
140 : end do
141 :
142 1 : close(iounit)
143 :
144 1 : if (num_ecapture_nuclei == 0) then
145 0 : ierr = -1
146 0 : write(*,*) 'failed trying to read special_weak_states_file -- no nuclei?'
147 0 : return
148 : end if
149 :
150 1 : if (num_ecapture_nuclei == max_num_ecapture_nuclei) then
151 0 : ierr = -1
152 0 : write(*,*) 'failed trying to read special_weak_states_file -- too many nuclei?'
153 0 : return
154 : end if
155 :
156 1 : call integer_dict_create_hash(ecapture_states_number_dict, ierr)
157 1 : if (ierr /= 0) return
158 :
159 1 : call integer_dict_create_hash(ecapture_states_offset_dict, ierr)
160 1 : if (ierr /= 0) return
161 :
162 1 : call realloc_double2(ecapture_states_data, num_ecapture_states, num_states_data, ierr)
163 1 : if (ierr /= 0) return
164 :
165 1 : call realloc_integer(ecapture_nuclide_id, num_ecapture_nuclei, ierr)
166 3 : if (ierr /= 0) return
167 :
168 : contains
169 :
170 2 : logical function failed(str)
171 : character (len=*) :: str
172 2 : failed = (ierr /= 0)
173 2 : if (failed) then
174 0 : write(*,*) 'failed: ' // trim(str)
175 : end if
176 1 : end function failed
177 :
178 :
179 : end subroutine load_ecapture_states_list
180 :
181 :
182 1 : subroutine load_ecapture_transitions_list(ierr)
183 : use utils_lib
184 : use chem_lib, only: chem_get_iso_id
185 : use chem_def, only: iso_name_length
186 :
187 : integer, intent(out) :: ierr
188 : integer :: iounit, i, j, id
189 : character (len=256) :: filename, string
190 : character(len=iso_name_length) :: lhs, rhs
191 : integer :: ntrans
192 : character(len=2*iso_name_length+1) :: key
193 :
194 : integer :: Si, Sf
195 : real(dp) :: logft
196 :
197 : logical, parameter :: dbg = .false.
198 :
199 : include 'formats'
200 :
201 1 : ierr = 0
202 :
203 1 : filename = ecapture_transitions_file
204 :
205 : if (dbg) then
206 : write(*,'(A)')
207 : write(*,*) 'ecapture transitions filename <' // trim(filename) // '>'
208 : write(*,'(A)')
209 : end if
210 :
211 : ierr = 0
212 1 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
213 1 : if (ierr /= 0) then
214 0 : write(*,*) 'failed to open special_weak_transitions_file ' // trim(filename)
215 0 : write(*,*) trim(string)
216 0 : return
217 : end if
218 :
219 : allocate(ecapture_lhs_nuclide_name(max_num_ecapture_reactions), &
220 : ecapture_rhs_nuclide_name(max_num_ecapture_reactions), &
221 : ecapture_lhs_nuclide_id(max_num_ecapture_reactions), &
222 1 : ecapture_rhs_nuclide_id(max_num_ecapture_reactions))
223 :
224 : allocate(ecapture_transitions_data(max_num_ecapture_reactions * &
225 : max_num_ecapture_transitions, &
226 1 : num_transitions_data))
227 : allocate(ecapture_logft_data(max_num_ecapture_reactions * &
228 1 : max_num_ecapture_transitions))
229 :
230 1 : nullify(ecapture_reactions_dict)
231 1 : nullify(ecapture_transitions_number_dict)
232 1 : nullify(ecapture_transitions_offset_dict)
233 1 : num_ecapture_reactions = 0
234 1 : num_ecapture_transitions = 0
235 :
236 3 : do i = 1, max_num_ecapture_reactions ! keep reading until end of file
237 :
238 : do ! skip any blank or comment lines
239 17 : read(iounit,'(a)',iostat=ierr) string
240 17 : if (ierr /= 0) exit
241 17 : if ((index(string,"!") /= 0) .or. (len_trim(string) == 0)) then
242 : cycle ! comment or blank line
243 : else
244 14 : exit ! good line
245 : end if
246 : end do
247 3 : if (ierr /= 0) then
248 1 : ierr = 0; exit
249 : end if
250 :
251 : ! string now holds the first good line
252 2 : read(string,fmt=*,iostat=ierr) lhs, rhs, ntrans
253 2 : if (ierr /= 0) then
254 0 : ierr = 0; exit
255 : end if
256 :
257 2 : id = chem_get_iso_id(lhs)
258 2 : if (id <= 0) then
259 0 : write(*,*) 'ecapture FATAL ERROR: unknown nuclide ' // lhs
260 0 : call mesa_error(__FILE__,__LINE__)
261 : end if
262 2 : ecapture_lhs_nuclide_id(i) = id
263 2 : id = chem_get_iso_id(rhs)
264 2 : if (id <= 0) then
265 0 : write(*,*) 'ecapture FATAL ERROR: unknown nuclide ' // rhs
266 0 : call mesa_error(__FILE__,__LINE__)
267 : end if
268 2 : ecapture_rhs_nuclide_id(i) = id
269 2 : ecapture_lhs_nuclide_name(i) = lhs
270 2 : ecapture_rhs_nuclide_name(i) = rhs
271 :
272 2 : call create_ecapture_dict_key(lhs, rhs, key)
273 :
274 : if (dbg) write(*,'(a)') 'ecapture transitions ' // trim(key)
275 :
276 2 : call integer_dict_define(ecapture_reactions_dict, key, i, ierr)
277 2 : if (failed('integer_dict_define')) return
278 2 : num_ecapture_reactions = i
279 :
280 2 : call integer_dict_define(ecapture_transitions_number_dict, key, ntrans, ierr)
281 2 : if (failed('integer_dict_define')) return
282 :
283 2 : call integer_dict_define(ecapture_transitions_offset_dict, key, num_ecapture_transitions, ierr)
284 2 : if (failed('integer_dict_define')) return
285 :
286 12 : do j = 1, ntrans
287 4 : read(iounit,fmt=*,iostat=ierr) Si, Sf, logft
288 4 : if (ierr /= 0) then
289 0 : ierr = 0; exit
290 : end if
291 4 : num_ecapture_transitions = num_ecapture_transitions + 1
292 :
293 : ! pack the data into the array
294 4 : ecapture_transitions_data(num_ecapture_transitions, i_Si) = Si
295 4 : ecapture_transitions_data(num_ecapture_transitions, i_Sf) = Sf
296 :
297 6 : ecapture_logft_data(num_ecapture_transitions) = logft
298 :
299 : end do
300 :
301 : end do
302 :
303 1 : close(iounit)
304 :
305 1 : if (num_ecapture_reactions == 0) then
306 0 : ierr = -1
307 0 : write(*,*) 'failed trying to read special_weak_transitions_file -- no reactions?'
308 0 : return
309 : end if
310 :
311 1 : if (num_ecapture_reactions == max_num_ecapture_reactions) then
312 0 : ierr = -1
313 0 : write(*,*) 'failed trying to read special_weak_transitions_file -- too many reactions?'
314 0 : return
315 : end if
316 :
317 1 : call integer_dict_create_hash(ecapture_reactions_dict, ierr)
318 1 : if (ierr /= 0) return
319 :
320 1 : call integer_dict_create_hash(ecapture_transitions_number_dict, ierr)
321 1 : if (ierr /= 0) return
322 :
323 1 : call integer_dict_create_hash(ecapture_transitions_offset_dict, ierr)
324 1 : if (ierr /= 0) return
325 :
326 1 : call realloc_integer2(ecapture_transitions_data, num_ecapture_transitions, num_transitions_data, ierr)
327 1 : if (ierr /= 0) return
328 :
329 1 : call realloc_double(ecapture_logft_data, num_ecapture_transitions, ierr)
330 1 : if (ierr /= 0) return
331 :
332 1 : call realloc_integer(ecapture_lhs_nuclide_id, num_ecapture_reactions, ierr)
333 1 : if (ierr /= 0) return
334 :
335 1 : call realloc_integer(ecapture_rhs_nuclide_id, num_ecapture_reactions, ierr)
336 6 : if (ierr /= 0) return
337 :
338 : contains
339 :
340 6 : logical function failed(str)
341 : character (len=*) :: str
342 6 : failed = (ierr /= 0)
343 6 : if (failed) then
344 0 : write(*,*) 'failed: ' // trim(str)
345 : end if
346 1 : end function failed
347 :
348 :
349 : end subroutine load_ecapture_transitions_list
350 :
351 :
352 :
353 : end module load_ecapture
|
