pro grey_atmoss,npe=npe

; Product of nlevels and emissivity
if (n_elements(e) eq 0) then npe=1.

; Do n=2,2^l
l=8

; 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 sfc and last temperature in
tsfc=fltarr(l)
tlast=fltarr(l)

; Frame for plotting in
!p.multi=[0,1,2]
plot,[200,350],[0,1.],/nodata, $
  xtitle='Temperature', $
  ytitle='Something like height or level', $
  title='Temperature through the atmosphere!C(nlevels)*(emiss-of-a-level) fixed at '+shtstr(npe)+'!CN levels=2^1...2^8', $
  ymargin=[5,5]

; Calculate ts's for each n, and add each to the plot
for m=0,l-1 do begin

; Number of levels...
  n=2^(m+1)
; Heights of the top of each level, scaled so that the TOA is at 1.
  lev_h=indgen(n)*2.^(l-m) & lev_h=lev_h/lev_h(n-1)

; Pass grey_atmos e so that e*(n-1)=npe. Note n-1 since there are only n-1
; atmospheric layers
  ts=grey_atmos(n=n,npe/(n-1.))

  plots,ts((indgen(n*2)+1)/2),lev_h(indgen(n*2)/2),col=cols(m),th=3
  plots,psym=2,ts(0),0,col=cols(m),th=3
  tsfc(m)=ts(0)
  tlast(m)=ts(n-1)

endfor

; Plot of sfc and last (top) temperatures
plot,2^(indgen(l)+1),tsfc, $
  yr=makerange([tsfc,tlast]), $
  xtitle='Number of levels', $
  ytitle='Temperature', $
  th=3,ymargin=[5,5], $
  title='Solid: Sfc T, Dashed: top layer T'
for i=0,m-1 do plots,psym=2,2^(i+1),tsfc(i),col=cols(i),th=3
oplot,2^(indgen(l)+1),tlast,th=3,li=2
for i=0,m-1 do plots,psym=4,2^(i+1),tlast(i),col=cols(i),th=3

end