Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Ed Brown & 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 chem_isos_io
21 :
22 : use chem_def
23 : use math_lib
24 : use const_def, only: dp, amu, clight, five_thirds, mev_to_ergs
25 :
26 : implicit none
27 :
28 : contains
29 :
30 16 : subroutine do_read_chem_isos(isotopes_filename, ierr)
31 : use utils_lib
32 : use const_def, only: mesa_data_dir
33 : character (len=*), intent(in) :: isotopes_filename
34 : integer, intent(out) :: ierr
35 : integer :: i, j, k, iounit, pass, nvec, Z, N
36 : character (len=256) :: filename, buf
37 4112 : real(dp), target :: vec_ary(256)
38 : real(dp), pointer :: vec(:)
39 :
40 16 : ierr = 0
41 16 : vec => vec_ary
42 :
43 16 : filename = trim(mesa_data_dir) // '/chem_data/' // trim(isotopes_filename)
44 16 : num_chem_isos = 0
45 :
46 48 : do pass = 1, 2
47 :
48 32 : open(newunit=iounit, file=trim(filename), iostat=ierr, status='old',action='read')
49 32 : if ( ierr /= 0 ) then
50 0 : write(*,*) 'unable to open '// trim(filename)
51 0 : return
52 : end if
53 32 : read(iounit,'(A)') buf ! skip line 1
54 :
55 32 : if (pass == 1) then
56 :
57 125696 : do ! 4 lines per nuclide
58 125712 : read(iounit, *, iostat=ierr)
59 125712 : if (ierr /= 0) exit
60 125696 : read(iounit, *, iostat=ierr)
61 125696 : if (ierr /= 0) exit
62 125696 : read(iounit, *, iostat=ierr)
63 125696 : if (ierr /= 0) exit
64 125696 : read(iounit, *, iostat=ierr)
65 125696 : if (ierr /= 0) exit
66 125696 : num_chem_isos = num_chem_isos+1
67 : end do
68 16 : if (num_chem_isos == 0) then
69 0 : write (*,*) 'unable to retrieve isotopes from '//trim(filename)
70 0 : return
71 : end if
72 : else
73 :
74 16 : call allocate_nuclide_data(chem_isos, num_chem_isos, ierr)
75 16 : if (ierr /= 0) then
76 0 : write(*,*) 'unable to allocate nuclide data'
77 0 : return
78 : end if
79 :
80 125712 : do i = 1, num_chem_isos
81 :
82 125696 : read(iounit, '(A)',iostat=ierr) buf
83 125696 : if (ierr /= 0) exit
84 125696 : call parse_line(buf,i,ierr)
85 125696 : if (ierr /= 0) exit
86 :
87 502784 : do k=1,3
88 377088 : read(iounit,'(a)',iostat=ierr) buf
89 377088 : if (ierr == 0) then
90 377088 : call str_to_vector(buf, vec, nvec, ierr)
91 377088 : if (nvec /= 8) ierr = -1
92 : end if
93 377088 : if (ierr /= 0) exit
94 3519488 : do j=1,8
95 3393792 : chem_isos% pfcn(j+(k-1)*8, i) = vec(j)
96 : end do
97 : end do
98 125696 : if (ierr /= 0) exit
99 :
100 125696 : chem_isos% chem_id(i) = i
101 125696 : chem_isos% nuclide(i) = i
102 125712 : chem_isos% isomeric_state(i) = get_isomeric_state(chem_isos% name(i), ierr)
103 :
104 : end do
105 16 : if (ierr /= 0) then
106 0 : write (*,*) 'something went wrong in read of '//trim(filename)
107 0 : return
108 : end if
109 :
110 : end if
111 :
112 48 : close(iounit)
113 :
114 : end do
115 :
116 16 : if (ierr /= 0) return
117 :
118 16 : call do_create_nuclides_dict(chem_isos, chem_isos_dict, ierr)
119 16 : if (ierr /= 0) return
120 :
121 : !Set mass excess of proton and neutron
122 125712 : do i = 1, num_chem_isos
123 125696 : Z = chem_isos% Z(i)
124 125696 : N = chem_isos% N(i)
125 125696 : if(Z==1 .and. N==0) del_Mp=chem_isos% mass_excess(i)
126 125712 : if(N==1 .and. Z==0) del_Mn=chem_isos% mass_excess(i)
127 : end do
128 :
129 125712 : chem_isos% Z_plus_N = chem_isos% Z + chem_isos% N
130 :
131 : ! pre-calculate Z^5/3
132 125712 : do i = 1, num_chem_isos
133 125712 : chem_isos% Z53(i) = pow(real(chem_isos% Z(i), dp), five_thirds)
134 : end do
135 :
136 125712 : chem_isos% binding_energy = chem_isos% Z*del_Mp + chem_isos% N*del_Mn - chem_isos% mass_excess
137 :
138 : ! Recompute Atomic masses for double precision consistency.
139 125712 : chem_isos% W = chem_isos% Z_plus_N + chem_isos% mass_excess*(mev_to_ergs/amu/(clight*clight))
140 :
141 1824 : element_min_N = 99999
142 1824 : element_max_N = -1
143 125712 : do i = 1, num_chem_isos
144 125696 : Z = chem_isos% Z(i)
145 125696 : N = chem_isos% N(i)
146 125696 : if (N < element_min_N(Z)) element_min_N(Z) = N
147 125712 : if (N > element_max_N(Z)) element_max_N(Z) = N
148 : end do
149 :
150 : contains
151 :
152 125696 : integer function get_isomeric_state(name, ierr)
153 : character (len=*), intent(in) :: name
154 : integer, intent(out) :: ierr
155 : integer :: i, len
156 125696 : ierr = 0
157 125696 : get_isomeric_state = 0
158 125696 : len = len_trim(name)
159 720560 : do i=1,len
160 720560 : if (name(i:i) == '-') then
161 32 : read(name(i+1:len),*,iostat=ierr) get_isomeric_state
162 32 : if (ierr /= 0 .or. get_isomeric_state < 0 .or. get_isomeric_state > 99) then
163 0 : write(*,*) 'ERROR: invalid name for iso ' // trim(name) // ' in ' // trim(filename)
164 0 : return
165 : end if
166 : return
167 : end if
168 : end do
169 : end function get_isomeric_state
170 :
171 :
172 502784 : subroutine parse_line(line,i,ierr)
173 : character(len=*),intent(in) :: line
174 : integer, intent(in) :: i
175 : integer, intent(inout) :: ierr
176 : integer :: j,k
177 : integer, parameter :: num_cols=6
178 : character(len=256),dimension(num_cols) :: tmp
179 : character(len=256) :: tmp2
180 :
181 125696 : k=1
182 125696 : tmp2=''
183 879872 : tmp=''
184 6410496 : do j=1,len(line)
185 6410496 : if(line(j:j)==' ') cycle
186 4512400 : tmp2=trim(tmp2)//line(j:j)
187 4512400 : if(j<len(line))then
188 4512400 : if(line(j+1:j+1)==' ') then
189 754176 : tmp(k)=tmp2(1:len_trim(tmp2))
190 754176 : k=k+1
191 754176 : tmp2=''
192 : end if
193 : end if
194 4512400 : if(k==num_cols+1) exit
195 : end do
196 :
197 125696 : chem_isos% name(i)=tmp(1)
198 125696 : call str_to_double(tmp(2),chem_isos% W(i),ierr)
199 125696 : if (ierr /= 0) return
200 125696 : read(tmp(3),*) chem_isos% Z(i)
201 125696 : if (ierr /= 0) return
202 125696 : read(tmp(4),*) chem_isos% N(i)
203 125696 : if (ierr /= 0) return
204 125696 : call str_to_double(tmp(5),chem_isos% spin(i),ierr)
205 125696 : if (ierr /= 0) return
206 125696 : call str_to_double(tmp(6),chem_isos% mass_excess(i),ierr)
207 125696 : if (ierr /= 0) return
208 :
209 : end subroutine parse_line
210 :
211 : end subroutine do_read_chem_isos
212 :
213 :
214 32 : subroutine do_create_nuclides_dict(nuclides, nuclides_dict, ierr)
215 : use utils_lib, only: integer_dict_define, integer_dict_create_hash, integer_dict_lookup
216 : type(nuclide_data), intent(in) :: nuclides
217 : type (integer_dict), pointer :: nuclides_dict ! will be allocated
218 : integer, intent(out) :: ierr
219 : integer :: i
220 :
221 : ierr = 0
222 16 : nullify(nuclides_dict)
223 125712 : do i=1,nuclides% nnuclides
224 125696 : call integer_dict_define(nuclides_dict, nuclides% name(i), i, ierr)
225 125712 : if (ierr /= 0) return
226 : end do
227 :
228 16 : call integer_dict_create_hash(nuclides_dict, ierr)
229 16 : if (ierr /= 0) return
230 :
231 : end subroutine do_create_nuclides_dict
232 :
233 : end module chem_isos_io
|
