pro grey_atmoss1,nl=nl,npe=npe,xr=xr,over=over

; N levels - default 100 (a random number)
if (n_elements(nl) eq 0) then n=100 else n=nl

; N*E's to use - can be a vector 
if (n_elements(npe) eq 0) then npe=[.1,1,10]
l=n_elements(npe)

; Colours
cols=[!p.color,1,2,3,10,40]
if (l gt 6) then cols=[cols,(wmc_guess_col(3*(l-6)))(2*(l-6)-1:3*(l-6)-2)]

; To keep temperatures in
tss=fltarr(l,n)

; Level heights
lev_h=indgen(n)/(n-1.)

; Calculate ts's for each npe
for m=0,l-1 do tss(m,*)=grey_atmos(n=n,npe(m)/n)

; Plot this
if (n_elements(xr) ne 2) then xr=makerange(tss)
if (not keyword_set(over)) then $
  plot,xr,[0,1.],/nodata, $
   xtitle='Temperature', $
   ytitle='Height', $
   title='Temperature through the atmosphere!C(nlevels)*(emiss-of-a-level) varying', $
   ymargin=[5,5]

for m=0,l-1 do begin
  plots,tss(m,(indgen(n*2)+1)/2),lev_h(indgen(n*2)/2),col=cols(m),th=3
  plots,psym=2,tss(m,0),0,col=cols(m),th=3
  xyouts,tss(m,0),0.5,align=0.5,'N*E='+shtstr(npe(m))
  plots,[1,1]*tss(m,0),[0.4,0.2]
endfor

end