clear
close all
% Simple climate model
% Set parameter in the model_fun.m 
sec=365*24*3600;
step=1*sec;
options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]); %,'MaxStep',1);
T0=[0 0];
t0=1850*sec;
tk=2006*sec;
[t,Temp]=ode45(@model_fun,[t0,tk],T0,options);
[a1,lambda,data]=model_fun(1,[1,1]);
t=t/sec;
plot(t,Temp)
hold on
year=data.year;
FOR=data.TOT;
dF=interp1(year,FOR,t);
dTeq=dF*lambda;
plot(t,dTeq,'g')
ylabel('Temp Anomaly')
xlabel('Year')
legend('Mixed layer temp','Deep layer temp','Equailibrium temp')
grid


[data,x]=forcingdata;

figure
bar(x(:,1),x(:,2:11),'stoced')
ylabel('Radiative forcing [W/m^2]')
xlabel('year')
legend('GH','trop O3','Strat O3','Sulfate','Soot','Org Carbon','Biomass Burn',...
'Indirect','Solar','Volcanic')


figure
subplot(2,1,1)
plot(x(:,1),x(:,2),'r')
hold on
plot(x(:,1),x(:,5),'b')
legend('GH','Sulfate')
ylabel('Radiative forcing [W/m^2]')
subplot(2,1,2)
plot(x(:,1),x(:,3),'r')
hold on
plot(x(:,1),x(:,4),'b')
legend('Trop. O3','Strat. O3')
xlabel('Year')
ylabel('Radiative forcing [W/m^2]')


figure
subplot(2,1,1)
plot(x(:,1),x(:,6),'r')
hold on
plot(x(:,1),x(:,7),'b')
legend('Soot','Org. carbon')
ylabel('Radiative forcing [W/m^2]')
subplot(2,1,2)
plot(x(:,1),x(:,10),'r')
hold on
plot(x(:,1),x(:,9),'b')
legend('Solar','Indirect')
xlabel('Year')
ylabel('Radiative forcing [W/m^2]')

figure
subplot(2,1,1)
plot(x(:,1),x(:,8),'r')
hold on
plot(x(:,1),x(:,11),'b')
legend('Biomass burn.','Volcanic')
ylabel('Radiative forcing [W/m^2]')
subplot(2,1,2)
plot(x(:,1),sum(x(:,2:11)'),'k')
legend('Total')
xlabel('Year')
ylabel('Radiative forcing [W/m^2]')



return

figure
for k=1:5
   subplot(5,1,k)
   bar(x(:,1),x(:,k+1)) 
end    

figure
for k=6:10
   subplot(5,1,k-5)
   bar(x(:,1),x(:,k+1)) 
end    
