%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 1a-d - Time Series %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script Path
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised/

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 = {'$\Sigma$ Phytoplankton Biomass',     '$\Sigma$ Zooplankton Biomass',         'NPP',                  'Grazing Losses',            'EP',                 ...
                 {'Phytotoplankton Specific';'Division Rate'},   {'Zooplankton Specific';'Grazing Rate'},          'GSP',                  '{Zooplankton Loss Rate}',               'Non-Grazing Losses', ...
                 'Specific Grazing Loss Rate',   {'Small Phytotoplankton';'Compositon'},      {'Small Zooplankton';'Compositon'},     'SST' ,                      'MLD' ,               ...
                 {'Phytotoplankton Specific';'Loss Rate'},           {'Phytoplankton Specific'; 'Non-grazing Loss Rate'}  'Export Efficeincy'     'Ecosystem Efficeincy'       '[Phyto]$_{\mathrm{Surf}}$' ,  ...
                 'DBM [Phyto]',               'DBM Depth',                   '$\overline{\textrm{[Phy]}_{\textrm{Zoo}}}$',       'Diaz Frac',                 'Diatom Frac',    'Non-Diatom Frac'...
                 'dgdT',                      'dgdVD',                       'dgdCC',                'dgdP',                      'dgdZ',...
                 'dngdT',                     'dngdVD',                      'dngdCC',               'dngdP',                     'dngdZ',...
                 'dldT',                      'dldVD',                       'dldCC',                'dldP',                      'dldZ',...
                 'Dilution Length Scale',    'Encounter Index' ,           'Clearance Rate',       'Clearance Rate (VD)',   '$\overline{\textrm{[Phy]}_{\textrm{MLD}}}$', 'Residual', 'Residual(m)',...
                 'Mesozooplantkton Assimilation', 'Trophic Transfer', 'Trophic Transfer w MicroZ' ,'Total Plankton Biomass','Trophic Transfer Efficiency'};


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}$)',                '($\%$)',               '($\%$)'  ,                  '($^{\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}$)'...
                      '(d$^{-1}$)',              '(d$^{-1}$)',             '(d$^{-1}$)',               '(d$^{-1}$)',                      '(d$^{-1}$)' ...
                      'm',                       '(m$^{-1}$)',   '(m$^{-3}$ mmol$^{-1}$ d$^{-1}$)',  '(m$^{-3}$ mmol$^{-1}$ d$^{-1}$)', '(mmol m$^{-3}$)', '(d$^{-1}$)', '(d$^{-1}$)', ...
                      '(mmol m$^{-2}$ d$^{-1}$)','(mmol m$^{-2}$ d$^{-1}$)','(mmol m$^{-2}$ d$^{-1}$)','(mmol m$^{-2}$)' ,'(frac)'}; 

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


end

for Variable_Catalogue = 1
%% Depth-Integrated Biomass Stocks
% 1  | Depth-Integrated Phytoplankton Biomass, 1–150 m (ΣP) | mmol m^-2
% 2  | Depth-Integrated Zooplankton Biomass, 1–150 m (ΣZ) | mmol m^-2
% 52 | Total Plankton Biomass (ΣBiomass = ΣP + ΣZ) | mmol m^-2

%% Bulk Carbon Fluxes
% 3  | Net Primary Production (NPP) | mmol m^-2 d^-1
% 5  | Export Production (EP) | mmol m^-2 d^-1
% 50 | Trophic Transfer (TT; mesozooplankton mortality to higher predation) | mmol m^-2 d^-1
% 4  | Grazing Losses (G) | mmol m^-2 d^-1
% 10 | Non-grazing Losses (NG) | mmol m^-2 d^-1
% 49 | Mesozooplankton Assimilation (GSPMeso; excludes microzooplankton) | mmol m^-2 d^-1
% 8  | Gross Secondary Production (GSP; all zooplankton, excludes assimilation inefficiencies) | mmol m^-2 d^-1
% 9  | Zooplankton Loss Rate (Z_m) | mmol m^-2 d^-1

%% Specific Rate Terms
% 6  | Phytoplankton-specific Division Rate (μ = NPP/ΣP) | d^-1
% 7  | Zooplankton-specific Grazing Rate (g = G/ΣZ) | d^-1
% 11 | Specific Grazing Loss Rate (l_g = G/ΣP) | d^-1
% 17 | Phytoplankton-specific Non-grazing Loss Rate (l_ng = NG/ΣP) | d^-1
% 16 | Phytoplankton-specific Loss Rate (l = l_g + l_ng) | d^-1

%% Community Composition
% 12 | Small Phytoplankton Composition (f_SP) | %
% 13 | Small Zooplankton Composition (f_SZ) | %
% 25 | Diatom Fraction | fraction
% 26 | Non-diatom Fraction | fraction

%% Physics
% 14 | Sea-surface Temperature (SST) | °C
% 15 | Mixed-layer Depth (MLD) | m

%% Carbon-Flux Transfer Efficiencies
% 18 | Export Efficiency (e-rat) | fraction
% 19 | Ecosystem Efficiency (eco-rat) | fraction
% 53 | Trophic Transfer Efficiency (TTE) | fraction

%% Phytoplankton Concentrations
% 23 | Zooplankton Biomass-weighted Phytoplankton Concentration ([Phy]_Zoo; concentration experienced by average zooplankton) | mmol m^-3
% 20 | Surface Phytoplankton Concentration ([Phyto]_Surf) | mmol m^-3
% 46 | Mean Phytoplankton Concentration within MLD ([Phyto]_MLD) | mmol m^-3

%% Vertical Profile
% 42 | Dilution Length Scale (DL) | m
% 21 | Deep Biomass Maximum Phytoplankton Concentration ([Phyto]_DBM) | mmol m^-3
% 22 | Deep Biomass Maximum Depth (Z_DBM) | m

%% Decomposed Loss-Rate Drivers: Grazing
% 27 | Direct Temperature Effect on Grazing Pressure (l_g,T) | d^-1
% 28 | Direct Vertical Redistribution Effect on Grazing Pressure (l_g,VD) | d^-1
% 29 | Direct Community Composition Effect on Grazing Pressure (l_g,CC) | d^-1
% 30 | Direct Phytoplankton Biomass Effect on Grazing Pressure (l_g,ΣP) | d^-1
% 31 | Direct Zooplankton Biomass Effect on Grazing Pressure (l_g,ΣZ) | d^-1

%% Decomposed Loss-Rate Drivers: Non-Grazing
% 32 | Direct Temperature Effect on Non-grazing Pressure (l_ng,T) | d^-1
% 33 | Direct Vertical Redistribution Effect on Non-grazing Pressure (l_ng,VD) | d^-1
% 34 | Direct Community Composition Effect on Non-grazing Pressure (l_ng,CC) | d^-1
% 35 | Direct Phytoplankton Biomass Effect on Non-grazing Pressure (l_ng,ΣP) | d^-1
% 36 | Direct Zooplankton Biomass Effect on Non-grazing Pressure (l_ng,ΣZ) | d^-1

%% Decomposed Loss-Rate Drivers: Total
% 37 | Direct Temperature Effect on Loss Rates (l_T) | d^-1
% 38 | Direct Vertical Redistribution Effect on Loss Rates (l_VD) | d^-1
% 39 | Direct Community Composition Effect on Loss Rates (l_CC) | d^-1
% 40 | Direct Phytoplankton Biomass Effect on Loss Rates (l_ΣP) | d^-1
% 41 | Direct Zooplankton Biomass Effect on Loss Rates (l_ΣZ) | d^-1
% 47 | Indirect Interaction Effects on Loss Rates (l_Int) | d^-1

%% Clearance Rates
% 44 | Clearance Rate (ε) | m^-3 mmol^-1 d^-1
% 45 | Direct Vertical Redistribution Effect on Clearance Rate (ε_VD) | m^-3 mmol^-1 d^-1
end

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


% === 2. 3x3 Time Series Subplots ===
vars = [14, 15,6, 1]; % See Variable_Catalogue above

% Variable Names
 % Dynamically name
  % ylabs = CoreVar_Inv_Names(vars);
 % Custom Names 
   yalbs = {'$\Delta_{\%} \Sigma \textrm{P} \; ($\%$)$','$\Delta_{\%} \mu  \;  ($\%$)$', '$\Delta_{\%}$ SST  ($\%$)', '$\Delta_{\%}$ MLD ($\%$)'};

label_char = ['a','b','c','d'];
subplot_w = 0.25;
subplot_h = 0.4;
start_x = 0.08; 
dx = 0.32;
start_y = 0.55;
dy = 0.5;
    
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.05, 1.07, ['\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'};

