mod_opal_opacity.t

Go to the documentation of this file.

!> This module reads in Rosseland-mean opacities from OPAL tables. Table opacity
!> values are given in base 10 logarithm and are a function of mass density (R)
!> and temperature (T), which are also both given in base 10 logarithm.
 
module mod_opal_opacity
  
  implicit none
 
  !> min and max indices for R,T-range in opacity table
  integer, parameter :: irmin = 2, irmax = 20
  integer, parameter :: itmin = 7, itmax = 76
 
  !> The opacity tables are read once and stored globally in Kappa_vals
  double precision, public :: kappa_vals(itmin:itmax,irmin:irmax)
  double precision, public :: logr_list(irmin:irmax), logt_list(itmin:itmax)
 
  public :: init_opal_table
  public :: set_opal_opacity
 
contains
 
  !> This subroutine is called when the FLD radiation module is initialised.
  !> The OPAL tables for different helium abundances are read and interpolated.
  subroutine init_opal_table(tablename)
    
    character(len=*), intent(in) :: tablename
 
    ! Local variables
    character(len=:), allocatable :: amrvac_dir, path_table_dir
    
    call get_environment_variable("AMRVAC_DIR", amrvac_dir)
        
    ! Absolute path to OPAL table location
    path_table_dir = trim(amrvac_dir)//"/src/rhd/OPAL_tables/"//trim(tablename)
    
    ! Read in the OPAL table
    call read_table(logr_list, logt_list, kappa_vals, path_table_dir)
 
  end subroutine init_opal_table
 
  !> This subroutine calculates the opacity for a given temperature-density
  !> structure. Opacities are read from a table with given metalicity.
  subroutine set_opal_opacity(rho, temp, kappa)
        
    double precision, intent(in)  :: rho, temp
    double precision, intent(out) :: kappa
 
    ! Local variables
    double precision :: logr_in, logt_in, logkappa_out
 
    ! Convert input density and temperature to OPAL table convention
    logr_in = log10(rho/(temp*1d-6)**3)
    logt_in = log10(temp)
 
    call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
 
    ! If the outcome is 9.999 (no table entry), look right in the table
    do while (logkappa_out .gt. 9.0d0)
        logr_in = logr_in + 0.5d0
        call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
    enddo
 
    ! If the outcome is NaN, look left in the table
    do while (logkappa_out .eq. 0.0d0)
        logr_in = logr_in - 0.5d0
        call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
    enddo
 
    ! Convert OPAL opacity for output
    kappa = 10d0**logkappa_out
 
  end subroutine set_opal_opacity
 
  !> This routine reads in 1-D (rho,T) values from an opacity table and gives
  !> back as output the 1-D (rho,T) values and 2-D opacity
  subroutine read_table(R, T, K, filename)
        
    character(*), intent(in)      :: filename
    double precision, intent(out) :: K(iTmin:iTmax,iRmin:iRmax), R(iRmin:iRmax), T(iTmin:iTmax)
 
    ! Local variables
    character :: dum
    integer   :: row, col
 
    open(unit=1, status='old', file=trim(filename))
 
    ! Skip first 4 lines
    do row = 1,4
        read(1,*)
    enddo
 
    ! Read rho
    read(1,*) dum, r(irmin:irmax)
    read(1,*)
 
    ! Read T and kappa
    do row = itmin,itmax
        read(1,'(f4.2,19f7.3)') t(row), k(row,irmin:irmax)
    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_kappa(Kappa_vals, logR_list, logT_list, R, T, K)
    
    double precision, intent(in)  :: Kappa_vals(iTmin:iTmax,iRmin:iRmax)
    double precision, intent(in)  :: logR_list(iRmin:iRmax), logT_list(iTmin:iTmax)
    double precision, intent(in)  :: R, T
    double precision, intent(out) :: K
 
    ! Local variables
    integer :: low_R_index, up_R_index
    integer :: low_T_index, up_T_index
 
    if (r .gt. maxval(logr_list)) then
        ! print*, 'Extrapolating in logR'
        low_r_index = irmax-1
        up_r_index = irmax
    elseif (r .lt. minval(logr_list)) then
        ! print*, 'Extrapolating in logR'
        low_r_index = irmin
        up_r_index = irmin+1
    else
        call get_low_up_index(r, logr_list, irmin, irmax, low_r_index, up_r_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_r_index, up_r_index, low_t_index, up_t_index, logr_list, logt_list, kappa_vals, r, t, k)
 
  end subroutine get_kappa
 
  !> 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, logR_list, logT_list, Kappa_vals, R, T, kappa_interp)
    
    integer, intent(in)           :: low_r, up_r, low_t, up_t
    double precision, intent(in)  :: Kappa_vals(iTmin:iTmax,iRmin:iRmax)
    double precision, intent(in)  :: logR_list(iRmin:iRmax), logT_list(iTmin:iTmax)
    double precision, intent(in)  :: R, 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 = logr_list(low_r)
    r2 = logr_list(up_r)
    t1 = logt_list(low_t)
    t2 = logt_list(up_t)
 
    k1 = kappa_vals(low_t, low_r)
    k2 = kappa_vals(low_t, up_r)
    k3 = kappa_vals(up_t, low_r)
    k4 = kappa_vals(up_t, up_r)
 
    call interpolate1d(r1,r2,r,k1,k2,ka)
    call interpolate1d(r1,r2,r,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_opal_opacity