Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 adjust_xyz
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10
24 : use chem_def
25 : use utils_lib
26 :
27 : implicit none
28 :
29 : private
30 : public :: change_net
31 : public :: change_small_net
32 : public :: get_xa_for_standard_metals
33 : public :: get_xa_for_accretion
34 : public :: set_z
35 : public :: set_y
36 : public :: set_uniform_composition
37 : public :: set_standard_composition
38 : public :: set_uniform_xa_from_file
39 : public :: set_composition
40 : public :: do_change_to_xa_for_accretion
41 : public :: set_abundance_ratio
42 : public :: do_replace
43 : public :: do_set_abundance
44 : public :: do_uniform_mix_section
45 : public :: do_uniform_mix_envelope_down_to_t
46 :
47 : logical, parameter :: dbg = .false.
48 :
49 : contains
50 :
51 0 : subroutine change_net( &
52 : id, adjust_abundances_for_new_isos, new_net_name, ierr)
53 : integer, intent(in) :: id
54 : logical, intent(in) :: adjust_abundances_for_new_isos
55 : character (len=*), intent(in) :: new_net_name
56 : integer, intent(out) :: ierr
57 : call do_composition_fixup( &
58 0 : id, adjust_abundances_for_new_isos, new_net_name, ierr)
59 0 : end subroutine change_net
60 :
61 :
62 0 : subroutine change_small_net( &
63 : id, adjust_abundances_for_new_isos, new_net_name, ierr)
64 : integer, intent(in) :: id
65 : logical, intent(in) :: adjust_abundances_for_new_isos
66 : character (len=*), intent(in) :: new_net_name
67 : integer, intent(out) :: ierr
68 : call do_composition_fixup( &
69 0 : id, adjust_abundances_for_new_isos, new_net_name, ierr)
70 0 : end subroutine change_small_net
71 :
72 :
73 0 : subroutine do_composition_fixup( &
74 : id, adjust_abundances_for_new_isos, new_net_name, ierr)
75 : use net, only: do_micro_change_net
76 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
77 : integer, intent(in) :: id
78 : logical, intent(in) :: adjust_abundances_for_new_isos
79 : character (len=*), intent(in) :: new_net_name
80 : integer, intent(out) :: ierr
81 :
82 : type (star_info), pointer :: s
83 : integer :: old_num_isos, old_num_reactions
84 0 : integer :: old_net_iso(num_chem_isos)
85 0 : integer :: old_chem_id(num_chem_isos)
86 : integer :: species, nz, num_reactions
87 0 : integer, pointer :: chem_id(:)
88 : character (len=net_name_len) :: old_net_name
89 :
90 0 : call get_star_ptr(id, s, ierr)
91 0 : if (ierr /= 0) return
92 :
93 0 : if (new_net_name == s% net_name) return
94 :
95 : old_net_name = s% net_name
96 0 : old_num_isos = s% species
97 0 : old_num_reactions = s% num_reactions
98 0 : old_net_iso = s% net_iso
99 0 : old_chem_id(1:old_num_isos) = s% chem_id(1:old_num_isos)
100 :
101 0 : call do_micro_change_net(s, new_net_name, ierr)
102 0 : if (ierr /= 0) then
103 0 : s% retry_message = 'micro_change_net failed'
104 0 : if (s% report_ierr) write(*,*) s% retry_message
105 0 : return
106 : end if
107 :
108 0 : species = s% species
109 0 : nz = max(s% nz, s% prev_mesh_nz)
110 0 : chem_id => s% chem_id
111 0 : num_reactions = s% num_reactions
112 :
113 0 : write(*,*) 'change to "' // trim(new_net_name)//'"'
114 0 : write(*,*) 'number of species', s% species
115 :
116 0 : call realloc(s% dxdt_nuc); if (ierr /= 0) return
117 0 : call realloc(s% dxdt_nuc_start); if (ierr /= 0) return
118 0 : call realloc(s% d_epsnuc_dX); if (ierr /= 0) return
119 0 : call realloc(s% d_dxdt_nuc_dRho); if (ierr /= 0) return
120 0 : call realloc(s% d_dxdt_nuc_dT); if (ierr /= 0) return
121 :
122 0 : call realloc(s% dxdt_mix); if (ierr /= 0) return
123 :
124 0 : call realloc(s% d_eps_grav_dX); if (ierr /= 0) return
125 0 : call realloc(s% dlnE_dxa_for_partials); if (ierr /= 0) return
126 0 : call realloc(s% dlnPeos_dxa_for_partials); if (ierr /= 0) return
127 :
128 0 : call realloc(s% extra_diffusion_factor); if (ierr /= 0) return
129 0 : call realloc(s% edv); if (ierr /= 0) return
130 0 : call realloc(s% v_rad); if (ierr /= 0) return
131 0 : call realloc(s% g_rad); if (ierr /= 0) return
132 0 : call realloc(s% typical_charge); if (ierr /= 0) return
133 0 : call realloc(s% diffusion_dX); if (ierr /= 0) return
134 0 : call realloc(s% diffusion_D_self); if (ierr /= 0) return
135 :
136 0 : if (associated(s% d_dXdt_nuc_dX)) deallocate(s% d_dXdt_nuc_dX)
137 0 : allocate(s% d_dXdt_nuc_dX(species, species, nz + nz_alloc_extra), stat=ierr)
138 0 : if (ierr /= 0) return
139 :
140 0 : if (associated(s% d_eos_dxa)) deallocate(s% d_eos_dxa)
141 0 : allocate(s% d_eos_dxa(num_eos_d_dxa_results, species, nz + nz_alloc_extra), stat=ierr)
142 0 : if (ierr /= 0) return
143 :
144 0 : call realloc(s% xa_sub_xa_start); if (ierr /= 0) return
145 0 : call realloc(s% xa_start); if (ierr /= 0) return
146 0 : call realloc(s% prev_mesh_xa); if (ierr /= 0) return
147 0 : call do_xa(s% nz, s% xh, s% xa)
148 0 : if (s% generations > 1) call do_xa(s% nz_old, s% xh_old, s% xa_old)
149 :
150 0 : call realloc_reactions(s% raw_rate)
151 0 : if(ierr/=0) return
152 0 : call realloc_reactions(s% screened_rate)
153 0 : if(ierr/=0) return
154 0 : call realloc_reactions(s% eps_nuc_rate)
155 0 : if(ierr/=0) return
156 0 : call realloc_reactions(s% eps_neu_rate)
157 0 : if(ierr/=0) return
158 :
159 0 : s% need_to_setvars = .true.
160 0 : s% prev_mesh_species_or_nvar_hydro_changed = .true.
161 :
162 : contains
163 :
164 0 : subroutine do_xa(nz, xh, xa_startv)
165 0 : use net_lib, only: clean_up_fractions
166 : integer, intent(in) :: nz
167 : real(dp), pointer :: xh(:,:)
168 : real(dp), pointer :: xa_startv(:,:)
169 : real(dp), pointer :: xa_new(:,:)
170 : real(dp), parameter :: max_sum_abs = 10d0
171 : real(dp), parameter :: xsum_tol = 1d-2
172 0 : allocate(xa_new(species, nz + nz_alloc_extra), stat=ierr)
173 0 : if (ierr /= 0) return
174 : call set_x_new( &
175 : s, adjust_abundances_for_new_isos, &
176 : nz, old_num_isos, species, xh, xa_startv, xa_new, &
177 0 : old_chem_id, old_net_iso, chem_id, ierr)
178 0 : if (associated(xa_startv)) deallocate(xa_startv)
179 0 : xa_startv => xa_new
180 0 : end subroutine do_xa
181 :
182 0 : subroutine realloc(ptr)
183 : real(dp), pointer :: ptr(:, :)
184 0 : if (associated(ptr)) deallocate(ptr)
185 0 : allocate(ptr(species, nz + nz_alloc_extra), stat=ierr)
186 0 : end subroutine realloc
187 :
188 0 : subroutine realloc_reactions(ptr)
189 : real(dp), pointer :: ptr(:, :)
190 0 : if (associated(ptr)) deallocate(ptr)
191 0 : allocate(ptr(num_reactions, nz + nz_alloc_extra), stat=ierr)
192 0 : end subroutine realloc_reactions
193 :
194 : subroutine realloc_integer(ptr)
195 : integer, pointer :: ptr(:, :)
196 : if (associated(ptr)) deallocate(ptr)
197 : allocate(ptr(species, nz + nz_alloc_extra), stat=ierr)
198 : end subroutine realloc_integer
199 :
200 : end subroutine do_composition_fixup
201 :
202 :
203 0 : subroutine set_x_new( &
204 : s, adjust_abundances_for_new_isos, &
205 : nz, old_num_isos, species, xh, xa_startv, xa_new, &
206 0 : old_chem_id, old_net_iso, chem_id, ierr)
207 : use net_lib, only: clean1
208 : type (star_info), pointer :: s
209 : logical, intent(in) :: adjust_abundances_for_new_isos
210 : integer, intent(in) :: nz, old_num_isos, species
211 : real(dp), pointer :: xh(:,:)
212 : real(dp), pointer :: xa_startv(:,:)
213 : real(dp), pointer :: xa_new(:,:)
214 : integer, pointer :: chem_id(:)
215 : integer, intent(in) :: old_chem_id(old_num_isos), old_net_iso(num_chem_isos)
216 : integer, intent(out) :: ierr
217 :
218 : integer :: k, op_err
219 : real(dp), parameter :: max_sum_abs = 10
220 : real(dp), parameter :: xsum_tol = 1d-2
221 : include 'formats'
222 :
223 0 : ierr = 0
224 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
225 : do k=1, nz
226 : call set_new_abundances( &
227 : s, adjust_abundances_for_new_isos, &
228 : nz, old_num_isos, species, xh, xa_startv, xa_new, k, &
229 : old_chem_id, old_net_iso, chem_id, op_err)
230 : if (op_err /= 0) ierr = op_err
231 : call clean1(species, xa_new(1:species,k), max_sum_abs, xsum_tol, op_err)
232 : if (op_err /= 0) ierr = op_err
233 : end do
234 : !$OMP END PARALLEL DO
235 0 : if (ierr /= 0) then
236 0 : s% retry_message = 'set_new_abundances failed in set_x_new'
237 0 : if (s% report_ierr) write(*, *) s% retry_message
238 : return
239 : end if
240 :
241 0 : s% need_to_setvars = .true.
242 :
243 0 : end subroutine set_x_new
244 :
245 :
246 0 : subroutine set_new_abundances( &
247 : s, adjust_abundances_for_new_isos, nz, &
248 : old_num_isos, species, xh, xa_startv, xa_new, &
249 0 : k, old_chem_id, old_net_iso, chem_id, ierr)
250 0 : use chem_lib, only: chem_Xsol
251 : type (star_info), pointer :: s
252 : logical, intent(in) :: adjust_abundances_for_new_isos
253 : integer, intent(in) :: nz, old_num_isos, species, k
254 : real(dp), pointer :: xh(:,:)
255 : real(dp), pointer :: xa_startv(:,:)
256 : real(dp), pointer :: xa_new(:,:)
257 : integer, pointer :: chem_id(:)
258 : integer, intent(in) :: &
259 : old_chem_id(old_num_isos), old_net_iso(num_chem_isos)
260 : integer, intent(out) :: ierr
261 :
262 : real(dp) :: &
263 0 : total_neut, total_h, total_he, total_c, total_n, total_o, &
264 0 : total_ne, total_mg, total_si, total_s, total_ar, total_ca, &
265 0 : total_fe, total_co, other, &
266 : old_total_neut, old_total_h, old_total_he, old_total_c, old_total_n, old_total_o, &
267 0 : old_total_ne, old_total_mg, old_total_si, old_total_s, old_total_ar, old_total_ca, &
268 0 : old_total_fe, old_total_co, old_total_ni, old_other, &
269 0 : lgT, zfrac, xsol, lgT_lo, lgT_hi
270 : integer :: i, j, cid, Z
271 : character(len=solnamelen) :: sol_name
272 : logical :: dbg, did_total_neut, did_total_h, did_total_he, &
273 : did_total_c, did_total_n, did_total_o, did_total_ne, &
274 : did_total_mg, did_total_si, did_total_s, did_total_ar, did_total_ca, &
275 : did_total_fe, did_total_co, did_total_ni, did_total_other
276 :
277 : include 'formats'
278 :
279 0 : dbg = .false. !(k == nz) .or. (k == nz-1)
280 :
281 0 : ierr = 0
282 0 : zfrac = s% initial_z / zsol
283 0 : lgT_lo = s% lgT_lo_for_set_new_abundances
284 0 : lgT_hi = s% lgT_hi_for_set_new_abundances
285 :
286 : if (dbg) then
287 : write(*,'(A)')
288 : write(*,2) 'set_new_abundances', k
289 : write(*,2) 'old_num_isos', old_num_isos
290 : do j=1,old_num_isos
291 : write(*,2) 'old ' // chem_isos% name(old_chem_id(j)), k, xa_startv(j,k)
292 : end do
293 : write(*,'(A)')
294 : end if
295 :
296 0 : old_total_neut = 0
297 0 : old_total_h = 0
298 0 : old_total_he = 0
299 0 : old_total_c = 0
300 0 : old_total_n = 0
301 0 : old_total_o = 0
302 0 : old_total_ne = 0
303 0 : old_total_mg = 0
304 0 : old_total_si = 0
305 0 : old_total_s = 0
306 0 : old_total_ar = 0
307 0 : old_total_ca = 0
308 0 : old_total_fe = 0
309 0 : old_total_co = 0
310 0 : old_total_ni = 0
311 0 : old_other = 0
312 0 : do j=1, old_num_isos
313 0 : Z = chem_isos% Z(old_chem_id(j))
314 0 : if (Z == 0) then
315 0 : old_total_neut = old_total_neut + xa_startv(j,k)
316 : else if (Z == 1) then
317 0 : old_total_h = old_total_h + xa_startv(j,k)
318 : else if (Z == 2) then
319 0 : old_total_he = old_total_he + xa_startv(j,k)
320 : else if (Z == 6) then
321 0 : old_total_c = old_total_c + xa_startv(j,k)
322 : else if (Z == 7) then
323 0 : old_total_n = old_total_n + xa_startv(j,k)
324 : else if (Z == 8) then
325 0 : old_total_o = old_total_o + xa_startv(j,k)
326 : else if (Z == 10) then
327 0 : old_total_ne = old_total_ne + xa_startv(j,k)
328 : else if (Z == 12) then
329 0 : old_total_mg = old_total_mg + xa_startv(j,k)
330 : else if (Z == 14) then
331 0 : old_total_si = old_total_si + xa_startv(j,k)
332 : else if (Z == 16) then
333 0 : old_total_s = old_total_s + xa_startv(j,k)
334 : else if (Z == 18) then
335 0 : old_total_ar = old_total_ar + xa_startv(j,k)
336 : else if (Z == 20) then
337 0 : old_total_ca = old_total_ca + xa_startv(j,k)
338 : else if (Z == 26) then
339 0 : old_total_fe = old_total_fe + xa_startv(j,k)
340 : else if (Z == 27) then
341 0 : old_total_co = old_total_co + xa_startv(j,k)
342 : else if (Z == 28) then
343 0 : old_total_fe = old_total_fe + xa_startv(j,k)
344 : else
345 0 : old_other = old_other + xa_startv(j,k)
346 : end if
347 : end do
348 :
349 0 : lgT = get_test_lgT(k)
350 : ! copy old isos and set abundances for new ones
351 0 : do j=1, species
352 0 : cid = chem_id(j)
353 0 : i = old_net_iso(cid) ! old index number for this species
354 0 : if (i /= 0) then
355 0 : xa_new(j,k) = xa_startv(i,k)
356 0 : else if (.not. adjust_abundances_for_new_isos) then
357 0 : xa_new(j,k) = 0d0
358 : else
359 0 : if (chem_isos% Z(cid) > 5) then
360 0 : sol_name = chem_isos% name(cid)
361 0 : xsol = chem_Xsol(sol_name)
362 0 : xa_new(j,k) = zfrac*xsol
363 : if (dbg) then
364 : write(*,2) 'xa new ' // trim(sol_name), j, xa_new(j,k), xsol, zfrac
365 : end if
366 : else
367 : ! cannot simply add light elements to high temperature region
368 : ! or will create an explosive situation that won't converge
369 0 : if (lgT >= lgT_hi) then ! don't add any
370 0 : xa_new(j,k) = 0
371 : else
372 0 : sol_name = chem_isos% name(cid)
373 0 : xsol = chem_Xsol(sol_name)
374 0 : if (lgT > lgT_lo) then ! interpolate
375 : xa_new(j, k) = &
376 0 : xsol*0.5d0*(1+cospi((lgT-lgT_lo)/(lgT_hi-lgT_lo)))
377 : else ! add solar
378 0 : xa_new(j,k) = xsol
379 : end if
380 : end if
381 : end if
382 : end if
383 : end do
384 :
385 : total_neut = 0
386 : total_h = 0
387 : total_he = 0
388 : total_c = 0
389 : total_n = 0
390 : total_o = 0
391 : total_ne = 0
392 : total_mg = 0
393 : total_si = 0
394 : total_s = 0
395 : total_ar = 0
396 : total_ca = 0
397 : total_fe = 0
398 : total_co = 0
399 : total_fe = 0
400 : other = 0
401 0 : do j=1, species
402 0 : select case(int(chem_isos% Z(chem_id(j))))
403 : case (0)
404 0 : total_neut = total_neut + xa_new(j,k)
405 : case (1)
406 0 : total_h = total_h + xa_new(j,k)
407 : case (2)
408 0 : total_he = total_he + xa_new(j,k)
409 : case (6)
410 0 : total_c = total_c + xa_new(j,k)
411 : case (7)
412 0 : total_n = total_n + xa_new(j,k)
413 : case (8)
414 0 : total_o = total_o + xa_new(j,k)
415 : case (10)
416 0 : total_ne = total_ne + xa_new(j,k)
417 : case (12)
418 0 : total_mg = total_mg + xa_new(j,k)
419 : case (14)
420 0 : total_si = total_si + xa_new(j,k)
421 : case (16)
422 0 : total_s = total_s + xa_new(j,k)
423 : case (18)
424 0 : total_ar = total_ar + xa_new(j,k)
425 : case (20)
426 0 : total_ca = total_ca + xa_new(j,k)
427 : case (26)
428 0 : total_fe = total_fe + xa_new(j,k)
429 : case (27)
430 0 : total_co = total_co + xa_new(j,k)
431 : case (28)
432 0 : total_fe = total_fe + xa_new(j,k)
433 : case default
434 0 : other = other + xa_new(j,k)
435 : end select
436 : end do
437 :
438 : if (dbg) then
439 : write(*,1) 'old_total_n', old_total_n
440 : write(*,1) 'total_n', total_n
441 : write(*,1) 'old_total_n/total_n', old_total_n/total_n
442 : end if
443 :
444 0 : did_total_neut = .false.
445 0 : did_total_h = .false.
446 0 : did_total_he = .false.
447 0 : did_total_c = .false.
448 0 : did_total_n = .false.
449 0 : did_total_o = .false.
450 0 : did_total_ne = .false.
451 0 : did_total_mg = .false.
452 0 : did_total_si = .false.
453 0 : did_total_s = .false.
454 0 : did_total_ar = .false.
455 0 : did_total_ca = .false.
456 : did_total_fe = .false.
457 0 : did_total_co = .false.
458 0 : did_total_fe = .false.
459 0 : did_total_ni = .false.
460 0 : did_total_other = .false.
461 0 : do j=1, species
462 0 : select case(int(chem_isos% Z(chem_id(j))))
463 : case (0)
464 0 : if (total_neut > 0) then
465 0 : xa_new(j,k) = xa_new(j,k)*old_total_neut/total_neut
466 0 : did_total_neut = .true.
467 : end if
468 : case (1)
469 0 : if (total_h > 0) then
470 0 : xa_new(j,k) = xa_new(j,k)*old_total_h/total_h
471 0 : did_total_h = .true.
472 : end if
473 : case (2)
474 0 : if (total_he > 0) then
475 0 : xa_new(j,k) = xa_new(j,k)*old_total_he/total_he
476 0 : did_total_he = .true.
477 : end if
478 : case (6)
479 0 : if (total_c > 0) then
480 0 : xa_new(j,k) = xa_new(j,k)*old_total_c/total_c
481 0 : did_total_c = .true.
482 : end if
483 : case (7)
484 0 : if (total_n > 0) then
485 0 : xa_new(j,k) = xa_new(j,k)*old_total_n/total_n
486 0 : did_total_n = .true.
487 : end if
488 : case (8)
489 0 : if (total_o > 0) then
490 0 : xa_new(j,k) = xa_new(j,k)*old_total_o/total_o
491 0 : did_total_o = .true.
492 : end if
493 : case (10)
494 0 : if (total_ne > 0) then
495 0 : xa_new(j,k) = xa_new(j,k)*old_total_ne/total_ne
496 0 : did_total_ne = .true.
497 : end if
498 : case (12)
499 0 : if (total_mg > 0) then
500 0 : xa_new(j,k) = xa_new(j,k)*old_total_mg/total_mg
501 0 : did_total_mg = .true.
502 : end if
503 : case (14)
504 0 : if (total_si > 0) then
505 0 : xa_new(j,k) = xa_new(j,k)*old_total_si/total_si
506 0 : did_total_si = .true.
507 : end if
508 : case (16)
509 0 : if (total_s > 0) then
510 0 : xa_new(j,k) = xa_new(j,k)*old_total_s/total_s
511 0 : did_total_s = .true.
512 : end if
513 : case (18)
514 0 : if (total_ar > 0) then
515 0 : xa_new(j,k) = xa_new(j,k)*old_total_ar/total_ar
516 0 : did_total_ar = .true.
517 : end if
518 : case (20)
519 0 : if (total_ca > 0) then
520 0 : xa_new(j,k) = xa_new(j,k)*old_total_ca/total_ca
521 0 : did_total_ca = .true.
522 : end if
523 : case (26)
524 0 : if (total_fe > 0) then
525 0 : xa_new(j,k) = xa_new(j,k)*old_total_fe/total_fe
526 0 : did_total_fe = .true.
527 : end if
528 : case (27)
529 0 : if (total_co > 0) then
530 0 : xa_new(j,k) = xa_new(j,k)*old_total_co/total_co
531 0 : did_total_co = .true.
532 : end if
533 : case (28)
534 0 : if (total_fe > 0) then
535 0 : xa_new(j,k) = xa_new(j,k)*old_total_fe/total_fe
536 0 : did_total_fe = .true.
537 : end if
538 : case default
539 0 : if (other > 0) then
540 0 : xa_new(j,k) = xa_new(j,k)*old_other/other
541 0 : did_total_other = .true.
542 : end if
543 : end select
544 : end do
545 :
546 : ! check for leftovers and dump them into the last iso
547 0 : j = species
548 0 : if (old_total_neut > 0 .and. .not. did_total_neut) &
549 0 : xa_new(j,k) = xa_new(j,k) + old_total_neut
550 0 : if (old_total_h > 0 .and. .not. did_total_h) &
551 0 : xa_new(j,k) = xa_new(j,k) + old_total_h
552 0 : if (old_total_he > 0 .and. .not. did_total_he) &
553 0 : xa_new(j,k) = xa_new(j,k) + old_total_he
554 0 : if (old_total_c > 0 .and. .not. did_total_c) &
555 0 : xa_new(j,k) = xa_new(j,k) + old_total_c
556 0 : if (old_total_n > 0 .and. .not. did_total_n) &
557 0 : xa_new(j,k) = xa_new(j,k) + old_total_n
558 0 : if (old_total_o > 0 .and. .not. did_total_o) &
559 0 : xa_new(j,k) = xa_new(j,k) + old_total_o
560 0 : if (old_total_ne > 0 .and. .not. did_total_ne) &
561 0 : xa_new(j,k) = xa_new(j,k) + old_total_ne
562 0 : if (old_total_mg > 0 .and. .not. did_total_mg) &
563 0 : xa_new(j,k) = xa_new(j,k) + old_total_mg
564 0 : if (old_total_si > 0 .and. .not. did_total_si) &
565 0 : xa_new(j,k) = xa_new(j,k) + old_total_si
566 :
567 0 : if (old_total_s > 0 .and. .not. did_total_s) &
568 0 : xa_new(j,k) = xa_new(j,k) + old_total_s
569 0 : if (old_total_ar > 0 .and. .not. did_total_ar) &
570 0 : xa_new(j,k) = xa_new(j,k) + old_total_ar
571 0 : if (old_total_ca > 0 .and. .not. did_total_ca) &
572 0 : xa_new(j,k) = xa_new(j,k) + old_total_ca
573 0 : if (old_total_fe > 0 .and. .not. did_total_fe) &
574 0 : xa_new(j,k) = xa_new(j,k) + old_total_fe
575 0 : if (old_total_co > 0 .and. .not. did_total_co) &
576 0 : xa_new(j,k) = xa_new(j,k) + old_total_co
577 : if (old_total_ni > 0 .and. .not. did_total_ni) &
578 : xa_new(j,k) = xa_new(j,k) + old_total_ni
579 0 : if (old_other > 0 .and. .not. did_total_other) &
580 0 : xa_new(j,k) = xa_new(j,k) + old_other
581 :
582 0 : if (abs(sum(xa_new(:,k)) - 1d0) > 0.1d0) then
583 0 : !$omp critical (new_abund)
584 0 : write(*,'(A)')
585 0 : write(*,2) 'bad sum: set_new_abundances', k, sum(xa_new(:,k))
586 0 : write(*,2) 'species', species
587 0 : do j=1,species
588 0 : write(*,2) trim(chem_isos% name(chem_id(j))), k, xa_new(j,k)
589 : end do
590 0 : write(*,'(A)')
591 0 : write(*,2) 'old_total_neut', k, old_total_neut
592 0 : write(*,2) 'old_total_h', k, old_total_h
593 0 : write(*,2) 'old_total_he', k, old_total_he
594 0 : write(*,2) 'old_total_c', k, old_total_c
595 0 : write(*,2) 'old_total_n', k, old_total_n
596 0 : write(*,2) 'old_total_o', k, old_total_o
597 0 : write(*,2) 'old_other', k, old_other
598 0 : write(*,'(A)')
599 0 : call mesa_error(__FILE__,__LINE__,'debug: set_new_abundances')
600 : !$omp end critical (new_abund)
601 :
602 : return
603 : end if
604 :
605 0 : if (dbg) then
606 : write(*,'(A)')
607 : write(*,2) 'new num species', species
608 : do j=1,species
609 : write(*,2) trim(chem_isos% name(chem_id(j))), j, xa_new(j,k)
610 : end do
611 : write(*,'(A)')
612 : call mesa_error(__FILE__,__LINE__,'debug: set_new_abundances')
613 : end if
614 :
615 : contains
616 :
617 0 : real(dp) function get_test_lgT(loc)
618 0 : use star_utils, only: get_lnT_from_xh
619 : integer, intent(in) :: loc
620 : integer :: k, kmax
621 0 : get_test_lgT = get_lnT_from_xh(s,loc)/ln10
622 0 : kmax = min(min(s%nz,nz),size(s% mlt_D,dim=1))
623 0 : do k=loc+1, kmax
624 0 : if (s% mlt_D(k) < 1d2) return
625 0 : get_test_lgT = get_lnT_from_xh(s,k)/ln10
626 : end do
627 0 : end function get_test_lgT
628 :
629 : end subroutine set_new_abundances
630 :
631 : subroutine get_ag89_composition(s, species, xa, ierr)
632 : use chem_lib, only: chem_Xsol
633 : type (star_info), pointer :: s
634 : integer, intent(in) :: species
635 : real(dp) :: xa(species)
636 : integer, intent(out) :: ierr
637 : integer :: i, mg24
638 : ierr = 0
639 : do i=1, species
640 : xa(i) = chem_Xsol(chem_isos% name(s% chem_id(i)))
641 : end do
642 : mg24 = s% net_iso(img24)
643 : xa(mg24) = xa(mg24) + (1 - sum(xa(1:species)))
644 : end subroutine get_ag89_composition
645 :
646 :
647 0 : subroutine set_y(s, y, nzlo, nzhi, ierr)
648 : use star_utils, only: eval_current_abundance
649 : type (star_info), pointer :: s
650 : integer, intent(in) :: nzlo, nzhi
651 : real(dp), intent(in) :: y
652 : integer, intent(out) :: ierr
653 :
654 0 : real(dp) :: xh1, xhe3, xhe4, z, ratio, desired_xh1, desired_xhe4
655 : include 'formats'
656 :
657 : ierr = 0
658 :
659 0 : if (y == 0) then ! convert all he4 to h1
660 0 : call set_abundance_ratio(s% id, ihe4, ih1, 0d0, nzlo, nzhi, ierr)
661 : end if
662 :
663 0 : xh1 = eval_current_abundance(s, s% net_iso(ih1), nzlo, nzhi, ierr)
664 0 : xhe3 = eval_current_abundance(s, s% net_iso(ihe3), nzlo, nzhi, ierr)
665 0 : xhe4 = eval_current_abundance(s, s% net_iso(ihe4), nzlo, nzhi, ierr)
666 0 : z = 1d0 - (xh1 + xhe3 + xhe4)
667 : ! keep xhe3 and z constant; change ratio of xh1 and xhe4
668 0 : desired_xhe4 = y - xhe3
669 0 : desired_xh1 = 1d0 - z - y
670 :
671 0 : if (desired_xh1 <= 0d0) then ! convert all h1 to he4
672 0 : ratio = 0d0
673 0 : else if (desired_xhe4 <= 0d0) then ! convert all he4 to h1
674 0 : ratio = 1d0
675 : else
676 0 : ratio = desired_xh1/desired_xhe4
677 : end if
678 :
679 0 : call set_abundance_ratio(s% id, ih1, ihe4, ratio, nzlo, nzhi, ierr)
680 :
681 0 : end subroutine set_y
682 :
683 :
684 0 : subroutine set_abundance_ratio(id, i1, i2, ratio, nzlo, nzhi, ierr)
685 : integer, intent(in) :: id, i1, i2, nzlo, nzhi
686 : real(dp), intent(in) :: ratio ! new x(i1)/x(i2)
687 : integer, intent(out) :: ierr
688 :
689 : type (star_info), pointer :: s
690 : integer :: k, j1, j2
691 0 : real(dp) :: xsum
692 :
693 : ierr = 0
694 :
695 0 : call get_star_ptr(id, s, ierr)
696 0 : if (ierr /= 0) return
697 :
698 0 : j1 = s% net_iso(i1)
699 0 : j2 = s% net_iso(i2)
700 0 : if (j1 == 0 .or. j2 == 0) then
701 0 : ierr = -1
702 0 : return
703 : end if
704 :
705 0 : do k=nzlo, nzhi
706 0 : xsum = s% xa(j1, k) + s% xa(j2, k)
707 0 : s% xa(j2, k) = xsum/(1+ratio)
708 0 : s% xa(j1, k) = xsum - s% xa(j2, k)
709 : end do
710 :
711 0 : s% need_to_setvars = .true.
712 :
713 0 : end subroutine set_abundance_ratio
714 :
715 :
716 0 : subroutine read_xa(s, species, xa, filename, ierr)
717 : use chem_def
718 : use chem_lib
719 : use utils_def
720 : use utils_lib
721 : type (star_info), pointer :: s
722 : integer, intent(in) :: species
723 : real(dp) :: xa(species)
724 : character (len=*), intent(in) :: filename
725 : integer, intent(out) :: ierr
726 :
727 : integer :: iounit, t, n, i, cid, j
728 : character (len=strlen) :: buffer, string
729 : logical, parameter :: dbg = .false.
730 :
731 0 : ierr = 0
732 :
733 0 : xa(1:species) = 0 ! default
734 :
735 0 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
736 0 : if (ierr /= 0) then
737 0 : write(*,*) 'failed to open abundance specification file "' // trim(filename)//'"'
738 0 : return
739 : end if
740 :
741 0 : n = 0
742 0 : i = 0
743 :
744 : do
745 0 : t = token(iounit, n, i, buffer, string)
746 0 : select case(t)
747 : case(name_token) ! iso name
748 0 : cid = chem_get_iso_id(string)
749 0 : if (cid <= 0) then
750 0 : write(*,*) 'reading ' // trim(filename)
751 0 : write(*,*) 'unknown chem name ' // trim(string)
752 0 : ierr = -1
753 0 : call cleanup
754 0 : return
755 : end if
756 0 : j = s% net_iso(cid)
757 0 : if (j == 0) then
758 0 : write(*,*) 'reading ' // trim(filename)
759 0 : write(*,*) 'iso not in current net ' // trim(string)
760 0 : ierr = -1
761 0 : call cleanup
762 0 : return
763 : end if
764 0 : t = token(iounit, n, i, buffer, string)
765 0 : if (t /= name_token) then
766 0 : write(*,*) 'reading ' // trim(filename)
767 : write(*,*) 'failed in reading abundance value for ' // &
768 0 : trim(chem_isos% name(cid))
769 0 : ierr = -1
770 0 : call cleanup
771 0 : return
772 : end if
773 :
774 : !read(string,fmt=*,iostat=ierr) xa(j)
775 0 : call str_to_double(string, xa(j), ierr)
776 0 : if (ierr /= 0 .or. xa(j) > 1 .or. xa(j) < 0) then
777 0 : write(*,*) 'reading ' // trim(filename)
778 : write(*,*) 'invalid number given as abundance value for ' // &
779 0 : trim(chem_isos% name(cid))
780 0 : ierr = -1
781 0 : call cleanup
782 0 : return
783 : end if
784 : case(eof_token)
785 0 : exit
786 : case default
787 0 : write(*,*) 'error in reading abundance file ' // trim(filename)
788 0 : ierr = -1
789 0 : call cleanup
790 0 : return
791 : end select
792 :
793 : end do
794 :
795 0 : call cleanup
796 :
797 0 : if (abs(sum(xa(1:species)) - 1d0) > 1d-6) then
798 0 : write(*,*) 'reading ' // trim(filename)
799 0 : write(*,*) 'abundances fail to add to 1.0: sum - 1 = ', sum(xa(:)) - 1d0
800 0 : ierr = -1
801 0 : return
802 : end if
803 :
804 0 : xa(1:species) = xa(1:species)/sum(xa)
805 :
806 : contains
807 :
808 0 : subroutine cleanup
809 0 : close(iounit)
810 0 : end subroutine cleanup
811 :
812 : end subroutine read_xa
813 :
814 :
815 0 : subroutine set_uniform_xa_from_file(id, file_for_uniform_xa, ierr)
816 : integer, intent(in) :: id
817 : character (len=*), intent(in) :: file_for_uniform_xa
818 : integer, intent(out) :: ierr
819 : type (star_info), pointer :: s
820 0 : call get_star_ptr(id, s, ierr)
821 0 : if (ierr /= 0) return
822 0 : call doit(s% species)
823 :
824 : contains
825 :
826 0 : subroutine doit(species)
827 : integer, intent(in) :: species
828 0 : real(dp) :: xa(species)
829 0 : call read_xa(s, species, xa, file_for_uniform_xa, ierr)
830 0 : if (ierr /= 0) return
831 0 : call set_uniform_composition(id, species, xa, ierr)
832 0 : end subroutine doit
833 :
834 : end subroutine set_uniform_xa_from_file
835 :
836 :
837 0 : subroutine set_uniform_composition(id, species, xa, ierr)
838 : integer, intent(in) :: species
839 : real(dp), intent(in) :: xa(species)
840 : integer, intent(in) :: id
841 : integer, intent(out) :: ierr
842 : type (star_info), pointer :: s
843 0 : call get_star_ptr(id, s, ierr)
844 0 : if (ierr /= 0) return
845 0 : call set_composition(id, 1, s% nz, species, xa, ierr)
846 : end subroutine set_uniform_composition
847 :
848 :
849 0 : subroutine set_composition(id, nzlo, nzhi, num_species, xa_new, ierr)
850 : integer, intent(in) :: nzlo, nzhi, num_species
851 : real(dp), intent(in) :: xa_new(num_species)
852 : integer, intent(in) :: id
853 : integer, intent(out) :: ierr
854 : integer :: j, k, nz, species
855 : type (star_info), pointer :: s
856 : include 'formats'
857 0 : call get_star_ptr(id, s, ierr)
858 0 : if (ierr /= 0) return
859 0 : nz = s% nz
860 0 : species = s% species
861 0 : if (num_species /= species) then
862 0 : ierr = -1
863 0 : s% retry_message = 'set_composition requires number of species to match current model'
864 0 : if (s% report_ierr) write(*, *) s% retry_message
865 0 : return
866 : end if
867 0 : if (nzlo < 1 .or. nzhi > nz) then
868 0 : ierr = -1
869 0 : s% retry_message = 'set_composition requires nzlo and nzhi to be within 1 to current num zones'
870 0 : if (s% report_ierr) write(*, *) s% retry_message
871 0 : return
872 : end if
873 0 : if (abs(1d0-sum(xa_new(1:species))) > 1d-6) then
874 0 : ierr = -1
875 0 : s% retry_message = 'set_composition requires new mass fractions to add to 1'
876 0 : if (s% report_ierr) write(*, *) s% retry_message
877 0 : return
878 : end if
879 0 : do k=nzlo,nzhi
880 0 : do j=1,species
881 0 : s% xa(j,k) = xa_new(j)
882 : end do
883 : end do
884 0 : ierr = 0
885 :
886 0 : s% need_to_setvars = .true.
887 :
888 : end subroutine set_composition
889 :
890 :
891 0 : subroutine set_standard_composition( &
892 : s, species, h1, h2, he3, he4, which_zfracs, &
893 : dump_missing_into_heaviest, ierr)
894 : type (star_info), pointer :: s
895 : integer, intent(in) :: species
896 : real(dp), intent(in) :: h1, h2, he3, he4 ! mass fractions
897 : integer, intent(in) :: which_zfracs ! defined in chem_def. e.g., GS98_zfracs
898 : logical, intent(in) :: dump_missing_into_heaviest
899 : integer, intent(out) :: ierr
900 0 : real(dp) :: xa(species)
901 : ierr = 0
902 : call get_xa_for_standard_metals(s, &
903 : species, s% chem_id, s% net_iso, h1, h2, he3, he4, which_zfracs, &
904 0 : dump_missing_into_heaviest, xa, ierr)
905 0 : if (ierr /= 0) return
906 0 : call set_uniform_composition(s% id, species, xa, ierr)
907 0 : end subroutine set_standard_composition
908 :
909 :
910 0 : subroutine get_xa_for_accretion(s, xa, ierr)
911 : use chem_lib, only: chem_get_iso_id
912 : type (star_info), pointer :: s
913 : real(dp) :: xa(:) ! (species)
914 : integer, intent(out) :: ierr
915 : integer :: j, i, species, cid
916 : include 'formats'
917 0 : ierr = 0
918 0 : species = s% species
919 0 : if (s% accrete_given_mass_fractions) then
920 0 : xa(1:species) = 0
921 0 : do j=1,s% num_accretion_species
922 0 : if (len_trim(s% accretion_species_id(j)) == 0) cycle
923 0 : cid = chem_get_iso_id(s% accretion_species_id(j))
924 0 : if (cid <= 0) cycle
925 0 : i = s% net_iso(cid)
926 0 : if (i == 0) cycle
927 0 : xa(i) = s% accretion_species_xa(j)
928 : end do
929 0 : if (abs(1d0 - sum(xa(1:species))) > 1d-2) then
930 : write(*,'(a)') &
931 0 : 'get_xa_for_accretion: accretion species mass fractions do not add to 1'
932 0 : write(*,1) 'sum(xa(1:species))', sum(xa(1:species))
933 0 : do j=1,s% num_accretion_species
934 0 : write(*,2) trim(s% accretion_species_id(j)), j, xa(j)
935 : end do
936 0 : ierr = -1
937 0 : return
938 : end if
939 0 : xa(1:species) = xa(1:species)/sum(xa(1:species))
940 : return
941 : end if
942 : call get_xa_for_standard_metals(s, &
943 : s% species, s% chem_id, s% net_iso, &
944 : s% accretion_h1, s% accretion_h2, s% accretion_he3, s% accretion_he4, &
945 : s% accretion_zfracs, &
946 0 : s% accretion_dump_missing_metals_into_heaviest, xa, ierr)
947 0 : end subroutine get_xa_for_accretion
948 :
949 :
950 0 : subroutine do_change_to_xa_for_accretion(id, nzlo_in, nzhi_in, ierr)
951 : integer, intent(in) :: id
952 : integer, intent(in) :: nzlo_in, nzhi_in
953 : integer, intent(out) :: ierr
954 : integer :: j, k, species
955 0 : real(dp), pointer :: xa(:) ! (species)
956 : type (star_info), pointer :: s
957 : integer :: nzlo, nzhi
958 : include 'formats'
959 : ierr = 0
960 0 : call get_star_ptr(id, s, ierr)
961 0 : if (ierr /= 0) return
962 0 : species = s% species
963 0 : nzlo = nzlo_in
964 0 : nzhi = nzhi_in
965 0 : if (nzlo < 1) nzlo = 1
966 0 : if (nzhi < 1 ) nzhi = s% nz
967 0 : allocate(xa(species))
968 0 : if (s% accrete_same_as_surface) then
969 0 : do j=1,species
970 0 : xa(j) = s% xa(j,1)
971 : end do
972 : else
973 0 : call get_xa_for_accretion(s, xa, ierr)
974 0 : if (ierr /= 0) then
975 0 : s% retry_message = 'get_xa_for_accretion failed in change_to_xa_for_accretion'
976 0 : if (s% report_ierr) write(*, *) s% retry_message
977 0 : deallocate(xa)
978 0 : return
979 : end if
980 : end if
981 0 : do k=nzlo,min(s% nz,nzhi)
982 0 : do j=1,species
983 0 : s% xa(j,k) = xa(j)
984 : end do
985 : end do
986 0 : deallocate(xa)
987 0 : s% need_to_setvars = .true.
988 0 : end subroutine do_change_to_xa_for_accretion
989 :
990 :
991 0 : subroutine get_xa_for_standard_metals( &
992 0 : s, species, chem_id, net_iso, h1_in, h2_in, he3_in, he4_in, which_zfracs, &
993 0 : dump_missing_into_heaviest, xa, ierr)
994 : use chem_def
995 : type (star_info), pointer :: s
996 : integer, intent(in) :: species, chem_id(:), net_iso(:), which_zfracs
997 : real(dp), intent(in) :: h1_in, h2_in, he3_in, he4_in
998 : logical, intent(in) :: dump_missing_into_heaviest
999 : real(dp) :: xa(:) ! (species)
1000 : integer, intent(out) :: ierr
1001 0 : real(dp) :: zfrac(num_chem_elements), Z, h1, h2, he3, he4
1002 : integer :: i
1003 : include 'formats'
1004 0 : ierr = 0
1005 0 : h1 = h1_in; h2 = h2_in; he3 = he3_in; he4 = he4_in
1006 0 : select case(which_zfracs)
1007 : case (AG89_zfracs)
1008 0 : zfrac(:) = AG89_element_zfrac(:)
1009 : case (GN93_zfracs)
1010 0 : zfrac(:) = GN93_element_zfrac(:)
1011 : case (GS98_zfracs)
1012 0 : zfrac(:) = GS98_element_zfrac(:)
1013 : case (L03_zfracs)
1014 0 : zfrac(:) = L03_element_zfrac(:)
1015 : case (AGS05_zfracs)
1016 0 : zfrac(:) = AGS05_element_zfrac(:)
1017 : case (AGSS09_zfracs)
1018 0 : zfrac(:) = AGSS09_element_zfrac(:)
1019 : case (L09_zfracs)
1020 0 : zfrac(:) = L09_element_zfrac(:)
1021 : case (A09_Prz_zfracs)
1022 0 : zfrac(:) = A09_Prz_zfrac(:)
1023 : case (MB22_photospheric_zfracs)
1024 0 : zfrac(:) = MB22_photospheric_element_zfrac(:)
1025 : case (AAG21_photospheric_zfracs)
1026 0 : zfrac(:) = AAG21_photospheric_element_zfrac(:)
1027 : case (Custom_zfracs) ! use non-standard values given in controls
1028 0 : zfrac(:) = 0
1029 0 : zfrac(e_li) = s% z_fraction_li
1030 0 : zfrac(e_be) = s% z_fraction_be
1031 0 : zfrac(e_b) = s% z_fraction_b
1032 0 : zfrac(e_c) = s% z_fraction_c
1033 0 : zfrac(e_n) = s% z_fraction_n
1034 0 : zfrac(e_o) = s% z_fraction_o
1035 0 : zfrac(e_f) = s% z_fraction_f
1036 0 : zfrac(e_ne) = s% z_fraction_ne
1037 0 : zfrac(e_na) = s% z_fraction_na
1038 0 : zfrac(e_mg) = s% z_fraction_mg
1039 0 : zfrac(e_al) = s% z_fraction_al
1040 0 : zfrac(e_si) = s% z_fraction_si
1041 0 : zfrac(e_p) = s% z_fraction_p
1042 0 : zfrac(e_s) = s% z_fraction_s
1043 0 : zfrac(e_cl) = s% z_fraction_cl
1044 0 : zfrac(e_ar) = s% z_fraction_ar
1045 0 : zfrac(e_k) = s% z_fraction_k
1046 0 : zfrac(e_ca) = s% z_fraction_ca
1047 0 : zfrac(e_sc) = s% z_fraction_sc
1048 0 : zfrac(e_ti) = s% z_fraction_ti
1049 0 : zfrac(e_v) = s% z_fraction_v
1050 0 : zfrac(e_cr) = s% z_fraction_cr
1051 0 : zfrac(e_mn) = s% z_fraction_mn
1052 0 : zfrac(e_fe) = s% z_fraction_fe
1053 0 : zfrac(e_co) = s% z_fraction_co
1054 0 : zfrac(e_ni) = s% z_fraction_ni
1055 0 : zfrac(e_cu) = s% z_fraction_cu
1056 0 : zfrac(e_zn) = s% z_fraction_zn
1057 0 : if (abs(sum(zfrac)-1) > 1d-6) then
1058 0 : write(*,*) 'bad sum(zfrac) for specified z fractions', sum(zfrac)
1059 0 : ierr = -1
1060 : end if
1061 0 : do i = 1, size(zfrac, dim=1)
1062 0 : if (zfrac(i) < 0) then
1063 0 : write(*,2) 'zfrac(i)', i, zfrac(i)
1064 0 : ierr = -1
1065 0 : return
1066 : end if
1067 : end do
1068 0 : zfrac(:) = zfrac(:) / sum(zfrac(:))
1069 : case default
1070 0 : if (abs(1d0 - (h1 + h2 + he3 + he4)) < 1d-8) then ! okay -- no metals
1071 0 : if (h1 > he4) then
1072 0 : h1 = max(0d0, min(1d0, 1d0 - (h2 + he3 + he4)))
1073 : else
1074 0 : he4 = max(0d0, min(1d0, 1d0 - (h1 + h2 + he3)))
1075 : end if
1076 0 : zfrac(:) = 0
1077 : else
1078 0 : ierr = -1
1079 : end if
1080 : end select
1081 0 : if (ierr /= 0) return
1082 0 : Z = max(0d0, min(1d0, 1d0 - (h1 + h2 + he3 + he4)))
1083 : call get_xa( &
1084 : species, chem_id, net_iso, h1, h2, he3, he4, zfrac, Z, &
1085 0 : dump_missing_into_heaviest, xa, ierr)
1086 0 : do i=1,species
1087 0 : if (xa(i) < 0 .or. is_bad(xa(i))) then
1088 0 : write(*,2) 'get_xa_for_standard_metals xa(i)', i, xa(i)
1089 0 : ierr = -1
1090 0 : return
1091 : end if
1092 : end do
1093 0 : end subroutine get_xa_for_standard_metals
1094 :
1095 :
1096 0 : subroutine get_xa( &
1097 0 : species, chem_id, net_iso, h1, h2, he3, he4, zfrac, Zinit, &
1098 0 : dump_missing_into_heaviest, xa, ierr)
1099 0 : use chem_def
1100 : integer, intent(in) :: species, chem_id(:), net_iso(:)
1101 : real(dp), intent(in) :: h1, h2, he3, he4, Zinit
1102 : real(dp), intent(in) :: zfrac(:) ! (num_chem_elements)
1103 : logical, intent(in) :: dump_missing_into_heaviest
1104 : real(dp) :: xa(:) ! (species)
1105 : real(dp), parameter :: tiny = 1d-15
1106 : integer, intent(out) :: ierr
1107 : include 'formats'
1108 : if (h1 < 0 .and. abs(h1)>tiny &
1109 : .or. h2 < 0 .and. abs(h2)>tiny &
1110 : .or. he3 < 0.and. abs(he3)>tiny &
1111 : .or. he4 < 0.and. abs(he4)>tiny &
1112 0 : .or.Zinit < 0.and.abs(Zinit)>tiny ) then
1113 0 : ierr = -1
1114 0 : write(*,*) 'when setting composition, need to provide values for h1, h2, he3, and he4'
1115 0 : write(*,*) 'H1=', h1
1116 0 : write(*,*) 'H2=', h2
1117 0 : write(*,*) 'He3=', he3
1118 0 : write(*,*) 'He4=', he4
1119 0 : write(*,*) 'Z =', Zinit
1120 0 : return
1121 0 : else if ( abs(1d0-(h1+h2+he3+he4+Zinit)) > tiny ) then
1122 0 : ierr = -2
1123 0 : write(*,1) 'h1', h1
1124 0 : write(*,1) 'h2', h2
1125 0 : write(*,1) 'he3', he3
1126 0 : write(*,1) 'he4', he4
1127 0 : write(*,1) 'Zinit', Zinit
1128 0 : write(*,*) 'sum of (H1+H2+He3+He4+Z) does not sum to 1: (1-sum)= ', &
1129 0 : 1d0-(h1+h2+he3+he4+Zinit)
1130 0 : call mesa_error(__FILE__,__LINE__,'get_xa')
1131 : end if
1132 0 : ierr = 0
1133 0 : xa(:) = 0
1134 0 : if (net_iso(ih2) /= 0) then
1135 0 : xa(net_iso(ih2)) = h2
1136 0 : xa(net_iso(ih1)) = h1
1137 0 : else if (net_iso(ih1) /= 0) then
1138 0 : xa(net_iso(ih1)) = h1 + h2
1139 : else
1140 0 : ierr = -1
1141 0 : write(*,*) 'require h1 to be in net'
1142 0 : return
1143 : end if
1144 0 : if (net_iso(ihe3) /= 0) xa(net_iso(ihe3)) = he3
1145 0 : if (net_iso(ihe4) /= 0) xa(net_iso(ihe4)) = he4
1146 : call adjust_z_fractions( &
1147 : species, chem_id, 1d0-(h1+h2+he3+he4), zfrac, &
1148 0 : dump_missing_into_heaviest, xa, ierr)
1149 0 : xa(:) = xa(:)/sum(xa)
1150 0 : end subroutine get_xa
1151 :
1152 :
1153 0 : subroutine adjust_z_fractions( &
1154 0 : species, chem_id, ztotal, zfrac, dump_missing_into_heaviest, xa, ierr)
1155 0 : use chem_def
1156 : use chem_lib, only: lodders03_element_atom_percent
1157 : integer, intent(in) :: species, chem_id(:) ! (species)
1158 : real(dp), intent(in) :: ztotal
1159 : real(dp), intent(in) :: zfrac(:) ! (num_chem_elements)
1160 : logical, intent(in) :: dump_missing_into_heaviest
1161 : real(dp) :: xa(:) ! (species) h & he not touched
1162 : integer, intent(out) :: ierr
1163 :
1164 0 : integer :: i, j, cid, cid2, element_id, jskip, Z(species)
1165 0 : real(dp) :: new_xa(species), frac, frac_sum, iso_frac, zsum
1166 :
1167 : include 'formats'
1168 :
1169 0 : ierr = 0
1170 0 : new_xa(:) = 0
1171 :
1172 0 : if (dump_missing_into_heaviest) then
1173 0 : jskip = maxloc(chem_isos% W(chem_id(:)),dim=1) ! find the heaviest
1174 : else
1175 : jskip = 0
1176 : end if
1177 :
1178 0 : do j=1,species
1179 0 : Z(j) = chem_isos% Z(chem_id(j))
1180 : end do
1181 :
1182 0 : do j=1,species
1183 0 : if (j == jskip) cycle
1184 0 : if (Z(j) <= 2) cycle
1185 0 : cid = chem_id(j)
1186 : ! multiply lodders' number fractions by chem_isos% A since want fractions by mass
1187 0 : frac = 1d-2*lodders03_element_atom_percent(chem_isos% name(cid))*chem_isos% W(cid)
1188 0 : if (frac == 0) cycle
1189 : frac_sum = 0
1190 0 : do i=1,species
1191 0 : if (Z(i) /= Z(j)) cycle
1192 0 : cid2 = chem_id(i)
1193 : frac_sum = frac_sum + &
1194 0 : 1d-2*lodders03_element_atom_percent(chem_isos% name(cid2))*chem_isos% W(cid2)
1195 : end do
1196 0 : element_id = Z(j) ! element index equals number of protons
1197 0 : iso_frac = frac/frac_sum ! this iso is iso_frac of the element abundance by mass
1198 0 : new_xa(j) = ztotal*zfrac(element_id)*iso_frac
1199 : end do
1200 :
1201 0 : if (jskip > 0) then
1202 : ! put the missing metals into jskip
1203 0 : new_xa(jskip) = max(0d0, ztotal-sum(new_xa(:)))
1204 0 : do j=1,species
1205 0 : if (j == jskip) cycle
1206 0 : if (Z(j) <= 2) cycle
1207 0 : xa(j) = new_xa(j)
1208 : end do
1209 0 : xa(jskip) = new_xa(jskip)
1210 : else ! renormalize all of the metals
1211 0 : zsum = sum(new_xa(:)) ! only have metals in new_xa
1212 0 : if (zsum > 0d0) then
1213 0 : do j=1,species
1214 0 : if (Z(j) <= 2) cycle
1215 0 : xa(j) = new_xa(j)*ztotal/zsum
1216 : end do
1217 : end if
1218 : end if
1219 :
1220 0 : end subroutine adjust_z_fractions
1221 :
1222 :
1223 0 : subroutine set_z(s, new_z, nzlo, nzhi, ierr)
1224 0 : use net_lib, only:clean_up_fractions
1225 : use star_utils, only: eval_current_z
1226 : type (star_info), pointer :: s
1227 : real(dp), intent(in) :: new_z
1228 : integer, intent(in) :: nzlo, nzhi
1229 : integer, intent(out) :: ierr
1230 :
1231 : integer :: nz, h1, he3, he4, species
1232 0 : real(dp) :: frac_z, current_z, initial_z
1233 : real(dp), parameter :: max_sum_abs = 10d0
1234 : real(dp), parameter :: xsum_tol = 1d-2
1235 :
1236 : include 'formats'
1237 :
1238 0 : ierr = 0
1239 :
1240 0 : if (nzlo > nzhi) then
1241 0 : ierr = -1; return
1242 : end if
1243 :
1244 0 : nz = s% nz
1245 0 : species = s% species
1246 0 : h1 = s% net_iso(ih1)
1247 0 : he3 = s% net_iso(ihe3)
1248 0 : he4 = s% net_iso(ihe4)
1249 :
1250 0 : call clean_up_fractions(1, nz, species, nz, s% xa, max_sum_abs, xsum_tol, ierr)
1251 0 : if (ierr /= 0) return
1252 :
1253 0 : current_z = eval_current_z(s, nzlo, nzhi, ierr)
1254 0 : initial_z = current_z
1255 :
1256 0 : if (current_z <= 0 .and. new_z > 0) then
1257 0 : ierr = -1
1258 0 : return
1259 : end if
1260 0 : frac_z = new_z/current_z
1261 :
1262 0 : if (abs(1-frac_z) < 1d-20) return
1263 :
1264 0 : call convert
1265 0 : if (ierr /= 0) return
1266 :
1267 : ! check
1268 0 : current_z = eval_current_z(s, nzlo, nzhi, ierr)
1269 :
1270 0 : if (abs(current_z-new_z) > 1d-8 .or. is_bad(current_z)) then
1271 0 : ierr = -1
1272 0 : s% retry_message = 'set_z failed'
1273 0 : if (s% report_ierr) then
1274 0 : write(*, *) 'set_z failed'
1275 0 : write(*, 1) 'initial_z', initial_z
1276 0 : write(*, 1) 'requested new_z', new_z
1277 0 : write(*, 1) 'actual current_z', current_z
1278 0 : write(*, *)
1279 0 : write(*, 1) 'requested/initial', new_z/initial_z
1280 0 : write(*, 1) 'requested/final', new_z/current_z
1281 0 : write(*, *)
1282 : end if
1283 0 : if (s% stop_for_bad_nums) then
1284 0 : write(*, 1) 'initial_z', initial_z
1285 0 : write(*, 1) 'requested new_z', new_z
1286 0 : write(*, 1) 'actual current_z', current_z
1287 0 : call mesa_error(__FILE__,__LINE__,'set_z')
1288 : end if
1289 : end if
1290 :
1291 0 : if (s% doing_first_model_of_run) then
1292 0 : s% initial_z = new_z
1293 0 : write(*,1) 's% initial_z =', new_z
1294 : end if
1295 :
1296 0 : s% need_to_setvars = .true.
1297 :
1298 :
1299 : contains
1300 :
1301 0 : subroutine convert
1302 : integer :: i, k
1303 0 : real(dp) :: old_z, old_xy, sumfracs, sum_z, frac_xy
1304 : include 'formats'
1305 0 : do k=nzlo, nzhi
1306 0 : old_z = 0; old_xy = 0
1307 0 : do i=species, 1, -1
1308 0 : if (i==h1 .or. i==he3 .or. i==he4) then
1309 0 : old_xy = old_xy + s% xa(i, k)
1310 : end if
1311 : end do
1312 0 : old_z = 1 - old_xy
1313 0 : sum_z = 0
1314 0 : do i=species, 1, -1
1315 0 : if (i==h1 .or. i==he3 .or. i==he4) cycle
1316 0 : s% xa(i, k) = s% xa(i, k)*frac_z
1317 0 : sum_z = sum_z + s% xa(i, k)
1318 : end do
1319 : ! sum_z is the new z for this cell
1320 : ! modify x and y to make mass fractions sum to 1
1321 0 : if (old_z < 1) then
1322 0 : frac_xy = (1 - old_z*frac_z)/old_xy
1323 0 : s% xa(h1, k) = s% xa(h1, k)*frac_xy
1324 0 : if (he3 /= 0) s% xa(he3, k) = s% xa(he3, k)*frac_xy
1325 0 : s% xa(he4, k) = s% xa(he4, k)*frac_xy
1326 : else
1327 0 : s% xa(h1, k) = (1-sum_z)/2
1328 0 : if (he3 /= 0) s% xa(he3, k) = 0
1329 0 : s% xa(he4, k) = (1-sum_z)/2
1330 : end if
1331 0 : sumfracs = 0
1332 0 : do i=species, 1, -1
1333 0 : sumfracs = sumfracs + s% xa(i, k)
1334 : end do
1335 0 : s% xa(1:species, k) = s% xa(1:species, k)/sumfracs
1336 0 : do i=1, species
1337 0 : if (is_bad(s% xa(i, k))) then
1338 0 : ierr = -1
1339 0 : s% retry_message = 'set_z failed - bad mass fraction'
1340 0 : if (s% report_ierr) then
1341 0 : write(*,2) 'set_z s% xa(i, k)', k, s% xa(i, k)
1342 : end if
1343 0 : if (s% stop_for_bad_nums) then
1344 0 : write(*,2) 'set_z s% xa(i, k)', k, s% xa(i, k)
1345 0 : call mesa_error(__FILE__,__LINE__,'set_z')
1346 : end if
1347 0 : return
1348 : end if
1349 : end do
1350 : end do
1351 0 : end subroutine convert
1352 :
1353 : end subroutine set_z
1354 :
1355 :
1356 0 : subroutine do_replace(s, id1, id2, nzlo, nzhi, ierr)
1357 : ! replaces species1 by species2
1358 : type (star_info), pointer :: s
1359 : integer, intent(in) :: id1, id2 ! values are chem_id's such as ihe4
1360 : integer, intent(in) :: nzlo, nzhi
1361 : integer, intent(out) :: ierr
1362 :
1363 : integer :: k, nz, species, species1, species2
1364 :
1365 0 : ierr = 0
1366 0 : nz = s% nz
1367 0 : species = s% species
1368 0 : species1 = s% net_iso(id1)
1369 0 : species2 = s% net_iso(id2)
1370 :
1371 0 : do k=nzlo, nzhi
1372 0 : s% xa(species2, k) = s% xa(species1, k) + s% xa(species2, k)
1373 0 : s% xa(species1, k) = 0
1374 : end do
1375 :
1376 0 : s% need_to_setvars = .true.
1377 :
1378 0 : end subroutine do_replace
1379 :
1380 :
1381 0 : subroutine do_uniform_mix_section(s, species, nzlo_in, nzhi_in, ierr)
1382 : type (star_info), pointer :: s
1383 : integer, intent(in) :: species, nzlo_in, nzhi_in
1384 : integer, intent(out) :: ierr
1385 : integer :: j, k, nzlo, nzhi
1386 0 : real(dp) :: dqsum, newxa(species), xsum
1387 : include 'formats'
1388 0 : ierr = 0
1389 0 : nzlo = max(1, nzlo_in)
1390 0 : nzhi = min(s% nz, nzhi_in)
1391 0 : dqsum = sum(s% dq(nzlo:nzhi))
1392 0 : do j=1,species
1393 0 : newxa(j) = dot_product(s% dq(nzlo:nzhi), s% xa(j,nzlo:nzhi))/dqsum
1394 : end do
1395 0 : xsum = sum(newxa(:))
1396 0 : do j=1,species
1397 0 : newxa(j) = newxa(j)/xsum
1398 : end do
1399 0 : do k=nzlo,nzhi
1400 0 : do j=1,species
1401 0 : s% xa(j,k) = newxa(j)
1402 : end do
1403 : end do
1404 0 : s% need_to_setvars = .true.
1405 0 : end subroutine do_uniform_mix_section
1406 :
1407 :
1408 0 : subroutine do_uniform_mix_envelope_down_to_T(s, T, ierr)
1409 : type (star_info), pointer :: s
1410 : real(dp), intent(in) :: T
1411 : integer, intent(out) :: ierr
1412 : integer :: k, kmix
1413 0 : ierr = 0
1414 0 : if (s% T(1) > T) return
1415 0 : kmix = s% nz
1416 0 : do k=2,s% nz
1417 0 : if (s% T(k) > T) then
1418 0 : kmix = k-1
1419 0 : exit
1420 : end if
1421 : end do
1422 0 : call do_uniform_mix_section(s, s% species, 1, kmix, ierr)
1423 : end subroutine do_uniform_mix_envelope_down_to_T
1424 :
1425 :
1426 0 : subroutine do_set_abundance(s, chem_id, new_frac, nzlo, nzhi, ierr)
1427 : ! set mass fraction of species to new_frac uniformly in cells nzlo to nzhi
1428 : type (star_info), pointer :: s
1429 : integer, intent(in) :: chem_id
1430 : real(dp), intent(in) :: new_frac
1431 : integer, intent(in) :: nzlo, nzhi
1432 : integer, intent(out) :: ierr
1433 :
1434 : integer :: k, j, i, nz, species
1435 0 : real(dp) :: old_frac, rescale
1436 :
1437 0 : ierr = 0
1438 0 : nz = s% nz
1439 0 : species = s% species
1440 0 : j = s% net_iso(chem_id)
1441 0 : if (j == 0) then
1442 0 : write(*,*) 'do_set_abundance: failed to find requested iso in current net'
1443 0 : ierr = -1
1444 0 : return
1445 : end if
1446 :
1447 0 : do k=nzlo, nzhi
1448 0 : old_frac = s% xa(j, k)
1449 0 : s% xa(j, k) = new_frac
1450 0 : if (1d0-old_frac > 1d-10) then
1451 0 : rescale = (1-new_frac)/(1-old_frac)
1452 0 : do i=1, species
1453 0 : if (i==j) cycle
1454 0 : s% xa(i, k) = rescale*s% xa(i, k)
1455 : end do
1456 : else
1457 0 : rescale = (1-new_frac)/(species-1)
1458 0 : do i=1, species
1459 0 : if (i==j) cycle
1460 0 : s% xa(i, k) = rescale
1461 : end do
1462 : end if
1463 : end do
1464 0 : s% need_to_setvars = .true.
1465 :
1466 : end subroutine do_set_abundance
1467 :
1468 : end module adjust_xyz
|
