Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-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 write_model
21 :
22 : use star_private_def
23 : use const_def, only: dp, lsun, rsun, four_thirds_pi
24 : use read_model
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_write_model
30 :
31 : contains
32 :
33 0 : subroutine do_write_model(id, filename, ierr)
34 : use utils_lib
35 : use chem_def
36 : use star_utils, only: get_scale_height_face_val
37 : integer, intent(in) :: id
38 : character (len=*), intent(in) :: filename
39 : integer, intent(out) :: ierr
40 :
41 : integer :: iounit, i, k, nvar_hydro, nz, species, file_type
42 : integer, pointer :: chem_id(:)
43 : type (star_info), pointer :: s
44 : logical :: v_flag, RTI_flag, &
45 : RSP2_flag, u_flag, prev_flag, rotation_flag, &
46 : write_mlt_vc, RSP_flag
47 :
48 : 1 format(a32, 2x, 1pd26.16)
49 : 11 format(a32, 2x, 1pd26.16, 2x, a, 2x, 99(1pd26.16))
50 : 2 format(a32, 2x, i30)
51 : 3 format(a32, 3x, a8)
52 : 4 format(a32, 3x, a)
53 :
54 : ierr = 0
55 0 : call get_star_ptr(id, s, ierr)
56 0 : if (ierr /= 0) return
57 :
58 0 : chem_id => s% chem_id
59 0 : nvar_hydro = s% nvar_hydro
60 0 : nz = s% nz
61 0 : v_flag = s% v_flag
62 0 : u_flag = s% u_flag
63 0 : RTI_flag = s% RTI_flag
64 0 : rotation_flag = s% rotation_flag
65 0 : RSP_flag = s% RSP_flag
66 0 : RSP2_flag = s% RSP2_flag
67 0 : write_mlt_vc = s% have_mlt_vc
68 :
69 0 : species = s% species
70 :
71 0 : open(newunit=iounit, file=trim(filename), action='write', status='replace')
72 0 : write(iounit,'(a)') '! note: initial lines of file can contain comments'
73 0 : write(iounit,'(a)') '!'
74 0 : prev_flag = (s% nz_old == s% nz .and. s% generations > 1)
75 0 : file_type = 0
76 0 : if (RTI_flag) file_type = file_type + 2**bit_for_RTI
77 0 : if (prev_flag) file_type = file_type + 2**bit_for_2models
78 0 : if (v_flag) file_type = file_type + 2**bit_for_velocity
79 0 : if (u_flag) file_type = file_type + 2**bit_for_u
80 0 : if (rotation_flag) file_type = file_type + 2**bit_for_rotation
81 0 : if (rotation_flag) file_type = file_type + 2**bit_for_j_rot
82 0 : if (RSP_flag) file_type = file_type + 2**bit_for_RSP
83 0 : if (RSP2_flag) file_type = file_type + 2**bit_for_RSP2
84 0 : if (write_mlt_vc) file_type = file_type + 2**bit_for_mlt_vc
85 :
86 0 : write(iounit, '(i14)', advance='no') file_type
87 0 : write(iounit,'(a)',advance='no') ' -- model for mesa/star'
88 0 : if (BTEST(file_type, bit_for_velocity)) &
89 0 : write(iounit,'(a)',advance='no') ', cell boundary velocities (v)'
90 0 : if (BTEST(file_type, bit_for_rotation)) &
91 0 : write(iounit,'(a)',advance='no') ', angular velocities (omega)'
92 0 : if (BTEST(file_type, bit_for_j_rot)) &
93 0 : write(iounit,'(a)',advance='no') ', specific angular momentum (j_rot)'
94 0 : if (BTEST(file_type, bit_for_D_omega)) &
95 0 : write(iounit,'(a)',advance='no') ', omega diffusion coefficients (D_omega)'
96 0 : if (BTEST(file_type, bit_for_am_nu_rot)) &
97 0 : write(iounit,'(a)',advance='no') ', am_nu_rot diffusion coefficients'
98 0 : if (BTEST(file_type, bit_for_u)) &
99 0 : write(iounit,'(a)',advance='no') ', cell center Riemann velocities (u)'
100 0 : if (BTEST(file_type, bit_for_RTI)) &
101 0 : write(iounit,'(a)',advance='no') ', Rayleigh-Taylor instabilities (alpha_RTI)'
102 0 : if (BTEST(file_type, bit_for_mlt_vc)) &
103 0 : write(iounit,'(a)',advance='no') ', mlt convection velocity (mlt_vc)'
104 0 : if (BTEST(file_type, bit_for_RSP)) &
105 0 : write(iounit,'(a)',advance='no') ', RSP luminosity (L), turbulent energy (et_RSP), and radiative flux (erad_RSP)'
106 0 : if (BTEST(file_type, bit_for_RSP2)) &
107 0 : write(iounit,'(a)',advance='no') ', RSP2 turbulent energy (w^2) and pressure scale height (Hp)'
108 : write(iounit,'(a)',advance='no') &
109 0 : '. cgs units. lnd=ln(density), lnT=ln(temperature), lnR=ln(radius)'
110 0 : if (s% i_lum /= 0) then
111 0 : write(iounit,'(a)',advance='no') ', L=luminosity'
112 : end if
113 0 : if (s% M_center /= 0) then
114 : write(iounit,'(a)',advance='no') &
115 0 : ', dq=fraction of xmstar=(mstar-mcenter) in cell; remaining cols are mass fractions.'
116 : else
117 : write(iounit,'(a)',advance='no') &
118 0 : ', dq=fraction of total mstar in cell; remaining cols are mass fractions.'
119 : end if
120 0 : write(iounit,'(A)')
121 : ! write property list
122 0 : write(iounit, '(a)') ! blank line before start of property list
123 0 : write(iounit, 4) 'version_number', "'" // trim(version_number) // "'"
124 0 : write(iounit, 1) 'M/Msun', s% star_mass
125 0 : write(iounit, 2) 'model_number', s% model_number
126 0 : write(iounit, 1) 'star_age', s% star_age
127 0 : write(iounit, 1) 'initial_z', s% initial_z
128 0 : write(iounit, 2) 'n_shells', nz
129 0 : write(iounit, 4) 'net_name', "'" // trim(s% net_name) // "'"
130 0 : write(iounit, 2) 'species', species
131 0 : if (s% M_center /= 0) then
132 0 : write(iounit, 11) 'xmstar', s% xmstar, &
133 0 : '! above core (g). core mass: Msun, grams:', &
134 0 : s% M_center/Msun, s% M_center
135 : end if
136 0 : if (s% R_center /= 0) then
137 0 : write(iounit, 11) 'R_center', s% R_center, &
138 0 : '! radius of core (cm). R/Rsun, avg core density (g/cm^3):', &
139 0 : s% R_center/Rsun, s% M_center/(four_thirds_pi*pow3(s% R_center))
140 : end if
141 0 : if (s% v_center /= 0) then
142 0 : write(iounit, 11) 'v_center', s% v_center, &
143 0 : '! velocity of outer edge of core (cm/s)'
144 : end if
145 0 : if (s% L_center /= 0) then
146 0 : write(iounit, 11) 'L_center', s% L_center, &
147 0 : '! luminosity of core (erg/s). L/Lsun, avg core eps (erg/g/s):', &
148 0 : s% L_center/Lsun, s% L_center/max(1d0,s% M_center)
149 : end if
150 0 : if (s% tau_factor /= 1) then
151 0 : write(iounit, 1) 'tau_factor', s% tau_factor
152 : end if
153 0 : if (s% Tsurf_factor /= 1) then
154 0 : write(iounit, 1) 'Tsurf_factor', s% Tsurf_factor
155 : end if
156 0 : if (s% opacity_factor /= 1) then
157 0 : write(iounit, 1) 'opacity_factor', s% opacity_factor
158 : end if
159 0 : if (s% crystal_core_boundary_mass > 0d0) then
160 0 : write(iounit, 1) 'crystal_core_boundary_mass', s% crystal_core_boundary_mass
161 : end if
162 0 : write(iounit, 1) 'Teff', s% Teff
163 0 : write(iounit, 1) 'power_nuc_burn', s% power_nuc_burn
164 0 : write(iounit, 1) 'power_h_burn', s% power_h_burn
165 0 : write(iounit, 1) 'power_he_burn', s% power_he_burn
166 0 : write(iounit, 1) 'power_z_burn', s% power_z_burn
167 0 : write(iounit, 1) 'power_photo', s% power_photo
168 0 : write(iounit, 1) 'total_energy', s% total_energy
169 0 : write(iounit, 1) 'cumulative_energy_error', s% cumulative_energy_error
170 0 : if (abs(s% total_energy) > 1d0) &
171 0 : write(iounit, 11) 'cumulative_error/total_energy', &
172 0 : s% cumulative_energy_error/s% total_energy, &
173 0 : 'log_rel_run_E_err', &
174 0 : safe_log10(abs(s% cumulative_energy_error/s% total_energy))
175 0 : write(iounit, 2) 'num_retries', s% num_retries
176 0 : write(iounit, '(a)') ! blank line for end of property list
177 :
178 0 : call header
179 0 : do k=1, nz
180 0 : if (ierr /= 0) exit
181 0 : write(iounit, fmt='(i5, 1x)', advance='no') k
182 0 : call write1(s% lnd(k),ierr); if (ierr /= 0) exit
183 0 : call write1(s% lnT(k),ierr); if (ierr /= 0) exit
184 0 : call write1(s% lnR(k),ierr); if (ierr /= 0) exit
185 0 : if (RSP_flag) then
186 0 : call write1(s% RSP_Et(k),ierr); if (ierr /= 0) exit
187 0 : call write1(s% erad(k),ierr); if (ierr /= 0) exit
188 0 : call write1(s% Fr(k),ierr); if (ierr /= 0) exit
189 0 : else if (RSP2_flag) then
190 0 : call write1(s% w(k),ierr); if (ierr /= 0) exit
191 0 : call write1(s% Hp_face(k),ierr)
192 0 : if (ierr /= 0) exit
193 : end if
194 0 : call write1(s% L(k),ierr); if (ierr /= 0) exit
195 0 : call write1(s% dq(k),ierr); if (ierr /= 0) exit
196 0 : if (v_flag) then
197 0 : call write1(s% v(k),ierr); if (ierr /= 0) exit
198 : end if
199 0 : if (rotation_flag) then
200 0 : call write1(s% omega(k),ierr); if (ierr /= 0) exit
201 0 : call write1(s% j_rot(k),ierr); if (ierr /= 0) exit
202 : end if
203 0 : if (u_flag) then
204 0 : call write1(s% u(k),ierr); if (ierr /= 0) exit
205 : end if
206 0 : if (RTI_flag) then
207 0 : call write1(s% alpha_RTI(k),ierr); if (ierr /= 0) exit
208 : end if
209 0 : if (write_mlt_vc) then
210 0 : call write1(s% mlt_vc(k),ierr); if (ierr /= 0) exit
211 : end if
212 0 : do i=1, species
213 0 : call write1(s% xa(i,k),ierr); if (ierr /= 0) exit
214 : end do
215 0 : write(iounit, *)
216 : end do
217 :
218 0 : if (prev_flag) then
219 0 : write(iounit, '(a)')
220 0 : write(iounit, '(a)') ' previous model'
221 : ! write prev properties
222 0 : write(iounit, '(a)')
223 0 : write(iounit, 2) 'previous n_shells', s% nz_old
224 0 : write(iounit, 1, advance='no') 'previous mass (grams)'
225 0 : call write1_eol(s% mstar_old,ierr)
226 0 : write(iounit, 1, advance='no') 'timestep (seconds)'
227 0 : call write1_eol(s% dt,ierr)
228 0 : write(iounit, 1, advance='no') 'dt_next (seconds)'
229 0 : call write1_eol(s% dt_next_unclipped,ierr)
230 0 : write(iounit, '(a)')
231 : end if
232 0 : close(iounit)
233 :
234 : contains
235 :
236 0 : subroutine write1_eol(val,ierr)
237 : real(dp), intent(in) :: val
238 : integer, intent(out) :: ierr
239 0 : call write1(val,ierr)
240 0 : write(iounit,'(A)') ''
241 0 : end subroutine write1_eol
242 :
243 0 : subroutine write1(val,ierr)
244 : real(dp), intent(in) :: val
245 : integer, intent(out) :: ierr
246 : integer, parameter :: str_len = 26
247 : character (len=str_len) :: string
248 0 : ierr = 0
249 0 : call double_to_str(val,string)
250 0 : write(iounit, fmt='(a26,1x)', advance='no') string
251 0 : end subroutine write1
252 :
253 0 : subroutine header
254 0 : write(iounit, fmt='(10x, a9, 1x, 99(a26, 1x))', advance='no') 'lnd', 'lnT', 'lnR'
255 0 : if (RSP_flag) then
256 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'et_RSP'
257 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'erad_RSP'
258 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'Fr_RSP'
259 0 : else if (RSP2_flag) then
260 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'w'
261 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'Hp'
262 : end if
263 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'L'
264 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'dq'
265 0 : if (v_flag) write(iounit, fmt='(a26, 1x)', advance='no') 'v'
266 0 : if (rotation_flag) write(iounit, fmt='(a26, 1x)', advance='no') 'omega'
267 0 : if (rotation_flag) write(iounit, fmt='(a26, 1x)', advance='no') 'j_rot'
268 0 : if (u_flag) write(iounit, fmt='(a26, 1x)', advance='no') 'u'
269 0 : if (RTI_flag) &
270 0 : write(iounit, fmt='(a26, 1x)', advance='no') 'alpha_RTI'
271 0 : if (write_mlt_vc) write(iounit, fmt='(a26, 1x)', advance='no') 'mlt_vc'
272 0 : do i=1, species
273 0 : write(iounit, fmt='(a26, 1x)', advance='no') chem_isos% name(chem_id(i))
274 : end do
275 0 : write(iounit, *)
276 0 : end subroutine header
277 :
278 : end subroutine do_write_model
279 :
280 :
281 : end module write_model
|
