% Plot monthly CO2 and HadCRUT4 from 1958 to ye
% Author: Vaughan R. Pratt, Stanford University
% September 21, 2020
% Must set ye before running this.  Can be any integer year in 1960..2020
%
% NB.  


figure(1); clf                              % The figure where it happened

% 1.  LOAD, CLEAN, AND PLOT CO2
websave('CO2.csv', 'http://clim8.stanford.edu/monthly_in_situ_co2_mlo.csv');
load CO2.csv
% Clean = Remove CO2 values that are undefined or beyond ye
defd = and((CO2(:, 4) < ye), (CO2(:, 5) > 0));
ddc = CO2(defd, 4);          % ddc - Domain of defined CO2
dCO2 = CO2(defd, 5);         % Defined CO2                                     % Defined CO2
yyaxis left
plot(ddc, dCO2); hold on
ylabel('Atmospheric CO2 (ppm)')

% 2.  LOAD, FILTER, AND PLOT HadCRUT4

options = weboptions;
options.Timeout = 30;
websave('HadCRUT4-gl.dat', ...
    'https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT4-gl.dat', ...
    options);
HC4mALL = importdata('HadCRUT4-gl.dat');
HC4m = HC4mALL(((1958:ye-1) - 1849)*2 - 1, 2:13)';
an = @(d)conv(d, ones(12, 1)/12, 'Valid');  % Annualize (12-mo moving avg)
% Annualizing is like annual mean but without sharp jump between Dec & Jan
yyaxis right
plot(an(1958+1/24:1/12:ye), an(HC4m(:)))            % 1/24 is mid-month
ylabel('Climate Anomaly (^oC)')

% 3. FINISH UP

set(gca, 'XLim', [1958 ye], 'XTick', 1800:5:2100)
title(sprintf('CO2 and Global Mean Surface Temperature 1958-%4d', ye))
legend('CO2 (Mauna Loa)', 'HadCRUT4', 'location', 'NorthWest')
saveas(1, sprintf('ClimCO2to%4d.jpg', ye))      % Save figure 1 as jpeg