%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=energy_balance_0d;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a zero dimensional global energy balance model including ice albedo feedback
% for APM 115
% Eli tziperman.
close all; clear all;

% integrate the time dependent equation:
figure
S0=1366/4; % watts/m^2
T_init=-10;
[time,T]=integrate_energy_balance_eqn(T_init,S0);
hold on
pause
T_init=-20;
[time,T]=integrate_energy_balance_eqn(T_init,S0);
pause

% plot albedo, terms in the energy balance, and RHS of energy balance equation:
S0=1366/4; % watts/m^2
plot_energy_balance_terms_and_rhs(S0)
pause

S0=1000/4; % watts/m^2
plot_energy_balance_terms_and_rhs(S0)
pause

S0=5000/4; % watts/m^2
plot_energy_balance_terms_and_rhs(S0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_energy_balance_terms_and_rhs(S0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ScreenSize is a four-element vector: [left, bottom, width, height]:
scrsz = get(0,'ScreenSize'); 
figure('Position',[scrsz(4)*0.1 scrsz(4)*0.35 scrsz(3)*0.8 scrsz(4)*0.5])

T_min=-150; T_max=200; dT=1;

i=0;
for T=T_min:dT:T_max
  i=i+1;
  alpha(i)=albedo(T);
  temperature(i)=T;
end
subplot(1,3,1)
plot(temperature,alpha);
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min,T_max,0, 1]);
h1=ylabel('albedo');
h2=xlabel('temperature');
h3=title(sprintf('S0=%gW/m^2',4*S0));
set([h1,h2,h3,gca],'fontsize',14)

% plot incoming*(1-albedo) vs outgoing radiation as function of temperature:
i=0;
for T=T_min:dT:T_max
  i=i+1;
  incoming(i)=S0*(1-albedo(T));
  outgoing(i)=longwave(T);
  temperature(i)=T;
end
subplot(1,3,2)
plot(temperature,incoming,'r');
hold on
plot(temperature,outgoing,'b');
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min,T_max,V(3),V(4)]);
h1=ylabel('radiation (watts/m^2)');
h2=xlabel('temperature');
h3=title('red: incoming; blue: outgoing');
set([h1,h2,h3,gca],'fontsize',14)


% plot incoming*(1-albedo)-outgoing radiation as function of temperature:
subplot(1,3,3)
plot(temperature,incoming-outgoing,'r');
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min,T_max,V(3),V(4)]);
hold on
plot(temperature,zeros(length(temperature),1),'k')
h1=ylabel('radiation (watts/m^2)');
h2=xlabel('temperature');
h3=title('dT/dt=incoming-outgoing');
set([h1,h2,h3,gca],'fontsize',14)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a=albedo(T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T in degree celsius
a_min=0.1;
a_max=0.9;
T_min=-40;
T_max=20;

if T>T_max 
  a=a_min;
elseif T<T_min
  a=a_max;
else
  a=a_max-(a_max-a_min)*(T-T_min)/(T_max-T_min);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function LW=longwave(T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% outgoing long wave radiation
% T in celsius
emissivity=0.6;
sigma=5.67e-8; % J/(sec m^2 K^4)

LW=emissivity*sigma*(T+273)^4;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [time,T]=integrate_energy_balance_eqn(T_init,S0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tspan=[0:0.011:10];
T0=T_init;  % -10 or -20 to get the two equilibria
[time,T] = ode45(@(t,X) myrhs(t,X,S0),tspan,T0);
plot(time,T)
h1=title('zero dimensional energy balance model');
h2=xlabel('time');
h3=ylabel('T');
set([h1,h2,h3,gca],'fontsize',14)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs=myrhs(t,T,S0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
incoming=S0*(1-albedo(T));
outgoing=longwave(T);
rhs=incoming-outgoing;

