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