From 304ab17de528906a687ad162f8d58b8dea483c1b Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Mon, 9 Dec 2024 17:54:57 +1100 Subject: [PATCH 1/2] (readwrite_mesa) Add a few more options for column names --- src/setup/readwrite_mesa.f90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index d69a8ff72..e72c313e8 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -101,8 +101,9 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, call read_column_labels(iu,nheaderlines,ncols,nlabels,header) if (nlabels /= ncols) print*,' WARNING: different number of labels compared to columns' - allocate(m(lines),r(lines),pres(lines),rho(lines),ene(lines), & - temp(lines),Xfrac(lines),Yfrac(lines)) + allocate(m(lines)) + m = -1d0 + allocate(r,pres,rho,ene,temp,Xfrac,Yfrac,source=m) over_directions: do idir=1,2 ! try backwards, then forwards if (idir==1) then @@ -152,12 +153,18 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, r = (10**dat(1:lines,i)) * solarr case('pressure','p') pres = dat(1:lines,i) + case('logP') + pres = 10**dat(1:lines,i) case('temperature','t') temp = dat(1:lines,i) - case('x_mass_fraction_h','xfrac') + case('x_mass_fraction_h','x_mass_fraction_H','x','xfrac','h1') Xfrac = dat(1:lines,i) - case('y_mass_fraction_he','yfrac') + case('log_x') + Xfrac = 10**dat(1:lines,i) + case('y_mass_fraction_he','y_mass_fraction_He','y','yfrac','he4') Yfrac = dat(1:lines,i) + case('log_y') + Yfrac = 10**dat(1:lines,i) case default got_column = .false. end select @@ -176,6 +183,13 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, enddo over_directions close(iu) + if(min(minval(m),minval(r),minval(pres),minval(rho),minval(ene))<0d0)ierr = 1 + + if (ierr /= 0) then + print "(a,/)",' ERROR reading MESA file [missing required columns]' + return + endif + if (.not. usecgs) then m = m / umass r = r / udist From 9c5df94dad81a2ad17e4360327e156295f74da7d Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Mon, 9 Dec 2024 18:18:13 +1100 Subject: [PATCH 2/2] (readwrite_mesa) Add a few more options for column names --- src/setup/readwrite_mesa.f90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index e72c313e8..8e02843a5 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -138,25 +138,36 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, case('mass','#mass','m') m = dat(1:lines,i) if (ismesafile .or. maxval(m) < 1.e-10*solarm) m = m * solarm ! If reading MESA profile, 'mass' is in units of Msun + case('logM','log_mass') + m = 10**dat(1:lines,i) + if (ismesafile .or. maxval(m) < 1.e-10*solarm) m = m * solarm ! If reading MESA profile, 'mass' is in units of Msun case('rho','density') rho = dat(1:lines,i) - case('logrho') + case('logrho','logRho') rho = 10**(dat(1:lines,i)) case('energy','e_int','e_internal') ene = dat(1:lines,i) + case('logE') + ene = 10**dat(1:lines,i) case('radius_cm') r = dat(1:lines,i) + case('radius_km') + r = dat(1:lines,i) * 1e5 case('radius','r') r = dat(1:lines,i) if (ismesafile .or. maxval(r) < 1e-10*solarr) r = r * solarr - case('logr') + case('logr','logR') r = (10**dat(1:lines,i)) * solarr + case('logR_cm') + r = 10**dat(1:lines,i) case('pressure','p') pres = dat(1:lines,i) case('logP') pres = 10**dat(1:lines,i) case('temperature','t') temp = dat(1:lines,i) + case('logT') + temp = 10**dat(1:lines,i) case('x_mass_fraction_h','x_mass_fraction_H','x','xfrac','h1') Xfrac = dat(1:lines,i) case('log_x')