% zero-dimensional energy balance model 
% planet with an Earth-like atmosphere and temperature-dependent albedo

clear
Se=1370;            % solar constant W/m^2, Earth=1370

A=202.1;              % radiative heat loss coefficient
B=1.9;             % radiative heat loss coefficient

%* constants for temperature-dependent albedo parameterization
albedo_i=0.6;     % albedo of ice & snow surface
albedo_h=0.3;     % albedo of land surface

T_i=263;         % mean surface temperature for snowball earth
T_h=283;          % mean surface temperature for small-ice earth
Tnot=273.15;      % convert from K to degrees C

albedo=0.38;         % initial value for albedo
albedo_est=0;        % array to track estimates of albedo
albedo_tol=0.05;     % tolerance for precision in albedo estimate

tally=1;

% temperature-dependent albedo

while abs((albedo-albedo_est(tally))/albedo)>albedo_tol
    albedo_est(tally+1)=albedo;

    T=Tnot + (0.25*Se*(1-albedo)-A)/B;
    T_est(tally)=T;

    albedo= (T<T_i)*albedo_i ...
        + (T>T_i & T<T_h)*(albedo_i + (albedo_h-albedo_i)*((T-T_i)/(T_h-T_i))) ...
        + (T>T_h)*albedo_h;
    
    tally=tally+1;
end

albedo_est=[albedo_est(2:tally) albedo];
T_est=[T_est T];


figure(2)
clf
plot(albedo_est,T_est,'b.:')
hold on
plot(albedo, T, 'ro')
xlabel('albedo')
ylabel('temperature (K)')
legend('guesses', 'steady-state')
title('iteration to steady-state temperature with temperature-dependent albedo')