function grey_atmos,e,s=s,n=n,stefan=stefan,tsfc=tsfc

; This pro simulates a grey atmosphere with n-1 levels and the surface (n in total).
; All energy transport is by radiation - there is no convection.
;
; Each atmospheric level has absorption/emission coefficient e (in the "infra red").
; The surface has coefficient 1 (in IR and "short wave").
; S w/m2 reach and are absorbed by the surface. The atmosphere does not absorb in the SW.
;
; Inputs:
;  e   - emissivity of each level
;  s=s - "solar constant" 
;  n=n - number of levels. n=1 is permitted (unless tsfc is set) but dull.
;  stefan=stefan - stefans constant in R=e*stefan*T^4
;  tsfc- if set, this is imposed as the surface temperature
;        (this is the case of "no sfc feedback" - the sfc radiation is out of balance)
;
; Outputs
;  t   - (vector of) temperatures. t(0)=sfc t, t(1)=lowest layer,...,t(n-2)=top layer
;
; See-also: http://geosci.uchicago.edu/~archer/PS134/LabManual/lab.layer.htm
;
; Author: W M Connolley, wmc@bas.ac.uk

;
; Setup
;
if (n_elements(s) eq 0) then s=500.			; "Solar radiation" in W/M^2
if (n_elements(n) eq 0) then n=3			; N levels (incl. sfc)
if (n lt 1 or keyword_set(tsfc) and n lt 2) then message,'Unsuitable n levels'
if (n_elements(stefan) eq 0) then stefan=5.7e-8		; Not really variable in reality!
; Check the e is in [0,1]
if (e lt 0 or e gt 1) then message,'Emissivity must be between 0 and 1, not: '+string(e)

; 
; Make the interaction matrix, a.
;
; a#[u_sfc,u_1,u_2,....,u_n-2] represents the energy balance in this model,
; where [u] is a vector of the LW radiances from the surface up to the TOA.
;
; Consider a layer (i):
;
;    | Downwelling: sum-of[U_j*(transmitted-through-(|i-j|-1)-layers))]
;    | Of this, (Downwelling)*(emissivity) is absorbed in the layer
;    v                            ^ Emits: stefan*emissivity*T_j^4
; -------------- layer j ---------|---------------
;    ^                            V Emits: stefan*emissivity*T_j^4
;    |
;    | Upwelling: sum-of[U_j*(transmitted-through-(|i-j|-1)-layers))]
;      Of this, (Upwelling)*(emissivity) is absorbed in the layer
;
; Sanity checks.
;  
;   If n=1, we have  |      ^ U=stefan*T^4    (irrespective of e, which does not appear)
;                    V S    |
;                  ----------------- sfc, temperature T ---
;   S=U, so T=(S/stefan)^0.25
;   With S=500, T should be (500/5.7e-8)^0.25 = 306.037
;
;   If n=2, we have  |      ^ U1  ^ U0*(1-e) is radiated      U1=stefan*e*T1^4
;                    |      |     |
;                 ---|------------+ U0*e is absorbed -------- L1, temperature T1 ---
;                    |      |     ^
;                    V S    V U1  | U0=stefan*T^4
;                 ------------------------------------------- sfc, temperature T ---
;
;   so S+U1=U0 and 2*U1=U0*e, hence U0=S+U0*e/2 => U0=S/(1-e/2) => T=(S/(1-e/2)/stefan)^0.25
;   If e=0, T is the same as for n=1 (as we should expect)
;   If e=1, T = (S/(1/2)/stefan)^0.25 = 363.941 if S=500.
;
a=dblarr(n,n)
for i=0,n-1 do begin
  if (i gt 0) then begin                          ; Layers above the surface
    for j=0,i-1 do   a(i,j)=e*(1.-e)^(i-j-1)      ;  Upwelling from layers below
    for j=i+1,n-1 do a(i,j)=e*(1.-e)^(j-i-1)      ;  Downwelling from layers above
    a(i,i)=-2                                     ;  This layer
  endif else begin                                ; The surface
    for j=1,n-1 do   a(i,j)=(1.-e)^(j-1)          ;  Downwelling from above
    a(0,0)=-1                                     ;  This layer. Note -1 not -2 since only radiate in one direction.
  endelse
endfor

;
; Solve for the long-wave (IR) fluxes u
; Note that "invert" is a bad way to solve this equation but its simple.
;
if (n_elements(tsfc) eq 0) then begin

; Standard case - equilibrium everywhere
  rhs=[-s] & if (n gt 1) then rhs=[rhs,replicate(0,n-1)]
  u=invert(a)#rhs

endif else begin

; Nonstandard - fixed sfc temperature
  usfc=stefan*tsfc^4
  u=[usfc,invert(a(1:n-1,1:n-1))#(-usfc*a(1:n-1,0))]

endelse

;
; Having found the u's, calculate the t's via U=(stefan)*(emissivity)*T^4
;
ee=[1] & if (n gt 1) then ee=[ee,replicate(e,n-1)]
t=(u/(stefan*ee))^0.25

;
; Return the temperatures
;
return,t

end