13 integer,
parameter :: iRmin = 2, irmax = 20
14 integer,
parameter :: iTmin = 7, itmax = 76
18 logical :: set_custom_tabledir = .false.
21 double precision,
public ::
kappa_vals(itmin:itmax,irmin:irmax)
33 character(len=*),
intent(in) :: tabledir
34 logical,
optional,
intent(in) :: set_custom_tabledir
37 character(len=256) :: amrvac_dir, path_table_dir
39 if (
present(set_custom_tabledir) .and. set_custom_tabledir)
then
40 path_table_dir = trim(tabledir)
42 call get_environment_variable(
"AMRVAC_DIR", amrvac_dir)
43 path_table_dir = trim(amrvac_dir)//
"/src/tables/OPAL_tables/"//trim(tabledir)
55 double precision,
intent(in) :: rho, temp
56 double precision,
intent(out) :: kappa
58 double precision :: logr_in, logt_in, logkappa_out
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.')
64 logr_in = log10(rho/(temp*1.0d-6)**3.0d0)
70 do while (logkappa_out > 9.0d0)
71 logr_in = logr_in + 0.5d0
76 do while (logkappa_out <= 1.0d-4)
77 logr_in = logr_in - 0.5d0
82 kappa = 10.0d0**logkappa_out
88 subroutine read_table(R, T, K, filename)
90 character(*),
intent(in) :: filename
91 double precision,
intent(out) :: k(itmin:itmax,irmin:irmax), r(irmin:irmax), t(itmin:itmax)
95 integer :: row, col, funit=98
98 inquire(file=trim(filename), exist=alive)
101 open(funit, file=trim(filename), status=
'unknown')
103 call mpistop(
'Table file you want to use cannot be found: '//filename)
112 read(funit,*) dum, r(irmin:irmax)
117 read(funit,
'(f4.2,19f7.3)') t(row), k(row,irmin:irmax)
122 end subroutine read_table
126 subroutine get_kappa(Kappa_vals, logR_list, logT_list, R, T, K)
128 double precision,
intent(in) ::
kappa_vals(itmin:itmax,irmin:irmax)
130 double precision,
intent(in) :: r, t
131 double precision,
intent(out) :: k
134 integer :: low_r_index, up_r_index
135 integer :: low_t_index, up_t_index
139 low_r_index = irmax-1
146 call get_low_up_index(r,
logr_list, irmin, irmax, low_r_index, up_r_index)
151 low_t_index = itmax-1
158 call get_low_up_index(t,
logt_list, itmin, itmax, low_t_index, up_t_index)
163 end subroutine get_kappa
167 subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
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
174 double precision :: low_val, up_val, res(imin:imax)
176 res = var_list - var_in
179 up_val = minval(res, mask = res >= 0) + var_in
180 low_val = maxval(res, mask = res <= 0) + var_in
183 up_i = minloc(abs(var_list - up_val),1) + imin -1
184 low_i = minloc(abs(var_list - low_val),1) + imin -1
186 if (up_i == low_i) low_i = low_i - 1
188 end subroutine get_low_up_index
191 subroutine interpolate_krt(low_r, up_r, low_t, up_t, logR_list, logT_list, Kappa_vals, R, T, kappa_interp)
193 integer,
intent(in) :: low_r, up_r, low_t, up_t
194 double precision,
intent(in) ::
kappa_vals(itmin:itmax,irmin:irmax)
196 double precision,
intent(in) :: r, t
197 double precision,
intent(out) :: kappa_interp
200 double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
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)
233 end subroutine interpolate_krt
236 subroutine interpolate1d(x1, x2, x, y1, y2, y)
238 double precision,
intent(in) :: x, x1, x2, y1, y2
239 double precision,
intent(out) :: y
241 y = y1 + (x-x1)*(y2-y1)/(x2-x1)
243 end subroutine interpolate1d
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