module coronal_impurities_m !! Module for reading and interpolating tabulated rate coefficients for !! radiation rates of impurities as a function of electron temperature !! Dependency on electron density is small and thus neglected use precision_grillix_m, only : GP, GP_NAN, GP_EPS use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST use variable_m, only: variable_t use screen_io_m, only : get_stdout use physical_constants_m, only : elementary_charge_SI use interpolation_m, only : interpol1d implicit none type, public :: coronal_impurity_t !! Class for radiation rate coefficients real(GP), allocatable, dimension(:) :: logte_tab !! Tabulated logarithm of the electron temperature in eV real(GP), allocatable, dimension(:) :: log_radiation_tab !! Tabulated logarithm of the normalised radiation rate coefficient integer, private :: tab_length !! Number of data points in the table real(GP), private :: te_ref !! Reference for electron temperature real(GP), private :: c_imp !! Impurity concentration (fraction of electron density) real(GP), private :: dlogte !! Spacing between logte entries contains procedure, public :: create => create_radrate_table procedure, public :: eval => eval_radrate end type contains subroutine create_radrate_table(self, file_name, tab_length, te_ref, & w_ref, c_imp) !! Create rate coefficient object from data file class(coronal_impurity_t), intent(inout) :: self !! Instance of class character(len=*), intent(in) :: file_name !! Name of atomic data file to read coefficients from integer, intent(in) :: tab_length !! Length of the data set real(GP), intent(in) :: te_ref !! Reference for electron temperature real(GP), intent(in) :: w_ref !! Reference for rate coefficient real(GP), intent(in) :: c_imp !! Impurity concentration (fraction of electron density) integer :: l, io_error character(len=256) :: io_errmsg real(GP) :: te, lz allocate(self%logte_tab(tab_length)) allocate(self%log_radiation_tab(tab_length)) open(unit = 40, file = file_name , status='old', action='read', & iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, & GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) read(40, *, iostat=io_error) ! Skip header read(40, *, iostat=io_error) ! Skip header read(40, *, iostat=io_error) ! Skip header read(40, *, iostat=io_error) ! Skip header read(40, *, iostat=io_error) ! Skip header if (io_error /= 0) call handle_error(io_errmsg, & GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) do l = 1, tab_length read(40, *, iostat=io_error) te, lz if (io_error /= 0) call handle_error(io_errmsg, & GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) self%logte_tab(l) = log(te) ! The input data is in [W*m3], w_ref is in [eV cm^3/s] self%log_radiation_tab(l) = log( lz * 1e6 / ( elementary_charge_SI * w_ref ) ) end do close(40) self%tab_length = tab_length self%te_ref = te_ref self%c_imp = c_imp self%dlogte = sum( self%logte_tab(2:tab_length) & - self%logte_tab(1:tab_length-1) ) / (tab_length-1) ! round to 4 digits (precision of logte table) self%dlogte = nint( self%dlogte * 1e4 ) * 1e-4 end subroutine real(GP) function eval_radrate(self, te, l) !! Evaluate the normalized rate coefficient through interpolation class(coronal_impurity_t), intent(in) :: self !! Instance of class real(GP), dimension(:), intent(in) :: te !! Electron temperature integer, intent(in) :: l !! Array index integer, parameter :: intorder=1 real(GP) :: logte_l, log_radiation logte_l = log(te(l) * self%te_ref) log_radiation = interpol1d(intorder, self%tab_length, self%dlogte, & self%log_radiation_tab, logte_l - self%logte_tab(1) ) eval_radrate = exp( log_radiation ) end function eval_radrate end module coronal_impurities_m