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 relax
21 :
22 : use star_private_def
23 : use const_def, only: dp, msun, rsun, secyer, one_third, four_thirds_pi
24 : use utils_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_relax_z
30 : public :: do_relax_mass
31 : public :: do_relax_to_limit
32 : public :: do_relax_mass_scale
33 : public :: do_relax_num_steps
34 : public :: do_relax_to_radiative_core
35 : public :: do_relax_composition
36 : public :: do_relax_angular_momentum
37 : public :: do_relax_entropy
38 : public :: do_relax_core
39 : public :: do_relax_m_center
40 : public :: do_relax_r_center
41 : public :: do_relax_v_center
42 : public :: do_relax_l_center
43 : public :: do_relax_dxdt_nuc_factor
44 : public :: do_relax_eps_nuc_factor
45 : public :: do_relax_opacity_max
46 : public :: do_relax_max_surf_dq
47 : public :: do_relax_to_xaccrete
48 : public :: do_relax_y
49 : public :: do_relax_tau_factor
50 : public :: do_relax_opacity_factor
51 : public :: do_relax_tsurf_factor
52 : public :: do_relax_uniform_omega
53 : public :: do_relax_irradiation
54 : public :: do_relax_mass_change
55 : public :: do_internal_evolve
56 :
57 : real(dp), parameter :: min_dlnz = -12
58 : real(dp), parameter :: min_z = 1d-12
59 :
60 : ! some relax routines depend on things such as other_energy and other_torque
61 : ! to which interpolation parameters cannot be passed directly. So for simplicity
62 : ! use these two global variables instead.
63 : integer :: relax_num_pts
64 : real(dp), pointer :: relax_work_array(:)
65 :
66 : contains
67 :
68 0 : subroutine do_relax_to_limit(id, restore_at_end, ierr)
69 : integer, intent(in) :: id
70 : logical, intent(in) :: restore_at_end
71 : integer, intent(out) :: ierr
72 :
73 : integer, parameter :: lipar=1, lrpar=1
74 :
75 : type (star_info), pointer :: s
76 :
77 : integer, target :: ipar_ary(lipar)
78 : integer, pointer :: ipar(:)
79 0 : real(dp), target :: rpar_ary(lrpar)
80 : real(dp), pointer :: rpar(:)
81 0 : rpar => rpar_ary
82 0 : ipar => ipar_ary
83 :
84 : ierr = 0
85 0 : call get_star_ptr(id, s, ierr)
86 0 : if (ierr /= 0) return
87 : call do_internal_evolve( &
88 : id, before_relax_to_limit, relax_to_limit_adjust_model, &
89 : relax_to_limit_check_model, null_finish_model, &
90 0 : restore_at_end, lipar, ipar, lrpar, rpar, ierr)
91 : end subroutine do_relax_to_limit
92 :
93 :
94 0 : subroutine before_relax_to_limit(s, id, lipar, ipar, lrpar, rpar, ierr)
95 : type (star_info), pointer :: s
96 : integer, intent(in) :: id, lipar, lrpar
97 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
98 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
99 : integer, intent(out) :: ierr
100 0 : ierr = 0
101 0 : end subroutine before_relax_to_limit
102 :
103 0 : integer function relax_to_limit_adjust_model(s, id, lipar, ipar, lrpar, rpar)
104 : type (star_info), pointer :: s
105 : integer, intent(in) :: id, lipar, lrpar
106 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
107 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
108 0 : relax_to_limit_adjust_model = keep_going
109 0 : end function relax_to_limit_adjust_model
110 :
111 0 : integer function relax_to_limit_check_model(s, id, lipar, ipar, lrpar, rpar)
112 : use do_one_utils, only: do_bare_bones_check_model, do_check_limits
113 : type (star_info), pointer :: s
114 : integer, intent(in) :: id, lipar, lrpar
115 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
116 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
117 0 : relax_to_limit_check_model = do_bare_bones_check_model(id)
118 0 : if (relax_to_limit_check_model == keep_going) then
119 0 : relax_to_limit_check_model = do_check_limits(id)
120 0 : if (relax_to_limit_check_model == terminate) &
121 0 : s% termination_code = t_relax_finished_okay
122 : end if
123 0 : end function relax_to_limit_check_model
124 :
125 :
126 0 : subroutine do_relax_mass(id, new_mass, lg_max_abs_mdot, ierr)
127 : integer, intent(in) :: id
128 : real(dp), intent(in) :: new_mass, lg_max_abs_mdot
129 : integer, intent(out) :: ierr
130 : integer, parameter :: lipar=1, lrpar=3
131 : integer :: max_model_number
132 : character (len=32) :: cool_wind_AGB_scheme, cool_wind_RGB_scheme
133 0 : real(dp) :: starting_dt_next, varcontrol_target, eps_mdot_factor
134 : logical :: adding_mass
135 : type (star_info), pointer :: s
136 : integer, target :: ipar_ary(lipar)
137 : integer, pointer :: ipar(:)
138 0 : real(dp), target :: rpar_ary(lrpar)
139 : real(dp), pointer :: rpar(:)
140 0 : rpar => rpar_ary
141 0 : ipar => ipar_ary
142 : include 'formats'
143 : ierr = 0
144 0 : call get_star_ptr(id, s, ierr)
145 0 : if (ierr /= 0) return
146 0 : if (abs(new_mass - s% star_mass) < 1d-12*new_mass) then
147 0 : s% star_mass = new_mass
148 0 : s% mstar = new_mass*Msun
149 0 : s% xmstar = s% mstar - s% M_center
150 0 : return
151 : end if
152 0 : write(*,'(A)')
153 0 : write(*,1) 'current mass', s% mstar/Msun
154 0 : write(*,1) 'relax to new_mass', new_mass
155 0 : write(*,1) 'lg_max_abs_mdot', lg_max_abs_mdot
156 0 : write(*,'(A)')
157 0 : if (new_mass <= 0) then
158 0 : ierr = -1
159 0 : write(*,*) 'invalid new mass'
160 0 : return
161 : end if
162 :
163 0 : rpar(1) = new_mass*Msun
164 0 : rpar(2) = lg_max_abs_mdot
165 0 : rpar(3) = s% mstar
166 :
167 0 : adding_mass = (new_mass > s% star_mass)
168 :
169 0 : cool_wind_AGB_scheme = s% cool_wind_AGB_scheme
170 0 : s% cool_wind_AGB_scheme = ''
171 :
172 0 : cool_wind_RGB_scheme = s% cool_wind_RGB_scheme
173 0 : s% cool_wind_RGB_scheme = ''
174 :
175 0 : varcontrol_target = s% varcontrol_target
176 0 : s% varcontrol_target = 1d-3
177 :
178 0 : eps_mdot_factor = s% eps_mdot_factor
179 0 : s% eps_mdot_factor = 0
180 :
181 0 : max_model_number = s% max_model_number
182 0 : s% max_model_number = -1111
183 0 : starting_dt_next = s% dt_next
184 : call do_internal_evolve( &
185 : id, before_evolve_relax_mass, relax_mass_adjust_model, relax_mass_check_model, &
186 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
187 0 : s% max_model_number = max_model_number
188 0 : s% dt_next = min(s% dt_next, starting_dt_next) * 1d-1
189 0 : s% initial_mass = new_mass
190 :
191 0 : s% cool_wind_AGB_scheme = cool_wind_AGB_scheme
192 0 : s% cool_wind_RGB_scheme = cool_wind_RGB_scheme
193 0 : s% varcontrol_target = varcontrol_target
194 0 : s% eps_mdot_factor = eps_mdot_factor
195 0 : s% star_mass = new_mass
196 0 : s% mstar = new_mass*Msun
197 0 : s% xmstar = s% mstar - s% M_center
198 :
199 0 : call error_check('relax mass',ierr)
200 :
201 0 : write(*,2) 's% max_model_number', s% max_model_number
202 :
203 :
204 0 : end subroutine do_relax_mass
205 :
206 :
207 0 : subroutine before_evolve_relax_mass(s, id, lipar, ipar, lrpar, rpar, ierr)
208 : type (star_info), pointer :: s
209 : integer, intent(in) :: id, lipar, lrpar
210 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
211 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
212 : integer, intent(out) :: ierr
213 0 : ierr = 0
214 0 : call setup_before_relax(s)
215 0 : s% mass_change = 0
216 0 : s% dt_next = min(s% dt_next, 1d4*secyer)
217 0 : end subroutine before_evolve_relax_mass
218 :
219 0 : integer function relax_mass_adjust_model(s, id, lipar, ipar, lrpar, rpar)
220 : type (star_info), pointer :: s
221 : integer, intent(in) :: id, lipar, lrpar
222 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
223 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
224 0 : relax_mass_adjust_model = keep_going
225 0 : end function relax_mass_adjust_model
226 :
227 0 : integer function relax_mass_check_model(s, id, lipar, ipar, lrpar, rpar)
228 : use do_one_utils, only:do_bare_bones_check_model
229 : type (star_info), pointer :: s
230 : integer, intent(in) :: id, lipar, lrpar
231 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
232 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
233 : integer :: ramp
234 : real(dp) :: &
235 0 : new_mass, old_mass, mdot, max_abs_mdot, abs_diff, lg_max_abs_mdot
236 :
237 : logical, parameter :: dbg = .false.
238 :
239 : include 'formats'
240 :
241 : if (dbg) write(*,*) 'relax_mass_check_model'
242 :
243 0 : relax_mass_check_model = do_bare_bones_check_model(id)
244 0 : if (relax_mass_check_model /= keep_going) then
245 0 : write(*,*) 'forced termination'
246 0 : return
247 : end if
248 :
249 0 : new_mass = rpar(1)
250 0 : lg_max_abs_mdot = rpar(2)
251 0 : if (lg_max_abs_mdot <= -100) then ! use default
252 0 : if (s% star_mass < 0.003d0) then
253 0 : lg_max_abs_mdot = -7d0
254 0 : else if (s% star_mass < 0.006d0) then
255 0 : lg_max_abs_mdot = -6.3d0
256 0 : else if (s% star_mass < 0.01d0) then
257 0 : lg_max_abs_mdot = -6d0
258 0 : else if (s% star_mass < 0.1d0) then
259 0 : lg_max_abs_mdot = -4d0
260 : else
261 0 : lg_max_abs_mdot = -3d0
262 : end if
263 : end if
264 0 : max_abs_mdot = exp10(lg_max_abs_mdot)*(Msun/secyer)
265 0 : old_mass = rpar(3)
266 :
267 0 : if (s% model_number >= s% max_model_number .and. s% max_model_number > 0) then
268 0 : write(*,3) 's% model_number >= s% max_model_number', &
269 0 : s% model_number, s% max_model_number
270 0 : relax_mass_check_model = terminate
271 0 : return
272 : end if
273 :
274 0 : abs_diff = abs(s% mstar - new_mass)
275 0 : mdot = (new_mass - s% mstar)/s% dt_next
276 0 : if (abs(mdot) > max_abs_mdot) then
277 0 : mdot = sign(max_abs_mdot, mdot)
278 : end if
279 :
280 0 : if (abs_diff < 1d-4*new_mass .or. abs(mdot) < 1d-50) then
281 0 : s% mass_change = 0
282 0 : s% star_mass = new_mass
283 0 : s% mstar = new_mass*Msun
284 0 : s% xmstar = s% mstar - s% M_center
285 : s% mass_change = 0
286 0 : write(*,1) 's% tau_base =', s% tau_base
287 0 : write(*,1) 's% tau_factor =', s% tau_factor
288 0 : write(*,1) 'atm_option = ' // trim(s% atm_option)
289 0 : relax_mass_check_model = terminate
290 0 : s% termination_code = t_relax_finished_okay
291 0 : return
292 : end if
293 :
294 0 : s% max_timestep = abs(new_mass - s% mstar)/mdot
295 :
296 : if (dbg) then
297 : write(*,1) 'new_mass/Msun', new_mass/Msun
298 : write(*,1) 'old_mass/Msun', old_mass/Msun
299 : write(*,1) 'abs_diff/Msun', abs_diff/Msun
300 : end if
301 :
302 0 : ramp = 12
303 0 : if (s% model_number < ramp) then
304 0 : mdot = mdot * pow(1.1d0,dble(s% model_number-ramp))
305 : end if
306 :
307 0 : if (abs(mdot)*s% dt > 0.05d0*s% mstar) mdot = sign(0.05d0*s% mstar/s% dt,mdot)
308 :
309 0 : s% mass_change = mdot/(Msun/secyer)
310 :
311 : if (dbg) write(*,1) 's% mass_change', s% mass_change
312 :
313 0 : end function relax_mass_check_model
314 :
315 :
316 0 : subroutine do_relax_composition( &
317 0 : id, num_steps_to_use, num_pts, species, xa, xq, ierr)
318 : integer, intent(in) :: id
319 : integer, intent(in) :: num_steps_to_use ! use this many steps to do conversion
320 : integer, intent(in) :: num_pts
321 : ! length of composition vector; need not equal nz for current model
322 : integer, intent(in) :: species
323 : ! per point; must = number of species for current model
324 : real(dp), intent(in) :: xa(species,num_pts) ! desired composition profile
325 : real(dp), intent(in) :: xq(num_pts)
326 : integer, intent(out) :: ierr
327 :
328 : integer, parameter :: lipar=3
329 : integer :: lrpar, max_model_number
330 0 : real(dp), pointer :: rpar(:)
331 0 : real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor
332 : logical :: do_element_diffusion, include_composition_in_eps_grav
333 : type (star_info), pointer :: s
334 0 : real(dp), pointer :: x(:), f1(:), f(:,:,:)
335 : integer, target :: ipar_ary(lipar)
336 : integer, pointer :: ipar(:)
337 0 : ipar => ipar_ary
338 :
339 : include 'formats'
340 :
341 : ierr = 0
342 :
343 0 : call get_star_ptr(id, s, ierr)
344 0 : if (ierr /= 0) return
345 :
346 0 : if (species /= s% species) then
347 0 : ierr = -1
348 0 : return
349 : end if
350 :
351 0 : ipar(1) = num_pts
352 0 : ipar(2) = num_steps_to_use
353 0 : ipar(3) = s% model_number
354 0 : lrpar = (1 + 4*species)*num_pts
355 0 : allocate(rpar(lrpar), stat=ierr)
356 0 : if (ierr /= 0) return
357 :
358 0 : x(1:num_pts) => rpar(1:num_pts)
359 0 : f1(1:4*num_pts*species) => rpar(num_pts+1:lrpar)
360 0 : f(1:4,1:num_pts,1:species) => f1(1:4*num_pts*species)
361 :
362 0 : call store_rpar(species, num_pts, ierr)
363 0 : if (ierr /= 0) return
364 :
365 0 : max_model_number = s% max_model_number
366 0 : s% max_model_number = num_steps_to_use + s% model_number + 1
367 0 : write(*,*) 'relax_composition: num_steps_to_use', num_steps_to_use
368 :
369 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
370 0 : s% dxdt_nuc_factor = 0 ! turn off composition change by nuclear burning
371 0 : mix_factor = s% mix_factor
372 0 : s% mix_factor = 0 ! turn off mixing
373 0 : do_element_diffusion = s% do_element_diffusion
374 0 : s% do_element_diffusion = .false. ! turn off diffusion
375 0 : include_composition_in_eps_grav = s% include_composition_in_eps_grav
376 0 : s% include_composition_in_eps_grav = .false. ! don't need energetic effects of artificial changes
377 0 : starting_dt_next = s% dt_next
378 : call do_internal_evolve( &
379 : id, before_evolve_relax_composition, &
380 : relax_composition_adjust_model, relax_composition_check_model, &
381 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
382 :
383 0 : s% max_model_number = max_model_number
384 0 : s% dt_next = starting_dt_next
385 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
386 0 : s% mix_factor = mix_factor
387 0 : s% do_element_diffusion = do_element_diffusion
388 0 : s% include_composition_in_eps_grav = include_composition_in_eps_grav
389 :
390 0 : call error_check('relax composition',ierr)
391 :
392 0 : deallocate(rpar)
393 :
394 :
395 : contains
396 :
397 :
398 0 : subroutine store_rpar(species, num_pts, ierr) ! get interpolation info
399 : use interp_1d_def, only: pm_work_size
400 : use interp_1d_lib, only: interp_pm
401 : integer, intent(in) :: species, num_pts
402 : integer, intent(out) :: ierr
403 : integer :: j, op_err
404 0 : real(dp), pointer :: interp_work(:), work(:), p(:)
405 0 : allocate(interp_work(num_pts*pm_work_size*species), stat=ierr)
406 0 : if (ierr /= 0) return
407 0 : x(:) = xq(:)
408 0 : do j=1, species ! make piecewise monotonic cubic interpolants
409 : op_err = 0
410 0 : f(1,1:num_pts,j) = xa(j,1:num_pts)
411 : work(1:num_pts*pm_work_size) => &
412 0 : interp_work(1+num_pts*pm_work_size*(j-1):num_pts*pm_work_size*j)
413 0 : p(1:4*num_pts) => f1(1+4*num_pts*(j-1):4*num_pts*j)
414 : call interp_pm(x, num_pts, p, pm_work_size, work, &
415 0 : 'do_relax_composition', op_err)
416 0 : if (op_err /= 0) ierr = op_err
417 : end do
418 0 : deallocate(interp_work)
419 0 : end subroutine store_rpar
420 :
421 : end subroutine do_relax_composition
422 :
423 :
424 0 : subroutine before_evolve_relax_composition(s, id, lipar, ipar, lrpar, rpar, ierr)
425 : type (star_info), pointer :: s
426 : integer, intent(in) :: id, lipar, lrpar
427 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
428 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
429 : integer, intent(out) :: ierr
430 0 : ierr = 0
431 0 : call setup_before_relax(s)
432 0 : end subroutine before_evolve_relax_composition
433 :
434 0 : integer function relax_composition_adjust_model(s, id, lipar, ipar, lrpar, rpar)
435 : type (star_info), pointer :: s
436 : integer, intent(in) :: id, lipar, lrpar
437 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
438 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
439 : integer :: num_pts, num_steps_to_use, species, starting_model_number, ierr
440 0 : real(dp) :: lambda, avg_err
441 :
442 0 : real(dp), pointer :: x(:) ! =(num_pts)
443 0 : real(dp), pointer :: f1(:) ! =(4, num_pts, species)
444 :
445 : include 'formats'
446 :
447 0 : num_pts = ipar(1)
448 0 : num_steps_to_use = ipar(2)
449 0 : starting_model_number = ipar(3)
450 0 : species = s% species
451 :
452 0 : if (lrpar /= (1 + 4*species)*num_pts) then
453 0 : write(*,*) 'bad lrpar for relax_composition_check_model'
454 0 : relax_composition_adjust_model = terminate
455 0 : return
456 : end if
457 :
458 0 : ierr = 0
459 0 : if (s% job% timescale_for_relax_composition < 0d0) then
460 : lambda = min(1d0, max(0d0, (s% model_number - starting_model_number - 1) &
461 0 : / max(1d0, dble(num_steps_to_use) - 1)))
462 : else
463 0 : lambda = min(1d0, s% dt/s% job% timescale_for_relax_composition)
464 : end if
465 :
466 0 : x(1:num_pts) => rpar(1:num_pts)
467 0 : f1(1:4*num_pts*species) => rpar(num_pts+1:lrpar)
468 0 : call adjust_xa(species, num_pts, avg_err, ierr)
469 0 : if (ierr /= 0) relax_composition_adjust_model = terminate
470 :
471 0 : write(*,1) 'avg remaining difference, lambda', avg_err, lambda
472 :
473 0 : relax_composition_adjust_model = keep_going
474 :
475 : contains
476 :
477 :
478 0 : subroutine adjust_xa(species, num_pts, avg_err, ierr)
479 : use interp_1d_lib, only: interp_values
480 : integer, intent(in) :: species, num_pts
481 : real(dp), intent(out) :: avg_err
482 : integer, intent(out) :: ierr
483 : integer :: j, k, nz, op_err
484 0 : real(dp), pointer :: vals(:,:), xq(:), f(:)
485 0 : ierr = 0
486 0 : nz = s% nz
487 0 : allocate(vals(species, nz), xq(nz), stat=ierr)
488 0 : if (ierr /= 0) return
489 0 : xq(1) = 0.5d0*s% dq(1) ! xq for cell center
490 0 : do k = 2, nz
491 0 : xq(k) = xq(k-1) + 0.5d0*(s% dq(k) + s% dq(k-1))
492 : end do
493 0 : !$OMP PARALLEL DO PRIVATE(j,op_err,f) SCHEDULE(dynamic,2)
494 : do j=1, species ! interpolate target composition
495 : f(1:4*num_pts) => f1(1+(j-1)*4*num_pts:j*4*num_pts)
496 : call interp_values(x, num_pts, f, nz, xq, vals(j,:), op_err)
497 : ! enforce non-negative mass fractions
498 : ! if the abundance switches back and forth between 0 and 1d-99,
499 : ! then small negative abundances ~ -1d-115 can be generated
500 : do k = 1, nz
501 : if (vals(j,k) < 0d0) vals(j,k) = 0d0
502 : end do
503 : if (op_err /= 0) ierr = op_err
504 : s% xa(j,1:nz) = (1d0-lambda)*s% xa(j,1:nz) + lambda*vals(j,1:nz)
505 : end do
506 : !$OMP END PARALLEL DO
507 0 : avg_err = sum(abs(s% xa(1:species,1:nz)-vals(1:species,1:nz)))/(nz*species)
508 0 : deallocate(vals, xq)
509 0 : end subroutine adjust_xa
510 :
511 : end function relax_composition_adjust_model
512 :
513 0 : integer function relax_composition_check_model(s, id, lipar, ipar, lrpar, rpar)
514 : use do_one_utils, only:do_bare_bones_check_model
515 : type (star_info), pointer :: s
516 : integer, intent(in) :: id, lipar, lrpar
517 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
518 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
519 :
520 : integer :: num_steps_to_use, starting_model_number
521 :
522 : include 'formats'
523 :
524 0 : relax_composition_check_model = do_bare_bones_check_model(id)
525 0 : if (relax_composition_check_model /= keep_going) return
526 :
527 0 : num_steps_to_use = ipar(2)
528 0 : starting_model_number = ipar(3)
529 :
530 0 : if (s% job% timescale_for_relax_composition < 0d0) then
531 0 : if (s% model_number - starting_model_number >= num_steps_to_use) then
532 0 : relax_composition_check_model = terminate
533 0 : s% termination_code = t_relax_finished_okay
534 0 : return
535 : end if
536 : else
537 : if (s% generations > 1 .and. &
538 0 : s% dt > s% job% timescale_for_relax_composition .and. &
539 : s% dt_old > s% job% timescale_for_relax_composition) then
540 0 : relax_composition_check_model = terminate
541 0 : s% termination_code = t_relax_finished_okay
542 0 : return
543 : end if
544 : end if
545 :
546 :
547 0 : end function relax_composition_check_model
548 :
549 :
550 0 : subroutine do_relax_to_xaccrete(id, num_steps_to_use, ierr)
551 0 : use adjust_xyz, only: get_xa_for_accretion
552 :
553 : integer, intent(in) :: id
554 : integer, intent(in) :: num_steps_to_use ! use this many steps to do conversion
555 : integer, intent(out) :: ierr
556 :
557 : integer, parameter :: lipar=2
558 : integer :: lrpar, max_model_number, species
559 0 : real(dp), pointer :: rpar(:)
560 0 : real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor
561 : logical :: do_element_diffusion
562 : type (star_info), pointer :: s
563 0 : real(dp), pointer :: xa(:)
564 : integer, target :: ipar_ary(lipar)
565 : integer, pointer :: ipar(:)
566 0 : ipar => ipar_ary
567 :
568 : include 'formats'
569 :
570 : ierr = 0
571 :
572 0 : call get_star_ptr(id, s, ierr)
573 0 : if (ierr /= 0) return
574 :
575 0 : species = s% species
576 0 : ipar(1) = s% model_number
577 0 : ipar(2) = num_steps_to_use
578 0 : if (num_steps_to_use <= 0) then
579 0 : ierr = -1
580 0 : write(*,2) 'invalid num_steps_to_use to relax_to_xaccrete', num_steps_to_use
581 0 : return
582 : end if
583 :
584 0 : lrpar = species
585 0 : allocate(rpar(lrpar), stat=ierr)
586 0 : if (ierr /= 0) return
587 :
588 0 : xa(1:species) => rpar(1:species)
589 :
590 0 : call get_xa_for_accretion(s, xa, ierr)
591 0 : if (ierr /= 0) then
592 0 : if (s% report_ierr) &
593 0 : write(*, *) 'get_xa_for_accretion failed in relax_to_xaccrete'
594 0 : deallocate(rpar)
595 0 : return
596 : end if
597 :
598 0 : max_model_number = s% max_model_number
599 0 : s% max_model_number = num_steps_to_use + 1
600 0 : write(*,*) 'num_steps_to_use', num_steps_to_use
601 :
602 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
603 0 : s% dxdt_nuc_factor = 0 ! turn off composition change by nuclear burning
604 0 : mix_factor = s% mix_factor
605 0 : s% mix_factor = 0 ! turn off mixing
606 0 : do_element_diffusion = s% do_element_diffusion
607 0 : do_element_diffusion = .false. ! turn off diffusion
608 0 : starting_dt_next = s% dt_next
609 :
610 : call do_internal_evolve( &
611 : id, before_evolve_relax_to_xaccrete, &
612 : relax_to_xaccrete_adjust_model, relax_to_xaccrete_check_model, &
613 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
614 :
615 0 : s% max_model_number = max_model_number
616 0 : s% dt_next = starting_dt_next
617 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
618 0 : s% mix_factor = mix_factor
619 0 : s% do_element_diffusion = do_element_diffusion
620 :
621 0 : call error_check('relax to xacrrete',ierr)
622 :
623 0 : deallocate(rpar)
624 :
625 0 : end subroutine do_relax_to_xaccrete
626 :
627 :
628 0 : subroutine before_evolve_relax_to_xaccrete(s, id, lipar, ipar, lrpar, rpar, ierr)
629 : type (star_info), pointer :: s
630 : integer, intent(in) :: id, lipar, lrpar
631 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
632 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
633 : integer, intent(out) :: ierr
634 0 : ierr = 0
635 0 : call setup_before_relax(s)
636 0 : end subroutine before_evolve_relax_to_xaccrete
637 :
638 0 : integer function relax_to_xaccrete_adjust_model(s, id, lipar, ipar, lrpar, rpar)
639 : type (star_info), pointer :: s
640 : integer, intent(in) :: id, lipar, lrpar
641 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
642 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
643 0 : relax_to_xaccrete_adjust_model = keep_going
644 0 : end function relax_to_xaccrete_adjust_model
645 :
646 0 : integer function relax_to_xaccrete_check_model(s, id, lipar, ipar, lrpar, rpar)
647 : use do_one_utils, only:do_bare_bones_check_model
648 : type (star_info), pointer :: s
649 : integer, intent(in) :: id, lipar, lrpar
650 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
651 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
652 :
653 : integer :: num_steps_to_use, starting_model_number, species, k, j
654 0 : real(dp), pointer :: xa(:)
655 : real(dp) :: frac
656 :
657 : include 'formats'
658 :
659 0 : relax_to_xaccrete_check_model = do_bare_bones_check_model(id)
660 0 : if (relax_to_xaccrete_check_model /= keep_going) return
661 :
662 0 : starting_model_number = ipar(1)
663 0 : num_steps_to_use = ipar(2)
664 0 : species = s% species
665 :
666 0 : frac = dble(s% model_number - starting_model_number)/dble(num_steps_to_use)
667 0 : frac = frac*frac
668 :
669 0 : if (lrpar /= species) then
670 0 : write(*,*) 'bad lrpar for relax_to_xaccrete_check_model'
671 0 : relax_to_xaccrete_check_model = terminate
672 : end if
673 :
674 0 : xa(1:species) => rpar(1:species)
675 :
676 0 : do k=1,s% nz
677 0 : do j=1,species
678 0 : s% xa(j,k) = (1d0-frac)*s% xa(j,k) + frac*xa(j)
679 : end do
680 : end do
681 :
682 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
683 0 : write(*,2) 'relax to xaccrete: fraction', s% model_number, frac
684 :
685 0 : if (s% model_number - starting_model_number >= num_steps_to_use) then
686 0 : relax_to_xaccrete_check_model = terminate
687 0 : s% termination_code = t_relax_finished_okay
688 0 : return
689 : end if
690 :
691 :
692 0 : end function relax_to_xaccrete_check_model
693 :
694 0 : subroutine do_relax_entropy( &
695 0 : id, max_steps_to_use, num_pts, entropy, xq, ierr)
696 0 : use alloc, only: alloc_extras, dealloc_extras
697 : integer, intent(in) :: id
698 : integer, intent(in) :: max_steps_to_use ! use this many steps to do conversion
699 : integer, intent(in) :: num_pts
700 : ! length of entropy vector; need not equal nz for current model
701 : real(dp), intent(in) :: entropy(num_pts) ! desired entropy profile
702 : real(dp), intent(in) :: xq(num_pts)
703 : integer, intent(out) :: ierr
704 :
705 : integer, parameter :: lipar=2
706 : integer :: lrpar, max_model_number
707 0 : real(dp), pointer :: rpar(:)
708 0 : real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor, max_years_for_timestep, time
709 : logical :: do_element_diffusion, use_other_energy
710 : type (star_info), pointer :: s
711 0 : real(dp), pointer :: x(:), f1(:), f(:,:)
712 : integer, target :: ipar_ary(lipar)
713 : integer, pointer :: ipar(:)
714 : procedure (other_energy_interface), pointer :: &
715 : other_energy => null()
716 :
717 0 : ipar => ipar_ary
718 :
719 : include 'formats'
720 :
721 : ierr = 0
722 :
723 0 : call get_star_ptr(id, s, ierr)
724 0 : if (ierr /= 0) return
725 :
726 0 : max_model_number = s% max_model_number
727 0 : s% max_model_number = max_steps_to_use + s% model_number + 1
728 0 : write(*,*) 'relax_entropy: max_steps_to_use', max_steps_to_use
729 :
730 0 : ipar(1) = num_pts
731 0 : ipar(2) = s% max_model_number
732 0 : lrpar = 5*num_pts
733 0 : allocate(rpar(lrpar), stat=ierr)
734 0 : if (ierr /= 0) return
735 :
736 0 : x(1:num_pts) => rpar(1:num_pts)
737 0 : f1(1:4*num_pts) => rpar(num_pts+1:lrpar)
738 0 : f(1:4,1:num_pts) => f1(1:4*num_pts)
739 :
740 0 : call store_rpar(num_pts, ierr)
741 0 : if (ierr /= 0) return
742 :
743 : ! need to use global variables, as relax_entropy uses
744 : ! the other_energy routine to which it can't pass rpar
745 0 : relax_num_pts = num_pts
746 0 : relax_work_array => rpar
747 :
748 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
749 0 : s% dxdt_nuc_factor = 0 ! turn off composition change by nuclear burning
750 0 : mix_factor = s% mix_factor
751 0 : s% mix_factor = 0 ! turn off mixing
752 0 : do_element_diffusion = s% do_element_diffusion
753 0 : s% do_element_diffusion = .false. ! turn off diffusion
754 0 : starting_dt_next = s% dt_next
755 0 : max_years_for_timestep = s% max_years_for_timestep
756 0 : s% max_years_for_timestep = s% job% max_dt_for_relax_entropy
757 0 : s% dt_next = min(s% dt_next, s% job% max_dt_for_relax_entropy * secyer)
758 0 : use_other_energy = s% use_other_energy
759 0 : s% use_other_energy = .true.
760 0 : other_energy => s% other_energy
761 0 : s% other_energy => entropy_relax_other_energy
762 0 : time = s% time
763 0 : s% time = 0d0
764 :
765 : call do_internal_evolve( &
766 : id, before_evolve_relax_entropy, &
767 : relax_entropy_adjust_model, relax_entropy_check_model, &
768 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
769 :
770 0 : s% max_model_number = max_model_number
771 0 : s% dt_next = starting_dt_next
772 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
773 0 : s% mix_factor = mix_factor
774 0 : s% do_element_diffusion = do_element_diffusion
775 0 : s% max_years_for_timestep = max_years_for_timestep
776 0 : s% use_other_energy = use_other_energy
777 0 : s% other_energy => other_energy
778 0 : s% time = time
779 :
780 0 : call error_check('relax entropy',ierr)
781 :
782 0 : deallocate(rpar)
783 :
784 :
785 : contains
786 :
787 :
788 0 : subroutine store_rpar(num_pts, ierr) ! get interpolation info
789 0 : use interp_1d_def, only: pm_work_size
790 : use interp_1d_lib, only: interp_pm
791 : integer, intent(in) :: num_pts
792 : integer, intent(out) :: ierr
793 : integer :: op_err
794 0 : real(dp), pointer :: interp_work(:), work(:), p(:)
795 0 : allocate(interp_work(num_pts*pm_work_size), stat=ierr)
796 0 : if (ierr /= 0) return
797 0 : x(:) = xq(:)
798 : op_err = 0
799 0 : f(1,1:num_pts) = entropy(1:num_pts)
800 : work(1:num_pts*pm_work_size) => &
801 0 : interp_work(1:num_pts*pm_work_size)
802 0 : p(1:4*num_pts) => f1(1:4*num_pts)
803 : call interp_pm(x, num_pts, p, pm_work_size, work, &
804 0 : 'do_relax_entropy', op_err)
805 0 : if (op_err /= 0) ierr = op_err
806 0 : deallocate(interp_work)
807 0 : end subroutine store_rpar
808 :
809 : end subroutine do_relax_entropy
810 :
811 0 : subroutine before_evolve_relax_entropy(s, id, lipar, ipar, lrpar, rpar, ierr)
812 : type (star_info), pointer :: s
813 : integer, intent(in) :: id, lipar, lrpar
814 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
815 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
816 : integer, intent(out) :: ierr
817 0 : ierr = 0
818 0 : call setup_before_relax(s)
819 0 : end subroutine before_evolve_relax_entropy
820 :
821 0 : integer function relax_entropy_adjust_model(s, id, lipar, ipar, lrpar, rpar)
822 : type (star_info), pointer :: s
823 : integer, intent(in) :: id, lipar, lrpar
824 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
825 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
826 0 : relax_entropy_adjust_model = keep_going
827 0 : end function relax_entropy_adjust_model
828 :
829 0 : integer function relax_entropy_check_model(s, id, lipar, ipar, lrpar, rpar)
830 : use do_one_utils, only: do_bare_bones_check_model
831 : type (star_info), pointer :: s
832 : integer, intent(in) :: id, lipar, lrpar
833 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
834 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
835 :
836 : integer :: num_pts, ierr, max_model_number
837 0 : real(dp) :: avg_err
838 0 : real(dp), pointer :: x(:) ! =(num_pts)
839 0 : real(dp), pointer :: f1(:) ! =(4, num_pts)
840 :
841 : include 'formats'
842 :
843 0 : relax_entropy_check_model = do_bare_bones_check_model(id)
844 0 : if (relax_entropy_check_model /= keep_going) return
845 :
846 0 : if (lipar /= 2) then
847 0 : write(*,*) 'bad lipar for relax_entropy_check_model'
848 0 : relax_entropy_check_model = terminate
849 0 : return
850 : end if
851 :
852 0 : num_pts = ipar(1)
853 0 : max_model_number = ipar(2)
854 :
855 0 : if (lrpar /= 5*num_pts) then
856 0 : write(*,*) 'bad lrpar for relax_entropy_check_model'
857 0 : relax_entropy_check_model = terminate
858 : end if
859 :
860 : ierr = 0
861 0 : x(1:num_pts) => rpar(1:num_pts)
862 0 : f1(1:4*num_pts) => rpar(num_pts+1:lrpar)
863 0 : call adjust_entropy(num_pts, avg_err, ierr)
864 0 : if (ierr /= 0) relax_entropy_check_model = terminate
865 :
866 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
867 0 : write(*,*) 'relax_entropy avg rel err, dt, model', avg_err, s% dt/secyer, s% model_number
868 :
869 :
870 0 : if (s% star_age >= s% job% timescale_for_relax_entropy*s% job% num_timescales_for_relax_entropy) then
871 0 : relax_entropy_check_model = terminate
872 0 : s% termination_code = t_relax_finished_okay
873 0 : return
874 : end if
875 :
876 0 : if (s% model_number >= max_model_number) then
877 0 : write(*,*) "Terminated relax because of max_model_number instead of relax_time"
878 0 : relax_entropy_check_model = terminate
879 0 : s% termination_code = t_relax_finished_okay
880 0 : return
881 : end if
882 :
883 :
884 : contains
885 :
886 :
887 0 : subroutine adjust_entropy(num_pts, avg_err, ierr)
888 0 : use interp_1d_lib, only: interp_values
889 : integer, intent(in) :: num_pts
890 : real(dp), intent(out) :: avg_err
891 : integer, intent(out) :: ierr
892 : integer :: k, nz, op_err
893 0 : real(dp) :: dentropy_sum
894 0 : real(dp), pointer :: vals(:), xq(:), f(:), entropy(:)
895 0 : ierr = 0
896 0 : nz = s% nz
897 0 : allocate(vals(nz), xq(nz), stat=ierr)
898 0 : if (ierr /= 0) return
899 0 : xq(1) = s% dq(1)/2 ! xq for cell center
900 0 : do k = 2, nz
901 0 : xq(k) = xq(k-1) + (s% dq(k) + s% dq(k-1))/2
902 : end do
903 0 : dentropy_sum = 0
904 0 : f(1:4*num_pts) => f1(1:4*num_pts)
905 0 : call interp_values(x, num_pts, f, nz, xq, vals(:), op_err)
906 0 : if (op_err /= 0) ierr = op_err
907 0 : allocate(entropy(1:nz))
908 0 : do k=1,nz
909 0 : entropy(k) = exp(s% lnS(k))
910 : end do
911 :
912 0 : dentropy_sum = sum(abs((entropy(1:nz)-vals(1:nz))/vals(1:nz)))
913 0 : avg_err = dentropy_sum/nz
914 0 : deallocate(vals, xq, entropy)
915 0 : end subroutine adjust_entropy
916 :
917 : end function relax_entropy_check_model
918 :
919 0 : subroutine entropy_relax_other_energy(id, ierr)
920 : use interp_1d_lib, only: interp_values
921 : use auto_diff_support
922 : integer, intent(in) :: id
923 : integer, intent(out) :: ierr
924 : type (star_info), pointer :: s
925 : integer :: k, nz, num_pts
926 0 : real(dp), pointer :: vals(:), xq(:), x(:), f(:)
927 : ierr = 0
928 0 : call star_ptr(id, s, ierr)
929 :
930 0 : nz = s% nz
931 0 : num_pts = relax_num_pts
932 0 : allocate(vals(nz), xq(nz), stat=ierr)
933 0 : f(1:4*num_pts) => relax_work_array(num_pts+1:5*num_pts)
934 0 : x(1:num_pts) => relax_work_array(1:num_pts)
935 0 : xq(1) = s% dq(1)/2 ! xq for cell center
936 0 : do k = 2, nz
937 0 : xq(k) = xq(k-1) + (s% dq(k) + s% dq(k-1))/2
938 : end do
939 0 : call interp_values(x, num_pts, f, nz, xq, vals(:), ierr)
940 0 : if (ierr /= 0) return
941 0 : do k = 1, s% nz
942 0 : s% extra_heat(k) = ( 1d0 - exp(s%lnS(k))/vals(k) ) * exp(s%lnE(k))
943 0 : s% extra_heat(k) = s% extra_heat(k) / (s% job% timescale_for_relax_entropy * secyer)
944 : end do
945 0 : end subroutine entropy_relax_other_energy
946 :
947 0 : subroutine do_relax_angular_momentum( &
948 0 : id, max_steps_to_use, num_pts, angular_momentum, xq, ierr)
949 0 : use alloc, only: alloc_extras, dealloc_extras
950 : integer, intent(in) :: id
951 : integer, intent(in) :: max_steps_to_use ! use this many steps to do conversion
952 : integer, intent(in) :: num_pts
953 : ! length of angular momentum vector; need not equal nz for current model
954 : real(dp), intent(in) :: angular_momentum(num_pts) ! desired angular momentum profile
955 : real(dp), intent(in) :: xq(num_pts)
956 : integer, intent(out) :: ierr
957 :
958 : integer, parameter :: lipar=2
959 : integer :: lrpar, max_model_number
960 0 : real(dp), pointer :: rpar(:)
961 0 : real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor, max_timestep, &
962 0 : am_D_mix_factor, time
963 : logical :: do_element_diffusion, use_other_torque
964 : type (star_info), pointer :: s
965 0 : real(dp), pointer :: x(:), f1(:), f(:,:)
966 : integer, target :: ipar_ary(lipar)
967 : integer, pointer :: ipar(:)
968 : procedure (other_torque_interface), pointer :: &
969 : other_torque => null()
970 :
971 0 : ipar => ipar_ary
972 :
973 : include 'formats'
974 :
975 : ierr = 0
976 :
977 0 : call get_star_ptr(id, s, ierr)
978 0 : if (ierr /= 0) return
979 :
980 0 : max_model_number = s% max_model_number
981 0 : s% max_model_number = max_steps_to_use + s% model_number + 1
982 0 : write(*,*) 'relax_angular_momentum: max_steps_to_use', max_steps_to_use
983 :
984 0 : ipar(1) = num_pts
985 0 : ipar(2) = s% max_model_number
986 0 : lrpar = 5*num_pts
987 0 : allocate(rpar(lrpar), stat=ierr)
988 0 : if (ierr /= 0) return
989 :
990 0 : x(1:num_pts) => rpar(1:num_pts)
991 0 : f1(1:4*num_pts) => rpar(num_pts+1:lrpar)
992 0 : f(1:4,1:num_pts) => f1(1:4*num_pts)
993 0 : call store_rpar(num_pts, ierr)
994 0 : if (ierr /= 0) return
995 :
996 : ! need to use global variables, as relax_angular_momentum uses
997 : ! the other_torque routine to which it can't pass rpar
998 0 : relax_num_pts = num_pts
999 0 : relax_work_array => rpar
1000 :
1001 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
1002 0 : s% dxdt_nuc_factor = 0 ! turn off composition change by nuclear burning
1003 0 : mix_factor = s% mix_factor
1004 0 : s% mix_factor = 0 ! turn off mixing
1005 0 : am_D_mix_factor = s% am_D_mix_factor
1006 0 : s% am_D_mix_factor = 0d0
1007 0 : do_element_diffusion = s% do_element_diffusion
1008 0 : s% do_element_diffusion = .false. ! turn off diffusion
1009 0 : starting_dt_next = s% dt_next
1010 0 : max_timestep = s% max_timestep
1011 0 : s% max_timestep = s% job% max_dt_for_relax_angular_momentum * secyer
1012 0 : s% dt_next = min(s% dt_next, s% job% max_dt_for_relax_angular_momentum * secyer)
1013 0 : use_other_torque = s% use_other_torque
1014 0 : s% use_other_torque = .true.
1015 0 : other_torque => s% other_torque
1016 0 : s% other_torque => angular_momentum_relax_other_torque
1017 0 : time = s% time
1018 0 : s% time = 0d0
1019 :
1020 : call do_internal_evolve( &
1021 : id, before_evolve_relax_angular_momentum, &
1022 : relax_angular_momentum_adjust_model, relax_angular_momentum_check_model, &
1023 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1024 :
1025 0 : s% max_model_number = max_model_number
1026 0 : s% dt_next = starting_dt_next
1027 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
1028 0 : s% mix_factor = mix_factor
1029 0 : s% am_D_mix_factor = am_D_mix_factor
1030 0 : s% do_element_diffusion = do_element_diffusion
1031 0 : s% max_timestep = max_timestep
1032 0 : s% use_other_torque = use_other_torque
1033 0 : s% other_torque => other_torque
1034 0 : s% time = time
1035 :
1036 0 : call error_check('relax angular momentum',ierr)
1037 :
1038 0 : deallocate(rpar)
1039 :
1040 :
1041 : contains
1042 :
1043 :
1044 0 : subroutine store_rpar(num_pts, ierr) ! get interpolation info
1045 0 : use interp_1d_def, only: pm_work_size
1046 : use interp_1d_lib, only: interp_pm
1047 : integer, intent(in) :: num_pts
1048 : integer, intent(out) :: ierr
1049 : integer :: op_err
1050 0 : real(dp), pointer :: interp_work(:), work(:), p(:)
1051 0 : allocate(interp_work(num_pts*pm_work_size), stat=ierr)
1052 0 : if (ierr /= 0) return
1053 0 : x(:) = xq(:)
1054 : op_err = 0
1055 0 : f(1,1:num_pts) = angular_momentum(1:num_pts)
1056 : work(1:num_pts*pm_work_size) => &
1057 0 : interp_work(1:num_pts*pm_work_size)
1058 0 : p(1:4*num_pts) => f1(1:4*num_pts)
1059 : call interp_pm(x, num_pts, p, pm_work_size, work, &
1060 0 : 'do_relax_entropy', op_err)
1061 0 : if (op_err /= 0) ierr = op_err
1062 0 : deallocate(interp_work)
1063 0 : end subroutine store_rpar
1064 :
1065 : end subroutine do_relax_angular_momentum
1066 :
1067 0 : subroutine before_evolve_relax_angular_momentum(s, id, lipar, ipar, lrpar, rpar, ierr)
1068 : type (star_info), pointer :: s
1069 : integer, intent(in) :: id, lipar, lrpar
1070 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1071 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1072 : integer, intent(out) :: ierr
1073 0 : ierr = 0
1074 0 : call setup_before_relax(s)
1075 0 : end subroutine before_evolve_relax_angular_momentum
1076 :
1077 0 : integer function relax_angular_momentum_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1078 : type (star_info), pointer :: s
1079 : integer, intent(in) :: id, lipar, lrpar
1080 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1081 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1082 0 : relax_angular_momentum_adjust_model = keep_going
1083 0 : end function relax_angular_momentum_adjust_model
1084 :
1085 0 : integer function relax_angular_momentum_check_model(s, id, lipar, ipar, lrpar, rpar)
1086 : use do_one_utils, only: do_bare_bones_check_model
1087 : type (star_info), pointer :: s
1088 : integer, intent(in) :: id, lipar, lrpar
1089 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1090 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1091 :
1092 : integer :: num_pts, ierr, max_model_number
1093 0 : real(dp) :: avg_err
1094 0 : real(dp), pointer :: x(:) ! =(num_pts)
1095 0 : real(dp), pointer :: f1(:) ! =(4, num_pts)
1096 :
1097 : include 'formats'
1098 :
1099 0 : relax_angular_momentum_check_model = do_bare_bones_check_model(id)
1100 0 : if (relax_angular_momentum_check_model /= keep_going) return
1101 :
1102 0 : if (lipar /= 2) then
1103 0 : write(*,*) 'bad lipar for relax_angular_momentum_check_model'
1104 0 : relax_angular_momentum_check_model = terminate
1105 0 : return
1106 : end if
1107 :
1108 0 : num_pts = ipar(1)
1109 0 : max_model_number = ipar(2)
1110 :
1111 0 : if (lrpar /= 5*num_pts) then
1112 0 : write(*,*) 'bad lrpar for relax_angular_momentum_check_model'
1113 0 : relax_angular_momentum_check_model = terminate
1114 : end if
1115 :
1116 : ierr = 0
1117 0 : x(1:num_pts) => rpar(1:num_pts)
1118 0 : f1(1:4*num_pts) => rpar(num_pts+1:lrpar)
1119 0 : call adjust_angular_momentum(num_pts, avg_err, ierr)
1120 0 : if (ierr /= 0) relax_angular_momentum_check_model = terminate
1121 :
1122 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1123 0 : write(*,*) 'relax_angular_momentum avg rel err, dt, model', avg_err, s% dt/secyer, s% model_number
1124 :
1125 0 : if (s% star_age >= &
1126 : s% job% timescale_for_relax_angular_momentum*&
1127 : s% job% num_timescales_for_relax_angular_momentum) then
1128 0 : relax_angular_momentum_check_model = terminate
1129 0 : s% termination_code = t_relax_finished_okay
1130 0 : return
1131 : end if
1132 :
1133 0 : if (s% model_number >= max_model_number) then
1134 0 : write(*,*) "Terminated relax because of max_model_number instead of relax_time"
1135 0 : relax_angular_momentum_check_model = terminate
1136 0 : s% termination_code = t_relax_finished_okay
1137 0 : return
1138 : end if
1139 :
1140 :
1141 : contains
1142 :
1143 :
1144 0 : subroutine adjust_angular_momentum(num_pts, avg_err, ierr)
1145 0 : use interp_1d_lib, only: interp_values
1146 : integer, intent(in) :: num_pts
1147 : real(dp), intent(out) :: avg_err
1148 : integer, intent(out) :: ierr
1149 : integer :: k, nz, op_err
1150 0 : real(dp) :: dangular_momentum_sum
1151 0 : real(dp), pointer :: vals(:), xq(:), f(:)
1152 0 : ierr = 0
1153 0 : nz = s% nz
1154 0 : allocate(vals(nz), xq(nz), stat=ierr)
1155 0 : if (ierr /= 0) return
1156 0 : xq(1) = s% dq(1)/2 ! xq for cell center
1157 0 : do k = 2, nz
1158 0 : xq(k) = xq(k-1) + (s% dq(k) + s% dq(k-1))/2
1159 : end do
1160 0 : dangular_momentum_sum = 0
1161 0 : f(1:4*num_pts) => f1(1:4*num_pts)
1162 0 : call interp_values(x, num_pts, f, nz, xq, vals(:), op_err)
1163 0 : if (op_err /= 0) ierr = op_err
1164 0 : dangular_momentum_sum = sum(abs((s% j_rot(1:nz)-vals(1:nz))/vals(1:nz)))
1165 0 : avg_err = dangular_momentum_sum/nz
1166 0 : deallocate(vals, xq)
1167 0 : end subroutine adjust_angular_momentum
1168 :
1169 : end function relax_angular_momentum_check_model
1170 :
1171 0 : subroutine angular_momentum_relax_other_torque(id, ierr)
1172 : use interp_1d_lib, only: interp_values
1173 : integer, intent(in) :: id
1174 : integer, intent(out) :: ierr
1175 : type (star_info), pointer :: s
1176 : integer :: k, k_inner, nz, num_pts, op_err
1177 0 : real(dp), pointer :: vals(:), xq(:), x(:), f(:)
1178 0 : real(dp) :: omega_target, j_rot_target
1179 : ierr = 0
1180 0 : call star_ptr(id, s, ierr)
1181 :
1182 0 : nz = s% nz
1183 0 : num_pts = relax_num_pts
1184 0 : allocate(vals(nz), xq(nz), stat=ierr)
1185 0 : f(1:4*num_pts) => relax_work_array(num_pts+1:5*num_pts)
1186 0 : x(1:num_pts) => relax_work_array(1:num_pts)
1187 0 : xq(1) = s% dq(1)/2 ! xq for cell center
1188 0 : do k = 2, nz
1189 0 : xq(k) = xq(k-1) + (s% dq(k) + s% dq(k-1))/2
1190 : end do
1191 0 : call interp_values(x, num_pts, f, nz, xq, vals(:), op_err)
1192 0 : if (op_err /= 0) ierr = op_err
1193 :
1194 0 : if (ierr /= 0) return
1195 0 : s% extra_jdot(:) = 0d0
1196 0 : do k = 1, s% nz
1197 0 : if (s% j_rot_flag) then
1198 : s% extra_jdot(k) = &
1199 : (1d0 - exp(-s% dt/(s% job% timescale_for_relax_angular_momentum*secyer))) * &
1200 0 : (vals(k) - s% j_rot_start(k))/s% dt
1201 : else
1202 : s% extra_jdot(k) = &
1203 : (1d0 - exp(-s% dt/(s% job% timescale_for_relax_angular_momentum*secyer))) * &
1204 0 : (vals(k) - s% j_rot(k))/s% dt
1205 : end if
1206 : end do
1207 : ! for cells near center use a constant omega, prevents extrapolating j_rot and causing artificially large core rotation
1208 : k_inner = -1
1209 0 : do k = 2, s% nz
1210 0 : if (xq(k) > x(num_pts)) then
1211 0 : k_inner = k-1
1212 : end if
1213 : end do
1214 0 : if (s% job% relax_angular_momentum_constant_omega_center) then
1215 0 : if (k_inner > 0) then
1216 0 : omega_target = vals(k_inner)/s% i_rot(k_inner)% val
1217 0 : do k=k_inner+1, s% nz
1218 0 : j_rot_target = omega_target*s% i_rot(k)% val
1219 0 : if (s% j_rot_flag) then
1220 : s% extra_jdot(k) = &
1221 : (1d0 - exp(-s% dt/(s% job% timescale_for_relax_angular_momentum*secyer))) * &
1222 0 : (j_rot_target - s% j_rot_start(k))/s% dt
1223 : else
1224 : s% extra_jdot(k) = &
1225 : (1d0 - exp(-s% dt/(s% job% timescale_for_relax_angular_momentum*secyer))) * &
1226 0 : (j_rot_target - s% j_rot(k))/s% dt
1227 : end if
1228 : end do
1229 : end if
1230 : end if
1231 0 : deallocate(vals, xq)
1232 0 : end subroutine angular_momentum_relax_other_torque
1233 :
1234 0 : subroutine do_relax_uniform_omega( &
1235 : id, kind_of_relax, target_value, num_steps_to_relax_rotation,&
1236 : relax_omega_max_yrs_dt, ierr)
1237 : integer, intent(in) :: id, kind_of_relax
1238 : real(dp), intent(in) :: target_value,relax_omega_max_yrs_dt
1239 : integer, intent(in) :: num_steps_to_relax_rotation
1240 : integer, intent(out) :: ierr
1241 :
1242 : integer, parameter :: lipar=3, lrpar=2
1243 : integer :: max_model_number, k
1244 0 : real(dp) :: max_years_for_timestep, mix_factor, dxdt_nuc_factor, am_D_mix_factor
1245 : type (star_info), pointer :: s
1246 0 : real(dp), target :: rpar_ary(lrpar)
1247 : real(dp), pointer :: rpar(:)
1248 : integer, target :: ipar_ary(lipar)
1249 : integer, pointer :: ipar(:)
1250 : logical :: okay
1251 :
1252 : include 'formats'
1253 :
1254 : ierr = 0
1255 :
1256 0 : rpar => rpar_ary
1257 0 : ipar => ipar_ary
1258 :
1259 0 : call get_star_ptr(id, s, ierr)
1260 0 : if (ierr /= 0) return
1261 :
1262 0 : if (.not. s% rotation_flag) return
1263 :
1264 0 : rpar(1) = target_value
1265 0 : okay = .true.
1266 0 : do k=1,s% nz
1267 0 : if (s% omega(k) /= target_value) then
1268 : okay = .false.
1269 : exit
1270 : end if
1271 : end do
1272 0 : if (okay) return
1273 :
1274 0 : rpar(2) = s% omega(1)
1275 0 : ipar(1) = s% model_number
1276 0 : ipar(2) = num_steps_to_relax_rotation
1277 0 : ipar(3) = kind_of_relax
1278 0 : if (num_steps_to_relax_rotation <= 0) then
1279 0 : ierr = -1
1280 0 : write(*,2) 'invalid num_steps_to_relax_rotation', num_steps_to_relax_rotation
1281 0 : return
1282 : end if
1283 :
1284 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
1285 0 : s% dxdt_nuc_factor = 0d0 ! turn off composition change by nuclear burning
1286 0 : mix_factor = s% mix_factor
1287 0 : am_D_mix_factor = s% am_D_mix_factor
1288 0 : s% mix_factor = 0d0
1289 0 : s% am_D_mix_factor = 0d0
1290 0 : max_model_number = s% max_model_number
1291 0 : s% max_model_number = num_steps_to_relax_rotation + 1
1292 0 : max_years_for_timestep = s% max_years_for_timestep
1293 0 : s% max_years_for_timestep = relax_omega_max_yrs_dt
1294 0 : write(*,*) 'num_steps_to_relax_rotation', num_steps_to_relax_rotation
1295 :
1296 : call do_internal_evolve( &
1297 : id, before_evolve_relax_omega, &
1298 : relax_omega_adjust_model, relax_omega_check_model, &
1299 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1300 :
1301 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
1302 0 : s% mix_factor = mix_factor
1303 0 : s% am_D_mix_factor = am_D_mix_factor
1304 0 : s% max_model_number = max_model_number
1305 0 : s% max_years_for_timestep = max_years_for_timestep
1306 0 : call error_check('relax uniform omega',ierr)
1307 0 : s% D_omega(1:s% nz) = 0
1308 0 : s% am_nu_rot(1:s% nz) = 0
1309 :
1310 0 : end subroutine do_relax_uniform_omega
1311 :
1312 :
1313 0 : subroutine before_evolve_relax_omega(s, id, lipar, ipar, lrpar, rpar, ierr)
1314 : type (star_info), pointer :: s
1315 : integer, intent(in) :: id, lipar, lrpar
1316 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1317 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1318 : integer, intent(out) :: ierr
1319 0 : ierr = 0
1320 0 : call setup_before_relax(s)
1321 0 : end subroutine before_evolve_relax_omega
1322 :
1323 0 : integer function relax_omega_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1324 : type (star_info), pointer :: s
1325 : integer, intent(in) :: id, lipar, lrpar
1326 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1327 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1328 0 : relax_omega_adjust_model = keep_going
1329 0 : end function relax_omega_adjust_model
1330 :
1331 0 : integer function relax_omega_check_model(s, id, lipar, ipar, lrpar, rpar)
1332 : use do_one_utils, only: do_bare_bones_check_model
1333 : use hydro_rotation, only: set_uniform_omega, set_i_rot, &
1334 : set_surf_avg_rotation_info
1335 : type (star_info), pointer :: s
1336 : integer, intent(in) :: id, lipar, lrpar
1337 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1338 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1339 :
1340 : integer :: num_steps_to_use, starting_model_number, kind_of_relax, ierr
1341 0 : real(dp) :: frac, target_value, starting_omega, new_omega, this_step_omega
1342 :
1343 : include 'formats'
1344 :
1345 0 : relax_omega_check_model = do_bare_bones_check_model(id)
1346 0 : if (relax_omega_check_model /= keep_going) return
1347 :
1348 0 : starting_model_number = ipar(1)
1349 0 : num_steps_to_use = ipar(2)
1350 0 : kind_of_relax = ipar(3)
1351 0 : target_value = rpar(1)
1352 0 : starting_omega = rpar(2)
1353 :
1354 : frac = max(0d0, min(1d0, &
1355 0 : dble(s% model_number - starting_model_number)/dble(num_steps_to_use)))
1356 :
1357 0 : if (kind_of_relax == relax_to_new_omega) then
1358 0 : new_omega = target_value
1359 0 : else if (kind_of_relax == relax_to_new_omega_div_omega_crit) then
1360 0 : call set_i_rot(s, .false.)
1361 0 : call set_surf_avg_rotation_info(s)
1362 0 : new_omega = target_value*s% omega_crit_avg_surf
1363 0 : else if (kind_of_relax == relax_to_new_surface_rotation_v) then
1364 0 : new_omega = target_value*1d5/(s% photosphere_r*Rsun)
1365 : else
1366 0 : write(*,2) 'bad value for kind_of_relax', kind_of_relax
1367 0 : call mesa_error(__FILE__,__LINE__,'relax_omega_check_model')
1368 : end if
1369 :
1370 0 : this_step_omega = frac*new_omega + (1 - frac)*starting_omega
1371 :
1372 0 : if (s% model_number > starting_model_number + num_steps_to_use + 50 .or. &
1373 : abs(s% omega(1) - new_omega) < 1d-4*max(1d-10,abs(new_omega))) then
1374 0 : write(*,2) 'final step: wanted-current, current, wanted', &
1375 0 : s% model_number, new_omega-s% omega(1), s% omega(1), new_omega
1376 0 : relax_omega_check_model = terminate
1377 0 : s% termination_code = t_relax_finished_okay
1378 0 : else if (mod(s% model_number, s% terminal_interval) == 0) then
1379 0 : write(*,2) 'relax to omega: wanted-current, current, wanted', &
1380 0 : s% model_number, new_omega-s% omega(1), s% omega(1), new_omega
1381 : end if
1382 :
1383 : ierr = 0
1384 0 : call set_uniform_omega(id, this_step_omega, ierr)
1385 0 : if (ierr /= 0) then
1386 0 : write(*,*) 'set_uniform_omega failed'
1387 0 : relax_omega_check_model = terminate
1388 0 : return
1389 : end if
1390 :
1391 0 : end function relax_omega_check_model
1392 :
1393 :
1394 0 : subroutine do_relax_tau_factor(id, new_tau_factor, dlogtau_factor, ierr)
1395 : integer, intent(in) :: id
1396 : real(dp), intent(in) :: new_tau_factor, dlogtau_factor
1397 : integer, intent(out) :: ierr
1398 : integer, parameter :: lipar=1, lrpar=2
1399 : integer :: max_model_number
1400 0 : real(dp) :: tau_factor
1401 : type (star_info), pointer :: s
1402 : integer, target :: ipar_ary(lipar)
1403 : integer, pointer :: ipar(:)
1404 0 : real(dp), target :: rpar_ary(lrpar)
1405 : real(dp), pointer :: rpar(:)
1406 0 : rpar => rpar_ary
1407 0 : ipar => ipar_ary
1408 : include 'formats'
1409 : ierr = 0
1410 0 : if (new_tau_factor <= 0) then
1411 0 : ierr = -1
1412 0 : write(*,*) 'invalid new_tau_factor', new_tau_factor
1413 0 : return
1414 : end if
1415 0 : call get_star_ptr(id, s, ierr)
1416 0 : if (ierr /= 0) return
1417 0 : tau_factor = s% tau_factor
1418 0 : if (abs(new_tau_factor - tau_factor) <= 1d-6) then
1419 0 : s% tau_factor = new_tau_factor
1420 0 : return
1421 : end if
1422 0 : write(*,'(A)')
1423 0 : write(*,1) 'current tau_factor', tau_factor
1424 0 : write(*,1) 'relax to new tau_factor', new_tau_factor
1425 0 : write(*,'(A)')
1426 0 : write(*,1) 'dlogtau_factor', dlogtau_factor
1427 0 : write(*,'(A)')
1428 0 : rpar(1) = new_tau_factor
1429 0 : rpar(2) = dlogtau_factor
1430 0 : max_model_number = s% max_model_number
1431 0 : s% max_model_number = -1111
1432 : call do_internal_evolve( &
1433 : id, before_evolve_relax_tau_factor, &
1434 : relax_tau_factor_adjust_model, relax_tau_factor_check_model, &
1435 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1436 0 : s% max_model_number = max_model_number
1437 0 : if (ierr == 0) s% force_tau_factor = s% tau_factor
1438 0 : call error_check('relax tau factor',ierr)
1439 0 : end subroutine do_relax_tau_factor
1440 :
1441 :
1442 0 : subroutine before_evolve_relax_tau_factor(s, id, lipar, ipar, lrpar, rpar, ierr)
1443 : type (star_info), pointer :: s
1444 : integer, intent(in) :: id, lipar, lrpar
1445 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1446 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1447 : integer, intent(out) :: ierr
1448 0 : ierr = 0
1449 0 : call turn_off_winds(s)
1450 0 : s% max_model_number = -111
1451 0 : end subroutine before_evolve_relax_tau_factor
1452 :
1453 0 : integer function relax_tau_factor_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1454 : type (star_info), pointer :: s
1455 : integer, intent(in) :: id, lipar, lrpar
1456 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1457 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1458 0 : relax_tau_factor_adjust_model = keep_going
1459 0 : end function relax_tau_factor_adjust_model
1460 :
1461 0 : integer function relax_tau_factor_check_model(s, id, lipar, ipar, lrpar, rpar)
1462 : use do_one_utils, only:do_bare_bones_check_model
1463 : type (star_info), pointer :: s
1464 : integer, intent(in) :: id, lipar, lrpar
1465 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1466 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1467 0 : real(dp) :: new_tau_factor, dlogtau_factor, current_tau_factor, next
1468 : logical, parameter :: dbg = .false.
1469 :
1470 : include 'formats'
1471 :
1472 0 : relax_tau_factor_check_model = do_bare_bones_check_model(id)
1473 0 : if (relax_tau_factor_check_model /= keep_going) return
1474 :
1475 0 : new_tau_factor = rpar(1)
1476 0 : dlogtau_factor = rpar(2)
1477 0 : current_tau_factor = s% tau_factor
1478 :
1479 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1480 0 : write(*,1) 'tau_factor target current', new_tau_factor, current_tau_factor
1481 :
1482 0 : if (abs(current_tau_factor-new_tau_factor) < 1d-15) then
1483 0 : s% tau_factor = new_tau_factor
1484 0 : s% termination_code = t_relax_finished_okay
1485 0 : relax_tau_factor_check_model = terminate
1486 0 : return
1487 : end if
1488 :
1489 0 : if (new_tau_factor < current_tau_factor) then
1490 0 : next = exp10(safe_log10(current_tau_factor) - dlogtau_factor)
1491 0 : if (next < new_tau_factor) next = new_tau_factor
1492 : else
1493 0 : next = exp10(safe_log10(current_tau_factor) + dlogtau_factor)
1494 0 : if (next > new_tau_factor) next = new_tau_factor
1495 : end if
1496 :
1497 : if (dbg) write(*,1) 'next', next, log10(next)
1498 :
1499 0 : s% tau_factor = next
1500 0 : s% max_timestep = secyer*s% time_step
1501 :
1502 0 : end function relax_tau_factor_check_model
1503 :
1504 :
1505 0 : subroutine do_relax_opacity_factor(id, new_opacity_factor, dopacity_factor, ierr)
1506 : integer, intent(in) :: id
1507 : real(dp), intent(in) :: new_opacity_factor, dopacity_factor
1508 : integer, intent(out) :: ierr
1509 : integer, parameter :: lipar=1, lrpar=2
1510 : integer :: max_model_number
1511 0 : real(dp) :: opacity_factor
1512 : type (star_info), pointer :: s
1513 : integer, target :: ipar_ary(lipar)
1514 : integer, pointer :: ipar(:)
1515 0 : real(dp), target :: rpar_ary(lrpar)
1516 : real(dp), pointer :: rpar(:)
1517 0 : rpar => rpar_ary
1518 0 : ipar => ipar_ary
1519 : include 'formats'
1520 : ierr = 0
1521 0 : if (new_opacity_factor <= 0) then
1522 0 : ierr = -1
1523 0 : write(*,*) 'invalid new_opacity_factor', new_opacity_factor
1524 0 : return
1525 : end if
1526 0 : call get_star_ptr(id, s, ierr)
1527 0 : if (ierr /= 0) return
1528 0 : opacity_factor = s% opacity_factor
1529 0 : if (abs(new_opacity_factor - opacity_factor) <= 1d-6) then
1530 0 : s% opacity_factor = new_opacity_factor
1531 0 : return
1532 : end if
1533 0 : write(*,'(A)')
1534 0 : write(*,1) 'current opacity_factor', opacity_factor
1535 0 : write(*,1) 'relax to new opacity_factor', new_opacity_factor
1536 0 : write(*,'(A)')
1537 0 : write(*,1) 'dopacity_factor', dopacity_factor
1538 0 : write(*,'(A)')
1539 0 : rpar(1) = new_opacity_factor
1540 0 : rpar(2) = dopacity_factor
1541 0 : max_model_number = s% max_model_number
1542 0 : s% max_model_number = -1111
1543 : call do_internal_evolve( &
1544 : id, before_evolve_relax_opacity_factor, &
1545 : relax_opacity_factor_adjust_model, relax_opacity_factor_check_model, &
1546 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1547 0 : s% max_model_number = max_model_number
1548 0 : if (ierr == 0) s% force_opacity_factor = s% opacity_factor
1549 0 : call error_check('relax opacity factor',ierr)
1550 0 : end subroutine do_relax_opacity_factor
1551 :
1552 :
1553 0 : subroutine before_evolve_relax_opacity_factor(s, id, lipar, ipar, lrpar, rpar, ierr)
1554 : type (star_info), pointer :: s
1555 : integer, intent(in) :: id, lipar, lrpar
1556 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1557 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1558 : integer, intent(out) :: ierr
1559 0 : ierr = 0
1560 0 : call turn_off_winds(s)
1561 0 : s% max_model_number = -111
1562 0 : end subroutine before_evolve_relax_opacity_factor
1563 :
1564 0 : integer function relax_opacity_factor_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1565 : type (star_info), pointer :: s
1566 : integer, intent(in) :: id, lipar, lrpar
1567 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1568 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1569 0 : relax_opacity_factor_adjust_model = keep_going
1570 0 : end function relax_opacity_factor_adjust_model
1571 :
1572 0 : integer function relax_opacity_factor_check_model(s, id, lipar, ipar, lrpar, rpar)
1573 : use do_one_utils, only:do_bare_bones_check_model
1574 : type (star_info), pointer :: s
1575 : integer, intent(in) :: id, lipar, lrpar
1576 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1577 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1578 0 : real(dp) :: new_opacity_factor, dopacity_factor, current_opacity_factor, next
1579 : logical, parameter :: dbg = .false.
1580 :
1581 : include 'formats'
1582 :
1583 0 : relax_opacity_factor_check_model = do_bare_bones_check_model(id)
1584 0 : if (relax_opacity_factor_check_model /= keep_going) return
1585 :
1586 0 : new_opacity_factor = rpar(1)
1587 0 : dopacity_factor = rpar(2)
1588 0 : current_opacity_factor = s% opacity_factor
1589 :
1590 : if (dbg) then
1591 : write(*,1) 'new_opacity_factor', new_opacity_factor
1592 : write(*,1) 'current_opacity_factor', current_opacity_factor
1593 : end if
1594 :
1595 0 : if (abs(current_opacity_factor-new_opacity_factor) < 1d-15) then
1596 0 : s% opacity_factor = new_opacity_factor
1597 0 : s% termination_code = t_relax_finished_okay
1598 0 : relax_opacity_factor_check_model = terminate
1599 0 : return
1600 : end if
1601 :
1602 0 : if (new_opacity_factor < current_opacity_factor) then
1603 0 : next = current_opacity_factor - dopacity_factor
1604 0 : if (next < new_opacity_factor) next = new_opacity_factor
1605 : else
1606 0 : next = current_opacity_factor + dopacity_factor
1607 0 : if (next > new_opacity_factor) next = new_opacity_factor
1608 : end if
1609 :
1610 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1611 0 : write(*,1) 'next opacity_factor', next ! provide terminal feedback to show working.
1612 :
1613 0 : s% opacity_factor = next
1614 0 : s% max_timestep = secyer*s% time_step
1615 :
1616 0 : end function relax_opacity_factor_check_model
1617 :
1618 :
1619 0 : subroutine do_relax_Tsurf_factor(id, new_Tsurf_factor, dlogTsurf_factor, ierr)
1620 : integer, intent(in) :: id
1621 : real(dp), intent(in) :: new_Tsurf_factor, dlogTsurf_factor
1622 : integer, intent(out) :: ierr
1623 : integer, parameter :: lipar=1, lrpar=2
1624 : integer :: max_model_number
1625 0 : real(dp) :: Tsurf_factor
1626 : type (star_info), pointer :: s
1627 : integer, target :: ipar_ary(lipar)
1628 : integer, pointer :: ipar(:)
1629 0 : real(dp), target :: rpar_ary(lrpar)
1630 : real(dp), pointer :: rpar(:)
1631 0 : rpar => rpar_ary
1632 0 : ipar => ipar_ary
1633 : include 'formats'
1634 : ierr = 0
1635 0 : if (new_Tsurf_factor <= 0) then
1636 0 : ierr = -1
1637 0 : write(*,*) 'invalid new_Tsurf_factor', new_Tsurf_factor
1638 0 : return
1639 : end if
1640 0 : call get_star_ptr(id, s, ierr)
1641 0 : if (ierr /= 0) return
1642 0 : Tsurf_factor = s% Tsurf_factor
1643 0 : if (abs(new_Tsurf_factor - Tsurf_factor) <= 1d-6) then
1644 0 : s% Tsurf_factor = new_Tsurf_factor
1645 0 : return
1646 : end if
1647 0 : write(*,'(A)')
1648 0 : write(*,1) 'current Tsurf_factor', Tsurf_factor
1649 0 : write(*,1) 'relax to new Tsurf_factor', new_Tsurf_factor
1650 0 : write(*,'(A)')
1651 0 : write(*,1) 'dlogTsurf_factor', dlogTsurf_factor
1652 0 : write(*,'(A)')
1653 0 : rpar(1) = new_Tsurf_factor
1654 0 : rpar(2) = dlogTsurf_factor
1655 0 : max_model_number = s% max_model_number
1656 0 : s% max_model_number = -1111
1657 : call do_internal_evolve( &
1658 : id, before_evolve_relax_Tsurf_factor, &
1659 : relax_Tsurf_factor_adjust_model, relax_Tsurf_factor_check_model, &
1660 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1661 0 : s% max_model_number = max_model_number
1662 0 : if (ierr /= 0) then
1663 0 : write(*,'(A)')
1664 0 : write(*,1) 'ERROR: failed doing relax Tsurf_factor', new_Tsurf_factor
1665 0 : write(*,'(A)')
1666 0 : call mesa_error(__FILE__,__LINE__,'do_relax_Tsurf_factor')
1667 : end if
1668 :
1669 0 : if (new_Tsurf_factor == 1d0) then
1670 0 : s% force_Tsurf_factor = 0d0
1671 : else
1672 0 : s% force_Tsurf_factor = s% Tsurf_factor
1673 : end if
1674 :
1675 0 : call error_check('relax tsurf factor',ierr)
1676 :
1677 0 : end subroutine do_relax_Tsurf_factor
1678 :
1679 :
1680 0 : subroutine before_evolve_relax_Tsurf_factor(s, id, lipar, ipar, lrpar, rpar, ierr)
1681 : type (star_info), pointer :: s
1682 : integer, intent(in) :: id, lipar, lrpar
1683 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1684 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1685 : integer, intent(out) :: ierr
1686 0 : ierr = 0
1687 0 : call turn_off_winds(s)
1688 0 : s% max_model_number = -111
1689 0 : end subroutine before_evolve_relax_Tsurf_factor
1690 :
1691 0 : integer function relax_Tsurf_factor_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1692 : type (star_info), pointer :: s
1693 : integer, intent(in) :: id, lipar, lrpar
1694 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1695 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1696 0 : relax_Tsurf_factor_adjust_model = keep_going
1697 0 : end function relax_Tsurf_factor_adjust_model
1698 :
1699 0 : integer function relax_Tsurf_factor_check_model(s, id, lipar, ipar, lrpar, rpar)
1700 : use do_one_utils, only:do_bare_bones_check_model
1701 : type (star_info), pointer :: s
1702 : integer, intent(in) :: id, lipar, lrpar
1703 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1704 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1705 0 : real(dp) :: new_Tsurf_factor, dlogTsurf_factor, current_Tsurf_factor, next
1706 : logical, parameter :: dbg = .false.
1707 :
1708 : include 'formats'
1709 :
1710 0 : relax_Tsurf_factor_check_model = do_bare_bones_check_model(id)
1711 0 : if (relax_Tsurf_factor_check_model /= keep_going) return
1712 :
1713 0 : new_Tsurf_factor = rpar(1)
1714 0 : dlogTsurf_factor = rpar(2)
1715 0 : current_Tsurf_factor = s% Tsurf_factor
1716 :
1717 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1718 0 : write(*,1) 'Tsurf_factor target current', new_Tsurf_factor, current_Tsurf_factor
1719 :
1720 0 : if (abs(current_Tsurf_factor-new_Tsurf_factor) < 1d-15) then
1721 0 : s% Tsurf_factor = new_Tsurf_factor
1722 0 : s% termination_code = t_relax_finished_okay
1723 0 : relax_Tsurf_factor_check_model = terminate
1724 0 : return
1725 : end if
1726 :
1727 0 : if (new_Tsurf_factor < current_Tsurf_factor) then
1728 0 : next = exp10(safe_log10(current_Tsurf_factor) - dlogTsurf_factor)
1729 0 : if (next < new_Tsurf_factor) next = new_Tsurf_factor
1730 : else
1731 0 : next = exp10(safe_log10(current_Tsurf_factor) + dlogTsurf_factor)
1732 0 : if (next > new_Tsurf_factor) next = new_Tsurf_factor
1733 : end if
1734 :
1735 : if (dbg) write(*,1) 'next Tsurf_factor', next, log10(next)
1736 :
1737 0 : s% Tsurf_factor = next
1738 0 : s% max_timestep = secyer*s% time_step
1739 :
1740 0 : end function relax_Tsurf_factor_check_model
1741 :
1742 :
1743 0 : subroutine do_relax_irradiation(id, &
1744 : min_steps, new_irrad_flux, new_irrad_col_depth, relax_irradiation_max_yrs_dt, ierr)
1745 : integer, intent(in) :: id, min_steps
1746 : real(dp), intent(in) :: &
1747 : new_irrad_flux, new_irrad_col_depth, relax_irradiation_max_yrs_dt
1748 : integer, intent(out) :: ierr
1749 :
1750 : integer, parameter :: lipar=2, lrpar=4
1751 : integer :: max_model_number, i
1752 0 : real(dp) :: max_years_for_timestep
1753 : type (star_info), pointer :: s
1754 : logical :: all_same
1755 : integer, target :: ipar_ary(lipar)
1756 : integer, pointer :: ipar(:)
1757 0 : real(dp), target :: rpar_ary(lrpar)
1758 : real(dp), pointer :: rpar(:)
1759 0 : rpar => rpar_ary
1760 0 : ipar => ipar_ary
1761 : include 'formats'
1762 : ierr = 0
1763 0 : call get_star_ptr(id, s, ierr)
1764 0 : if (ierr /= 0) return
1765 :
1766 0 : ipar(1) = s% model_number
1767 0 : ipar(2) = min_steps
1768 :
1769 0 : rpar(1) = new_irrad_flux
1770 0 : rpar(2) = new_irrad_col_depth
1771 0 : rpar(3) = s% irradiation_flux
1772 0 : rpar(4) = s% column_depth_for_irradiation
1773 :
1774 0 : all_same = .true.
1775 0 : do i = 1, 2
1776 0 : if (abs(rpar(i)-rpar(i+2)) > 1d-10) then
1777 : all_same = .false.; exit
1778 : end if
1779 : end do
1780 0 : if (all_same) return
1781 :
1782 0 : write(*,'(A)')
1783 0 : write(*,2) 'relax to new irradiation -- min steps', min_steps
1784 0 : write(*,'(A)')
1785 :
1786 0 : max_model_number = s% max_model_number
1787 0 : s% max_model_number = -1111
1788 0 : max_years_for_timestep = s% max_years_for_timestep
1789 0 : s% max_years_for_timestep = relax_irradiation_max_yrs_dt
1790 : call do_internal_evolve( &
1791 : id, before_evolve_relax_irradiation, &
1792 : relax_irradiation_adjust_model, relax_irradiation_check_model, &
1793 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1794 0 : s% max_model_number = max_model_number
1795 0 : s% max_years_for_timestep = max_years_for_timestep
1796 0 : call error_check('relax irradiation',ierr)
1797 0 : end subroutine do_relax_irradiation
1798 :
1799 :
1800 0 : subroutine before_evolve_relax_irradiation(s, id, lipar, ipar, lrpar, rpar, ierr)
1801 : type (star_info), pointer :: s
1802 : integer, intent(in) :: id, lipar, lrpar
1803 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1804 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1805 : integer, intent(out) :: ierr
1806 0 : ierr = 0
1807 0 : call turn_off_winds(s)
1808 0 : s% max_model_number = -111
1809 0 : end subroutine before_evolve_relax_irradiation
1810 :
1811 0 : integer function relax_irradiation_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1812 : type (star_info), pointer :: s
1813 : integer, intent(in) :: id, lipar, lrpar
1814 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1815 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1816 0 : relax_irradiation_adjust_model = keep_going
1817 0 : end function relax_irradiation_adjust_model
1818 :
1819 0 : integer function relax_irradiation_check_model(s, id, lipar, ipar, lrpar, rpar)
1820 : use do_one_utils, only:do_bare_bones_check_model
1821 : type (star_info), pointer :: s
1822 : integer, intent(in) :: id, lipar, lrpar
1823 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1824 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1825 :
1826 : integer :: adjust_model, max_num_steps, num_steps
1827 0 : real(dp) :: old_irrad_flux, old_irrad_col_depth
1828 0 : real(dp) :: new_irrad_flux, new_irrad_col_depth, frac
1829 : logical, parameter :: dbg = .false.
1830 :
1831 : include 'formats'
1832 :
1833 0 : relax_irradiation_check_model = do_bare_bones_check_model(id)
1834 0 : if (relax_irradiation_check_model /= keep_going) return
1835 :
1836 0 : adjust_model = ipar(1)
1837 0 : max_num_steps = ipar(2)
1838 :
1839 0 : new_irrad_flux = rpar(1)
1840 0 : new_irrad_col_depth = rpar(2)
1841 0 : old_irrad_flux = rpar(3)
1842 0 : old_irrad_col_depth = rpar(4)
1843 0 : num_steps = s% model_number - adjust_model
1844 0 : frac = dble(num_steps)/dble(max_num_steps)
1845 :
1846 0 : if (s% dt < s% max_years_for_timestep*secyer) then
1847 0 : ipar(1) = adjust_model + 1
1848 0 : write(*,'(a60,2i6,3x,99e12.3)') 'relax irradiation, model, step, frac, flux, wait for dt', &
1849 0 : s% model_number, num_steps, frac, s% irradiation_flux
1850 0 : return
1851 : end if
1852 :
1853 0 : if (frac >= 1d0) then
1854 0 : s% irradiation_flux = new_irrad_flux
1855 0 : s% column_depth_for_irradiation = new_irrad_col_depth
1856 0 : relax_irradiation_check_model = terminate
1857 0 : s% termination_code = t_relax_finished_okay
1858 : write(*,'(a60,2i6,3x,99e12.3)') &
1859 0 : 'DONE: relax irradiation, model, step, fraction done, flux', &
1860 0 : s% model_number, num_steps, frac, s% irradiation_flux
1861 0 : return
1862 : end if
1863 :
1864 : s% irradiation_flux = &
1865 0 : frac*new_irrad_flux + (1-frac)*old_irrad_flux
1866 : s% column_depth_for_irradiation = &
1867 0 : frac*new_irrad_col_depth + (1-frac)*old_irrad_col_depth
1868 :
1869 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1870 0 : write(*,'(a60,2i6,3x,99e12.3)') 'relax irradiation, model, step, fraction done, flux', &
1871 0 : s% model_number, num_steps, frac, s% irradiation_flux
1872 :
1873 0 : end function relax_irradiation_check_model
1874 :
1875 :
1876 0 : subroutine do_relax_mass_change( &
1877 : id, min_steps, initial_mass_change, final_mass_change, relax_mass_change_max_yrs_dt, ierr)
1878 : integer, intent(in) :: id, min_steps
1879 : real(dp), intent(in) :: initial_mass_change, final_mass_change, relax_mass_change_max_yrs_dt
1880 : integer, intent(out) :: ierr
1881 :
1882 : integer, parameter :: lipar=2, lrpar=2
1883 : integer :: max_model_number
1884 0 : real(dp) :: max_years_for_timestep
1885 : type (star_info), pointer :: s
1886 : integer, target :: ipar_ary(lipar)
1887 : integer, pointer :: ipar(:)
1888 0 : real(dp), target :: rpar_ary(lrpar)
1889 : real(dp), pointer :: rpar(:)
1890 :
1891 0 : rpar => rpar_ary
1892 0 : ipar => ipar_ary
1893 : include 'formats'
1894 :
1895 0 : ierr = 0
1896 0 : if (abs(initial_mass_change - final_mass_change) < 1d-10) return
1897 :
1898 0 : call get_star_ptr(id, s, ierr)
1899 0 : if (ierr /= 0) return
1900 :
1901 :
1902 0 : ipar(1) = s% model_number
1903 0 : ipar(2) = min_steps
1904 0 : rpar(1) = initial_mass_change
1905 0 : rpar(2) = final_mass_change
1906 :
1907 0 : write(*,'(A)')
1908 0 : write(*,2) 'relax_mass_change -- min steps, init, final, max dt', &
1909 0 : min_steps, initial_mass_change, final_mass_change, relax_mass_change_max_yrs_dt
1910 0 : write(*,'(A)')
1911 :
1912 0 : max_model_number = s% max_model_number
1913 0 : s% max_model_number = -1111
1914 0 : max_years_for_timestep = s% max_years_for_timestep
1915 0 : s% max_years_for_timestep = relax_mass_change_max_yrs_dt
1916 :
1917 : call do_internal_evolve( &
1918 : id, before_evolve_relax_mass_change, &
1919 : relax_mass_change_adjust_model, relax_mass_change_check_model, &
1920 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
1921 :
1922 0 : s% max_model_number = max_model_number
1923 0 : s% max_years_for_timestep = max_years_for_timestep
1924 :
1925 0 : if (ierr == 0) s% mass_change = final_mass_change
1926 :
1927 0 : call error_check('relax mass change',ierr)
1928 :
1929 0 : end subroutine do_relax_mass_change
1930 :
1931 :
1932 0 : subroutine before_evolve_relax_mass_change(s, id, lipar, ipar, lrpar, rpar, ierr)
1933 : type (star_info), pointer :: s
1934 : integer, intent(in) :: id, lipar, lrpar
1935 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1936 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1937 : integer, intent(out) :: ierr
1938 0 : ierr = 0
1939 0 : call turn_off_winds(s)
1940 0 : s% max_model_number = -111
1941 0 : end subroutine before_evolve_relax_mass_change
1942 :
1943 0 : integer function relax_mass_change_adjust_model(s, id, lipar, ipar, lrpar, rpar)
1944 : type (star_info), pointer :: s
1945 : integer, intent(in) :: id, lipar, lrpar
1946 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1947 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1948 0 : relax_mass_change_adjust_model = keep_going
1949 0 : end function relax_mass_change_adjust_model
1950 :
1951 0 : integer function relax_mass_change_check_model(s, id, lipar, ipar, lrpar, rpar)
1952 : use do_one_utils, only:do_bare_bones_check_model
1953 : type (star_info), pointer :: s
1954 : integer, intent(in) :: id, lipar, lrpar
1955 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
1956 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
1957 :
1958 : integer :: adjust_model, max_num_steps, num_steps
1959 0 : real(dp) :: init_mass_change, final_mass_change, frac
1960 : logical, parameter :: dbg = .false.
1961 :
1962 : include 'formats'
1963 :
1964 0 : relax_mass_change_check_model = do_bare_bones_check_model(id)
1965 0 : if (relax_mass_change_check_model /= keep_going) return
1966 :
1967 0 : adjust_model = ipar(1)
1968 0 : max_num_steps = ipar(2)
1969 0 : num_steps = s% model_number - adjust_model
1970 :
1971 0 : init_mass_change = rpar(1)
1972 0 : final_mass_change = rpar(2)
1973 0 : frac = dble(num_steps)/dble(max_num_steps)
1974 :
1975 0 : if (s% dt < s% max_years_for_timestep*secyer) then
1976 0 : ipar(1) = adjust_model + 1 ! don't count this one
1977 0 : write(*,'(a60,2i6,3x,99e12.3)') 'relax_mass_change wait for dt: model, step, frac', &
1978 0 : s% model_number, num_steps, frac
1979 0 : return
1980 : end if
1981 :
1982 0 : if (frac >= 1d0) then
1983 0 : s% mass_change = final_mass_change
1984 0 : relax_mass_change_check_model = terminate
1985 0 : s% termination_code = t_relax_finished_okay
1986 0 : write(*,'(a60,2i6,3x,99e12.3)') 'DONE: relax_mass_change'
1987 0 : return
1988 : end if
1989 :
1990 0 : s% mass_change = frac*final_mass_change + (1-frac)*init_mass_change
1991 :
1992 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
1993 0 : write(*,'(a60,2i6,3x,99e12.3)') 'relax_mass_change, model, step, fraction done, mass_change', &
1994 0 : s% model_number, num_steps, frac, s% mass_change
1995 :
1996 0 : end function relax_mass_change_check_model
1997 :
1998 :
1999 0 : subroutine do_relax_core( &
2000 : id, new_core_mass, dlg_core_mass_per_step, &
2001 : relax_core_years_for_dt, core_avg_rho, core_avg_eps, ierr)
2002 : integer, intent(in) :: id
2003 : real(dp), intent(in) :: new_core_mass ! in Msun units
2004 : real(dp), intent(in) :: dlg_core_mass_per_step, relax_core_years_for_dt
2005 : real(dp), intent(in) :: core_avg_rho, core_avg_eps
2006 : ! adjust R_center according to core_avg_rho (g cm^-3)
2007 : ! adjust L_center according to core_avg_eps (erg g^-1 s^-1)
2008 : integer, intent(out) :: ierr
2009 :
2010 : integer, parameter :: lipar=1, lrpar=5
2011 : integer :: max_model_number
2012 0 : real(dp) :: max_years_for_timestep
2013 0 : real(dp) :: current_core_mass
2014 : type (star_info), pointer :: s
2015 : integer, target :: ipar_ary(lipar)
2016 : integer, pointer :: ipar(:)
2017 0 : real(dp), target :: rpar_ary(lrpar)
2018 : real(dp), pointer :: rpar(:)
2019 0 : rpar => rpar_ary
2020 0 : ipar => ipar_ary
2021 : include 'formats'
2022 : ierr = 0
2023 0 : if (new_core_mass <= 0) then
2024 0 : ierr = -1
2025 0 : write(*,*) 'invalid new_core_mass', new_core_mass
2026 0 : return
2027 : end if
2028 0 : call get_star_ptr(id, s, ierr)
2029 0 : if (ierr /= 0) return
2030 0 : current_core_mass = s% M_center/Msun
2031 0 : if (abs(new_core_mass - current_core_mass) <= 1d-12*new_core_mass) then
2032 0 : call do1_relax_core(s, new_core_mass, core_avg_rho, core_avg_eps, ierr)
2033 0 : if (ierr /= 0) then
2034 0 : write(*,*) 'failed in do1_relax_core'
2035 : end if
2036 0 : return
2037 : end if
2038 0 : write(*,'(A)')
2039 0 : write(*,1) 'current core mass', current_core_mass
2040 0 : write(*,1) 'relax to new_core_mass', new_core_mass
2041 0 : write(*,'(A)')
2042 0 : write(*,1) 'dlg_core_mass_per_step', dlg_core_mass_per_step
2043 0 : write(*,'(A)')
2044 0 : rpar(1) = new_core_mass
2045 0 : rpar(2) = dlg_core_mass_per_step
2046 0 : rpar(3) = relax_core_years_for_dt
2047 0 : rpar(4) = core_avg_rho
2048 0 : rpar(5) = core_avg_eps
2049 0 : max_model_number = s% max_model_number
2050 0 : s% max_model_number = -1111
2051 0 : max_years_for_timestep = s% max_years_for_timestep
2052 0 : s% max_years_for_timestep = relax_core_years_for_dt
2053 : call do_internal_evolve( &
2054 : id, before_evolve_relax_core, relax_core_adjust_model, &
2055 : relax_core_check_model, null_finish_model, &
2056 0 : .true., lipar, ipar, lrpar, rpar, ierr)
2057 0 : s% max_model_number = max_model_number
2058 0 : s% max_years_for_timestep = max_years_for_timestep
2059 0 : call error_check('relax core',ierr)
2060 0 : end subroutine do_relax_core
2061 :
2062 :
2063 0 : subroutine before_evolve_relax_core(s, id, lipar, ipar, lrpar, rpar, ierr)
2064 : type (star_info), pointer :: s
2065 : integer, intent(in) :: id, lipar, lrpar
2066 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2067 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2068 : integer, intent(out) :: ierr
2069 0 : real(dp) :: relax_core_years_for_dt
2070 0 : ierr = 0
2071 0 : call setup_before_relax(s)
2072 0 : s% max_model_number = -111
2073 0 : relax_core_years_for_dt = rpar(3)
2074 0 : s% dt_next = secyer*relax_core_years_for_dt
2075 0 : end subroutine before_evolve_relax_core
2076 :
2077 0 : integer function relax_core_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2078 : type (star_info), pointer :: s
2079 : integer, intent(in) :: id, lipar, lrpar
2080 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2081 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2082 0 : relax_core_adjust_model = keep_going
2083 0 : end function relax_core_adjust_model
2084 :
2085 0 : integer function relax_core_check_model(s, id, lipar, ipar, lrpar, rpar)
2086 : use do_one_utils, only:do_bare_bones_check_model
2087 : type (star_info), pointer :: s
2088 : integer, intent(in) :: id, lipar, lrpar
2089 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2090 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2091 0 : real(dp) :: new_core_mass, dlg_core_mass_per_step, next
2092 0 : real(dp) :: relax_core_dt, relax_core_years_for_dt
2093 0 : real(dp) :: core_avg_rho, core_avg_eps, current_core_mass
2094 : integer :: ierr
2095 : logical :: end_now
2096 :
2097 : logical, parameter :: dbg = .false.
2098 :
2099 : include 'formats'
2100 :
2101 0 : relax_core_check_model = do_bare_bones_check_model(id)
2102 0 : if (relax_core_check_model /= keep_going) return
2103 :
2104 0 : new_core_mass = rpar(1)
2105 0 : dlg_core_mass_per_step = rpar(2)
2106 0 : relax_core_years_for_dt = rpar(3)
2107 0 : core_avg_rho = rpar(4)
2108 0 : core_avg_eps = rpar(5)
2109 0 : ierr = 0
2110 :
2111 0 : current_core_mass = s% M_center/Msun
2112 0 : if (abs(new_core_mass - current_core_mass) <= 1d-12*new_core_mass) then
2113 0 : call do1_relax_core(s, new_core_mass, core_avg_rho, core_avg_eps, ierr)
2114 0 : if (ierr /= 0) then
2115 0 : write(*,*) 'failed in do1_relax_core'
2116 : end if
2117 0 : relax_core_check_model = terminate
2118 0 : s% termination_code = t_relax_finished_okay
2119 0 : return
2120 : end if
2121 :
2122 0 : if (mod(s% model_number, s% terminal_interval) == 0) then
2123 0 : write(*,1) 'current & target core masses', &
2124 0 : current_core_mass, new_core_mass, &
2125 0 : (new_core_mass - current_core_mass)/new_core_mass
2126 : end if
2127 :
2128 0 : relax_core_dt = secyer*relax_core_years_for_dt
2129 :
2130 0 : if (s% dt < relax_core_dt*0.9d0) then
2131 0 : write(*,1) 's% dt < relax_core_dt*0.9d0', s% dt, relax_core_dt*0.9d0
2132 0 : write(*,1) 's% max_timestep', s% max_timestep
2133 0 : write(*,1) 's% max_years_for_timestep*secyer', s% max_years_for_timestep*secyer
2134 0 : write(*,'(A)')
2135 0 : return ! give a chance to stabilize
2136 : end if
2137 :
2138 0 : end_now=.false.
2139 0 : if (new_core_mass < current_core_mass) then
2140 0 : next = exp10(safe_log10(current_core_mass) - dlg_core_mass_per_step)
2141 0 : if (next < new_core_mass) then
2142 0 : next = new_core_mass
2143 0 : end_now = .true.
2144 : end if
2145 : else
2146 0 : next = exp10(log10(max(1d-8,current_core_mass)) + dlg_core_mass_per_step)
2147 0 : if (next > new_core_mass) then
2148 0 : next = new_core_mass
2149 0 : end_now = .true.
2150 : end if
2151 : end if
2152 :
2153 : if (dbg) write(*,1) 'next', next, log10(next)
2154 :
2155 0 : call do1_relax_core(s, next, core_avg_rho, core_avg_eps, ierr)
2156 0 : if (ierr /= 0) then
2157 0 : write(*,*) 'failed in do1_relax_core'
2158 0 : relax_core_check_model = terminate
2159 : end if
2160 :
2161 : if (.true.) then
2162 0 : write(*,1) 's% M_center', s% M_center
2163 0 : write(*,1) 's% L_center', s% L_center
2164 0 : write(*,1) 's% R_center', s% R_center
2165 0 : write(*,1) 's% xmstar', s% xmstar
2166 0 : write(*,'(A)')
2167 : end if
2168 :
2169 0 : if(end_now) then
2170 0 : relax_core_check_model = terminate
2171 0 : s% termination_code = t_relax_finished_okay
2172 0 : return
2173 : end if
2174 :
2175 0 : end function relax_core_check_model
2176 :
2177 :
2178 0 : subroutine do1_relax_core(s, next, core_avg_rho, core_avg_eps, ierr)
2179 : type (star_info), pointer :: s
2180 : real(dp), intent(in) :: next, core_avg_rho, core_avg_eps
2181 : integer, intent(out) :: ierr
2182 0 : real(dp) :: next_M_center, next_R_center, next_L_center
2183 0 : ierr = 0
2184 0 : next_M_center = next*Msun
2185 0 : s% M_center = next_M_center
2186 0 : s% xmstar = s% mstar - s% M_center
2187 0 : next_R_center = pow(s% M_center/(core_avg_rho*four_thirds_pi),one_third)
2188 0 : call do1_relax_R_center(s, next_R_center, ierr)
2189 0 : if (ierr /= 0) return
2190 0 : next_L_center = s% M_center*core_avg_eps
2191 0 : call do1_relax_L_center(s, next_L_center, ierr)
2192 0 : end subroutine do1_relax_core
2193 :
2194 :
2195 0 : subroutine do_relax_mass_scale( &
2196 : id, new_mass, dlgm_per_step, change_mass_years_for_dt, ierr)
2197 : integer, intent(in) :: id
2198 : real(dp), intent(in) :: new_mass, dlgm_per_step, change_mass_years_for_dt
2199 : integer, intent(out) :: ierr
2200 : integer, parameter :: lipar=1, lrpar=3
2201 : integer :: max_model_number
2202 0 : real(dp) :: max_years_for_timestep, relax_mass_scale_dt, eps_mdot_factor
2203 : logical :: adding_mass
2204 : type (star_info), pointer :: s
2205 : integer, target :: ipar_ary(lipar)
2206 : integer, pointer :: ipar(:)
2207 0 : real(dp), target :: rpar_ary(lrpar)
2208 : real(dp), pointer :: rpar(:)
2209 0 : rpar => rpar_ary
2210 0 : ipar => ipar_ary
2211 : include 'formats'
2212 : ierr = 0
2213 0 : if (new_mass <= 0) then
2214 0 : ierr = -1
2215 0 : write(*,*) 'invalid new_mass', new_mass
2216 0 : return
2217 : end if
2218 0 : call get_star_ptr(id, s, ierr)
2219 0 : if (ierr /= 0) return
2220 0 : if (abs(new_mass - s% star_mass) <= 1d-12*new_mass) then
2221 0 : s% star_mass = new_mass
2222 0 : s% mstar = new_mass*Msun
2223 0 : s% xmstar = s% mstar - s% M_center
2224 0 : return
2225 : end if
2226 0 : write(*,'(A)')
2227 0 : write(*,1) 'relax_mass_scale'
2228 0 : write(*,1) 'current star_mass', s% star_mass
2229 0 : write(*,1) 'relax to new_mass', new_mass
2230 0 : write(*,1) 'dlgm_per_step', dlgm_per_step
2231 0 : write(*,'(A)')
2232 0 : rpar(1) = new_mass
2233 0 : rpar(2) = dlgm_per_step
2234 0 : rpar(3) = change_mass_years_for_dt
2235 0 : max_model_number = s% max_model_number
2236 0 : s% max_model_number = -1111
2237 :
2238 0 : adding_mass = (new_mass > s% star_mass)
2239 :
2240 0 : eps_mdot_factor = s% eps_mdot_factor
2241 0 : s% eps_mdot_factor = 0
2242 :
2243 0 : max_years_for_timestep = s% max_years_for_timestep
2244 0 : relax_mass_scale_dt = secyer*change_mass_years_for_dt
2245 0 : s% max_years_for_timestep = relax_mass_scale_dt/secyer
2246 : call do_internal_evolve( &
2247 : id, before_evolve_relax_mass_scale, &
2248 : relax_mass_scale_adjust_model, relax_mass_scale_check_model, null_finish_model, &
2249 0 : .true., lipar, ipar, lrpar, rpar, ierr)
2250 0 : s% max_model_number = max_model_number
2251 0 : s% star_mass = new_mass
2252 0 : s% mstar = new_mass*Msun
2253 0 : s% xmstar = s% mstar - s% M_center
2254 0 : s% eps_mdot_factor = eps_mdot_factor
2255 0 : s% max_years_for_timestep = max_years_for_timestep
2256 :
2257 0 : call error_check('relax mass scale',ierr)
2258 : end subroutine do_relax_mass_scale
2259 :
2260 :
2261 0 : subroutine before_evolve_relax_mass_scale(s, id, lipar, ipar, lrpar, rpar, ierr)
2262 : type (star_info), pointer :: s
2263 : integer, intent(in) :: id, lipar, lrpar
2264 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2265 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2266 : integer, intent(out) :: ierr
2267 0 : real(dp) :: change_mass_years_for_dt
2268 0 : ierr = 0
2269 0 : call setup_before_relax(s)
2270 0 : s% max_model_number = -111
2271 0 : change_mass_years_for_dt = rpar(3)
2272 0 : s% dt_next = secyer*change_mass_years_for_dt
2273 0 : end subroutine before_evolve_relax_mass_scale
2274 :
2275 0 : integer function relax_mass_scale_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2276 : type (star_info), pointer :: s
2277 : integer, intent(in) :: id, lipar, lrpar
2278 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2279 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2280 0 : relax_mass_scale_adjust_model = keep_going
2281 0 : end function relax_mass_scale_adjust_model
2282 :
2283 0 : integer function relax_mass_scale_check_model(s, id, lipar, ipar, lrpar, rpar)
2284 : use do_one_utils, only:do_bare_bones_check_model
2285 : type (star_info), pointer :: s
2286 : integer, intent(in) :: id, lipar, lrpar
2287 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2288 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2289 0 : real(dp) :: new_mass, dlgm_per_step, next
2290 0 : real(dp) :: relax_mass_scale_dt, change_mass_years_for_dt
2291 : logical, parameter :: dbg = .false.
2292 :
2293 : include 'formats'
2294 :
2295 0 : relax_mass_scale_check_model = do_bare_bones_check_model(id)
2296 0 : if (relax_mass_scale_check_model /= keep_going) return
2297 :
2298 0 : new_mass = rpar(1)
2299 0 : dlgm_per_step = rpar(2)
2300 0 : change_mass_years_for_dt = rpar(3)
2301 0 : if (s% star_mass < 0.01d0) dlgm_per_step = dlgm_per_step*0.1d0
2302 : !if (s% star_mass < 0.001d0) dlgm_per_step = dlgm_per_step*0.1d0
2303 :
2304 : if (dbg) then
2305 : write(*,1) 'new_mass', new_mass
2306 : write(*,1) 'current mass', s% star_mass
2307 : end if
2308 :
2309 0 : if (abs(s% star_mass-new_mass) < 1d-12*new_mass) then
2310 0 : s% star_mass = new_mass
2311 0 : s% mstar = new_mass*Msun
2312 0 : s% xmstar = s% mstar - s% M_center
2313 0 : s% termination_code = t_relax_finished_okay
2314 0 : relax_mass_scale_check_model = terminate
2315 0 : return
2316 : end if
2317 :
2318 0 : relax_mass_scale_dt = secyer*change_mass_years_for_dt
2319 :
2320 0 : if (s% dt < relax_mass_scale_dt*0.9d0) return ! give a chance to stabilize
2321 :
2322 0 : if (new_mass < s% star_mass) then
2323 0 : next = exp10(safe_log10(s% star_mass) - dlgm_per_step)
2324 0 : if (next < new_mass) next = new_mass
2325 : else
2326 0 : next = exp10(safe_log10(s% star_mass) + dlgm_per_step)
2327 0 : if (next > new_mass) next = new_mass
2328 : end if
2329 :
2330 : if (dbg) write(*,1) 'next', next, log10(next)
2331 :
2332 0 : s% star_mass = next
2333 0 : s% mstar = next*Msun
2334 0 : s% xmstar = s% mstar - s% M_center
2335 :
2336 0 : end function relax_mass_scale_check_model
2337 :
2338 :
2339 0 : subroutine do_relax_M_center(id, new_mass, dlgm_per_step, relax_M_center_dt, ierr)
2340 : integer, intent(in) :: id
2341 : real(dp), intent(in) :: new_mass, dlgm_per_step, relax_M_center_dt
2342 : integer, intent(out) :: ierr
2343 : integer, parameter :: lipar=1, lrpar=3
2344 : integer :: max_model_number
2345 0 : real(dp) :: max_years_for_timestep
2346 : type (star_info), pointer :: s
2347 : integer, target :: ipar_ary(lipar)
2348 : integer, pointer :: ipar(:)
2349 0 : real(dp), target :: rpar_ary(lrpar)
2350 : real(dp), pointer :: rpar(:)
2351 0 : rpar => rpar_ary
2352 0 : ipar => ipar_ary
2353 : include 'formats'
2354 : ierr = 0
2355 0 : if (new_mass <= 0) then
2356 0 : ierr = -1
2357 0 : write(*,*) 'invalid new_mass', new_mass
2358 0 : return
2359 : end if
2360 0 : call get_star_ptr(id, s, ierr)
2361 0 : if (ierr /= 0) return
2362 0 : if (abs(new_mass - s% star_mass) <= 1d-6) then
2363 0 : s% star_mass = new_mass
2364 0 : s% mstar = new_mass*Msun
2365 0 : s% xmstar = s% mstar - s% M_center
2366 0 : return
2367 : end if
2368 0 : write(*,'(A)')
2369 0 : write(*,1) 'current star_mass', s% star_mass
2370 0 : write(*,1) 'relax to new_mass', new_mass
2371 0 : write(*,'(A)')
2372 0 : write(*,1) 'dlgm_per_step', dlgm_per_step
2373 0 : write(*,'(A)')
2374 0 : rpar(1) = new_mass
2375 0 : rpar(2) = dlgm_per_step
2376 0 : rpar(3) = relax_M_center_dt
2377 0 : max_model_number = s% max_model_number
2378 0 : max_years_for_timestep = s% max_years_for_timestep
2379 0 : s% max_model_number = -1111
2380 0 : s% max_years_for_timestep = relax_M_center_dt/secyer
2381 : call do_internal_evolve( &
2382 : id, before_evolve_relax_M_center, &
2383 : relax_M_center_adjust_model, relax_M_center_check_model, &
2384 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
2385 :
2386 0 : s% max_model_number = max_model_number
2387 0 : s% max_years_for_timestep = max_years_for_timestep
2388 0 : call error_check('relax M center',ierr)
2389 0 : end subroutine do_relax_M_center
2390 :
2391 :
2392 0 : subroutine before_evolve_relax_M_center(s, id, lipar, ipar, lrpar, rpar, ierr)
2393 : type (star_info), pointer :: s
2394 : integer, intent(in) :: id, lipar, lrpar
2395 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2396 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2397 : integer, intent(out) :: ierr
2398 0 : ierr = 0
2399 0 : call setup_before_relax(s)
2400 0 : s% max_model_number = -111
2401 0 : s% dt_next = rpar(3) ! relax_M_center_dt
2402 0 : end subroutine before_evolve_relax_M_center
2403 :
2404 0 : integer function relax_M_center_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2405 : type (star_info), pointer :: s
2406 : integer, intent(in) :: id, lipar, lrpar
2407 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2408 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2409 0 : relax_M_center_adjust_model = keep_going
2410 0 : end function relax_M_center_adjust_model
2411 :
2412 0 : integer function relax_M_center_check_model(s, id, lipar, ipar, lrpar, rpar)
2413 : use do_one_utils, only:do_bare_bones_check_model
2414 : type (star_info), pointer :: s
2415 : integer, intent(in) :: id, lipar, lrpar
2416 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2417 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2418 : integer :: ierr
2419 0 : real(dp) :: new_mass, dlgm_per_step, relax_M_center_dt, next
2420 : logical, parameter :: dbg = .false.
2421 : logical :: end_now
2422 :
2423 : include 'formats'
2424 :
2425 0 : relax_M_center_check_model = do_bare_bones_check_model(id)
2426 0 : if (relax_M_center_check_model /= keep_going) return
2427 :
2428 0 : new_mass = rpar(1)
2429 0 : dlgm_per_step = rpar(2)
2430 0 : relax_M_center_dt = rpar(3)
2431 :
2432 0 : if (mod(s% model_number, s% terminal_interval) == 0 .and. s% M_center>0.0d0) &
2433 0 : write(*,1) 'relax_M_center target/current', new_mass/(s% M_center/Msun)
2434 :
2435 0 : end_now=.false.
2436 0 : if (new_mass < s% star_mass) then
2437 0 : next = exp10(safe_log10(s% star_mass) - dlgm_per_step)
2438 0 : if (next < new_mass) then
2439 0 : next = new_mass
2440 0 : end_now=.true.
2441 : end if
2442 : else
2443 0 : next = exp10(safe_log10(s% star_mass) + dlgm_per_step)
2444 0 : if (next > new_mass) then
2445 0 : next = new_mass
2446 0 : end_now=.true.
2447 : end if
2448 : end if
2449 :
2450 0 : if (abs(s% star_mass - new_mass) < 1d-15 .or. end_now) then
2451 0 : call set_new_mass_for_relax_M_center(s, new_mass, ierr)
2452 0 : if (ierr /= 0) then
2453 0 : write(*,*) 'failed to set mass for relax_M_center'
2454 0 : relax_M_center_check_model = terminate
2455 0 : return
2456 : end if
2457 0 : write(*,1) 'final mass', s% star_mass, s% mstar, s% M_center, s% xmstar
2458 0 : relax_M_center_check_model = terminate
2459 0 : s% termination_code = t_relax_finished_okay
2460 0 : return
2461 : end if
2462 :
2463 0 : if (s% dt < relax_M_center_dt*0.9d0) return ! give a chance to stabilize
2464 :
2465 : if (dbg) write(*,1) 'next', next, log10(next)
2466 :
2467 0 : call set_new_mass_for_relax_M_center(s, next, ierr)
2468 0 : if (ierr /= 0) then
2469 0 : write(*,*) 'failed to set mass for relax_M_center'
2470 0 : relax_M_center_check_model = terminate
2471 0 : return
2472 : end if
2473 :
2474 0 : end function relax_M_center_check_model
2475 :
2476 :
2477 0 : subroutine set_new_mass_for_relax_M_center(s, new_mass, ierr)
2478 0 : use star_utils, only: set_qs
2479 : type (star_info), pointer :: s
2480 : real(dp), intent(in) :: new_mass ! Msun
2481 : integer, intent(out) :: ierr
2482 : include 'formats'
2483 0 : ierr = 0
2484 0 : s% star_mass = new_mass
2485 0 : s% mstar = new_mass*Msun
2486 0 : s% M_center = s% mstar - s% xmstar
2487 0 : end subroutine set_new_mass_for_relax_M_center
2488 :
2489 :
2490 0 : subroutine do_relax_R_center(id, new_R_center, dlgR_per_step, relax_R_center_dt, ierr)
2491 : integer, intent(in) :: id
2492 : real(dp), intent(in) :: new_R_center, dlgR_per_step, relax_R_center_dt
2493 : integer, intent(out) :: ierr
2494 : integer, parameter :: lipar=1, lrpar=3
2495 : integer :: max_model_number
2496 0 : real(dp) :: max_years_for_timestep
2497 : type (star_info), pointer :: s
2498 : integer, target :: ipar_ary(lipar)
2499 : integer, pointer :: ipar(:)
2500 0 : real(dp), target :: rpar_ary(lrpar)
2501 : real(dp), pointer :: rpar(:)
2502 0 : rpar => rpar_ary
2503 0 : ipar => ipar_ary
2504 : include 'formats'
2505 : ierr = 0
2506 0 : if (new_R_center < 0) then
2507 0 : ierr = -1
2508 0 : write(*,*) 'invalid new_R_center', new_R_center
2509 0 : return
2510 : end if
2511 0 : call get_star_ptr(id, s, ierr)
2512 0 : if (ierr /= 0) return
2513 0 : if (abs(new_R_center - s% R_center) <= 1d-6) then
2514 0 : s% R_center = new_R_center
2515 0 : return
2516 : end if
2517 0 : write(*,'(A)')
2518 0 : write(*,1) 'current R_center', s% R_center
2519 0 : write(*,1) 'relax to new_R_center', new_R_center
2520 0 : write(*,'(A)')
2521 0 : write(*,1) 'dlgR_per_step', dlgR_per_step
2522 0 : write(*,'(A)')
2523 0 : rpar(1) = new_R_center
2524 0 : rpar(2) = dlgR_per_step
2525 0 : rpar(3) = relax_R_center_dt
2526 0 : max_model_number = s% max_model_number
2527 0 : max_years_for_timestep = s% max_years_for_timestep
2528 0 : s% max_model_number = -1111
2529 0 : s% max_years_for_timestep = relax_R_center_dt/secyer
2530 : call do_internal_evolve( &
2531 : id, before_evolve_relax_R_center, &
2532 : relax_R_center_adjust_model, relax_R_center_check_model, &
2533 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
2534 :
2535 0 : s% max_model_number = max_model_number
2536 0 : s% max_years_for_timestep = max_years_for_timestep
2537 0 : call error_check('relax R center',ierr)
2538 0 : end subroutine do_relax_R_center
2539 :
2540 :
2541 0 : subroutine before_evolve_relax_R_center(s, id, lipar, ipar, lrpar, rpar, ierr)
2542 : type (star_info), pointer :: s
2543 : integer, intent(in) :: id, lipar, lrpar
2544 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2545 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2546 : integer, intent(out) :: ierr
2547 0 : ierr = 0
2548 0 : call setup_before_relax(s)
2549 0 : s% max_model_number = -111
2550 0 : s% dt_next = rpar(3) ! relax_R_center_dt
2551 0 : end subroutine before_evolve_relax_R_center
2552 :
2553 0 : integer function relax_R_center_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2554 : type (star_info), pointer :: s
2555 : integer, intent(in) :: id, lipar, lrpar
2556 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2557 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2558 0 : relax_R_center_adjust_model = keep_going
2559 0 : end function relax_R_center_adjust_model
2560 :
2561 0 : integer function relax_R_center_check_model(s, id, lipar, ipar, lrpar, rpar)
2562 : use do_one_utils, only:do_bare_bones_check_model
2563 : type (star_info), pointer :: s
2564 : integer, intent(in) :: id, lipar, lrpar
2565 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2566 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2567 : integer :: ierr
2568 0 : real(dp) :: new_R_center, dlgR_per_step, relax_R_center_dt, next
2569 : logical, parameter :: dbg = .false.
2570 :
2571 : include 'formats'
2572 :
2573 0 : relax_R_center_check_model = do_bare_bones_check_model(id)
2574 0 : if (relax_R_center_check_model /= keep_going) return
2575 :
2576 0 : new_R_center = rpar(1)
2577 0 : dlgR_per_step = rpar(2)
2578 0 : relax_R_center_dt = rpar(3)
2579 :
2580 0 : if (mod(s% model_number, s% terminal_interval) == 0 .and. s% R_center > 0d0) &
2581 0 : write(*,1) 'relax_R_center target/current', new_R_center/s% R_center
2582 :
2583 0 : if (abs(s% R_center - new_R_center) < 1d-15) then
2584 0 : call do1_relax_R_center(s, new_R_center, ierr)
2585 0 : if (ierr /= 0) then
2586 0 : write(*,*) 'failed in relax_R_center'
2587 : end if
2588 0 : relax_R_center_check_model = terminate
2589 0 : s% termination_code = t_relax_finished_okay
2590 0 : return
2591 : end if
2592 :
2593 0 : if (s% dt < relax_R_center_dt*0.9d0) return ! give a chance to stabilize
2594 :
2595 0 : if (new_R_center < s% R_center) then
2596 0 : next = exp10(safe_log10(s% R_center) - dlgR_per_step)
2597 0 : if (next < new_R_center) next = new_R_center
2598 0 : else if (s% R_center < 1) then
2599 0 : next = 1
2600 : else
2601 0 : next = exp10(safe_log10(s% R_center) + dlgR_per_step)
2602 0 : if (next > new_R_center) next = new_R_center
2603 : end if
2604 :
2605 : if (dbg) write(*,1) 'next', next
2606 :
2607 0 : call do1_relax_R_center(s, next, ierr)
2608 0 : if (ierr /= 0) then
2609 0 : write(*,*) 'failed in relax_R_center'
2610 0 : relax_R_center_check_model = terminate
2611 : end if
2612 :
2613 0 : end function relax_R_center_check_model
2614 :
2615 :
2616 0 : subroutine do1_relax_R_center(s, new_Rcenter, ierr)
2617 : ! adjust all lnR's to keep same density for each cell as 1st guess for next model
2618 0 : use star_utils, only: set_qs, store_lnR_in_xh
2619 : type (star_info), pointer :: s
2620 : real(dp), intent(in) :: new_Rcenter ! cm
2621 : integer, intent(out) :: ierr
2622 0 : real(dp) :: dm, rho, dr3, rp13
2623 : integer :: k
2624 : include 'formats'
2625 0 : ierr = 0
2626 0 : s% R_center = new_Rcenter
2627 : ! adjust lnR's
2628 0 : rp13 = s% R_center*s% R_center*s% R_center
2629 0 : do k = s% nz, 1, -1
2630 0 : dm = s% dm(k)
2631 0 : rho = s% rho(k)
2632 0 : dr3 = dm/(rho*four_thirds_pi) ! dm/rho is cell volume
2633 0 : call store_lnR_in_xh(s, k, log(rp13 + dr3)*one_third)
2634 0 : rp13 = rp13 + dr3
2635 : end do
2636 0 : end subroutine do1_relax_R_center
2637 :
2638 :
2639 0 : subroutine do_relax_v_center(id, new_v_center, dv_per_step, relax_v_center_dt, ierr)
2640 : integer, intent(in) :: id
2641 : real(dp), intent(in) :: new_v_center, dv_per_step, relax_v_center_dt
2642 : integer, intent(out) :: ierr
2643 : integer, parameter :: lipar=1, lrpar=3
2644 : integer :: max_model_number
2645 0 : real(dp) :: max_years_for_timestep
2646 : type (star_info), pointer :: s
2647 : integer, target :: ipar_ary(lipar)
2648 : integer, pointer :: ipar(:)
2649 0 : real(dp), target :: rpar_ary(lrpar)
2650 : real(dp), pointer :: rpar(:)
2651 0 : rpar => rpar_ary
2652 0 : ipar => ipar_ary
2653 : include 'formats'
2654 : ierr = 0
2655 0 : if (new_v_center < 0) then
2656 0 : ierr = -1
2657 0 : write(*,*) 'invalid new_v_center', new_v_center
2658 0 : return
2659 : end if
2660 0 : call get_star_ptr(id, s, ierr)
2661 0 : if (ierr /= 0) return
2662 0 : if (abs(s% v_center - new_v_center) < &
2663 : 1d-6*max(1d-6,abs(s% v_center),abs(new_v_center))) then
2664 0 : s% v_center = new_v_center
2665 0 : return
2666 : end if
2667 0 : write(*,'(A)')
2668 0 : write(*,1) 'current v_center', s% v_center
2669 0 : write(*,1) 'relax to new_v_center', new_v_center
2670 0 : write(*,'(A)')
2671 0 : write(*,1) 'dv_per_step', dv_per_step
2672 0 : write(*,'(A)')
2673 0 : rpar(1) = new_v_center
2674 0 : rpar(2) = dv_per_step
2675 0 : rpar(3) = relax_v_center_dt
2676 0 : max_model_number = s% max_model_number
2677 0 : max_years_for_timestep = s% max_years_for_timestep
2678 0 : s% max_model_number = -1111
2679 0 : s% max_years_for_timestep = relax_v_center_dt/secyer
2680 : call do_internal_evolve( &
2681 : id, before_evolve_relax_v_center, &
2682 : relax_v_center_adjust_model, relax_v_center_check_model, &
2683 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
2684 :
2685 0 : s% max_model_number = max_model_number
2686 0 : s% max_years_for_timestep = max_years_for_timestep
2687 0 : call error_check('relax V center',ierr)
2688 0 : end subroutine do_relax_v_center
2689 :
2690 :
2691 0 : subroutine before_evolve_relax_v_center(s, id, lipar, ipar, lrpar, rpar, ierr)
2692 : type (star_info), pointer :: s
2693 : integer, intent(in) :: id, lipar, lrpar
2694 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2695 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2696 : integer, intent(out) :: ierr
2697 0 : ierr = 0
2698 0 : call setup_before_relax(s)
2699 0 : s% max_model_number = -111
2700 0 : s% dt_next = rpar(3) ! relax_v_center_dt
2701 0 : end subroutine before_evolve_relax_v_center
2702 :
2703 0 : integer function relax_v_center_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2704 : type (star_info), pointer :: s
2705 : integer, intent(in) :: id, lipar, lrpar
2706 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2707 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2708 0 : relax_v_center_adjust_model = keep_going
2709 0 : end function relax_v_center_adjust_model
2710 :
2711 0 : integer function relax_v_center_check_model(s, id, lipar, ipar, lrpar, rpar)
2712 : use do_one_utils, only:do_bare_bones_check_model
2713 : type (star_info), pointer :: s
2714 : integer, intent(in) :: id, lipar, lrpar
2715 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2716 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2717 0 : real(dp) :: new_v_center, dv_per_step, relax_v_center_dt, next
2718 : logical, parameter :: dbg = .false.
2719 :
2720 : include 'formats'
2721 :
2722 0 : relax_v_center_check_model = do_bare_bones_check_model(id)
2723 0 : if (relax_v_center_check_model /= keep_going) return
2724 :
2725 0 : new_v_center = rpar(1)
2726 0 : dv_per_step = rpar(2)
2727 0 : relax_v_center_dt = rpar(3)
2728 :
2729 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
2730 0 : write(*,1) 'target v_center current', new_v_center, s% v_center
2731 :
2732 0 : if (abs(s% v_center - new_v_center) < &
2733 : 1d-6*max(1d-6,abs(s% v_center),abs(new_v_center))) then
2734 0 : s% v_center = new_v_center
2735 0 : relax_v_center_check_model = terminate
2736 0 : s% termination_code = t_relax_finished_okay
2737 0 : return
2738 : end if
2739 :
2740 0 : if (s% dt < relax_v_center_dt*0.9d0) return ! give a chance to stabilize
2741 :
2742 0 : if (new_v_center < s% v_center) then
2743 0 : next = s% v_center - dv_per_step
2744 0 : if (next < new_v_center) next = new_v_center
2745 : else
2746 0 : next = s% v_center + dv_per_step
2747 0 : if (next > new_v_center) next = new_v_center
2748 : end if
2749 :
2750 : if (dbg) write(*,1) 'next', next
2751 :
2752 0 : s% v_center = next
2753 :
2754 0 : end function relax_v_center_check_model
2755 :
2756 :
2757 0 : subroutine do_relax_L_center(id, new_L_center, dlgL_per_step, relax_L_center_dt, ierr)
2758 : integer, intent(in) :: id
2759 : real(dp), intent(in) :: new_L_center, dlgL_per_step, relax_L_center_dt
2760 : integer, intent(out) :: ierr
2761 : integer, parameter :: lipar=1, lrpar=3
2762 : integer :: max_model_number
2763 0 : real(dp) :: max_years_for_timestep
2764 : type (star_info), pointer :: s
2765 : integer, target :: ipar_ary(lipar)
2766 : integer, pointer :: ipar(:)
2767 0 : real(dp), target :: rpar_ary(lrpar)
2768 : real(dp), pointer :: rpar(:)
2769 0 : rpar => rpar_ary
2770 0 : ipar => ipar_ary
2771 : include 'formats'
2772 : ierr = 0
2773 0 : call get_star_ptr(id, s, ierr)
2774 0 : if (ierr /= 0) return
2775 0 : if (abs(new_L_center - s% L_center) <= &
2776 : 1d-10*max(abs(new_L_center),abs(s% L_center),1d0)) then
2777 0 : s% L_center = new_L_center
2778 0 : return
2779 : end if
2780 0 : write(*,'(A)')
2781 0 : write(*,1) 'current L_center', s% L_center
2782 0 : write(*,1) 'relax to new_L_center', new_L_center
2783 0 : write(*,'(A)')
2784 0 : write(*,1) 'dlgL_per_step', dlgL_per_step
2785 0 : write(*,'(A)')
2786 0 : rpar(1) = new_L_center
2787 0 : rpar(2) = dlgL_per_step*(new_L_center - s% L_center)
2788 0 : rpar(3) = relax_L_center_dt
2789 0 : max_model_number = s% max_model_number
2790 0 : max_years_for_timestep = s% max_years_for_timestep
2791 0 : s% max_model_number = -1111
2792 0 : s% max_years_for_timestep = relax_L_center_dt/secyer
2793 : call do_internal_evolve( &
2794 : id, before_evolve_relax_L_center, &
2795 : relax_L_center_adjust_model, relax_L_center_check_model, &
2796 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
2797 :
2798 0 : s% max_model_number = max_model_number
2799 0 : s% max_years_for_timestep = max_years_for_timestep
2800 0 : call error_check('relax L center',ierr)
2801 0 : end subroutine do_relax_L_center
2802 :
2803 :
2804 0 : subroutine before_evolve_relax_L_center(s, id, lipar, ipar, lrpar, rpar, ierr)
2805 : type (star_info), pointer :: s
2806 : integer, intent(in) :: id, lipar, lrpar
2807 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2808 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2809 : integer, intent(out) :: ierr
2810 0 : ierr = 0
2811 0 : call setup_before_relax(s)
2812 0 : s% max_model_number = -111
2813 0 : s% dt_next = rpar(3) ! relax_L_center_dt
2814 0 : end subroutine before_evolve_relax_L_center
2815 :
2816 0 : integer function relax_L_center_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2817 : type (star_info), pointer :: s
2818 : integer, intent(in) :: id, lipar, lrpar
2819 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2820 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2821 0 : relax_L_center_adjust_model = keep_going
2822 0 : end function relax_L_center_adjust_model
2823 :
2824 0 : integer function relax_L_center_check_model(s, id, lipar, ipar, lrpar, rpar)
2825 : use do_one_utils, only:do_bare_bones_check_model
2826 : type (star_info), pointer :: s
2827 : integer, intent(in) :: id, lipar, lrpar
2828 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2829 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2830 : integer :: ierr
2831 0 : real(dp) :: new_L_center, dL, relax_L_center_dt, next
2832 : logical, parameter :: dbg = .false.
2833 :
2834 : include 'formats'
2835 :
2836 0 : relax_L_center_check_model = do_bare_bones_check_model(id)
2837 0 : if (relax_L_center_check_model /= keep_going) return
2838 :
2839 0 : new_L_center = rpar(1)
2840 0 : dL = rpar(2)
2841 0 : relax_L_center_dt = rpar(3)
2842 :
2843 :
2844 0 : if (mod(s% model_number, s% terminal_interval) == 0 .and. s% L_center>0.d0) &
2845 0 : write(*,1) 'relax_L_center target/current', new_L_center/s% L_center
2846 :
2847 0 : if (abs(new_L_center - s% L_center) < abs(dL)) then
2848 0 : call do1_relax_L_center(s, new_L_center, ierr)
2849 0 : if (ierr /= 0) then
2850 0 : write(*,*) 'failed in relax_L_center'
2851 : end if
2852 0 : relax_L_center_check_model = terminate
2853 0 : s% termination_code = t_relax_finished_okay
2854 0 : return
2855 : end if
2856 :
2857 0 : if (s% dt < relax_L_center_dt*0.9d0) return ! give a chance to stabilize
2858 :
2859 0 : next = s% L_center + dL
2860 : if (dbg) write(*,1) 'next', next
2861 :
2862 0 : call do1_relax_L_center(s, next, ierr)
2863 0 : if (ierr /= 0) then
2864 0 : write(*,*) 'failed in relax_L_center'
2865 0 : relax_L_center_check_model = terminate
2866 : end if
2867 :
2868 0 : end function relax_L_center_check_model
2869 :
2870 :
2871 0 : subroutine do1_relax_L_center(s, new_Lcenter, ierr)
2872 : type (star_info), pointer :: s
2873 : real(dp), intent(in) :: new_Lcenter
2874 : integer, intent(out) :: ierr
2875 0 : real(dp) :: L_center_prev, dL
2876 : integer :: i_lum, nz, k
2877 : include 'formats'
2878 0 : ierr = 0
2879 0 : nz = s% nz
2880 0 : L_center_prev = s% L_center
2881 0 : s% L_center = new_Lcenter
2882 0 : i_lum = s% i_lum
2883 0 : if (i_lum == 0) return
2884 0 : dL = new_Lcenter - L_center_prev
2885 0 : do k=1,nz
2886 0 : s% xh(i_lum,k) = s% xh(i_lum,k) + dL
2887 : end do
2888 0 : end subroutine do1_relax_L_center
2889 :
2890 :
2891 0 : subroutine do_relax_dxdt_nuc_factor(id, new_value, per_step_multiplier, ierr)
2892 : integer, intent(in) :: id
2893 : real(dp), intent(in) :: new_value, per_step_multiplier
2894 : integer, intent(out) :: ierr
2895 : integer, parameter :: lipar=1, lrpar=2
2896 : integer :: max_model_number
2897 : type (star_info), pointer :: s
2898 : integer, target :: ipar_ary(lipar)
2899 : integer, pointer :: ipar(:)
2900 0 : real(dp), target :: rpar_ary(lrpar)
2901 : real(dp), pointer :: rpar(:)
2902 0 : rpar => rpar_ary
2903 0 : ipar => ipar_ary
2904 : include 'formats'
2905 : ierr = 0
2906 0 : if (new_value <= 0) then
2907 0 : ierr = -1
2908 0 : write(*,*) 'invalid new_value', new_value
2909 0 : return
2910 : end if
2911 0 : if (per_step_multiplier <= 0 .or. per_step_multiplier == 1) then
2912 0 : ierr = -1
2913 0 : write(*,*) 'invalid per_step_multiplier', per_step_multiplier
2914 0 : return
2915 : end if
2916 0 : call get_star_ptr(id, s, ierr)
2917 0 : if (ierr /= 0) return
2918 0 : if (abs(new_value - s% dxdt_nuc_factor) <= 1d-6) then
2919 0 : s% dxdt_nuc_factor = new_value
2920 0 : return
2921 : end if
2922 0 : write(*,'(A)')
2923 0 : write(*,1) 'current dxdt_nuc_factor', s% dxdt_nuc_factor
2924 0 : write(*,1) 'relax to new_value', new_value
2925 0 : write(*,'(A)')
2926 0 : write(*,1) 'per_step_multiplier', per_step_multiplier
2927 0 : write(*,'(A)')
2928 0 : rpar(1) = new_value
2929 0 : rpar(2) = per_step_multiplier
2930 0 : max_model_number = s% max_model_number
2931 0 : s% max_model_number = -1111
2932 : call do_internal_evolve( &
2933 : id, before_evolve_relax_dxdt_nuc_factor, &
2934 : relax_dxdt_nuc_factor_adjust_model, relax_dxdt_nuc_factor_check_model, &
2935 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
2936 0 : s% max_model_number = max_model_number
2937 0 : s% dxdt_nuc_factor = new_value
2938 0 : call error_check('relax dxdt nuc factor',ierr)
2939 : end subroutine do_relax_dxdt_nuc_factor
2940 :
2941 :
2942 0 : subroutine before_evolve_relax_dxdt_nuc_factor(s, id, lipar, ipar, lrpar, rpar, ierr)
2943 : type (star_info), pointer :: s
2944 : integer, intent(in) :: id, lipar, lrpar
2945 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2946 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2947 : integer, intent(out) :: ierr
2948 0 : ierr = 0
2949 0 : call turn_off_winds(s)
2950 0 : s% max_model_number = -111
2951 0 : s% dt_next = secyer*1d-3
2952 0 : end subroutine before_evolve_relax_dxdt_nuc_factor
2953 :
2954 0 : integer function relax_dxdt_nuc_factor_adjust_model(s, id, lipar, ipar, lrpar, rpar)
2955 : type (star_info), pointer :: s
2956 : integer, intent(in) :: id, lipar, lrpar
2957 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2958 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2959 0 : relax_dxdt_nuc_factor_adjust_model = keep_going
2960 0 : end function relax_dxdt_nuc_factor_adjust_model
2961 :
2962 0 : integer function relax_dxdt_nuc_factor_check_model(s, id, lipar, ipar, lrpar, rpar)
2963 : use do_one_utils, only:do_bare_bones_check_model
2964 : type (star_info), pointer :: s
2965 : integer, intent(in) :: id, lipar, lrpar
2966 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
2967 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
2968 0 : real(dp) :: new_value, per_step_multiplier
2969 : logical, parameter :: dbg = .false.
2970 :
2971 : include 'formats'
2972 :
2973 0 : relax_dxdt_nuc_factor_check_model = do_bare_bones_check_model(id)
2974 0 : if (relax_dxdt_nuc_factor_check_model /= keep_going) return
2975 :
2976 0 : new_value = rpar(1)
2977 0 : per_step_multiplier = rpar(2)
2978 :
2979 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
2980 0 : write(*,1) 'new_value current', new_value, s% dxdt_nuc_factor
2981 :
2982 0 : s% dxdt_nuc_factor = s% dxdt_nuc_factor * per_step_multiplier
2983 :
2984 0 : if ((per_step_multiplier < 1 .and. s% dxdt_nuc_factor < new_value) .or. &
2985 : (per_step_multiplier > 1 .and. s% dxdt_nuc_factor > new_value)) then
2986 0 : s% dxdt_nuc_factor = new_value
2987 0 : relax_dxdt_nuc_factor_check_model = terminate
2988 0 : s% termination_code = t_relax_finished_okay
2989 0 : return
2990 : end if
2991 :
2992 0 : end function relax_dxdt_nuc_factor_check_model
2993 :
2994 :
2995 0 : subroutine do_relax_eps_nuc_factor(id, new_value, per_step_multiplier, ierr)
2996 : integer, intent(in) :: id
2997 : real(dp), intent(in) :: new_value, per_step_multiplier
2998 : integer, intent(out) :: ierr
2999 : integer, parameter :: lipar=1, lrpar=2
3000 : integer :: max_model_number
3001 : type (star_info), pointer :: s
3002 : integer, target :: ipar_ary(lipar)
3003 : integer, pointer :: ipar(:)
3004 0 : real(dp), target :: rpar_ary(lrpar)
3005 : real(dp), pointer :: rpar(:)
3006 0 : rpar => rpar_ary
3007 0 : ipar => ipar_ary
3008 : include 'formats'
3009 : ierr = 0
3010 0 : if (new_value <= 0) then
3011 0 : ierr = -1
3012 0 : write(*,*) 'invalid new_value', new_value
3013 0 : return
3014 : end if
3015 0 : if (per_step_multiplier <= 0 .or. per_step_multiplier == 1) then
3016 0 : ierr = -1
3017 0 : write(*,*) 'invalid per_step_multiplier', per_step_multiplier
3018 0 : return
3019 : end if
3020 0 : call get_star_ptr(id, s, ierr)
3021 0 : if (ierr /= 0) return
3022 0 : if (abs(new_value - s% eps_nuc_factor) <= 1d-6) then
3023 0 : s% eps_nuc_factor = new_value
3024 0 : return
3025 : end if
3026 0 : write(*,'(A)')
3027 0 : write(*,1) 'current eps_nuc_factor', s% eps_nuc_factor
3028 0 : write(*,1) 'relax to new_value', new_value
3029 0 : write(*,'(A)')
3030 0 : write(*,1) 'per_step_multiplier', per_step_multiplier
3031 0 : write(*,'(A)')
3032 0 : rpar(1) = new_value
3033 0 : rpar(2) = per_step_multiplier
3034 0 : max_model_number = s% max_model_number
3035 0 : s% max_model_number = -1111
3036 : call do_internal_evolve( &
3037 : id, before_evolve_relax_eps_nuc_factor, &
3038 : relax_eps_nuc_factor_adjust_model, relax_eps_nuc_factor_check_model, &
3039 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3040 0 : s% max_model_number = max_model_number
3041 0 : s% eps_nuc_factor = new_value
3042 0 : call error_check('relax eps nuc factor',ierr)
3043 0 : end subroutine do_relax_eps_nuc_factor
3044 :
3045 :
3046 0 : subroutine before_evolve_relax_eps_nuc_factor(s, id, lipar, ipar, lrpar, rpar, ierr)
3047 : type (star_info), pointer :: s
3048 : integer, intent(in) :: id, lipar, lrpar
3049 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3050 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3051 : integer, intent(out) :: ierr
3052 0 : ierr = 0
3053 0 : call turn_off_winds(s)
3054 0 : s% max_model_number = -111
3055 0 : s% dt_next = secyer*1d-3
3056 0 : end subroutine before_evolve_relax_eps_nuc_factor
3057 :
3058 0 : integer function relax_eps_nuc_factor_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3059 : type (star_info), pointer :: s
3060 : integer, intent(in) :: id, lipar, lrpar
3061 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3062 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3063 0 : relax_eps_nuc_factor_adjust_model = keep_going
3064 0 : end function relax_eps_nuc_factor_adjust_model
3065 :
3066 0 : integer function relax_eps_nuc_factor_check_model(s, id, lipar, ipar, lrpar, rpar)
3067 : use do_one_utils, only:do_bare_bones_check_model
3068 : type (star_info), pointer :: s
3069 : integer, intent(in) :: id, lipar, lrpar
3070 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3071 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3072 0 : real(dp) :: new_value, per_step_multiplier
3073 : logical, parameter :: dbg = .false.
3074 :
3075 : include 'formats'
3076 :
3077 0 : relax_eps_nuc_factor_check_model = do_bare_bones_check_model(id)
3078 0 : if (relax_eps_nuc_factor_check_model /= keep_going) return
3079 :
3080 0 : new_value = rpar(1)
3081 0 : per_step_multiplier = rpar(2)
3082 :
3083 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
3084 0 : write(*,1) 'new_value, current', new_value, s% eps_nuc_factor
3085 :
3086 0 : s% eps_nuc_factor = s% eps_nuc_factor * per_step_multiplier
3087 :
3088 0 : if ((per_step_multiplier < 1 .and. s% eps_nuc_factor < new_value) .or. &
3089 : (per_step_multiplier > 1 .and. s% eps_nuc_factor > new_value)) then
3090 0 : s% eps_nuc_factor = new_value
3091 0 : relax_eps_nuc_factor_check_model = terminate
3092 0 : s% termination_code = t_relax_finished_okay
3093 0 : return
3094 : end if
3095 :
3096 0 : end function relax_eps_nuc_factor_check_model
3097 :
3098 :
3099 0 : subroutine do_relax_opacity_max(id, new_value, per_step_multiplier, ierr)
3100 : integer, intent(in) :: id
3101 : real(dp), intent(in) :: new_value, per_step_multiplier
3102 : integer, intent(out) :: ierr
3103 : integer, parameter :: lipar=1, lrpar=2
3104 : integer :: max_model_number
3105 : type (star_info), pointer :: s
3106 : integer, target :: ipar_ary(lipar)
3107 : integer, pointer :: ipar(:)
3108 0 : real(dp), target :: rpar_ary(lrpar)
3109 : real(dp), pointer :: rpar(:)
3110 0 : rpar => rpar_ary
3111 0 : ipar => ipar_ary
3112 : include 'formats'
3113 : ierr = 0
3114 0 : if (new_value <= 0) then
3115 0 : ierr = -1
3116 0 : write(*,*) 'invalid new_value', new_value
3117 0 : return
3118 : end if
3119 0 : if (per_step_multiplier <= 0 .or. per_step_multiplier == 1) then
3120 0 : ierr = -1
3121 0 : write(*,*) 'invalid per_step_multiplier', per_step_multiplier
3122 0 : return
3123 : end if
3124 0 : call get_star_ptr(id, s, ierr)
3125 0 : if (ierr /= 0) return
3126 0 : if (s% opacity_max <= 0) then
3127 0 : ierr = -1
3128 0 : write(*,*) 'invalid opacity_max', s% opacity_max
3129 0 : return
3130 : end if
3131 0 : if (abs(new_value - s% opacity_max) <= 1d-6) then
3132 0 : s% opacity_max = new_value
3133 0 : return
3134 : end if
3135 0 : write(*,'(A)')
3136 0 : write(*,1) 'current opacity_max', s% opacity_max
3137 0 : write(*,1) 'relax to new_value', new_value
3138 0 : write(*,'(A)')
3139 0 : write(*,1) 'per_step_multiplier', per_step_multiplier
3140 0 : write(*,'(A)')
3141 0 : rpar(1) = new_value
3142 0 : rpar(2) = per_step_multiplier
3143 0 : max_model_number = s% max_model_number
3144 0 : s% max_model_number = -1111
3145 : call do_internal_evolve( &
3146 : id, before_evolve_relax_opacity_max, &
3147 : relax_opacity_max_adjust_model, relax_opacity_max_check_model, &
3148 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3149 :
3150 0 : s% max_model_number = max_model_number
3151 0 : s% opacity_max = new_value
3152 0 : s% dt_next = rpar(1) ! keep dt from relax
3153 0 : call error_check('relax opacity max',ierr)
3154 0 : end subroutine do_relax_opacity_max
3155 :
3156 :
3157 0 : subroutine before_evolve_relax_opacity_max(s, id, lipar, ipar, lrpar, rpar, ierr)
3158 : type (star_info), pointer :: s
3159 : integer, intent(in) :: id, lipar, lrpar
3160 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3161 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3162 : integer, intent(out) :: ierr
3163 0 : ierr = 0
3164 0 : s% max_model_number = -111
3165 0 : s% dt_next = secyer*1d-3
3166 0 : call turn_off_winds(s)
3167 0 : end subroutine before_evolve_relax_opacity_max
3168 :
3169 0 : integer function relax_opacity_max_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3170 : type (star_info), pointer :: s
3171 : integer, intent(in) :: id, lipar, lrpar
3172 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3173 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3174 0 : relax_opacity_max_adjust_model = keep_going
3175 0 : end function relax_opacity_max_adjust_model
3176 :
3177 0 : integer function relax_opacity_max_check_model(s, id, lipar, ipar, lrpar, rpar)
3178 : use do_one_utils, only:do_bare_bones_check_model
3179 : type (star_info), pointer :: s
3180 : integer, intent(in) :: id, lipar, lrpar
3181 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3182 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3183 0 : real(dp) :: new_value, per_step_multiplier
3184 : logical, parameter :: dbg = .false.
3185 :
3186 : include 'formats'
3187 :
3188 0 : relax_opacity_max_check_model = do_bare_bones_check_model(id)
3189 0 : if (relax_opacity_max_check_model /= keep_going) return
3190 :
3191 0 : new_value = rpar(1)
3192 0 : per_step_multiplier = rpar(2)
3193 :
3194 0 : s% opacity_max = s% opacity_max * per_step_multiplier
3195 :
3196 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
3197 0 : write(*,1) 'relax opacity', s% opacity_max, new_value
3198 :
3199 0 : if ((per_step_multiplier < 1 .and. s% opacity_max < new_value) .or. &
3200 : (per_step_multiplier > 1 .and. s% opacity_max > new_value)) then
3201 0 : s% opacity_max = new_value
3202 0 : relax_opacity_max_check_model = terminate
3203 0 : s% termination_code = t_relax_finished_okay
3204 0 : rpar(1) = s% dt
3205 0 : return
3206 : end if
3207 :
3208 0 : end function relax_opacity_max_check_model
3209 :
3210 :
3211 0 : subroutine do_relax_max_surf_dq(id, new_value, per_step_multiplier, ierr)
3212 : integer, intent(in) :: id
3213 : real(dp), intent(in) :: new_value, per_step_multiplier
3214 : integer, intent(out) :: ierr
3215 : integer, parameter :: lipar=1, lrpar=2
3216 : integer :: max_model_number
3217 : type (star_info), pointer :: s
3218 : integer, target :: ipar_ary(lipar)
3219 : integer, pointer :: ipar(:)
3220 0 : real(dp), target :: rpar_ary(lrpar)
3221 : real(dp), pointer :: rpar(:)
3222 0 : rpar => rpar_ary
3223 0 : ipar => ipar_ary
3224 : include 'formats'
3225 : ierr = 0
3226 0 : if (new_value <= 0) then
3227 0 : ierr = -1
3228 0 : write(*,*) 'invalid new_value', new_value
3229 0 : return
3230 : end if
3231 0 : if (per_step_multiplier <= 0 .or. per_step_multiplier == 1) then
3232 0 : ierr = -1
3233 0 : write(*,*) 'invalid per_step_multiplier', per_step_multiplier
3234 0 : return
3235 : end if
3236 0 : call get_star_ptr(id, s, ierr)
3237 0 : if (ierr /= 0) return
3238 0 : if (s% max_surface_cell_dq <= 0) then
3239 0 : ierr = -1
3240 0 : write(*,*) 'invalid max_surf_dq', s% max_surface_cell_dq
3241 0 : return
3242 : end if
3243 0 : if (abs(new_value - s% max_surface_cell_dq) <= &
3244 : 1d-6*min(new_value,s% max_surface_cell_dq)) then
3245 0 : s% max_surface_cell_dq = new_value
3246 0 : return
3247 : end if
3248 0 : write(*,'(A)')
3249 0 : write(*,1) 'current max_surf_dq', s% max_surface_cell_dq
3250 0 : write(*,1) 'relax to new_value', new_value
3251 0 : write(*,'(A)')
3252 0 : write(*,1) 'per_step_multiplier', per_step_multiplier
3253 0 : write(*,'(A)')
3254 0 : rpar(1) = new_value
3255 0 : rpar(2) = per_step_multiplier
3256 0 : max_model_number = s% max_model_number
3257 0 : s% max_model_number = -1111
3258 : call do_internal_evolve( &
3259 : id, before_evolve_relax_max_surf_dq, &
3260 : relax_max_surf_dq_adjust_model, relax_max_surf_dq_check_model, &
3261 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3262 0 : s% max_model_number = max_model_number
3263 0 : s% max_surface_cell_dq = new_value
3264 0 : s% dt_next = rpar(1) ! keep dt from relax
3265 0 : call error_check('relax max surf dq',ierr)
3266 0 : end subroutine do_relax_max_surf_dq
3267 :
3268 :
3269 0 : subroutine before_evolve_relax_max_surf_dq(s, id, lipar, ipar, lrpar, rpar, ierr)
3270 : type (star_info), pointer :: s
3271 : integer, intent(in) :: id, lipar, lrpar
3272 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3273 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3274 : integer, intent(out) :: ierr
3275 0 : ierr = 0
3276 0 : s% max_model_number = -111
3277 0 : s% dt_next = secyer*1d-3
3278 0 : call turn_off_winds(s)
3279 0 : end subroutine before_evolve_relax_max_surf_dq
3280 :
3281 0 : integer function relax_max_surf_dq_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3282 : type (star_info), pointer :: s
3283 : integer, intent(in) :: id, lipar, lrpar
3284 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3285 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3286 0 : relax_max_surf_dq_adjust_model = keep_going
3287 0 : end function relax_max_surf_dq_adjust_model
3288 :
3289 0 : integer function relax_max_surf_dq_check_model(s, id, lipar, ipar, lrpar, rpar)
3290 : use do_one_utils, only:do_bare_bones_check_model
3291 : type (star_info), pointer :: s
3292 : integer, intent(in) :: id, lipar, lrpar
3293 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3294 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3295 0 : real(dp) :: new_value, per_step_multiplier
3296 : logical, parameter :: dbg = .false.
3297 :
3298 : include 'formats'
3299 :
3300 0 : relax_max_surf_dq_check_model = do_bare_bones_check_model(id)
3301 0 : if (relax_max_surf_dq_check_model /= keep_going) return
3302 :
3303 0 : new_value = rpar(1)
3304 0 : per_step_multiplier = rpar(2)
3305 :
3306 0 : s% max_surface_cell_dq = s% max_surface_cell_dq * per_step_multiplier
3307 :
3308 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
3309 0 : write(*,1) 'relax max_surface_cell_dq', s% max_surface_cell_dq, new_value
3310 :
3311 0 : if ((per_step_multiplier < 1 .and. s% max_surface_cell_dq < new_value) .or. &
3312 : (per_step_multiplier > 1 .and. s% max_surface_cell_dq > new_value)) then
3313 0 : s% max_surface_cell_dq = new_value
3314 0 : relax_max_surf_dq_check_model = terminate
3315 0 : s% termination_code = t_relax_finished_okay
3316 0 : rpar(1) = s% dt
3317 0 : return
3318 : end if
3319 :
3320 0 : end function relax_max_surf_dq_check_model
3321 :
3322 :
3323 0 : subroutine do_relax_num_steps(id, num_steps, max_timestep, ierr)
3324 : integer, intent(in) :: id, num_steps
3325 : real(dp), intent(in) :: max_timestep
3326 : integer, intent(out) :: ierr
3327 : integer, parameter :: lipar=1, lrpar=1
3328 : integer :: max_model_number, model_number
3329 0 : real(dp) :: save_max_timestep, save_max_years_for_timestep
3330 : logical :: save_use_gradT
3331 : type (star_info), pointer :: s
3332 : integer, target :: ipar_ary(lipar)
3333 : integer, pointer :: ipar(:)
3334 0 : real(dp), target :: rpar_ary(lrpar)
3335 : real(dp), pointer :: rpar(:)
3336 0 : rpar => rpar_ary
3337 0 : ipar => ipar_ary
3338 :
3339 : include 'formats'
3340 0 : ierr = 0
3341 0 : if (num_steps <= 0) return
3342 :
3343 :
3344 0 : call get_star_ptr(id, s, ierr)
3345 0 : if (ierr /= 0) return
3346 :
3347 0 : write(*,'(A)')
3348 0 : write(*,2) 'relax_num_steps', num_steps
3349 0 : write(*,'(A)')
3350 0 : ipar(1) = num_steps
3351 0 : if (max_timestep <= 0) then
3352 0 : rpar(1) = secyer
3353 : else
3354 0 : rpar(1) = max_timestep
3355 : end if
3356 0 : max_model_number = s% max_model_number
3357 0 : model_number = s% model_number
3358 0 : save_max_timestep = s% max_timestep
3359 0 : save_max_years_for_timestep = s% max_years_for_timestep
3360 0 : save_use_gradT = s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn
3361 0 : s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.
3362 0 : s% model_number = 0
3363 : call do_internal_evolve( &
3364 : id, before_evolve_relax_num_steps, &
3365 : relax_num_steps_adjust_model, relax_num_steps_check_model, &
3366 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3367 0 : s% max_model_number = max_model_number
3368 0 : s% model_number = model_number
3369 0 : s% max_timestep = save_max_timestep
3370 0 : s% max_years_for_timestep = save_max_years_for_timestep
3371 0 : s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = save_use_gradT
3372 0 : call error_check('relax num steps',ierr)
3373 :
3374 0 : end subroutine do_relax_num_steps
3375 :
3376 :
3377 0 : subroutine before_evolve_relax_num_steps(s, id, lipar, ipar, lrpar, rpar, ierr)
3378 : type (star_info), pointer :: s
3379 : integer, intent(in) :: id, lipar, lrpar
3380 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3381 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3382 : integer, intent(out) :: ierr
3383 0 : ierr = 0
3384 0 : call setup_before_relax(s)
3385 0 : s% max_timestep = rpar(1)
3386 0 : s% max_years_for_timestep = s% max_timestep/secyer
3387 0 : s% dt_next = s% max_timestep
3388 0 : s% max_model_number = ipar(1)
3389 0 : end subroutine before_evolve_relax_num_steps
3390 :
3391 0 : integer function relax_num_steps_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3392 : type (star_info), pointer :: s
3393 : integer, intent(in) :: id, lipar, lrpar
3394 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3395 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3396 0 : relax_num_steps_adjust_model = keep_going
3397 0 : end function relax_num_steps_adjust_model
3398 :
3399 0 : integer function relax_num_steps_check_model(s, id, lipar, ipar, lrpar, rpar)
3400 : use do_one_utils, only:do_bare_bones_check_model
3401 :
3402 : type (star_info), pointer :: s
3403 : integer, intent(in) :: id, lipar, lrpar
3404 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3405 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3406 :
3407 : logical, parameter :: dbg = .false.
3408 :
3409 : include 'formats'
3410 :
3411 0 : relax_num_steps_check_model = do_bare_bones_check_model(id)
3412 0 : if (relax_num_steps_check_model /= keep_going) return
3413 0 : if (s% model_number >= ipar(1)) then
3414 0 : relax_num_steps_check_model = terminate
3415 0 : s% termination_code = t_relax_finished_okay
3416 : end if
3417 :
3418 0 : end function relax_num_steps_check_model
3419 :
3420 :
3421 0 : subroutine do_relax_to_radiative_core(id, ierr)
3422 : integer, intent(in) :: id
3423 : integer, intent(out) :: ierr
3424 : integer, parameter :: lipar=1, lrpar=2
3425 : integer :: max_model_number, model_number
3426 0 : real(dp) :: save_max_timestep, save_max_years_for_timestep
3427 : type (star_info), pointer :: s
3428 : integer, target :: ipar_ary(lipar)
3429 : integer, pointer :: ipar(:)
3430 0 : real(dp), target :: rpar_ary(lrpar)
3431 : real(dp), pointer :: rpar(:)
3432 0 : real(dp) :: max_timestep
3433 0 : rpar => rpar_ary
3434 0 : ipar => ipar_ary
3435 :
3436 : include 'formats'
3437 : ierr = 0
3438 :
3439 0 : call get_star_ptr(id, s, ierr)
3440 0 : if (ierr /= 0) return
3441 :
3442 0 : if (s% star_mass < s% job% pre_ms_check_radiative_core_min_mass) then
3443 0 : write(*,*) 'stop relax to begin radiative core because star_mass < pre_ms_check_radiative_core_min_mass'
3444 0 : return
3445 : end if
3446 :
3447 0 : max_timestep = 1d3*secyer ! can provide a parameter for this if necessary
3448 :
3449 0 : write(*,'(A)')
3450 0 : write(*,1) 'relax_to_radiative_core'
3451 0 : write(*,'(A)')
3452 : if (max_timestep <= 0) then
3453 : rpar(1) = secyer
3454 : else
3455 0 : rpar(1) = max_timestep
3456 : end if
3457 0 : rpar(2) = 1d99 ! min_conv_mx1_bot
3458 0 : max_model_number = s% max_model_number
3459 0 : model_number = s% model_number
3460 0 : save_max_timestep = s% max_timestep
3461 0 : save_max_years_for_timestep = s% max_years_for_timestep
3462 0 : s% model_number = 0
3463 : call do_internal_evolve( &
3464 : id, before_evolve_relax_to_radiative_core, &
3465 : relax_to_radiative_core_adjust_model, relax_to_radiative_core_check_model, &
3466 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3467 0 : s% max_model_number = max_model_number
3468 0 : s% model_number = model_number
3469 0 : s% max_timestep = save_max_timestep
3470 0 : s% max_years_for_timestep = save_max_years_for_timestep
3471 0 : call error_check('relax_to_radiative_core',ierr)
3472 :
3473 0 : end subroutine do_relax_to_radiative_core
3474 :
3475 0 : subroutine before_evolve_relax_to_radiative_core(s, id, lipar, ipar, lrpar, rpar, ierr)
3476 : type (star_info), pointer :: s
3477 : integer, intent(in) :: id, lipar, lrpar
3478 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3479 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3480 : integer, intent(out) :: ierr
3481 0 : ierr = 0
3482 0 : call setup_before_relax(s)
3483 : !s% max_timestep = rpar(1)
3484 : !s% max_years_for_timestep = s% max_timestep/secyer
3485 : !s% dt_next = s% max_timestep
3486 0 : end subroutine before_evolve_relax_to_radiative_core
3487 :
3488 0 : integer function relax_to_radiative_core_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3489 : type (star_info), pointer :: s
3490 : integer, intent(in) :: id, lipar, lrpar
3491 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3492 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3493 0 : relax_to_radiative_core_adjust_model = keep_going
3494 0 : end function relax_to_radiative_core_adjust_model
3495 :
3496 0 : integer function relax_to_radiative_core_check_model(s, id, lipar, ipar, lrpar, rpar)
3497 : use do_one_utils, only:do_bare_bones_check_model
3498 : use report, only: set_power_info
3499 : type (star_info), pointer :: s
3500 : integer, intent(in) :: id, lipar, lrpar
3501 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3502 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3503 : logical, parameter :: dbg = .false.
3504 : integer :: ierr
3505 0 : real(dp) :: min_conv_mx1_bot
3506 : include 'formats'
3507 0 : relax_to_radiative_core_check_model = do_bare_bones_check_model(id)
3508 0 : if (relax_to_radiative_core_check_model /= keep_going) return
3509 0 : ierr = 0
3510 0 : call set_power_info(s)
3511 0 : if (s% L_nuc_burn_total/s% L_phot > s% job% pre_ms_check_radiative_core_Lnuc_div_L_limit) then
3512 0 : write(*,*) 'reached pre_ms_check_radiative_core_Lnuc_div_L_limit in relax to begin radiative core'
3513 0 : relax_to_radiative_core_check_model = terminate
3514 0 : s% termination_code = t_relax_finished_okay
3515 0 : return
3516 : end if
3517 0 : min_conv_mx1_bot = rpar(2)
3518 0 : if (s% conv_mx1_bot < min_conv_mx1_bot) then
3519 0 : min_conv_mx1_bot = s% conv_mx1_bot
3520 0 : rpar(2) = min_conv_mx1_bot
3521 : end if
3522 0 : if (min_conv_mx1_bot < s% job% pre_ms_check_radiative_core_start) then
3523 0 : if (s% conv_mx1_bot > s% job% pre_ms_check_radiative_core_stop) then
3524 0 : write(*,2) 'finished relax to begin radiative core', s% model_number, s% conv_mx1_bot
3525 0 : relax_to_radiative_core_check_model = terminate
3526 0 : s% termination_code = t_relax_finished_okay
3527 0 : else if (mod(s% model_number, s% terminal_interval) == 0) then
3528 0 : write(*,2) 'relative mass of radiative core still tiny', s% model_number, s% conv_mx1_bot
3529 : end if
3530 0 : else if (mod(s% model_number, s% terminal_interval) == 0) then
3531 0 : write(*,2) 'waiting for fully convective core to develop', s% model_number, s% conv_mx1_bot
3532 : end if
3533 0 : end function relax_to_radiative_core_check_model
3534 :
3535 :
3536 0 : subroutine do_relax_Z(id, new_z, dlnz, minq, maxq, ierr)
3537 0 : use star_utils, only: eval_current_z
3538 : use adjust_xyz, only:set_z
3539 : integer, intent(in) :: id
3540 : real(dp), intent(in) :: new_z, dlnz, minq, maxq
3541 : integer, intent(out) :: ierr
3542 : integer, parameter :: lipar=1, lrpar=5
3543 : integer :: max_model_number
3544 0 : real(dp) :: z
3545 : type (star_info), pointer :: s
3546 : integer, target :: ipar_ary(lipar)
3547 : integer, pointer :: ipar(:)
3548 0 : real(dp), target :: rpar_ary(lrpar)
3549 : real(dp), pointer :: rpar(:)
3550 0 : rpar => rpar_ary
3551 0 : ipar => ipar_ary
3552 : include 'formats'
3553 : ierr = 0
3554 0 : if (new_Z < 0 .or. new_Z > 1) then
3555 0 : ierr = -1
3556 0 : write(*,*) 'invalid new_Z', new_Z
3557 0 : return
3558 : end if
3559 0 : call get_star_ptr(id, s, ierr)
3560 0 : if (ierr /= 0) return
3561 0 : z = eval_current_z(s, 1, s% nz, ierr)
3562 0 : if (ierr /= 0) return
3563 0 : if (abs(new_z - z) <= 1d-6*z) return
3564 0 : if (max(new_z, z) > 1d-6) then
3565 0 : if (abs(new_z - z) <= 1d-3*new_z) then
3566 0 : call set_z(s, new_z, 1, s% nz, ierr)
3567 0 : return
3568 : end if
3569 : end if
3570 0 : write(*,'(A)')
3571 0 : write(*,1) 'current Z', z
3572 0 : write(*,1) 'relax to new Z', new_z
3573 0 : write(*,1) '(new - current) / current', (new_z - z) / z
3574 0 : write(*,'(A)')
3575 0 : write(*,1) 'dlnz per step', dlnz
3576 0 : write(*,'(A)')
3577 0 : rpar(1) = log(max(min_z,new_z))
3578 0 : rpar(2) = new_z
3579 0 : rpar(3) = abs(dlnz)
3580 0 : rpar(4) = minq
3581 0 : rpar(5) = maxq
3582 0 : max_model_number = s% max_model_number
3583 0 : s% max_model_number = -1111
3584 0 : s% initial_z = z
3585 : call do_internal_evolve( &
3586 : id, before_evolve_relax_Z, &
3587 : relax_Z_adjust_model, relax_Z_check_model, &
3588 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3589 0 : s% max_model_number = max_model_number
3590 0 : call error_check('relax Z',ierr)
3591 0 : end subroutine do_relax_Z
3592 :
3593 :
3594 0 : subroutine before_evolve_relax_Z(s, id, lipar, ipar, lrpar, rpar, ierr)
3595 : type (star_info), pointer :: s
3596 : integer, intent(in) :: id, lipar, lrpar
3597 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3598 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3599 : integer, intent(out) :: ierr
3600 0 : ierr = 0
3601 0 : call setup_before_relax(s)
3602 0 : s% max_model_number = -111
3603 0 : s% max_timestep = secyer
3604 0 : s% dt_next = s% max_timestep
3605 0 : end subroutine before_evolve_relax_Z
3606 :
3607 0 : integer function relax_Z_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3608 : type (star_info), pointer :: s
3609 : integer, intent(in) :: id, lipar, lrpar
3610 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3611 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3612 0 : relax_Z_adjust_model = keep_going
3613 0 : end function relax_Z_adjust_model
3614 :
3615 0 : integer function relax_Z_check_model(s, id, lipar, ipar, lrpar, rpar)
3616 : use adjust_xyz, only: set_z
3617 : use star_utils, only: k_for_q, eval_current_z
3618 : use do_one_utils, only:do_bare_bones_check_model
3619 : type (star_info), pointer :: s
3620 : integer, intent(in) :: id, lipar, lrpar
3621 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3622 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3623 : integer :: ierr, klo, khi
3624 0 : real(dp) :: lnz_target, new_z, new_lnz, dlnz, lnz, current_z, next_z, &
3625 0 : min_q_for_relax_Z, max_q_for_relax_Z
3626 :
3627 : logical, parameter :: zdbg = .true.
3628 :
3629 : include 'formats'
3630 :
3631 0 : relax_Z_check_model = do_bare_bones_check_model(id)
3632 0 : if (relax_Z_check_model /= keep_going) return
3633 :
3634 0 : lnz_target = rpar(1)
3635 0 : new_z = rpar(2)
3636 0 : dlnz = rpar(3)
3637 0 : min_q_for_relax_Z = rpar(4)
3638 0 : max_q_for_relax_Z = rpar(5)
3639 :
3640 0 : khi = k_for_q(s, min_q_for_relax_Z)
3641 0 : klo = k_for_q(s, max_q_for_relax_Z)
3642 0 : if (zdbg) write(*,2) 'klo', klo, max_q_for_relax_Z
3643 0 : if (zdbg) write(*,2) 'khi', khi, min_q_for_relax_Z
3644 :
3645 0 : current_z = eval_current_z(s, klo, khi, ierr)
3646 0 : if (ierr /= 0) return
3647 :
3648 0 : if (mod(s% model_number, s% terminal_interval) == 0) &
3649 0 : write(*,1) 'new_z, current', new_z, current_z
3650 :
3651 0 : if (abs(current_z-new_z) <= 1d-6*new_z) then
3652 0 : relax_Z_check_model = terminate
3653 0 : s% termination_code = t_relax_finished_okay
3654 0 : return
3655 : end if
3656 :
3657 0 : lnz = log(max(min_z,current_z))
3658 :
3659 : if (zdbg) then
3660 0 : write(*,1) 'lnz_target', lnz_target
3661 0 : write(*,1) 'lnz', lnz
3662 0 : write(*,1) 'lnz - lnz_target', lnz - lnz_target
3663 0 : write(*,1) 'dlnz', dlnz
3664 : end if
3665 :
3666 0 : if (abs(lnz - lnz_target) < dlnz) then
3667 0 : dlnz = abs(lnz - lnz_target)
3668 0 : if (zdbg) write(*,1) 'reduced dlnz', dlnz
3669 : end if
3670 :
3671 0 : if (lnz_target < lnz) then
3672 0 : new_lnz = lnz - dlnz
3673 : else
3674 0 : new_lnz = lnz + dlnz
3675 : end if
3676 :
3677 0 : if (zdbg) write(*,1) 'new_lnz', new_lnz
3678 :
3679 0 : if (new_lnz >= min_dlnz) then
3680 0 : next_z = exp(new_lnz)
3681 : else
3682 0 : next_z = new_z
3683 : end if
3684 :
3685 0 : if (zdbg) write(*,1) 'next_z', next_z
3686 :
3687 : ierr = 0
3688 0 : call set_z(s, next_z, klo, khi, ierr)
3689 0 : if (ierr /= 0) then
3690 0 : if (s% report_ierr) write(*, *) 'relax_Z_check_model ierr', ierr
3691 0 : relax_Z_check_model = terminate
3692 0 : s% result_reason = nonzero_ierr
3693 0 : return
3694 : end if
3695 :
3696 0 : write(*,1) 'relax Z, z diff, new, current', new_z - current_z, new_z, current_z
3697 :
3698 0 : if (klo == 1 .and. khi == s% nz) s% initial_z = next_z
3699 0 : s% max_timestep = secyer*s% time_step
3700 :
3701 0 : end function relax_Z_check_model
3702 :
3703 :
3704 0 : subroutine do_relax_Y(id, new_Y, dY, minq, maxq, ierr)
3705 0 : use star_utils, only: k_for_q, eval_current_y
3706 : integer, intent(in) :: id
3707 : real(dp), intent(in) :: new_Y, dY, minq, maxq
3708 : integer, intent(out) :: ierr
3709 : integer, parameter :: lipar=1, lrpar=4
3710 : integer :: max_model_number, khi, klo
3711 0 : real(dp) :: y
3712 : type (star_info), pointer :: s
3713 : integer, target :: ipar_ary(lipar)
3714 : integer, pointer :: ipar(:)
3715 0 : real(dp), target :: rpar_ary(lrpar)
3716 : real(dp), pointer :: rpar(:)
3717 0 : rpar => rpar_ary
3718 0 : ipar => ipar_ary
3719 : include 'formats'
3720 : ierr = 0
3721 0 : call get_star_ptr(id, s, ierr)
3722 0 : if (ierr /= 0) return
3723 :
3724 0 : khi = k_for_q(s, minq)
3725 0 : klo = k_for_q(s, maxq)
3726 :
3727 0 : y = eval_current_y(s, klo, khi, ierr)
3728 0 : if (ierr /= 0) return
3729 0 : if (is_bad(y)) then
3730 0 : write(*,1) 'y', y
3731 0 : call mesa_error(__FILE__,__LINE__,'do_relax_Y')
3732 : end if
3733 :
3734 0 : if (abs(new_Y - y) <= 1d-6*new_Y) return
3735 0 : if (new_Y < 0 .or. new_Y > 1) then
3736 0 : ierr = -1
3737 0 : write(*,*) 'invalid new_Y', new_Y
3738 0 : return
3739 : end if
3740 0 : write(*,'(A)')
3741 0 : write(*,1) 'current Y', Y
3742 0 : write(*,1) 'relax to new_Y', new_Y
3743 0 : write(*,1) 'dY per step', dY
3744 0 : write(*,'(A)')
3745 0 : rpar(1) = new_Y
3746 0 : rpar(2) = abs(dY)
3747 0 : rpar(3) = minq
3748 0 : rpar(4) = maxq
3749 0 : max_model_number = s% max_model_number
3750 0 : s% max_model_number = -1111
3751 0 : s% initial_y = y
3752 : call do_internal_evolve( &
3753 : id, before_evolve_relax_Y, &
3754 : relax_Y_adjust_model, relax_Y_check_model, &
3755 0 : null_finish_model, .true., lipar, ipar, lrpar, rpar, ierr)
3756 0 : s% max_model_number = max_model_number
3757 0 : call error_check('relax Y',ierr)
3758 0 : end subroutine do_relax_Y
3759 :
3760 :
3761 0 : subroutine before_evolve_relax_Y(s, id, lipar, ipar, lrpar, rpar, ierr)
3762 : type (star_info), pointer :: s
3763 : integer, intent(in) :: id, lipar, lrpar
3764 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3765 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3766 : integer, intent(out) :: ierr
3767 0 : ierr = 0
3768 0 : call setup_before_relax(s)
3769 0 : s% max_model_number = -111
3770 0 : s% max_timestep = secyer
3771 0 : s% dt_next = s% max_timestep
3772 0 : end subroutine before_evolve_relax_Y
3773 :
3774 0 : integer function relax_Y_adjust_model(s, id, lipar, ipar, lrpar, rpar)
3775 : type (star_info), pointer :: s
3776 : integer, intent(in) :: id, lipar, lrpar
3777 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3778 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3779 0 : relax_Y_adjust_model = keep_going
3780 0 : end function relax_Y_adjust_model
3781 :
3782 0 : integer function relax_Y_check_model(s, id, lipar, ipar, lrpar, rpar)
3783 : use star_utils, only: k_for_q, eval_current_y
3784 : use adjust_xyz, only: set_y
3785 : use do_one_utils, only: do_bare_bones_check_model
3786 : type (star_info), pointer :: s
3787 : integer, intent(in) :: id, lipar, lrpar
3788 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3789 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3790 : integer :: ierr, klo, khi
3791 0 : real(dp) :: new_y, dy, current_y, next_y, minq, maxq, actual_next_y
3792 : logical, parameter :: ydbg = .true.
3793 : logical :: end_now
3794 :
3795 : include 'formats'
3796 :
3797 0 : end_now=.false.
3798 :
3799 0 : relax_Y_check_model = do_bare_bones_check_model(id)
3800 0 : if (relax_Y_check_model /= keep_going) return
3801 :
3802 0 : new_y = rpar(1)
3803 0 : dy = rpar(2)
3804 0 : minq = rpar(3)
3805 0 : maxq = rpar(4)
3806 :
3807 0 : khi = k_for_q(s, minq)
3808 0 : klo = k_for_q(s, maxq)
3809 :
3810 : if (ydbg) then
3811 0 : write(*,4) 'klo, khi nz', klo, khi, s% nz
3812 : end if
3813 :
3814 0 : current_y = eval_current_y(s, klo, khi, ierr)
3815 0 : if (is_bad(current_y)) then
3816 0 : write(*,1) 'current_y', current_y
3817 0 : call mesa_error(__FILE__,__LINE__,'relax_y_check_model')
3818 : end if
3819 0 : if (ierr /= 0) return
3820 :
3821 0 : if (mod(s% model_number, s% terminal_interval) == 0) then
3822 0 : write(*,1) 'new_y', new_y
3823 0 : write(*,1) 'dy', dy
3824 0 : write(*,1) 'current_y', current_y
3825 0 : write(*,1) 'current_y - new_y', current_y - new_y
3826 : end if
3827 :
3828 0 : if (abs(current_y - new_y) < 1d-10) then
3829 0 : relax_Y_check_model = terminate
3830 0 : s% termination_code = t_relax_finished_okay
3831 0 : return
3832 : end if
3833 :
3834 0 : if (abs(current_y - new_y) < dY) then
3835 0 : dY = abs(current_y - new_y)
3836 0 : end_now = .true.
3837 0 : if (ydbg) write(*,1) 'reduced dY', dY
3838 : end if
3839 :
3840 0 : if (new_y >= current_y) then
3841 0 : next_y = current_y + dy
3842 : else
3843 0 : next_y = current_y - dy
3844 : end if
3845 :
3846 0 : if (ydbg) write(*,1) 'next_y', next_y
3847 :
3848 : ierr = 0
3849 :
3850 0 : call set_y(s, next_y, klo, khi, ierr)
3851 0 : if (ierr /= 0) then
3852 0 : if (s% report_ierr) write(*, *) 'relax_Y_check_model ierr', ierr
3853 0 : relax_Y_check_model = terminate
3854 0 : s% result_reason = nonzero_ierr
3855 0 : return
3856 : end if
3857 :
3858 0 : actual_next_y = eval_current_y(s, klo, khi, ierr)
3859 :
3860 0 : write(*,1) 'relax Y, y diff, new, current', new_y - current_y, new_y, actual_next_y
3861 :
3862 0 : if (ydbg) write(*,1) 'actual_next_y', actual_next_y
3863 0 : if (ydbg) write(*,1) 'actual_next_y - next_y', actual_next_y - next_y
3864 0 : if (ydbg) write(*,1) 'y diff', actual_next_y - new_y
3865 :
3866 0 : if (abs(actual_next_y - new_y) < 1d-10 .or. end_now) then
3867 0 : relax_Y_check_model = terminate
3868 0 : s% termination_code = t_relax_finished_okay
3869 0 : return
3870 : end if
3871 :
3872 0 : if (klo == 1 .and. khi == s% nz) s% initial_y = next_y
3873 0 : s% max_timestep = secyer*s% time_step
3874 :
3875 0 : end function relax_Y_check_model
3876 :
3877 :
3878 0 : subroutine setup_before_relax(s)
3879 : type (star_info), pointer :: s
3880 0 : s% dxdt_nuc_factor = 0
3881 0 : s% max_age = 1d50
3882 0 : s% max_age_in_days = 1d50
3883 0 : s% max_age_in_seconds = 1d50
3884 0 : s% max_timestep_factor = 2
3885 0 : s% max_model_number = -1111
3886 0 : call turn_off_winds(s)
3887 0 : end subroutine setup_before_relax
3888 :
3889 :
3890 0 : subroutine turn_off_winds(s)
3891 : type (star_info), pointer :: s
3892 0 : s% mass_change = 0
3893 0 : s% Reimers_scaling_factor = 0d0
3894 0 : s% Blocker_scaling_factor = 0d0
3895 0 : s% de_Jager_scaling_factor = 0d0
3896 0 : s% van_Loon_scaling_factor = 0d0
3897 0 : s% Nieuwenhuijzen_scaling_factor = 0d0
3898 0 : s% Vink_scaling_factor = 0d0
3899 0 : s% Dutch_scaling_factor = 0d0
3900 0 : s% Bjorklund_scaling_factor = 0d0
3901 0 : s% use_other_wind = .false.
3902 0 : end subroutine turn_off_winds
3903 :
3904 :
3905 0 : subroutine do_internal_evolve( &
3906 : id, before_evolve, adjust_model, check_model,&
3907 : finish_model, restore_at_end, &
3908 : lipar, ipar, lrpar, rpar, ierr)
3909 : use evolve
3910 : use star_utils, only: yrs_for_init_timestep
3911 : use pgstar
3912 : use history_specs, only: set_history_columns
3913 : use profile, only: set_profile_columns
3914 : integer, intent(in) :: id, lipar, lrpar
3915 : logical, intent(in) :: restore_at_end
3916 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3917 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3918 : interface
3919 : subroutine before_evolve(s, id, lipar, ipar, lrpar, rpar, ierr)
3920 : use const_def, only: dp
3921 : use star_private_def, only:star_info
3922 : implicit none
3923 : type (star_info), pointer :: s
3924 : integer, intent(in) :: id, lipar, lrpar
3925 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3926 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3927 : integer, intent(out) :: ierr
3928 : end subroutine before_evolve
3929 : integer function adjust_model(s, id, lipar, ipar, lrpar, rpar)
3930 : use const_def, only: dp
3931 : use star_private_def, only:star_info
3932 : implicit none
3933 : type (star_info), pointer :: s
3934 : integer, intent(in) :: id, lipar, lrpar
3935 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3936 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3937 : end function adjust_model
3938 : integer function check_model(s, id, lipar, ipar, lrpar, rpar)
3939 : use const_def, only: dp
3940 : use star_private_def, only:star_info
3941 : implicit none
3942 : type (star_info), pointer :: s
3943 : integer, intent(in) :: id, lipar, lrpar
3944 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
3945 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
3946 : end function check_model
3947 : integer function finish_model(s)
3948 : use star_def, only:star_info
3949 : implicit none
3950 : type (star_info), pointer :: s
3951 : end function finish_model
3952 : end interface
3953 : integer, intent(out) :: ierr
3954 : type (star_info), pointer :: s
3955 : character (len=10) :: MLT_option
3956 : integer :: result, model_number, model_number_for_last_retry, &
3957 : recent_log_header, num_retries, &
3958 : photo_interval, profile_interval, priority_profile_interval, &
3959 : model_number_old, max_number_retries, &
3960 : solver_iters_timestep_limit, iter_for_resid_tol2, iter_for_resid_tol3, &
3961 : steps_before_use_gold_tolerances, steps_before_use_gold2_tolerances
3962 0 : real(dp) :: star_age, time, max_age, max_age_in_days, max_age_in_seconds, max_timestep, &
3963 0 : Reimers_scaling_factor, Blocker_scaling_factor, de_Jager_scaling_factor, Dutch_scaling_factor, &
3964 0 : van_Loon_scaling_factor, Nieuwenhuijzen_scaling_factor, Vink_scaling_factor, Bjorklund_scaling_factor,&
3965 0 : dxdt_nuc_factor, tol_correction_norm, tol_max_correction, warning_limit_for_max_residual, &
3966 0 : tol_residual_norm1, tol_max_residual1, &
3967 0 : tol_residual_norm2, tol_max_residual2, &
3968 0 : tol_residual_norm3, tol_max_residual3, &
3969 0 : max_timestep_factor, mass_change, varcontrol_target, dt_next, &
3970 0 : time_old, maxT_for_gold_tolerances
3971 : logical :: do_history_file, write_profiles_flag, first_try, use_other_wind, &
3972 : use_gold_tolerances, use_gold2_tolerances
3973 :
3974 : procedure(integer), pointer :: tmp_ptr1 => null(), tmp_ptr3 => null()
3975 : procedure(), pointer :: tmp_ptr2 => null(), tmp_ptr4 => null()
3976 : logical, parameter :: dbg = .false.
3977 :
3978 : include 'formats'
3979 :
3980 0 : ierr = 0
3981 0 : call get_star_ptr(id, s, ierr)
3982 0 : if (ierr /= 0) return
3983 :
3984 0 : call save_stuff
3985 :
3986 0 : s% do_history_file = .false.
3987 0 : s% write_profiles_flag = .false.
3988 0 : s% warning_limit_for_max_residual = 1d99
3989 0 : s% recent_log_header = -1
3990 0 : s% max_number_retries = s% relax_max_number_retries
3991 0 : s% use_gold_tolerances = s% relax_use_gold_tolerances
3992 0 : s% steps_before_use_gold_tolerances = -1
3993 0 : s% use_gold2_tolerances = .false.
3994 0 : s% steps_before_use_gold2_tolerances = -1
3995 0 : if (s% MLT_option == 'TDC') then
3996 0 : s% MLT_option = 'Cox'
3997 : end if
3998 :
3999 0 : if (s% relax_solver_iters_timestep_limit /= 0) &
4000 0 : s% solver_iters_timestep_limit = s% relax_solver_iters_timestep_limit
4001 :
4002 0 : if (s% relax_tol_correction_norm /= 0) &
4003 0 : s% tol_correction_norm = s% relax_tol_correction_norm
4004 0 : if (s% relax_tol_max_correction /= 0) &
4005 0 : s% tol_max_correction = s% relax_tol_max_correction
4006 :
4007 0 : if (s% relax_iter_for_resid_tol2 /= 0) &
4008 0 : s% iter_for_resid_tol2 = s% relax_iter_for_resid_tol2
4009 0 : if (s% relax_tol_residual_norm1 /= 0) &
4010 0 : s% tol_residual_norm1 = s% relax_tol_residual_norm1
4011 0 : if (s% relax_tol_max_residual1 /= 0) &
4012 0 : s% tol_max_residual1 = s% relax_tol_max_residual1
4013 :
4014 0 : if (s% relax_iter_for_resid_tol3 /= 0) &
4015 0 : s% iter_for_resid_tol3 = s% relax_iter_for_resid_tol3
4016 0 : if (s% relax_tol_residual_norm2 /= 0) &
4017 0 : s% tol_residual_norm2 = s% relax_tol_residual_norm2
4018 0 : if (s% relax_tol_max_residual2 /= 0) &
4019 0 : s% tol_max_residual2 = s% relax_tol_max_residual2
4020 :
4021 0 : if (s% relax_tol_residual_norm3 /= 0) &
4022 0 : s% tol_residual_norm3 = s% relax_tol_residual_norm3
4023 0 : if (s% relax_tol_max_residual3 /= 0) &
4024 0 : s% tol_max_residual3 = s% relax_tol_max_residual3
4025 :
4026 0 : if (s% relax_maxT_for_gold_tolerances /= 0) &
4027 0 : s% maxT_for_gold_tolerances = s% relax_maxT_for_gold_tolerances
4028 :
4029 0 : if (s% doing_first_model_of_run) then
4030 0 : s% num_retries = 0
4031 0 : s% time = 0
4032 0 : s% star_age = 0
4033 0 : s% model_number_for_last_retry = 0
4034 0 : s% photo_interval = 0
4035 0 : s% profile_interval = 0
4036 0 : s% priority_profile_interval = 0
4037 : end if
4038 :
4039 0 : if( s% job% pgstar_flag) then
4040 : ! Can't use the star_lib versions otherwise we have a circular dependency in the makefile
4041 0 : call set_history_columns(id,s% job% history_columns_file, .true., ierr)
4042 0 : if (ierr /= 0) return
4043 0 : call set_profile_columns(id, s% job% profile_columns_file, .true., ierr)
4044 0 : if (ierr /= 0) return
4045 : end if
4046 :
4047 0 : if (s% doing_first_model_of_run .and. s% job% pgstar_flag) then
4048 0 : call do_start_new_run_for_pgstar(s, ierr)
4049 0 : if (ierr /= 0) return
4050 : end if
4051 :
4052 0 : call before_evolve(s, id, lipar, ipar, lrpar, rpar, ierr)
4053 0 : if (ierr /= 0) then
4054 0 : write(*,*) 'failed in before_evolve'
4055 0 : return
4056 : end if
4057 :
4058 0 : s% termination_code = -1
4059 0 : s% doing_relax = .true.
4060 0 : s% need_to_setvars = .true.
4061 :
4062 : evolve_loop: do ! evolve one step per loop
4063 :
4064 0 : first_try = .true.
4065 :
4066 0 : step_loop: do ! may need to repeat this loop for retry
4067 :
4068 0 : result = do_evolve_step_part1(id, first_try)
4069 0 : if (result == keep_going) &
4070 0 : result = adjust_model(s, id, lipar, ipar, lrpar, rpar)
4071 0 : if (result == keep_going) &
4072 0 : result = do_evolve_step_part2(id, first_try)
4073 :
4074 0 : if (result == keep_going) result = check_model(s, id, lipar, ipar, lrpar, rpar)
4075 0 : if (result == keep_going) result = pick_next_timestep(id)
4076 0 : if (result == keep_going) exit step_loop
4077 :
4078 0 : if (result == retry) then
4079 0 : else if (s% result_reason /= result_reason_normal) then
4080 0 : write(*, *) model_number, 'terminate reason: ' // trim(result_reason_str(s% result_reason))
4081 : end if
4082 :
4083 0 : if (result == redo) result = prepare_to_redo(id)
4084 0 : if (result == retry) result = prepare_to_retry(id)
4085 0 : if (result == terminate) exit evolve_loop
4086 0 : first_try = .false.
4087 :
4088 : end do step_loop
4089 :
4090 0 : if (.not. s% job% disable_pgstar_during_relax_flag .and. s% job% pgstar_flag) then
4091 : ! Can't use the star_lib versions otherwise we have a circular dependency in the makefile
4092 : if(dbg) write(*,2) 'after step_loop: call update_pgstar_data', s% model_number
4093 0 : call update_pgstar_data(s, ierr)
4094 : if (failed()) return
4095 0 : call do_read_pgstar_controls(s, s% inlist_fname, ierr)
4096 : if (failed()) return
4097 0 : call do_pgstar_plots( s, .false., ierr)
4098 : if (failed()) return
4099 : end if
4100 :
4101 0 : result = finish_model(s)
4102 0 : if (result /= keep_going) exit evolve_loop
4103 :
4104 0 : tmp_ptr1 => s% how_many_extra_history_columns
4105 0 : tmp_ptr2 => s% data_for_extra_history_columns
4106 0 : tmp_ptr3 => s% how_many_extra_profile_columns
4107 0 : tmp_ptr4 => s% data_for_extra_profile_columns
4108 :
4109 0 : s% how_many_extra_profile_columns => no_extra_profile_columns
4110 0 : s% data_for_extra_profile_columns => none_for_extra_profile_columns
4111 0 : s% how_many_extra_history_columns => no_extra_history_columns
4112 0 : s% data_for_extra_history_columns => none_for_extra_history_columns
4113 :
4114 0 : result = finish_step(id, ierr)
4115 0 : s% how_many_extra_history_columns => tmp_ptr1
4116 0 : s% data_for_extra_history_columns => tmp_ptr2
4117 0 : s% how_many_extra_profile_columns => tmp_ptr3
4118 0 : s% data_for_extra_profile_columns => tmp_ptr4
4119 0 : nullify(tmp_ptr1,tmp_ptr2,tmp_ptr3,tmp_ptr4)
4120 :
4121 0 : if (result /= keep_going) exit evolve_loop
4122 :
4123 0 : if (associated(s% finish_relax_step)) then
4124 0 : result = s% finish_relax_step(id)
4125 0 : if (result /= keep_going) exit evolve_loop
4126 : end if
4127 :
4128 : end do evolve_loop
4129 :
4130 0 : if (.not. s% job% disable_pgstar_during_relax_flag .and. s% job% pgstar_flag) then
4131 : ! Can't use the star_lib versions otherwise we have a circular dependency in the makefile
4132 0 : call update_pgstar_data(s, ierr)
4133 0 : if (ierr /= 0) return
4134 0 : call do_read_pgstar_controls(s, s% inlist_fname, ierr)
4135 0 : if (ierr /= 0) return
4136 0 : call do_pgstar_plots(s, s% job% save_pgstar_files_when_terminate, ierr)
4137 0 : if (ierr /= 0) return
4138 : end if
4139 :
4140 0 : s% doing_relax = .false.
4141 0 : s% need_to_setvars = .true. ! just to be safe
4142 0 : s% have_ST_start_info = .false. ! for ST time smoothing
4143 :
4144 0 : if (.not. (s% termination_code == t_relax_finished_okay .or. &
4145 0 : s% termination_code == t_extras_check_model)) ierr = -1
4146 :
4147 0 : s% termination_code = -1
4148 :
4149 0 : if (associated(s% finished_relax)) call s% finished_relax(id)
4150 :
4151 0 : if (restore_at_end) call restore_stuff
4152 :
4153 0 : if (s% job% set_cumulative_energy_error_each_relax) &
4154 0 : s% cumulative_energy_error = s% job% new_cumulative_energy_error
4155 :
4156 0 : s% dt = 0
4157 0 : s% dt_old = 0
4158 :
4159 0 : s% timestep_hold = -100
4160 0 : s% model_number_for_last_retry = -100
4161 :
4162 : contains
4163 :
4164 0 : logical function failed()
4165 0 : failed = .false.
4166 0 : if (ierr == 0) return
4167 0 : failed = .true.
4168 0 : end function failed
4169 :
4170 0 : subroutine save_stuff
4171 :
4172 0 : warning_limit_for_max_residual = s% warning_limit_for_max_residual
4173 0 : do_history_file = s% do_history_file
4174 0 : write_profiles_flag = s% write_profiles_flag
4175 0 : recent_log_header = s% recent_log_header
4176 0 : mass_change = s% mass_change
4177 0 : Reimers_scaling_factor = s% Reimers_scaling_factor
4178 0 : Blocker_scaling_factor = s% Blocker_scaling_factor
4179 0 : de_Jager_scaling_factor = s% de_Jager_scaling_factor
4180 0 : van_Loon_scaling_factor = s% van_Loon_scaling_factor
4181 0 : Nieuwenhuijzen_scaling_factor = s% Nieuwenhuijzen_scaling_factor
4182 0 : Vink_scaling_factor = s% Vink_scaling_factor
4183 0 : Dutch_scaling_factor = s% Dutch_scaling_factor
4184 0 : Bjorklund_scaling_factor = s% Bjorklund_scaling_factor
4185 0 : use_other_wind = s% use_other_wind
4186 :
4187 0 : num_retries = s% num_retries
4188 0 : star_age = s% star_age
4189 0 : time = s% time
4190 0 : model_number = s% model_number
4191 0 : dxdt_nuc_factor = s% dxdt_nuc_factor
4192 0 : max_age = s% max_age
4193 0 : max_age_in_days = s% max_age_in_days
4194 0 : max_age_in_seconds = s% max_age_in_seconds
4195 0 : max_timestep_factor = s% max_timestep_factor
4196 0 : varcontrol_target = s% varcontrol_target
4197 0 : max_timestep = s% max_timestep
4198 0 : model_number_for_last_retry = s% model_number_for_last_retry
4199 0 : photo_interval = s% photo_interval
4200 0 : profile_interval = s% profile_interval
4201 0 : priority_profile_interval = s% priority_profile_interval
4202 0 : dt_next = s% dt_next
4203 0 : max_number_retries = s% max_number_retries
4204 0 : MLT_option = s% MLT_option
4205 :
4206 0 : use_gold2_tolerances = s% use_gold2_tolerances
4207 0 : steps_before_use_gold2_tolerances = s% steps_before_use_gold2_tolerances
4208 0 : use_gold_tolerances = s% use_gold_tolerances
4209 0 : steps_before_use_gold_tolerances = s% steps_before_use_gold_tolerances
4210 0 : solver_iters_timestep_limit = s% solver_iters_timestep_limit
4211 0 : tol_correction_norm = s% tol_correction_norm
4212 0 : tol_max_correction = s% tol_max_correction
4213 0 : iter_for_resid_tol2 = s% iter_for_resid_tol2
4214 0 : tol_residual_norm1 = s% tol_residual_norm1
4215 0 : tol_max_residual1 = s% tol_max_residual1
4216 0 : iter_for_resid_tol3 = s% iter_for_resid_tol3
4217 0 : tol_residual_norm2 = s% tol_residual_norm2
4218 0 : tol_max_residual2 = s% tol_max_residual2
4219 0 : tol_residual_norm3 = s% tol_residual_norm3
4220 0 : tol_max_residual3 = s% tol_max_residual3
4221 0 : maxT_for_gold_tolerances = s% maxT_for_gold_tolerances
4222 :
4223 : ! selected history
4224 0 : time_old = s% time_old
4225 0 : model_number_old = s% model_number_old
4226 :
4227 0 : end subroutine save_stuff
4228 :
4229 0 : subroutine restore_stuff
4230 0 : s% warning_limit_for_max_residual = warning_limit_for_max_residual
4231 0 : s% do_history_file = do_history_file
4232 0 : s% write_profiles_flag = write_profiles_flag
4233 0 : s% recent_log_header = recent_log_header
4234 0 : s% mass_change = mass_change
4235 0 : s% Reimers_scaling_factor = Reimers_scaling_factor
4236 0 : s% Blocker_scaling_factor = Blocker_scaling_factor
4237 0 : s% de_Jager_scaling_factor = de_Jager_scaling_factor
4238 0 : s% van_Loon_scaling_factor = van_Loon_scaling_factor
4239 0 : s% Nieuwenhuijzen_scaling_factor = Nieuwenhuijzen_scaling_factor
4240 0 : s% Vink_scaling_factor = Vink_scaling_factor
4241 0 : s% Dutch_scaling_factor = Dutch_scaling_factor
4242 0 : s% Bjorklund_scaling_factor = Bjorklund_scaling_factor
4243 0 : s% use_other_wind = use_other_wind
4244 0 : s% num_retries = num_retries
4245 0 : s% star_age = star_age
4246 0 : s% time = time
4247 0 : s% model_number = model_number
4248 0 : s% dxdt_nuc_factor = dxdt_nuc_factor
4249 0 : s% max_age = max_age
4250 0 : s% max_age_in_days = max_age_in_days
4251 0 : s% max_age_in_seconds = max_age_in_seconds
4252 0 : s% max_timestep_factor = max_timestep_factor
4253 0 : s% varcontrol_target = varcontrol_target
4254 0 : s% max_timestep = max_timestep
4255 0 : s% model_number_for_last_retry = model_number_for_last_retry
4256 0 : s% photo_interval = photo_interval
4257 0 : s% profile_interval = profile_interval
4258 0 : s% priority_profile_interval = priority_profile_interval
4259 0 : s% dt_next = dt_next
4260 0 : s% max_number_retries = max_number_retries
4261 0 : s% MLT_option = MLT_option
4262 :
4263 0 : s% use_gold2_tolerances = use_gold2_tolerances
4264 0 : s% steps_before_use_gold2_tolerances = steps_before_use_gold2_tolerances
4265 0 : s% use_gold_tolerances = use_gold_tolerances
4266 0 : s% steps_before_use_gold_tolerances = steps_before_use_gold_tolerances
4267 0 : s% solver_iters_timestep_limit = solver_iters_timestep_limit
4268 0 : s% tol_correction_norm = tol_correction_norm
4269 0 : s% tol_max_correction = tol_max_correction
4270 0 : s% iter_for_resid_tol2 = iter_for_resid_tol2
4271 0 : s% tol_residual_norm1 = tol_residual_norm1
4272 0 : s% tol_max_residual1 = tol_max_residual1
4273 0 : s% iter_for_resid_tol3 = iter_for_resid_tol3
4274 0 : s% tol_residual_norm2 = tol_residual_norm2
4275 0 : s% tol_max_residual2 = tol_max_residual2
4276 0 : s% tol_residual_norm3 = tol_residual_norm3
4277 0 : s% tol_max_residual3 = tol_max_residual3
4278 0 : s% maxT_for_gold_tolerances = maxT_for_gold_tolerances
4279 :
4280 : ! selected history
4281 0 : s% time_old = time_old
4282 0 : s% model_number_old = model_number_old
4283 :
4284 0 : end subroutine restore_stuff
4285 :
4286 : end subroutine do_internal_evolve
4287 :
4288 :
4289 0 : integer function null_finish_model(s)
4290 : use star_def, only:star_info
4291 : type (star_info), pointer :: s
4292 0 : null_finish_model = keep_going
4293 0 : end function null_finish_model
4294 :
4295 :
4296 0 : integer function no_extra_history_columns(id)
4297 : integer, intent(in) :: id
4298 0 : no_extra_history_columns = 0
4299 0 : end function no_extra_history_columns
4300 :
4301 :
4302 0 : subroutine none_for_extra_history_columns(id, n, names, vals, ierr)
4303 : integer, intent(in) :: id, n
4304 : character (len=maxlen_history_column_name) :: names(n)
4305 : real(dp) :: vals(n)
4306 : integer, intent(out) :: ierr
4307 0 : ierr = 0
4308 0 : end subroutine none_for_extra_history_columns
4309 :
4310 :
4311 0 : integer function no_extra_profile_columns(id)
4312 : integer, intent(in) :: id
4313 0 : no_extra_profile_columns = 0
4314 0 : end function no_extra_profile_columns
4315 :
4316 :
4317 0 : subroutine none_for_extra_profile_columns(id, n, nz, names, vals, ierr)
4318 : integer, intent(in) :: id, n, nz
4319 : character (len=maxlen_profile_column_name) :: names(n)
4320 : real(dp) :: vals(nz,n)
4321 : integer, intent(out) :: ierr
4322 0 : ierr = 0
4323 0 : end subroutine none_for_extra_profile_columns
4324 :
4325 0 : subroutine error_check(name,ierr)
4326 : character(len=*), intent(in) :: name
4327 : integer, intent(in) :: ierr
4328 : include 'formats'
4329 :
4330 0 : if (ierr /= 0) then
4331 0 : write(*,*) 'failed in ', name
4332 0 : return
4333 : end if
4334 :
4335 0 : write(*,'(A)')
4336 0 : write(*,*) 'finished doing ', name
4337 0 : write(*,'(A)')
4338 : end subroutine error_check
4339 :
4340 : end module relax
|
