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 2 : subroutine load_ecapture_data(ierr)
32 : integer, intent(out) :: ierr
33 2 : ierr = 0
34 2 : if (do_ecapture) then
35 2 : call load_ecapture_states_list(ierr)
36 2 : if (ierr /= 0) return
37 2 : call load_ecapture_transitions_list(ierr)
38 : end if
39 : end subroutine load_ecapture_data
40 :
41 :
42 2 : 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 2 : real(dp) :: Ei, Ji
54 :
55 : logical, parameter :: dbg = .false.
56 :
57 : include 'formats'
58 :
59 2 : ierr = 0
60 :
61 2 : 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 2 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
71 2 : 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 2 : ecapture_nuclide_id(max_num_ecapture_nuclei))
78 :
79 : allocate(ecapture_states_data(max_num_ecapture_nuclei * &
80 2 : max_num_ecapture_states, num_states_data))
81 :
82 2 : nullify(ecapture_states_number_dict)
83 2 : nullify(ecapture_states_offset_dict)
84 :
85 2 : num_ecapture_nuclei = 0
86 2 : num_ecapture_states = 0
87 6 : 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 18 : read(iounit,'(a)',iostat=ierr) string
91 18 : if (ierr /= 0) exit
92 18 : if ((index(string,"!") /= 0) .or. (len_trim(string) == 0)) then
93 : cycle ! comment or blank line
94 : else
95 12 : exit ! good line
96 : end if
97 : end do
98 6 : if (ierr /= 0) then
99 2 : ierr = 0; exit
100 : end if
101 :
102 : ! string now holds the first good line
103 4 : read(string,fmt=*,iostat=ierr) iso, nstates
104 4 : if (ierr /= 0) then
105 0 : ierr = 0; exit
106 : end if
107 4 : num_ecapture_nuclei = i
108 :
109 4 : if (nstates > max_num_ecapture_states) stop "ecapture: too many states"
110 4 : call integer_dict_define(ecapture_states_number_dict, iso, nstates, ierr)
111 :
112 4 : id = chem_get_iso_id(iso)
113 4 : 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 4 : ecapture_nuclide_id(i) = id
119 4 : 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 4 : call integer_dict_define(ecapture_states_offset_dict, iso, num_ecapture_states, ierr)
125 4 : if (failed('integer_dict_define')) return
126 :
127 28 : do j = 1, nstates
128 16 : read(iounit,fmt=*,iostat=ierr) Ei, Ji
129 16 : if (ierr /= 0) then
130 0 : ierr = 0; exit
131 : end if
132 16 : num_ecapture_states = num_ecapture_states + 1
133 :
134 : ! pack the data into the array
135 16 : ecapture_states_data(num_ecapture_states, i_E) = Ei
136 20 : ecapture_states_data(num_ecapture_states, i_J) = Ji
137 :
138 : end do
139 :
140 : end do
141 :
142 2 : close(iounit)
143 :
144 2 : 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 2 : 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 2 : call integer_dict_create_hash(ecapture_states_number_dict, ierr)
157 2 : if (ierr /= 0) return
158 :
159 2 : call integer_dict_create_hash(ecapture_states_offset_dict, ierr)
160 2 : if (ierr /= 0) return
161 :
162 2 : call realloc_double2(ecapture_states_data, num_ecapture_states, num_states_data, ierr)
163 2 : if (ierr /= 0) return
164 :
165 2 : call realloc_integer(ecapture_nuclide_id, num_ecapture_nuclei, ierr)
166 6 : if (ierr /= 0) return
167 :
168 : contains
169 :
170 4 : logical function failed(str)
171 : character (len=*) :: str
172 4 : failed = (ierr /= 0)
173 4 : if (failed) then
174 0 : write(*,*) 'failed: ' // trim(str)
175 : end if
176 2 : end function failed
177 :
178 :
179 : end subroutine load_ecapture_states_list
180 :
181 :
182 2 : 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 2 : real(dp) :: logft
196 :
197 : logical, parameter :: dbg = .false.
198 :
199 : include 'formats'
200 :
201 2 : ierr = 0
202 :
203 2 : 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 2 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
213 2 : 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 2 : ecapture_rhs_nuclide_id(max_num_ecapture_reactions))
223 :
224 : allocate(ecapture_transitions_data(max_num_ecapture_reactions * &
225 : max_num_ecapture_transitions, &
226 2 : num_transitions_data))
227 : allocate(ecapture_logft_data(max_num_ecapture_reactions * &
228 2 : max_num_ecapture_transitions))
229 :
230 2 : nullify(ecapture_reactions_dict)
231 2 : nullify(ecapture_transitions_number_dict)
232 2 : nullify(ecapture_transitions_offset_dict)
233 2 : num_ecapture_reactions = 0
234 2 : num_ecapture_transitions = 0
235 :
236 6 : do i = 1, max_num_ecapture_reactions ! keep reading until end of file
237 :
238 : do ! skip any blank or comment lines
239 34 : read(iounit,'(a)',iostat=ierr) string
240 34 : if (ierr /= 0) exit
241 34 : if ((index(string,"!") /= 0) .or. (len_trim(string) == 0)) then
242 : cycle ! comment or blank line
243 : else
244 28 : exit ! good line
245 : end if
246 : end do
247 6 : if (ierr /= 0) then
248 2 : ierr = 0; exit
249 : end if
250 :
251 : ! string now holds the first good line
252 4 : read(string,fmt=*,iostat=ierr) lhs, rhs, ntrans
253 4 : if (ierr /= 0) then
254 0 : ierr = 0; exit
255 : end if
256 :
257 4 : id = chem_get_iso_id(lhs)
258 4 : if (id <= 0) then
259 0 : write(*,*) 'ecapture FATAL ERROR: unknown nuclide ' // lhs
260 0 : call mesa_error(__FILE__,__LINE__)
261 : end if
262 4 : ecapture_lhs_nuclide_id(i) = id
263 4 : id = chem_get_iso_id(rhs)
264 4 : if (id <= 0) then
265 0 : write(*,*) 'ecapture FATAL ERROR: unknown nuclide ' // rhs
266 0 : call mesa_error(__FILE__,__LINE__)
267 : end if
268 4 : ecapture_rhs_nuclide_id(i) = id
269 4 : ecapture_lhs_nuclide_name(i) = lhs
270 4 : ecapture_rhs_nuclide_name(i) = rhs
271 :
272 4 : call create_ecapture_dict_key(lhs, rhs, key)
273 :
274 : if (dbg) write(*,'(a)') 'ecapture transitions ' // trim(key)
275 :
276 4 : call integer_dict_define(ecapture_reactions_dict, key, i, ierr)
277 4 : if (failed('integer_dict_define')) return
278 4 : num_ecapture_reactions = i
279 :
280 4 : call integer_dict_define(ecapture_transitions_number_dict, key, ntrans, ierr)
281 4 : if (failed('integer_dict_define')) return
282 :
283 4 : call integer_dict_define(ecapture_transitions_offset_dict, key, num_ecapture_transitions, ierr)
284 4 : if (failed('integer_dict_define')) return
285 :
286 24 : do j = 1, ntrans
287 8 : read(iounit,fmt=*,iostat=ierr) Si, Sf, logft
288 8 : if (ierr /= 0) then
289 0 : ierr = 0; exit
290 : end if
291 8 : num_ecapture_transitions = num_ecapture_transitions + 1
292 :
293 : ! pack the data into the array
294 8 : ecapture_transitions_data(num_ecapture_transitions, i_Si) = Si
295 8 : ecapture_transitions_data(num_ecapture_transitions, i_Sf) = Sf
296 :
297 12 : ecapture_logft_data(num_ecapture_transitions) = logft
298 :
299 : end do
300 :
301 : end do
302 :
303 2 : close(iounit)
304 :
305 2 : 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 2 : 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 2 : call integer_dict_create_hash(ecapture_reactions_dict, ierr)
318 2 : if (ierr /= 0) return
319 :
320 2 : call integer_dict_create_hash(ecapture_transitions_number_dict, ierr)
321 2 : if (ierr /= 0) return
322 :
323 2 : call integer_dict_create_hash(ecapture_transitions_offset_dict, ierr)
324 2 : if (ierr /= 0) return
325 :
326 2 : call realloc_integer2(ecapture_transitions_data, num_ecapture_transitions, num_transitions_data, ierr)
327 2 : if (ierr /= 0) return
328 :
329 2 : call realloc_double(ecapture_logft_data, num_ecapture_transitions, ierr)
330 2 : if (ierr /= 0) return
331 :
332 2 : call realloc_integer(ecapture_lhs_nuclide_id, num_ecapture_reactions, ierr)
333 2 : if (ierr /= 0) return
334 :
335 2 : call realloc_integer(ecapture_rhs_nuclide_id, num_ecapture_reactions, ierr)
336 12 : if (ierr /= 0) return
337 :
338 : contains
339 :
340 12 : logical function failed(str)
341 : character (len=*) :: str
342 12 : failed = (ierr /= 0)
343 12 : if (failed) then
344 0 : write(*,*) 'failed: ' // trim(str)
345 : end if
346 2 : end function failed
347 :
348 :
349 : end subroutine load_ecapture_transitions_list
350 :
351 :
352 :
353 : end module load_ecapture
|
