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 kap_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10
24 : use auto_diff
25 :
26 : implicit none
27 :
28 : private
29 : public :: prepare_kap
30 : public :: setup_for_op_mono
31 : public :: fraction_of_op_mono
32 : public :: frac_op_mono
33 : public :: get_kap
34 :
35 : contains
36 :
37 80060 : subroutine prepare_kap (s, ierr)
38 :
39 : type(star_info), pointer :: s
40 : integer, intent(out) :: ierr
41 :
42 : ! Set opacity parameters and monochromatic stuff
43 :
44 66 : call setup_for_op_mono(s, .true., ierr)
45 : if (ierr /= 0) return
46 :
47 : end subroutine prepare_kap
48 :
49 :
50 66 : subroutine setup_for_op_mono(s, check_if_need, ierr)
51 :
52 : use kap_lib, only: load_op_mono_data, get_op_mono_params
53 : use utils_lib, only: utils_OMP_GET_MAX_THREADS
54 :
55 : type (star_info), pointer :: s
56 : logical, intent(in) :: check_if_need
57 : integer, intent(out) :: ierr
58 :
59 : integer :: k, nptot, ipe, nrad, n
60 : logical :: need_op
61 : type(auto_diff_real_2var_order1) :: beta
62 :
63 : include 'formats'
64 :
65 66 : ierr = 0
66 66 : if (s% op_mono_n > 0) return ! already setup
67 :
68 66 : if (check_if_need) then
69 66 : if (s% high_logT_op_mono_full_off < 0d0 .or. &
70 : s% low_logT_op_mono_full_off > 99d0) return
71 0 : need_op = .false.
72 0 : do k=1,s% nz
73 0 : beta = fraction_of_op_mono(s, k)
74 0 : if (beta > 0d0) then
75 : need_op = .true.
76 : exit
77 : end if
78 : end do
79 0 : if (.not. need_op) return
80 : end if
81 :
82 : call load_op_mono_data( &
83 0 : s% op_mono_data_path, s% op_mono_data_cache_filename, ierr)
84 0 : if (ierr /= 0) then
85 0 : write(*,*) 'error while loading OP data, ierr = ',ierr
86 0 : return
87 : end if
88 :
89 0 : call get_op_mono_params(nptot, ipe, nrad)
90 :
91 0 : n = utils_OMP_GET_MAX_THREADS()
92 :
93 : if (n /= s% op_mono_n .or. &
94 : nptot /= s% op_mono_nptot .or. &
95 0 : ipe /= s% op_mono_ipe .or. &
96 : nrad /= s% op_mono_nrad) then
97 0 : if (associated(s% op_mono_umesh1)) deallocate(s% op_mono_umesh1)
98 0 : if (associated(s% op_mono_semesh1)) deallocate(s% op_mono_semesh1)
99 0 : if (associated(s% op_mono_ff1)) deallocate(s% op_mono_ff1)
100 0 : if (associated(s% op_mono_rs1)) deallocate(s% op_mono_rs1)
101 0 : if (s% use_op_mono_alt_get_kap) then
102 : allocate( &
103 : s% op_mono_umesh1(nptot*n), s% op_mono_semesh1(nptot*n), s% op_mono_ff1(nptot*ipe*6*6*n), &
104 0 : s% op_mono_rs1(nptot*6*6*n), stat=ierr)
105 : else
106 : allocate( &
107 : s% op_mono_umesh1(nptot*n), s% op_mono_semesh1(nptot*n), s% op_mono_ff1(nptot*ipe*4*4*n), &
108 0 : s% op_mono_rs1(nptot*4*4*n), stat=ierr)
109 : end if
110 0 : if (ierr /= 0) return
111 0 : s% op_mono_n = n
112 0 : s% op_mono_nptot = nptot
113 0 : s% op_mono_ipe = ipe
114 0 : s% op_mono_nrad = nrad
115 : end if
116 :
117 66 : end subroutine setup_for_op_mono
118 :
119 :
120 0 : real(dp) function fraction_of_op_mono(s, k) result(beta)
121 : ! returns the real(dp) value of the blend function for cell k
122 : type (star_info), pointer :: s
123 : integer, intent(in) :: k
124 :
125 : type(auto_diff_real_2var_order1) :: frac
126 :
127 0 : if (s% op_mono_method == 'mombarg' .and. k < 2) then
128 0 : beta = 0d0 ! don't use 'mombarg' op mono method for atmospheres
129 : else
130 0 : frac = frac_op_mono(s, s% lnd(k)/ln10, s% lnT(k)/ln10)
131 : end if
132 :
133 0 : beta = frac% val
134 :
135 66 : end function fraction_of_op_mono
136 :
137 :
138 79994 : type(auto_diff_real_2var_order1) function frac_op_mono(s, logRho, logT) result(beta)
139 : ! returns an auto_diff type: var1 = lnd, var2 = lnT (derivs w.r.t. ln *not* log)
140 : use utils_lib, only: mesa_error
141 : type (star_info), pointer :: s
142 : real(dp), intent(in) :: logRho, logT
143 :
144 : real(dp) :: high_full_off, high_full_on, low_full_off, low_full_on
145 : type(auto_diff_real_2var_order1) :: alfa, log10_T
146 :
147 : include 'formats'
148 :
149 : ! default to 0
150 79994 : beta = 0
151 :
152 79994 : high_full_off = s% high_logT_op_mono_full_off
153 79994 : high_full_on = s% high_logT_op_mono_full_on
154 79994 : low_full_off = s% low_logT_op_mono_full_off
155 79994 : low_full_on = s% low_logT_op_mono_full_on
156 :
157 : if (high_full_off < high_full_on &
158 : .or. high_full_on < low_full_on &
159 79994 : .or. low_full_on < low_full_off) then
160 0 : write(*,*) "ERROR: OP mono controls must satisfy"
161 0 : write(*,*) "high_logT_op_mono_full_off >= high_logT_op_mono_full_on"
162 0 : write(*,*) "high_logT_op_mono_full_on >= low_logT_op_mono_full_on"
163 0 : write(*,*) "low_logT_op_mono_full_on >= high_logT_op_mono_full_off"
164 0 : call mesa_error(__FILE__,__LINE__,'fraction_of_op_mono')
165 : end if
166 :
167 79994 : if (high_full_off <= 0d0 .or. high_full_on <= 0d0) then
168 79994 : beta = 0._dp
169 79994 : return
170 : end if
171 :
172 : ! auto_diff version of inputs
173 0 : log10_T = logT
174 0 : log10_T % d1val2 = 1d0/ln10
175 :
176 : ! alfa is fraction standard opacity
177 :
178 0 : if (log10_T >= high_full_off .or. log10_T <= low_full_off) then
179 0 : alfa = 1d0
180 0 : else if (log10_T <= high_full_on .and. log10_T >= low_full_on) then
181 0 : alfa = 0d0
182 0 : else if (log10_T > high_full_on) then ! between high_on and high_off
183 0 : if (high_full_off - high_full_on > 1d-10) then
184 0 : alfa = (log10_T - high_full_on) / (high_full_off - high_full_on)
185 : else
186 0 : alfa = 1d0
187 : end if
188 : else ! between low_off and low_on
189 0 : if (low_full_on - low_full_off > 1d-10) then
190 0 : alfa = (log10_T - low_full_on) / (low_full_off - low_full_on)
191 : else
192 0 : alfa = 1d0
193 : end if
194 : end if
195 :
196 : ! beta is fraction of op mono
197 : ! transform linear blend to smooth quintic blend
198 0 : alfa = -alfa*alfa*alfa*(-10d0 + alfa*(15d0 - 6d0*alfa))
199 0 : beta = 1d0 - alfa
200 :
201 0 : end function frac_op_mono
202 :
203 :
204 79994 : subroutine get_kap( &
205 79994 : s, k, zbar, xa, logRho, logT, &
206 : lnfree_e, dlnfree_e_dlnRho, dlnfree_e_dlnT, &
207 : eta, d_eta_dlnRho, d_eta_dlnT, &
208 : kap_fracs, kap, dlnkap_dlnd, dlnkap_dlnT, ierr)
209 :
210 : use utils_lib
211 : use num_lib
212 : use kap_def, only: kap_test_partials, &
213 : kap_test_partials_val, kap_test_partials_dval_dx, &
214 : num_kap_fracs
215 : use kap_lib, only: &
216 : load_op_mono_data, get_op_mono_params, &
217 : get_op_mono_args, kap_get_op_mono, kap_get, &
218 : call_compute_kappa_mombarg
219 :
220 : use chem_def, only: chem_isos
221 : use star_utils, only: get_XYZ, lookup_nameofvar
222 :
223 : type (star_info), pointer :: s
224 : integer, intent(in) :: k
225 : real(dp), intent(in) :: &
226 : zbar, logRho, logT
227 : real(dp), intent(in) :: lnfree_e, dlnfree_e_dlnRho, dlnfree_e_dlnT
228 : real(dp), intent(in) :: eta, d_eta_dlnRho, d_eta_dlnT
229 : real(dp), intent(in) :: xa(:)
230 : real(dp), intent(out) :: kap_fracs(num_kap_fracs)
231 : real(dp), intent(out) :: kap, dlnkap_dlnd, dlnkap_dlnT
232 : integer, intent(out) :: ierr
233 :
234 719946 : real(dp) :: dlnkap_dxa(s% species)
235 :
236 : integer :: i, iz, nptot, ipe, nrad, thread_num, sz, offset
237 : type(auto_diff_real_2var_order1) :: beta, lnkap, lnkap_op
238 79994 : real(dp) :: kap_op, dlnkap_op_dlnRho, dlnkap_op_dlnT, &
239 79994 : Z, xh, xhe, xc, xn, xo, xne, &
240 1439892 : log_kap_rad, fk(17)
241 :
242 : character(len=4) :: e_name
243 :
244 79994 : integer, pointer :: net_iso(:)
245 : real, pointer :: &
246 79994 : umesh(:), semesh(:), ff(:,:,:,:), rs(:,:,:)
247 159988 : integer :: nel, izzp(s% species)
248 1359898 : real(dp) :: fap(s% species), fac(s% species)
249 : logical :: screening
250 : real(dp), parameter :: &
251 : eps = 1d-6, minus_eps = -eps, one_plus_eps = 1d0 + eps
252 :
253 : include 'formats'
254 :
255 79994 : if (s% doing_timing) s% timing_num_get_kap_calls = s% timing_num_get_kap_calls + 1
256 :
257 79994 : kap = -1d99
258 79994 : if (k > 0 .and. k <= s% nz) s% kap_frac_op_mono(k) = 0
259 :
260 79994 : net_iso => s% net_iso
261 79994 : xc = 0; xn = 0; xo = 0; xne = 0
262 719946 : do i=1, s% species
263 639952 : iz = chem_isos% Z(s% chem_id(i))
264 79994 : select case(iz)
265 : case (6)
266 79994 : xc = xc + xa(i)
267 : case (7)
268 79994 : xn = xn + xa(i)
269 : case (8)
270 79994 : xo = xo + xa(i)
271 : case (10)
272 639952 : xne = xne + xa(i)
273 : end select
274 : end do
275 :
276 : if (xc < minus_eps .or. xn < minus_eps .or. &
277 : xo < minus_eps .or. xne < minus_eps .or. &
278 : xc > one_plus_eps .or. xn > one_plus_eps .or. &
279 79994 : xo > one_plus_eps .or. xne > one_plus_eps) then
280 0 : ierr = -1
281 0 : if (xc < 0d0 .or. xc > 1d0) write(s% retry_message,2) 'get_kap: xc', k, xc
282 0 : if (xn < 0d0 .or. xn > 1d0) write(s% retry_message,2) 'get_kap: xn', k, xn
283 0 : if (xo < 0d0 .or. xo > 1d0) write(s% retry_message,2) 'get_kap: xo', k, xo
284 0 : if (xne < 0d0 .or. xne > 1d0) write(s% retry_message,2) 'get_kap: xne', k, xne
285 0 : if (s% report_ierr) write(*, *) s% retry_message
286 0 : return
287 : call mesa_error(__FILE__,__LINE__,'bad mass fraction: get_kap')
288 : end if
289 :
290 79994 : call get_XYZ(s, xa, xh, xhe, Z)
291 :
292 79994 : if (xh < 0d0 .or. xhe < 0d0 .or. Z < 0d0) then
293 0 : ierr = -1
294 0 : if (xh < 0d0) write(s% retry_message,2) 'xh', k, xh
295 0 : if (xhe < 0d0) write(s% retry_message,2) 'xhe', k, xhe
296 0 : if (Z < 0d0) write(s% retry_message,2) 'Z', k, Z
297 0 : if (s% report_ierr) write(*, *) s% retry_message
298 0 : return
299 : call mesa_error(__FILE__,__LINE__,'negative mass fraction: get_kap')
300 : end if
301 :
302 79994 : if (s% use_simple_es_for_kap) then
303 0 : kap = 0.2d0*(1 + xh)
304 0 : dlnkap_dlnd = 0
305 0 : dlnkap_dlnT = 0
306 0 : return
307 : end if
308 :
309 79994 : if (s% op_mono_method == 'mombarg' .and. k < 2) then
310 0 : beta = 0d0 ! don't use 'mombarg' op mono method for atmospheres
311 : else
312 79994 : beta = frac_op_mono(s, logRho, logT)
313 : end if
314 :
315 79994 : if (k > 0 .and. k <= s% nz) s% kap_frac_op_mono(k) = beta % val
316 :
317 79994 : if (beta > 0d0) then
318 :
319 0 : if (s% op_mono_method == 'hu') then
320 : call get_op_mono_args( &
321 : s% species, xa, s% op_mono_min_X_to_include, s% chem_id, &
322 0 : s% op_mono_factors, nel, izzp, fap, fac, ierr)
323 0 : if (ierr /= 0) then
324 0 : write(*,*) 'error in get_op_mono_args, ierr = ',ierr
325 0 : return
326 : end if
327 :
328 0 : if (associated(s% op_mono_umesh1)) then
329 :
330 0 : thread_num = utils_OMP_GET_THREAD_NUM() ! in range 0 to op_mono_n-1
331 0 : if (thread_num < 0) then
332 0 : write(*,3) 'thread_num < 0', thread_num, s% op_mono_n
333 0 : ierr = -1
334 0 : return
335 : end if
336 0 : if (thread_num >= s% op_mono_n) then
337 0 : write(*,3) 'thread_num >= s% op_mono_n', thread_num, s% op_mono_n
338 0 : ierr = -1
339 0 : return
340 : end if
341 0 : nptot = s% op_mono_nptot
342 0 : ipe = s% op_mono_ipe
343 0 : nrad = s% op_mono_nrad
344 :
345 0 : sz = nptot; offset = thread_num*sz
346 0 : umesh(1:nptot) => s% op_mono_umesh1(offset+1:offset+sz)
347 0 : semesh(1:nptot) => s% op_mono_semesh1(offset+1:offset+sz)
348 0 : if (s% use_op_mono_alt_get_kap) then
349 0 : sz = nptot*ipe*6*6; offset = thread_num*sz
350 0 : ff(1:nptot,1:ipe,1:6,1:6) => s% op_mono_ff1(offset+1:offset+sz)
351 0 : sz = nptot*6*6; offset = thread_num*sz
352 0 : rs(1:nptot,1:6,1:6) => s% op_mono_rs1(offset+1:offset+sz)
353 0 : sz = nptot*nrad*6*6; offset = thread_num*sz
354 : else
355 0 : sz = nptot*ipe*4*4; offset = thread_num*sz
356 0 : ff(1:nptot,1:ipe,1:4,1:4) => s% op_mono_ff1(offset+1:offset+sz)
357 0 : sz = nptot*4*4; offset = thread_num*sz
358 0 : rs(1:nptot,1:4,1:4) => s% op_mono_rs1(offset+1:offset+sz)
359 0 : sz = nptot*nrad*4*4; offset = thread_num*sz
360 : end if
361 :
362 : else
363 :
364 : call load_op_mono_data( &
365 0 : s% op_mono_data_path, s% op_mono_data_cache_filename, ierr)
366 0 : if (ierr /= 0) then
367 0 : write(*,*) 'error while loading OP data, ierr = ',ierr
368 0 : return
369 : end if
370 :
371 0 : call get_op_mono_params(nptot, ipe, nrad)
372 0 : if (s% use_op_mono_alt_get_kap) then
373 : allocate( &
374 : umesh(nptot), semesh(nptot), ff(nptot,ipe,6,6), &
375 0 : rs(nptot,6,6), stat=ierr)
376 : else
377 : allocate( &
378 : umesh(nptot), semesh(nptot), ff(nptot,ipe,4,4), &
379 0 : rs(nptot,4,4), stat=ierr)
380 : end if
381 0 : if (ierr /= 0) return
382 :
383 : end if
384 :
385 0 : if (s% solver_test_kap_partials) then
386 : kap_test_partials = (k == s% solver_test_partials_k .and. &
387 : s% solver_call_number == s% solver_test_partials_call_number .and. &
388 0 : s% solver_iter == s% solver_test_partials_iter_number )
389 : end if
390 :
391 0 : screening = .true.
392 0 : if (s% use_other_kap) then
393 : call s% other_kap_get_op_mono( &
394 : s% kap_handle, zbar, logRho, logT, &
395 : s% use_op_mono_alt_get_kap, &
396 : nel, izzp, fap, fac, screening, umesh, semesh, ff, rs, &
397 0 : kap_op, dlnkap_op_dlnRho, dlnkap_op_dlnT, ierr)
398 : else
399 : call kap_get_op_mono( &
400 : s% kap_handle, zbar, logRho, logT, &
401 : s% use_op_mono_alt_get_kap, &
402 : nel, izzp, fap, fac, screening, umesh, semesh, ff, rs, &
403 0 : kap_op, dlnkap_op_dlnRho, dlnkap_op_dlnT, ierr)
404 : end if
405 :
406 0 : if (s% solver_test_kap_partials .and. kap_test_partials) then
407 0 : s% solver_test_partials_val = kap_test_partials_val
408 0 : s% solver_test_partials_dval_dx = kap_test_partials_dval_dx
409 : end if
410 :
411 0 : if (.not. associated(s% op_mono_umesh1)) deallocate(umesh, semesh, ff, rs)
412 :
413 0 : else if (s% op_mono_method == 'mombarg') then
414 0 : fk = 0
415 0 : if (logT > 3.5d0 .and. logT < 8.0d0) then
416 0 : do i=1, s% species
417 0 : e_name = chem_isos% name(s% chem_id(i))
418 0 : if (e_name == 'h1') fk(1) = xa(i)/ chem_isos% W(s% chem_id(i))
419 0 : if (e_name == 'he4') fk(2) = xa(i)/ chem_isos% W(s% chem_id(i))
420 0 : if (e_name == 'c12') fk(3) = xa(i)/ chem_isos% W(s% chem_id(i))
421 0 : if (e_name == 'n14') fk(4) = xa(i)/ chem_isos% W(s% chem_id(i))
422 0 : if (e_name == 'o16') fk(5) = xa(i)/ chem_isos% W(s% chem_id(i))
423 0 : if (e_name == 'ne20')fk(6) = xa(i)/ chem_isos% W(s% chem_id(i))
424 0 : if (e_name == 'na23')fk(7) = xa(i)/ chem_isos% W(s% chem_id(i))
425 0 : if (e_name == 'mg24')fk(8) = xa(i)/ chem_isos% W(s% chem_id(i))
426 0 : if (e_name == 'al27')fk(9) = xa(i)/ chem_isos% W(s% chem_id(i))
427 0 : if (e_name == 'si28')fk(10) = xa(i)/ chem_isos% W(s% chem_id(i))
428 0 : if (e_name == 's32') fk(11) = xa(i)/ chem_isos% W(s% chem_id(i))
429 0 : if (e_name == 'ar40')fk(12) = xa(i)/ chem_isos% W(s% chem_id(i))
430 0 : if (e_name == 'ca40')fk(13) = xa(i)/ chem_isos% W(s% chem_id(i))
431 0 : if (e_name == 'cr52')fk(14) = xa(i)/ chem_isos% W(s% chem_id(i))
432 0 : if (e_name == 'mn55')fk(15) = xa(i)/ chem_isos% W(s% chem_id(i))
433 0 : if (e_name == 'fe56')fk(16) = xa(i)/ chem_isos% W(s% chem_id(i))
434 0 : if (e_name == 'ni58')fk(17) = xa(i)/ chem_isos% W(s% chem_id(i))
435 : end do
436 0 : fk = fk / sum(fk)
437 : call call_compute_kappa_mombarg(s% kap_handle, k, &
438 : fk, logT, logRho, &
439 : zbar, lnfree_e, dlnfree_e_dlnRho, dlnfree_e_dlnT, &
440 0 : kap_op, dlnkap_op_dlnT, dlnkap_op_dlnRho, log_kap_rad, ierr)
441 : end if
442 : else
443 0 : write(*,*) 'Invalid argument for op_mono_method.'
444 0 : stop
445 : end if
446 :
447 0 : if (ierr /= 0) then
448 : ! Fall back to standard opacities, not OP
449 0 : beta = 0
450 0 : if (k > 0 .and. k <= s% nz) s% kap_frac_op_mono(k) = 0
451 0 : ierr = 0
452 0 : else if (beta == 1d0) then
453 0 : kap = kap_op
454 0 : dlnkap_dlnT = dlnkap_op_dlnT
455 0 : dlnkap_dlnd = dlnkap_op_dlnRho
456 0 : return
457 : end if
458 :
459 0 : if (is_bad(kap_op)) then
460 0 : ierr = 1
461 0 : return
462 : end if
463 :
464 0 : if (kap_op <= 0d0) then
465 0 : ierr = 1
466 0 : return
467 : end if
468 :
469 0 : lnkap_op = log(kap_op)
470 0 : lnkap_op% d1val1 = dlnkap_op_dlnRho
471 0 : lnkap_op% d1val2 = dlnkap_op_dlnT
472 :
473 : end if
474 :
475 79994 : if (s% solver_test_kap_partials) then
476 : kap_test_partials = (k == s% solver_test_partials_k .and. &
477 : s% solver_call_number == s% solver_test_partials_call_number .and. &
478 0 : s% solver_iter == s% solver_test_partials_iter_number )
479 : end if
480 :
481 79994 : if (s% use_other_kap) then
482 : call s% other_kap_get( &
483 : s% id, k, s% kap_handle, s% species, s% chem_id, s% net_iso, xa, &
484 : logRho, logT, &
485 : lnfree_e, dlnfree_e_dlnRho, dlnfree_e_dlnT, &
486 : eta, d_eta_dlnRho, d_eta_dlnT, &
487 0 : kap_fracs, kap, dlnkap_dlnd, dlnkap_dlnT, dlnkap_dxa, ierr)
488 : else
489 : call kap_get( &
490 : s% kap_handle, s% species, s% chem_id, s% net_iso, xa, &
491 : logRho, logT, &
492 : lnfree_e, dlnfree_e_dlnRho, dlnfree_e_dlnT, &
493 : eta, d_eta_dlnRho, d_eta_dlnT, &
494 79994 : kap_fracs, kap, dlnkap_dlnd, dlnkap_dlnT, dlnkap_dxa, ierr)
495 : end if
496 :
497 79994 : if (ierr /= 0) then
498 : return
499 : end if
500 :
501 79994 : if (s% solver_test_kap_partials .and. kap_test_partials) then
502 0 : s% solver_test_partials_val = kap_test_partials_val
503 0 : s% solver_test_partials_dval_dx = kap_test_partials_dval_dx
504 : end if
505 :
506 79994 : if (beta == 0d0) then
507 : return
508 : end if
509 :
510 : ! OP mono already packed in lnkap_op
511 : ! pack lnkap with standard opacities
512 0 : lnkap = log(kap)
513 0 : lnkap% d1val1 = dlnkap_dlnd
514 0 : lnkap% d1val2 = dlnkap_dlnT
515 :
516 : ! Blend opacities
517 0 : lnkap = (1d0-beta)*lnkap + beta*lnkap_op
518 :
519 : ! unpack autodiff type
520 0 : kap = exp(lnkap% val)
521 0 : dlnkap_dlnd = lnkap% d1val1
522 0 : dlnkap_dlnT = lnkap% d1val2
523 :
524 159988 : end subroutine get_kap
525 :
526 : end module kap_support
|
