%**  create plots of sea surface temperature from sea-floor sediment cores
%    compare with computed insolation and deep-sea del-O18
 
%** data sets
%  insolation_1Ma.mat  
%  INSOL_Jun    1001 x 8 array: June insolation at 8 latitudes, 
%                               1 ka interval
%                               calculated by Berger, A., 1992
%               col1: age, ka,  origin (0) is 1950 A.D.
%               col2: insolation at 90 N W m^-2
%               col2: insolation at 60 N W m^-2
%               col2: insolation at 30 N W m^-2
%               col2: insolation at 0    W m^-2
%               col2: insolation at 30 S W m^-2
%               col2: insolation at 60 S W m^-2
%               col2: insolation at 90 S W m^-2
%
%  INSOL_Dec    1001 x 8 array: December insolation, as above
%
%
%  Berger, A., 1992, Orbital Variations and Insolation Database. IGBP 
%  PAGES/World Data Center-A for Paleoclimatology Data Contribution Series 
%  # 92-007.  NOAA/NGDC Paleoclimatology Program, Boulder CO, USA.
%
%  Berger A. and M.F. Loutre, 1991, Insolation values for the climate of 
%  the last 10 million of years. Quaternary Sciences Review, 10 (4), 
%  297-317.
% 
%
%  specmap_delO18_SST.mat       
%  RC11_SST     294 x 3   sea surface temperature: radiocarbon age ka, 
%                         warm SST deg C, cold SST deg C 
%  V27_SST      123 x 3   sea surface temperature: radiocarbon age ka, 
%                         warm SST deg C, cold SST deg C     
%  V30_SST      326 x 3   sea surface temperature: radiocarbon age ka, 
%                         warm SST deg C, cold SST deg C
%  SPEC_delO18  783 x 2   stacked benthic oxygen isotope record 
%                         radiocarbon age ka, normalized delO18 PDB per mil
%
%  Imbrie,J., A. McIntyre, and A. C. Mix, 1989, Oceanic Response to Orbital
%  Forcing in the Late Quaternary: Observational and Experimental Strategies, 
%  A. Berger, S. H. Schneider, and J.-C. Duplessy, (Eds.), Climate and 
%  Geosciences, A Challenge for Science and Society in the 21st Century,  
%  D. Reidel Publishing Company
%
%
%  PliPleHol_s_delO18.mat
%  PliPleHol_s_delO18 2116 x 3  stacked benthic oxygen isotopes and errors
%                               age ka, delta O18 PDB per mil, error per mil
%
%  Lisiecki, L. E., and M. E. Raymo (2004), A Plio-Pleistocene stack of 
%  57 globally distributed benthic d18O records, Paleoceanography, 
%  doi:10.1029/2004PA001071.
%

%clear
load insolation_1Ma.mat
load specmap_delO18_SST.mat
 
% plot core locations on a map
load topo
 
figure(20)
%clma                       %clear map axes if we've been fiddling already
axesm robinson
graticule=[180 360];
h=meshm(topo,topolegend,graticule);
demcmap(topo)
shading interp
plotm(52.836, -30.34, 'r+') %location of V27-007
textm(52.836, -30.34+1.5, 'V27')        %location of V27-007
plotm(-43.517, 79.867, 'r+')            %   "   of RC11-120
textm(-43.517, 79.867+1.5, 'RC')        %   "   of RC11-120
plotm(18.43, -21.08, 'r+')              %   "   of V30-49
textm(18.43, -21.08+1.5, 'V30')         %   "   of V30-49
title('three seafloor sediment cores')



figure(21)
clf
subplot(3,1,1)
plot(V27_SST(:,1), V27_SST(:,2), 'b-')
hold on
plot(V27_SST(:,1), V27_SST(:,3), 'g-')
ylabel('SST, \\circ C')
title('SST reconstructed from N. Atlantic core V27-116, b=summer, g=winter')
set(gca,'XGrid','on')
set(gca,'XTick',[0:25:250])
set(gca,'YTick',[0:5:20])
 
subplot(3,1,2)
[AX,H1,H2]=plotyy(INSOL_Jun(1:251,1), INSOL_Jun(1:251,3),V27_SST(:,1),V27_SST(:,2)-V27_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XTick',[0:25:250])
set(AX(1),'YLim',[400 600])
set(AX(1),'YTick',[400:50:600])
set(AX(2),'YLim',[2 10])
set(AX(2),'YTick',[2:2:10])
title('June insolation at 60 N and seasonality from V27-116')
 
subplot(3,1,3)
[AX,H1,H2]=plotyy(INSOL_Dec(1:251,1), INSOL_Dec(1:251,3),V27_SST(:,1),V27_SST(:,2)-V27_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XTick',[0:25:250])
set(AX(1),'YLim',[0 50])
set(AX(1),'YTick',[0:10:50])
set(AX(2),'YLim',[2 10])
set(AX(2),'YTick',[2:2:10])
title('December insolation at 60 N and seasonality from V27-116')
xlabel('thousands of years ago')
 
 
figure(22)
clf
subplot(3,1,1)
plot(RC11_SST(:,1), RC11_SST(:,2), 'b-')
hold on
plot(RC11_SST(:,1), RC11_SST(:,3), 'g-')
ylabel('SST, \circ C')
title('SST reconstructed from S. Indian Ocean core RC11-120, b=summer, g=winter')
set(gca,'XGrid','on')
set(gca,'XLim',[0 300])
set(gca,'XTick',[0:25:300])
set(gca,'YLim',[0 20])
set(gca,'YTick',[0:5:20])
 
subplot(3,1,2)
[AX,H1,H2]=plotyy(INSOL_Jun(1:301,1), INSOL_Jun(1:301,4),RC11_SST(:,1),RC11_SST(:,2)-RC11_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 300])
set(AX(2),'XLim',[0 300])
set(AX(1),'XTick',[0:25:300])
set(AX(1),'YLim',[400 600])
set(AX(1),'YTick',[400:50:600])
set(AX(2),'YLim',[0 8])
set(AX(2),'YTick',[0:2:8])
title('June insolation at 30 N and seasonality from RC11-120')
 
subplot(3,1,3)
[AX,H1,H2]=plotyy(INSOL_Dec(1:301,1), INSOL_Dec(1:301,4),RC11_SST(:,1),RC11_SST(:,2)-RC11_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 300])
set(AX(2),'XLim',[0 300])
set(AX(1),'XTick',[0:25:300])
set(AX(1),'YLim',[100 300])
set(AX(1),'YTick',[100:50:300])
set(AX(2),'YLim',[0 8])
set(AX(2),'YTick',[0:2:8])
title('December insolation at 30 N and seasonality from RC11-120')
xlabel('thousands of years ago')
 
 
figure(23)
clf
subplot(3,1,1)
plot(V30_SST(:,1), V30_SST(:,2), 'b-')
hold on
plot(V30_SST(:,1), V30_SST(:,3), 'g-')
ylabel('SST, \circ C')
title('SST reconstructed from Equatorial Atlantic core V30-49, b=summer, g=winter')
set(gca,'XGrid','on')
set(gca,'XLim',[0 325])
set(gca,'XTick',[0:25:325])
set(gca,'YLim',[10 30])
set(gca,'YTick',[10:5:30])
 
subplot(3,1,2)
[AX,H1,H2]=plotyy(INSOL_Jun(1:326,1), INSOL_Jun(1:326,3),V30_SST(:,1),V30_SST(:,2)-V30_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 325])
set(AX(2),'XLim',[0 325])
set(AX(1),'XTick',[0:25:325])
set(AX(1),'YLim',[400 600])
set(AX(1),'YTick',[400:50:600])
set(AX(2),'YLim',[2 12])
set(AX(2),'YTick',[2:2:12])
title('June insolation at 60 N and seasonality from site V30-49')
 
 
subplot(3,1,3)
[AX,H1,H2]=plotyy(INSOL_Jun(1:326,1), INSOL_Jun(1:326,5),V30_SST(:,1),V30_SST(:,2)-V30_SST(:,3), 'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','difference, \circ C')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 325])
set(AX(2),'XLim',[0 325])
set(AX(1),'XTick',[0:25:325])
set(AX(1),'YLim',[300 500])
set(AX(1),'YTick',[300:50:500])
set(AX(2),'YLim',[2 12])
set(AX(2),'YTick',[2:2:12])
title('June insolation at equator and seasonality from site V30-49')
xlabel('thousands of years ago')
 
 
figure(24)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(V30_SST(:,1),V30_SST(:,2)-V30_SST(:,3),SPEC_delO18(1:326,1), SPEC_delO18(1:326,2), 'plot');
set(get(AX(1),'Ylabel'),'String','difference, \circ C')
set(get(AX(2),'Ylabel'),'String','\delta O^{18} per mil')
hold on
plot(RC11_SST(:,1),RC11_SST(:,2)-RC11_SST(:,3),'r-')
plot(V27_SST(:,1),V27_SST(:,2)-V27_SST(:,3), 'g-')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 325])
set(AX(2),'XLim',[0 325])
set(AX(1),'XTick',[0:25:325])
set(AX(1),'YLim',[0 12])
set(AX(1),'YTick',[0:2:12])
set(AX(2),'YLim',[-2.5 2])
set(AX(2),'YTick',[-3:1:2])
title({'seasonality b=equatorial Atlantic, g=N. Atlantic, r=S. Indian ocean'; 'and global benthic \delta O^{18}'})
 
subplot(2,1,2)
[AX,H1,H2]=plotyy(V30_SST(:,1),V30_SST(:,2)-V30_SST(:,3),INSOL_Jun(1:326,1), INSOL_Jun(1:326,3), 'plot');
set(get(AX(1),'Ylabel'),'String','difference, \circ C')
set(get(AX(2),'Ylabel'),'String','W m^-^2')
hold on
plot(RC11_SST(:,1),RC11_SST(:,2)-RC11_SST(:,3),'r-')
plot(V27_SST(:,1),V27_SST(:,2)-V27_SST(:,3), 'g-')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 325])
set(AX(2),'XLim',[0 325])
set(AX(1),'XTick',[0:25:325])
set(AX(1),'YLim',[0 12])
set(AX(1),'YTick',[0:2:12])
set(AX(2),'YLim',[400 600])
set(AX(2),'YTick',[400:50:600])
title({'seasonality b=equatorial Atlantic, g=N. Atlantic, r=S. Indian ocean';'and June insolation 60 N'})
xlabel('thousands of years ago')