Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015-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 remove_shells
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi4, ln10, boltz_sigma, avo, kerg, msun
24 : use utils_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_remove_center_at_cell_k
30 : public :: do_remove_center_by_temperature
31 : public :: do_remove_center_by_radius_cm
32 : public :: do_remove_inner_fraction_q
33 : public :: do_remove_center_by_he4
34 : public :: do_remove_center_by_c12_o16
35 : public :: do_remove_center_by_si28
36 : public :: do_remove_center_to_reduce_co56_ni56
37 : public :: do_remove_center_by_ye
38 : public :: do_remove_center_by_entropy
39 : public :: do_remove_center_by_infall_kms
40 : public :: do_remove_center_at_inner_max_abs_v
41 : public :: do_remove_center_by_mass_gm
42 : public :: do_remove_fe_core
43 : public :: do_zero_inner_v_by_mass_gm
44 : public :: do_relax_to_star_cut
45 : public :: do_remove_surface_by_v_surf_km_s
46 : public :: do_remove_surface_by_v_surf_div_cs
47 : public :: do_remove_surface_by_v_surf_div_v_escape
48 : public :: do_remove_surface_at_cell_k
49 : public :: do_remove_surface_at_he_core_boundary
50 : public :: do_remove_surface_by_optical_depth
51 : public :: do_remove_surface_by_density
52 : public :: do_remove_surface_by_pressure
53 : public :: do_remove_surface_by_radius_cm
54 : public :: do_remove_surface_by_q
55 : public :: do_remove_surface_by_mass_gm
56 : public :: do_limit_center_logp
57 : public :: do_remove_center_by_logrho
58 : public :: do_remove_fallback
59 :
60 : contains
61 :
62 0 : subroutine do_remove_center_at_cell_k(id, k, ierr)
63 : integer, intent(in) :: id, k
64 : integer, intent(out) :: ierr
65 : type (star_info), pointer :: s
66 0 : call get_star_ptr(id, s, ierr)
67 0 : if (ierr /= 0) then
68 0 : write(*,*) 'do_remove_center_at_cell_k: get_star_ptr ierr', ierr
69 0 : return
70 : end if
71 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
72 : end subroutine do_remove_center_at_cell_k
73 :
74 :
75 0 : subroutine do_remove_center_by_temperature(id, temperature, ierr)
76 : integer, intent(in) :: id
77 : real(dp), intent(in) :: temperature
78 : integer, intent(out) :: ierr
79 : type (star_info), pointer :: s
80 : integer :: k
81 0 : call get_star_ptr(id, s, ierr)
82 0 : if (ierr /= 0) then
83 0 : write(*,*) 'do_remove_center_by_temperature: get_star_ptr ierr', ierr
84 0 : return
85 : end if
86 0 : do k=1,s% nz
87 0 : if (s% T(k) >= temperature) then
88 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
89 0 : return
90 : end if
91 : end do
92 0 : ierr = -1
93 : end subroutine do_remove_center_by_temperature
94 :
95 :
96 0 : subroutine do_remove_center_by_he4(id, x, ierr)
97 : use chem_def, only: ihe4
98 : integer, intent(in) :: id
99 : real(dp), intent(in) :: x
100 : integer, intent(out) :: ierr
101 : type (star_info), pointer :: s
102 : integer :: k, he4
103 0 : call get_star_ptr(id, s, ierr)
104 0 : if (ierr /= 0) then
105 0 : write(*,*) 'do_remove_center_by_he4: get_star_ptr ierr', ierr
106 0 : return
107 : end if
108 0 : he4 = s% net_iso(ihe4)
109 0 : if (he4 <= 0) then
110 0 : ierr = -1
111 0 : write(*,*) 'do_remove_center_by_he4: no he4 in current net'
112 0 : return
113 : end if
114 0 : do k=1,s% nz
115 0 : if (s% xa(he4,k) >= x) then
116 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
117 0 : return
118 : end if
119 : end do
120 0 : ierr = -1
121 0 : end subroutine do_remove_center_by_he4
122 :
123 :
124 0 : subroutine do_remove_center_by_c12_o16(id, x, ierr)
125 0 : use chem_def, only: ic12, io16
126 : integer, intent(in) :: id
127 : real(dp), intent(in) :: x
128 : integer, intent(out) :: ierr
129 : type (star_info), pointer :: s
130 : integer :: k, c12, o16
131 0 : call get_star_ptr(id, s, ierr)
132 0 : if (ierr /= 0) then
133 0 : write(*,*) 'do_remove_center_by_c12_o16: get_star_ptr ierr', ierr
134 0 : return
135 : end if
136 0 : c12 = s% net_iso(ic12)
137 0 : if (c12 <= 0) then
138 0 : ierr = -1
139 0 : write(*,*) 'do_remove_center_by_c12_o16: no c12 in current net'
140 0 : return
141 : end if
142 0 : o16 = s% net_iso(io16)
143 0 : if (o16 <= 0) then
144 0 : ierr = -1
145 0 : write(*,*) 'do_remove_center_by_c12_o16: no o16 in current net'
146 0 : return
147 : end if
148 0 : do k=1,s% nz
149 0 : if (s% xa(c12,k) + s% xa(o16,k) >= x) then
150 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
151 0 : return
152 : end if
153 : end do
154 0 : ierr = -1
155 0 : end subroutine do_remove_center_by_c12_o16
156 :
157 :
158 0 : subroutine do_remove_center_by_si28(id, x, ierr)
159 0 : use chem_def, only: isi28
160 : integer, intent(in) :: id
161 : real(dp), intent(in) :: x
162 : integer, intent(out) :: ierr
163 : type (star_info), pointer :: s
164 : integer :: k, si28
165 0 : call get_star_ptr(id, s, ierr)
166 0 : if (ierr /= 0) then
167 0 : write(*,*) 'do_remove_center_by_si28: get_star_ptr ierr', ierr
168 0 : return
169 : end if
170 0 : si28 = s% net_iso(isi28)
171 0 : if (si28 <= 0) then
172 0 : ierr = -1
173 0 : write(*,*) 'do_remove_center_by_si28: no si28 in current net'
174 0 : return
175 : end if
176 0 : do k=1,s% nz
177 0 : if (s% xa(si28,k) >= x) then
178 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
179 0 : return
180 : end if
181 : end do
182 0 : ierr = -1
183 0 : end subroutine do_remove_center_by_si28
184 :
185 :
186 0 : subroutine do_remove_center_to_reduce_co56_ni56(id, x, ierr)
187 0 : use chem_def, only: ico56, ini56
188 : integer, intent(in) :: id
189 : real(dp), intent(in) :: x
190 : integer, intent(out) :: ierr
191 : type (star_info), pointer :: s
192 : integer :: j, jj, k, species, nz, co56, ni56
193 0 : real(dp) :: mtotal, dm56, alfa_co56, dm56_new, dm56_old, x56_new, x56_old
194 : include 'formats'
195 0 : call get_star_ptr(id, s, ierr)
196 0 : if (ierr /= 0) then
197 0 : write(*,*) 'do_remove_center_to_reduce_co56_ni56: get_star_ptr ierr', ierr
198 0 : return
199 : end if
200 0 : nz = s% nz
201 0 : species = s% species
202 0 : co56 = s% net_iso(ico56)
203 0 : ni56 = s% net_iso(ini56)
204 0 : if (co56 <= 0 .or. ni56 <= 0) then
205 0 : ierr = -1
206 0 : write(*,*) 'do_remove_center_to_reduce_co56_ni56: must have both in current net'
207 0 : return
208 : end if
209 0 : mtotal = 0d0
210 0 : do k=1,nz
211 0 : dm56 = s% dm(k)*(s% xa(co56,k) + s% xa(ni56,k))/Msun
212 0 : if (mtotal+dm56 >= x) then
213 0 : write(*,2) 'mtotal+dm56 >= x', k, mtotal+dm56, x, mtotal
214 0 : if (mtotal+dm56 - x < x - mtotal) then
215 0 : if (k < nz) call do_remove_inner_fraction_q(id, s% q(k+1), ierr)
216 : else
217 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
218 : end if
219 0 : if (ierr == 0) then ! adjust ni+co in center zone
220 0 : nz = s% nz
221 : mtotal = dot_product(s% dm(1:nz), &
222 0 : s% xa(co56,1:nz) + s% xa(ni56,1:nz))/Msun
223 0 : write(*,1) 'mtotal after remove', mtotal
224 0 : write(*,2) 'nz after remove', nz
225 0 : dm56 = x - mtotal ! change Ni+Co in nz by this much
226 0 : x56_old = s% xa(co56,nz) + s% xa(ni56,nz)
227 0 : dm56_old = s% dm(nz)*x56_old/Msun
228 0 : if (x56_old <= 0d0) then
229 : alfa_co56 = 0d0
230 : else
231 0 : alfa_co56 = s% xa(co56,nz)/x56_old
232 : end if
233 0 : dm56_new = dm56_old + dm56
234 0 : x56_new = dm56_new/(s% dm(nz)/Msun)
235 0 : write(*,3) 'old x co56', co56, species, s% xa(co56,nz)
236 0 : write(*,3) 'old x ni56', ni56, species, s% xa(ni56,nz)
237 :
238 0 : s% xa(co56,nz) = 0d0
239 0 : s% xa(ni56,nz) = 0d0
240 0 : j = 1
241 0 : do jj=2,species
242 0 : if (s% xa(jj,nz) > s% xa(j,nz)) j = jj
243 : end do
244 :
245 0 : s% xa(co56,nz) = x56_new*alfa_co56
246 0 : s% xa(ni56,nz) = x56_new*(1d0 - alfa_co56)
247 0 : s% xa(j,nz) = s% xa(j,nz) - (x56_new - x56_old)
248 :
249 0 : write(*,1) 'new s% xa(co56,nz)', s% xa(co56,nz)
250 0 : write(*,1) 'new s% xa(ni56,nz)', s% xa(ni56,nz)
251 : mtotal = dot_product(s% dm(1:nz), &
252 0 : s% xa(co56,1:nz) + s% xa(ni56,1:nz))/Msun
253 0 : write(*,1) 'final mass for Ni56+Co56', mtotal
254 0 : write(*,'(A)')
255 0 : call mesa_error(__FILE__,__LINE__,'do_remove_center_to_reduce_co56_ni56')
256 : end if
257 0 : return
258 : end if
259 0 : mtotal = mtotal + dm56
260 0 : write(*,2) 'mtotal', k, mtotal
261 : end do
262 0 : write(*,1) 'failed in do_remove_center_to_reduce_co56_ni56'
263 0 : write(*,1) 'not enough Ni56+Co56 in model', mtotal, x
264 0 : ierr = -1
265 0 : end subroutine do_remove_center_to_reduce_co56_ni56
266 :
267 :
268 0 : subroutine do_remove_fallback(id, ierr)
269 : integer, intent(in) :: id
270 : integer, intent(out) :: ierr
271 : type (star_info), pointer :: s
272 : integer :: k, k0, nz
273 0 : real(dp) :: ie, ke, pe, rR, rL, rC, m_cntr, &
274 0 : sum_total_energy, speed_limit
275 0 : real(dp), pointer :: v(:)
276 :
277 : include 'formats'
278 :
279 0 : call get_star_ptr(id, s, ierr)
280 0 : if (ierr /= 0) then
281 0 : write(*,*) 'do_remove_fallback: get_star_ptr ierr', ierr
282 0 : return
283 : end if
284 :
285 0 : if (s% u_flag) then
286 0 : v => s% u
287 0 : else if (s% v_flag) then
288 0 : v => s% v
289 : else
290 : return
291 : end if
292 :
293 0 : nz = s% nz
294 :
295 : ! check to see how far extend fallback above innermost cell
296 0 : k0 = nz
297 0 : if (s% job% fallback_check_total_energy) then ! remove_bound_inner_region
298 : ! integrate total energy outward looking for sign going negative.
299 : ! if find, then continue until reach minimum integral and cut there.
300 0 : sum_total_energy = 0d0
301 0 : do k = nz,1,-1
302 0 : ie = s% energy(k)*s% dm(k)
303 0 : ke = 0.5d0*v(k)*v(k)*s% dm(k)
304 0 : if (k == s% nz) then
305 0 : rL = s% R_center
306 : else
307 0 : rL = s% r(k+1)
308 : end if
309 0 : rR = s% r(k)
310 0 : rC = 0.5d0*(rR + rL)
311 0 : m_cntr = s% m(k) - 0.5d0*s% dm(k)
312 0 : pe = -s% cgrav(k)*m_cntr*s% dm(k)/rC
313 0 : if (is_bad(ie + ke + pe)) then
314 0 : write(*,2) 'ie', k, ie
315 0 : write(*,2) 'ke', k, ke
316 0 : write(*,2) 'pe', k, pe
317 0 : call mesa_error(__FILE__,__LINE__,'do_remove_fallback')
318 : end if
319 0 : sum_total_energy = sum_total_energy + ie + ke + pe
320 0 : if (is_bad(sum_total_energy)) then
321 0 : write(*,2) 'sum_total_energy', k, sum_total_energy
322 0 : write(*,2) 'ie', k, ie
323 0 : write(*,2) 'ke', k, ke
324 0 : write(*,2) 'pe', k, pe
325 0 : call mesa_error(__FILE__,__LINE__,'do_remove_fallback')
326 : end if
327 : !write(*,3) 'sum_total_energy', k, nz, sum_total_energy, ie + ke + pe
328 0 : if (sum_total_energy < 0d0) exit
329 0 : k0 = k
330 : end do
331 0 : if (sum_total_energy >= 0d0) then
332 : !write(*,1) 'no bound inner region', sum_total_energy
333 : return ! no bound inner region
334 : end if
335 0 : do k=k0-1,1,-1
336 0 : ie = s% energy(k)*s% dm(k)
337 0 : ke = 0.5d0*v(k)*v(k)*s% dm(k)
338 0 : rL = s% r(k+1)
339 0 : rR = s% r(k)
340 0 : rC = 0.5d0*(rR + rL)
341 0 : m_cntr = s% m(k) - 0.5d0*s% dm(k)
342 0 : pe = -s% cgrav(k)*m_cntr*s% dm(k)/rC
343 0 : if (ie + ke + pe > 0d0) then ! now back to unbound cell
344 : k0 = k+1
345 : !write(*,2) 'top', k0, ie + ke + pe
346 : exit
347 : end if
348 : end do
349 : else
350 0 : speed_limit = s% job% remove_fallback_speed_limit
351 0 : if (-v(nz) <= speed_limit*s% csound(nz)) return
352 0 : do k = nz-1,1,-1
353 0 : k0 = k+1
354 0 : if (-v(k) < speed_limit*s% csound(k)) exit
355 : end do
356 : end if
357 :
358 : !call mesa_error(__FILE__,__LINE__,'do_remove_fallback')
359 :
360 : !write(*,3) 'k0 old nz', k0, s% nz, s% m(k0)/Msun
361 :
362 : ! remove cells k0..nz
363 0 : call do_remove_inner_fraction_q(id, s% q(k0), ierr)
364 :
365 0 : if (s% job% remove_center_set_zero_v_center) s% v_center = 0d0
366 :
367 0 : end subroutine do_remove_fallback
368 :
369 :
370 0 : subroutine do_remove_center_by_logRho(id, logRho_limit, ierr)
371 : integer, intent(in) :: id
372 : real(dp), intent(in) :: logRho_limit
373 : integer, intent(out) :: ierr
374 : type (star_info), pointer :: s
375 : integer :: k, k0
376 0 : real(dp) :: lnd_limit
377 : include 'formats'
378 :
379 0 : call get_star_ptr(id, s, ierr)
380 0 : if (ierr /= 0) then
381 0 : write(*,*) 'do_remove_center_by_logRho: get_star_ptr ierr', ierr
382 0 : return
383 : end if
384 :
385 0 : lnd_limit = logRho_limit*ln10
386 0 : k0 = 0
387 0 : do k = s% nz,1,-1
388 0 : if (s% lnd(k) < lnd_limit) then
389 : k0 = k
390 : exit
391 : end if
392 0 : if (s% q(k) > 0.01d0) return
393 : end do
394 :
395 : ! k0 is innermost cell with density below limit
396 : ! search out from there for outermost with density too low
397 :
398 0 : do k = k0,1,-1
399 0 : if (s% lnd(k) < lnd_limit) cycle
400 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
401 0 : return
402 : end do
403 :
404 : end subroutine do_remove_center_by_logRho
405 :
406 :
407 0 : subroutine do_limit_center_logP(id, logP_limit, ierr)
408 : integer, intent(in) :: id
409 : real(dp), intent(in) :: logP_limit
410 : integer, intent(out) :: ierr
411 : type (star_info), pointer :: s
412 : integer :: k
413 0 : real(dp) :: lnP_limit
414 : include 'formats'
415 0 : call get_star_ptr(id, s, ierr)
416 0 : if (ierr /= 0) then
417 0 : write(*,*) 'do_limit_center_logP: get_star_ptr ierr', ierr
418 0 : return
419 : end if
420 0 : lnP_limit = logP_limit*ln10
421 0 : k = s% nz
422 0 : if (s% lnPeos(k) > lnP_limit) then
423 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
424 0 : if (ierr == 0) &
425 0 : write(*,3)' remove center for logP_limit', &
426 0 : k, s% nz, s% m(k)/Msun, s% lnPeos(k)/ln10, logP_limit
427 : end if
428 : end subroutine do_limit_center_logP
429 :
430 :
431 0 : subroutine do_remove_center_by_ye(id, ye, ierr)
432 : integer, intent(in) :: id
433 : real(dp), intent(in) :: ye
434 : integer, intent(out) :: ierr
435 : type (star_info), pointer :: s
436 : integer :: k
437 0 : call get_star_ptr(id, s, ierr)
438 0 : if (ierr /= 0) then
439 0 : write(*,*) 'do_remove_center_by_ye: get_star_ptr ierr', ierr
440 0 : return
441 : end if
442 0 : do k=1,s% nz
443 0 : if (s% ye(k) <= ye) then
444 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
445 0 : return
446 : end if
447 : end do
448 0 : ierr = -1
449 : end subroutine do_remove_center_by_ye
450 :
451 :
452 0 : subroutine do_remove_center_by_entropy(id, entropy, ierr)
453 : integer, intent(in) :: id
454 : real(dp), intent(in) :: entropy
455 : integer, intent(out) :: ierr
456 : type (star_info), pointer :: s
457 : integer :: k
458 0 : call get_star_ptr(id, s, ierr)
459 0 : if (ierr /= 0) then
460 0 : write(*,*) 'do_remove_center_by_entropy: get_star_ptr ierr', ierr
461 0 : return
462 : end if
463 0 : do k=1,s% nz
464 0 : if (s% entropy(k) <= entropy) then
465 0 : call do_remove_inner_fraction_q(id, s% q(k), ierr)
466 0 : return
467 : end if
468 : end do
469 0 : ierr = -1
470 : end subroutine do_remove_center_by_entropy
471 :
472 :
473 0 : subroutine do_remove_center_by_infall_kms(id, infall_kms, ierr)
474 : integer, intent(in) :: id
475 : real(dp), intent(in) :: infall_kms
476 : integer, intent(out) :: ierr
477 : type (star_info), pointer :: s
478 : integer :: k, k_infall
479 0 : real(dp) :: v_infall, v
480 : include 'formats'
481 0 : call get_star_ptr(id, s, ierr)
482 0 : if (ierr /= 0) then
483 0 : write(*,*) 'do_remove_center_by_infall_kms: get_star_ptr ierr', ierr
484 0 : return
485 : end if
486 0 : v_infall = -abs(infall_kms)*1d5
487 0 : k_infall = 0
488 0 : do k=s% nz,1,-1
489 0 : if (s% v_flag) then
490 0 : v = s% v(k)
491 0 : else if (s% u_flag) then
492 0 : v = s% u(k)
493 : else
494 0 : write(*,*) 'must have v or u for do_remove_center_by_infall_kms'
495 0 : ierr = -1
496 0 : return
497 : end if
498 0 : if (v > v_infall) exit ! not falling fast enough
499 0 : k_infall = k
500 : end do
501 0 : if (k_infall == 0) return ! no infall
502 0 : call do_remove_inner_fraction_q(id, s% q(k_infall), ierr)
503 0 : write(*,1) 'new inner boundary mass', s% m_center/Msun
504 : end subroutine do_remove_center_by_infall_kms
505 :
506 :
507 0 : subroutine do_remove_center_by_radius_cm(id, r, ierr)
508 : integer, intent(in) :: id
509 : real(dp), intent(in) :: r
510 : integer, intent(out) :: ierr
511 : type (star_info), pointer :: s
512 0 : real(dp) :: q_r, rp1, r00, qp1
513 : integer :: k
514 :
515 0 : call get_star_ptr(id, s, ierr)
516 0 : if (ierr /= 0) then
517 0 : write(*,*) 'do_remove_center_by_radius_cm: get_star_ptr ierr', ierr
518 0 : return
519 : end if
520 0 : rp1 = s% R_center
521 0 : if (rp1 > r) then
522 0 : ierr = -1
523 0 : s% retry_message = 'error in remove center by radius: r < R_center'
524 0 : if (s% report_ierr) write(*, *) s% retry_message
525 : end if
526 0 : if (s% r(1) <= r) then
527 0 : ierr = -1
528 0 : s% retry_message = 'error in remove center by radius: r >= R_surface'
529 0 : if (s% report_ierr) write(*, *) s% retry_message
530 : end if
531 0 : if (rp1 == r) return
532 0 : qp1 = 0d0
533 0 : do k=s% nz, 1, -1
534 0 : r00 = s% r(k)
535 0 : if (r00 > r .and. r >= rp1) then
536 : q_r = qp1 + s% dq(k)* &
537 0 : (r*r*r - rp1*rp1*rp1)/(r00*r00*r00 - rp1*rp1*rp1)
538 0 : exit
539 : end if
540 0 : rp1 = r00
541 0 : qp1 = s% q(k)
542 : end do
543 0 : call do_remove_inner_fraction_q(id, q_r, ierr)
544 :
545 : end subroutine do_remove_center_by_radius_cm
546 :
547 :
548 0 : subroutine do_remove_center_at_inner_max_abs_v(id, ierr)
549 : integer, intent(in) :: id
550 : integer, intent(out) :: ierr
551 : type (star_info), pointer :: s
552 0 : real(dp) :: q_max, abs_v_max
553 : integer :: k_max
554 0 : real(dp), pointer :: v(:)
555 : include 'formats'
556 0 : call get_star_ptr(id, s, ierr)
557 0 : if (ierr /= 0) then
558 0 : write(*,*) 'do_remove_center_at_inner_max_abs_v: get_star_ptr ierr', ierr
559 0 : return
560 : end if
561 0 : if (s% u_flag) then
562 0 : v => s% u
563 0 : else if (s% v_flag) then
564 0 : v => s% v
565 : else
566 0 : call mesa_error(__FILE__,__LINE__,'no u or v for do_remove_center_at_inner_max_abs_v?')
567 0 : return
568 : end if
569 0 : k_max = minloc(v(1:s% nz),dim=1)
570 0 : q_max = s% q(k_max)
571 0 : abs_v_max = abs(v(k_max))
572 : !write(*,2) 'v abs_v_max q_max', k_max, v(k_max), abs_v_max, q_max
573 0 : call do_remove_center(id, k_max, ierr)
574 0 : end subroutine do_remove_center_at_inner_max_abs_v
575 :
576 :
577 0 : subroutine do_remove_fe_core(id, ierr)
578 : integer, intent(in) :: id
579 : integer, intent(out) :: ierr
580 : type (star_info), pointer :: s
581 0 : call get_star_ptr(id, s, ierr)
582 0 : if (ierr /= 0) then
583 0 : write(*,*) 'do_remove_fe_core: get_star_ptr ierr', ierr
584 0 : return
585 : end if
586 0 : if (s% fe_core_k <= 0 .or. s% fe_core_k > s% nz) return
587 0 : call do_remove_inner_fraction_q(id, s% q(s% fe_core_k), ierr)
588 : end subroutine do_remove_fe_core
589 :
590 :
591 0 : subroutine do_remove_center_by_mass_gm(id, m, ierr)
592 : integer, intent(in) :: id
593 : real(dp), intent(in) :: m
594 : integer, intent(out) :: ierr
595 : type (star_info), pointer :: s
596 : real(dp) :: q_m
597 : include 'formats'
598 0 : call get_star_ptr(id, s, ierr)
599 0 : if (ierr /= 0) then
600 0 : write(*,*) 'do_remove_center_by_mass_gm: get_star_ptr ierr', ierr
601 0 : return
602 : end if
603 0 : q_m = (m - s% M_center)/s% xmstar
604 0 : if (q_m <= 0d0) return
605 0 : call do_remove_inner_fraction_q(id, q_m, ierr)
606 : end subroutine do_remove_center_by_mass_gm
607 :
608 :
609 0 : subroutine do_remove_inner_fraction_q(id, q, ierr)
610 : ! note: we have not implemented fractional cell removal.
611 : ! it adds a lot of complexity and, so far, we don't need it.
612 : use alloc, only: prune_star_info_arrays
613 : use star_utils, only: set_qs
614 : integer, intent(in) :: id
615 : real(dp), intent(in) :: q
616 : integer, intent(out) :: ierr
617 : type (star_info), pointer :: s
618 : integer :: k
619 : include 'formats'
620 0 : call get_star_ptr(id, s, ierr)
621 0 : if (ierr /= 0) then
622 0 : write(*,*) 'do_remove_inner_fraction_q: get_star_ptr ierr', ierr
623 0 : return
624 : end if
625 0 : if (q < 0d0 .or. q > 1d0) then
626 0 : ierr = -1
627 0 : s% retry_message = 'error in remove center: invalid location q'
628 0 : if (s% report_ierr) write(*, *) s% retry_message
629 0 : return
630 : end if
631 0 : do k = 1, s% nz
632 0 : if (s% q(k) <= q) exit
633 : end do
634 0 : call do_remove_center(id, k, ierr)
635 0 : end subroutine do_remove_inner_fraction_q
636 :
637 :
638 0 : subroutine do_remove_center(id, k, ierr)
639 0 : use read_model, only: finish_load_model
640 : use alloc, only: prune_star_info_arrays
641 : use star_utils, only: normalize_dqs, set_qs
642 : integer, intent(in) :: id, k
643 : integer, intent(out) :: ierr
644 : type (star_info), pointer :: s
645 0 : real(dp) :: old_xmstar, new_xmstar
646 : logical :: restart
647 : integer :: kk
648 : include 'formats'
649 0 : call get_star_ptr(id, s, ierr)
650 0 : if (ierr /= 0) then
651 0 : write(*,*) 'do_remove_center: get_star_ptr ierr', ierr
652 0 : return
653 : end if
654 0 : if (k <= 1 .or. k > s% nz) return
655 0 : old_xmstar = s% xmstar
656 0 : s% M_center = s% m(k)
657 0 : new_xmstar = s% m(1) - s% M_center
658 0 : s% xmstar = new_xmstar
659 0 : s% R_center = s% r(k)
660 0 : if (s% job% remove_center_adjust_L_center) s% L_center = s% L(k)
661 0 : if (s% u_flag) then
662 0 : kk = minloc(s% u(1:s% nz),dim=1)
663 0 : s% v_center = s% u(k)
664 0 : if (is_bad(s% v_center)) then
665 0 : write(*,2) 'center s% u(k)', k, s% u(k)
666 0 : call mesa_error(__FILE__,__LINE__,'do_remove_center')
667 : end if
668 0 : else if (s% v_flag) then
669 0 : s% v_center = s% v(k)
670 0 : if (is_bad(s% v_center)) then
671 0 : write(*,2) 'center s% v(k)', k, s% v(k)
672 0 : call mesa_error(__FILE__,__LINE__,'do_remove_center')
673 : end if
674 : else
675 0 : s% v_center = 0d0
676 : end if
677 0 : if (s% job% remove_center_set_zero_v_center) s% v_center = 0d0
678 0 : s% nz = k-1
679 0 : do kk=1,k-1
680 0 : s% dq(kk) = s% dm(kk)/new_xmstar
681 : end do
682 0 : if (.not. s% do_normalize_dqs_as_part_of_set_qs) then
683 0 : call normalize_dqs(s, s% nz, s% dq, ierr)
684 0 : if (ierr /= 0) then
685 0 : if (s% report_ierr) write(*,*) 'normalize_dqs failed in do_remove_center'
686 0 : return
687 : end if
688 : end if
689 0 : call set_qs(s, s% nz, s% q, s% dq, ierr)
690 0 : if (ierr /= 0) return
691 0 : s% generations = 1 ! memory leak, but hopefully not necessary to fix
692 : ! assuming remove center is a rare operation
693 0 : call prune_star_info_arrays(s, ierr)
694 0 : if (ierr /= 0) return
695 0 : s% need_to_setvars = .true.
696 0 : restart = .false.
697 0 : call finish_load_model(s, restart, ierr)
698 0 : end subroutine do_remove_center
699 :
700 :
701 0 : subroutine do_zero_inner_v_by_mass_gm(id, m, ierr)
702 : integer, intent(in) :: id
703 : real(dp), intent(in) :: m
704 : integer, intent(out) :: ierr
705 : type (star_info), pointer :: s
706 : integer :: k, i_u, i_v
707 : include 'formats'
708 0 : call get_star_ptr(id, s, ierr)
709 0 : if (ierr /= 0) then
710 0 : write(*,*) 'do_zero_inner_v_by_mass_gm: get_star_ptr ierr', ierr
711 0 : return
712 : end if
713 0 : i_u = s% i_u
714 0 : i_v = s% i_v
715 0 : do k=s% nz, 1, -1
716 0 : if (s% m(k) > m) exit
717 0 : if (i_u > 0) then
718 0 : s% xh(i_u, k) = 0d0
719 0 : s% u(k) = 0d0
720 : end if
721 0 : if (i_v > 0) then
722 0 : s% xh(i_v, k) = 0d0
723 0 : s% v(k) = 0d0
724 : end if
725 : end do
726 0 : end subroutine do_zero_inner_v_by_mass_gm
727 :
728 :
729 0 : subroutine do_remove_surface_at_cell_k(id, k, ierr)
730 : integer, intent(in) :: id, k
731 : integer, intent(out) :: ierr
732 : type (star_info), pointer :: s
733 0 : call get_star_ptr(id, s, ierr)
734 0 : if (ierr /= 0) then
735 0 : write(*,*) 'do_remove_surface_at_cell_k: get_star_ptr ierr', ierr
736 0 : return
737 : end if
738 0 : if (k <= 1) return
739 0 : call do_remove_surface(id, k, ierr)
740 : end subroutine do_remove_surface_at_cell_k
741 :
742 :
743 0 : subroutine do_remove_surface_at_he_core_boundary(id, h1_fraction, ierr)
744 : integer, intent(in) :: id
745 : real(dp), intent(in) :: h1_fraction
746 : integer, intent(out) :: ierr
747 : type (star_info), pointer :: s
748 : integer :: k
749 0 : call get_star_ptr(id, s, ierr)
750 0 : if (ierr /= 0) then
751 0 : write(*,*) 'do_remove_surface_at_he_core_boundary: get_star_ptr ierr', ierr
752 0 : return
753 : end if
754 0 : if (s% X(1) <= h1_fraction) return
755 0 : do k=2,s% nz
756 0 : if (s% X(k) <= h1_fraction) then
757 0 : call do_remove_surface(id, k-1, ierr)
758 0 : return
759 : end if
760 : end do
761 0 : ierr = -1
762 : end subroutine do_remove_surface_at_he_core_boundary
763 :
764 :
765 0 : subroutine do_remove_surface_by_optical_depth(id, optical_depth, ierr)
766 : integer, intent(in) :: id
767 : real(dp), intent(in) :: optical_depth
768 : integer, intent(out) :: ierr
769 : type (star_info), pointer :: s
770 : integer :: k
771 0 : call get_star_ptr(id, s, ierr)
772 0 : if (ierr /= 0) then
773 0 : write(*,*) 'do_remove_surface_by_optical_depth: get_star_ptr ierr', ierr
774 0 : return
775 : end if
776 0 : if (optical_depth <= s% tau(1)) return
777 0 : do k=1,s% nz
778 0 : if (s% tau(k) >= optical_depth) then
779 0 : call do_remove_surface(id, k, ierr)
780 0 : return
781 : end if
782 : end do
783 0 : ierr = -1
784 : end subroutine do_remove_surface_by_optical_depth
785 :
786 :
787 0 : subroutine do_remove_surface_by_v_surf_km_s(id, v_surf_km_s, ierr)
788 : integer, intent(in) :: id
789 : real(dp), intent(in) :: v_surf_km_s
790 : integer, intent(out) :: ierr
791 : type (star_info), pointer :: s
792 : integer :: k
793 0 : real(dp), dimension(:), pointer :: v
794 0 : real(dp) :: v_max
795 : include 'formats'
796 0 : call get_star_ptr(id, s, ierr)
797 0 : if (ierr /= 0) then
798 0 : write(*,*) 'do_remove_surface_by_v_surf_km_s: get_star_ptr ierr', ierr
799 0 : return
800 : end if
801 0 : if (s% u_flag) then
802 0 : v => s% u
803 0 : else if (s% v_flag) then
804 0 : v => s% v
805 : else
806 : return
807 : end if
808 0 : v_max = 1d5*v_surf_km_s
809 0 : if (v(1) < v_max) return
810 0 : do k=2,3 ! s% nz
811 0 : if (v(k) < v_max) exit
812 0 : write(*,2) 'v', k-1, v(k-1)/1d5, v_surf_km_s
813 : end do
814 0 : write(*,2) 'do_remove_surface_by_v_surf_km_s', k-1, v(k-1)/1d5, v_surf_km_s
815 0 : call do_remove_surface(id, k-1, ierr)
816 0 : return
817 0 : end subroutine do_remove_surface_by_v_surf_km_s
818 :
819 :
820 0 : subroutine do_remove_surface_by_v_surf_div_cs(id, v_surf_div_cs, ierr)
821 : integer, intent(in) :: id
822 : real(dp), intent(in) :: v_surf_div_cs
823 : integer, intent(out) :: ierr
824 : type (star_info), pointer :: s
825 : integer :: k
826 0 : real(dp), dimension(:), pointer :: v
827 : include 'formats'
828 0 : call get_star_ptr(id, s, ierr)
829 0 : if (ierr /= 0) then
830 0 : write(*,*) 'do_remove_surface_by_v_surf_div_cs: get_star_ptr ierr', ierr
831 0 : return
832 : end if
833 0 : if (s% u_flag) then
834 0 : v => s% u
835 0 : else if (s% v_flag) then
836 0 : v => s% v
837 : else
838 : return
839 : end if
840 : !write(*,1) 'v(1)/cs', v(1)/s% csound(1), v_surf_div_cs
841 0 : if (v(1) < s% csound(1)*v_surf_div_cs) return
842 0 : do k=2,30 ! s% nz
843 0 : if (v(k) < s% csound(k)*v_surf_div_cs) exit
844 0 : write(*,2) 'v/cs', k-1, v(k-1)/s% csound(k-1)
845 : end do
846 0 : write(*,2) 'do_remove_surface_by_v_surf_div_cs', k-1, v(k-1)/s% csound(k-1)
847 0 : call do_remove_surface(id, k-1, ierr)
848 : !call mesa_error(__FILE__,__LINE__,'do_remove_surface_by_v_surf_div_cs')
849 0 : return
850 0 : end subroutine do_remove_surface_by_v_surf_div_cs
851 :
852 :
853 0 : subroutine do_remove_surface_by_v_surf_div_v_escape(id, v_surf_div_v_escape, ierr)
854 : integer, intent(in) :: id
855 : real(dp), intent(in) :: v_surf_div_v_escape
856 : integer, intent(out) :: ierr
857 : type (star_info), pointer :: s
858 : integer :: k, k_vesc
859 0 : real(dp) :: vesc
860 0 : real(dp), dimension(:), pointer :: v
861 : include 'formats'
862 0 : call get_star_ptr(id, s, ierr)
863 0 : if (ierr /= 0) then
864 0 : write(*,*) 'do_remove_surface_by_v_surf_div_v_escape: get_star_ptr ierr', ierr
865 0 : return
866 : end if
867 0 : if (s% u_flag) then
868 0 : v => s% u
869 0 : else if (s% v_flag) then
870 0 : v => s% v
871 : else
872 : return
873 : end if
874 0 : k_vesc = 0
875 0 : do k=2, s% nz
876 0 : if (s% q(k) > s% job% max_q_for_remove_surface_by_v_surf_div_v_escape) cycle
877 0 : if (s% q(k) < s% job% min_q_for_remove_surface_by_v_surf_div_v_escape) exit
878 0 : vesc = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k)))
879 0 : if (v(k) >= vesc*v_surf_div_v_escape) k_vesc = k
880 : end do
881 0 : if (k_vesc == 0) return
882 0 : write(*,2) 'do_remove_surface_by_v_surf_div_v_escape q', k_vesc, s% q(k_vesc)
883 0 : call do_remove_surface(id, k_vesc, ierr)
884 0 : end subroutine do_remove_surface_by_v_surf_div_v_escape
885 :
886 :
887 0 : subroutine do_remove_surface_by_pressure(id, pressure, ierr)
888 : integer, intent(in) :: id
889 : real(dp), intent(in) :: pressure
890 : integer, intent(out) :: ierr
891 : type (star_info), pointer :: s
892 : integer :: k
893 0 : call get_star_ptr(id, s, ierr)
894 0 : if (ierr /= 0) then
895 0 : write(*,*) 'do_remove_surface_by_pressure: get_star_ptr ierr', ierr
896 0 : return
897 : end if
898 0 : if (pressure <= s% Peos(1)) return
899 0 : do k=1,s% nz
900 0 : if (s% Peos(k) >= pressure) then
901 0 : call do_remove_surface(id, k, ierr)
902 0 : return
903 : end if
904 : end do
905 0 : ierr = -1
906 : end subroutine do_remove_surface_by_pressure
907 :
908 :
909 0 : subroutine do_remove_surface_by_density(id, density, ierr)
910 : integer, intent(in) :: id
911 : real(dp), intent(in) :: density
912 : integer, intent(out) :: ierr
913 : type (star_info), pointer :: s
914 : ierr = 0
915 : include 'formats'
916 0 : call get_star_ptr(id, s, ierr)
917 0 : if (ierr /= 0) then
918 0 : write(*,*) 'do_remove_surface_by_density: get_star_ptr ierr', ierr
919 0 : return
920 : end if
921 0 : if (s% rho(1) < density) then
922 0 : write(*,2) 'do_remove_surface', 1, s% rho(1), density
923 0 : call do_remove_surface(id, 2, ierr)
924 0 : return
925 : end if
926 : end subroutine do_remove_surface_by_density
927 :
928 :
929 0 : subroutine do_remove_surface_by_radius_cm(id, r, ierr)
930 : integer, intent(in) :: id
931 : real(dp), intent(in) :: r
932 : integer, intent(out) :: ierr
933 : type (star_info), pointer :: s
934 : integer :: k
935 0 : call get_star_ptr(id, s, ierr)
936 0 : if (ierr /= 0) then
937 0 : write(*,*) 'do_remove_surface_by_radius_cm: get_star_ptr ierr', ierr
938 0 : return
939 : end if
940 0 : if (r >= s% r(1)) return
941 0 : do k=1,s% nz
942 0 : if (s% r(k) <= r) then
943 0 : call do_remove_surface(id, k, ierr)
944 0 : return
945 : end if
946 : end do
947 0 : ierr = -1
948 : end subroutine do_remove_surface_by_radius_cm
949 :
950 :
951 0 : subroutine do_remove_surface_by_mass_gm(id, m, ierr)
952 : integer, intent(in) :: id
953 : real(dp), intent(in) :: m
954 : integer, intent(out) :: ierr
955 : type (star_info), pointer :: s
956 : integer :: k
957 : include 'formats'
958 0 : call get_star_ptr(id, s, ierr)
959 0 : if (ierr /= 0) then
960 0 : write(*,*) 'do_remove_surface_by_mass_gm: get_star_ptr ierr', ierr
961 0 : return
962 : end if
963 0 : if (m >= s% m(1)) return
964 0 : do k=1,s% nz
965 0 : if (s% m(k) <= m) then
966 0 : call do_remove_surface(id, k, ierr)
967 0 : return
968 : end if
969 : end do
970 0 : ierr = -1
971 : end subroutine do_remove_surface_by_mass_gm
972 :
973 :
974 0 : subroutine do_remove_surface_by_q(id, q, ierr)
975 : integer, intent(in) :: id
976 : real(dp), intent(in) :: q
977 : integer, intent(out) :: ierr
978 : type (star_info), pointer :: s
979 : integer :: k
980 0 : call get_star_ptr(id, s, ierr)
981 0 : if (ierr /= 0) then
982 0 : write(*,*) 'do_remove_surface_by_q: get_star_ptr ierr', ierr
983 0 : return
984 : end if
985 0 : if (q >= 1d0) return
986 0 : do k=1,s% nz
987 0 : if (s% q(k) <= q) then
988 0 : call do_remove_surface(id, k, ierr)
989 0 : return
990 : end if
991 : end do
992 0 : ierr = -1
993 : end subroutine do_remove_surface_by_q
994 :
995 :
996 0 : subroutine do_remove_surface(id, surface_k, ierr)
997 : use read_model, only: finish_load_model
998 : use mesh_adjust, only: do_prune_mesh_surface
999 : use alloc, only: resize_star_info_arrays
1000 : use star_utils, only: tau_eff
1001 : integer, intent(in) :: id, surface_k
1002 : integer, intent(out) :: ierr
1003 : type (star_info), pointer :: s, c, prv
1004 0 : type (star_info), target :: copy_info
1005 0 : real(dp) :: tau_surf_new, tau_factor_new, Lmid, Rmid, T, P, T_black_body
1006 : integer :: k, k_old, nz, nz_old, skip
1007 :
1008 : logical, parameter :: dbg = .false., restart = .false.
1009 :
1010 : include 'formats'
1011 :
1012 0 : if (surface_k == 1) then
1013 0 : ierr = 0
1014 0 : return
1015 : end if
1016 :
1017 0 : call get_star_ptr(id, s, ierr)
1018 0 : if (ierr /= 0) then
1019 0 : if (s% report_ierr) &
1020 0 : write(*,*) 'do_remove_surface: get_star_ptr ierr'
1021 0 : return
1022 : end if
1023 :
1024 0 : if (s% job% remove_surface_by_relax_to_star_cut) then
1025 0 : if (s% R_center /= 0d0) then
1026 0 : write(*,*) 'remove surface currently requires model with inner boundary at true center of star'
1027 0 : ierr = -1
1028 0 : call mesa_error(__FILE__,__LINE__,'do_remove_surface')
1029 : end if
1030 : call do_relax_to_star_cut( &
1031 : id, surface_k, s% job% remove_surface_do_jrot, &
1032 : s% job% remove_surface_do_entropy, &
1033 0 : s% job% remove_surface_turn_off_energy_sources_and_sinks, ierr)
1034 0 : return
1035 : end if
1036 :
1037 0 : nz_old = s% nz
1038 0 : skip = surface_k - 1
1039 :
1040 : if (dbg) write(*,2) 'do remove surface skip', skip
1041 0 : if (skip < 1 .or. skip >= nz_old) return
1042 :
1043 0 : tau_surf_new = tau_eff(s,1+skip)
1044 0 : tau_factor_new = s% tau_factor*tau_surf_new/s% tau(1)
1045 :
1046 : if (dbg) write(*,1) 'tau_surf_old', s% tau(1)
1047 : if (dbg) write(*,1) 'tau_factor_old', s% tau_factor
1048 : if (dbg) write(*,1) 'tau_surf_new', tau_surf_new
1049 : if (dbg) write(*,1) 'tau_factor_new', tau_factor_new
1050 :
1051 0 : rmid = s% rmid(1+skip)
1052 0 : Lmid = (s% L(1+skip) + s% L(2+skip))/2
1053 0 : T = s% T(1+skip)
1054 0 : P = s% Peos(1+skip)
1055 :
1056 0 : if (.not. associated(s% other_star_info)) then
1057 0 : allocate(s% other_star_info)
1058 0 : prv => s% other_star_info
1059 0 : c => null()
1060 : if (dbg) write(*,1) 'c is null'
1061 : else
1062 0 : prv => s% other_star_info
1063 0 : c => copy_info
1064 0 : c = prv
1065 : end if
1066 :
1067 0 : prv = s ! this makes copies of pointers and scalars
1068 :
1069 0 : nz = nz_old - skip
1070 0 : s% nz = nz
1071 : if (dbg) write(*,2) 'nz_old', nz_old
1072 : if (dbg) write(*,2) 'nz', nz
1073 :
1074 : if (dbg) write(*,1) 'call resize_star_info_arrays'
1075 0 : call resize_star_info_arrays(s, c, ierr)
1076 0 : if (ierr /= 0) then
1077 0 : if (s% report_ierr) &
1078 0 : write(*,*) 'resize_star_info_arrays failed in do_remove_surface'
1079 0 : return
1080 : end if
1081 :
1082 : if (dbg) write(*,2) 'prv% dm(1)/Msun', 1, prv% dm(1)/Msun
1083 : if (dbg) write(*,2) 'prv% m(1)/msun', 1, prv% m(1)/Msun
1084 : if (dbg) write(*,2) 'prv% dm(nz_old)/Msun', nz_old, prv% dm(nz_old)/Msun
1085 : if (dbg) write(*,2) 'prv% m(nz_old)/msun', nz_old, prv% m(nz_old)/Msun
1086 :
1087 0 : s% mstar = prv% m(1 + skip)
1088 0 : s% xmstar = s% mstar - prv% M_center
1089 0 : s% q(1) = 1d0
1090 0 : do k = 1, nz
1091 0 : k_old = k + skip
1092 0 : s% dq(k) = prv% dm(k_old)/s% xmstar
1093 0 : if (k > 1) s% q(k) = s% q(k-1) - s% dq(k-1)
1094 : end do
1095 0 : s% dq(nz) = s% q(nz)
1096 0 : do k = 1, nz
1097 0 : s% dm(k) = s% dq(k)*s% xmstar
1098 0 : s% m(k) = s% q(k)*s% xmstar + s% M_center
1099 : end do
1100 : if (dbg) write(*,2) 's% dq(1)', 1, s% dq(1)
1101 : if (dbg) write(*,2) 's% dm(1)/Msun', 1, s% dm(1)/Msun
1102 : if (dbg) write(*,2) 's% m(1)/msun', 1, s% m(1)/Msun
1103 : if (dbg) write(*,2) 's% dm(nz)/Msun', nz, s% dm(nz)/Msun
1104 : if (dbg) write(*,2) 's% m(nz)/msun', nz, s% m(nz)/Msun
1105 :
1106 : if (dbg) write(*,1) 'call do_prune_mesh_surface'
1107 : call do_prune_mesh_surface( &
1108 : s, nz, nz_old, prv% xh, prv% xa, &
1109 : prv% j_rot, prv% i_rot, &
1110 : prv% omega, prv% D_omega, prv% am_nu_rot, &
1111 : prv% conv_vel, prv% lnT, &
1112 : prv% dPdr_dRhodr_info, prv% nu_ST, prv% D_ST, prv% D_DSI, prv% D_SH, &
1113 : prv% D_SSI, prv% D_ES, prv% D_GSF, prv% D_mix, &
1114 0 : s% xh, s% xa, ierr)
1115 0 : if (ierr /= 0) then
1116 : return
1117 : end if
1118 : if (dbg) write(*,2) 's% dq(1)', 1, s% dq(1)
1119 : if (dbg) write(*,2) 's% dm(1)/Msun', 1, s% dm(1)/Msun
1120 : if (dbg) write(*,2) 's% m(1)/msun', 1, s% m(1)/Msun
1121 : if (dbg) write(*,2) 's% dm(nz)/Msun', nz, s% dm(nz)/Msun
1122 : if (dbg) write(*,2) 's% m(nz)/msun', nz, s% m(nz)/Msun
1123 :
1124 0 : if (Lmid > 0d0) then
1125 0 : T_black_body = pow(Lmid/(pi4*rmid*rmid*boltz_sigma), 0.25d0)
1126 0 : s% Tsurf_factor = T/T_black_body
1127 : else
1128 0 : s% Tsurf_factor = 1d0
1129 : end if
1130 0 : s% force_Tsurf_factor = s% Tsurf_factor
1131 :
1132 0 : if (s% use_momentum_outer_bc) then
1133 0 : s% tau_factor = tau_factor_new
1134 0 : s% force_tau_factor = s% tau_factor
1135 : end if
1136 :
1137 0 : s% need_to_setvars = .true.
1138 : if (dbg) write(*,1) 'call finish_load_model'
1139 0 : call finish_load_model(s, restart, ierr)
1140 0 : if (ierr /= 0) then
1141 0 : if (s% report_ierr) &
1142 0 : write(*,*) 'finish_load_model failed in do_remove_surface'
1143 0 : return
1144 : end if
1145 :
1146 : if (dbg) write(*,1) 'do_remove_surface tau_factor, Tsurf_factor', &
1147 : s% tau_factor, s% Tsurf_factor
1148 :
1149 : if (dbg) call mesa_error(__FILE__,__LINE__,'do_remove_surface')
1150 :
1151 0 : end subroutine do_remove_surface
1152 :
1153 :
1154 : ! Relax to a trimmed stellar model with surface cells removed down to k_remove
1155 : ! (the cell k_remove will be the outermost in the new model).
1156 0 : subroutine do_relax_to_star_cut( &
1157 : id, k_remove, do_jrot, do_entropy, turn_off_energy_sources_and_sinks, ierr)
1158 :
1159 0 : use interp_1d_def, only: pm_work_size
1160 : use interp_1d_lib, only: interp_pm, interp_values, interp_value
1161 : use adjust_xyz, only: change_net
1162 : use set_flags, only: set_v_flag, set_u_flag, set_RTI_flag, set_rotation_flag
1163 : use rotation_mix_info, only: set_rotation_mixing_info
1164 : use hydro_rotation, only: set_i_rot, set_rotation_info
1165 : use relax, only: do_relax_composition, do_relax_angular_momentum, do_relax_entropy
1166 : use init, only: load_zams_model
1167 :
1168 : integer, intent(in) :: id, k_remove
1169 : logical, intent(in) :: do_jrot, do_entropy
1170 : logical, intent(in) :: turn_off_energy_sources_and_sinks
1171 : ! determines if we turn off non_nuc_neu and eps_nuc for entropy relax
1172 : integer, intent(out) :: ierr
1173 :
1174 : logical :: v_flag, u_flag, RTI_flag, rotation_flag
1175 : type (star_info), pointer :: s
1176 : character (len=net_name_len) :: net_name
1177 : integer :: model_number, num_trace_history_values, photo_interval
1178 0 : real(dp) :: eps_nuc_factor, non_nuc_neu_factor, &
1179 0 : initial_z, initial_y, initial_mass, &
1180 0 : cumulative_energy_error, cumulative_extra_heating
1181 :
1182 0 : real(dp), pointer :: q(:), xq(:), xa(:,:), j_rot(:), entropy(:)
1183 0 : real(dp) :: time
1184 : integer :: num_pts, k0, species
1185 : logical :: save_have_mlt_vc
1186 : logical, parameter :: dbg = .false.
1187 :
1188 : ierr = 0
1189 0 : call get_star_ptr(id, s, ierr)
1190 0 : if (ierr /= 0) return
1191 :
1192 0 : eps_nuc_factor = s% eps_nuc_factor
1193 0 : non_nuc_neu_factor = s% non_nuc_neu_factor
1194 0 : net_name = s% net_name
1195 0 : num_trace_history_values = s% num_trace_history_values
1196 :
1197 0 : time = s% time
1198 0 : model_number = s% model_number
1199 0 : num_trace_history_values = s% num_trace_history_values
1200 0 : cumulative_energy_error = s% cumulative_energy_error
1201 0 : cumulative_extra_heating = s% cumulative_extra_heating
1202 :
1203 : ! zero model_number and time (will restore later)
1204 0 : s% model_number = 0
1205 0 : s% time = 0
1206 :
1207 0 : species = s% species
1208 0 : num_pts = s% nz - k_remove + 1
1209 0 : allocate(q(num_pts), xq(num_pts), xa(species, num_pts))
1210 0 : rotation_flag = .false.
1211 0 : if (do_jrot .and. s% rotation_flag) then
1212 0 : allocate(j_rot(num_pts))
1213 0 : rotation_flag = .true.
1214 : end if
1215 0 : if (do_entropy) then
1216 0 : allocate(entropy(num_pts))
1217 : end if
1218 : !need to compute cell-centered q for remaining mass
1219 0 : xq(1) = s% dq(k_remove)/2/s% q(k_remove)
1220 0 : do k0 = 1, num_pts-1
1221 0 : xq(1+k0) = xq(1+k0-1) + (s% dq(k_remove+k0) + s% dq(k_remove+k0-1))/s% q(k_remove)/2
1222 : end do
1223 :
1224 : !!create interpolant for convective velocities
1225 : ! TODO: adapt this for TDC
1226 : !conv_vel_flag = .false.
1227 : !if (s% conv_vel_flag) then
1228 : ! conv_vel_flag = .true.
1229 : ! allocate(interp_work((num_pts)*pm_work_size), &
1230 : ! conv_vel_interp(4*(num_pts)), stat=ierr)
1231 : ! do k0 = 1, num_pts
1232 : ! conv_vel_interp(4*k0-3) = s% conv_vel(k0+k_remove-1)
1233 : ! q(k0) = s% q(k0+k_remove-1)/s% q(k_remove)
1234 : ! end do
1235 : ! call interp_pm(q, num_pts, conv_vel_interp,&
1236 : ! pm_work_size, interp_work, 'conv_vel interpolant', ierr)
1237 :
1238 : ! ! turn off conv vel flag to load model
1239 : ! call set_conv_vel_flag(id, .false., ierr)
1240 : ! if (dbg) write(*,*) "set_conv_vel_flag ierr", ierr
1241 : !end if
1242 :
1243 : ! save have_mlt_vc and set to false (to load ZAMS model)
1244 0 : save_have_mlt_vc = s% have_mlt_vc
1245 0 : s% have_mlt_vc = .false.
1246 :
1247 : !save composition and entropy profiles
1248 0 : xa(:,:) = s% xa(:,k_remove:s% nz)
1249 0 : if (rotation_flag) then
1250 0 : j_rot(:) = s% j_rot(k_remove:s% nz)
1251 : end if
1252 0 : if (do_entropy) then
1253 0 : entropy(:) = s% entropy(k_remove:s% nz)*avo*kerg
1254 : end if
1255 :
1256 : ! various flags need to be turned off for the ZAMS model to load
1257 0 : v_flag = .false.
1258 0 : if (s% v_flag) then
1259 0 : call set_v_flag(id, .false., ierr)
1260 : if (dbg) write(*,*) "set_v_flag ierr", ierr
1261 0 : v_flag = .true.
1262 : end if
1263 :
1264 0 : u_flag = .false.
1265 0 : if (s% u_flag) then
1266 0 : call set_u_flag(id, .false., ierr)
1267 : if (dbg) write(*,*) "set_u_flag ierr", ierr
1268 0 : u_flag = .true.
1269 : end if
1270 :
1271 0 : RTI_flag = .false.
1272 0 : if (s% RTI_flag) then
1273 0 : call set_RTI_flag(id, .false., ierr)
1274 : if (dbg) write(*,*) "set_RTI_flag ierr", ierr
1275 0 : RTI_flag = .true.
1276 : end if
1277 :
1278 0 : if (s% rotation_flag) then
1279 0 : call set_rotation_flag(id, .false., ierr)
1280 : if (dbg) write(*,*) "set_rotation_flag ierr", ierr
1281 : end if
1282 :
1283 : ! avoid making photos
1284 0 : photo_interval = s% photo_interval
1285 0 : s% photo_interval = 10000000
1286 0 : s% have_previous_conv_vel = .false.
1287 0 : s% have_j_rot = .false.
1288 : ! WARNING, might need to add stuff here to actually get the ZAMS model to load.
1289 : ! otherwise can get an error of the form "error in reading model data j+species > nvec"
1290 : ! if you happen to run into these problem, check for flags being checked in read1_model in read_model.f90
1291 : ! and be sure they're turned off.
1292 :
1293 : ! set values used to load the starting model that will be relaxed
1294 0 : initial_z = s% initial_z
1295 0 : initial_y = s% initial_y
1296 0 : initial_mass = s% initial_mass
1297 0 : s% initial_z = 0.02d0
1298 0 : s% initial_y = 0.28d0
1299 0 : s% initial_mass = s% m(k_remove)/Msun
1300 :
1301 0 : s% prev_mesh_nz = 0
1302 :
1303 0 : call change_net(id, .true., 'basic.net', ierr) ! TODO:need to allow specification of different net
1304 : if (dbg) write(*,*) "check change_net ierr", ierr
1305 0 : if (ierr /= 0) return
1306 0 : call load_zams_model(id, ierr)
1307 : if (dbg) write(*,*) "check load_zams ierr", ierr
1308 0 : if (ierr /= 0) return
1309 0 : call change_net(id, .true., net_name, ierr)
1310 : if (dbg) write(*,*) "check ierr", ierr
1311 0 : if (ierr /= 0) return
1312 :
1313 : ! restore have_mlt_vc
1314 0 : s% have_mlt_vc = save_have_mlt_vc
1315 :
1316 0 : if (turn_off_energy_sources_and_sinks) then
1317 0 : s% non_nuc_neu_factor = 0d0
1318 0 : s% eps_nuc_factor = 0d0
1319 : end if
1320 :
1321 0 : s% num_trace_history_values = 0
1322 : call do_relax_composition( &
1323 0 : id, s% job% num_steps_to_relax_composition, num_pts, species, xa, xq, ierr)
1324 : if (dbg) write(*,*) "check ierr", ierr
1325 0 : if (ierr /= 0) return
1326 0 : deallocate(xa)
1327 :
1328 0 : if (rotation_flag) then
1329 0 : call set_rotation_flag(id, .true., ierr)
1330 : if (dbg) write(*,*) "set_rotation_flag true ierr", ierr
1331 0 : if (ierr /= 0) return
1332 0 : call set_rotation_info(s, .false., ierr)
1333 : if (dbg) write(*,*) "set_rotation_info ierr", ierr
1334 0 : if (ierr /= 0) return
1335 0 : call set_rotation_mixing_info(s, ierr)
1336 : if (dbg) write(*,*) "set_rotation_mixing_info ierr", ierr
1337 0 : if (ierr /= 0) return
1338 : call do_relax_angular_momentum( &
1339 0 : id, s% job% max_steps_to_relax_angular_momentum, num_pts, j_rot, xq, ierr)
1340 : if (dbg) write(*,*) "check ierr", ierr
1341 0 : if (ierr /= 0) return
1342 0 : deallocate(j_rot)
1343 : end if
1344 :
1345 0 : if (do_entropy) then
1346 : call do_relax_entropy( &
1347 0 : id, s% job% max_steps_to_relax_entropy, num_pts, entropy, xq, ierr)
1348 : if (dbg) write(*,*) "check ierr", ierr
1349 0 : if (ierr /= 0) return
1350 0 : deallocate(entropy)
1351 : end if
1352 :
1353 : !!take care of convective velocities
1354 : ! TODO: adapt for TDC
1355 : !if (s% conv_vel_flag) then
1356 : ! do k0=1, s% nz
1357 : ! call interp_value(q, num_pts, conv_vel_interp, s% q(k0), s% conv_vel(k0), ierr)
1358 : ! !avoid extending regions with non-zero conv vel
1359 : ! do k=2, num_pts-1
1360 : ! if(s% q(k0) < q(k) .and. s% q(k0) > q(k+1) &
1361 : ! .and. (conv_vel_interp(4*k-3)<1d-5 .or. conv_vel_interp(4*(k+1)-3)<1d-5)) then
1362 : ! s% conv_vel(k0) = 0d0
1363 : ! exit
1364 : ! end if
1365 : ! end do
1366 : ! s% xh(s% i_ln_cvpv0, k0) = log(s% conv_vel(k0)+s% conv_vel_v0)
1367 : ! end do
1368 : ! write(*,*) 'need to rewrite some things here in do_relax_to_star_cut'
1369 : ! call mesa_error(__FILE__,__LINE__,'do_relax_to_star_cut')
1370 : ! deallocate(conv_vel_interp, interp_work)
1371 : !end if
1372 :
1373 0 : s% generations = 1
1374 :
1375 : ! restore v_flag and u_flag
1376 0 : if (v_flag) then
1377 0 : call set_v_flag(id, .true., ierr)
1378 : end if
1379 0 : if (u_flag) then
1380 0 : call set_u_flag(id, .true., ierr)
1381 : end if
1382 :
1383 : ! this avoids the history file from being rewritten
1384 0 : s% doing_first_model_of_run = .false.
1385 :
1386 0 : s% time = time
1387 0 : s% model_number = model_number
1388 0 : s% num_trace_history_values = num_trace_history_values
1389 0 : s% cumulative_energy_error = cumulative_energy_error
1390 0 : s% cumulative_extra_heating = cumulative_extra_heating
1391 :
1392 0 : s% non_nuc_neu_factor = non_nuc_neu_factor
1393 0 : s% eps_nuc_factor = eps_nuc_factor
1394 :
1395 0 : s% initial_z = initial_z
1396 0 : s% initial_y = initial_y
1397 0 : s% initial_mass = initial_mass
1398 0 : s% photo_interval = photo_interval
1399 :
1400 0 : deallocate(q, xq)
1401 :
1402 0 : s% need_to_setvars = .true.
1403 :
1404 0 : end subroutine do_relax_to_star_cut
1405 :
1406 : end module remove_shells
|
