%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. S5 — Zonal Drivers Decomposition 
%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))

for startup =1  
% Load regridded lat lon
cd /Users/twrohr/Desktop/data.nosync/SSP585
         load('Longhurst.mat');
          Lon =  LonGrid; Lat = LatGrid;

% Constants
Ra = 6.371e6;                  % Earth radius [m]
deg2rad = pi / 180;
dphi = deg2rad;               % 1 degree in radians
dlambda = deg2rad;

% Latitude vector (center of cells)
lat = linspace(-89.5, 89.5, 180);  % 180 latitudes from -89.5 to +89.5
lat_rad = lat * deg2rad;

% Area per 1x1 degree cell as a function of latitude
area_per_lat = Ra^2 * dphi * dlambda * cos(lat_rad);  % [m^2]

% Repeat across 360 longitudes to make full 180x360 grid
area_map = repmat(area_per_lat', 1, 360);  % [180x360]

% Convert lon for warppign --- Shift longitudes to [0, 360]
        lonvec = linspace(-179.5, 179.5, 360);
        latvec = linspace(-89.5, 89.5, 180);
        lonvec_wrapped = lonvec;
        lonvec_wrapped(lonvec_wrapped < 0) = lonvec_wrapped(lonvec_wrapped < 0) + 360;
        [lonvec_sorted, lon_order] = sort(lonvec_wrapped);

% Colormaps
    % Get full bam colormap (e.g., 256 colors)
    n_colors = 256;

    % --- Define BU Map  ---
    blue_rgb   = [0.35, 0.55, 0.85];  % muted slate blue
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green
        % --- Interpolation ---
        n_half = floor(n_colors/2);
        blue_to_white = [linspace(blue_rgb(1), white_rgb(1), n_half)', ...
                         linspace(blue_rgb(2), white_rgb(2), n_half)', ...
                         linspace(blue_rgb(3), white_rgb(3), n_half)'];
        white_to_green = [linspace(white_rgb(1), green_rgb(1), n_colors - n_half)', ...
                          linspace(white_rgb(2), green_rgb(2), n_colors - n_half)', ...
                          linspace(white_rgb(3), green_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapBU = [blue_to_white; white_to_green];
    
    % --- Define TD map  ---
    red_rgb     = [0.78, 0.25, 0.25];  % soft classy red
    mustard_rgb = [0.90, 0.7, 0.25];   % strong yellow
        % --- Interpolation ---
        red_to_white = [linspace(red_rgb(1), white_rgb(1), n_half)', ...
                        linspace(red_rgb(2), white_rgb(2), n_half)', ...
                        linspace(red_rgb(3), white_rgb(3), n_half)'];
        white_to_mustard = [linspace(white_rgb(1), mustard_rgb(1), n_colors - n_half)', ...
                            linspace(white_rgb(2), mustard_rgb(2), n_colors - n_half)', ...
                            linspace(white_rgb(3), mustard_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapTD = [red_to_white; white_to_mustard];

    %%% Discrete
    green_start = [0.40, 0.60, 0.40];
    green_end   = [0.85, 0.95, 0.85];
    greens = [linspace(green_start(1), green_end(1), 6)', ...
              linspace(green_start(2), green_end(2), 6)', ...
              linspace(green_start(3), green_end(3), 6)'];
    red_start = [0.90, 0.70, 0.70];
    red_end   = [0.70, 0.00, 0.00];
    reds = [linspace(red_start(1), red_end(1), 6)', ...
            linspace(red_start(2), red_end(2), 6)', ...
            linspace(red_start(3), red_end(3), 6)'];
    custom_cmap = [greens; reds];

    % Model Colors
    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',                 ...
                     'Spec. Division Rate',       'Spec. Grazing Rate',          'GSP',                        'Z Loss rate',               'Non-Grazing Losses', ...
                     'Spec. Grazing Loss Rate',   'Fraction Small P',            'Fraction Small Z',           'SST' ,                      'MLD' ,               ...
                     'Spec. Loss Rate',           'Spec. Non-grazing Loss Rate'  'Export Efficeincy'           'Ecosystem Efficeincy'       'Surface [Phyto]' ,  ...
                     'DBM [Phyto]',               'DBM Depth',                   'Dilution Index',             'Diaz Frac',                 'Diatom Frac',           'Non- Diatom Frac'...
                     'Temperature (g)',           'Vertical Dist (g) ',          'Community Composition (g)',  'Phyto. Biomass (g)',        'Zooplankton Biomass (g)',...
                     'Temperature (n)',           'Vertical Dist (n) ',          'Community Composition (n)',  'Phyto. Biomass (n)',        'Zooplankton Biomass (n)',...
                     'Temperature (T)',           'Vertical Dist (T) ',          'Community Composition (T)',  'Phyto. Biomass (T)',        'Zooplankton Biomass (T)',...
                     'Dilution Index',};

     decomp_names = {{'$\Delta l_{T}$'; 'Contribution of' ; 'Temperture'}, {'$\Delta l_{VD}$'; 'Contribution of' ; 'Vertical Distribution'},...
                     {'$\Delta l_{CC}$';'Contribution of' ; 'Community Composition'}, {'$\Delta l_{\Sigma Phy}$';'Contribution of' ; 'Phytoplankton Biomass'},...
                     {'$\Delta l_{\Sigma Zoo}$';'Contribution of';'Zooplankton Biomass'}};
     decomp_symbols = {{'$\Delta l_{T}$ (d$^{-1}$)'}, {'$\Delta l_{VD}$ (d$^{-1}$)'},...
                     {'$\Delta l_{CC}$ (d$^{-1}$)'}, {'$\Delta l_{\Sigma Phy}$ (d$^{-1}$)'},...
                     {'$\Delta l_{\Sigma Zoo}$ (d$^{-1}$)'}};
    
     decomp_names2 =  { '$\frac{\Delta \; l_{T}}{\Delta\;  \mathrm{SST}}$' , '$\frac{\Delta \; l_{VD} }{\Delta \; \mathrm{D}_L}$' ,...
                    '$\frac{\Delta \; l_{CC} }{\Delta\;  f_{\mathrm{SP}}}$' , '$\frac{\Delta \; l_{\Sigma Phy}}{\Delta \; \Sigma\,{\mathrm{Phy}}}$' ,...
                     '$\frac{\Delta \; l_{\Sigma Zoo} }{\Delta\;  \Sigma \mathrm{Zoo} }$' };

     decomp_units = {'(d$^{-1}$ $^{\circ}$C warmer $^{-1}$)', '(d$^{-1}$ m $^{-1}$)', '(d$^{-1}$ $\%$ inc. small phy. $^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)'};

     Driver_label = {'$\Delta$ SST ($^\circ$C)',' $\Delta$ Dilution Length Scale ($D_{L}$; m)', ' $\Delta$ Fraction Sm. Phytoplankton ($f_{sp}$; $\%$)', '$\Delta$ $\Sigma$ Phytoplankton (mmolC m$^{-2}$)', '$\Delta$ $\Sigma$ Zooplankton (mmolC m$^{-2}$)'};

     decomp_names =decomp_names([4,5,1,3,2]);
     decomp_names2 =decomp_names2([4,5,1,3,2]);
     Driver_label =Driver_label([4,5,1,3,2]);
     decomp_units =decomp_units([4,5,1,3,2]);
     decomp_symbols =decomp_symbols([4,5,1,3,2]); 
     CoreVar_Names{37} = '-';

     wedge_colors = [
    0.90 0.55 0.25;   % SST     → soft orange
    0.35 0.60 0.85;   % MLD     → clearer blue
    0.25 0.75 0.65;   % Diatoms → more greenish teal
    0.60 0.75 0.45;   % Phyto   → sage green
    0.85 0.55 0.75    % Zoo     → soft magenta
];
lighter_colors = wedge_colors + 0.6 .* (1 - wedge_colors);  % shift toward white, preserve hue
   
end





 %% Just Decomposition
  % Zonal Average — 12-panel "Loss Mechanism" only
clc; close all

latims = [-70, 70];
lat     =  latims(1):latims(2);
latlim  =  90+latims(1) : 90+latims(2);   % index into 180-lat zonal arrays

figure('Units','normalized','Position',[0 0 .5 1],'Color','w');

% --- Load zonal stats and names
cd /Users/twrohr/Desktop/data.nosync/SSP585
load('All_Stats_Raw.mat')                   % provides zon_stats
NM           = resolve_names_from_stats(zon_stats);
modelNames   = NM.modelNames(model_id);
modelFields  = NM.modelFields(model_id);

% --- indices we need (as in your script)
idx_vars = [27:36, 47, 17, 16, 11, 1, 2];   % (keep same order you used)

% --- colormaps you already have in workspace; if not, define placeholders
if ~exist('wedge_colors','var') || size(wedge_colors,1)<5
    wedge_colors = repmat([0.30 0.50 0.85],5,1);   % fallback blue-ish
end
if ~exist('lighter_colors','var') || size(lighter_colors,1)<5
    lighter_colors = repmat([0.85 0.30 0.30],5,1); % fallback red-ish (dotted)
end

% panel counter
pcount = 0;

% ======================== Loop over 12 panels ============================
for s = 1:12
    pcount = pcount + 1;
    ax = subplot(4,3,pcount); hold(ax,'on');
    set(ax,'TickLabelInterpreter','latex','FontSize',12,'Box','on','YGrid','on');

    % -------- build Plotvar for this panel (model or ensemble) ----------
    clear Plotvar Ens_Vars
    if s==12
        % ---- Ensemble across the 11 models ----
        for m = 1:11
            mfld = modelFields{m};
            for f = 1:numel(idx_vars)
                Plotvar_tmp = reshape(zon_stats.(NM.varFields{idx_vars(f)}).(mfld).median, [], 1);
                Ens_Vars(m, idx_vars(f), :) = Plotvar_tmp;  
            end
        end
        Plotvar = squeeze(mean(Ens_Vars, 1, 'omitnan'));  
        panel_title = 'Ensemble';
    else
        % ---- Single model ----
        mfld = modelFields{s};
        for f = 1:numel(idx_vars)
            Plotvar(idx_vars(f),:) = reshape(zon_stats.(NM.varFields{idx_vars(f)}).(mfld).median, [], 1);
        end
        panel_title = model_names{s};
    end

    % ====================== Plot Loss Mechanism only ======================
    lw=2;
    % Grazing pressure (27:31)
    for ii = 27:31
        plot(ax, squeeze(Plotvar(ii,latlim)), lat, ...
            'Color', wedge_colors(ii-26,:), 'LineWidth', lw);
    end

    % Non-grazing pressure (32:36)
    for ii = 32:36
        plot(ax, squeeze(Plotvar(ii,latlim)), lat, ...
            'Color', lighter_colors(ii-31,:), 'LineStyle', ':', 'LineWidth', lw);
    end

    % Residual (47)
    plot(ax, squeeze(Plotvar(47,latlim)), lat, 'k', 'LineWidth', 1);
    %plot(ax, squeeze(Plotvar(16,latlim)), lat, 'r', 'LineWidth', lw);

    % Axes cosmetics
    xline(ax, 0, 'k-', 'LineWidth', 1.5);
    ylim(ax, latims); yticks(ax, -60:20:60);
    dum =Plotvar([27:36,47],latlim);  dum=dum(:);
    xlims = [ min(dum(:)) , max(dum(:)) ];
    xlim([-max(abs(xlims)) - .001 , max(abs(xlims))+.001 ])
    
    % Titles / labels
      yl = ylim(ax);xl = xlim(ax);
      text(ax,  mean(xl),yl(2) + 0.15*diff(yl), panel_title, ...
                'HorizontalAlignment','center','VerticalAlignment','top', ...
                'FontSize',18,'Interpreter','latex','Color','k'); %cmap_mod(s,:))
    if ismember(pcount, [1 4 7 10])   % first column → y label
        ylabel(ax, 'Latitude ($^\circ$)', 'Interpreter','latex','FontSize',16);
    end
    if ismember(pcount, [10:12])                   % bottom row → x label
        xlabel(ax, '$\Delta l_i$ (d$^{-1}$)', 'Interpreter','latex','FontSize',16);
    end
end

% Tighten layout a bit
set(gcf,'Renderer','painters');



