MPI-AMRVAC 3.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
Loading...
Searching...
No Matches
mod_cak_opacity.t
Go to the documentation of this file.
1!> This module reads in CAK line opacities in the Gayley (1995) notation
2!> (alpha, Qbar, Q0, kappae) from corresponding tables. Tabulated values assume
3!> LTE conditions and are a function of mass density (D) and temperature (T),
4!> which are both given in base 10 logarithm.
5!> The construction of the tables is outlined in Poniatowski+ (2021), A&A, 667.
6
8
9 use mod_comm_lib, only: mpistop
10
11 implicit none
12 private
13
14 !> min and max indices for R,T-range in opacity table
15 integer, parameter :: iDmin = 2, idmax = 51
16 integer, parameter :: iTmin = 2, itmax = 21
17
18 !> If user wants to read tables from different location than AMRVAC source
19 !> NOTE: the tables should have same size as AMRVAC tables
20 logical :: set_custom_tabledir = .false.
21
22 !> The opacity tables are read once and stored globally
23 double precision, public :: alpha_vals(idmin:idmax,itmin:itmax)
24 double precision, public :: qbar_vals(idmin:idmax,itmin:itmax)
25 double precision, public :: q0_vals(idmin:idmax,itmin:itmax)
26 double precision, public :: kappae_vals(idmin:idmax,itmin:itmax)
27 double precision, public :: logd_list(idmin:idmax), logt_list(itmin:itmax)
28
29 public :: init_cak_table
30 public :: set_cak_opacity
31
32contains
33
34 !> This routine is called when the FLD radiation module is initialised.
35 subroutine init_cak_table(tabledir,set_custom_tabledir)
36
37 character(len=*), intent(in) :: tabledir
38 logical, optional, intent(in) :: set_custom_tabledir
39
40 ! Local variables
41 character(len=256) :: amrvac_dir, path_table_dir
42
43 if (present(set_custom_tabledir) .and. set_custom_tabledir) then
44 path_table_dir = trim(tabledir)
45 else
46 call get_environment_variable("AMRVAC_DIR", amrvac_dir)
47 path_table_dir = trim(amrvac_dir)//"/src/tables/CAK_tables/"//trim(tabledir)
48 endif
49
50 call read_table(logd_list, logt_list, alpha_vals, trim(path_table_dir)//"/al_TD")
51 call read_table(logd_list, logt_list, qbar_vals, trim(path_table_dir)//"/Qb_TD")
52 call read_table(logd_list, logt_list, q0_vals, trim(path_table_dir)//"/Q0_TD")
53 call read_table(logd_list, logt_list, kappae_vals, trim(path_table_dir)//"/Ke_TD")
54
55 end subroutine init_cak_table
56
57 !> This subroutine calculates the opacity for a given temperature-density
58 !> structure. Opacities are read from a table with given metalicity.
59 subroutine set_cak_opacity(rho, temp, alpha_output, Qbar_output, Q0_output, kappae_output)
60
61 double precision, intent(in) :: rho, temp
62 double precision, intent(out) :: alpha_output, qbar_output, q0_output, kappae_output
63
64 ! Local variables
65 double precision :: d_input, t_input, kappa_cak
66
67 if (temp <= 0.0d0) call mpistop('Input temperature < 0 in set_cak_opacity.')
68 if (rho <= 0.0d0) call mpistop('Input density < 0 in set_cak_opacity.')
69
70 d_input = log10(rho)
71 t_input = log10(temp)
72
73 call get_val_comb(alpha_vals,qbar_vals,q0_vals,kappae_vals, &
74 logd_list, logt_list, d_input, t_input, &
75 alpha_output, qbar_output, q0_output, kappae_output)
76
77 end subroutine set_cak_opacity
78
79 !> This routine reads in 1-D (T,rho) values from a line opacity table and gives
80 !> back as output the 1-D (T,rho) values and 2-D line opacity value
81 subroutine read_table(D, T, K, filename)
82
83 character(*), intent(in) :: filename
84 double precision, intent(out) :: k(idmin:idmax,itmin:itmax), d(idmin:idmax), t(itmin:itmax)
85
86 ! Local variables
87 character :: dum
88 integer :: row, col, funit=99
89 logical :: alive
90
91 inquire(file=trim(filename), exist=alive)
92
93 if (alive) then
94 open(funit, file=trim(filename), status='unknown')
95 else
96 call mpistop('Table file you want to use cannot be found: '//filename)
97 endif
98
99 ! Read temperature
100 read(funit,*) dum, t(itmin:itmax)
101
102 ! Read rho and kappa
103 do row = idmin,idmax
104 read(funit,*) d(row), k(row,itmin:itmax)
105 enddo
106
107 close(funit)
108
109 end subroutine read_table
110
111 !> This subroutine looks in the table for the four couples (T,rho) surrounding
112 !> a given input T and rho
113 subroutine get_val_comb(K1_vals,K2_vals,K3_vals,K4_vals, &
114 logD_list, logT_list, D, T, K1, K2, K3, K4)
115
116 double precision, intent(in) :: k1_vals(idmin:idmax,itmin:itmax)
117 double precision, intent(in) :: k2_vals(idmin:idmax,itmin:itmax)
118 double precision, intent(in) :: k3_vals(idmin:idmax,itmin:itmax)
119 double precision, intent(in) :: k4_vals(idmin:idmax,itmin:itmax)
120 double precision, intent(in) :: logd_list(idmin:idmax), logt_list(itmin:itmax)
121 double precision, intent(in) :: d, t
122 double precision, intent(out) :: k1, k2, k3, k4
123
124 ! Local variables
125 integer :: low_d_index, up_d_index, low_t_index, up_t_index
126
127 if (d > maxval(logd_list)) then
128 ! print*, 'Extrapolating in logR'
129 low_d_index = idmax-1
130 up_d_index = idmax
131 elseif (d < minval(logd_list)) then
132 ! print*, 'Extrapolating in logR'
133 low_d_index = idmin
134 up_d_index = idmin+1
135 else
136 call get_low_up_index(d, logd_list, idmin, idmax, low_d_index, up_d_index)
137 endif
138
139 if (t > maxval(logt_list)) then
140 ! print*, 'Extrapolating in logT'
141 low_t_index = itmax-1
142 up_t_index = itmax
143 elseif (t < minval(logt_list)) then
144 ! print*, 'Extrapolating in logT'
145 low_t_index = itmin
146 up_t_index = itmin+1
147 else
148 call get_low_up_index(t, logt_list, itmin, itmax, low_t_index, up_t_index)
149 endif
150
151 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
152 logd_list, logt_list, k1_vals, d, t, k1)
153
154 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
155 logd_list, logt_list, k2_vals, d, t, k2)
156
157 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
158 logd_list, logt_list, k3_vals, d, t, k3)
159
160 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
161 logd_list, logt_list, k4_vals, d, t, k4)
162
163 end subroutine get_val_comb
164
165 !> This subroutine finds the indices in rho and T arrays of the two values
166 !> surrounding the input rho and T
167 subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
168
169 integer, intent(in) :: imin, imax
170 double precision, intent(in) :: var_in, var_list(imin:imax)
171 integer, intent(out) :: low_i, up_i
172
173 ! Local variables
174 double precision :: low_val, up_val, res(imin:imax)
175
176 res = var_list - var_in
177
178 ! Find all bounding values for given input
179 up_val = minval(res, mask = res >= 0) + var_in
180 low_val = maxval(res, mask = res <= 0) + var_in
181
182 ! Find all bounding indices
183 up_i = minloc(abs(var_list - up_val),1) + imin -1
184 low_i = minloc(abs(var_list - low_val),1) + imin -1
185
186 if (up_i == low_i) low_i = low_i - 1
187
188 end subroutine get_low_up_index
189
190 !> This subroutine does a bilinear interpolation in the R,T-plane
191 subroutine interpolate_krt(low_r, up_r, low_t, up_t, logD_list, logT_list, Kappa_vals, D, T, kappa_interp)
192
193 integer, intent(in) :: low_r, up_r, low_t, up_t
194 double precision, intent(in) :: kappa_vals(idmin:idmax,itmin:itmax)
195 double precision, intent(in) :: logd_list(idmin:idmax), logt_list(itmin:itmax)
196 double precision, intent(in) :: d, t
197 double precision, intent(out) :: kappa_interp
198
199 ! Local variables
200 double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
201
202 ! First interpolate twice in the T coord to get ka and kb, then interpolate
203 ! in the R coord to get ki
204
205 ! r_1 rho r_2
206 ! | |
207 ! | |
208 ! ----k1--------ka-----k2----- t_1
209 ! | | |
210 ! | | |
211 ! T | | |
212 ! | | |
213 ! | ki |
214 ! | | |
215 ! ----k3--------kb-----k4----- t_2
216 ! | |
217 ! | |
218
219 r1 = logd_list(low_r)
220 r2 = logd_list(up_r)
221 t1 = logt_list(low_t)
222 t2 = logt_list(up_t)
223
224 k1 = kappa_vals(low_r, low_t)
225 k2 = kappa_vals(low_r, up_t)
226 k3 = kappa_vals(up_r, low_t)
227 k4 = kappa_vals(up_r, up_t)
228
229 call interpolate1d(r1,r2,d,k1,k2,ka)
230 call interpolate1d(r1,r2,d,k3,k4,kb)
231 call interpolate1d(t1,t2,t,ka,kb,kappa_interp)
232
233 end subroutine interpolate_krt
234
235 !> Interpolation in one dimension
236 subroutine interpolate1d(x1, x2, x, y1, y2, y)
237
238 double precision, intent(in) :: x, x1, x2, y1, y2
239 double precision, intent(out) :: y
240
241 y = y1 + (x-x1)*(y2-y1)/(x2-x1)
242
243 end subroutine interpolate1d
244
245end module mod_cak_opacity
This module reads in CAK line opacities in the Gayley (1995) notation (alpha, Qbar,...
double precision, dimension(itmin:itmax), public logt_list
double precision, dimension(idmin:idmax, itmin:itmax), public q0_vals
double precision, dimension(idmin:idmax, itmin:itmax), public qbar_vals
subroutine, public init_cak_table(tabledir, set_custom_tabledir)
This routine is called when the FLD radiation module is initialised.
double precision, dimension(idmin:idmax, itmin:itmax), public alpha_vals
The opacity tables are read once and stored globally.
double precision, dimension(idmin:idmax), public logd_list
subroutine, public set_cak_opacity(rho, temp, alpha_output, qbar_output, q0_output, kappae_output)
This subroutine calculates the opacity for a given temperature-density structure. Opacities are read ...
double precision, dimension(idmin:idmax, itmin:itmax), public kappae_vals
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.