mod_cak_opacity.t

Go to the documentation of this file.

!> This module reads in CAK line opacities in the Gayley (1995) notation
!> (alpha, Qbar, Q0, kappae) from corresponding tables. Tabulated values assume
!> LTE conditions and are a function of mass density (D) and temperature (T),
!> which are both given in base 10 logarithm.
!> The construction of the tables is outlined in Poniatowski+ (2021), A&A, 667.
 
module mod_cak_opacity
 
  implicit none
 
  !> min and max indices for R,T-range in opacity table
  integer, parameter :: idmin = 2, idmax = 16
  integer, parameter :: itmin = 2, itmax = 51
 
  !> The opacity tables are read once and stored globally
  double precision, public :: alpha_vals(idmin:idmax,itmin:itmax)
  double precision, public :: qbar_vals(idmin:idmax,itmin:itmax)
  double precision, public :: q0_vals(idmin:idmax,itmin:itmax)
  double precision, public :: kappae_vals(idmin:idmax,itmin:itmax)
  double precision, public :: logd_list(idmin:idmax), logt_list(itmin:itmax)
 
  public :: init_cak_table
  public :: set_cak_opacity
 
contains
 
  !> This routine is called when the FLD radiation module is initialised.
  subroutine init_cak_table(tabledir)
 
    character(len=*), intent(in) :: tabledir
 
    ! Local variables
    character(len=:), allocatable :: amrvac_dir, path_table_dir
    
    call get_environment_variable("AMRVAC_DIR", amrvac_dir)
    path_table_dir = trim(amrvac_dir)//"/src/rhd/CAK_tables/"//trim(tabledir)
 
    call read_table(logd_list, logt_list, alpha_vals, path_table_dir//"/al_TD")
    call read_table(logd_list, logt_list, qbar_vals, path_table_dir//"/Qb_TD")
    call read_table(logd_list, logt_list, q0_vals, path_table_dir//"/Q0_TD")
    call read_table(logd_list, logt_list, kappae_vals, path_table_dir//"/Ke_TD")
 
  end subroutine init_cak_table
 
  !> This subroutine calculates the opacity for a given temperature-density
  !> structure. Opacities are read from a table with given metalicity.
  subroutine set_cak_opacity(rho, temp, alpha_output, Qbar_output, Q0_output, kappae_output)
        
    double precision, intent(in)  :: rho, temp
    double precision, intent(out) :: alpha_output, qbar_output, q0_output, kappae_output
 
    ! Local variables
    double precision :: d_input, t_input, kappa_cak
 
    d_input = log10(rho)
    t_input = log10(temp)
 
    d_input = min(-7.d0-1.d-5, d_input)
    d_input = max(-16.d0+1.d-5, d_input)
    t_input = min(5d0-1.d-5, t_input)
    t_input = max(4d0+1.d-5, t_input)
 
    call get_val_comb(alpha_vals,qbar_vals,q0_vals,kappae_vals, &
                      logd_list, logt_list, d_input, t_input, &
                      alpha_output, qbar_output, q0_output, kappae_output)
 
  end subroutine set_cak_opacity
 
  !> This routine reads in 1-D (T,rho) values from a line opacity table and gives
  !> back as output the 1-D (T,rho) values and 2-D line opacity value
  subroutine read_table(D, T, K, filename)
 
    character(*), intent(in)      :: filename
    double precision, intent(out) :: K(iDmin:iDmax,iTmin:iTmax), D(iDmin:iDmax), T(iTmin:iTmax)
 
    ! Local variables
    character :: dum
    integer   :: row, col
 
    open(unit=1, status='old', file=trim(filename))
 
    ! Read temperature
    read(1,*) dum, t(itmin:itmax)
 
    ! Read rho and kappa
    do row = idmin,idmax
      read(1,*) d(row), k(row,itmin:itmax)
    enddo
 
    close(1)
 
  end subroutine read_table
 
  !> This subroutine looks in the table for the four couples (T,rho) surrounding
  !> a given input T and rho
  subroutine get_val_comb(K1_vals,K2_vals,K3_vals,K4_vals, &
                     logD_list, logT_list, D, T, &
                     K1, K2, K3, K4)
 
    double precision, intent(in)  :: K1_vals(iDmin:iDmax,iTmin:iTmax)
    double precision, intent(in)  :: K2_vals(iDmin:iDmax,iTmin:iTmax)
    double precision, intent(in)  :: K3_vals(iDmin:iDmax,iTmin:iTmax)
    double precision, intent(in)  :: K4_vals(iDmin:iDmax,iTmin:iTmax)
    double precision, intent(in)  :: logD_list(iDmin:iDmax), logT_list(iTmin:iTmax)
    double precision, intent(in)  :: D, T
    double precision, intent(out) :: K1, K2, K3, K4
 
    ! Local variables
    integer :: low_D_index, up_D_index, low_T_index, up_T_index
 
    if (d .gt. maxval(logd_list)) then
        ! print*, 'Extrapolating in logR'
        low_d_index = idmax-1
        up_d_index = idmax
    elseif (d .lt. minval(logd_list)) then
        ! print*, 'Extrapolating in logR'
        low_d_index = idmin
        up_d_index = idmin+1
    else
        call get_low_up_index(d, logd_list, idmin, idmax, low_d_index, up_d_index)
    endif
 
    if (t .gt. maxval(logt_list)) then
        ! print*, 'Extrapolating in logT'
        low_t_index = itmax-1
        up_t_index = itmax
    elseif (t .lt. minval(logt_list)) then
        ! print*, 'Extrapolating in logT'
        low_t_index = itmin
        up_t_index = itmin+1
    else
        call get_low_up_index(t, logt_list, itmin, itmax, low_t_index, up_t_index)
    endif
 
    call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
                       logd_list, logt_list, k1_vals, d, t, k1)
 
    call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
                       logd_list, logt_list, k2_vals, d, t, k2)
 
    call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
                       logd_list, logt_list, k3_vals, d, t, k3)
 
    call interpolate_krt(low_d_index, up_d_index, low_t_index, up_t_index, &
                       logd_list, logt_list, k4_vals, d, t, k4)
                       
  end subroutine get_val_comb
 
  !> This subroutine finds the indices in rho and T arrays of the two values
  !> surrounding the input rho and T
  subroutine get_low_up_index(var_in, var_list, imin, imax, low_i, up_i)
    
    integer, intent(in)          :: imin, imax
    double precision, intent(in) :: var_in, var_list(imin:imax)
    integer, intent(out)         :: low_i, up_i
 
    ! Local variables
    double precision :: low_val, up_val, res(imin:imax)
 
    res = var_list - var_in
    
    ! Find all bounding values for given input
    up_val = minval(res, mask = res .ge. 0) + var_in
    low_val = maxval(res, mask = res .le. 0) + var_in
 
    ! Find all bounding indices
    up_i = minloc(abs(var_list - up_val),1) + imin -1
    low_i = minloc(abs(var_list - low_val),1) + imin -1
 
    if (up_i .eq. low_i) low_i = low_i - 1
 
  end subroutine get_low_up_index
 
  !> This subroutine does a bilinear interpolation in the R,T-plane
  subroutine interpolate_krt(low_r, up_r, low_t, up_t, logD_list, logT_list, Kappa_vals, D, T, kappa_interp)
 
    integer, intent(in)           :: low_r, up_r, low_t, up_t
    double precision, intent(in)  :: Kappa_vals(iDmin:iDmax,iTmin:iTmax)
    double precision, intent(in)  :: logD_list(iDmin:iDmax), logT_list(iTmin:iTmax)
    double precision, intent(in)  :: D, T
    double precision, intent(out) :: kappa_interp
 
    ! Local variables
    double precision :: r1, r2, t1, t2, k1, k2, k3, k4, ka, kb
 
    ! First interpolate twice in the T coord to get ka and kb, then interpolate
    ! in the R coord to get ki
 
    !    r_1    rho       r_2
    !     |                |
    !     |                |
    ! ----k1--------ka-----k2----- t_1
    !     |          |     |
    !     |          |     |
    !   T |          |     |
    !     |          |     |
    !     |          ki    |
    !     |          |     |
    ! ----k3--------kb-----k4----- t_2
    !     |                |
    !     |                |
 
    r1 = logd_list(low_r)
    r2 = logd_list(up_r)
    t1 = logt_list(low_t)
    t2 = logt_list(up_t)
 
    k1 = kappa_vals(low_r, low_t)
    k2 = kappa_vals(low_r, up_t)
    k3 = kappa_vals(up_r, low_t)
    k4 = kappa_vals(up_r, up_t)
 
    call interpolate1d(r1,r2,d,k1,k2,ka)
    call interpolate1d(r1,r2,d,k3,k4,kb)
    call interpolate1d(t1,t2,t,ka,kb,kappa_interp)
 
  end subroutine interpolate_krt
 
  !> Interpolation in one dimension
  subroutine interpolate1d(x1, x2, x, y1, y2, y)
        
    double precision, intent(in)  :: x, x1, x2, y1, y2
    double precision, intent(out) :: y
 
    y = y1 + (x-x1)*(y2-y1)/(x2-x1)
 
  end subroutine interpolate1d
 
end module mod_cak_opacity