%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 2 - Global TD_BU Map %%
% NOTE plotvar and ctrlvar (lines 110-11) can be changed
%      to plot the % change in a diffent varibale (plotvar) 
%      or to change the definition of top-down/bottom up control with respect to a vairbale other than divisions rates (ctrlvar)
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs


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 -- 0.45, 0.60, 0.80
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green -- 0.40, 0.60, 0.40
        % --- 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 -- 0.75, 0.30, 0.30
    mustard_rgb = [0.90, 0.7, 0.25];  % strong end member yellow -- 0.80, 0.65, 0.20
        % --- 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];



    %%% Discete
    % Interpolate 6 steps from strong to weak BU
    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)'];
    % Interpolate 6 steps from weak to strong TD
    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)'];
    % Combine
    custom_cmap = [greens; reds];

end

%% === Settings ===
cmins = [-40, 40]; % Color axis limits
hcmins = [-50, 50]; % Color axis limits
nbins = 100;
model_names = {'iHAMMOCC','CMOC','WOMBAT','OECOv2','MARBL','CanOE','MEDUSA2.1',...
               'PISCESv2-IPSL','PISCESv2-CNRM','BFM5.2','COBALTv2','Ensemble'};

% Initialize storage for ensemble
TD_all_vals = []; BU_all_vals = []; TD_stack = []; BU_stack = [];

% === Load data ===
cd /Users/twrohr/Desktop/Data.nosync/SSP585/Abridged/regridded2
mods = dir('*.mat');
mod_order = [10,5,1,9,2,6,11,8,4,3,7];

figure('Units','normalized','Position',[0 0 .85 1.25]); set(gcf,'color','w');

for s = 1:12
    % Load or average data
    if s <= 11
        load(mods(mod_order(s)).name, 'CoreVar_Inv_pc','Lat_t','Lon_t','Bio_Weights')
        area_map = Bio_Weights;  % replace area map with total biomass (inv * area)
        % Pick variable to plot and mask with    
            plotvar = squeeze(CoreVar_Inv_pc(1,:,:)); % 1. biomass
            ctrlvar = squeeze(CoreVar_Inv_pc(6,:,:)); % 6. Div rates; 8. GSP; 3. NPP 5. EP
            maskvar = sign(plotvar .* ctrlvar);
        % Adjust to plot specified geogrpahic range
            %plotvar([1:70,110:180],:)= nan;
            %plotvar([70:110],:)= nan;
        if s == 1
            plot_stack = plotvar;
            ctrl_stack = ctrlvar;
            mask_stack = maskvar;
        else
            plot_stack = cat(3, plot_stack, plotvar);
            ctrl_stack = cat(3, ctrl_stack, ctrlvar);
            mask_stack = cat(3, mask_stack, maskvar);
        end
    else
        mask_stack(mask_stack>0) = 0; 
        plotvar = -100*(sum(mask_stack, 3)./11);
        % plotvar = mean(plot_stack, 3, 'omitnan');
        % ctrlvar = mean(ctrl_stack, 3, 'omitnan');
        % maskvar = sign(plotvar .* ctrlvar);
    end

    BUmask = maskvar; BUmask(BUmask<0) = nan;
    TDmask = maskvar; TDmask(TDmask>0) = nan;

        % ============================================================
        % Median biomass change (all cells / TD-controlled / BU-controlled)
        % Print per model as loop runs
        % ============================================================
        if s < 12
            % All valid biomass-change cells
            all_vals = plotvar(~isnan(plotvar));
        
            % TD-controlled cells (where TDmask is valid)
            td_vals  = plotvar(~isnan(TDmask) & ~isnan(plotvar));
        
            % BU-controlled cells (where BUmask is valid)
            bu_vals  = plotvar(~isnan(BUmask) & ~isnan(plotvar));
        
            med_all = mean(all_vals, 'omitnan');
            med_TD  = mean(td_vals,  'omitnan');
            med_BU  = mean(bu_vals,  'omitnan');
        
            fprintf('%-12s | median ΔBio (all/TD/BU) = %+7.2f / %+7.2f / %+7.2f %%\n', ...
                model_names{s}, med_all, med_TD, med_BU);
        else
            fprintf('%-12s | (ensemble panel) medians not printed here\n', model_names{s});
        end


   % Histogram Values

    % Area Weight each bin.
     bin_edges = linspace(hcmins(1), hcmins(2), nbins+1);
     bin_centers = 0.5 * (bin_edges(1:end-1) + bin_edges(2:end));

        % Count and weight by area
        countsTD = zeros(1, nbins);
        countsBU = zeros(1, nbins);

    % Preallocate storage for per-model normalized histograms
    if s == 1
        TD_frac_all = nan(11, nbins);
        BU_frac_all = nan(11, nbins);
    end
    
    % Initialize
    countsTD_raw = zeros(1, nbins);
    countsBU_raw = zeros(1, nbins);

    if s < 12
        % --- Raw biomass-weighted counts per bin (NOT normalized yet) ---
        for b = 1:nbins
            % TD: only use grid cells where TDmask is valid
            bin_mask_TD = (plotvar >= bin_edges(b)) & (plotvar < bin_edges(b+1)) & ~isnan(TDmask);
            countsTD_raw(b) = nansum( (area_map .* bin_mask_TD), 'all' );
    
            % BU: only use grid cells where BUmask is valid
            bin_mask_BU = (plotvar >= bin_edges(b)) & (plotvar < bin_edges(b+1)) & ~isnan(BUmask);
            countsBU_raw(b) = nansum( (area_map .* bin_mask_BU), 'all' );
        end
    
        % --- Normalize within-model by that model's total biomass weight over valid ocean ---
        totW = nansum(area_map(~isnan(plotvar)), 'all');  % NOTE: biomass-weight total, not geometric area
    
        if totW > 0
            TD_frac_all(s,:) = countsTD_raw / totW;  % sums to ~fraction of model biomass in TD bins
            BU_frac_all(s,:) = countsBU_raw / totW;
        else
            TD_frac_all(s,:) = nan(1, nbins);
            BU_frac_all(s,:) = nan(1, nbins);
        end
    
        % --- What you plot for this model panel (percent) ---
        countsTD = 100 * TD_frac_all(s,:);
        countsBU = 100 * BU_frac_all(s,:);
    
    else
        % --- Ensemble: equal model weighting = mean of per-model normalized histograms ---
        countsTD = 100 * mean(TD_frac_all, 1, 'omitnan');
        countsBU = 100 * mean(BU_frac_all, 1, 'omitnan');
    end

    % Colormaps
    colorsTD = interp1(linspace(cmins(1), cmins(2), size(cmapTD,1)), cmapTD, bin_centers, 'linear', 'extrap');
    colorsBU = interp1(linspace(cmins(1), cmins(2), size(cmapBU,1)), cmapBU, bin_centers, 'linear', 'extrap');

    % === Dynamic position adjustments ===
    h_map = subplot(3,4,s);
    pos = get(h_map, 'Position'); delete(h_map);
    posnew = pos;

    if     s<=4;   posnew(2) = 0.73-.03; 
    elseif s<=8;   posnew(2) = 0.47-.04; 
    elseif s<=12;  posnew(2) = 0.21-.05;
    end

    if     any(s==[2,6,10]); posnew(1) = posnew(1) - 0.02; 
    elseif any(s==[3,7,11]); posnew(1) = posnew(1) - 0.04; 
    elseif any(s==[4,8,12]); posnew(1) = posnew(1) - 0.06;
    end

    map_pos = posnew; hist_ht =.045;
    hist_top = [posnew(1), posnew(2)+map_pos(4)-.035+hist_ht, posnew(3), hist_ht];
    hist_bot = [posnew(1), posnew(2)+map_pos(4)-.035, posnew(3), hist_ht];


    % Plot map with TD overlay
    Ra = georefcells([min( Lat_t(:,1)), max( Lat_t(:,1))], [min(Lon_t(1,:)), max(Lon_t(1,:))], size(plotvar));
    if s<12
    h1 = axes('Position', map_pos); axesm robinson; axis off; tightmap
    axesm('robinson', 'MapLonLimit', [0 360], 'MapLatLimit', [-80 80], 'Origin', [0 180 0]); 
    % gi1 = geoshow(Lat,Lon,plotvar,'DisplayType','texturemap','CData',plotvar,'ZData',Z);
    gi1 = geoshow(plotvar, Ra, 'DisplayType', 'texturemap');
    gi1.AlphaDataMapping = 'none';
    gi1.FaceAlpha = 'texturemap';
    alpha(gi1,double(~isnan(TDmask)))
    colormap(h1,cmapTD); clim(cmins);

    % BU overlay on second axis
    h2 = axes('Position', map_pos); axesm robinson; axis off; tightmap
    axesm('robinson', 'MapLonLimit', [0 360], 'MapLatLimit', [-80 80], 'Origin', [0 180 0]); 
    gi2 = geoshow(plotvar, Ra, 'DisplayType', 'texturemap');
    gi2.AlphaDataMapping = 'none';
    gi2.FaceAlpha = 'texturemap';
    alpha(gi2,double(~isnan(BUmask)))
    colormap(h2,cmapBU); clim(cmins);
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7 ]);
    else
    % 
    % % percent top down cells, ensemble
    % h3 = axes('Position', map_pos); axesm robinson; axis off; tightmap
    % axesm('robinson', 'MapLonLimit', [0 360], '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) + map_pos(3)-.12, map_pos(2)+.043, .15, map_pos(4)-.095];
    % cb_ax = axes('Position', cb_pos); colormap(cb_ax, custom_cmap); caxis([1 12]);
    % cb = colorbar(cb_ax, 'Ticks', 1:12,  'TickLabelInterpreter', 'latex',  'Location', 'eastoutside');
    % bin_centers1 = (1:12) - 0.5; bin_labels = [compose('%d', round(100 * (0:11) / 11))];
    % cb = colorbar(cb_ax,   'Ticks', bin_centers1, 'TickLabels', bin_labels, 'TickLabelInterpreter', 'latex', 'Location', 'eastoutside');
    %  cb.Label.String = 'Models with TD Control ($\%$)'; cb.Label.Interpreter = 'latex'; cb.Label.Rotation = 90;
    % axis off

    end

    % Quadrant Text (percent of non-NaN plotvar)
        % Total ocean area (excluding NaNs)
        total_ocean_area = nansum(area_map(~isnan(plotvar)));
        if s<12
        % Top-down: < 0 and > 0
        p_TD_neg = 100 * nansum(area_map((TDmask < 0) & (plotvar < 0))) / total_ocean_area;
        p_TD_pos = 100 * nansum(area_map((TDmask < 0) & (plotvar > 0))) / total_ocean_area;
        p_TD_neg_all(s)=p_TD_neg;  p_TD_pos_all(s)=p_TD_pos; 
        % Bottom-up: < 0 and > 0
        p_BU_neg = 100 * nansum(area_map((BUmask > 0) & (plotvar < 0))) / total_ocean_area;
        p_BU_pos = 100 * nansum(area_map((BUmask > 0) & (plotvar > 0))) / total_ocean_area;
        p_BU_neg_all(s)=p_BU_neg;  p_BU_pos_all(s)=p_BU_pos;  
        else
         p_TD_neg=mean(p_TD_neg_all); p_TD_pos=mean(p_TD_pos_all);
         p_BU_neg=mean(p_BU_neg_all); p_BU_pos=mean(p_BU_pos_all);
        end

% === Histogram settings ===
if s<12
hist_ht = 0.05;
hist_top = [map_pos(1), map_pos(2) + map_pos(4) - 0.04 + hist_ht, map_pos(3), hist_ht];
hist_bot = [map_pos(1), map_pos(2) + map_pos(4) - 0.04,            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])]);
        axis off
        arrayfun(@(x) xline(x, 'Color', [.7 .7 .7], 'LineWidth', 0.5), -40:20:40);
        title(model_names{s}, 'Interpreter','latex','FontSize',16);

        % Summary stats
        text(-39, 0.82 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_BU_neg), ...
            'Parent', axBU, 'Interpreter','latex', 'FontSize', 16, 'HorizontalAlignment', 'left');
        text(40,  0.82 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_BU_pos), ...
            'Parent', axBU, 'Interpreter','latex', 'FontSize', 16, 'HorizontalAlignment', 'right');
        text(-40, 1.25 * max([countsTD,countsBU]), sprintf('\\textbf{%s}', char('a' + s - 1)), ...
            'Parent', axBU, 'Interpreter','latex', 'FontSize', 20, 'HorizontalAlignment','left');      
        
        % ===========================
        % TD histogram (now BOTTOM, pointing DOWN)
        % ===========================
        axTD = axes('Position', hist_bot); 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), -countsTD(b), 1.0, ...
                'FaceColor', colorsTD(b,:), 'EdgeColor', [.5 .5 .5], 'LineWidth', 0.2);
        end
        
        xlim([-40, 40]);
        ylim([-max([countsTD,countsBU]) 0]);
        axis off
        arrayfun(@(x) xline(x, 'Color', [.7 .7 .7], 'LineWidth', 0.5), -40:20:40);
        
        % Summary stats
        text(-39, -0.85 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_TD_neg), ...
            'Parent', axTD, 'Interpreter','latex', 'FontSize', 16, 'HorizontalAlignment', 'left');
        text(40,  -0.85 * max([countsTD,countsBU]), sprintf('%.0f\\%%', p_TD_pos), ...
            'Parent', axTD, 'Interpreter','latex', 'FontSize', 16, 'HorizontalAlignment', 'right');
else
end

end

%% === Side-by-side Colorbars (Horizontal Only, Taller) ===
% cb_w = 0.277; cb_h = 0.55; base_y = 0.11;
% %cb_w = 0.13; cb_h = 0.55; base_y = 0.11;
% 
% % Top-down colorbar
% cb1_ax = axes('Position', [0.16, base_y, cb_w, cb_h]); axis off
% colormap(cb1_ax, cmapTD); caxis(cmins);
% cb1 = colorbar(cb1_ax, 'southoutside', 'TickLabelInterpreter','latex', 'Orientation','horizontal','fontsize',20);
% set(get(cb1,'Title'), 'String', '$\Delta$ Phytoplankton Biomass driven from top-down ($\%$)', ...
%     'Interpreter','latex','FontSize',18);
% 
% % Bottom-up colorbar
% cb2_ax = axes('Position', [0.54, base_y, cb_w, cb_h]); axis off
% colormap(cb2_ax, cmapBU); caxis(cmins);
% cb2 = colorbar(cb2_ax, 'southoutside', 'TickLabelInterpreter','latex', 'Orientation','horizontal','fontsize',20);
% set(get(cb2,'Title'), 'String', '$\Delta$ Phytoplankton Biomass driven from bottom-up ($\%$)', ...
%     'Interpreter','latex','FontSize',18);

% === Side-by-side Colorbars  ===
% cb_w = 0.27; cb_h = 0.225; base_y = 0.21;
% %cb_w = 0.13; cb_h = 0.55; base_y = 0.11;
% 
% % Top-down colorbar
% cb1_ax = axes('Position', [0.45, base_y, cb_w, cb_h]); axis off
% colormap(cb1_ax, cmapTD); caxis(cmins);
% cb1 = colorbar(cb1_ax, 'eastoutside', 'TickLabelInterpreter','latex','fontsize',20);
% set(get(cb1,'ylabel'), 'String', {'$\Delta$ Phytoplankton Biomass';'driven from top-down ($\%$)'}, ...
%     'Interpreter','latex','FontSize',18);
% 
% % Bottom-up colorbar
% cb2_ax = axes('Position', [0.54, base_y, cb_w, cb_h]); axis off
% colormap(cb2_ax, cmapBU); caxis(cmins);
% cb2 = colorbar(cb2_ax, 'eastoutside', 'TickLabelInterpreter','latex','fontsize',20);
% set(get(cb2,'ylabel'), 'String', {'$\Delta$ Phytoplankton Biomass';'driven from bottom-up ($\%$)'}, ...
%     'Interpreter','latex','FontSize',18);

% === Side-by-side Colorbars  ===
cb_w = 0.17; cb_h = 0.225; base_y = 0.19;
%cb_w = 0.13; cb_h = 0.55; base_y = 0.11;

% Top-down colorbar
cb1_ax = axes('Position', [0.68, base_y+.05, cb_w, cb_h]); axis off
colormap(cb1_ax, cmapTD); caxis(cmins);
cb1 = colorbar(cb1_ax, 'southoutside', 'TickLabelInterpreter','latex','fontsize',14);
set(get(cb1,'title'), 'String', {'$\Delta$ Phytoplankton Biomass ($\%$)'}, ...
    'Interpreter','latex','FontSize',14);

% Bottom-up colorbar
cb2_ax = axes('Position', [0.68, base_y+.135+.05, cb_w, cb_h]); axis off
colormap(cb2_ax, cmapBU); caxis(cmins);
cb2 = colorbar(cb2_ax, 'southoutside', 'TickLabelInterpreter','latex','fontsize',14);
set(get(cb2,'title'), 'String', {'$\Delta$ Phytoplankton Biomass ($\%$)'}, ...
    'Interpreter','latex','FontSize',14);



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


%% Stats on Ensemble (must run with ensembel created last)
% % Logical mask where plotvar equals exactly zero (and is not NaN)
% zero_mask = (Ens_dist >50 ) & ~isnan(Ens_dist);
% 
% % Area where plotvar == 0
% zero_area = nansum(area_map(zero_mask));
% 
% % Percentage of total ocean area
% percent_zero = 100 * zero_area / total_ocean_area;
% 
% fprintf('%.2f%% of the ocean area has plotvar == 0\n', percent_zero);
