mod_cak_opacity.t

Go to the documentation of this file.

module mod_cak_opacity
    implicit NONE
 
    !> min and max indices for R,T-range in opacity table
    integer, parameter :: dmin = 2
    integer, parameter :: dmax = 21
    integer, parameter :: tmin = 2
    integer, parameter :: tmax = 21
 
    !> The opacity tables are read once and stored globally in Kappa_vals
    double precision, public :: alpha_vals(dmin:dmax,tmin:tmax)
    double precision, public :: qbar_vals(dmin:dmax,tmin:tmax)
    double precision, public :: q0_vals(dmin:dmax,tmin:tmax)
    double precision, public :: kappa_e_vals(dmin:dmax,tmin:tmax)
 
    double precision, public :: log_d_list(dmin:dmax)
    double precision, public :: log_t_list(tmin:tmax)
 
    character(255), public :: amrvac_dir
    character(255), public :: fileplace
 
    public :: init_cak
    public :: set_cak_opacity
 
  contains
 
!> This routine is called when the fld radiation module is initialised.
!> Here, the tables for different He Abndcs are read and interpolated
subroutine init_cak(tablename)
  ! use mod_global_parameters
  character(6), intent(in) :: tablename
 
  CALL get_environment_variable("AMRVAC_DIR", amrvac_dir)
  fileplace = trim(amrvac_dir)//"/src/rhd/CAK_tables/"
 
  call read_table(log_d_list, log_t_list, alpha_vals, trim(tablename)//"/al_TD")
  call read_table(log_d_list, log_t_list, qbar_vals, trim(tablename)//"/Qb_TD")
  call read_table(log_d_list, log_t_list, q0_vals, trim(tablename)//"/Q0_TD")
  call read_table(log_d_list, log_t_list, kappa_e_vals, trim(tablename)//"/Ke_TD")
 
  ! print*, "Read Luka's tables"
 
end subroutine init_cak
 
!> This subroutine calculates the opacity for
!> a given temperature-density structure.
!> The opacities are read from a table that has the initialised metalicity
subroutine set_cak_opacity(rho,temp, gradv,alpha_output, Qbar_output, Q0_output, kappa_e_output)
  double precision, intent(in) :: rho, temp, gradv
  double precision, intent(out) :: alpha_output, qbar_output, q0_output, kappa_e_output
 
  double precision, PARAMETER :: const_c     = 2.99792458d10   ! cm S^-1           ; Speed of light
 
  double precision :: d_input, t_input
 
  double precision :: kappa_cak
  double precision :: tau, m_t
 
  d_input = dlog10(rho)
  t_input = dlog10(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,kappa_e_vals, &
                    log_d_list, log_t_list, d_input, t_input, &
                    alpha_output, qbar_output, q0_output, kappa_e_output)
 
  ! alpha = alpha_output
  ! Qbar = Qbar_output
  ! Q0 = Q0_output
  ! kappa_e = kappa_e_output
 
end subroutine set_cak_opacity
 
!> This routine reads out values and arguments from a table
subroutine read_table(D, T, K, filename)
    !> This routine reads in the the values for log kappa, and the values for log T and log R on the x and y axis
 
    double precision, intent(out) :: K(Dmin:Dmax,Tmin:Tmax), D(Dmin:Dmax), T(Tmin:Tmax)
    character(*), intent(in) :: filename
 
    character :: dum
    integer :: row, col
 
    OPEN(1,status = 'old', file=trim(fileplace)//filename)
 
    !> Read logT
    READ(1,*) dum,t(tmin:tmax)
 
    !> Read T and K
    do row = dmin,dmax !> NOT READING ENTIRE TABLE
      READ(1,*) d(row), k(row,tmin:tmax)
    enddo
 
    CLOSE(1)
 
end subroutine read_table
 
!>This subroutine looks in the table for the four couples (T,R)
!surrounding a given input for T and R
subroutine get_val_comb(K1_vals,K2_vals,K3_vals,K4_vals, &
                   Log_D_list, Log_T_list, D, T, &
                   K1, K2, K3, K4)
 
    double precision, intent(in) :: K1_vals(Dmin:Dmax,Tmin:Tmax)
    double precision, intent(in) :: K2_vals(Dmin:Dmax,Tmin:Tmax)
    double precision, intent(in) :: K3_vals(Dmin:Dmax,Tmin:Tmax)
    double precision, intent(in) :: K4_vals(Dmin:Dmax,Tmin:Tmax)
 
    double precision, intent(in) :: Log_D_list(Dmin:Dmax)
    double precision, intent(in) :: Log_T_list(Tmin:Tmax)
    double precision, intent(in) :: D, T
 
    double precision, intent(out) :: K1
    double precision, intent(out) :: K2
    double precision, intent(out) :: K3
    double precision, intent(out) :: K4
 
    integer :: low_r_index, up_r_index
    integer :: low_t_index, up_t_index
 
    if (d .gt. maxval(log_d_list)) then
        ! print*, 'Extrapolating in logR'
        low_r_index = dmax-1
        up_r_index = dmax
    elseif (d .lt. minval(log_d_list)) then
        ! print*, 'Extrapolating in logR'
        low_r_index = dmin
        up_r_index = dmin+1
    else
        call get_low_up_index(d, log_d_list, dmin, dmax, low_r_index, up_r_index)
    endif
 
    if (t .gt. maxval(log_t_list)) then
        ! print*, 'Extrapolating in logT'
        low_t_index = tmax-1
        up_t_index = tmax
    elseif ( t .lt. minval(log_t_list)) then
        ! print*, 'Extrapolating in logT'
        low_t_index = tmin
        up_t_index = tmin+1
    else
        call get_low_up_index(t, log_t_list, tmin, tmax, low_t_index, up_t_index)
    endif
 
    call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, &
                       log_d_list, log_t_list, k1_vals, d, t, k1)
 
    call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, &
                       log_d_list, log_t_list, k2_vals, d, t, k2)
 
    call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, &
                       log_d_list, log_t_list, k3_vals, d, t, k3)
 
    call interpolate_krt(low_r_index, up_r_index, low_t_index, up_t_index, &
                       log_d_list, log_t_list, k4_vals, d, t, k4)
end subroutine get_val_comb
 
 
!> this subroutine finds the indexes in R and T arrays of the two values surrounding the input R and T
subroutine get_low_up_index(x, x_list, imin, imax, low_i, up_i)
 
    integer, intent(in) :: imin, imax
    double precision, intent(in) :: x
    double precision, intent(in) :: x_list(imin:imax)
 
    integer, intent(out) :: low_i, up_i
 
    double precision :: low_val, up_val
    double precision :: res(imin:imax)
 
    res = x_list - x
 
    up_val = minval(res, mask = res .ge. 0) + x
    low_val = maxval(res, mask = res .le. 0) + x
 
    up_i = minloc(abs(x_list - up_val),1) + imin -1
    low_i = minloc(abs(x_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, Log_D_list, Log_T_list, Kappa_vals, D, T, k_interp)
 
    integer, intent(in) :: low_r, up_r, low_t, up_t
    double precision, intent(in) :: Kappa_vals(Dmin:Dmax,Tmin:Tmax)
    double precision, intent(in) :: Log_D_list(Dmin:Dmax)
    double precision, intent(in) :: Log_T_list(Tmin:Tmax)
    double precision, intent(in) :: D,T
    double precision, intent(out) :: k_interp
 
    double precision :: r1,r2,t1,t2
    double precision :: k1, k2, k3, k4
    double precision :: ka, kb
 
    !Cool ascii drawing of interpolation scheme: first interpolate twice in the T coord to get
    !ka and kb, then interpolate in the R coord to get ki
 
!   r_1    R        r_2
!     |                |
!     |                |
! ----k1--------ka-----k2----- t_1
!     |          |     |
!     |          |     |
!   T |          |     |
!     |          |     |
!     |          ki    |
!     |          |     |
! ----k3--------kb-----k4----- t_2
!     |                |
!     |                |
 
 
    r1 = log_d_list(low_r)
    r2 = log_d_list(up_r)
    t1 = log_t_list(low_t)
    t2 = log_t_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)
 
    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)
 
    ! print*, 'surounding indexes'
    ! print*, low_r, up_r, low_t, up_t
    ! print*, 'surounding input values'
    ! print*, r1,r2,t1,t2
    ! print*, 'surounding table values'
    ! print*, k1, k2, k3, k4
 
    call interpolate1d(r1,r2,d,k1,k2,ka)
    call interpolate1d(r1,r2,d,k3,k4,kb)
    call interpolate1d(t1,t2,t,ka,kb,k_interp)
 
    ! print*, 'interpolated value'
    ! print*, ka, kb, k_interp
 
end subroutine interpolate_krt
 
 
!> Interpolation in one dimension
subroutine interpolate1d(x1, x2, x, y1, y2, y)
 
    double precision, intent(in) :: x, x1, x2
    double precision, intent(in) :: y1, y2
    double precision, intent(out) :: y
 
    y = y1 + (x-x1)*(y2-y1)/(x2-x1)
 
end subroutine interpolate1d
 
end module mod_cak_opacity
 
! !> Interpolation on logarithmic scale
! subroutine log_interpolate1D(x1, x2, x, y1, y2, y)
!
!    double precision, intent(in) :: x, x1, x2
!    double precision, intent(in) :: y1, y2
!    double precision, intent(out) :: y
!
!    double precision :: expx, expx1, expx2
!    double precision :: expy1, expy2
!
!    expx = 10**x
!    expx1 = 10**x1
!    expx2 = 10**x2
!
!    expy1 = 10**y1
!    expy2 = 10**y2
!
!    y = expy1 + (expx-expx1)*(expy2-expy1)/(expx2-expx1)
!    y = log10(y)
!
! end subroutine log_interpolate1D