MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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  implicit none
8 
9  !> min and max indices for R,T-range in opacity table
10  integer, parameter :: irmin = 2, irmax = 20
11  integer, parameter :: itmin = 7, itmax = 76
12 
13  !> The opacity tables are read once and stored globally in Kappa_vals
14  double precision, public :: kappa_vals(itmin:itmax,irmin:irmax)
15  double precision, public :: logr_list(irmin:irmax), logt_list(itmin:itmax)
16 
17  public :: init_opal_table
18  public :: set_opal_opacity
19 
20 contains
21 
22  !> This subroutine is called when the FLD radiation module is initialised.
23  !> The OPAL tables for different helium abundances are read and interpolated.
24  subroutine init_opal_table(tablename)
25 
26  character(len=*), intent(in) :: tablename
27 
28  ! Local variables
29  character(len=:), allocatable :: amrvac_dir, path_table_dir
30 
31  call get_environment_variable("AMRVAC_DIR", amrvac_dir)
32 
33  ! Absolute path to OPAL table location
34  path_table_dir = trim(amrvac_dir)//"/src/rhd/OPAL_tables/"//trim(tablename)
35 
36  ! Read in the OPAL table
37  call read_table(logr_list, logt_list, kappa_vals, path_table_dir)
38 
39  end subroutine init_opal_table
40 
41  !> This subroutine calculates the opacity for a given temperature-density
42  !> structure. Opacities are read from a table with given metalicity.
43  subroutine set_opal_opacity(rho, temp, kappa)
44 
45  double precision, intent(in) :: rho, temp
46  double precision, intent(out) :: kappa
47 
48  ! Local variables
49  double precision :: logr_in, logt_in, logkappa_out
50 
51  ! Convert input density and temperature to OPAL table convention
52  logr_in = log10(rho/(temp*1d-6)**3)
53  logt_in = log10(temp)
54 
55  call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
56 
57  ! If the outcome is 9.999 (no table entry), look right in the table
58  do while (logkappa_out .gt. 9.0d0)
59  logr_in = logr_in + 0.5d0
60  call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
61  enddo
62 
63  ! If the outcome is NaN, look left in the table
64  do while (logkappa_out .eq. 0.0d0)
65  logr_in = logr_in - 0.5d0
66  call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
67  enddo
68 
69  ! Convert OPAL opacity for output
70  kappa = 10d0**logkappa_out
71 
72  end subroutine set_opal_opacity
73 
74  !> This routine reads in 1-D (rho,T) values from an opacity table and gives
75  !> back as output the 1-D (rho,T) values and 2-D opacity
76  subroutine read_table(R, T, K, filename)
77 
78  character(*), intent(in) :: filename
79  double precision, intent(out) :: K(iTmin:iTmax,iRmin:iRmax), R(iRmin:iRmax), T(iTmin:iTmax)
80 
81  ! Local variables
82  character :: dum
83  integer :: row, col
84 
85  open(unit=1, status='old', file=trim(filename))
86 
87  ! Skip first 4 lines
88  do row = 1,4
89  read(1,*)
90  enddo
91 
92  ! Read rho
93  read(1,*) dum, r(irmin:irmax)
94  read(1,*)
95 
96  ! Read T and kappa
97  do row = itmin,itmax
98  read(1,'(f4.2,19f7.3)') t(row), k(row,irmin:irmax)
99  enddo
100 
101  close(1)
102 
103  end subroutine read_table
104 
105  !> This subroutine looks in the table for the four couples (T,rho) surrounding
106  !> a given input T and rho
107  subroutine get_kappa(Kappa_vals, logR_list, logT_list, R, T, K)
108 
109  double precision, intent(in) :: Kappa_vals(iTmin:iTmax,iRmin:iRmax)
110  double precision, intent(in) :: logR_list(iRmin:iRmax), logT_list(iTmin:iTmax)
111  double precision, intent(in) :: R, T
112  double precision, intent(out) :: K
113 
114  ! Local variables
115  integer :: low_R_index, up_R_index
116  integer :: low_T_index, up_T_index
117 
118  if (r .gt. maxval(logr_list)) then
119  ! print*, 'Extrapolating in logR'
120  low_r_index = irmax-1
121  up_r_index = irmax
122  elseif (r .lt. minval(logr_list)) then
123  ! print*, 'Extrapolating in logR'
124  low_r_index = irmin
125  up_r_index = irmin+1
126  else
127  call get_low_up_index(r, logr_list, irmin, irmax, low_r_index, up_r_index)
128  endif
129 
130  if (t .gt. maxval(logt_list)) then
131  ! print*, 'Extrapolating in logT'
132  low_t_index = itmax-1
133  up_t_index = itmax
134  elseif ( t .lt. minval(logt_list)) then
135  ! print*, 'Extrapolating in logT'
136  low_t_index = itmin
137  up_t_index = itmin+1
138  else
139  call get_low_up_index(t, logt_list, itmin, itmax, low_t_index, up_t_index)
140  endif
141 
142  call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, logr_list, logt_list, kappa_vals, r, t, k)
143 
144  end subroutine get_kappa
145 
146  !> This subroutine finds the indices in rho and T arrays of the two values
147  !> surrounding the input rho and T
148  subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
149 
150  integer, intent(in) :: imin, imax
151  double precision, intent(in) :: var_in, var_list(imin:imax)
152  integer, intent(out) :: low_i, up_i
153 
154  ! Local variables
155  double precision :: low_val, up_val, res(imin:imax)
156 
157  res = var_list - var_in
158 
159  ! Find all bounding values for given input
160  up_val = minval(res, mask = res .ge. 0) + var_in
161  low_val = maxval(res, mask = res .le. 0) + var_in
162 
163  ! Find all bounding indices
164  up_i = minloc(abs(var_list - up_val),1) + imin -1
165  low_i = minloc(abs(var_list - low_val),1) + imin -1
166 
167  if (up_i .eq. low_i) low_i = low_i - 1
168 
169  end subroutine get_low_up_index
170 
171  !> This subroutine does a bilinear interpolation in the R,T-plane
172  subroutine interpolate_krt(low_r, up_r, low_t, up_t, logR_list, logT_list, Kappa_vals, R, T, kappa_interp)
173 
174  integer, intent(in) :: low_r, up_r, low_t, up_t
175  double precision, intent(in) :: Kappa_vals(iTmin:iTmax,iRmin:iRmax)
176  double precision, intent(in) :: logR_list(iRmin:iRmax), logT_list(iTmin:iTmax)
177  double precision, intent(in) :: R, T
178  double precision, intent(out) :: kappa_interp
179 
180  ! Local variables
181  double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
182 
183  ! First interpolate twice in the T coord to get ka and kb, then interpolate
184  ! in the R coord to get ki
185 
186  ! r_1 rho r_2
187  ! | |
188  ! | |
189  ! ----k1--------ka-----k2----- t_1
190  ! | | |
191  ! | | |
192  ! T | | |
193  ! | | |
194  ! | ki |
195  ! | | |
196  ! ----k3--------kb-----k4----- t_2
197  ! | |
198  ! | |
199 
200  r1 = logr_list(low_r)
201  r2 = logr_list(up_r)
202  t1 = logt_list(low_t)
203  t2 = logt_list(up_t)
204 
205  k1 = kappa_vals(low_t, low_r)
206  k2 = kappa_vals(low_t, up_r)
207  k3 = kappa_vals(up_t, low_r)
208  k4 = kappa_vals(up_t, up_r)
209 
210  call interpolate1d(r1,r2,r,k1,k2,ka)
211  call interpolate1d(r1,r2,r,k3,k4,kb)
212  call interpolate1d(t1,t2,t,ka,kb,kappa_interp)
213 
214  end subroutine interpolate_krt
215 
216  !> Interpolation in one dimension
217  subroutine interpolate1d(x1, x2, x, y1, y2, y)
218 
219  double precision, intent(in) :: x, x1, x2, y1, y2
220  double precision, intent(out) :: y
221 
222  y = y1 + (x-x1)*(y2-y1)/(x2-x1)
223 
224  end subroutine interpolate1d
225 
226 end module mod_opal_opacity
This module reads in Rosseland-mean opacities from OPAL tables. Table opacity values are given in bas...
integer, parameter irmin
min and max indices for R,T-range in opacity table
double precision, dimension(itmin:itmax), public logt_list
subroutine, public init_opal_table(tablename)
This subroutine is called when the FLD radiation module is initialised. The OPAL tables for different...
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 ...
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 interpolate1d(x1, x2, x, y1, y2, y)
Interpolation in one dimension.
subroutine get_kappa(Kappa_vals, logR_list, logT_list, R, T, K)
This subroutine looks in the table for the four couples (T,rho) surrounding a given input T and rho.
integer, parameter irmax
integer, parameter itmin
subroutine interpolate_krt(low_r, up_r, low_t, up_t, logR_list, logT_list, Kappa_vals, R, T, kappa_interp)
This subroutine does a bilinear interpolation in the R,T-plane.
double precision, dimension(irmin:irmax), public logr_list
integer, parameter itmax
subroutine read_table(R, T, K, filename)
This routine reads in 1-D (rho,T) values from an opacity table and gives back as output the 1-D (rho,...