15 integer,
parameter :: iDmin = 2, idmax = 51
16 integer,
parameter :: iTmin = 2, itmax = 21
20 logical :: set_custom_tabledir = .false.
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)
37 character(len=*),
intent(in) :: tabledir
38 logical,
optional,
intent(in) :: set_custom_tabledir
41 character(len=256) :: amrvac_dir, path_table_dir
43 if (
present(set_custom_tabledir) .and. set_custom_tabledir)
then
44 path_table_dir = trim(tabledir)
46 call get_environment_variable(
"AMRVAC_DIR", amrvac_dir)
47 path_table_dir = trim(amrvac_dir)//
"/src/tables/CAK_tables/"//trim(tabledir)
59 subroutine set_cak_opacity(rho, temp, alpha_output, Qbar_output, Q0_output, kappae_output)
61 double precision,
intent(in) :: rho, temp
62 double precision,
intent(out) :: alpha_output, qbar_output, q0_output, kappae_output
65 double precision :: d_input, t_input, kappa_cak
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.')
75 alpha_output, qbar_output, q0_output, kappae_output)
81 subroutine read_table(D, T, K, filename)
83 character(*),
intent(in) :: filename
84 double precision,
intent(out) :: k(idmin:idmax,itmin:itmax), d(idmin:idmax), t(itmin:itmax)
88 integer :: row, col, funit=99
91 inquire(file=trim(filename), exist=alive)
94 open(funit, file=trim(filename), status=
'unknown')
96 call mpistop(
'Table file you want to use cannot be found: '//filename)
100 read(funit,*) dum, t(itmin:itmax)
104 read(funit,*) d(row), k(row,itmin:itmax)
109 end subroutine read_table
113 subroutine get_val_comb(K1_vals,K2_vals,K3_vals,K4_vals, &
114 logD_list, logT_list, D, T, K1, K2, K3, K4)
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)
121 double precision,
intent(in) :: d, t
122 double precision,
intent(out) :: k1, k2, k3, k4
125 integer :: low_d_index, up_d_index, low_t_index, up_t_index
129 low_d_index = idmax-1
136 call get_low_up_index(d,
logd_list, idmin, idmax, low_d_index, up_d_index)
141 low_t_index = itmax-1
148 call get_low_up_index(t,
logt_list, itmin, itmax, low_t_index, up_t_index)
151 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
154 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
157 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
160 call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
163 end subroutine get_val_comb
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, logD_list, logT_list, Kappa_vals, D, T, kappa_interp)
193 integer,
intent(in) :: low_r, up_r, low_t, up_t
194 double precision,
intent(in) :: kappa_vals(idmin:idmax,itmin:itmax)
196 double precision,
intent(in) :: d, t
197 double precision,
intent(out) :: kappa_interp
200 double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
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)
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)
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
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.