%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Supplemntary Figure 1 - Time Series %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_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');

% ==== Boxplot Axis (LEFT) ====
for Make_BoxPlot=[1]
box_pos = [0.05, 0.1, 0.28, 0.84];

% ==== Time Series Subplot Grid (RIGHT) ====
subplot_w = 0.18;
subplot_h = 0.25;
start_x = 0.38;
dx = 0.205;
start_y = 0.69;
dy = 0.295;

% === 1. Wide Left Boxplot Axis ===
ax_box = axes('Position', box_pos);

CoreVar_Labels_bp = fliplr({'SST', 'MLD', '$\mu$', '$l$', 'Phyto', 'Zoo', 'NPP', 'GSP', 'EP'});
vars_bp = fliplr([14, 15, 6, 16, 1, 2, 3 ,8 ,5]);

xL = [-25 25]; yL = [0 length(vars_bp)+1];
hold(ax_box, 'on');
fill([0 xL(2) xL(2) 0], [yL(1) yL(1) yL(2) yL(2)], ...
     [0.6 0.8 1.0], 'FaceAlpha', 0.1, 'EdgeColor', 'none');
fill([xL(1) 0 0 xL(1)], [yL(1) yL(1) yL(2) yL(2)], ...
     [1.0 0.6 0.6], 'FaceAlpha', 0.1, 'EdgeColor', 'none');
line([0 0], [0 length(vars_bp)+1], 'Parent', ax_box, ...
     'Color', [0.6 0.6 0.6], 'LineStyle', '--', 'LineWidth', 1.2);

data = squeeze(CoreVar_Int_pc(:,1,vars_bp));
% --- Call boxplot (but don't restyle yet)
boxplot(ax_box, data, 'orientation','horizontal', 'colors','k', ...
    'symbol','k+', 'Whisker',1.5, 'Notch','off', 'BoxStyle','outline');
hold(ax_box, 'on');

% --- Redraw patches (white fill, blue edge)
h = findobj(ax_box,'Tag','Box');
for j = 1:length(h)
    patch(get(h(j),'XData'), get(h(j),'YData'), 'w', ...
          'EdgeColor','b', 'LineWidth', 1.2);
end

% --- Manually re-plot median lines (so they appear on top of the patch)
for i = 1:length(vars_bp)
    med = median(data(:,i), 'omitnan');
    plot(ax_box, [med, med], [i-0.25, i+0.25], 'r--', 'LineWidth', 1.2);
end

% --- Add red mean marker
for i = 1:length(vars_bp)
    mu = mean(data(:,i), 'omitnan');
    plot(ax_box, mu, i, 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');
end

set(ax_box, 'XLim', [-25, 25], 'FontSize', 16, 'TickLabelInterpreter','latex');
yticks(ax_box, 1:length(vars_bp));
yticklabels(ax_box, CoreVar_Sym(vars_bp));
xlabel(ax_box, {'Percent change under SSP585 forcing', ...
    '$\overline{2015 - 2020} \rightarrow \overline{2095 - 2100}$'}, ...
    'Interpreter','latex', 'fontsize',18);

% Add subplot label
text(ax_box, 0.03, 1.016, '\textbf{a}', 'Units','normalized', ...
    'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
    'Interpreter','latex', 'FontSize', 18);
end

% === 2. 3x3 Time Series Subplots ===
vars = [14, 6, 1, 15, 16, 2, 3, 8, 5]; % Us efor Sup figure
vars = [14, 14, 15, 6, 16, 1, 3, 8, 5]; % Us efor Sup figure
label_char = 'b':'j';

td_mods = [];  bu_mods = [];
for i = 1:9
    row = ceil(i / 3);
    col = mod(i-1, 3) + 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);
    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:120:1000,1000]);
        xticklabels(2015 + [60:120:1032]/12);
    else
        xticklabels([]);
    end

    if i ==1 || i ==4 || i==7
        ylabel('Percent Change', 'FontSize', 18, 'Interpreter','latex') 
    else
    end

    ylim([-25 25]);
    xlim([30 1000]);
    set(gca, 'FontSize', 14, 'TickLabelInterpreter','latex');
    title(CoreVar_Names{vars(i)}, 'FontSize', 18, 'Interpreter','latex');

    % Add subplot label
    text(ax, 0.03, 1.05, ['\textbf{' label_char(i) '}'], ...
        'Units','normalized', 'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', 'Interpreter','latex', 'FontSize', 18);
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'};



% BU Legend 
legend_ax = axes('Position', [.23, 0.015, 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','north', 'Box','off','NumColumns', 3);
set(lgd, 'Interpreter','latex', 'FontSize', 16);

% TD Legend 
legend_ax = axes('Position', [0.48, 0.015, 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','north', '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')