pro rd_basel,fname,lambda,nmod,mnum,teff,logg,logz,flux
; This routine reads the BaSeL model file fname and returns the following values
;  lambda = wavelengths for flux samples (nm)
;  nmod = number of models in this file
;  mnum(nmod) = model number of each model
;  teff(nmod) = effective temp of each model (K)
;  logg(nmod) = log(g) for each model (cgs)
;  logz(nmod) = [Fe/H] for each model
;  flux(nlam,nmod) = surface flux in ergs/cm^2-s-Hz

; constants
c=2.997925e17                  ; speed of light, nm/s
h=6.626e-27                    ; Planck's constant, cgs

; open file, read wavelengths, make some derived quantities
nlam=1221
openr,iun,fname,/get_lun
lambda=fltarr(nlam)
f1='(8f10.1)'
readf,iun,lambda,format=f1
lam2=lambda^2
nu=c/lambda

; count lines to the end, reposition
point_lun,-iun,posn
ss=''
nl=0L
while(not eof(iun)) do begin
  readf,iun,ss
  nl=nl+1
endwhile
nmod=round(float(nl)/176.)
point_lun,iun,posn

; make output arrays
mnum=lonarr(nmod)
teff=fltarr(nmod)
logg=fltarr(nmod)
logz=fltarr(nmod)
flux=fltarr(nlam,nmod)
tflux=fltarr(nlam)

; read the data
f2='(i6,i6,2f6.2,2f6.2)'
f3='(7e11.4)'
for i=0,nmod-1 do begin
  readf,iun,mn,itef,lg,lm,vt,xh,format=f2
  mnum(i)=mn
  teff(i)=float(itef)
  logg(i)=lg
  logz(i)=lm
  readf,iun,tflux,format=f3
  flux(*,i)=tflux
endfor

; convert flux into desired units
for i=0,nmod-1 do begin
  ;flux(*,i)=flux(*,i)*0.4*c/lam2             ; ergs/cm^2-s-A
  ;flux(*,i)=flux(*,i)/(h*nu)
endfor

close,iun
free_lun,iun

end