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 mesh_functions
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, msun, rsun, secyer
24 : use num_lib
25 : use utils_lib
26 : use chem_def
27 :
28 : implicit none
29 :
30 : private
31 : public :: num_mesh_functions, set_mesh_function_data, max_allowed_gvals
32 :
33 : integer, parameter :: max_allowed_gvals = 50
34 :
35 : contains
36 :
37 30 : integer function get_net_iso(s, species)
38 : use chem_lib, only: chem_get_iso_id
39 : type (star_info), pointer :: s
40 : character (len=*) :: species
41 : integer :: j
42 30 : get_net_iso = 0
43 30 : if (len_trim(species) == 0) return
44 30 : j = chem_get_iso_id(species)
45 30 : if (j <= 0) then
46 0 : write(*,*) 'unknown species name for mesh function: ' // trim(species)
47 0 : write(*,*) 'len_trim(species)', len_trim(species)
48 0 : return
49 : end if
50 30 : get_net_iso = s% net_iso(j) ! 0 if species not in current net
51 30 : end function get_net_iso
52 :
53 :
54 180 : logical function do_mass_function(s,species,weight,j)
55 : type (star_info), pointer :: s
56 : character (len=iso_name_length) :: species
57 : real(dp), intent(in) :: weight
58 : integer, intent(out) :: j
59 180 : j = 0
60 20 : if (weight > 0) j = get_net_iso(s, species)
61 200 : do_mass_function = (j > 0)
62 30 : end function do_mass_function
63 :
64 :
65 10 : integer function num_mesh_functions(s)
66 : type (star_info), pointer :: s
67 : integer :: i, j, k
68 10 : i = 0
69 10 : if (s% use_other_mesh_functions) &
70 0 : call s% how_many_other_mesh_fcns(s% id, i)
71 10 : if (s% convective_bdy_weight > 0) i=i+1
72 10 : if (s% E_function_weight > 0) i=i+1
73 10 : if (s% P_function_weight > 0) i=i+1
74 10 : if (s% T_function1_weight > 0) i=i+1
75 10 : if (s% T_function2_weight > 0) i=i+1
76 10 : if (s% R_function_weight > 0) i=i+1
77 10 : if (s% R_function2_weight > 0) i=i+1
78 10 : if (s% R_function3_weight > 0) i=i+1
79 10 : if (s% M_function_weight > 0) i=i+1
80 10 : if (s% gradT_function_weight > 0) i=i+1
81 10 : if (s% log_tau_function_weight > 0) i=i+1
82 10 : if (s% log_kap_function_weight > 0) i=i+1
83 10 : if (s% gam_function_weight > 0) i=i+1
84 10 : if (s% omega_function_weight > 0 .and. s% rotation_flag) i=i+1
85 100 : do k=1,num_xa_function
86 110 : if (do_mass_function( &
87 : s, s% xa_function_species(k), s% xa_function_weight(k), j)) &
88 20 : i=i+1
89 : end do
90 10 : num_mesh_functions = i
91 10 : end function num_mesh_functions
92 :
93 :
94 10 : subroutine set_mesh_function_data( &
95 : s, nfcns, names, gval_is_xa_function, gval_is_logT_function, vals1, ierr)
96 : type (star_info), pointer :: s
97 : integer, intent(in) :: nfcns
98 : character (len=32), intent(out) :: names(max_allowed_gvals)
99 : logical, intent(out), dimension(max_allowed_gvals) :: &
100 : gval_is_xa_function, gval_is_logT_function
101 : real(dp), pointer :: vals1(:) ! =(nz, nfcns)
102 : integer, intent(out) :: ierr
103 :
104 : integer :: i, nz, j, k, i_other
105 : logical, parameter :: dbg = .false.
106 : real(dp), dimension(:,:), pointer :: vals
107 :
108 10 : ierr = 0
109 10 : nz = s% nz
110 :
111 10 : vals(1:nz,1:nfcns) => vals1(1:nz*nfcns)
112 10 : gval_is_xa_function = .false.
113 10 : gval_is_logT_function = .false.
114 :
115 10 : i_other = 0
116 10 : if (s% use_other_mesh_functions) then
117 0 : call s% how_many_other_mesh_fcns(s% id, i_other)
118 0 : if (i_other > 0) then
119 0 : if (i_other > nfcns) then
120 0 : ierr = -1
121 0 : return
122 : end if
123 : call s% other_mesh_fcn_data( &
124 0 : s% id, i_other, names, gval_is_xa_function, vals1, ierr)
125 0 : if (ierr /= 0) return
126 : end if
127 : end if
128 :
129 10 : i = i_other
130 10 : if (s% E_function_weight > 0) then
131 0 : i = i+1; names(i) = 'E_function'
132 : end if
133 10 : if (s% P_function_weight > 0) then
134 10 : i = i+1; names(i) = 'P_function'
135 : end if
136 10 : if (s% T_function1_weight > 0) then
137 10 : i = i+1; names(i) = 'T_function1'
138 10 : gval_is_logT_function(i) = .true.
139 : end if
140 10 : if (s% T_function2_weight > 0) then
141 0 : i = i+1; names(i) = 'T_function2'
142 0 : gval_is_logT_function(i) = .true.
143 : end if
144 10 : if (s% R_function_weight > 0) then
145 0 : i = i+1; names(i) = 'R_function'
146 : end if
147 10 : if (s% R_function2_weight > 0) then
148 0 : i = i+1; names(i) = 'R_function2'
149 : end if
150 10 : if (s% R_function3_weight > 0) then
151 0 : i = i+1; names(i) = 'R_function3'
152 : end if
153 10 : if (s% M_function_weight > 0) then
154 0 : i = i+1; names(i) = 'M_function'
155 : end if
156 10 : if (s% gradT_function_weight > 0) then
157 0 : i = i+1; names(i) = 'gradT_function'
158 : end if
159 10 : if (s% log_tau_function_weight > 0) then
160 0 : i = i+1; names(i) = 'log_tau_function'
161 : end if
162 10 : if (s% log_kap_function_weight > 0) then
163 0 : i = i+1; names(i) = 'log_kap_function'
164 : end if
165 10 : if (s% gam_function_weight > 0) then
166 0 : i = i+1; names(i) = 'gam_function'
167 : end if
168 10 : if (s% omega_function_weight > 0 .and. s% rotation_flag) then
169 0 : i = i+1; names(i) = 'omega_function'
170 : end if
171 10 : if (s% convective_bdy_weight > 0) then
172 0 : i = i+1; names(i) = 'convective_bdy'
173 : end if
174 100 : do k=1,num_xa_function
175 120 : if (do_mass_function(s, s% xa_function_species(k), s% xa_function_weight(k), j)) then
176 10 : i = i+1; names(i) = trim(s% xa_function_species(k))
177 10 : gval_is_xa_function(i) = .true.
178 : !write(names(i),'(a,i1)') 'xa_function_', k
179 : end if
180 : end do
181 10 : if (i /= nfcns) then
182 0 : write(*,*) 'error in set_mesh_function_names: incorrect nfcns'
183 0 : ierr = -1
184 : end if
185 :
186 40 : do i=i_other+1,nfcns
187 :
188 30 : if (ierr /= 0) cycle
189 : if (dbg) write(*,*) trim(names(i))
190 :
191 40 : if (names(i) == 'E_function') then
192 0 : do k=1,nz
193 0 : vals(k,i) = s% E_function_weight * max(s% E_function_param, s% lnE(k)/ln10)
194 : end do
195 :
196 30 : else if (names(i) == 'P_function') then
197 12105 : do k=1,nz
198 12105 : vals(k,i) = s% P_function_weight * s% lnPeos(k)/ln10
199 : end do
200 :
201 20 : else if (names(i) == 'T_function1') then
202 12105 : do k=1,nz
203 12105 : vals(k,i) = s% T_function1_weight*s% lnT(k)/ln10
204 : end do
205 :
206 10 : else if (names(i) == 'T_function2') then
207 0 : do k=1,nz
208 : vals(k,i) = &
209 0 : s% T_function2_weight*log10(s% T(k) / (s% T(k) + s% T_function2_param))
210 : end do
211 :
212 10 : else if (names(i) == 'R_function') then
213 0 : do k=1,nz
214 : vals(k,i) = &
215 0 : s% R_function_weight*log10(1 + (s% r(k)/Rsun)/s% R_function_param)
216 : end do
217 :
218 10 : else if (names(i) == 'R_function2') then
219 0 : do k=1,nz
220 : vals(k,i) = &
221 : s% R_function2_weight * &
222 0 : min(s% R_function2_param1, max(s% R_function2_param2,s% r(k)/s% r(1)))
223 : end do
224 :
225 10 : else if (names(i) == 'R_function3') then
226 0 : do k=1,nz
227 0 : vals(k,i) = s% R_function3_weight*s% r(k)/s% r(1)
228 : end do
229 :
230 10 : else if (names(i) == 'M_function') then
231 0 : do k=1,nz
232 : vals(k,i) = &
233 0 : s% M_function_weight*log10(1 + (s% xmstar*s% q(k)/Msun)/s% M_function_param)
234 : end do
235 :
236 10 : else if (names(i) == 'gradT_function') then
237 0 : do k=1,nz
238 0 : vals(k,i) = s% gradT_function_weight*s% gradT(k)
239 : end do
240 :
241 10 : else if (names(i) == 'log_tau_function') then
242 0 : do k=1,nz
243 0 : vals(k,i) = s% log_tau_function_weight*log10(s% tau(k))
244 : end do
245 :
246 10 : else if (names(i) == 'log_kap_function') then
247 0 : do k=1,nz
248 0 : vals(k,i) = s% log_kap_function_weight*log10(s% opacity(k))
249 : end do
250 :
251 10 : else if (names(i) == 'gam_function') then
252 0 : do k=1,nz
253 : vals(k,i) = s% gam_function_weight* &
254 0 : tanh((s% gam(k) - s% gam_function_param1)/s% gam_function_param2)
255 : end do
256 :
257 10 : else if (names(i) == 'omega_function') then
258 0 : do k=1,nz
259 0 : vals(k,i) = s% omega_function_weight*log10(max(1d-99,abs(s% omega(k))))
260 : end do
261 :
262 10 : else if (names(i) == 'convective_bdy') then
263 0 : call do_conv_bdy(i)
264 :
265 : else
266 100 : do k=1,num_xa_function
267 100 : call do1_xa_function(k,i)
268 : end do
269 :
270 : end if
271 :
272 : end do
273 :
274 : contains
275 :
276 :
277 0 : subroutine do_conv_bdy(i)
278 : integer, intent(in) :: i
279 : integer :: k, j
280 :
281 0 : vals(1:nz,i) = 0
282 :
283 0 : if (s% dt < s% convective_bdy_min_dt_yrs*secyer) return
284 :
285 : !do k=1,nz
286 : ! if (s% dq(k) < s% convective_bdy_dq_limit) cycle
287 : ! if (s% cz_bdy_dq(k) /= 0) vals(k,i) = s% convective_bdy_weight
288 : !end do
289 :
290 0 : do j = 1, s% num_conv_boundaries
291 0 : k = s% conv_bdy_loc(j)
292 0 : if (s% dq(k) < s% convective_bdy_dq_limit) cycle
293 0 : vals(k,i) = s% convective_bdy_weight
294 0 : if (k < nz .and. s% top_conv_bdy(i)) then
295 0 : vals(k+1,i) = s% convective_bdy_weight
296 0 : else if (k > 1 .and. .not. s% top_conv_bdy(i)) then
297 0 : vals(k-1,i) = s% convective_bdy_weight
298 : end if
299 : end do
300 :
301 0 : do k=2,nz
302 0 : vals(k,i) = vals(k,i) + vals(k-1,i)
303 : end do
304 :
305 : end subroutine do_conv_bdy
306 :
307 90 : subroutine do1_xa_function(k,i)
308 : integer, intent(in) :: k,i
309 90 : real(dp) :: weight, param
310 : integer :: j, m
311 90 : if (len_trim(s% xa_function_species(k))==0) return
312 10 : if (trim(names(i)) /= trim(s% xa_function_species(k))) return
313 10 : j = get_net_iso(s, s% xa_function_species(k))
314 10 : if (j <= 0) then
315 0 : ierr = -1
316 : else
317 10 : weight = s% xa_function_weight(k)
318 10 : param = s% xa_function_param(k)
319 12105 : do m=1,s% nz
320 12105 : vals(m,i) = weight*log10(s% xa(j,m) + param)
321 : end do
322 : end if
323 : end subroutine do1_xa_function
324 :
325 : end subroutine set_mesh_function_data
326 :
327 : end module mesh_functions
|
