Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2017-2018 Rich Townsend & 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 conv_premix
21 :
22 : use const_def, only: dp, i8, msun, convective_mixing
23 : use star_private_def
24 : use chem_def
25 : use chem_lib
26 : use num_lib
27 :
28 : implicit none
29 :
30 : private
31 : public :: do_conv_premix
32 :
33 : integer, parameter :: BURN_TYPE_NONE = 1
34 : integer, parameter :: BURN_TYPE_H = 2
35 : integer, parameter :: BURN_TYPE_HE = 3
36 : integer, parameter :: BURN_TYPE_Z = 4
37 :
38 : integer, parameter :: FIXED_PT_MODE = 5
39 : integer, parameter :: FIXED_DT_MODE = 6
40 :
41 : logical, parameter :: TRACE_ALL = .FALSE.
42 :
43 : ! The zone_info type stores information about the extent &
44 : ! properties of each convective zone
45 :
46 : type zone_info
47 : integer :: kc_t = 0 ! Cell index of top boundary
48 : integer :: kc_b = 0 ! Cell index of bot boundary
49 : real(dp) :: vc_t = 0._dp ! Convective velocity near top boundary
50 : real(dp) :: vc_b = 0._dp ! Convective velocity near bot boundary
51 : real(dp) :: vp_t = 0._dp ! Penetration velocity at top boundary
52 : real(dp) :: vp_b = 0._dp ! Penetration velocity at bot boundary
53 : real(dp) :: dt_t = 0._dp ! Clock for top boundary
54 : real(dp) :: dt_b = 0._dp ! Clock for bottom boundary
55 : integer :: burn_type = BURN_TYPE_NONE ! Type of burning in zone
56 : logical :: sel_t = .FALSE. ! Flag to select top boundary for modification
57 : logical :: sel_b = .FALSE. ! Flag to select bot boundary for modification
58 : logical :: split_retreat = .FALSE. ! Flag to indicate the zone has split, or the trailing bdy retreated
59 : real(dp), allocatable :: avg_xa(:) ! Initial average abundances
60 : real(dp), allocatable :: davg_xa_dt(:) ! Initial rate of change of average abundances (due to burning)
61 : integer :: avoid_inc_iso = 0 ! Isotope id for which increases should be avoided
62 : end type zone_info
63 :
64 : ! The saved_data type saves cell and face data from star_info, to
65 : ! enable us to restore the latter to an earlier state. Only data in
66 : ! cells kc_t:kc_b are used, but the arrays are full-size
67 :
68 : type saved_data
69 : real(dp), allocatable :: xa(:,:)
70 : real(dp), allocatable :: lnd(:)
71 : real(dp), allocatable :: rho(:)
72 : real(dp), allocatable :: lnPgas(:)
73 : real(dp), allocatable :: Pgas(:)
74 : real(dp), allocatable :: gradL_composition_term(:)
75 : integer, allocatable :: update_mode(:)
76 : integer :: kc_t = 0
77 : integer :: kc_b = 0
78 : end type saved_data
79 :
80 : contains
81 :
82 0 : subroutine do_conv_premix (s, ierr)
83 :
84 : use star_utils, only: start_time, update_time
85 :
86 : type(star_info), pointer :: s
87 : integer, intent(out) :: ierr
88 :
89 : logical, parameter :: TRACE_CONV_PREMIX = TRACE_ALL
90 :
91 0 : integer :: update_mode(s%nz)
92 0 : type(saved_data) :: sd
93 0 : type(zone_info), allocatable :: zi(:)
94 : integer :: j_it
95 0 : real(dp), allocatable :: dr(:)
96 0 : real(dp), allocatable :: dt_t(:)
97 0 : real(dp), allocatable :: dt_b(:)
98 : integer :: i_bdy
99 : logical :: t_bdy
100 :
101 : integer(i8) :: time0
102 0 : real(dp) :: total
103 :
104 0 : ierr = 0
105 :
106 0 : if (s% doing_timing) call start_time(s, time0, total)
107 :
108 : ! Initialize the update_mode values. These control how cells will
109 : ! be updated (i.e., how the thermodynamic state will be
110 : ! re-evaluated after abundances have changed). For untouched
111 : ! cells, update_mode is determined by the lnPgas_flag setting
112 :
113 : !if (s%lnPgas_flag) then
114 : ! update_mode = FIXED_PT_MODE
115 : !else
116 0 : update_mode = FIXED_DT_MODE
117 : !end if
118 :
119 : ! Allocate the saved data arrays
120 :
121 0 : call alloc_saved_data_(s, sd)
122 :
123 : ! Initialize the zone info table
124 :
125 0 : call init_zone_info_(s, zi)
126 :
127 : ! Loop until there are no boundaries left to advance
128 :
129 0 : j_it = 0
130 :
131 0 : iter_loop : do
132 :
133 0 : j_it = j_it + 1
134 :
135 : ! Validate the current zone info table
136 :
137 0 : call validate_zone_info_(s, zi)
138 :
139 : ! Check for completion
140 :
141 0 : if (.NOT. (ANY(zi%sel_t .OR. zi%sel_b))) exit iter_loop
142 :
143 : ! Calculate mixing timescales for each boundary (the explicit
144 : ! allocations of dt_t and dt_b are required to workaround a
145 : ! gfortran bug, where allocate-on-assign will raise bogus array
146 : ! bounds violations)
147 :
148 0 : dr = s%r(zi%kc_t+1) - s%r(zi%kc_b)
149 :
150 0 : allocate(dt_t(SIZE(zi)))
151 0 : allocate(dt_b(SIZE(zi)))
152 :
153 0 : dt_t = MERGE(dr/zi%vc_t, HUGE(0._dp), MASK=zi%sel_t)
154 0 : dt_b = MERGE(dr/zi%vc_b, HUGE(0._dp), MASK=zi%sel_b)
155 :
156 : ! Select the boundary with the shortest timescale
157 :
158 0 : i_bdy = MINLOC(MIN(dt_t, dt_b), DIM=1)
159 0 : t_bdy = dt_t(i_bdy) < dt_b(i_bdy)
160 :
161 : if (TRACE_CONV_PREMIX) then
162 : write(*,*) 'Selected i_bdy, t_bdy:', i_bdy, t_bdy
163 : end if
164 :
165 : ! Advance the selected boundary, modifying the model and the
166 : ! zone info table
167 :
168 0 : call advance_bdy_(s, update_mode, zi, i_bdy, t_bdy, sd, j_it)
169 :
170 : ! Loop around
171 :
172 0 : deallocate(dt_t)
173 0 : deallocate(dt_b)
174 :
175 : end do iter_loop
176 0 : s% need_to_setvars = .true.
177 :
178 0 : if (s% doing_timing) call update_time(s, time0, total, s% time_conv_premix)
179 :
180 0 : return
181 :
182 0 : end subroutine do_conv_premix
183 :
184 :
185 0 : subroutine advance_bdy_ (s, update_mode, zi, i_bdy, t_bdy, sd, j_it)
186 :
187 : type(star_info), pointer :: s
188 : integer, intent(inout) :: update_mode(:)
189 : type(zone_info), allocatable, intent(inout) :: zi(:)
190 : integer, intent(inout) :: i_bdy
191 : logical, intent(in) :: t_bdy
192 : type(saved_data), intent(inout) :: sd
193 : integer, intent(in) :: j_it
194 :
195 : logical, parameter :: TRACE_ADVANCE_BDY = TRACE_ALL
196 :
197 : integer :: j_ad
198 : character(256) :: filename
199 :
200 : ! Advance the top (bottom) boundary of the zone zi(i_bdy). t_bdy
201 : ! is .true. for the top boundary, .false. for the bottom. Because
202 : ! the zone info table can potentially change during the process
203 : ! (as zones split/merge), so can i_bdy
204 :
205 : if (TRACE_ADVANCE_BDY) then
206 : write(*,*) 'Advancing zone boundary: i_bdy, t_bdy, j_it=', i_bdy, t_bdy, j_it
207 : end if
208 :
209 0 : if (s% conv_premix_dump_snapshots) then
210 :
211 0 : write(filename, 100) 'cpm-snap.', s%model_number, '.it_', j_it, '.ad_', 0, '.dat'
212 : 100 format(A,I5.5,A,I5.5,A,I5.5,A)
213 :
214 0 : call dump_snapshot_(s, filename)
215 :
216 : end if
217 :
218 0 : j_ad = 0
219 :
220 : advance_loop : do
221 :
222 0 : j_ad = j_ad + 1
223 :
224 : ! Check whether we've reached the edge of the grid
225 :
226 0 : if (t_bdy .AND. zi(i_bdy)%kc_t == 1) then
227 0 : zi(i_bdy)%sel_t = .FALSE.
228 0 : elseif (.NOT. t_bdy .AND. zi(i_bdy)%kc_b == s%nz) then
229 0 : zi(i_bdy)%sel_b = .FALSE.
230 : end if
231 :
232 : ! Check whether the advancing face is still selected
233 :
234 0 : if ((t_bdy .AND. .NOT. zi(i_bdy)%sel_t) .OR. &
235 : (.NOT. t_bdy .AND. .NOT. zi(i_bdy)%sel_b)) exit advance_loop
236 :
237 : ! Advance the boundary by one cell
238 :
239 0 : call advance_bdy_by_one_cell_(s, update_mode, zi, i_bdy, t_bdy, sd)
240 :
241 : ! If necessary, dump a snapshot
242 :
243 0 : if (s% conv_premix_dump_snapshots .AND. MOD(j_ad, 50) == 0) then
244 :
245 0 : write(filename, 100) 'cpm-snap.', s%model_number, '.it_', j_it, '.ad_', j_ad, '.dat'
246 :
247 0 : call dump_snapshot_(s, filename)
248 :
249 : end if
250 :
251 : ! Loop around
252 :
253 : end do advance_loop
254 :
255 0 : if (s% conv_premix_dump_snapshots) then
256 :
257 0 : write(filename, 100) 'cpm-snap.', s%model_number, '.it_', j_it, '.ad_', j_ad-1, '.dat'
258 :
259 0 : call dump_snapshot_(s, filename)
260 :
261 : end if
262 :
263 : if (TRACE_ADVANCE_BDY) then
264 : write(*,*) 'Done advancing zone boundary'
265 : end if
266 :
267 0 : return
268 :
269 0 : end subroutine advance_bdy_
270 :
271 :
272 0 : subroutine advance_bdy_by_one_cell_ (s, update_mode, zi, i_bdy, t_bdy, sd)
273 :
274 : type(star_info), pointer :: s
275 : integer, intent(inout) :: update_mode(:)
276 : type(zone_info), allocatable, intent(inout) :: zi(:)
277 : integer, intent(inout) :: i_bdy
278 : logical, intent(in) :: t_bdy
279 : type(saved_data), intent(inout) :: sd
280 :
281 : integer, parameter :: N_SUB_INI = 1
282 : integer, parameter :: N_SUB_MAX = 1024
283 : logical, parameter :: TRACE_MIX_CELL = TRACE_ALL
284 : logical, parameter :: TRACE_MIX_SUBCELL = TRACE_ALL
285 :
286 0 : real(dp) :: dr
287 0 : real(dp) :: vp
288 0 : real(dp) :: delta_dt
289 0 : real(dp) :: dt_limit
290 : integer :: kc_new
291 : integer :: n_sub
292 : integer :: kc_t
293 : integer :: kc_b
294 : integer :: kf_t
295 : integer :: kf_b
296 0 : real(dp), allocatable :: xa_sub(:,:)
297 0 : real(dp) :: m_sub
298 : logical :: has_split
299 : integer :: kc_sub
300 0 : real(dp) :: m_zone
301 : integer :: l
302 0 : real(dp) :: avg_xa(s%species)
303 : integer :: n_nc
304 : integer :: kf
305 0 : real(dp) :: davg_xa_dt(s%species)
306 :
307 : ! Advance the top (bottom) boundary of the zone zi(i_bdy) by one
308 : ! cell
309 :
310 : ! Calculate the clock increment for adding the new cell
311 :
312 0 : if (t_bdy) then
313 0 : kc_t = zi(i_bdy)%kc_t
314 0 : dr = s%r(kc_t) - s%r(kc_t+1)
315 0 : vp = zi(i_bdy)%vp_t
316 : else
317 0 : kc_b = zi(i_bdy)%kc_b
318 0 : if (kc_b < s%nz) then
319 0 : dr = s%r(kc_b) - s%r(kc_b+1)
320 : else
321 0 : dr = s%r(kc_b)
322 : end if
323 0 : vp = zi(i_bdy)%vp_b
324 : end if
325 :
326 0 : delta_dt = dr/vp
327 :
328 : ! Check whether we have enough time; if not, return
329 :
330 0 : if (s%conv_premix_time_factor > 0._dp) then
331 0 : dt_limit = s%conv_premix_time_factor*s%dt
332 : else
333 : dt_limit = HUGE(0._dp)
334 : end if
335 :
336 0 : if (t_bdy) then
337 :
338 0 : if (zi(i_bdy)%dt_t + delta_dt > dt_limit) then
339 :
340 0 : zi(i_bdy)%sel_t = .FALSE.
341 :
342 : if (TRACE_MIX_CELL) then
343 : write(*,*) ' Outcome: out of time', zi(i_bdy)%dt_t, delta_dt, dt_limit
344 : end if
345 :
346 0 : return
347 :
348 : else
349 :
350 0 : zi(i_bdy)%dt_t = zi(i_bdy)%dt_t + delta_dt
351 :
352 : end if
353 :
354 : else
355 :
356 0 : if (zi(i_bdy)%dt_b + delta_dt > dt_limit) then
357 :
358 0 : zi(i_bdy)%sel_b = .FALSE.
359 :
360 : if (TRACE_MIX_CELL) then
361 : write(*,*) ' Outcome: out of time', zi(i_bdy)%dt_b, delta_dt, dt_limit
362 : end if
363 :
364 0 : return
365 :
366 : else
367 :
368 0 : zi(i_bdy)%dt_b = zi(i_bdy)%dt_b + delta_dt
369 :
370 : end if
371 :
372 : end if
373 :
374 : ! Determine the index of the new cell we're going to add
375 :
376 0 : if (t_bdy) then
377 0 : kc_new = zi(i_bdy)%kc_t - 1
378 : else
379 0 : kc_new = zi(i_bdy)%kc_b + 1
380 : end if
381 :
382 : ! Save the current model state over all cells we *might* modify
383 :
384 0 : if (t_bdy) then
385 0 : call save_model_(s, update_mode, kc_new, zi(i_bdy)%kc_b, sd)
386 : else
387 0 : call save_model_(s, update_mode, zi(i_bdy)%kc_t, kc_new, sd)
388 : end if
389 :
390 : ! Loop over subcell refinement levels
391 :
392 : n_sub = N_SUB_INI
393 :
394 0 : refine_loop : do
395 :
396 : ! Set initial cell and face indices
397 :
398 0 : kc_t = zi(i_bdy)%kc_t
399 0 : kc_b = zi(i_bdy)%kc_b
400 :
401 0 : kf_t = kc_t
402 0 : kf_b = kc_b + 1
403 :
404 : if (TRACE_MIX_CELL) then
405 : write(*,*) 'Attempting to add cell via n_sub subcells:', n_sub
406 : write(*,*) ' kc_t :', kc_t
407 : write(*,*) ' kc_b :', kc_b
408 : write(*,*) ' kc_new :', kc_new
409 : end if
410 :
411 : ! Set update_mode over the cells being mixed
412 :
413 0 : if (s%conv_premix_fix_pgas) then
414 0 : update_mode(kc_t:kc_b) = FIXED_PT_MODE
415 : else
416 0 : update_mode(kc_t:kc_b) = FIXED_DT_MODE
417 : end if
418 :
419 : ! Set gradL_composition_term to zero on interior faces plus the
420 : ! advancing face
421 :
422 0 : if (t_bdy) then
423 0 : s%gradL_composition_term(kf_t:kf_b-1) = 0._dp
424 : else
425 0 : s%gradL_composition_term(kf_t+1:kf_b) = 0._dp
426 : end if
427 :
428 : ! Divide the new cell into n_sub virtual subcells, and store the
429 : ! (initially uniform) abundances of these subcells
430 :
431 0 : if (ALLOCATED(xa_sub)) deallocate(xa_sub)
432 0 : allocate(xa_sub(s%species,n_sub))
433 :
434 0 : do l = 1, s%species
435 0 : xa_sub(l,:) = s%xa(l,kc_new)
436 : end do
437 :
438 0 : m_sub = s%dm(kc_new)/n_sub
439 :
440 : ! Mix subcells into the zone, one by one
441 :
442 0 : has_split = .FALSE.
443 :
444 0 : subcell_loop : do kc_sub = 1, n_sub
445 :
446 : if (TRACE_MIX_SUBCELL) then
447 : write(*,*) 'In subcell loop:', kc_sub, n_sub
448 : end if
449 :
450 : ! Determine average abundances across the zone and the 1:k_s
451 : ! subcells
452 :
453 0 : m_zone = SUM(s%dm(kc_t:kc_b))
454 :
455 0 : do l = 1, s%species
456 : avg_xa(l) = (SUM(s%dm(kc_t:kc_b)*s%xa(l,kc_t:kc_b)) + &
457 0 : SUM(m_sub*xa_sub(l,1:kc_sub))) / (m_zone + m_sub*kc_sub)
458 : end do
459 :
460 0 : avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp)
461 0 : avg_xa = avg_xa/SUM(avg_xa)
462 :
463 : ! Update abundances using the average
464 :
465 0 : do l = 1, s%species
466 0 : s%xa(l,kc_t:kc_b) = avg_xa(l)
467 0 : xa_sub(l,1:kc_sub) = avg_xa(l)
468 : end do
469 :
470 0 : call update_model_(s, update_mode, kc_t, kc_b)
471 :
472 : ! Determine how many interior faces are now non-convective
473 :
474 0 : n_nc = COUNT(s%mlt_mixing_type(kf_t+1:kf_b-1) /= convective_mixing)
475 :
476 0 : if (n_nc > 1 .AND. n_sub < N_SUB_MAX) then
477 :
478 : ! More than a single one: double the number of subcells
479 : ! and restart (as long as we're below the subcell limit)
480 :
481 0 : n_sub = 2*n_sub
482 :
483 0 : call restore_model_(s, update_mode, sd)
484 :
485 : cycle refine_loop
486 :
487 0 : elseif (n_nc > 0) then
488 :
489 : ! A single one (or more, if we're over the subcell
490 : ! limit): move the trailing boundary to its new position
491 : ! (defined as the closest non-convective face to the
492 : ! advancing boundary)
493 :
494 0 : has_split = .TRUE.
495 :
496 0 : if (t_bdy) then
497 :
498 0 : search_down_loop : do kf = kf_t+1, kf_b-1
499 0 : if (s%mlt_mixing_type(kf) /= convective_mixing) exit search_down_loop
500 : end do search_down_loop
501 :
502 0 : kc_b = kf - 1
503 0 : kf_b = kc_b + 1
504 :
505 : if (TRACE_MIX_SUBCELL) then
506 : write(*,*) ' Moved lower boundary to', kc_b
507 : end if
508 :
509 : else
510 :
511 0 : search_up_loop : do kf = kf_b-1, kf_t+1, -1
512 0 : if (s%mlt_mixing_type(kf) /= convective_mixing) exit search_up_loop
513 : end do search_up_loop
514 :
515 0 : kc_t = kf
516 0 : kf_t = kc_t
517 :
518 : if (TRACE_MIX_SUBCELL) then
519 : write(*,*) ' Moved upper boundary to', kc_t
520 : end if
521 :
522 : end if
523 :
524 : end if
525 :
526 : end do subcell_loop
527 :
528 : ! If we reach this point, all went well; exit
529 :
530 : if (TRACE_MIX_SUBCELL) then
531 : write(*,*) ' Completed mixing subcell'
532 : end if
533 :
534 : exit refine_loop
535 :
536 : end do refine_loop
537 :
538 : ! Update the abundances in the new cell, and adjust the cell/face
539 : ! indices
540 :
541 0 : s%xa(:,kc_new) = avg_xa
542 :
543 0 : if (s%conv_premix_fix_pgas) then
544 0 : update_mode(kc_new) = FIXED_PT_MODE
545 : else
546 0 : update_mode(kc_new) = FIXED_DT_MODE
547 : end if
548 :
549 0 : if (t_bdy) then
550 0 : kc_t = kc_new
551 0 : kf_t = kc_t
552 : else
553 0 : kc_b = kc_new
554 0 : kf_b = kc_b + 1
555 : end if
556 :
557 0 : call update_model_(s, update_mode, kc_t, kc_b)
558 :
559 : ! Check whether an abundance increase has occurred; if so, revert
560 : ! back to the starting model and return
561 :
562 0 : if (zi(i_bdy)%avoid_inc_iso /= 0) then
563 :
564 : ! Evaluate the rate of change of the zone-average abundances
565 :
566 0 : davg_xa_dt = (avg_xa - zi(i_bdy)%avg_xa)/s%dt
567 :
568 : ! Check whether an increase will occur (davg_xa_dt is positive,
569 : ! and greater in magnitude than the rate of decrease due to
570 : ! burning)
571 :
572 : associate (iso => zi(i_bdy)%avoid_inc_iso)
573 :
574 0 : if (davg_xa_dt(iso) > MAX(-zi(i_bdy)%davg_xa_dt(iso), 0._dp)) then
575 :
576 0 : call restore_model_(s, update_mode, sd)
577 :
578 0 : if (t_bdy) then
579 0 : zi(i_bdy)%sel_t = .FALSE.
580 : else
581 0 : zi(i_bdy)%sel_b = .FALSE.
582 : end if
583 :
584 : if (TRACE_MIX_CELL) then
585 : write(*,*) ' Outcome: abundance reversal for isotope', iso
586 : end if
587 :
588 0 : return
589 :
590 : end if
591 :
592 : end associate
593 :
594 : end if
595 :
596 : ! Determine whether the face just inside the advancing boundary
597 : ! has become/remained convective; if not, revert back to the
598 : ! starting model and return
599 :
600 0 : if (t_bdy) then
601 :
602 0 : if (s%mlt_mixing_type(kf_t+1) /= convective_mixing) then
603 :
604 0 : call restore_model_(s, update_mode, sd)
605 :
606 0 : zi(i_bdy)%sel_t = .FALSE.
607 :
608 : if (TRACE_MIX_CELL) then
609 : write(*,*) ' Outcome: non-convective at top'
610 : end if
611 :
612 0 : return
613 :
614 : end if
615 :
616 : else
617 :
618 0 : if (s%mlt_mixing_type(kf_b-1) /= convective_mixing) then
619 :
620 0 : call restore_model_(s, update_mode, sd)
621 :
622 0 : zi(i_bdy)%sel_b = .FALSE.
623 :
624 : if (TRACE_MIX_CELL) then
625 : write(*,*) ' Outcome: non-convective at bottom'
626 : end if
627 :
628 0 : return
629 :
630 : end if
631 :
632 : end if
633 :
634 : ! Update the zone info table to reflect all the changes we made
635 :
636 0 : call update_zone_info_(s, i_bdy, t_bdy, has_split, zi)
637 :
638 0 : return
639 :
640 0 : end subroutine advance_bdy_by_one_cell_
641 :
642 :
643 0 : subroutine init_zone_info_ (s, zi)
644 :
645 : type(star_info), pointer :: s
646 : type(zone_info), allocatable, intent(out) :: zi(:)
647 :
648 : integer :: i
649 :
650 : ! Initialize the zone info table for the whole star
651 :
652 0 : call create_zone_info_(s, 1, s%nz, zi)
653 :
654 : ! Perform additional initializations
655 :
656 0 : zone_loop : do i = 1, SIZE(zi)
657 :
658 : ! Set penetration velocities
659 :
660 0 : call set_penetration_velocities_(s, zi(i))
661 :
662 : ! Set burn types and initial average abundances
663 :
664 0 : call set_burn_data_(s, zi(i))
665 0 : call set_abund_data_(s, zi(i))
666 :
667 : ! Initialize clocks
668 :
669 0 : zi(i)%dt_t = 0._dp
670 0 : zi(i)%dt_b = 0._dp
671 :
672 : ! Set flags
673 :
674 0 : zi(i)%sel_t = zi(i)%kc_t /= 1
675 0 : zi(i)%sel_b = zi(i)%kc_b /= s%nz
676 :
677 0 : zi(i)%split_retreat = .FALSE.
678 :
679 : end do zone_loop
680 :
681 0 : return
682 :
683 : end subroutine init_zone_info_
684 :
685 :
686 0 : subroutine create_zone_info_ (s, kc_t, kc_b, zi)
687 :
688 : type(star_info), pointer :: s
689 : integer, intent(in) :: kc_t
690 : integer, intent(in) :: kc_b
691 : type(zone_info), allocatable, intent(out) :: zi(:)
692 :
693 : logical :: in_conv
694 0 : type(zone_info) :: zi_new
695 : integer :: kf
696 :
697 : ! Create a zone info table spanning the cells kc_t:kc_b. For the
698 : ! purposes of convective premixing, a zone is defined as a regions
699 : ! of two or more cells with (i) convective faces inside the
700 : ! region, (ii) non-convective faces at the top and bottom
701 : ! boundaries of the region. The possible exceptions to these
702 : ! rules apply to zones which begin (end) at kc_t (kc_b); their
703 : ! upper (lower) boundary faces are not checked for being
704 : ! non-convective.
705 :
706 0 : if (kc_t >= kc_b) then
707 0 : write(*,*) 'conv_premix: Invalid cell range in create_zone_info_'
708 0 : stop
709 : end if
710 :
711 : ! Create the empty zone info table
712 :
713 0 : allocate(zi(0))
714 :
715 : ! Loop through faces, looking for zone boundaries and adding zones
716 : ! (cf. locate_convection_boundaries in mix_info.f90)
717 :
718 0 : kf = kc_b
719 :
720 0 : in_conv = s%mlt_mixing_type(kf) == convective_mixing
721 :
722 0 : if (in_conv) then
723 :
724 0 : zi_new%kc_b = kf
725 0 : zi_new%vc_b = s%mlt_vc(kf)
726 :
727 : end if
728 :
729 0 : face_loop : do kf = kc_b-1, kc_t+1, -1
730 :
731 0 : if (in_conv .AND. s%mlt_mixing_type(kf) /= convective_mixing) then
732 :
733 : ! Transitioning out of a convective zone; complete
734 : ! definition of new zone info and add to to the table
735 :
736 0 : zi_new%kc_t = kf
737 0 : zi_new%vc_t = s%mlt_vc(kf+1)
738 :
739 0 : zi = [zi,zi_new]
740 :
741 0 : in_conv = .FALSE.
742 :
743 0 : elseif (.NOT. in_conv .AND. s%mlt_mixing_type(kf) == convective_mixing) then
744 :
745 : ! Transitioning into a convective zone; begin definition of
746 : ! new zone info
747 :
748 0 : zi_new%kc_b = kf
749 0 : zi_new%vc_b = s%mlt_vc(kf)
750 :
751 0 : in_conv = .TRUE.
752 :
753 : end if
754 :
755 : end do face_loop
756 :
757 0 : kf = kc_t
758 :
759 0 : if (in_conv) then
760 :
761 0 : zi_new%kc_t = kf
762 0 : zi_new%vc_t = s%mlt_vc(kf+1)
763 :
764 0 : zi = [zi,zi_new]
765 :
766 : end if
767 :
768 0 : return
769 :
770 : end subroutine create_zone_info_
771 :
772 :
773 0 : subroutine update_zone_info_ (s, i_bdy, t_bdy, has_split, zi)
774 :
775 : type(star_info), pointer :: s
776 : integer, intent(inout) :: i_bdy
777 : logical, intent(in) :: t_bdy
778 : logical, intent(in) :: has_split
779 : type(zone_info), allocatable, intent(inout) :: zi(:)
780 :
781 : logical, parameter :: TRACE_UPDATE_ZONE = TRACE_ALL
782 :
783 0 : type(zone_info), allocatable :: zi_new(:)
784 : integer :: i
785 :
786 : ! Update the zone info table to account for a single-cell
787 : ! advancement in zone i_bdy, direction t_bdy
788 :
789 : ! Advance the zone by one cell
790 :
791 0 : if (t_bdy) then
792 0 : zi(i_bdy)%kc_t = zi(i_bdy)%kc_t - 1
793 0 : zi(i_bdy)%vc_t = s%mlt_vc(zi(i_bdy)%kc_t+1)
794 : else
795 0 : zi(i_bdy)%kc_b = zi(i_bdy)%kc_b + 1
796 0 : zi(i_bdy)%vc_b = s%mlt_vc(zi(i_bdy)%kc_b)
797 : end if
798 :
799 0 : call set_penetration_velocities_(s, zi(i_bdy))
800 :
801 : ! If necessary, truncate (or even delete) the adjacent zone we've
802 : ! encroached into
803 :
804 0 : if (t_bdy .AND. i_bdy < SIZE(zi)) then
805 :
806 0 : if (zi(i_bdy)%kc_t == zi(i_bdy+1)%kc_b) then
807 :
808 : ! Encroached into zone above
809 :
810 0 : if (zi(i_bdy+1)%kc_b-zi(i_bdy+1)%kc_t > 1) then
811 :
812 : ! Truncate the zone
813 :
814 0 : zi(i_bdy+1)%kc_b = zi(i_bdy+1)%kc_b - 1
815 0 : zi(i_bdy+1)%vc_b = s%mlt_vc(zi(i_bdy+1)%kc_b)
816 :
817 : if (TRACE_UPDATE_ZONE) then
818 : write(*,*) 'Truncated zone above to', zi(i_bdy+1)%kc_b, zi(i_bdy+1)%vc_b, &
819 : s%mlt_mixing_type(zi(i_bdy+1)%kc_b), zi(i_bdy+1)%kc_b-zi(i_bdy+1)%kc_t+1
820 : end if
821 :
822 : else
823 :
824 : ! Delete the zone
825 :
826 0 : zi = [zi(:i_bdy),zi(i_bdy+2:)]
827 :
828 : if (TRACE_UPDATE_ZONE) then
829 : write(*,*) 'Deleted zone above'
830 : end if
831 :
832 : end if
833 :
834 : end if
835 :
836 0 : elseif (.NOT. t_bdy .AND. i_bdy > 1) then
837 :
838 0 : if (zi(i_bdy)%kc_b == zi(i_bdy-1)%kc_t) then
839 :
840 : ! Encroached into zone below
841 :
842 0 : if (zi(i_bdy-1)%kc_b-zi(i_bdy-1)%kc_t > 1) then
843 :
844 : ! Truncate the zone
845 :
846 0 : zi(i_bdy-1)%kc_t = zi(i_bdy-1)%kc_t + 1
847 0 : zi(i_bdy-1)%vc_t = s%mlt_vc(zi(i_bdy-1)%kc_t+1)
848 :
849 : if (TRACE_UPDATE_ZONE) then
850 : write(*,*) 'Truncated zone below to', zi(i_bdy-1)%kc_t
851 : end if
852 :
853 : else
854 :
855 : ! Delete the zone
856 :
857 0 : zi = [zi(:i_bdy-2),zi(i_bdy:)]
858 0 : i_bdy = i_bdy - 1
859 :
860 : if (TRACE_UPDATE_ZONE) then
861 : write(*,*) 'Deleted zone below'
862 : end if
863 :
864 : end if
865 :
866 : end if
867 :
868 : end if
869 :
870 : ! If the zone has split (or its trailing boundary has retreated),
871 : ! then replace the zone with as many new zones as necessary
872 :
873 0 : if (has_split) then
874 :
875 : ! Create new zones to replace the zone
876 :
877 0 : call create_zone_info_(s, zi(i_bdy)%kc_t, zi(i_bdy)%kc_b, zi_new)
878 :
879 : ! Set up parameters for the new zones
880 :
881 0 : new_zone_loop : do i = 1, SIZE(zi_new)
882 :
883 : ! Clocks & selection flags
884 :
885 0 : if (zi_new(i)%kc_t == zi(i_bdy)%kc_t) then
886 0 : zi_new(i)%dt_t = zi(i_bdy)%dt_t
887 0 : zi_new(i)%sel_t = zi(i_bdy)%sel_t
888 : else
889 0 : if (t_bdy) then
890 0 : zi_new(i)%dt_t = zi(i_bdy)%dt_t
891 0 : zi_new(i)%sel_t = .FALSE.
892 : else
893 0 : zi_new(i)%dt_t = zi(i_bdy)%dt_b
894 0 : zi_new(i)%sel_t = .FALSE.
895 : end if
896 : end if
897 :
898 0 : if (zi_new(i)%kc_b == zi(i_bdy)%kc_b) then
899 0 : zi_new(i)%dt_b = zi(i_bdy)%dt_b
900 0 : zi_new(i)%sel_b = zi(i_bdy)%sel_b
901 : else
902 0 : if (t_bdy) then
903 0 : zi_new(i)%dt_b = zi(i_bdy)%dt_t
904 0 : zi_new(i)%sel_b = .FALSE.
905 : else
906 0 : zi_new(i)%dt_b = zi(i_bdy)%dt_b
907 0 : zi_new(i)%sel_b = .FALSE.
908 : end if
909 : end if
910 :
911 : ! Penetration velocities
912 :
913 0 : call set_penetration_velocities_(s, zi_new(i))
914 :
915 : ! Burn data
916 :
917 0 : call set_burn_data_(s, zi_new(i))
918 :
919 : ! Initial abundances
920 :
921 0 : zi_new(i)%avg_xa = zi(i_bdy)%avg_xa
922 0 : zi_new(i)%davg_xa_dt = zi(i_bdy)%davg_xa_dt
923 :
924 : end do new_zone_loop
925 :
926 : ! Replace the zone with the new zones
927 :
928 0 : zi = [zi(:i_bdy-1),zi_new,zi(i_bdy+1:)]
929 :
930 : if (TRACE_UPDATE_ZONE) then
931 : write(*,*) 'Replaced zone with new zones', SIZE(zi), SIZE(zi_new)
932 : end if
933 :
934 : else
935 :
936 : if (TRACE_UPDATE_ZONE) then
937 : write(*,*) 'Zone did not split, updating indices',&
938 : COUNT(s%mlt_mixing_type(zi(i_bdy)%kc_t+1:zi(i_bdy)%kc_b) /= convective_mixing)
939 : end if
940 :
941 : end if
942 :
943 0 : return
944 :
945 0 : end subroutine update_zone_info_
946 :
947 :
948 0 : subroutine set_penetration_velocities_ (s, zi)
949 :
950 : type(star_info), pointer :: s
951 : type(zone_info), intent(inout) :: zi
952 :
953 0 : real(dp) :: mu_i
954 0 : real(dp) :: mu_e
955 : integer :: kf
956 0 : real(dp) :: z_s
957 :
958 : ! Set the penetration velocities for the zone info. These are
959 : ! calculated by evaluating the molecular weight contrast between
960 : ! the boundary cell and the adjacent (radiative) cell, and
961 : ! calculating a characteristic penetration velocity using the
962 : ! approach by Castellani (1971)
963 :
964 : ! Bottom boundary
965 :
966 0 : if (zi%kc_b < s%nz) then
967 :
968 0 : mu_i = s%mu(zi%kc_b)
969 0 : mu_e = s%mu(zi%kc_b+1)
970 :
971 0 : kf = zi%kc_b
972 :
973 0 : if (mu_e - mu_i > EPSILON(0._dp)) then
974 0 : z_s = MIN(zi%vc_b*zi%vc_b/(2._dp*s%grav(kf)*(1._dp - mu_i/mu_e)), s%mlt_mixing_length(kf))
975 : else
976 0 : z_s = s%mlt_mixing_length(kf)
977 : end if
978 :
979 0 : zi%vp_b = zi%vc_b*z_s/s%mlt_mixing_length(kf)
980 :
981 : else
982 :
983 0 : zi%vp_b = zi%vc_b
984 :
985 : end if
986 :
987 0 : if (zi%vp_b == 0._dp) print *,'Bottom bdy has zero vp', zi%kc_b, zi%vp_b, zi%vc_b, z_s
988 :
989 : ! Top boundary
990 :
991 0 : if (zi%kc_t > 1) then
992 :
993 0 : mu_i = s%mu(zi%kc_t)
994 0 : mu_e = s%mu(zi%kc_t-1)
995 :
996 0 : kf = zi%kc_t + 1
997 :
998 0 : if (mu_i - mu_e > EPSILON(0._dp)) then
999 0 : z_s = MIN(zi%vc_t*zi%vc_t/(2._dp*s%grav(kf)*(1._dp - mu_e/mu_i)), s%mlt_mixing_length(kf))
1000 : else
1001 0 : z_s = s%mlt_mixing_length(kf)
1002 : end if
1003 :
1004 0 : zi%vp_t = zi%vc_t*z_s/s%mlt_mixing_length(kf)
1005 :
1006 : else
1007 :
1008 0 : zi%vp_t = zi%vc_t
1009 :
1010 : end if
1011 :
1012 0 : if (zi%vp_t == 0._dp) print *,'Top bdy has zero vp', zi%kc_t, zi%vp_t, zi%vc_t, z_s
1013 :
1014 0 : return
1015 :
1016 : end subroutine set_penetration_velocities_
1017 :
1018 :
1019 0 : subroutine set_burn_data_ (s, zi)
1020 :
1021 : type(star_info), pointer :: s
1022 : type(zone_info), intent(inout) :: zi
1023 :
1024 : real(dp), parameter :: EPS_THRESHOLD = 1._dp ! Threshold eps for a zone to be classified as burning
1025 :
1026 : integer :: iso_h1
1027 : integer :: iso_he4
1028 : integer :: iso_c12
1029 0 : real(dp) :: eps_max
1030 0 : real(dp) :: eps_h_max
1031 0 : real(dp) :: eps_he_max
1032 : integer :: kc
1033 :
1034 : ! Set burning data for the zone info
1035 :
1036 0 : iso_h1 = s%net_iso(chem_get_iso_id('h1'))
1037 0 : iso_he4 = s%net_iso(chem_get_iso_id('he4'))
1038 0 : iso_c12 = s%net_iso(chem_get_iso_id('c12'))
1039 :
1040 : associate (kc_t => zi%kc_t, &
1041 : kc_b => zi%kc_b)
1042 :
1043 : ! Find the maximum burning rate, and the H/He burning rate at
1044 : ! the same location
1045 :
1046 0 : eps_max = -HUGE(0._dp)
1047 0 : eps_h_max = -HUGE(0._dp)
1048 0 : eps_he_max = -HUGE(0._dp)
1049 :
1050 0 : cell_loop : do kc = kc_t, kc_b
1051 :
1052 0 : if (s%eps_nuc(kc) > eps_max) then
1053 :
1054 0 : eps_max = s%eps_nuc(kc)
1055 :
1056 : eps_h_max = s%eps_nuc_categories(ipp, kc) + &
1057 0 : s%eps_nuc_categories(icno, kc)
1058 : eps_he_max = s%eps_nuc_categories(i3alf, kc) + &
1059 0 : s%eps_nuc_categories(i_burn_c, kc)
1060 :
1061 : end if
1062 :
1063 : end do cell_loop
1064 :
1065 : ! Classify the principal burning type
1066 :
1067 0 : if (eps_max > EPS_THRESHOLD) then
1068 :
1069 0 : if (eps_h_max > 0.5_dp*eps_max) then
1070 :
1071 : ! Hydrogen burning zone
1072 :
1073 0 : zi%burn_type = BURN_TYPE_H
1074 :
1075 0 : elseif (eps_he_max > 0.5_dp*eps_max) then
1076 :
1077 : ! Helium burning zone
1078 :
1079 0 : zi%burn_type = BURN_TYPE_HE
1080 :
1081 : else
1082 :
1083 : ! Metal burning zone
1084 :
1085 0 : zi%burn_type = BURN_TYPE_Z
1086 :
1087 : end if
1088 :
1089 : else
1090 :
1091 : ! Non-burning zone
1092 :
1093 0 : zi%burn_type = BURN_TYPE_NONE
1094 :
1095 : end if
1096 :
1097 : end associate
1098 :
1099 : ! Based on the burning type, decide which isotope (if any)
1100 : ! should be monitored for abundance reversal avoidance
1101 :
1102 0 : if (s%conv_premix_avoid_increase) then
1103 :
1104 0 : select case (zi%burn_type)
1105 : case (BURN_TYPE_H)
1106 0 : zi%avoid_inc_iso = iso_h1
1107 : case (BURN_TYPE_HE)
1108 0 : zi%avoid_inc_iso = iso_he4
1109 : case default
1110 0 : zi%avoid_inc_iso = 0
1111 : end select
1112 :
1113 : end if
1114 :
1115 0 : return
1116 :
1117 : end subroutine set_burn_data_
1118 :
1119 :
1120 0 : subroutine set_abund_data_ (s, zi)
1121 :
1122 : type(star_info), pointer :: s
1123 : type(zone_info), intent(inout) :: zi
1124 :
1125 0 : real(dp) :: m_zone
1126 : integer :: l
1127 :
1128 : ! Set abundance data for the zone info
1129 :
1130 0 : allocate(zi%avg_xa(s%species))
1131 0 : allocate(zi%davg_xa_dt(s%species))
1132 :
1133 : associate (kc_t => zi%kc_t, &
1134 : kc_b => zi%kc_b, &
1135 : avg_xa => zi%avg_xa, &
1136 : davg_xa_dt => zi%davg_xa_dt)
1137 :
1138 : ! Calculate average abundances and their rates of change
1139 :
1140 0 : m_zone = SUM(s%dm(kc_t:kc_b))
1141 :
1142 0 : do l = 1, s%species
1143 0 : avg_xa(l) = SUM(s%dm(kc_t:kc_b)*s%xa(l,kc_t:kc_b))/m_zone
1144 0 : davg_xa_dt(l) = SUM(s%dm(kc_t:kc_b)*s%dxdt_nuc(l,kc_t:kc_b))/m_zone
1145 : end do
1146 :
1147 0 : avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp)
1148 0 : avg_xa = avg_xa/SUM(avg_xa)
1149 :
1150 : end associate
1151 :
1152 0 : return
1153 :
1154 : end subroutine set_abund_data_
1155 :
1156 :
1157 0 : subroutine validate_zone_info_ (s, zi)
1158 :
1159 : type(star_info), pointer :: s
1160 : type(zone_info), intent(in) :: zi(:)
1161 :
1162 : logical, parameter :: TRACE_VALIDATE_INFO = TRACE_ALL
1163 :
1164 : logical :: valid
1165 : integer :: n_zi
1166 : integer :: i
1167 :
1168 : integer :: kf
1169 :
1170 : ! Validate zone info data
1171 :
1172 0 : n_zi = SIZE(zi)
1173 :
1174 0 : valid = .TRUE.
1175 :
1176 0 : zone_info_loop : do i = 1, n_zi
1177 :
1178 : if (TRACE_VALIDATE_INFO) then
1179 : write(*,*) 'Validating info for zone', i
1180 : call print_zone_info_(s, zi(i))
1181 : end if
1182 :
1183 : ! (i) Degenerate zone (just one cell)
1184 :
1185 0 : if (zi(i)%kc_t == zi(i)%kc_b) then
1186 0 : write(*,*) 'conv_premix: Degenerate zone'
1187 0 : valid = .FALSE.
1188 : end if
1189 :
1190 : ! (iii) Out of bounds zone (cells outside grid)
1191 :
1192 0 : if (zi(i)%kc_t < 1) then
1193 0 : write(*,*) 'conv_premix: Zone top outside grid'
1194 0 : valid = .FALSE.
1195 : end if
1196 :
1197 0 : if (zi(i)%kc_b > s%nz) then
1198 0 : write(*,*) 'conv_premix: Zone bottom outside grid'
1199 0 : valid = .FALSE.
1200 : end if
1201 :
1202 : ! (iv) Overlapping zone (cells inside previous/next zone)
1203 :
1204 0 : if (i > 1) then
1205 0 : if (zi(i)%kc_b >= zi(max(1,i-1))%kc_t) then
1206 : ! bp: max(1,i-1) to prevent bogus warning from gfortran
1207 0 : write(*,*) 'conv_premix: Zone bottom inside previous zone'
1208 0 : valid = .FALSE.
1209 : end if
1210 : end if
1211 :
1212 0 : if (i < n_zi) then
1213 0 : if (zi(i)%kc_t <= zi(i+1)%kc_b) then
1214 0 : write(*,*) 'conv_premix: Zone top inside next zone'
1215 0 : valid = .FALSE.
1216 : end if
1217 : end if
1218 :
1219 : ! (vi) Convective velocities
1220 :
1221 0 : if (zi(i)%vc_t <= 0._dp) then
1222 0 : write(*,*) 'conv_premix: Convective velocity is not greater than zero at zone top'
1223 0 : valid = .FALSE.
1224 : end if
1225 :
1226 0 : if (zi(i)%vc_b <= 0._dp) then
1227 0 : write(*,*) 'conv_premix: Convective velocity is not greater than zero at zone bottom'
1228 0 : valid = .FALSE.
1229 : end if
1230 :
1231 : ! (vii) Mixing type
1232 :
1233 0 : if (.NOT. ALL(s%mlt_mixing_type(zi(i)%kc_t+1:zi(i)%kc_b) == convective_mixing)) then
1234 0 : do kf = zi(i)%kc_t+1, zi(i)%kc_b
1235 0 : if (s%mlt_mixing_type(kf) /= convective_mixing) then
1236 0 : write(*,*) 'Mix type at kf=',kf, s%mlt_mixing_type(kf)
1237 : end if
1238 : end do
1239 0 : write(*,*) 'conv_premix: Zone contains cells with mixing_type /= convective_mixing'
1240 0 : valid = .FALSE.
1241 : end if
1242 :
1243 : end do zone_info_loop
1244 :
1245 0 : if (.NOT. valid) stop
1246 :
1247 0 : return
1248 :
1249 : end subroutine validate_zone_info_
1250 :
1251 :
1252 : subroutine print_zone_info_ (s, zi)
1253 :
1254 : type(star_info), pointer :: s
1255 : type(zone_info), intent(in) :: zi
1256 :
1257 : ! Print out zone info
1258 :
1259 : write(*,*) ' nz :', s%nz
1260 : write(*,*) ' kc_t :', zi%kc_t
1261 : write(*,*) ' kc_b :', zi%kc_b
1262 : write(*,*) ' mass(kf_t) :', (s%M_center + s%xmstar*s%q(zi%kc_t))/Msun
1263 : if (zi%kc_b < s%nz) then
1264 : write(*,*) ' mass(kf_b) :', (s%M_center + s%xmstar*s%q(zi%kc_b+1))/Msun
1265 : else
1266 : write(*,*) ' mass(kf_b) :', (s%M_center)/Msun
1267 : end if
1268 : write(*,*) ' vc_t :', zi%vc_t
1269 : write(*,*) ' vc_b :', zi%vc_b
1270 : write(*,*) ' dt_t :', zi%dt_t
1271 : write(*,*) ' dt_b :', zi%dt_b
1272 : write(*,*) ' burn_type :', zi%burn_type
1273 : write(*,*) ' sel_t :', zi%sel_t
1274 : write(*,*) ' sel_b :', zi%sel_b
1275 : write(*,*) ' avoid_inc_iso :', zi%avoid_inc_iso
1276 : write(*,'(A)')
1277 :
1278 : return
1279 :
1280 : end subroutine print_zone_info_
1281 :
1282 :
1283 0 : subroutine update_model_ (s, update_mode, kc_t, kc_b)
1284 :
1285 : use turb_info, only: set_mlt_vars
1286 : use micro
1287 :
1288 : type(star_info), pointer :: s
1289 : integer, intent(in) :: update_mode(:)
1290 : integer, intent(in) :: kc_t
1291 : integer, intent(in) :: kc_b
1292 :
1293 : logical, parameter :: TRACE_UPDATE_MODEL = .FALSE.
1294 :
1295 : integer :: ierr
1296 : integer :: kf_t
1297 : integer :: kf_b
1298 :
1299 : ! Update the model to reflect changes in the abundances across
1300 : ! cells kc_t:kc_b
1301 :
1302 : if (TRACE_UPDATE_MODEL) then
1303 : write(*,*) ' Updating model'
1304 : end if
1305 :
1306 : ! Perform the EOS updates, in accordance with the update_mode
1307 : ! values. We make two passes -- the first to handle cells with
1308 : ! update_mode == FIXED_PT_MODE, the second to handle cells with
1309 : ! update_mode == FIXED_DT_MODE.
1310 :
1311 0 : s%fix_Pgas = .TRUE.
1312 :
1313 0 : call set_eos_with_mask(s, kc_t, kc_b, update_mode==FIXED_PT_MODE, ierr)
1314 0 : if (ierr /= 0) then
1315 0 : write(*,*) 'conv_premix: error from call to set_eos_with_mask'
1316 0 : stop
1317 : end if
1318 :
1319 0 : s%fix_Pgas = .FALSE.
1320 :
1321 0 : call set_eos_with_mask(s, kc_t, kc_b, update_mode==FIXED_DT_MODE, ierr)
1322 0 : if (ierr /= 0) then
1323 0 : write(*,*) 'conv_premix: error from call to set_eos_with_mask'
1324 0 : stop
1325 : end if
1326 :
1327 : ! Update opacities across cells kc_t:kc_b (this also sets rho_face
1328 : ! and related quantities on faces kc_t:kc_b)
1329 :
1330 : call set_micro_vars(s, kc_t, kc_b, &
1331 0 : skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.FALSE., ierr=ierr)
1332 0 : if (ierr /= 0) then
1333 0 : write(*,*) 'conv_premix: error from call to set_micro_vars'
1334 0 : stop
1335 : end if
1336 :
1337 : ! Finally update MLT for interior faces
1338 :
1339 0 : kf_t = kc_t
1340 0 : kf_b = kc_b + 1
1341 :
1342 0 : call set_mlt_vars(s, kf_t+1, kf_b-1, ierr)
1343 0 : if (ierr /= 0) then
1344 0 : write(*,*) 'conv_premix: failed in call to set_mlt_vars during update_model_'
1345 0 : stop
1346 : end if
1347 :
1348 0 : return
1349 :
1350 0 : end subroutine update_model_
1351 :
1352 :
1353 0 : subroutine dump_snapshot_ (s, filename)
1354 :
1355 0 : use brunt
1356 : use hydro_vars, only: set_grads
1357 :
1358 : type(star_info), pointer :: s
1359 : character(*), intent(in) :: filename
1360 :
1361 0 : real(dp), allocatable :: gradL_composition_term(:)
1362 : integer :: ierr
1363 : integer :: unit
1364 : integer :: iso_h1
1365 : integer :: iso_he4
1366 : integer :: iso_c12
1367 : integer :: iso_o16
1368 : integer :: k
1369 :
1370 : ! Dump a snapshot of the current model state to file
1371 :
1372 : ! First, re-calculate gradL_composition_term
1373 0 : allocate(gradL_composition_term(1:s%nz))
1374 :
1375 0 : gradL_composition_term = s%gradL_composition_term(1:s%nz)
1376 :
1377 0 : call do_brunt_B(s, 1, s%nz, ierr)
1378 0 : if (ierr /= 0) then
1379 0 : write(*,*) 'conv_premix: error from call to do_brunt_B'
1380 0 : stop
1381 : end if
1382 :
1383 0 : call set_grads(s, ierr)
1384 0 : if (ierr /= 0) then
1385 0 : write(*,*) 'conv_premix: error from call to set_grads'
1386 0 : stop
1387 : end if
1388 :
1389 : ! Dump the snapshot
1390 :
1391 0 : open(NEWUNIT=unit, FILE=filename, STATUS='REPLACE')
1392 :
1393 0 : iso_h1 = s%net_iso(chem_get_iso_id('h1'))
1394 0 : iso_he4 = s%net_iso(chem_get_iso_id('he4'))
1395 0 : iso_c12 = s%net_iso(chem_get_iso_id('c12'))
1396 0 : iso_o16 = s%net_iso(chem_get_iso_id('o16'))
1397 :
1398 0 : write(*,*) ' Dumping to file:', TRIM(filename)
1399 :
1400 0 : write(unit, *) '1 2'
1401 0 : write(unit, *) 'star_age model_number'
1402 0 : write(unit, *) s%star_age, s%model_number
1403 0 : write(unit, *) ''
1404 0 : write(unit, *) '1 2 3 4 5 6 7 8'
1405 0 : write(unit, *) 'k mass gradr grada gradL_composition_term mixing_type x y'
1406 :
1407 0 : do k = 1, s%nz
1408 0 : write(unit, *) k, s%m(k)/msun, s%gradr(k), s%grada_face(k), &
1409 0 : s%gradL_composition_term(k), s%mlt_mixing_type(k), s%xa(iso_h1,k), s%xa(iso_he4,k)
1410 : end do
1411 :
1412 0 : close(unit)
1413 :
1414 : ! Restore gradL_composition_term
1415 :
1416 0 : s%gradL_composition_term(1:s%nz) = gradL_composition_term
1417 :
1418 0 : return
1419 :
1420 0 : end subroutine dump_snapshot_
1421 :
1422 :
1423 0 : subroutine alloc_saved_data_ (s, sd)
1424 :
1425 : type(star_info), pointer :: s
1426 : type(saved_data) :: sd
1427 :
1428 : ! Allocate cell data arrays
1429 :
1430 0 : allocate(sd%xa(s%species,s%nz))
1431 :
1432 0 : allocate(sd%lnd(s%nz))
1433 0 : allocate(sd%rho(s%nz))
1434 :
1435 0 : allocate(sd%lnPgas(s%nz))
1436 0 : allocate(sd%Pgas(s%nz))
1437 :
1438 0 : allocate(sd%update_mode(s%nz))
1439 :
1440 0 : allocate(sd%gradL_composition_term(s%nz))
1441 :
1442 0 : end subroutine alloc_saved_data_
1443 :
1444 :
1445 0 : subroutine save_model_ (s, update_mode, kc_t, kc_b, sd)
1446 :
1447 : type(star_info), pointer :: s
1448 : integer, intent(in) :: update_mode(:)
1449 : integer, intent(in) :: kc_t
1450 : integer, intent(in) :: kc_b
1451 : type(saved_data), intent(inout) :: sd
1452 :
1453 : logical, parameter :: TRACE_SAVE_MODEL = .FALSE.
1454 :
1455 : integer :: kf_t
1456 : integer :: kf_b
1457 :
1458 : ! Save data for cells kc_t:kc_b and interior faces
1459 :
1460 : if (TRACE_SAVE_MODEL) then
1461 : write(*,*) ' Saving model'
1462 : end if
1463 :
1464 : ! Save cell data
1465 :
1466 : if (TRACE_SAVE_MODEL) then
1467 : write(*,*) ' kc_t:', kc_t
1468 : write(*,*) ' kc_b:', kc_b
1469 : end if
1470 :
1471 0 : sd%update_mode(kc_t:kc_b) = update_mode(kc_t:kc_b)
1472 :
1473 0 : sd%xa(:,kc_t:kc_b) = s%xa(:,kc_t:kc_b)
1474 :
1475 0 : sd%lnd(kc_t:kc_b) = s%lnd(kc_t:kc_b)
1476 0 : sd%Rho(kc_t:kc_b) = s%Rho(kc_t:kc_b)
1477 :
1478 0 : sd%lnPgas(kc_t:kc_b) = s%lnPgas(kc_t:kc_b)
1479 0 : sd%Pgas(kc_t:kc_b) = s%Pgas(kc_t:kc_b)
1480 :
1481 : ! Save face data
1482 :
1483 0 : kf_t = kc_t
1484 0 : kf_b = kc_b + 1
1485 :
1486 : if (TRACE_SAVE_MODEL) then
1487 : write(*,*) ' kf_t:', kf_t
1488 : write(*,*) ' kf_b:', kf_b
1489 : end if
1490 :
1491 0 : sd%gradL_composition_term(kf_t+1:kf_b-1) = s%gradL_composition_term(kf_t+1:kf_b-1)
1492 :
1493 : ! Store indices (used when we call restore_model_)
1494 :
1495 0 : sd%kc_t = kc_t
1496 0 : sd%kc_b = kc_b
1497 :
1498 0 : return
1499 :
1500 : end subroutine save_model_
1501 :
1502 :
1503 0 : subroutine restore_model_ (s, update_mode, sd)
1504 :
1505 : type(star_info), pointer :: s
1506 : integer, intent(inout) :: update_mode(:)
1507 : type(saved_data), intent(inout) :: sd
1508 :
1509 : logical, parameter :: TRACE_RESTORE_MODEL = .FALSE.
1510 :
1511 : integer :: kc_t
1512 : integer :: kc_b
1513 : integer :: kf_t
1514 : integer :: kf_b
1515 :
1516 : ! Restore data from previous call to save_model_
1517 :
1518 : if (TRACE_RESTORE_MODEL) then
1519 : write(*,*) ' Restoring model'
1520 : end if
1521 :
1522 : ! Restore cell data
1523 :
1524 0 : kc_t = sd%kc_t
1525 0 : kc_b = sd%kc_b
1526 :
1527 : if (TRACE_RESTORE_MODEL) then
1528 : write(*,*) ' kc_t:', kc_t
1529 : write(*,*) ' kc_b:', kc_b
1530 : end if
1531 :
1532 0 : update_mode(kc_t:kc_b) = sd%update_mode(kc_t:kc_b)
1533 :
1534 0 : s%xa(:,kc_t:kc_b) = sd%xa(:,kc_t:kc_b)
1535 :
1536 0 : s%lnd(kc_t:kc_b) = sd%lnd(kc_t:kc_b)
1537 0 : s%rho(kc_t:kc_b) = sd%rho(kc_t:kc_b)
1538 :
1539 0 : s%lnPgas(kc_t:kc_b) = sd%lnPgas(kc_t:kc_b)
1540 0 : s%Pgas(kc_t:kc_b) = sd%Pgas(kc_t:kc_b)
1541 :
1542 : ! Restore face data
1543 :
1544 0 : kf_t = kc_t
1545 0 : kf_b = kc_b + 1
1546 :
1547 : if (TRACE_RESTORE_MODEL) then
1548 : write(*,*) ' kf_t:', kf_t
1549 : write(*,*) ' kf_b:', kf_b
1550 : end if
1551 :
1552 0 : s%gradL_composition_term(kf_t+1:kf_b-1) = sd%gradL_composition_term(kf_t+1:kf_b-1)
1553 :
1554 : ! Update the model set those quantities that are not stored
1555 : ! explicitly in sd
1556 :
1557 0 : call update_model_(s, update_mode, kc_t, kc_b)
1558 :
1559 0 : end subroutine restore_model_
1560 :
1561 0 : end module conv_premix
|
