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