clear
close all
% climate sensitivity	0,60 	K/(W/sqm)
% mixed-depth	100,0 	metres
% deep-depth	900,0 	metres
% diffusivity	0,00010 	m^2/sec
% time step	31536000,0 	secs
%	cs_mixed	421800000,0 	J/kg/sqm
%	cs_deep	3796200000,0 	J/kg/sqm
%	Diffusivity*density*specific ht	421,8 	J/K/m/sec
%	init mixed layer temp	0,0 	K
%	init deep layer temp	0,0 	K

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);
t=t/sec;
plot(t,Temp)
hold on
[data,x]=forcingdata;
year=data.year;
FOR=data.TOT;
dF=interp1(year,FOR,t);
[a,lambda]=model_fun(1900,[0,0]);
dTeq=dF*lambda;
plot(t,dTeq,'g')
ylabel('Temp Anomaly')
xlabel('Year')
legend('Mixed layer temp','Deep layer temp','Equailibrium temp')
grid

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)
bar(x(:,1),x(:,2),'r')
hold on
bar(x(:,1),x(:,5),'b')
legend('GH','Sulfate')
ylabel('Radiative forcing [W/m^2]')
subplot(2,1,2)
bar(x(:,1),x(:,3),'r')
hold on
bar(x(:,1),x(:,4),'b')
legend('Trop. O3','Strat. O3')
xlabel('Year')
ylabel('Radiative forcing [W/m^2]')
%subplotfix

figure
subplot(2,1,1)
bar(x(:,1),x(:,6),'r')
hold on
bar(x(:,1),x(:,5),'b')
legend('GH','Sulfate')
ylabel('Radiative forcing [W/m^2]')
subplot(2,1,2)
bar(x(:,1),x(:,3),'r')
hold on
bar(x(:,1),x(:,4),'b')
legend('Trop. O3','Strat. O3')
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