%** Vostok ice core records of climate variables
%
%  V_CH4              457x2		methane: age (ka ago), gas concentration (ppmv)           
%  V_CO2              283x2  	carbon dioxide: age (ka ago),  gas conc. (ppmv)  
%  V_agedepth        3311x2  	age-depth model for this core; age (ka ago) depth (meters)
%  V_delduet_delT    3311x3  	delta-D and conversion to relative temperature: age (ka ago), ratio (per mil), relative temperature (C)
%  V_dustflux         500x2  	dust flux


load insolation_1Ma.mat
load PliPleHol_s_delO18.mat
load vostok_icecore_records.mat

%plot Vostok location on a map
load topo

%flat stereographic map
figure(30)
axesm('MapProjection','stereo','origin',[-90 0], 'FLatLimit',[-Inf 60]) ;
graticule=[180 360];
h=meshm(topo,topolegend,graticule);
demcmap(topo)
shading interp
plotm(-78.463, 106.826, 'r*')	%Vostok location
title('location of Vostok station in Antarctica (polar stereographic map)')



figure(31)
clf
subplot(3,1,1)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(1:423,1), PliPleHol_s_delO18(1:423,2), V_delduet_delT(:,1),V_delduet_delT(:,3),'plot');
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(get(AX(2),'Ylabel'),'String','\circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 425])
set(AX(2),'XLim',[0 425])
set(AX(1),'XTick',[0:25:425])
set(AX(1),'YLim',[2.4 4.5])
set(AX(1),'YTick',[2.5:0.5:4.5])
set(AX(2),'YLim',[-10 5])
set(AX(2),'YTick',[-10:5:5])
title('global benthic \delta O^{18} and Vostok temperature anomaly from \delta D')

subplot(3,1,2)
[AXa,H1a,H2a]=plotyy(V_CO2(:,1), V_CO2(:,2), V_delduet_delT(:,1), V_delduet_delT(:,3), 'plot');
set(get(AXa(1),'Ylabel'),'String','CO2 ppmv')
set(get(AXa(2),'Ylabel'),'String','\circ C')
set(AXa(1),'XGrid','on')
set(AXa(1),'XLim',[0 425])
set(AXa(2),'XLim',[0 425])
set(AXa(1),'XTick',[0:25:425])
set(AXa(1),'YLim',[100 300])
set(AXa(1),'YTick',[100:50:300])
set(AXa(2),'YLim',[-10 5])
set(AXa(2),'YTick',[-10:5:5])
title('Vostok CO2 and Vostok temperature anomaly from \delta D')

subplot(3,1,3)
[AX,H1,H2]=plotyy( V_CH4(:,1), V_CH4(:,2), V_delduet_delT(:,1), V_delduet_delT(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','CH4 ppmv')
set(get(AX(2),'Ylabel'),'String','\circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 425])
set(AX(2),'XLim',[0 425])
set(AX(1),'XTick',[0:25:425])
set(AX(1),'YLim',[100 800])
set(AX(1),'YTick',[200:200:800])
set(AX(2),'YLim',[-10 5])
set(AX(2),'YTick',[-10:5:5])
title('Vostok CH4 and Vostok temperature anomaly from \delta D')
xlabel('thousands of years ago')


figure(32)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(INSOL_Jun(1:426,1), INSOL_Jun(1:426,3), V_delduet_delT(:,1),V_delduet_delT(:,3),'plot');
set(get(AX(1),'Ylabel'),'String','W m^{-2}')
set(get(AX(2),'Ylabel'),'String','\circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 425])
set(AX(2),'XLim',[0 425])
set(AX(1),'XTick',[0:25:425])
set(AX(1),'YLim',[400 550])
set(AX(1),'YTick',[400:50:550])
set(AX(2),'YLim',[-10 5])
set(AX(2),'YTick',[-10:5:5])
title('June insolation at 60 N and Vostok temperature anomaly from \delta D')

subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Dec(1:426,1), INSOL_Dec(1:426,7), V_delduet_delT(:,1), V_delduet_delT(:,3), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^{-2}')
set(get(AXa(2),'Ylabel'),'String','\circ C')
set(AXa(1),'XGrid','on')
set(AXa(1),'XLim',[0 425])
set(AXa(2),'XLim',[0 425])
set(AXa(1),'XTick',[0:25:425])
set(AXa(1),'YLim',[400 550])
set(AXa(1),'YTick',[400:50:550])
set(AXa(2),'YLim',[-10 5])
set(AXa(2),'YTick',[-10:5:5])
title('December insolation at 60 S and Vostok temperature anomaly from \delta D')
xlabel('thousands of years ago')




figure(33)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(INSOL_Jun(1:426,1), INSOL_Jun(1:426,4), V_CH4(:,1), V_CH4(:,2),'plot');
set(get(AX(1),'Ylabel'),'String','W m^{-2}')
set(get(AX(2),'Ylabel'),'String','CH4 ppmv')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 425])
set(AX(2),'XLim',[0 425])
set(AX(1),'XTick',[0:25:425])
set(AX(1),'YLim',[425 575])
set(AX(1),'YTick',[400:50:600])
set(AX(2),'YLim',[300 800])
set(AX(2),'YTick',[200:200:800])
title('June insolation at 30 N and Vostok CH4')

subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Jun(1:426,1), INSOL_Jun(1:426,4), V_CO2(:,1), V_CO2(:,2), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^{-2}')
set(get(AXa(2),'Ylabel'),'String','CH4 ppmv')
set(AXa(1),'XGrid','on')
set(AXa(1),'XLim',[0 425])
set(AXa(2),'XLim',[0 425])
set(AXa(1),'XTick',[0:25:425])
set(AXa(1),'YLim',[425 575])
set(AXa(1),'YTick',[400:50:600])
set(AXa(2),'YLim',[100 300])
set(AXa(2),'YTick',[100:50:300])
title('June insolation at 30 N and Vostok CO2')
xlabel('thousands of years ago')


figure(34)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(1:418,1), PliPleHol_s_delO18(1:418,2), V_dustflux(:,1),V_dustflux(:,2),'plot');
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(get(AX(2),'Ylabel'),'String','dust flux mg m^{-2} a^-^1')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 425])
set(AX(2),'XLim',[0 425])
set(AX(1),'XTick',[0:25:425])
set(AX(1),'YLim',[2.4 4.5])
set(AX(1),'YTick',[2.5:0.5:4.5])
set(AX(2),'YLim',[0 15])
set(AX(2),'YTick',[0:5:15])
title('global benthic \delta O^{18} and Vostok dust flux')

subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Jun(1:426,1), INSOL_Jun(1:426,4), V_dustflux(:,1),V_dustflux(:,2),'plot');
set(get(AXa(1),'Ylabel'),'String','W m^{-2}')
set(get(AXa(2),'Ylabel'),'String','dust flux mg m^{-2} a^-^1')
set(AXa(1),'XGrid','on')
set(AXa(1),'XLim',[0 425])
set(AXa(2),'XLim',[0 425])
set(AXa(1),'XTick',[0:25:425])
set(AXa(1),'YLim',[425 575])
set(AXa(1),'YTick',[400:50:600])
set(AXa(2),'YLim',[0 15])
set(AXa(2),'YTick',[0:5:15])
title('June insolation at 30 N and Vostok dust flux')
xlabel('thousands of years ago')