%% Choose some years, plot the fraction of ensemble members for which it was
%warmest or coldest in the 1400-2011 interval, plot some
%histograms, and plot the block average according to the process level. 

%% LOAD the arrays of max/min years, chooisng between the four time frames


selector=2011;


if selector==2011;
    
    load MaxYears_1400_2011 %MaxYears_1400_2011
    load MinYears_1400_2011 %MinYears_1400_2011
    timeline=1400:1:2011;

    Years_Array_Max=MaxYears_1400_2011; clear MaxYears_1400_2011
    Years_Array_Min=MinYears_1400_2011; clear MinYears_1400_2011
    
elseif selector==2003
    
    load MaxYears_1400_2003 %MaxYears_1400_2003
    load MinYears_1400_2003 %MinYears_1400_2003
    timeline=1400:1:2003;

    Years_Array_Max=MaxYears_1400_2003; clear MaxYears_1400_2003
    Years_Array_Min=MinYears_1400_2003; clear MinYears_1400_2003
    
elseif selector==1994
    
    load MaxYears_1400_1994 %MaxYears_1400_1994
    load MinYears_1400_1994 %MinYears_1400_1994
    timeline=1400:1:1994;

    Years_Array_Max=MaxYears_1400_1994; clear MaxYears_1400_1994
    Years_Array_Min=MinYears_1400_1994; clear MinYears_1400_1994
    
elseif selector==1959
    
    load MaxYears_1400_1959 %MaxYears_1400_1959
    load MinYears_1400_1959 %MinYears_1400_1959
    timeline=1400:1:1959;

    Years_Array_Max=MaxYears_1400_1959; clear MaxYears_1400_1959
    Years_Array_Min=MinYears_1400_1959; clear MinYears_1400_1959
    
    
end

%%

% we can also truncate off the last digit in these arrays. 
% gives the decade in which the warmest/coldest year occurred in each draw
% IN EACH CASE, THESE ARE XXX0-XXX9. (save for the last, which is
% variable)
Years_Array_Max_By_D=floor(Years_Array_Max/10)*10;
Years_Array_Min_By_D=floor(Years_Array_Min/10)*10;
timeline_By_D=unique(floor(timeline/10)*10);


%Load the input to BARCAST
load Inputs/BARCAST_INPUT_vHNL_Jan2012

Locs=BARCAST_INPUT.Master_Locs(BARCAST_INPUT.Inds_Central,:);
N_Locs=length(BARCAST_INPUT.Inds_Central);
N_Samps=length(Years_Array_Max(1,:));

%load the colormap
load HotCoolColormap_LONG


%% Decadal bins into a figure (cf Fig.2, thrid row):

fig=figure(606); clf
set(fig, 'position', [10,10,1250,400], 'renderer', 'painters');
plot([1,2]*1000, ones(1,2)*1000, 'r-', 'linewidth', 11), hold on
plot([1,2]*1000, ones(1,2)*1000, 'b-', 'linewidth', 7), hold on
legend('Max.', 'Min.', 'location', 'NorthWest') 

set(gca, 'fontsize', 15)
[cs, vs]=hist(Years_Array_Max_By_D(:), [timeline_By_D]);
ins=find(cs>0);
for kk=1:1:length(ins)
        %fill([vs(ins(kk))-4, vs(ins(kk))+4, vs(ins(kk))+4, vs(ins(kk))-4], [0, 0, cs(ins(kk))/length(Years_Array_Max_By_D(:)), cs(ins(kk))/length(Years_Array_Max_By_D(:))], 'r-', 'edgecolor', 'r'), hold on
        fill([vs(ins(kk))+1, vs(ins(kk))+9, vs(ins(kk))+9, vs(ins(kk))+1], [0, 0, cs(ins(kk))/length(Years_Array_Max_By_D(:)), cs(ins(kk))/length(Years_Array_Max_By_D(:))], 'r-', 'edgecolor', 'r'), hold on
end

[cs, vs]=hist(Years_Array_Min_By_D(:), [timeline_By_D]);
ins=find(cs>0);
for kk=1:1:length(ins)
        %fill([vs(ins(kk))-2, vs(ins(kk))+2, vs(ins(kk))+2, vs(ins(kk))-2], [0, 0, cs(ins(kk))/length(Years_Array_Min_By_D(:)), cs(ins(kk))/length(Years_Array_Min_By_D(:))], 'b-', 'edgecolor', 'b'), hold on
        fill([vs(ins(kk))+3, vs(ins(kk))+7, vs(ins(kk))+7, vs(ins(kk))+3], [0, 0, cs(ins(kk))/length(Years_Array_Min_By_D(:)), cs(ins(kk))/length(Years_Array_Min_By_D(:))], 'b-', 'edgecolor', 'b'), hold on
end
axis([1398 2020 0 0.25])
ylabel('Fraction')
xlabel('Year')
title('All data')

%% HISTOGRAMS FOR ICE-VARVE ENDING IN 1994

% also load the ICE_VARVE min/max (1400--1994)
load MaxYears_1400_1994_ICE_VARVE %MaxYears_1400_1994_ICE_VARVE
load MinYears_1400_1994_ICE_VARVE %MinYears_1400_1994_ICE_VARVE
timeline_IV=1400:1:1994;

Years_Array_Max_IV_By_D=floor(MaxYears_1400_1994_ICE_VARVE/10)*10;
Years_Array_Min_IV_By_D=floor(MinYears_1400_1994_ICE_VARVE/10)*10;
timeline_IV_By_D=unique(floor(timeline_IV/10)*10);

fig=figure(606); clf
set(fig, 'position', [10,10,1250,400], 'renderer', 'painters');
plot([1,2]*1000, ones(1,2)*1000, 'r-', 'linewidth', 11), hold on
plot([1,2]*1000, ones(1,2)*1000, 'b-', 'linewidth', 7), hold on
legend('Max.', 'Min.', 'location', 'NorthWest') 

set(gca, 'fontsize', 15)
[cs, vs]=hist(Years_Array_Max_IV_By_D(:), [timeline_IV_By_D]);
ins=find(cs>0);
for kk=1:1:length(ins)
        fill([vs(ins(kk))+1, vs(ins(kk))+9, vs(ins(kk))+9, vs(ins(kk))+1], [0, 0, cs(ins(kk))/length(Years_Array_Max_IV_By_D(:)), cs(ins(kk))/length(Years_Array_Max_IV_By_D(:))], 'r-', 'edgecolor', 'r'), hold on
end

[cs, vs]=hist(Years_Array_Min_IV_By_D(:), [timeline_IV_By_D]);
ins=find(cs>0);
for kk=1:1:length(ins)
        fill([vs(ins(kk))+3, vs(ins(kk))+7, vs(ins(kk))+7, vs(ins(kk))+3], [0, 0, cs(ins(kk))/length(Years_Array_Min_IV_By_D(:)), cs(ins(kk))/length(Years_Array_Min_IV_By_D(:))], 'b-', 'edgecolor', 'b'), hold on
end
axis([1398 2015 0 0.25])
ylabel('Fraction')
xlabel('Year')
title('Ice cores and varves, 1400-1994')

%% MAKE THE MAPS,2 ROWS FOR COLUMNS

% SELECT YEARS AND DECADES:
% ROW 1: Prob Hottest for two given years, prob hottest fell in two given
% decades
% ROW 2: As above, but for coldest. 
Timers_W=[2003, 2011, 1990, 2000];
Timers_C=[1453, 1816, 1690, 1880];

Fracs=NaN(N_Locs, 8);
for KK=1:1:N_Locs
    Fracs(KK,1)=length(find(Years_Array_Max(KK,:) == Timers_W(1)))/N_Samps;
    Fracs(KK,2)=length(find(Years_Array_Max(KK,:) == Timers_W(2)))/N_Samps;
    
    Fracs(KK,3)=length(find(Years_Array_Max_By_D(KK,:) == Timers_W(3)))/N_Samps;
    Fracs(KK,4)=length(find(Years_Array_Max_By_D(KK,:) == Timers_W(4)))/N_Samps;    
    
    Fracs(KK,5)=length(find(Years_Array_Min(KK,:) == Timers_C(1)))/N_Samps;
    Fracs(KK,6)=length(find(Years_Array_Min(KK,:) == Timers_C(2)))/N_Samps;
    
    Fracs(KK,7)=length(find(Years_Array_Min_By_D(KK,:) == Timers_C(3)))/N_Samps;
    Fracs(KK,8)=length(find(Years_Array_Min_By_D(KK,:) == Timers_C(4)))/N_Samps;    
end

% Plot 
Lat_pred_range=[45, 90];
Lon_pred_range=[-180, 180];

%find the locations that will be shown that will always be NaNs:
Pos_locs=[kron([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]', ones(size([Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]'))) , kron(ones(size([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]')), [Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]')];
Empty_locs=setdiff(Pos_locs, Locs, 'rows');

ColLims_T=[0,1];

figure(101), clf
fig=figure(101);
set(fig,'DoubleBuffer','on', 'color',  [1.0 1.0 1.0], 'position', [10,10,1350,950], 'renderer', 'painters');
    
ccmap=HotCoolColormap_LONG(16:1:end, :);
cN=length(ccmap);
colormap(ccmap)

for LL=1:1:4



    subplot(2, 4, LL)
    set(gca, 'fontSize', 15)
    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('patch',[.75 .75 .75],'edgecolor','k', 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

    
    
    %Find the NON NaN entries for the year:
    g_inds=find(isnan(Fracs(:, LL))==0);    
    for kk=1:1:length(g_inds);
            indi=max(1, ceil((Fracs(g_inds(kk), LL)-ColLims_T(1))/diff(ColLims_T)*cN));
            lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
            lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        if Fracs(g_inds(kk), LL)~=0
            m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
            m_plot([lon_mat, lon_mat], [lat_mat, lat_mat], 'k-', 'Color', ccmap(indi,:)), hold on
        elseif Fracs(g_inds(kk), LL)==0
            m_patch(fliplr(lon_mat), fliplr(lat_mat), [1,1,1], 'edgecolor', [1,1,1]), hold on
        end
        
    end

    %and the NaN entries
     for kk=1:1:length(Empty_locs(:,1));
         lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
         lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
         m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
     end


    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on
end

subplot(2,4,1)
ht1=text(-0.25,   1.1, ['Year'], 'fontsize', 18);
ht2=text(-0.25, -1.2, [num2str(Timers_W(1))], 'fontsize', 18);
tt=text(-1.25, -0.625,'Prob. Warmest', 'rotation', 90, 'fontsize', 18);

subplot(2,4,2)
%hh=title(['Decade: (',num2str(floor(Timers_W(2))-4), ', ',num2str(floor(Timers_W(2))+5), ')'], 'Position', [-0.15, +1.05, 1], 'fontsize', 18);
ht1=text(-0.25,   1.1, ['Year'], 'fontsize', 18);
ht2=text(-0.25, -1.2, [num2str(Timers_W(2))], 'fontsize', 18);

subplot(2,4,3)
ht1=text(-0.3,   1.1, ['Decade'], 'fontsize', 18);
ht2=text(-0.45, -1.2, [num2str(Timers_W(3)), ' to ',num2str(min(Timers_W(3)+9, selector))], 'fontsize', 18);

subplot(2,4,4)
ht1=text(-0.3, +1.1, ['Decade'], 'fontsize', 18);
ht2=text(-0.45, -1.2, [num2str(Timers_W(4)), ' to ',num2str(min(Timers_W(4)+9, selector))], 'fontsize', 18);

caxis('manual')
caxis(ColLims_T)
pp(LL)=colorbar('vert', 'Position', [0.9326 0.5976 0.0215 0.3392], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');    
cbfreeze  


ccmap=flipud(HotCoolColormap_LONG(1:1:13, :));    
cN=length(ccmap);
colormap(ccmap)

for LL=5:1:8

    %Max
    subplot(2, 4, LL)
    set(gca, 'fontSize', 15)
    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('patch',[.75 .75 .75],'edgecolor','k', 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on    
    
    %Find the NON NaN entries for the year:
    g_inds=find(isnan(Fracs(:, LL))==0);    
    for kk=1:1:length(g_inds);
            indi=max(1, ceil((Fracs(g_inds(kk), LL)-ColLims_T(1))/diff(ColLims_T)*cN));
            lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
            lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        if Fracs(g_inds(kk), LL)~=0
            m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
            m_plot([lon_mat, lon_mat], [lat_mat, lat_mat], 'k-', 'Color', ccmap(indi,:)), hold on
        %white for zero
        elseif Fracs(g_inds(kk), LL)==0
            m_patch(fliplr(lon_mat), fliplr(lat_mat), [1,1,1], 'edgecolor', [1,1,1]), hold on
        end
        
    end

    %and the NaN entries
     for kk=1:1:length(Empty_locs(:,1));
         lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
         lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
         m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
     end

    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

end
subplot(2,4,5)
ht2=text(-0.25, -1.2, [num2str(Timers_C(1))], 'fontsize', 18);
text(-1.25, -0.625,'Prob. Coldest', 'rotation', 90, 'fontsize', 18);
subplot(2,4,6)
ht2=text(-0.25, -1.2, [num2str(Timers_C(2))], 'fontsize', 18);
subplot(2,4,7)
ht2=text(-0.45, -1.2, [num2str(Timers_C(3)), ' to ',num2str(min(Timers_C(3)+9, selector))], 'fontsize', 18);
subplot(2,4,8)
ht2=text(-0.45, -1.2, [num2str(Timers_C(4)), ' to ',num2str(min(Timers_C(4)+9, selector))], 'fontsize', 18);

caxis('manual')
caxis(ColLims_T)
pp(LL)=colorbar('vert', 'Position', [0.9326 0.13 0.0215 0.3392], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');    
cbfreeze  




