Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2019 Radek Smolec & 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 rsp
21 : use star_def, only: star_ptr, star_info
22 : use rsp_eval_eos_and_kap, only: init_for_rsp_eos_and_kap
23 : use rsp_def
24 : use rsp_step, only: calculate_energies, init_HYD, HYD
25 : use utils_lib, only: mesa_error
26 :
27 : implicit none
28 :
29 : private
30 : public :: rsp_setup_part1
31 : public :: rsp_setup_part2
32 : public :: rsp_one_step
33 : public :: build_rsp_model
34 : public :: rsp_total_energy_integrals
35 : public :: do1_rsp_build
36 :
37 : contains
38 :
39 0 : subroutine do1_rsp_build(s,ierr)
40 : ! call from other_rsp_build_model after changing params.
41 : ! can change rsp_* params; but cannot change nz or net.
42 : ! multiple calls are ok to search.
43 : use rsp_build, only: do_rsp_build
44 : use hydro_vars, only: set_vars, set_Teff
45 : type (star_info), pointer :: s
46 : integer, intent(out) :: ierr
47 : integer :: k
48 : include 'formats'
49 0 : call init_def(s)
50 0 : call init_for_rsp_eos_and_kap(s)
51 0 : s% rsp_period = 0d0
52 0 : call do_rsp_build(s,ierr)
53 0 : if (ierr /= 0) return
54 0 : do k=1,NZN
55 0 : s% v(k) = 0d0
56 0 : s% xh(s% i_v,k) = s% v(k)
57 : end do
58 : ierr = 0
59 0 : call finish_build_rsp_model(s,ierr)
60 0 : if (ierr /= 0) return
61 0 : s% doing_finish_load_model = .true.
62 0 : call set_vars(s, 0d0, ierr)
63 0 : if (ierr /= 0) return
64 0 : s% doing_finish_load_model = .false.
65 0 : call set_Teff(s, ierr)
66 0 : if (ierr /= 0) return
67 0 : end subroutine do1_rsp_build
68 :
69 :
70 0 : subroutine build_rsp_model(s,ierr)
71 : ! called by model_builder in place of loading a model
72 0 : use alloc, only: allocate_star_info_arrays
73 : use set_flags, only: set_RSP_flag
74 : use rsp_build, only: do_rsp_build
75 : use net, only: do_micro_change_net
76 : use const_def, only: Lsun, Rsun, Msun
77 : type (star_info), pointer :: s
78 : integer, intent(out) :: ierr
79 : include 'formats'
80 0 : NSTART = 1
81 0 : s% nz = s% RSP_nz
82 0 : if (s% job% change_initial_net) then
83 0 : call do_micro_change_net(s, s% job% new_net_name, ierr)
84 : else
85 0 : call do_micro_change_net(s, 'o18_and_ne22.net', ierr)
86 : end if
87 0 : if (ierr /= 0) then
88 0 : write(*,*) 'failed in do_micro_change_net'
89 0 : return
90 : end if
91 0 : s% tau_factor = s% RSP_surface_tau/s% tau_base
92 0 : call init_def(s)
93 0 : call init_allocate(s,s% nz)
94 0 : call allocate_star_info_arrays(s, ierr)
95 0 : if (ierr /= 0) then
96 0 : write(*,*) 'failed in allocate_star_info_arrays'
97 0 : return
98 : end if
99 0 : call set_RSP_flag(s% id, .true., ierr)
100 0 : if (ierr /= 0) then
101 0 : write(*,*) 'failed in set_RSP_flag'
102 0 : return
103 : end if
104 0 : call init_for_rsp_eos_and_kap(s)
105 0 : s% rsp_period = 0d0
106 0 : if (s% RSP_use_atm_grey_with_kap_for_Psurf) then
107 0 : s% tau_factor = s% RSP_tau_surf_for_atm_grey_with_kap/s% tau_base
108 0 : write(*,1) 'set s% tau_factor using RSP_tau_surf', s% tau_factor
109 : end if
110 0 : if (s% use_other_rsp_build_model) then
111 0 : call s% other_rsp_build_model(s% id,ierr)
112 0 : if (ierr /= 0) then
113 0 : write(*,*) 'failed in other_rsp_build_model'
114 0 : return
115 : end if
116 : else
117 0 : call do_rsp_build(s,ierr)
118 0 : if (ierr /= 0) then
119 0 : write(*,*) 'failed in do_rsp_build'
120 0 : return
121 : end if
122 : end if
123 0 : if (.not. s% use_RSP_new_start_scheme) &
124 0 : call set_random_velocities(s)
125 0 : rsp_tau_factor = s% tau_factor
126 0 : PERIODLIN = s% rsp_period
127 0 : if (s% rsp_period > 0d0) then
128 0 : s% rsp_dt = s% rsp_period/dble(s% rsp_target_steps_per_cycle)
129 : else
130 0 : s% rsp_dt = 1d0
131 : end if
132 0 : s% dt_next = s% rsp_dt
133 0 : s% dt = s% rsp_dt
134 0 : if (is_bad(s% dt)) then
135 0 : write(*,1) 'dt', s% dt
136 0 : return
137 : end if
138 0 : call finish_build_rsp_model(s, ierr)
139 0 : write(*,2) 'nz', s%nz
140 0 : write(*,1) 'T(nz)', s% T(s%nz)
141 0 : write(*,1) 'L_center/Lsun', s% L_center/Lsun
142 0 : write(*,1) 'R_center/Rsun', s% R_center/Rsun
143 0 : write(*,1) 'M_center/Msun', s% M_center/Msun
144 0 : write(*,1) 'L(1)/Lsun', s% L(1)/Lsun
145 0 : write(*,1) 'R(1)/Rsun', s% r(1)/Rsun
146 0 : write(*,1) 'M(1)/Msun', s% m(1)/Msun
147 0 : write(*,1) 'v(1)/1d5', s% v(1)/1d5
148 0 : write(*,1) 'tau_factor', s% tau_factor
149 0 : write(*,1) 'tau_base', s% tau_base
150 0 : write(*,*)
151 0 : end subroutine build_rsp_model
152 :
153 :
154 0 : subroutine finish_build_rsp_model(s,ierr)
155 : use star_utils, only: &
156 : normalize_dqs, set_qs, set_m_and_dm, set_dm_bar, set_m_grav_and_grav, &
157 0 : store_rho_in_xh, store_T_in_xh, store_r_in_xh
158 : type (star_info), pointer :: s
159 : integer, intent(out) :: ierr
160 : integer :: i, k, j
161 : include 'formats'
162 0 : do i=1,NZN
163 0 : k = NZN+1 - i
164 0 : s% Prad(k) = f_Edd_isotropic*s% erad(k)/s% Vol(k)
165 0 : if (is_bad(s% Prad(k))) then
166 0 : write(*,2) 's% Prad(k)', k, s% Prad(k)
167 0 : write(*,2) 's% erad(k)', k, s% erad(k)
168 0 : write(*,2) 's% Vol(k)', k, s% Vol(k)
169 0 : call mesa_error(__FILE__,__LINE__,'build_rsp_model')
170 : end if
171 0 : s% dq(k) = s% dm(k)/s% xmstar
172 0 : if (is_bad(s% Vol(k)) .or. s% Vol(k) <= 0d0) then
173 0 : write(*,2) 's% Vol(k)', I, s% Vol(k)
174 0 : call mesa_error(__FILE__,__LINE__,'build_rsp_model')
175 : end if
176 0 : call store_rho_in_xh(s, k, 1d0/s% Vol(k))
177 0 : call store_T_in_xh(s, k, s% T(k))
178 0 : call store_r_in_xh(s, k, s% r(k))
179 0 : s% xh(s% i_Et_RSP,k) = s% RSP_w(k)*s% RSP_w(k)
180 0 : do j=1,s% species
181 0 : s% xa(j,k) = xa(j)
182 : end do
183 0 : s% xh(s% i_v,k) = s% v(k)
184 : end do
185 0 : s% dq(s% nz) = (s% m(NZN) - s% M_center)/s% xmstar
186 0 : if (.not. s% do_normalize_dqs_as_part_of_set_qs) then
187 0 : call normalize_dqs(s, NZN, s% dq, ierr)
188 0 : if (ierr /= 0) then
189 0 : if (s% report_ierr) write(*,*) 'normalize_dqs failed in finish_build_rsp_model'
190 : return
191 : end if
192 : end if
193 0 : call set_qs(s, NZN, s% q, s% dq, ierr)
194 0 : if (ierr /= 0) then
195 0 : write(*,*) 'failed in set_qs'
196 0 : call mesa_error(__FILE__,__LINE__,'build_rsp_model')
197 : end if
198 0 : call set_m_and_dm(s)
199 0 : call set_m_grav_and_grav(s)
200 0 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
201 0 : end subroutine finish_build_rsp_model
202 :
203 :
204 0 : subroutine set_random_velocities(s)
205 0 : use star_utils, only: rand
206 : type (star_info), pointer :: s
207 : integer :: i, k
208 0 : if (s% RSP_Avel /= 0d0 .or. s% RSP_Arnd /= 0d0) then
209 0 : do i=1,NZN
210 0 : k = NZN+1-i
211 : s% v(k) = s% v(k) + &
212 0 : 1d5*dble(i)/dble(NZN)*(s% RSP_Avel + s% RSP_Arnd*(2d0*rand(s) - 1d0))
213 : end do
214 0 : write(*,*) 'set random velocities'
215 0 : s% RSP_have_set_velocities = .true.
216 : end if
217 0 : end subroutine set_random_velocities
218 :
219 :
220 0 : subroutine rsp_setup_part1(s,restart,ierr)
221 : ! called by finish_load_model before set_vars
222 0 : use const_def, only: crad
223 : use rsp_eval_eos_and_kap, only: &
224 : restart_rsp_eos_and_kap, get_surf_P_T_kap
225 : use alloc, only: resize_star_info_arrays
226 : use star_utils, only: get_XYZ
227 : type (star_info), pointer :: s
228 : logical, intent(in) :: restart
229 : integer, intent(out) :: ierr
230 0 : real(dp) :: tau_surf, kap_guess, T_surf, Psurf, kap_surf, Teff_atm, Y
231 : integer :: k
232 : include 'formats'
233 0 : ierr = 0
234 0 : if (s% RSP_use_atm_grey_with_kap_for_Psurf .and. &
235 : (s% atm_option /= 'T_tau' .or. &
236 : s% atm_T_tau_relation /= 'Eddington' .or. &
237 : s% atm_T_tau_opacity /= 'iterated')) then
238 0 : write(*,'(A)')
239 0 : write(*,*) 'when RSP_use_atm_grey_with_kap_for_Psurf, must set:'
240 0 : write(*,*) ' atm_option = ''T_tau'''
241 0 : write(*,*) ' atm_T_tau_relation = ''Eddington'''
242 0 : write(*,*) ' atm_T_tau_opacity = ''iterated'''
243 0 : ierr = -1
244 0 : call mesa_error(__FILE__,__LINE__,'rsp_setup_part1')
245 : return
246 : end if
247 0 : if (restart) then
248 0 : NSTART = 2
249 0 : call init_def(s)
250 0 : call restart_rsp_eos_and_kap(s)
251 : else
252 0 : prev_cycle_run_num_steps = 0
253 0 : run_num_iters_prev_period = 0
254 0 : run_num_retries_prev_period = 0
255 0 : NSTART = 1
256 0 : if (.not. s% job% create_RSP_model) then
257 0 : call init_def(s)
258 0 : call init_allocate(s,s% nz)
259 0 : call get_XYZ(s, s% xa(:,1), s% RSP_X, Y, s% RSP_Z)
260 0 : call init_for_rsp_eos_and_kap(s)
261 0 : IWORK=0
262 0 : NZN = s% nz
263 0 : ELSTA = s% L(1)
264 0 : RSTA = s% r(1)
265 0 : s% rsp_dt = s% dt_next
266 0 : if (s% max_timestep > 0d0 .and. s% rsp_dt > s% max_timestep) &
267 0 : s% rsp_dt = s% max_timestep
268 0 : rsp_tau_factor = s% tau_factor
269 0 : s% rsp_period = s% rsp_dt*dble(s% RSP_target_steps_per_cycle)
270 0 : s% RSP_have_set_velocities = .true.
271 0 : call copy_from_xh_to_rsp(s,-1)
272 0 : do k=1,NZN
273 0 : s% L_start(k) = 0d0
274 : end do
275 0 : if (s% RSP_use_atm_grey_with_kap_for_Psurf) then
276 0 : tau_surf = s% RSP_tau_surf_for_atm_grey_with_kap
277 0 : kap_guess = 1d-2
278 : call get_surf_P_T_kap(s, &
279 : s% m(1), s% r(1), s% L(1), tau_surf, kap_guess, &
280 0 : T_surf, Psurf, kap_surf, Teff_atm, ierr)
281 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed in get_surf_P_T_kap')
282 0 : else if (s% RSP_use_Prad_for_Psurf) then
283 0 : Psurf = crad*s% T(1)**4/3d0
284 : else
285 0 : Psurf = 0d0
286 : end if
287 0 : Psurf_from_atm = Psurf
288 : else
289 0 : s% dt_next = s% rsp_dt
290 0 : s% dt = s% rsp_dt
291 : end if
292 0 : rsp_min_dr_div_cs = 1d99
293 0 : rsp_min_rad_diff_time = 1d99
294 0 : call begin_calculation(s,restart,ierr)
295 : end if
296 0 : s% tau_factor = rsp_tau_factor
297 0 : end subroutine rsp_setup_part1
298 :
299 :
300 0 : subroutine rsp_setup_part2(s, restart, ierr)
301 0 : use hydro_vars, only: set_Teff
302 : use rsp_step
303 : ! called by finish_load_model after set_vars
304 : type (star_info), pointer :: s
305 : logical, intent(in) :: restart
306 : integer, intent(out) :: ierr
307 : integer :: nz
308 : include 'formats'
309 0 : if (restart) then
310 0 : call finish_rsp_photo_in(s)
311 0 : return
312 : end if
313 0 : ierr = 0
314 0 : nz = s% nz
315 0 : call finish_after_build_model(s)
316 0 : call copy_results(s)
317 0 : call set_Teff(s, ierr)
318 0 : if (ierr /= 0) then
319 0 : if (s% report_ierr) &
320 0 : write(*,*) 'rsp_setup_part2: set_Teff returned ierr', ierr
321 0 : return
322 : end if
323 0 : TEFF = s% Teff
324 0 : if (s% rsp_period > 0d0) then
325 0 : s% rsp_dt = s% rsp_period/dble(s% rsp_target_steps_per_cycle)
326 0 : s% dt_next = s% rsp_dt*s% RSP_initial_dt_factor
327 0 : s% dt = s% dt_next
328 : end if
329 0 : if (s% use_other_rsp_build_model .and. &
330 : s% set_RSP_Psurf_to_multiple_of_initial_P1 > 0d0) then
331 0 : s% RSP_Psurf = s% Peos(1)*s% set_RSP_Psurf_to_multiple_of_initial_P1
332 0 : write(*,1) 'rsp_setup_part2 set RSP_Psurf', s% RSP_Psurf
333 : end if
334 0 : end subroutine rsp_setup_part2
335 :
336 :
337 0 : subroutine get_LINA_info(s,ierr)
338 0 : use rsp_lina, only: do_LINA
339 : type (star_info), pointer :: s
340 : integer, intent(out) :: ierr
341 :
342 0 : real(dp), allocatable :: VEL(:,:)
343 : real(dp), allocatable, dimension(:) :: &
344 0 : M, DM, DM_BAR, R, Vol, T, RSP_Et, Lr
345 : integer :: NMODES, I, k, sz
346 0 : real(dp) :: amix1, amix2, velkm
347 :
348 : include 'formats'
349 0 : ierr = 0
350 :
351 0 : if (s% RSP_kick_vsurf_km_per_sec == 0d0) then
352 0 : write(*,*) 'skip calling LINA since RSP_kick_vsurf_km_per_sec = 0'
353 0 : return
354 : end if
355 :
356 0 : sz = NZN+1
357 :
358 : allocate(VEL(sz,15), &
359 0 : M(sz), DM(sz), DM_BAR(sz), R(sz), Vol(sz), T(sz), RSP_Et(sz), Lr(sz))
360 :
361 0 : do i=1,NZN
362 0 : k = NZN+1-i
363 0 : M(i) = s% m(k)
364 0 : DM(i) = s% dm(k)
365 0 : DM_BAR(i) = s% dm_bar(k)
366 0 : R(i) = s% r(k)
367 0 : Vol(i) = s% Vol(k)
368 0 : T(i) = s% T(k)
369 0 : RSP_Et(i) = s% RSP_Et(k)
370 0 : Lr(i) = 4d0*pi*s% r(k)**2*s% Fr(k)
371 : end do
372 :
373 0 : NMODES = s% RSP_nmodes
374 :
375 : call do_LINA(s, s% RSP_L*SUNL, NZN, NMODES, VEL, &
376 : s% rsp_LINA_periods, s% rsp_LINA_growth_rates, &
377 0 : M, DM, DM_BAR, R, Vol, T, RSP_Et, Lr, ierr)
378 0 : if (ierr /= 0) return
379 :
380 0 : write(*,'(a)') ' P(days) growth'
381 0 : do I=1,NMODES
382 0 : write(*,'(I3,2X,99e16.5)') I-1, &
383 0 : s% rsp_LINA_periods(I)/86400.d0, &
384 0 : s% rsp_LINA_growth_rates(I)
385 : end do
386 :
387 : s% rsp_period = &
388 0 : s% rsp_LINA_periods(s% RSP_mode_for_setting_PERIODLIN + 1)
389 :
390 0 : amix1 = s% RSP_fraction_1st_overtone
391 0 : amix2 = s% RSP_fraction_2nd_overtone
392 0 : velkm = s% RSP_kick_vsurf_km_per_sec
393 0 : s% v_center = 0d0
394 0 : do i=1,NZN
395 0 : k = NZN+1-i
396 : s% v(k)=1.0d5*VELKM* &
397 0 : ((1.0d0-AMIX1-AMIX2)*vel(i,1)+AMIX1*vel(i,2)+AMIX2*vel(i,3))
398 : end do
399 :
400 0 : s% RSP_have_set_velocities = .true.
401 :
402 0 : end subroutine get_LINA_info
403 :
404 :
405 0 : subroutine rsp_one_step(s,ierr)
406 0 : use brunt, only: do_brunt_N2
407 : use rsp_step, only: rsp_set_Teff, &
408 : turn_off_time_weighting, turn_on_time_weighting
409 : type (star_info), pointer :: s
410 : integer, intent(out) :: ierr
411 : integer :: k, j, k_max_abs_rel_hse_err
412 0 : real(dp) :: hse_err, max_abs_rel_hse_err
413 : logical :: restart
414 :
415 : include 'formats'
416 :
417 0 : ierr = 0
418 0 : s% RSP_just_set_velocities = .false.
419 0 : if (.not. s% RSP_have_set_velocities) then
420 :
421 0 : max_abs_rel_hse_err = 0d0
422 0 : k_max_abs_rel_hse_err = 0
423 0 : do k=2,s% nz
424 : hse_err = &
425 0 : (s% Peos(k-1) - s% Peos(k))/(-s% cgrav(k)*s% m(k)*s% dm_bar(k)/(4d0*pi*s% r(k)**4)) - 1d0
426 0 : if (abs(hse_err) >= max_abs_rel_hse_err) then
427 0 : max_abs_rel_hse_err = abs(hse_err)
428 0 : k_max_abs_rel_hse_err = k
429 : end if
430 : end do
431 :
432 0 : s% need_to_save_profiles_now = .true.
433 0 : s% RSP_just_set_velocities = .true.
434 :
435 0 : write(*,3) 'relaxation max_abs_rel_hse_err, days', s% model_number, k_max_abs_rel_hse_err, &
436 0 : max_abs_rel_hse_err, s% time/(24*3600)
437 :
438 0 : if (.not. s% use_other_RSP_linear_analysis) then
439 0 : call get_LINA_info(s,ierr)
440 : else
441 : ! must set gradT before calling since gyre needs it.
442 : ! Y_face is superadiabatic gradient
443 0 : do k=1,NZN
444 0 : s% gradT_sub_grada(k) = s% Y_face(k)
445 0 : if (k > 1) s% gradT(k) = &
446 0 : s% Y_face(k) + 0.5d0*(s% grada(k-1) + s% grada(k))
447 0 : s% abar(k) = abar
448 0 : s% zbar(k) = zbar
449 0 : s% X(k) = X
450 0 : s% Z(k) = Z
451 0 : do j=1,s% species
452 0 : s% xa(j,k) = xa(j)
453 : end do
454 : end do
455 0 : s% gradT(1) = s% gradT(2)
456 0 : s% calculate_Brunt_N2 = .true.
457 0 : call do_brunt_N2(s, 1, NZN, ierr)
458 0 : if (ierr /= 0) return
459 0 : restart = .false.
460 0 : call s% other_rsp_linear_analysis(s% id, restart, ierr)
461 0 : if (ierr /= 0) then
462 0 : if (s% report_ierr) &
463 0 : write(*,*) 'other_rsp_linear_analysis ierr', ierr
464 0 : return
465 : end if
466 0 : s% RSP_have_set_velocities = .true.
467 : end if
468 :
469 0 : PERIODLIN = s% rsp_period
470 0 : s% rsp_dt = s% rsp_period/dble(s% rsp_target_steps_per_cycle)
471 0 : s% dt = s% rsp_dt
472 :
473 0 : s% cumulative_energy_error_old = 0d0
474 0 : s% time = 0d0
475 0 : s% time_old = 0d0
476 0 : write(*,*) 'automatically resets age and cumulative energy error info when sets velocities'
477 0 : s% need_to_save_profiles_now = .true.
478 :
479 0 : call set_random_velocities(s)
480 :
481 : end if
482 :
483 0 : s% do_history_file = s% RSP_have_set_velocities ! don't write history entries until set velocities
484 : !call turn_on_time_weighting(s)
485 :
486 0 : if (s% dt > s% RSP_max_dt .and. s% RSP_max_dt > 0d0) then
487 0 : s% dt = s% RSP_max_dt
488 : end if
489 :
490 0 : call do1_step(s,ierr)
491 0 : if (ierr /= 0) return
492 :
493 0 : call copy_results(s)
494 0 : call rsp_set_Teff(s)
495 0 : if (s% RSP_write_map) call add_to_map
496 :
497 :
498 : contains
499 :
500 0 : subroutine add_to_map
501 0 : use profile_getval, only: get_profile_val
502 : integer :: i, k, NPCH1, NPCH2, IP, n, io
503 0 : real(dp) :: FASE
504 : character (len=256) :: fname
505 : include 'formats'
506 0 : NPCH1 = s% RSP_map_first_period
507 0 : NPCH2 = s% RSP_map_last_period
508 0 : IP = s% RSP_num_periods
509 0 : io = 77
510 0 : if (IP+1>=NPCH1.and.IP+1<=NPCH2) then
511 0 : if(.not. writing_map) then
512 0 : call read_map_specs(s,ierr)
513 0 : if (ierr /= 0) then
514 0 : write(*,*) 'failed in read_map_specs'
515 : return
516 : end if
517 0 : if (len_trim(s% RSP_map_filename) == 0) &
518 0 : s% RSP_map_filename = 'map.data'
519 0 : fname = trim(s% log_directory) // '/' // trim(s% RSP_map_filename)
520 0 : open(io,file=trim(fname),status='unknown')
521 0 : write(*,*) 'writing ' // trim(fname)
522 0 : write(io,'(i18,1x,i4)',advance='no') 1, 2
523 0 : do n=1,num_map_cols
524 0 : write(io,'(1x,i18)',advance='no') n+2
525 : end do
526 0 : write(io,'(A)')
527 0 : write(io,'(a18,1x,a4)',advance='no') 'phase', 'zone'
528 0 : do n=1,num_map_cols
529 0 : write(io,'(1x,a18)',advance='no') trim(map_col_names(n))
530 : end do
531 0 : write(io,'(A)')
532 0 : writing_map = .true.
533 0 : done_writing_map = .false.
534 0 : s% need_to_set_history_names_etc = .true.
535 0 : s% star_history_name = s% RSP_map_history_filename
536 0 : FASE0 = s% time
537 : end if
538 0 : FASE=(s% time-FASE0)/s% rsp_period
539 : !write(*,4) 'add to map', s% model_number, IP, NPCH2, FASE
540 0 : do k=1,NZN,s% RSP_map_zone_interval ! gnuplot pm3d map
541 0 : I = NZN+1 - k
542 0 : if(I>IBOTOM.and.I<NZN) then
543 0 : write(io,'(d18.10,1x,i4)',advance='no') FASE, k
544 0 : do n=1,num_map_cols
545 : write(io,'(1x,d18.10)',advance='no') &
546 0 : get_profile_val(s,map_ids(n),k)
547 : end do
548 0 : write(io,'(A)')
549 : !write(io,778) FASE,I,s% T(k), &
550 : ! s% Hp_face(k),s% Y_face(k),s% PII(k),s% Lc(k)/s% L(k), &
551 : ! 4d0*pi*s% r(k)**2*s% Fr(k)/s% L(k),s% Lt(k)/s% L(k), &
552 : ! s% RSP_w(k)**2,s% egas(k)+s% erad(k),s% csound(k), &
553 : ! s% Pt(k),s% Pgas(k)+s% Prad(k),s% Eq(k), &
554 : ! s% COUPL(k)
555 : end if
556 : end do
557 : !write(io,*)
558 : end if
559 0 : if(IP==NPCH2 .and. .not. done_writing_map) then
560 0 : close(io)
561 0 : fname = trim(s% log_directory) // '/' // trim(s% RSP_map_filename)
562 0 : write(*,*) ' close ' // trim(fname)
563 0 : done_writing_map = .true.
564 : end if
565 : 778 format(d16.10,1x,i3,14(1x,d16.10))
566 0 : end subroutine add_to_map
567 :
568 : end subroutine rsp_one_step
569 :
570 :
571 0 : subroutine read_map_specs(s,ierr)
572 : use utils_lib
573 : use utils_def
574 : use profile_getval, only: get_profile_id
575 : type (star_info), pointer :: s
576 : integer, intent(out) :: ierr
577 :
578 : integer :: iounit, n, i, t, id, col
579 : character (len=strlen) :: buffer, string, filename
580 :
581 : include 'formats'
582 :
583 0 : ierr = 0
584 :
585 0 : filename = s% RSP_map_columns_filename
586 0 : if (len_trim(filename) == 0) filename = 'map_columns.list'
587 0 : open(newunit=iounit, file=trim(filename), action='read', status='old', iostat=ierr)
588 0 : if (ierr /= 0) then
589 0 : write(*,*) 'failed to open ' // trim(filename)
590 0 : return
591 : end if
592 :
593 0 : n = 0
594 0 : i = 0
595 0 : col = 0
596 :
597 0 : num_map_cols = 0
598 :
599 0 : do
600 0 : t = token(iounit, n, i, buffer, string)
601 0 : if (t == eof_token) exit
602 0 : if (t /= name_token) then
603 0 : write(*,*) 'error reading map columns list ' // trim(filename)
604 0 : call mesa_error(__FILE__,__LINE__,'read_map_specs')
605 0 : call error; return
606 : end if
607 0 : id = get_profile_id(s, string)
608 0 : if (id < 0) then
609 0 : write(*,*) 'error: <' // trim(string) // '> in map columns is not in your profile list'
610 0 : write(*,*) 'please add it to your profile columns list and try again'
611 0 : write(*,*) 'also, replace any TAB characters by spaces in ' // trim(filename)
612 0 : call mesa_error(__FILE__,__LINE__,'read_map_specs')
613 0 : call error; return
614 : end if
615 0 : col = col+1
616 0 : if (col > max_map_cols) then
617 0 : write(*,*) 'error: ' // trim(filename) // ' has too many map columns'
618 0 : write(*,*) 'the max is currently fixed at ', max_map_cols
619 0 : write(*,*) 'please delete some and try again'
620 0 : call error; return
621 : end if
622 0 : map_col_names(col) = trim(string)
623 0 : map_ids(col) = id
624 : end do
625 :
626 0 : num_map_cols = col
627 :
628 0 : close(iounit)
629 :
630 : contains
631 :
632 0 : subroutine error
633 0 : ierr = -1
634 0 : close(iounit)
635 0 : end subroutine error
636 :
637 : end subroutine read_map_specs
638 :
639 :
640 0 : subroutine do1_step(s,ierr)
641 : type (star_info), pointer :: s
642 : integer, intent(out) :: ierr
643 : integer :: i, k
644 0 : real(dp) :: dr_div_cs, r_in, r_00, max_dt, target_dt, total_radiation
645 :
646 : include 'formats'
647 :
648 0 : ID=ID+1
649 :
650 : target_dt = min( &
651 : s% rsp_period/dble(s% RSP_target_steps_per_cycle), &
652 0 : s% dt*s% max_timestep_factor)
653 0 : if (s% rsp_max_dt > 0) target_dt = min(target_dt, s% rsp_max_dt)
654 0 : if (s% max_timestep > 0) target_dt = min(target_dt, s% max_timestep)
655 0 : s% dt = target_dt
656 :
657 0 : if (is_bad(s% dt)) then
658 0 : write(*,1) 'do1_step dt', s% dt
659 0 : write(*,1) 'rsp_period', s% rsp_period
660 0 : write(*,2) 'RSP_target_steps_per_cycle', s% RSP_target_steps_per_cycle
661 0 : write(*,1) 'max_timestep_factor', s% max_timestep_factor
662 0 : write(*,1) 'rsp_period/RSP_target_steps_per_cycle/', s% rsp_period/s% RSP_target_steps_per_cycle
663 0 : call mesa_error(__FILE__,__LINE__,'do1_step 1')
664 : end if
665 :
666 0 : max_dt = rsp_min_dr_div_cs*s% RSP_max_dt_times_min_dr_div_cs
667 : !if (s% RSP_max_dt_times_min_rad_diff_time > 0d0 .and. rsp_min_rad_diff_time > 0d0) then
668 : ! if (rsp_min_rad_diff_time*s% RSP_max_dt_times_min_rad_diff_time < max_dt) then
669 : ! max_dt = rsp_min_rad_diff_time*s% RSP_max_dt_times_min_rad_diff_time
670 : ! if (s% dt > max_dt) then
671 : ! write(*,3) 'dt limited by rad diff time', NZN+1-i_min_rad_diff_time, s% model_number, &
672 : ! s% dt, rsp_min_rad_diff_time, s% RSP_max_dt_times_min_rad_diff_time
673 : ! !call mesa_error(__FILE__,__LINE__,'rsp')
674 : ! end if
675 : ! end if
676 : !end if
677 0 : if (s% dt > max_dt) then
678 0 : if (s% RSP_report_limit_dt) then
679 0 : write(*,4) 'limit to RSP_max_dt_times_min_dr_div_cs', s% model_number, max_dt, s% dt
680 : end if
681 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_step 1')
682 0 : s% dt = max_dt
683 : end if
684 :
685 0 : if (s% force_timestep > 0d0) s% dt = s% force_timestep ! overrides everything else
686 :
687 0 : if (is_bad(s% dt) .or. s% dt <= 0d0) then
688 0 : write(*,1) 'dt', s% dt
689 0 : write(*,1) 'RSP_max_dt_times_min_dr_div_cs', s% RSP_max_dt_times_min_dr_div_cs
690 0 : write(*,1) 'rsp_min_dr_div_cs', rsp_min_dr_div_cs
691 0 : write(*,1) 'rsp_min_rad_diff_time', rsp_min_rad_diff_time
692 0 : call mesa_error(__FILE__,__LINE__,'do1_step 2')
693 : end if
694 :
695 : ierr = 0
696 0 : call HYD(s,ierr)
697 0 : if (ierr /= 0) return
698 : ! s% dt might have been reduced by retries in HYD
699 0 : s% time = s% time_old + s% dt
700 0 : s% rsp_dt = s% dt ! will be used to set dt for next step
701 :
702 : ! set this here for use in next step. to avoid restart problems.
703 0 : rsp_min_dr_div_cs = 1d99
704 0 : i_min_dr_div_cs = -1
705 0 : r_00 = s% R_center
706 0 : do i = 1,nzn
707 0 : k = NZN+1-i
708 0 : r_in = r_00
709 0 : r_00 = s% r(k)
710 0 : if (abs(s% v(k))/s% csound(k) < s% RSP_v_div_cs_threshold_for_dt_limit) cycle
711 0 : k = nzn+1-i
712 0 : dr_div_cs = (r_00 - r_in)/s% csound(k)
713 0 : if (dr_div_cs < rsp_min_dr_div_cs) then
714 0 : rsp_min_dr_div_cs = dr_div_cs
715 0 : i_min_dr_div_cs = i
716 : end if
717 : end do
718 :
719 0 : rsp_min_rad_diff_time = 1d99
720 0 : i_min_rad_diff_time = -1
721 0 : if (s% RSP_max_dt_times_min_rad_diff_time > 0d0) then
722 0 : rsp_min_rad_diff_time = dt_for_radiative_diffusion(i_min_rad_diff_time)
723 : end if
724 :
725 0 : call calculate_work_integrals(s)
726 0 : call calculate_energies(s,total_radiation)
727 0 : call gather_pulse_statistics(s)
728 0 : if (s% RSP_max_num_periods < 0 .or. &
729 : s% rsp_num_periods < s% RSP_max_num_periods) return
730 0 : call get_GRPDV(s)
731 :
732 : contains
733 :
734 0 : real(dp) function dt_for_radiative_diffusion(i_min_rad_diff_time)
735 : integer, intent(out) :: i_min_rad_diff_time
736 0 : real(dp) :: min_dt, dt_cell, l, D, dr
737 : integer :: k, k_min_dt, nz
738 : include 'formats'
739 0 : min_dt = 1d99
740 0 : nz = s% nz
741 0 : k_min_dt = -1
742 : i_min_rad_diff_time = -1
743 0 : do k=1,nz
744 0 : l = s% V(k)/s% opacity(k) ! photon mean free path
745 0 : D = clight*l/3d0 ! diffusion coefficient, clight/(3*opacity*rho)
746 0 : if (k < nz) then
747 0 : dr = s% r(k) - s% r(k+1)
748 : else
749 0 : dr = s% r(k) - s% r_center
750 : end if
751 : ! if curious, ask Lars about the Pgas/Prad term
752 0 : dt_cell = dr**2/D*(1d0 + s% Pgas(k)/s% Prad(k))
753 0 : if (dt_cell < min_dt) then
754 0 : min_dt = dt_cell
755 0 : k_min_dt = k
756 : end if
757 : end do
758 0 : i_min_rad_diff_time = NZN-k_min_dt+1
759 0 : dt_for_radiative_diffusion = min_dt
760 0 : end function dt_for_radiative_diffusion
761 :
762 : end subroutine do1_step
763 :
764 :
765 0 : subroutine gather_pulse_statistics(s) ! assumes have set EKMAX and EKMIN
766 : ! updates LMAX, LMIN, RMAX, RMIN,
767 : ! s% rsp_GREKM, s% rsp_GREKM_avg_abs, s% rsp_DeltaR, s% rsp_DeltaMAG
768 : type (star_info), pointer :: s
769 : logical :: cycle_complete
770 : include 'formats'
771 0 : if(s% L(1)/SUNL>LMAX) LMAX=s% L(1)/SUNL
772 0 : if(s% L(1)/SUNL<LMIN) LMIN=s% L(1)/SUNL
773 0 : INSIDE=0
774 0 : call check_cycle_completed(s,cycle_complete)
775 0 : ULL=UN
776 0 : TE_start=s% time
777 0 : if (cycle_complete) then
778 0 : if (s% rsp_num_periods > 1) then
779 0 : s% rsp_GREKM = (EKMAX-EKMAXL)/(EKMAX+EKMAXL)*2.d0
780 : s% rsp_GREKM_avg_abs = s% rsp_GREKM_avg_abs_frac_new*abs(s% rsp_GREKM) + &
781 0 : (1d0 - s% rsp_GREKM_avg_abs_frac_new)*s% rsp_GREKM_avg_abs
782 : end if
783 0 : EKDEL = EKMAX-EKMIN
784 0 : EKMAXL = EKMAX
785 0 : EKMAX =-10.d50
786 0 : EKMIN = -EKMAX
787 0 : s% rsp_DeltaR = RMAX-RMIN
788 0 : RMAX = -SUNR
789 0 : RMIN = -RMAX
790 0 : s% rsp_DeltaMAG = 2.5d0*log10(LMAX/LMIN)
791 0 : LMIN = 10.d10
792 0 : LMAX = -LMIN
793 0 : VMAX = 0
794 : end if
795 0 : end subroutine gather_pulse_statistics
796 :
797 :
798 0 : subroutine get_GRPDV(s)
799 : type (star_info), pointer :: s
800 : integer :: I, k
801 0 : PDVWORK=0.d0
802 0 : do I=1,NZN
803 0 : k = NZN+1-i
804 : WORK(I)= WORK(I)+(VV0(I)-s% Vol(k))*s% dm(k)* &
805 : (THETA*(PPP0(I)+PPQ0(I))+ &
806 : THETA1*((s% Pgas(k)+s% Prad(k))+s% Pvsc(k))) &
807 0 : -s% dt*s% dm(k)*s% Eq(k)
808 : WORKQ(I)= WORKQ(I)+(VV0(I)-s% Vol(k))*s% dm(k)* &
809 0 : (THETA*PPQ0(I)+THETA1*s% Pvsc(k))
810 0 : PDVWORK=PDVWORK+WORK(i)
811 : end do
812 0 : s% rsp_GRPDV=PDVWORK/EKDEL
813 0 : if (is_bad(s% rsp_GRPDV)) s% rsp_GRPDV=0d0
814 0 : end subroutine get_GRPDV
815 :
816 :
817 0 : subroutine begin_calculation(s,restart,ierr)
818 : type (star_info), pointer :: s
819 : logical, intent(in) :: restart
820 : integer, intent(out) :: ierr
821 : real(dp) :: total_radiation
822 : include 'formats'
823 : ierr = 0
824 0 : FIRST = 0
825 0 : TT1 = 0.d0
826 0 : EKMAX = -10.d50
827 0 : EKMIN = -EKMAX
828 0 : EKMAXL = EKMAX
829 0 : s% rsp_num_periods =-1
830 0 : ID = 0 !number of timesteps done in one period
831 0 : INSIDE = 1 ! for initial call
832 : !s% mstar = M(1)
833 0 : call set_star_vars(s,ierr)
834 0 : if(s% rsp_num_periods==1)s% rsp_GREKM=0.d0
835 0 : EKDEL = EKMAX-EKMIN
836 0 : EKMAXL = EKMAX
837 0 : EKMAX =-10.d50
838 0 : EKMIN =-EKMAX
839 0 : LMIN = 10.d10
840 0 : LMAX =-LMIN
841 0 : RMAX = -SUNR
842 0 : RMIN = -RMAX
843 0 : VMAX = 0
844 0 : s% rsp_GREKM = 0
845 0 : s% rsp_GREKM_avg_abs = 0
846 0 : s% rsp_GRPDV = 0
847 0 : s% rsp_DeltaR = 0
848 0 : s% rsp_DeltaMAG = 0
849 0 : ID = 0
850 0 : E0 = 0.d0
851 0 : call init_HYD(s,ierr)
852 0 : if (ierr /= 0) return
853 0 : s% rsp_num_periods = 0
854 0 : call calculate_energies(s,total_radiation)
855 0 : E0 = EDE_start
856 0 : call calculate_work_integrals(s)
857 : end subroutine begin_calculation
858 :
859 :
860 0 : subroutine calculate_work_integrals(s)
861 : type (star_info), pointer :: s
862 : integer :: i, k
863 0 : real(dp) :: dt, dm, dVol, P_tw, Pvsc_tw, Ptrb_tw
864 : character (len=256) :: fname
865 0 : dt = s% dt
866 : ! LAST STEP OF PdV
867 0 : if(INSIDE==1.and.IWORK==1) then
868 0 : IWORK=0
869 0 : do I=1,NZN
870 0 : k = NZN+1-i
871 0 : dm = s% dm(k)
872 0 : dVol = VV0(I) - s% Vol_start(k)
873 : P_tw = THETA*PPP0(I) + &
874 0 : THETA1*(s% Pgas_start(k) + s% Prad_start(k))
875 0 : Pvsc_tw = THETA*PPQ0(I) + THETA1*s% Pvsc_start(k)
876 0 : Ptrb_tw = THETAT*PPT0(I) + THETAT1*s% Ptrb_start(k)
877 : WORK(I) = WORK(I) + &
878 : dVol*s% dm(k)*(P_tw + Pvsc_tw + Ptrb_tw) &
879 0 : - dt*dm*s% Eq(k)
880 0 : WORKQ(I) = WORKQ(I) + dVol*dm*Pvsc_tw
881 0 : WORKT(I) = WORKT(I) + dVol*dm*Ptrb_tw
882 0 : WORKC(I) = WORKC(I) - dt*dm*s% Eq(k)
883 : end do
884 0 : if (s% rsp_num_periods == s% RSP_work_period) then
885 0 : fname = trim(s% log_directory) // '/' // trim(s% RSP_work_filename)
886 0 : write(*,*) 'write work integral data to ' // trim(fname)
887 0 : open(78,file=trim(fname),status='unknown')
888 0 : do I=1,NZN
889 0 : k = NZN+1-i
890 : write(78,'(I3,tr1,F7.4,4(tr1,d16.10))') &
891 0 : I, log10(sum(s% dm(1:k))), &
892 0 : WORK(I)/EKDEL, WORKQ(I)/EKDEL, &
893 0 : WORKT(I)/EKDEL, WORKC(I)/EKDEL
894 : end do
895 0 : close(78)
896 : end if
897 0 : PDVWORK=0.d0
898 0 : do I=1,NZN
899 0 : k = NZN+1-i
900 0 : PDVWORK=PDVWORK+WORK(I)
901 0 : WORK(I)=0.d0
902 0 : WORKQ(I)=0.d0
903 0 : WORKT(I)=0.d0
904 0 : WORKC(I)=0.d0
905 : end do
906 0 : s% rsp_GRPDV=PDVWORK/EKDEL
907 : end if
908 :
909 : ! INITIAL STEP OF PdV:
910 0 : if((INSIDE==1.and.IWORK==0).or. &
911 : (s% rsp_num_periods==0.and.IWORK==0))then
912 0 : IWORK=1
913 0 : do I=1,NZN
914 0 : k = NZN+1-i
915 0 : VV0(I) = s% Vol_start(k)
916 0 : PPP0(I) = s% Pgas_start(k) + s% Prad_start(k)
917 0 : PPQ0(I) = s% Pvsc_start(k)
918 0 : PPT0(I) = s% Ptrb_start(k)
919 0 : PPC0(I) = s% Chi_start(k)
920 : end do
921 : end if
922 :
923 : ! FIRST AND NEXT STEPS of PdV:
924 0 : if(IWORK==1)then
925 0 : do I=1,NZN
926 0 : k = NZN+1-i
927 0 : dm = s% dm(k)
928 0 : dVol = s% Vol(k) - s% Vol_start(k)
929 : P_tw = THETA*(s% Pgas(k) + s% Prad(k)) &
930 0 : + THETA1*(s% Pgas_start(k) + s% Prad_start(k))
931 0 : Pvsc_tw = THETA*s% Pvsc(k) + THETA1*s% Pvsc_start(k)
932 0 : Ptrb_tw = THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k)
933 : WORK(I) = WORK(I) + &
934 0 : dm*(dVol*(P_tw + Pvsc_tw + Ptrb_tw) - dt*s% Eq(k))
935 0 : WORKQ(I)= WORKQ(I) + dm*dVol*Pvsc_tw
936 0 : WORKT(I)= WORKT(I) + dm*dVol*Ptrb_tw
937 0 : WORKC(I)= WORKC(I) - dt*dm*s% Eq(k)
938 : end do
939 : end if
940 0 : end subroutine calculate_work_integrals
941 :
942 :
943 0 : subroutine rsp_total_energy_integrals(s, &
944 : total_internal_energy, total_gravitational_energy, &
945 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
946 : total_turbulent_energy, sum_total, total_radiation)
947 : type (star_info), pointer :: s
948 : real(dp), intent(out) :: &
949 : total_internal_energy, total_gravitational_energy, &
950 : total_radial_kinetic_energy, total_rotational_kinetic_energy, &
951 : total_turbulent_energy, sum_total, total_radiation
952 : include 'formats'
953 0 : call calculate_energies(s,total_radiation)
954 0 : total_internal_energy = ETHE
955 0 : total_gravitational_energy = EGRV
956 0 : total_radial_kinetic_energy = EKIN
957 0 : total_rotational_kinetic_energy = 0d0
958 0 : total_turbulent_energy = ECON
959 0 : sum_total = ETOT
960 0 : end subroutine rsp_total_energy_integrals
961 :
962 :
963 0 : subroutine check_cycle_completed(s,cycle_complete)
964 : ! uses ULL, FIRST, s% RSP_min_max_R_for_periods, PERIODLIN, s% RSP_min_PERIOD_div_PERIODLIN
965 : ! depends on s% RSP_have_set_velocities = .true.
966 : ! sets s% rsp_num_periods, s% rsp_period
967 :
968 : type (star_info), pointer :: s
969 : logical, intent(out) :: cycle_complete
970 0 : real(dp) :: TET, min_PERIOD
971 : include 'formats'
972 0 : TET = s% time
973 0 : cycle_complete = .false.
974 0 : UN=s% v(1)
975 0 : if(UN>0.d0.and.ULL<=0.d0) then
976 0 : RMIN=s% r(1)/SUNR
977 : end if
978 0 : if (s% model_number==1) return
979 0 : if (.not. s% RSP_have_set_velocities) return
980 0 : if (s% r(1)/SUNR < s% RSP_min_max_R_for_periods) return
981 0 : if (UN/s% csound(1) > VMAX) then
982 0 : VMAX = UN/s% csound(1)
983 : end if
984 0 : if(UN*ULL>0.0d0.or.UN>0.d0) return
985 0 : T0=TET
986 0 : min_PERIOD = PERIODLIN*s% RSP_min_PERIOD_div_PERIODLIN
987 0 : if (abs(UN-ULL)>1.0d-10) T0=TE_start-(TE_start-TET)*ULL/(ULL-UN)
988 0 : if (min_PERIOD > 0d0 .and. T0-TT1 < min_PERIOD) return
989 0 : if (s% r(1)/SUNR - RMIN < s% RSP_min_deltaR_for_periods) return
990 0 : if(FIRST==1)then
991 0 : cycle_complete = .true.
992 0 : s% rsp_num_periods=s% rsp_num_periods+1
993 0 : s% rsp_period=T0-TT1
994 0 : if (is_bad(s% rsp_period)) then
995 0 : write(*,2) 'PERIOD', s% model_number, s% rsp_period
996 0 : write(*,2) 'TET', s% model_number, TET
997 0 : write(*,2) 'T0', s% model_number, T0
998 0 : write(*,2) 'TT1', s% model_number, TT1
999 0 : call mesa_error(__FILE__,__LINE__,'check_cycle_completed')
1000 : end if
1001 0 : RMAX=s% r(1)/SUNR !(NOT INTERPOLATED)
1002 : write(*,'(a7,i7,f11.5,a9,f11.5,a14,f9.5,a9,i3,a7,i6,a16,f9.5,a6,i10,a6,f10.3)') &
1003 0 : 'period', s% rsp_num_periods, s% rsp_period/(24*3600), &
1004 0 : 'delta R', RMAX - RMIN, &
1005 0 : 'max vsurf/cs', VMAX, &
1006 0 : 'retries', &
1007 0 : s% num_retries - run_num_retries_prev_period, &
1008 0 : 'steps', &
1009 0 : s% model_number - prev_cycle_run_num_steps, &
1010 0 : 'iters/step', &
1011 : dble(s% total_num_solver_iterations - run_num_iters_prev_period)/ &
1012 0 : dble(s% model_number - prev_cycle_run_num_steps), &
1013 0 : 'model', s% model_number, 'days', s% time/(24*3600)
1014 0 : prev_cycle_run_num_steps = s% model_number
1015 0 : run_num_iters_prev_period = s% total_num_solver_iterations
1016 0 : run_num_retries_prev_period = s% num_retries
1017 0 : TT1=T0
1018 0 : INSIDE=1
1019 0 : VMAX = 0d0
1020 : else
1021 0 : write(*,*) 'first maximum radius, period calculations start at model, day', &
1022 0 : s% model_number, s% time/(24*3600)
1023 0 : TT1=T0
1024 0 : FIRST=1
1025 0 : ID=0
1026 : end if
1027 : end subroutine check_cycle_completed
1028 :
1029 : end module rsp
|
