MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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  implicit none
10 
11  !> min and max indices for R,T-range in opacity table
12  integer, parameter :: idmin = 2, idmax = 16
13  integer, parameter :: itmin = 2, itmax = 51
14 
15  !> The opacity tables are read once and stored globally
16  double precision, public :: alpha_vals(idmin:idmax,itmin:itmax)
17  double precision, public :: qbar_vals(idmin:idmax,itmin:itmax)
18  double precision, public :: q0_vals(idmin:idmax,itmin:itmax)
19  double precision, public :: kappae_vals(idmin:idmax,itmin:itmax)
20  double precision, public :: logd_list(idmin:idmax), logt_list(itmin:itmax)
21 
22  public :: init_cak_table
23  public :: set_cak_opacity
24 
25 contains
26 
27  !> This routine is called when the FLD radiation module is initialised.
28  subroutine init_cak_table(tabledir)
29 
30  character(len=*), intent(in) :: tabledir
31 
32  ! Local variables
33  character(len=:), allocatable :: amrvac_dir, path_table_dir
34 
35  call get_environment_variable("AMRVAC_DIR", amrvac_dir)
36  path_table_dir = trim(amrvac_dir)//"/src/rhd/CAK_tables/"//trim(tabledir)
37 
38  call read_table(logd_list, logt_list, alpha_vals, path_table_dir//"/al_TD")
39  call read_table(logd_list, logt_list, qbar_vals, path_table_dir//"/Qb_TD")
40  call read_table(logd_list, logt_list, q0_vals, path_table_dir//"/Q0_TD")
41  call read_table(logd_list, logt_list, kappae_vals, path_table_dir//"/Ke_TD")
42 
43  end subroutine init_cak_table
44 
45  !> This subroutine calculates the opacity for a given temperature-density
46  !> structure. Opacities are read from a table with given metalicity.
47  subroutine set_cak_opacity(rho, temp, alpha_output, Qbar_output, Q0_output, kappae_output)
48 
49  double precision, intent(in) :: rho, temp
50  double precision, intent(out) :: alpha_output, qbar_output, q0_output, kappae_output
51 
52  ! Local variables
53  double precision :: d_input, t_input, kappa_cak
54 
55  d_input = log10(rho)
56  t_input = log10(temp)
57 
58  d_input = min(-7.d0-1.d-5, d_input)
59  d_input = max(-16.d0+1.d-5, d_input)
60  t_input = min(5d0-1.d-5, t_input)
61  t_input = max(4d0+1.d-5, t_input)
62 
64  logd_list, logt_list, d_input, t_input, &
65  alpha_output, qbar_output, q0_output, kappae_output)
66 
67  end subroutine set_cak_opacity
68 
69  !> This routine reads in 1-D (T,rho) values from a line opacity table and gives
70  !> back as output the 1-D (T,rho) values and 2-D line opacity value
71  subroutine read_table(D, T, K, filename)
72 
73  character(*), intent(in) :: filename
74  double precision, intent(out) :: K(iDmin:iDmax,iTmin:iTmax), D(iDmin:iDmax), T(iTmin:iTmax)
75 
76  ! Local variables
77  character :: dum
78  integer :: row, col
79 
80  open(unit=1, status='old', file=trim(filename))
81 
82  ! Read temperature
83  read(1,*) dum, t(itmin:itmax)
84 
85  ! Read rho and kappa
86  do row = idmin,idmax
87  read(1,*) d(row), k(row,itmin:itmax)
88  enddo
89 
90  close(1)
91 
92  end subroutine read_table
93 
94  !> This subroutine looks in the table for the four couples (T,rho) surrounding
95  !> a given input T and rho
96  subroutine get_val_comb(K1_vals,K2_vals,K3_vals,K4_vals, &
97  logD_list, logT_list, D, T, &
98  K1, K2, K3, K4)
99 
100  double precision, intent(in) :: K1_vals(iDmin:iDmax,iTmin:iTmax)
101  double precision, intent(in) :: K2_vals(iDmin:iDmax,iTmin:iTmax)
102  double precision, intent(in) :: K3_vals(iDmin:iDmax,iTmin:iTmax)
103  double precision, intent(in) :: K4_vals(iDmin:iDmax,iTmin:iTmax)
104  double precision, intent(in) :: logD_list(iDmin:iDmax), logT_list(iTmin:iTmax)
105  double precision, intent(in) :: D, T
106  double precision, intent(out) :: K1, K2, K3, K4
107 
108  ! Local variables
109  integer :: low_D_index, up_D_index, low_T_index, up_T_index
110 
111  if (d .gt. maxval(logd_list)) then
112  ! print*, 'Extrapolating in logR'
113  low_d_index = idmax-1
114  up_d_index = idmax
115  elseif (d .lt. minval(logd_list)) then
116  ! print*, 'Extrapolating in logR'
117  low_d_index = idmin
118  up_d_index = idmin+1
119  else
120  call get_low_up_index(d, logd_list, idmin, idmax, low_d_index, up_d_index)
121  endif
122 
123  if (t .gt. maxval(logt_list)) then
124  ! print*, 'Extrapolating in logT'
125  low_t_index = itmax-1
126  up_t_index = itmax
127  elseif (t .lt. minval(logt_list)) then
128  ! print*, 'Extrapolating in logT'
129  low_t_index = itmin
130  up_t_index = itmin+1
131  else
132  call get_low_up_index(t, logt_list, itmin, itmax, low_t_index, up_t_index)
133  endif
134 
135  call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
136  logd_list, logt_list, k1_vals, d, t, k1)
137 
138  call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
139  logd_list, logt_list, k2_vals, d, t, k2)
140 
141  call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
142  logd_list, logt_list, k3_vals, d, t, k3)
143 
144  call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
145  logd_list, logt_list, k4_vals, d, t, k4)
146 
147  end subroutine get_val_comb
148 
149  !> This subroutine finds the indices in rho and T arrays of the two values
150  !> surrounding the input rho and T
151  subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
152 
153  integer, intent(in) :: imin, imax
154  double precision, intent(in) :: var_in, var_list(imin:imax)
155  integer, intent(out) :: low_i, up_i
156 
157  ! Local variables
158  double precision :: low_val, up_val, res(imin:imax)
159 
160  res = var_list - var_in
161 
162  ! Find all bounding values for given input
163  up_val = minval(res, mask = res .ge. 0) + var_in
164  low_val = maxval(res, mask = res .le. 0) + var_in
165 
166  ! Find all bounding indices
167  up_i = minloc(abs(var_list - up_val),1) + imin -1
168  low_i = minloc(abs(var_list - low_val),1) + imin -1
169 
170  if (up_i .eq. low_i) low_i = low_i - 1
171 
172  end subroutine get_low_up_index
173 
174  !> This subroutine does a bilinear interpolation in the R,T-plane
175  subroutine interpolate_krt(low_r, up_r, low_t, up_t, logD_list, logT_list, Kappa_vals, D, T, kappa_interp)
176 
177  integer, intent(in) :: low_r, up_r, low_t, up_t
178  double precision, intent(in) :: Kappa_vals(iDmin:iDmax,iTmin:iTmax)
179  double precision, intent(in) :: logD_list(iDmin:iDmax), logT_list(iTmin:iTmax)
180  double precision, intent(in) :: D, T
181  double precision, intent(out) :: kappa_interp
182 
183  ! Local variables
184  double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
185 
186  ! First interpolate twice in the T coord to get ka and kb, then interpolate
187  ! in the R coord to get ki
188 
189  ! r_1 rho r_2
190  ! | |
191  ! | |
192  ! ----k1--------ka-----k2----- t_1
193  ! | | |
194  ! | | |
195  ! T | | |
196  ! | | |
197  ! | ki |
198  ! | | |
199  ! ----k3--------kb-----k4----- t_2
200  ! | |
201  ! | |
202 
203  r1 = logd_list(low_r)
204  r2 = logd_list(up_r)
205  t1 = logt_list(low_t)
206  t2 = logt_list(up_t)
207 
208  k1 = kappa_vals(low_r, low_t)
209  k2 = kappa_vals(low_r, up_t)
210  k3 = kappa_vals(up_r, low_t)
211  k4 = kappa_vals(up_r, up_t)
212 
213  call interpolate1d(r1,r2,d,k1,k2,ka)
214  call interpolate1d(r1,r2,d,k3,k4,kb)
215  call interpolate1d(t1,t2,t,ka,kb,kappa_interp)
216 
217  end subroutine interpolate_krt
218 
219  !> Interpolation in one dimension
220  subroutine interpolate1d(x1, x2, x, y1, y2, y)
221 
222  double precision, intent(in) :: x, x1, x2, y1, y2
223  double precision, intent(out) :: y
224 
225  y = y1 + (x-x1)*(y2-y1)/(x2-x1)
226 
227  end subroutine interpolate1d
228 
229 end 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
integer, parameter idmin
min and max indices for R,T-range in opacity table
double precision, dimension(idmin:idmax, itmin:itmax), public alpha_vals
The opacity tables are read once and stored globally.
integer, parameter itmin
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 ...
subroutine get_val_comb(K1_vals, K2_vals, K3_vals, K4_vals, logD_list, logT_list, D, T, K1, K2, K3, K4)
This subroutine looks in the table for the four couples (T,rho) surrounding a given input T and rho.
integer, parameter idmax
double precision, dimension(idmin:idmax), public logd_list
subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
This subroutine finds the indices in rho and T arrays of the two values surrounding the input rho and...
subroutine interpolate_krt(low_r, up_r, low_t, up_t, logD_list, logT_list, Kappa_vals, D, T, kappa_interp)
This subroutine does a bilinear interpolation in the R,T-plane.
subroutine interpolate1d(x1, x2, x, y1, y2, y)
Interpolation in one dimension.
subroutine, public init_cak_table(tabledir)
This routine is called when the FLD radiation module is initialised.
subroutine read_table(D, T, K, filename)
This routine reads in 1-D (T,rho) values from a line opacity table and gives back as output the 1-D (...
double precision, dimension(idmin:idmax, itmin:itmax), public kappae_vals
integer, parameter itmax