MPI-AMRVAC 3.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
Loading...
Searching...
No Matches
mod_opal_opacity.t
Go to the documentation of this file.
1!> This module reads in Rosseland-mean opacities from OPAL tables. Table opacity
2!> values are given in base 10 logarithm and are a function of mass density (R)
3!> and temperature (T), which are also both given in base 10 logarithm.
4
6
7 use mod_comm_lib, only: mpistop
8
9 implicit none
10 private
11
12 !> min and max indices for R,T-range in opacity table
13 integer, parameter :: iRmin = 2, irmax = 20
14 integer, parameter :: iTmin = 7, itmax = 76
15
16 !> If user wants to read tables from different location than AMRVAC source
17 !> NOTE: the tables should have same size as AMRVAC tables
18 logical :: set_custom_tabledir = .false.
19
20 !> The opacity tables are read once and stored globally in Kappa_vals
21 double precision, public :: kappa_vals(itmin:itmax,irmin:irmax)
22 double precision, public :: logr_list(irmin:irmax), logt_list(itmin:itmax)
23
24 public :: init_opal_table
25 public :: set_opal_opacity
26
27contains
28
29 !> This subroutine is called when the FLD radiation module is initialised.
30 !> The OPAL tables for different helium abundances are read and interpolated.
31 subroutine init_opal_table(tabledir,set_custom_tabledir)
32
33 character(len=*), intent(in) :: tabledir
34 logical, optional, intent(in) :: set_custom_tabledir
35
36 ! Local variables
37 character(len=256) :: amrvac_dir, path_table_dir
38
39 if (present(set_custom_tabledir) .and. set_custom_tabledir) then
40 path_table_dir = trim(tabledir)
41 else
42 call get_environment_variable("AMRVAC_DIR", amrvac_dir)
43 path_table_dir = trim(amrvac_dir)//"/src/tables/OPAL_tables/"//trim(tabledir)
44 endif
45
46 ! Read in the OPAL table
47 call read_table(logr_list, logt_list, kappa_vals, path_table_dir)
48
49 end subroutine init_opal_table
50
51 !> This subroutine calculates the opacity for a given temperature-density
52 !> structure. Opacities are read from a table with given metalicity.
53 subroutine set_opal_opacity(rho, temp, kappa)
54
55 double precision, intent(in) :: rho, temp
56 double precision, intent(out) :: kappa
57 ! Local variables
58 double precision :: logr_in, logt_in, logkappa_out
59
60 if (temp <= 0.0d0) call mpistop('Input temperature < 0 in set_opal_opacity.')
61 if (rho <= 0.0d0) call mpistop('Input density < 0 in set_opal_opacity.')
62
63 ! Convert input density and temperature to OPAL table convention
64 logr_in = log10(rho/(temp*1.0d-6)**3.0d0)
65 logt_in = log10(temp)
66
67 call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
68
69 ! If the outcome is 9.999 (no table entry), look right in the table
70 do while (logkappa_out > 9.0d0)
71 logr_in = logr_in + 0.5d0
72 call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
73 enddo
74
75 ! If the outcome is NaN, look left in the table
76 do while (logkappa_out <= 1.0d-4)
77 logr_in = logr_in - 0.5d0
78 call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
79 enddo
80
81 ! Convert OPAL opacity for output
82 kappa = 10.0d0**logkappa_out
83
84 end subroutine set_opal_opacity
85
86 !> This routine reads in 1-D (rho,T) values from an opacity table and gives
87 !> back as output the 1-D (rho,T) values and 2-D opacity
88 subroutine read_table(R, T, K, filename)
89
90 character(*), intent(in) :: filename
91 double precision, intent(out) :: k(itmin:itmax,irmin:irmax), r(irmin:irmax), t(itmin:itmax)
92
93 ! Local variables
94 character :: dum
95 integer :: row, col, funit=98
96 logical :: alive
97
98 inquire(file=trim(filename), exist=alive)
99
100 if (alive) then
101 open(funit, file=trim(filename), status='unknown')
102 else
103 call mpistop('Table file you want to use cannot be found: '//filename)
104 endif
105
106 ! Skip first 4 lines
107 do row = 1,4
108 read(funit,*)
109 enddo
110
111 ! Read rho
112 read(funit,*) dum, r(irmin:irmax)
113 read(funit,*)
114
115 ! Read T and kappa
116 do row = itmin,itmax
117 read(funit,'(f4.2,19f7.3)') t(row), k(row,irmin:irmax)
118 enddo
119
120 close(funit)
121
122 end subroutine read_table
123
124 !> This subroutine looks in the table for the four couples (T,rho) surrounding
125 !> a given input T and rho
126 subroutine get_kappa(Kappa_vals, logR_list, logT_list, R, T, K)
127
128 double precision, intent(in) :: kappa_vals(itmin:itmax,irmin:irmax)
129 double precision, intent(in) :: logr_list(irmin:irmax), logt_list(itmin:itmax)
130 double precision, intent(in) :: r, t
131 double precision, intent(out) :: k
132
133 ! Local variables
134 integer :: low_r_index, up_r_index
135 integer :: low_t_index, up_t_index
136
137 if (r > maxval(logr_list)) then
138 ! print*, 'Extrapolating in logR'
139 low_r_index = irmax-1
140 up_r_index = irmax
141 elseif (r < minval(logr_list)) then
142 ! print*, 'Extrapolating in logR'
143 low_r_index = irmin
144 up_r_index = irmin+1
145 else
146 call get_low_up_index(r, logr_list, irmin, irmax, low_r_index, up_r_index)
147 endif
148
149 if (t > maxval(logt_list)) then
150 ! print*, 'Extrapolating in logT'
151 low_t_index = itmax-1
152 up_t_index = itmax
153 elseif (t < minval(logt_list)) then
154 ! print*, 'Extrapolating in logT'
155 low_t_index = itmin
156 up_t_index = itmin+1
157 else
158 call get_low_up_index(t, logt_list, itmin, itmax, low_t_index, up_t_index)
159 endif
160
161 call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, logr_list, logt_list, kappa_vals, r, t, k)
162
163 end subroutine get_kappa
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, logR_list, logT_list, Kappa_vals, R, T, kappa_interp)
192
193 integer, intent(in) :: low_r, up_r, low_t, up_t
194 double precision, intent(in) :: kappa_vals(itmin:itmax,irmin:irmax)
195 double precision, intent(in) :: logr_list(irmin:irmax), logt_list(itmin:itmax)
196 double precision, intent(in) :: r, 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 = logr_list(low_r)
220 r2 = logr_list(up_r)
221 t1 = logt_list(low_t)
222 t2 = logt_list(up_t)
223
224 k1 = kappa_vals(low_t, low_r)
225 k2 = kappa_vals(low_t, up_r)
226 k3 = kappa_vals(up_t, low_r)
227 k4 = kappa_vals(up_t, up_r)
228
229 call interpolate1d(r1,r2,r,k1,k2,ka)
230 call interpolate1d(r1,r2,r,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_opal_opacity
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
This module reads in Rosseland-mean opacities from OPAL tables. Table opacity values are given in bas...
subroutine, public init_opal_table(tabledir, set_custom_tabledir)
This subroutine is called when the FLD radiation module is initialised. The OPAL tables for different...
double precision, dimension(itmin:itmax), public logt_list
double precision, dimension(itmin:itmax, irmin:irmax), public kappa_vals
The opacity tables are read once and stored globally in Kappa_vals.
subroutine, public set_opal_opacity(rho, temp, kappa)
This subroutine calculates the opacity for a given temperature-density structure. Opacities are read ...
double precision, dimension(irmin:irmax), public logr_list