Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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 evolve_support
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 :
25 : implicit none
26 :
27 : private
28 : public :: set_current_to_old, new_generation, output, output_to_file
29 :
30 : contains
31 :
32 23 : subroutine new_generation(s, ierr)
33 : use utils_lib
34 : type (star_info), pointer :: s
35 : integer, intent(out) :: ierr
36 :
37 : integer :: nz, j, k
38 :
39 : include 'formats'
40 :
41 11 : ierr = 0
42 11 : nz = s% nz
43 :
44 11 : if (.not. s% rsp_flag) then
45 :
46 11 : call copy_to_old(s% dq, s% dq_old, ierr)
47 11 : if (ierr /= 0) return
48 :
49 11 : call copy_to_old(s% q, s% q_old, ierr)
50 11 : if (ierr /= 0) return
51 :
52 11 : call copy_to_old(s% omega, s% omega_old, ierr)
53 11 : if (ierr /= 0) return
54 :
55 11 : call copy_to_old(s% j_rot, s% j_rot_old, ierr)
56 11 : if (ierr /= 0) return
57 :
58 11 : call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr)
59 11 : if (ierr /= 0) return
60 :
61 11 : call enlarge_if_needed_2(s% xh_old,s% nvar_hydro,nz,nz_alloc_extra,ierr)
62 11 : if (ierr /= 0) return
63 11 : if (s% fill_arrays_with_NaNs) call fill_with_NaNs_2d(s% xh_old)
64 :
65 11 : call enlarge_if_needed_2(s% xa_old,s% species,nz,nz_alloc_extra,ierr)
66 11 : if (ierr /= 0) return
67 11 : if (s% fill_arrays_with_NaNs) call fill_with_NaNs_2d(s% xh_old)
68 :
69 13098 : do k = 1, s% nz
70 65435 : do j=1, s% nvar_hydro
71 65435 : s% xh_old(j,k) = s% xh(j,k)
72 : end do
73 117794 : do j=1,s% species
74 117783 : s% xa_old(j,k) = s% xa(j,k)
75 : end do
76 : end do
77 :
78 : end if
79 :
80 11 : s% model_number_old = s% model_number
81 11 : s% nz_old = s% nz
82 11 : s% time_old = s% time
83 11 : s% dt_old = s% dt
84 11 : s% mstar_old = s% mstar
85 11 : s% xmstar_old = s% xmstar
86 11 : s% M_center_old = s% M_center
87 11 : s% R_center_old = s% R_center
88 11 : s% L_center_old = s% L_center
89 11 : s% v_center_old = s% v_center
90 11 : s% cumulative_energy_error_old = s% cumulative_energy_error
91 11 : s% total_energy_old = s% total_energy
92 11 : s% total_internal_energy_old = s% total_internal_energy
93 11 : s% total_angular_momentum_old = s% total_angular_momentum
94 11 : s% Teff_old = s% Teff
95 11 : s% power_nuc_burn_old = s% power_nuc_burn
96 11 : s% power_h_burn_old = s% power_h_burn
97 11 : s% power_he_burn_old = s% power_he_burn
98 11 : s% power_z_burn_old = s% power_z_burn
99 11 : s% power_photo_old = s% power_photo
100 11 : s% mstar_dot_old = s% mstar_dot
101 11 : s% L_phot_old = s% L_phot
102 11 : s% L_surf_old = s% L_surf
103 11 : s% dt_limit_ratio_old = s% dt_limit_ratio
104 11 : s% gradT_excess_alpha_old = s% gradT_excess_alpha
105 11 : s% crystal_core_boundary_mass_old = s% crystal_core_boundary_mass
106 :
107 11 : do j = 1, s% len_extra_work
108 11 : s% extra_work_old(j) = s% extra_work(j)
109 : end do
110 :
111 11 : do j = 1, s% len_extra_iwork
112 11 : s% extra_iwork_old(j) = s% extra_iwork(j)
113 : end do
114 :
115 341 : s% ixtra_old = s% ixtra
116 341 : s% xtra_old = s% xtra
117 341 : s% lxtra_old = s% lxtra
118 :
119 11 : call s% other_new_generation(s% id, ierr)
120 :
121 11 : s% need_to_setvars = .true.
122 :
123 : contains
124 :
125 55 : subroutine copy_to_old(ptr, ptr_old, ierr)
126 : real(dp), pointer, dimension(:) :: ptr, ptr_old
127 : integer, intent(out) :: ierr
128 : logical :: first_time
129 : ierr = 0
130 55 : first_time = (.not. associated(ptr_old))
131 55 : call realloc_if_needed_1(ptr_old,nz,nz_alloc_extra,ierr)
132 55 : if (ierr /= 0) return
133 55 : if (s% fill_arrays_with_NaNs) then
134 0 : call fill_with_NaNs(ptr_old)
135 55 : else if (s% zero_when_allocate) then
136 0 : ptr_old(:) = 0
137 55 : else if (first_time) then
138 12325 : ptr_old(1:nz) = -9d99
139 : end if
140 65490 : do j = 1, s% nz
141 65490 : ptr_old(j) = ptr(j)
142 : end do
143 : end subroutine copy_to_old
144 :
145 : end subroutine new_generation
146 :
147 11 : subroutine set_current_to_old(s)
148 : use star_utils, only: total_angular_momentum, set_m_and_dm, set_dm_bar, set_qs
149 : use hydro_rotation, only: use_xh_to_update_i_rot
150 : use utils_lib
151 : type (star_info), pointer :: s
152 : integer :: j, k, ierr
153 :
154 : include 'formats'
155 :
156 11 : s% nz = s% nz_old
157 11 : s% mstar = s% mstar_old
158 11 : s% xmstar = s% xmstar_old
159 11 : s% M_center = s% M_center_old
160 11 : s% R_center = s% R_center_old
161 11 : s% L_center = s% L_center_old
162 11 : s% v_center = s% v_center_old
163 11 : s% cumulative_energy_error = s% cumulative_energy_error_old
164 11 : s% total_energy = s% total_energy_old
165 11 : s% total_internal_energy = s% total_internal_energy_old
166 11 : s% total_angular_momentum = s% total_angular_momentum_old
167 11 : s% Teff = s% Teff_old
168 11 : s% power_nuc_burn = s% power_nuc_burn_old
169 11 : s% power_h_burn = s% power_h_burn_old
170 11 : s% power_he_burn = s% power_he_burn_old
171 11 : s% power_z_burn = s% power_z_burn_old
172 11 : s% power_photo = s% power_photo_old
173 11 : s% mstar_dot = s% mstar_dot_old
174 11 : s% L_phot = s% L_phot_old
175 11 : s% L_surf = s% L_surf_old
176 11 : s% dt_limit_ratio = s% dt_limit_ratio_old
177 11 : s% gradT_excess_alpha = s% gradT_excess_alpha_old
178 11 : s% crystal_core_boundary_mass = s% crystal_core_boundary_mass_old
179 :
180 11 : if (.not. s% RSP_flag) then
181 11 : if (s% fill_arrays_with_NaNs) then
182 0 : call fill_with_NaNs_2d(s% xh)
183 0 : call fill_with_NaNs_2d(s% xa)
184 : end if
185 13098 : do k = 1, s% nz
186 65435 : do j=1, s% nvar_hydro
187 65435 : s% xh(j,k) = s% xh_old(j,k)
188 : end do
189 117783 : do j=1,s% species
190 117783 : s% xa(j,k) = s% xa_old(j,k)
191 : end do
192 13087 : s% dq(k) = s% dq_old(k)
193 13098 : s% mlt_vc(k) = s% mlt_vc_old(k)
194 : end do
195 11 : s% okay_to_set_mlt_vc = .true.
196 :
197 11 : call set_qs(s, s% nz, s% q, s% dq, ierr)
198 11 : if (ierr /= 0) then
199 0 : write(*,*) 'set_current_to_old failed in set_qs'
200 0 : stop
201 : end if
202 11 : call set_m_and_dm(s)
203 11 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
204 :
205 11 : if (s% rotation_flag) then
206 0 : do k=1,s% nz
207 0 : s% j_rot(k) = s% j_rot_old(k)
208 : end do
209 0 : call use_xh_to_update_i_rot(s)
210 0 : do k=1,s% nz
211 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
212 0 : if (is_bad_num(s% omega(k)) .or. abs(s% omega(k)) > 1d50) then
213 0 : if (s% stop_for_bad_nums) then
214 0 : write(*,2) 's% omega(k)', k, s% omega(k)
215 0 : stop 'set_current_to_old'
216 : end if
217 : end if
218 : end do
219 0 : s% total_angular_momentum = total_angular_momentum(s)
220 : end if
221 :
222 : end if
223 :
224 11 : do j = 1, s% len_extra_work
225 11 : s% extra_work(j) = s% extra_work_old(j)
226 : end do
227 :
228 11 : do j = 1, s% len_extra_iwork
229 11 : s% extra_iwork(j) = s% extra_iwork_old(j)
230 : end do
231 :
232 341 : s% ixtra = s% ixtra_old
233 341 : s% xtra = s% xtra_old
234 341 : s% lxtra = s% lxtra_old
235 :
236 11 : call s% other_set_current_to_old(s% id)
237 :
238 11 : s% need_to_setvars = .true.
239 :
240 11 : end subroutine set_current_to_old
241 :
242 :
243 1 : subroutine output(id, ierr)
244 11 : use star_utils, only: get_name_for_restart_file
245 : interface
246 : subroutine save_restart_info(iounit, id, ierr)
247 : implicit none
248 : integer, intent(in) :: iounit
249 : integer, intent(in) :: id
250 : integer, intent(out) :: ierr
251 : end subroutine save_restart_info
252 : end interface
253 : integer, intent(in) :: id
254 : integer, intent(out) :: ierr
255 :
256 : character (len=strlen) :: filename, num_str, fstring
257 : type (star_info), pointer :: s
258 : integer :: num_digits
259 :
260 1 : call get_star_ptr(id, s, ierr)
261 1 : if (ierr /= 0) return
262 1 : call get_name_for_restart_file(s% model_number, s% photo_digits, num_str)
263 1 : filename = trim(s% photo_directory) // '/' // trim(num_str)
264 1 : call output_to_file(filename, id, ierr)
265 1 : if (ierr /= 0) return
266 :
267 1 : write(*, '(a)', advance='no') 'save ' // trim(filename)
268 1 : num_digits = 1 + log10(dble(max(1,s% model_number)))
269 1 : write(fstring,'( "(a,i",i2.2,".",i2.2,")" )') num_digits, num_digits
270 1 : write(*,fstring) ' for model ', s% model_number
271 :
272 1 : end subroutine output
273 :
274 :
275 1 : subroutine output_to_file(filename, id, ierr)
276 1 : use photo_out, only: output_star_photo
277 : use utils_lib, only: folder_exists, mkdir
278 : character (len=*) :: filename
279 : integer, intent(in) :: id
280 : integer, intent(out) :: ierr
281 :
282 : integer :: iounit
283 : type (star_info), pointer :: s
284 : character(len=strlen) :: iomsg
285 :
286 : include 'formats'
287 :
288 1 : call get_star_ptr(id, s, ierr)
289 1 : if (ierr /= 0) return
290 :
291 1 : if(.not. folder_exists(trim(s% photo_directory))) call mkdir(trim(s% photo_directory))
292 :
293 : open(newunit=iounit, file=trim(filename), action='write', &
294 1 : status='replace', iostat=ierr, iomsg=iomsg, form='unformatted')
295 1 : if (ierr == 0) then
296 1 : s% most_recent_photo_name = trim(filename)
297 1 : call output_star_photo(s, iounit, ierr)
298 1 : close(iounit)
299 : else
300 0 : write(*,*) trim(iomsg)
301 : end if
302 :
303 2 : end subroutine output_to_file
304 :
305 : end module evolve_support
|
