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
 
  use mod_comm_lib, only: mpistop
  use mod_global_parameters, only: mype
 
  implicit none
  private
 
  !> min and max indices for R,T-range in opacity table
  integer, parameter :: iRmin = 2, irmax = 20
  integer, parameter :: iTmin = 7, itmax = 76
 
  !> If user wants to read tables from different location than AMRVAC source
  !> NOTE: the tables should have same size as AMRVAC tables
  logical :: set_custom_tabledir = .false.
 
  !> 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(tabledir,set_custom_tabledir)
    
    character(len=*),  intent(in) :: tabledir
    logical, optional, intent(in) :: set_custom_tabledir
 
    ! Local variables
    character(len=256) :: amrvac_dir, path_table_dir
 
    if (present(set_custom_tabledir) .and. set_custom_tabledir) then
      path_table_dir = trim(tabledir)
    else
      call get_environment_variable("AMRVAC_DIR", amrvac_dir)
      path_table_dir = trim(amrvac_dir)//"/src/tables/OPAL_tables/"//trim(tabledir)
    endif
 
    ! Read in the OPAL table
    call read_table(logr_list, logt_list, kappa_vals, path_table_dir)
 
    if(mype==0)then
      write(*,*)' Initialized OPAL tables: read in the table from', path_table_dir
    endif
 
  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
 
    if (temp <= 0.0d0) call mpistop('Input temperature < 0 in set_opal_opacity.')
    if (rho <= 0.0d0) call mpistop('Input density < 0 in set_opal_opacity.')
 
    ! Convert input density and temperature to OPAL table convention
    logr_in = log10(rho/(temp*1.0d-6)**3.0d0)
    logt_in = log10(temp)
 
    !!print *,'enter get_kappa with logR_in and logT_in=',logR_in,logT_in
    call get_kappa(kappa_vals, logr_list, logt_list, logr_in, logt_in, logkappa_out)
    !!print *,'got logKappa_out=',logKappa_out
 
    ! If the outcome is 9.999 (no table entry), look right in the table
    !!do while (logKappa_out > 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 <= 1.0d-4)
    !!    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 = 10.0d0**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, funit=98
    logical   :: alive
 
    inquire(file=trim(filename), exist=alive)
 
    if (alive) then
      open(funit, file=trim(filename), status='unknown')
    else
      call mpistop('Table file you want to use cannot be found: '//filename)
    endif
 
    ! Skip first 4 lines
    do row = 1,4
        read(funit,*)
    enddo
 
    ! Read rho
    read(funit,*) dum, r(irmin:irmax)
    read(funit,*)
 
    ! Read T and kappa
    do row = itmin,itmax
        read(funit,'(f4.2,19f7.3)') t(row), k(row,irmin:irmax)
    enddo
 
    close(funit)
 
  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
 
    !print *,'enter with R and T=',R,T
    !print *,'check max-min(logR_list)=',maxval(logR_list),minval(logR_list)
    if (r .ge. maxval(logr_list)) then
      low_r_index = irmax-1
      up_r_index = irmax
    elseif (r .le. minval(logr_list)) then
      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 .ge. maxval(logt_list)) then
      low_t_index = itmax-1
      up_t_index = itmax
    elseif (t .le. minval(logt_list)) then
      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 >= 0) + var_in
    low_val = maxval(res, mask = res <= 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 == 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