% A simple radiative equilibrium energy balance model with an
% arbitrary number of layers (at least two).
%
% Ben Kravitz (ben.kravitz.work@gmail.com) 12 October 2011
clear all;
%% inputs
% L is the number of layers
% The model must have at least a surface and a TOA
% i.e., L must be at least 2
L=15;
% epsilon is an array of emissivities/absorptivities for each layer
% you can modify all of the layers
% individually or simultaneously as you see fit
% epsilon=1 is a completely opaque layer
% epsilon=0 is a completely transparent layer
epsilon=1.0*ones(1,L);
% Solar constant (in W/m2)
S0=1366;
% Planetary albedo (fraction)
alpha=0.3;
% Accuracy of the model
% If temperature changes are below this amount on successive
% iterations, then the model quits iterating. I've only tested
% this at 0.01, but it finishes in less than a second or two.
accuracy=0.01;
% Output temperature field
% If this is 1, it will output the equilibrium temperature at
% each layer when it is done iterating.
output=1;
%%%%%%%%%%%%%% model begins here %%%%%%%%%%%%%%
%%% modify anything below at your own peril %%%
%% constants
sigma=5.67e-8; % Stefan-Boltzmann constant
T=zeros(1,L); % initializing temperatures array
surf_sw=(S0/4)*(1-alpha);
if T(1)==0;
T(1)=(surf_sw/sigma)^(1/4); % initial surface temperature (shortwave only)
% the model is relatively insensitive to the choice of initial surface temperature
% so I just picked something reasonable
end
T2=zeros(1,L); % initializing comparison array
%% iterate over everything and stop when temperature change is small
while norm(T-T2)>accuracy;
T2=T;
rads=epsilon.*sigma.*T.^4;
%% surface (layer 1)
surf_sw=(S0/4)*(1-alpha);
terms1=zeros(1,L); % contributions of each layer to layer 1
terms1(2:L)=rads(2:L);
for n=2:L;
counter=n;
while counter>2;
terms1(n)=terms1(n)*(1-epsilon(counter-1));
counter=counter-1;
end
end
surf_lw=sum(terms1)*epsilon(1);
T(1)=((surf_sw+surf_lw)/sigma)^(1/4); % adjusting surface temperature to include longwave
%% layer n (2 to L-1) if these layers exist (i.e., if L>2)
if L>2;
for n=2:L-1;
termsn=zeros(1,L); % contributions of each layer to layer n
termsn(n+1:L)=rads(n+1:L); % layers above
for k=n+1:L;
counter=k;
while counter>n+1;
termsn(k)=termsn(k)*(1-epsilon(counter-1));
counter=counter-1;
end
end
termsn(1:n-1)=rads(1:n-1); % layers below
for k=1:n-1;
counter=k;
while counter