%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 1 - Time Series %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTANT!! -- MUST run Figure 2 first to create Ensemble data
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs

for startup =1
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/SOLACE'))
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))

% Model Info and Colormaps
mods = dir;
model_names = {'iHAMMOCC','CMOC','WOMBAT','OECOv2','MARBL','CanOE','MEDUSA2.1','PISCESv2-IPSL','PISCESv2-CNRM','BFM5.2','COBALTv2'};
mod_order=[10,5,1,9,2,6,11,8,4,3,7,12];

cmap1 = brewermap(13,'Oranges');  cmap2 = brewermap(15,'Greens'); cmap3 = brewermap(11,'Purples');
cmap_mod = [cmap2([5,7,11,13,15],:);cmap1([5,7,11,13],:);cmap3([7,11],:);[0,0,0]];

CoreVar_Names = {'Phytoplankton Biomass',     'Zooplankton Biomass',         'NPP',                  'Grazing Losses',            'EP',                 ...
                 'Specific Division Rate',       'Specific Grazing Rate',          'GSP',                  'Z Loss rate',               'Non-Grazing Losses', ...
                 'Specific Grazing Loss Rate',   'Fraction Small P',            'Fraction Small Z',     'SST' ,                      'MLD' ,               ...
                 'Specific Loss Rate',           'Specific Non-grazing Loss Rate'  'Export Efficeincy'     'Ecosystem Efficeincy'       'Surface [Phyto]' ,  ...
                 'DBM [Phyto]',               'DBM Depth',                   'Dilution Index',       'Diaz Frac',                 'Diatom Frac',    'Non- Diatom Frac'...
                 'dgdT',                      'dgdVD',                       'dgdCC',                'dgdP',                      'dgdZ',...
                 'dngdT',                     'dngdVD',                      'dngdCC',               'dngdP',                     'dngdZ'};


CoreVar_Sym    = {'$\Sigma$P',                '$\Sigma$Z',                   'NPP',                  'G',                         'EP',   ...
                  '$\mu$',                    'g',                           'GSP',                  'Z$_M$',                     'NGL', ...
                  'gp',                       'f$_{SP}$',                    'f_{SZ}',               'SST'                        'MLD', ...
                  'l',                        'ngl'                          'e-rat'                 'eco-rat'};


CoreVar_Inv_Units = {'(mmol $m^{-2}$)',            '(mmol $m^{-2}$)',     '(mmol $m^{-2}$ $d^{-1}$)',    '(mmol $m^{-2}$ $d^{-1}$)',     '(mmol $m^{-2}$ $d^{-1}$)'...
                     '($d^{-1}$)',                 '($d^{-1}$)',          '(mmol $m^{-2}$ $d^{-1}$)',    '(mmol $m^{-2}$ $d^{-1}$)',     '(mmol $m^{-2}$ $d^{-1}$)' ...
                      '($d^{-1}$)',                '(frac)',               '(frac)'  ,                  '($^{\circ}$C)'                     'm'...
                      '($d^{-1}$)',             '($d^{-1}$)',              '(frac)'                        '(frac)'   ,                  '(mmol $m^{-3}$)',...
                      '(mmol $m^{-3}$)',           '(m)',                  '(mmol $m^{-3}$)'            ,  '(frac)',                       '(frac)', '(frac)',...
                      '($d^{-1}$)',              '($d^{-1}$)',             '($d^{-1}$)',               '($d^{-1}$)',                      '($d^{-1}$)',... 
                       '($d^{-1}$)',              '($d^{-1}$)',             '($d^{-1}$)',               '($d^{-1}$)',                      '($d^{-1}$)'};

% Load Model compiles data (Made in SSP_Top_Down_Plots_All
cd /Users/twrohr/Desktop/data.nosync/SSP585
load('All_Model_Sturctures.mat')
end

%% Plot Figure: Global Change Time Series with Boxplot Summary
figure('Units','normalized','Position',[2 2 1 1]); set(gcf,'Color','w');


% === 2. 3x3 Time Series Subplots ===
yalbs = {'$\Delta_{\%} \Sigma \textrm{P} \; ($\%$)$','$\Delta_{\%} \mu  \;  ($\%$)$', '$\Delta_{\%}$ SST  ($\%$)', '$\Delta_{\%}$ MLD ($\%$)'};
vars = [1,6, 14, 15];
label_char = ['a','b','c','d'];
subplot_w = 0.22;
subplot_h = 0.3;
start_x = 0.033; dx = 0.27;
start_y = 0.6;
dy = 0.4;
    
td_mods = [];  bu_mods = [];
for i = 1:4
    % Positioning
    row = ceil(i / 2);
    col = mod(i-1, 2) + 1;
    xpos = start_x + (col-1) * dx;
    ypos = start_y - (row-1) * dy;
    ax = axes('Position', [xpos, ypos, subplot_w, subplot_h]);
    
    for m = 1:11
    kern = 12.*5; % Monthly moving mean
            divdum = movmean(squeeze(CoreVar_Int_pc_ts(m,1,6,:)),12);
            phydum = movmean(squeeze(CoreVar_Int_pc_ts(m,1,1,:)),12);  
            if phydum(end)*divdum(end) < 0
              plot(movmean(squeeze(CoreVar_Int_pc_ts(m,1,vars(i),:)),kern), ...
             'color', cmap_mod(m,:), 'LineWidth', 2, 'LineStyle', ':'); hold on
              if vars(i) == 6; td_mods = [td_mods,m]; end
            else
              plot(movmean(squeeze(CoreVar_Int_pc_ts(m,1,vars(i),:)),kern), ...
             'color', cmap_mod(m,:), 'LineWidth', 2); hold on
              if vars(i) == 6; bu_mods = [bu_mods,m]; end
            end
    end

    mod_smooth = movmean(squeeze(CoreVar_Int_pc_ts(:,1,vars(i),:)),kern,2);
    mod_mean = mean(mod_smooth,1);
    mod_std  = std(mod_smooth,[],1);
    all_std(i)= mean(mod_std(end-5:end));
    all_mean(i)= mean(mod_mean(end-5:end));
    plot(mod_mean, 'k', 'LineWidth', 5);
    fill([1:length(mod_mean), fliplr(1:length(mod_mean))], ...
         [mod_mean+mod_std, fliplr(mod_mean-mod_std)], ...
          [.5 .5 .5], 'FaceAlpha', 0.15, 'EdgeColor', 'none');
    yline(0,'r', 'LineWidth', 1); hold on 
    
    
    %if i > 6
        xticks([60:240:1000,1000]);
        xticklabels(2015 + [60:240:1032]/12);
    %else
    %    xticks([]);
    %end
    if vars(i) == 1 
    ylim([-15 15]);
    elseif vars(i) == 6
    ylim([-20 20]);
    else
    ylim([-22 22]); 
    end
    xlim([13 1000]);
    set(gca, 'FontSize', 18, 'TickLabelInterpreter','latex');
    title(CoreVar_Names{vars(i)}, 'FontSize', 25, 'Interpreter','latex');
    ylabel(yalbs{(i)}, 'FontSize', 20, 'Interpreter','latex');


    % Add subplot label
    text(ax, 0.03, 1.05, ['\textbf{' label_char(i) '}'], ...
        'Units','normalized', 'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', 'Interpreter','latex', 'FontSize', 30);
end

% === 3. Add Legend Below Subplots (Two Rows) ===
model_names = {'iHAMMOCC', 'CMOC', 'WOMBAT', 'OECOv2', 'MARBL', ...
               'CanOE', 'MEDUSA2.1', 'PISCESv2-IPSL', 'PISCESv2-CNRM', 'BFMv5.2', 'COBALTv2'};


%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Add Map
%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % percent top down cells, ensemble
    map_pos= [.57, ypos, .4, subplot_h];
    h3 = axes('Position', map_pos); axesm robinson; axis off; tightmap
    axesm('robinson', 'MapLonLimit', [0 360], 'MapLatLimit', [-80 80], 'Origin', [0 180 0]); 
    gi3 = geoshow(plotvar, Ra, 'DisplayType', 'texturemap');
    colormap(h3,custom_cmap);
    load coastlines; patchm(coastlat,coastlon,	[.7 .7 .7 ]);
    %colormap(cmocean('dense')); colorbar
    clim([0,100]);
    set(gca, 'Color', [0.85 0.85 0.85]);
    Ens_dist = plotvar; 

    % % === Colorbar in its own axis to the right ===
    cb_pos = [map_pos(1)+.02 , map_pos(2)-.06, map_pos(3)-.04, .25];
    cb_ax = axes('Position', cb_pos); colormap(cb_ax, custom_cmap); caxis([1 13]);
    cb = colorbar(cb_ax, 'Ticks', 1:13,  'TickLabelInterpreter', 'latex',  'Location', 'southoutside','orientation','horizontal');

    bin_centers1 = (1:12)+12/13/2; bin_labels = [compose('%d', round(100 * (0:11) / 11))];
    cb = colorbar(cb_ax,   'Ticks', bin_centers1, 'TickLabels', bin_labels, 'TickLabelInterpreter','latex', ...
                 'Location', 'southoutside','orientation','horizontal', 'fontsize', 20);
    cb.Label.String = 'Percent of models with top-down dominance'; cb.Label.Interpreter = 'latex';  cb.Label.FontSize = 25;

    axis off

text(h3, 0.0, 1.05, '\textbf{f}', ...
        'Units','normalized', 'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', 'Interpreter','latex', 'FontSize', 30);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Add Histogram
%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % === Histogram settings ===
hist_ht = subplot_h/2;
hist_top = [map_pos(1), map_pos(2) + map_pos(4) + 0.1 + hist_ht, map_pos(3), hist_ht];
hist_bot = [map_pos(1), map_pos(2) + map_pos(4) + 0.1,            map_pos(3), hist_ht];

% ===========================
% BU histogram (now TOP, pointing UP)
% ===========================
axBU = axes('Position', hist_top); hold on
fill([0 hcmins(2) hcmins(2) 0], [0 0 max([countsTD,countsBU]) max([countsTD,countsBU])], ...
    [0.95 0.95 0.95], 'EdgeColor', 'none');

for b = 1:nbins
    bar(bin_centers(b), countsBU(b), 1.0, ...
        'FaceColor', colorsBU(b,:), 'EdgeColor', [.5 .5 .5], 'LineWidth', 0.2);
end

xlim([-40, 40]);
ylim([0 max([countsTD,countsBU])]);
yticks([]); xticks([]);
%axis off
box on
arrayfun(@(x) xline(x, 'Color', [.7 .7 .7], 'LineWidth', 0.5), -40:20:40);
title({'Ensemble distribution of biomass change'}, 'Interpreter','latex', 'FontSize', 25)

% Summary stats
text(-39, 0.88 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_BU_neg), ...
    'Parent', axBU, 'Interpreter','latex', 'FontSize', 30, 'HorizontalAlignment', 'left');
text(39,  0.88 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_BU_pos), ...
    'Parent', axBU, 'Interpreter','latex', 'FontSize', 30, 'HorizontalAlignment', 'right');
text(axBU, 0.0, 1.1, '\textbf{e}', ...
        'Units','normalized', 'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', 'Interpreter','latex', 'FontSize', 30);


% ===========================
% TD histogram (now BOTTOM, pointing DOWN)
% ===========================
axTD = axes('Position', hist_bot); hold on
fill([0 hcmins(2) hcmins(2) 0], [0 0 -1.01*max([countsTD,countsBU]) -1.01*max([countsTD,countsBU])], ...
    [0.95 0.95 0.95], 'EdgeColor', 'none');

for b = 1:nbins
    bar(bin_centers(b), -countsTD(b), 1.0, ...
        'FaceColor', colorsTD(b,:), 'EdgeColor', [.5 .5 .5], 'LineWidth', 0.2);
end

xlim([-40, 40]);
ylim([-1.02*max([countsTD,countsBU]) 0]);
xticks(-40:20:40);
xticklabels(compose('%d', -40:20:40));
yticks([]);
set(gca, 'TickLabelInterpreter', 'latex', 'FontSize', 20);
xlabel({'Change in phytoplankton biomass', ...
    '$\overline{2015 - 2020} \rightarrow \overline{2095 - 2100}$'}, 'Interpreter','latex', 'FontSize', 25)

arrayfun(@(x) xline(x, 'Color', [.7 .7 .7], 'LineWidth', 0.5), -40:20:40);

% Summary stats
text(-39, -0.9 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_TD_neg), ...
    'Parent', axTD, 'Interpreter','latex', 'FontSize', 30, 'HorizontalAlignment', 'left');
text(39,  -0.9 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_TD_pos), ...
    'Parent', axTD, 'Interpreter','latex', 'FontSize', 30, 'HorizontalAlignment', 'right');


% LEGEND
% BU Legend 
legend_ax = axes('Position', [-.16, 0.13, 0.6, 0.035]);  % Taller axis for 2 rows
axis off; hold on;

% Dummy plots for legend
for m = bu_mods
    plot(0,0, '-', 'Color', cmap_mod(m,:), 'LineWidth', 4, 'DisplayName', model_names{m});
end

% Create legend
lgd = legend('Orientation','horizontal', 'Location','southoutside', 'Box','off','NumColumns', 3);
set(lgd, 'Interpreter','latex', 'FontSize', 16);

% TD Legend 
legend_ax = axes('Position', [0.12, 0.13, 0.6, 0.035]);  % Taller axis for 2 rows
axis off; hold on;

% Dummy plots for legend
for m = td_mods
    plot(0,0, ':', 'Color', cmap_mod(m,:), 'LineWidth', 4, 'DisplayName', model_names{m});
end

% Create  legend
lgd = legend('Orientation','horizontal', 'Location','southoutside', 'Box','off','NumColumns', 3);
set(lgd, 'Interpreter','latex', 'FontSize', 16);


% Save
%cd /Users/twrohr/Desktop/Manuscripts/Rohr_NCC/Figures
%print(gcf, '-dpng', '-r300', 'Fig1.png')