|
From: Corentin C. <cor...@en...> - 2016-08-10 20:42:43
|
Hi,
I am trying to add some extra columns to the output of MESA (kappa_T and
kappa_mu). I wrote the following piece of code:
subroutine data_for_extra_profile_columns(id, id_extra, n, nz, names,
vals, ierr)
use star_def, only: star_info, maxlen_profile_column_name
use const_def, only: dp
integer, intent(in) :: id, id_extra, n, nz
character (len=maxlen_profile_column_name) :: names(n)
real(dp) :: vals(nz,n)
real(dp), allocatable :: logLambdaHHe(:)
integer, intent(out) :: ierr
type (star_info), pointer :: s
integer :: k
ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
!note: do NOT add the extra names to profile_columns.list
! the profile_columns.list is only for the built-in profile column
options.
! it must not include the new column names you are adding here.
! here is an example for adding a profile column
!if (n /= 1) stop 'data_for_extra_profile_columns'
!names(1) = 'beta'
!do k = 1, nz
! vals(k,1) = s% Pgas(k)/s% P(k)
!end do
allocate(logLambdaHHe(s%nz))
logLambdaHHe = -19.26d0 - 0.5d0*log(s%rho) + 1.5d0*log(s%T) &
& - 0.5d0*log((s%X + 3d0)/2d0) - log(2d0)
names(1) = 'kappa_T'
vals(:, 1) = (16*boltz_sigma*s%T**3) / (3*s%opacity*s%rho**2*s%cp)
names(2) = 'kappa_mu'
vals(:, 2) = (15d0)/(16d0 * s%rho * logLambdaHHe) &
&* sqrt(2d0*mp/(5*pi)) * (boltzm*s%T)**(5d0/2d0) / qe**4 &
&* (3d0+s%X) / ((1d0+s%X)*(3d0+5d0*s%X)*(0.7d0 + 0.3d0*s%X))
end subroutine data_for_extra_profile_columns
However, if I run it, the code crashes because the size of s%rho is not
s%nz (nb, s%nz == nz). Is there a reason for that?
Thank you,
Corentin
|