Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2020 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 kap_ctrls_io
21 :
22 : use const_def, only: dp
23 : use utils_namelist, only: max_extra_inlists
24 : use kap_def
25 : use math_lib
26 :
27 : implicit none
28 :
29 : public :: read_kap_namelist, write_namelist, get_kap_controls, set_kap_controls
30 : private
31 :
32 : real(dp) :: Zbase
33 :
34 : integer :: kap_option, kap_CO_option, kap_lowT_option
35 :
36 : character(len=strlen) :: kap_file_prefix, kap_CO_prefix, kap_lowT_prefix
37 :
38 : ! user table info
39 : integer :: user_num_kap_Xs = 0
40 : real(dp), dimension(kap_max_dim) :: user_kap_Xs = -1d0
41 : integer :: user_num_kap_Zs = 0
42 : real(dp), dimension(kap_max_dim) :: user_kap_Zs= -1d0
43 : integer, dimension(kap_max_dim) :: user_num_kap_Xs_for_this_Z = 0
44 :
45 : integer :: user_num_kap_CO_Xs = 0
46 : real(dp), dimension(kap_max_dim) :: user_kap_CO_Xs = -1d0
47 : integer :: user_num_kap_CO_Zs = 0
48 : real(dp), dimension(kap_max_dim) :: user_kap_CO_Zs= -1d0
49 : integer, dimension(kap_max_dim) :: user_num_kap_CO_Xs_for_this_Z = 0
50 :
51 : integer :: user_num_kap_lowT_Xs = 0
52 : real(dp), dimension(kap_max_dim) :: user_kap_lowT_Xs = -1d0
53 : integer :: user_num_kap_lowT_Zs = 0
54 : real(dp), dimension(kap_max_dim) :: user_kap_lowT_Zs= -1d0
55 : integer, dimension(kap_max_dim) :: user_num_kap_lowT_Xs_for_this_Z = 0
56 :
57 :
58 : real(dp) :: kap_blend_logT_upper_bdy, kap_blend_logT_lower_bdy
59 :
60 : logical :: cubic_interpolation_in_X
61 : logical :: cubic_interpolation_in_Z
62 : logical :: include_electron_conduction
63 : logical :: use_blouin_conductive_opacities
64 :
65 : logical :: use_Zbase_for_Type1
66 : logical :: use_Type2_opacities
67 :
68 : real(dp) :: kap_Type2_full_off_X, kap_Type2_full_on_X
69 : real(dp) :: kap_Type2_full_off_dZ, kap_Type2_full_on_dZ
70 :
71 : logical :: show_info
72 :
73 : ! hooks
74 : logical :: use_other_elect_cond_opacity, &
75 : use_other_compton_opacity, use_other_radiative_opacity
76 :
77 : ! debugging
78 : logical :: dbg
79 : real(dp) :: logT_lo, logT_hi
80 : real(dp) :: logRho_lo, logRho_hi
81 : real(dp) :: X_lo, X_hi
82 : real(dp) :: Z_lo, Z_hi
83 :
84 : logical, dimension(max_extra_inlists) :: read_extra_kap_inlist
85 : character (len=strlen), dimension(max_extra_inlists) :: extra_kap_inlist_name
86 :
87 : ! User supplied inputs
88 : real(dp) :: kap_ctrl(10)
89 : integer :: kap_integer_ctrl(10)
90 : logical :: kap_logical_ctrl(10)
91 : character(len=strlen) :: kap_character_ctrl(10)
92 :
93 : namelist /kap/ &
94 :
95 : Zbase, &
96 :
97 : kap_file_prefix, kap_CO_prefix, kap_lowT_prefix, aesopus_filename, &
98 :
99 : user_num_kap_Xs, user_kap_Xs, &
100 : user_num_kap_Zs, user_kap_Zs, user_num_kap_Xs_for_this_Z, &
101 :
102 : user_num_kap_CO_Xs, user_kap_CO_Xs, &
103 : user_num_kap_CO_Zs, user_kap_CO_Zs, user_num_kap_CO_Xs_for_this_Z, &
104 :
105 : user_num_kap_lowT_Xs, user_kap_lowT_Xs, &
106 : user_num_kap_lowT_Zs, user_kap_lowT_Zs, user_num_kap_lowT_Xs_for_this_Z, &
107 :
108 : kap_blend_logT_upper_bdy, kap_blend_logT_lower_bdy, &
109 :
110 : cubic_interpolation_in_X, cubic_interpolation_in_Z, &
111 :
112 : include_electron_conduction, use_blouin_conductive_opacities, &
113 :
114 : use_Zbase_for_Type1, use_Type2_opacities, &
115 :
116 : kap_Type2_full_off_X, kap_Type2_full_on_X, &
117 : kap_Type2_full_off_dZ, kap_Type2_full_on_dZ, &
118 :
119 : show_info, &
120 :
121 : use_other_elect_cond_opacity, &
122 : use_other_compton_opacity, &
123 : use_other_radiative_opacity, &
124 :
125 : ! User supplied inputs
126 : kap_ctrl, &
127 : kap_integer_ctrl, &
128 : kap_logical_ctrl, &
129 : kap_character_ctrl,&
130 :
131 : read_extra_kap_inlist, extra_kap_inlist_name
132 :
133 : contains
134 :
135 :
136 : ! read a "namelist" file and set parameters
137 5 : subroutine read_kap_namelist(handle, inlist, ierr)
138 : use utils_namelist, only: read_namelist, missing_namelist_warning
139 : integer, intent(in) :: handle
140 : character (len=*), intent(in) :: inlist
141 : integer, intent(out) :: ierr ! 0 means AOK.
142 : type (Kap_General_Info), pointer :: rq
143 :
144 5 : call get_kap_ptr(handle,rq,ierr)
145 :
146 5 : if (ierr /= 0) return
147 :
148 5 : call set_default_controls
149 :
150 5 : if (inlist /= '') then
151 4 : call read_namelist(inlist, read_kap_file, "kap", ierr, missing_namelist_warning)
152 : end if
153 :
154 5 : if (ierr /= 0) return
155 :
156 5 : call store_controls(rq, ierr)
157 : end subroutine read_kap_namelist
158 :
159 4 : subroutine read_kap_file(unit, iostat, iomsg, extra_inlists, extra_inlists_mask)
160 : use const_def, only: strlen
161 : use utils_namelist, only: max_extra_inlists
162 :
163 : integer, intent(in) :: unit
164 : integer, intent(out) :: iostat
165 : character(len=strlen), intent(out) :: iomsg
166 : character(len=strlen), dimension(max_extra_inlists), intent(out) :: extra_inlists
167 : logical, dimension(max_extra_inlists), intent(out) :: extra_inlists_mask
168 :
169 : integer :: i
170 :
171 4 : read_extra_kap_inlist(:) = .false.
172 :
173 4 : read(unit, nml=kap, iostat=iostat, iomsg=iomsg)
174 :
175 4 : if (iostat /= 0) then
176 : return
177 : end if
178 :
179 24 : do i=1, max_extra_inlists
180 20 : extra_inlists(i) = extra_kap_inlist_name(i)
181 24 : extra_inlists_mask(i) = read_extra_kap_inlist(i)
182 : end do
183 :
184 : end subroutine read_kap_file
185 :
186 5 : subroutine set_default_controls
187 : include 'kap.defaults'
188 5 : end subroutine set_default_controls
189 :
190 :
191 5 : subroutine store_controls(rq, ierr)
192 : type (Kap_General_Info), pointer :: rq
193 :
194 : integer :: i, ierr
195 :
196 5 : rq% Zbase = Zbase
197 :
198 5 : rq% cubic_interpolation_in_X = cubic_interpolation_in_X
199 5 : rq% cubic_interpolation_in_Z = cubic_interpolation_in_Z
200 5 : rq% include_electron_conduction = include_electron_conduction
201 5 : rq% use_blouin_conductive_opacities = use_blouin_conductive_opacities
202 5 : rq% use_Zbase_for_Type1 = use_Zbase_for_Type1
203 5 : rq% use_Type2_opacities = use_Type2_opacities
204 :
205 : ! check for limits on full_off/on options
206 5 : if (kap_Type2_full_off_X > 0.71d0) then
207 0 : write(*,*) "kap_Type2_full_off_X must be smaller than 0.71"
208 0 : ierr = -1
209 0 : return
210 : end if
211 5 : if (kap_Type2_full_on_X > 0.71d0) then
212 0 : write(*,*) "kap_Type2_full_on_X must be smaller than 0.71"
213 0 : ierr = -1
214 0 : return
215 : end if
216 5 : if (kap_Type2_full_off_X < kap_Type2_full_on_X) then
217 0 : write(*,*) "kap_Type2_full_off_X has to be bigger than kap_Type2_full_on_X"
218 0 : ierr = -1
219 0 : return
220 : end if
221 5 : if (kap_Type2_full_off_dZ > kap_Type2_full_on_dZ) then
222 0 : write(*,*) "kap_Type2_full_off_dZ has to be smaller than kap_Type2_full_on_dZ"
223 0 : ierr = -1
224 0 : return
225 : end if
226 :
227 5 : rq% kap_Type2_full_off_X = kap_Type2_full_off_X
228 5 : rq% kap_Type2_full_on_X = kap_Type2_full_on_X
229 5 : rq% kap_Type2_full_off_dZ = kap_Type2_full_off_dZ
230 5 : rq% kap_Type2_full_on_dZ = kap_Type2_full_on_dZ
231 :
232 :
233 5 : if (kap_blend_logT_upper_bdy > 0) rq% kap_blend_logT_upper_bdy = kap_blend_logT_upper_bdy
234 5 : if (kap_blend_logT_lower_bdy > 0) rq% kap_blend_logT_lower_bdy = kap_blend_logT_lower_bdy
235 :
236 :
237 5 : kap_option = 0
238 11 : do i=1,kap_options_max
239 11 : if (kap_file_prefix == kap_option_str(i)) then
240 5 : kap_option = i
241 5 : exit
242 : end if
243 : end do
244 5 : if (kap_option == 0) then
245 0 : write(*,*) 'WARNING: unknown kap_file_prefix (assuming user table): ' // trim(kap_file_prefix)
246 0 : kap_option = kap_user
247 0 : kap_option_str(kap_user) = trim(kap_file_prefix)
248 :
249 0 : if (user_num_kap_Xs == 0 .or. user_num_kap_Zs == 0) then
250 0 : write(*,*) 'ERROR: must set user_num_kap_Xs, user_num_kap_Zs, and related variables'
251 0 : ierr = -1
252 0 : return
253 : end if
254 :
255 0 : if (user_num_kap_Xs > kap_max_dim .or. user_num_kap_Zs > kap_max_dim) then
256 0 : write(0,*) ' failed in kap_read_config_file: maximum X or Z dimensions exceeded'
257 0 : write(0,*) ' maximum dimension is ', kap_max_dim
258 0 : write(0,*) ' num_kap_Xs = ', num_kap_Xs
259 0 : write(0,*) ' num_kap_Zs = ', num_kap_Zs
260 0 : ierr = -1
261 0 : return
262 : end if
263 :
264 0 : num_kap_Xs(kap_user) = user_num_kap_Xs
265 0 : kap_Xs(:, kap_user) = user_kap_Xs
266 :
267 0 : num_kap_Zs(kap_user) = user_num_kap_Zs
268 0 : kap_Zs(:, kap_user) = user_kap_Zs
269 :
270 0 : num_kap_Xs_for_this_Z(:, kap_user) = user_num_kap_Xs_for_this_Z
271 :
272 : end if
273 5 : rq% kap_option = kap_option
274 :
275 :
276 5 : kap_CO_option = 0
277 10 : do i=1,kap_CO_options_max
278 10 : if (kap_CO_prefix == kap_CO_option_str(i)) then
279 5 : kap_CO_option = i
280 5 : exit
281 : end if
282 : end do
283 5 : if (kap_CO_option == 0) then
284 0 : write(*,*) 'WARNING: unknown kap_CO_prefix (assuming user table): ' // trim(kap_CO_prefix)
285 0 : kap_CO_option = kap_CO_user
286 0 : kap_CO_option_str(kap_CO_user) = trim(kap_CO_prefix)
287 :
288 0 : if (user_num_kap_CO_Xs == 0 .or. user_num_kap_CO_Zs == 0) then
289 0 : write(*,*) 'ERROR: must set user_num_kap_CO_Xs, user_num_kap_CO_Zs, and related variables'
290 0 : ierr = -1
291 0 : return
292 : end if
293 :
294 0 : num_kap_CO_Xs(kap_CO_user) = user_num_kap_CO_Xs
295 0 : kap_CO_Xs(:, kap_CO_user) = user_kap_CO_Xs
296 :
297 0 : num_kap_CO_Zs(kap_CO_user) = user_num_kap_CO_Zs
298 0 : kap_CO_Zs(:, kap_CO_user) = user_kap_CO_Zs
299 :
300 0 : num_kap_CO_Xs_for_this_Z(:, kap_CO_user) = user_num_kap_CO_Xs_for_this_Z
301 :
302 : end if
303 5 : rq% kap_CO_option = kap_CO_option
304 :
305 :
306 5 : kap_lowT_option = 0
307 26 : do i=1,kap_lowT_options_max
308 26 : if (kap_lowT_prefix == kap_lowT_option_str(i)) then
309 5 : kap_lowT_option = i
310 5 : exit
311 : end if
312 : end do
313 5 : if (kap_lowT_option == 0) then
314 0 : write(*,*) 'WARNING: unknown kap_lowT_prefix (assuming user table): ' // trim(kap_lowT_prefix)
315 0 : kap_lowT_option = kap_lowT_user
316 0 : kap_lowT_option_str(kap_lowT_user) = trim(kap_lowT_prefix)
317 :
318 0 : if (user_num_kap_lowT_Xs == 0 .or. user_num_kap_lowT_Zs == 0) then
319 0 : write(*,*) 'ERROR: must set user_num_kap_lowT_Xs, user_num_kap_lowT_Zs, and related variables'
320 0 : ierr = -1
321 0 : return
322 : end if
323 :
324 0 : if (user_num_kap_lowT_Xs > kap_max_dim .or. user_num_kap_lowT_Zs > kap_max_dim) then
325 0 : write(0,*) ' failed in kap_read_config_file: maximum X or Z dimensions exceeded'
326 0 : write(0,*) ' maximum dimension is ', kap_max_dim
327 0 : write(0,*) ' num_kap_lowT_Xs = ', num_kap_lowT_Xs
328 0 : write(0,*) ' num_kap_lowT_Zs = ', num_kap_lowT_Zs
329 0 : ierr = -1
330 0 : return
331 : end if
332 :
333 0 : num_kap_lowT_Xs(kap_lowT_user) = user_num_kap_lowT_Xs
334 0 : kap_lowT_Xs(:, kap_lowT_user) = user_kap_lowT_Xs
335 :
336 0 : num_kap_lowT_Zs(kap_lowT_user) = user_num_kap_lowT_Zs
337 0 : kap_lowT_Zs(:, kap_lowT_user) = user_kap_lowT_Zs
338 :
339 0 : num_kap_lowT_Xs_for_this_Z(:, kap_lowT_user) = user_num_kap_lowT_Xs_for_this_Z
340 :
341 : end if
342 5 : rq% kap_lowT_option = kap_lowT_option
343 :
344 5 : rq% show_info = show_info
345 :
346 5 : rq% use_other_elect_cond_opacity = use_other_elect_cond_opacity
347 5 : rq% use_other_compton_opacity = use_other_compton_opacity
348 5 : rq% use_other_radiative_opacity = use_other_radiative_opacity
349 :
350 : ! user inputs
351 55 : rq% kap_ctrl = kap_ctrl
352 55 : rq% kap_integer_ctrl = kap_integer_ctrl
353 55 : rq% kap_logical_ctrl = kap_logical_ctrl
354 55 : rq% kap_character_ctrl = kap_character_ctrl
355 :
356 : end subroutine store_controls
357 :
358 :
359 0 : subroutine write_namelist(handle, filename, ierr)
360 : integer, intent(in) :: handle
361 : character(*), intent(in) :: filename
362 : integer, intent(out) :: ierr
363 : type (Kap_General_Info), pointer :: rq
364 : integer :: iounit
365 : open(newunit=iounit, file=trim(filename), &
366 0 : action='write', status='replace', iostat=ierr)
367 0 : if (ierr /= 0) then
368 0 : write(*,*) 'failed to open ' // trim(filename)
369 0 : return
370 : end if
371 0 : call get_kap_ptr(handle,rq,ierr)
372 0 : if (ierr /= 0) then
373 0 : close(iounit)
374 0 : return
375 : end if
376 0 : call set_controls_for_writing(rq)
377 0 : write(iounit, nml=kap, iostat=ierr)
378 0 : close(iounit)
379 : end subroutine write_namelist
380 :
381 :
382 0 : subroutine set_controls_for_writing(rq)
383 : type (Kap_General_Info), pointer :: rq
384 :
385 0 : Zbase = rq% Zbase
386 :
387 0 : kap_blend_logT_upper_bdy = rq% kap_blend_logT_upper_bdy
388 0 : kap_blend_logT_lower_bdy = rq% kap_blend_logT_lower_bdy
389 :
390 0 : cubic_interpolation_in_X = rq% cubic_interpolation_in_X
391 0 : cubic_interpolation_in_Z = rq% cubic_interpolation_in_Z
392 :
393 0 : include_electron_conduction = rq% include_electron_conduction
394 0 : use_blouin_conductive_opacities = rq% use_blouin_conductive_opacities
395 :
396 0 : use_Zbase_for_Type1 = rq% use_Zbase_for_Type1
397 0 : use_Type2_opacities = rq% use_Type2_opacities
398 :
399 0 : kap_Type2_full_off_X = rq% kap_Type2_full_off_X
400 0 : kap_Type2_full_on_X = rq% kap_Type2_full_on_X
401 0 : kap_Type2_full_off_dZ = rq% kap_Type2_full_off_dZ
402 0 : kap_Type2_full_on_dZ = rq% kap_Type2_full_on_dZ
403 :
404 0 : show_info = rq% show_info
405 :
406 0 : use_other_elect_cond_opacity = rq% use_other_elect_cond_opacity
407 0 : use_other_compton_opacity = rq% use_other_compton_opacity
408 0 : use_other_radiative_opacity = rq% use_other_radiative_opacity
409 :
410 : ! user inputs
411 0 : kap_ctrl = rq% kap_ctrl
412 0 : kap_integer_ctrl = rq% kap_integer_ctrl
413 0 : kap_logical_ctrl = rq% kap_logical_ctrl
414 0 : kap_character_ctrl = rq% kap_character_ctrl
415 :
416 :
417 0 : end subroutine set_controls_for_writing
418 :
419 :
420 :
421 0 : subroutine get_kap_controls(rq, name, val, ierr)
422 : use utils_lib, only: StrUpCase
423 : type (kap_General_Info), pointer :: rq
424 : character(len=*),intent(in) :: name
425 : character(len=*), intent(out) :: val
426 : integer, intent(out) :: ierr
427 :
428 0 : character(len(name)+1) :: upper_name
429 : character(len=512) :: str
430 : integer :: iounit,iostat,ind,i
431 :
432 0 : ierr = 0
433 :
434 :
435 : ! First save current controls
436 0 : call set_controls_for_writing(rq)
437 :
438 : ! Write namelist to temporary file
439 0 : open(newunit=iounit,status='scratch')
440 0 : write(iounit,nml=kap)
441 0 : rewind(iounit)
442 :
443 : ! Namelists get written in capitals
444 0 : upper_name = trim(StrUpCase(name))//'='
445 0 : val = ''
446 : ! Search for name inside namelist
447 : do
448 0 : read(iounit,'(A)',iostat=iostat) str
449 0 : ind = index(trim(str),trim(upper_name))
450 0 : if( ind /= 0 ) then
451 0 : val = str(ind+len_trim(upper_name):len_trim(str)-1) ! Remove final comma and starting =
452 0 : do i=1,len(val)
453 0 : if(val(i:i)=='"') val(i:i) = ' '
454 : end do
455 : exit
456 : end if
457 0 : if(is_iostat_end(iostat)) exit
458 : end do
459 :
460 0 : if(len_trim(val) == 0 .and. ind==0 ) ierr = -1
461 :
462 0 : close(iounit)
463 :
464 0 : end subroutine get_kap_controls
465 :
466 0 : subroutine set_kap_controls(rq, name, val, ierr)
467 : type (kap_General_Info), pointer :: rq
468 : character(len=*), intent(in) :: name, val
469 0 : character(len=len(name)+len(val)+8) :: tmp
470 : integer, intent(out) :: ierr
471 :
472 0 : ierr = 0
473 :
474 : ! First save current kap_controls
475 0 : call set_controls_for_writing(rq)
476 :
477 0 : tmp=''
478 0 : tmp = '&kap '//trim(name)//'='//trim(val)//' /'
479 :
480 : ! Load into namelist
481 0 : read(tmp, nml=kap)
482 :
483 : ! Add to kap
484 0 : call store_controls(rq,ierr)
485 0 : if(ierr/=0) return
486 :
487 : end subroutine set_kap_controls
488 :
489 :
490 :
491 : end module kap_ctrls_io
|
