program readMARCS * Sample programme to call subroutine readmodel to store data in arrays * Version 1.1, 21-Dec-2004 /BE * Version 1.2, 3-Apr-2012 /BE * * Parameters are described in subroutine readmodel * * Compile by e.g.: "f90 -o readmarcs readmarcs.f" * * Run by either "cat model_file | readmarcs" if unzipped * or perhaps "gzcat model_file | readmarcs" if gzipped * implicit none integer ndp,nmolmax,k parameter(ndp=100) parameter (nmolmax=30) character fileID*72 *** real teff,flux,g,xite,rmass,rmetal,alpha,radius,rluminosity real palfa,pny,py,pbeta,abund(92),xx,yy,zz common /global/ teff,flux,g,xite,rmass,rmetal,alpha,radius, & rluminosity,palfa,pny,py,pbeta,abund,xx,yy,zz *** integer ndepth real tau(ndp),taus(ndp),z(ndp),t(ndp),pe(ndp),pg(ndp),prad(ndp) real pturb(ndp),xkapr(ndp),ro(ndp),emu(ndp),v(ndp),frconv(ndp) real rhox(ndp) common /structure/ ndepth,tau,taus,z,t,pe,pg,prad,pturb,xkapr, & ro,emu,v,frconv,rhox *** real hnic(ndp),other(ndp),presmo(nmolmax,ndp) common /pressures/ hnic,presmo,other *** call readmodel(fileID) * *** List part of the data: * print *,'ID: ',fileID print *,'Teff=',teff,', log g=',alog10(g),', [Fe/H]=',rmetal print *,' X, Y and Z=',xx,yy,zz write(*,1000) 1000 format(' k log(TauRoss) T column mass log P(CO)') do k=1,ndepth write(*,1010) k,tau(k),t(k),rhox(k),presmo(7,k) 1010 format(i3,f13.3,f9.2,f15.5,f10.3) enddo END * * *...v....1....v....2....v....3....v....4....v....5....v....6....v....7.. * * subroutine readmodel(fileID) * reads a MARCS model and stores data into arrays * in 3 common blocks: 'global', 'structure' and 'pressures' * ********************************************************************** * fileID returns the model ID string (a72) ********************************************************************** * Global properties of the model: common "global" * teff = effective temperature [K] * flux = physical flux [erg/cm2/s] * g = surface acceleration of gravity at tau(Rosseland)=1.0 [cm/s2] * xite = microturbulence parameter [km/s] * rmass = stellar mass [Msun] (0.0 for plane-parallel models) * rmetal = [Fe/H] * alpha = [alpha/Fe] * radius = Stellar radius [cm] at Tau(Rosseland) = 1.0 * (1.0 cm for plane-parallel models) * rluminosity = luminosity [Lsun] for the radius given above * palfa pny py pbeta = convection parameters * Ref: Henyey L., Vardaya M.S., Bodenheimer P. 1965, ApJ 142, 841 * abund() = logarithmic number abundances of the 92 first chemical * elements in the periodic table * XX = Hydrogen mass fraction * YY = Helium mass fraction * ZZ = Metal mass fraction ************** * Model structure and thermodynamics, depth index = k: common "structure" * ndepth = number of depth points * tau(k) = log(tau(Rosseland)) * taus(k) = log(tau(5000 A)) * z(k) = depth [cm], depth=0.0 @ tau(Rosseland)=1.0 * t(k) = temperature [K] * pe(k) = electron pressure [dyn/cm2] * pg(k) = gas pressure [dyn/cm2] * prad(k) = radiation pressure [dyn/cm2] * pturb(k) = turbulence pressure [dyn/cm2] * xkapr(k) = Rosseland opacity [cm2/g] * ro(k) = density [g/cm3] * emu(k) = mean molecular weight [amu] * v(k) = convective velocity [cm/s] * frconv(k) = Fractional convective flux Fconv/F * rhox(k) = Column mass above point k [g/cm2] ************** * Partial pressures: common "pressures" * hnic(k) = log(H I pressure/[1 dyn/cm2]) * presmo(j,k) = log(XXX pressure/[1 dyn/cm2]) for j:XXX * 1:H-, 2:H2, 3:H2+, 4:H2O, 5:OH, 6:CH, 7:CO, * 8:CN, 9:C2, 10:N2, 11:O2, 12:NO, 13:NH, * 14:TiO, 15:C2H2, 16:HCN, 17:C2H, 18:HS, * 19:SiH, 20:C3H, 21:C3, 22:CS, 23:SiC, * 24:SiC2, 25:NS, 26:SiN, 27:SiO, 28:SO, 29:S2, 30:SiS * other(k) = log(pressure of not listed gas constituents/[dyn/cm2]) ********************************************************************** implicit none integer ndp,nmolmax parameter(ndp=100) parameter (nmolmax=30) character fileID*72 real taudum,pgas integer k,kk,j *** real teff,flux,g,xite,rmass,rmetal,alpha,radius,rluminosity real palfa,pny,py,pbeta,abund(92),xx,yy,zz common /global/ teff,flux,g,xite,rmass,rmetal,alpha,radius, & rluminosity,palfa,pny,py,pbeta,abund,xx,yy,zz *** integer ndepth real tau(ndp),taus(ndp),z(ndp),t(ndp),pe(ndp),pg(ndp),prad(ndp) real pturb(ndp),xkapr(ndp),ro(ndp),emu(ndp),v(ndp),frconv(ndp) real rhox(ndp) common /structure/ ndepth,tau,taus,z,t,pe,pg,prad,pturb,xkapr, & ro,emu,v,frconv,rhox *** real hnic(ndp),other(ndp),presmo(nmolmax,ndp) common /pressures/ hnic,presmo,other *** read(*,'(a)') fileID read(*,1170) TEFF,FLUX,G,xite,rmass,rmetal,alpha,radius, & rluminosity 1170 format(f7.0,/,e12.4,/,e12.4,/,f5.1,/,f5.1,/,2f6.2,/,e12.4,/,e12.4) read(*,1176) PALFA,PNY,PY,PBETA 1176 format(f6.2,f5.2,f6.3,f5.2) read(*,*) xx,yy,zz 1178 format(x,3f8.5) read(*,*) read(*,1180) abund 1180 format(10f7.2) * * structure * read(*,2105) ndepth 2105 format(i4) read(*,*) read(*,*) do kk=1,ndepth read(*,*) k,tau(k),taus(k),z(k),t(k), & pe(k),pg(k),prad(k),pturb(k) if(kk.ne.k) stop 'model and loop indices differ' enddo * * thermodynamics * read(*,*) do kk=1,ndepth read(*,*) k,taudum,xkapr(k),ro(k),emu(k),v(k),frconv(k), & rhox(k) if(kk.ne.k) stop 'model and loop indices differ' enddo * * molecular partial pressures * read(*,*) read(*,*) do kk=1,ndepth read(*,*) k,pgas,hnic(k),(presmo(j,k),j=1,9) if(kk.ne.k) stop 'model and loop indices differ' enddo read(*,*) do kk=1,ndepth read(*,*) k,(presmo(j,k),j=10,20) if(kk.ne.k) stop 'model and loop indices differ' enddo read(*,*) do kk=1,ndepth read(*,*) k,(presmo(j,k),j=21,30),other(k) if(kk.ne.k) stop 'model and loop indices differ' enddo * return END * * *...v....1....v....2....v....3....v....4....v....5....v....6....v....7.. * *