%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. 3bcde — Decomposed Lose Rate Distributions
% NOTE 18 Aug. Ensemble correlations use the ensemble-mean spatial Bio_Weights map
%            - These are slightly higher than those that appear in Manusript figure, which doesnt bias correct
%            - Absolute difference annotations need to be vertically adjusted
%            - Line width for zonal plot made slightly smaller 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Paths
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig_dir  = '/Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised';
util_dir = '/Users/twrohr/Desktop/Matlab/Matlab/Utilities';
mod_dir  = '/Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2';

cd(fig_dir)
addpath(genpath(util_dir))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% User controls
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Variables
var_g     = 11;   % grazing pressure, l_g or Delta l_g
var_ng    = 17;   % non-grazing pressure, l_ng or Delta l_ng
var_ltot  = 16;   % total specific loss rate, l or Delta l
var_mu    = 6;    % division rate, mu or Delta mu

% Depth range for contemporary distributions
depth_range = 1:60;

% Optional top-down mask for CHANGE non-grazing data
TDmask2 = false;

% Model selection
model_id = 1:11;
include_ensemble = true;

% -------------------------
% Panel widths
% -------------------------
% Set either zonal width to 0 to skip that zonal column.
contemp_hist_w  = 0.095;
contemp_zonal_w = 0.04;

change_hist_w   = 0.095;
change_zonal_w  = 0.04;

% Horizontal spacing
left_margin       = 0.020;
gap_hist_to_zonal = 0.010;
gap_between_pairs = 0.012;

% Vertical layout
panel_bottom = 0.080;
panel_height = 0.860;
row_gap      = 0.002;

% Figure size
fig_pos = [0.06 0.05 0.86 0.90];

% -------------------------
% Contemporary x-limits
% -------------------------
contemp_xmin = 0;
contemp_xmax = 1;
contemp_xt   = 0:0.2:1.0;

contemp_zxmin = 0;
contemp_zxmax = 0.65;
contemp_zxt   = 0:0.25:0.5;

% CESM-specific contemporary scaling
cesm_model_index = 5;
cesm_contemp_xscale = 1;

contemp_xmax_cesm  = contemp_xmax  * cesm_contemp_xscale;
contemp_zxmax_cesm = contemp_zxmax * cesm_contemp_xscale;

contemp_xt_cesm  = contemp_xt;
contemp_zxt_cesm = 0:0.5:contemp_zxmax_cesm;

% -------------------------
% Change x-limits
% -------------------------
change_xmin = -0.10;
change_xmax =  0.15;
change_xt   = -0.10:0.05:0.15;

change_zxmin = -0.08;
change_zxmax =  0.13;
change_zxt   = -0.05:0.05:0.10;

% KDE grid
nbins = 250;

% Distribution shape controls
rowHalfHeight = 0.42;
cap           = 0.045;
lw            = 1.2;

% Show percentage labels on change histograms
show_change_percentages = true;

% Usually cleaner to label models only once, on far-left contemporary panels
show_contemp_model_labels = false;
show_change_model_labels  = false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model names, order, colors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model_names = {'iHAMMOCC', ...
               'CMOC', ...
               'WOMBAT', ...
               'OECOv2', ...
               'MARBL', ...
               'CanOE', ...
               'MEDUSA2.1', ...
               {'PISCESv2';'IPSL'}, ...
               {'PISCESv2';'CNRM'}, ...
               'BFM5.2', ...
               'COBALTv2'};

% Order of files in /Abridged/regridded2
mod_order = [10,5,1,9,2,6,11,8,4,3,7,12];

% Model colors
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];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Files, rows, latitude and area weights
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd(mod_dir)
mods = dir('*.mat');

numMods = numel(model_id);

if include_ensemble
    nRows = numMods + 1;
else
    nRows = numMods;
end

% Assumes 1 x 1 degree gridded output with [lat x lon] = [180 x 360]
latvec = linspace(-89.5, 89.5, 180)';

Ra = 6.371e6;
deg2rad = pi/180;
dphi = deg2rad;
dlambda = deg2rad;

lat_rad = latvec * deg2rad;
area_per_lat = Ra^2 * dphi * dlambda * cos(lat_rad);
area_map = repmat(area_per_lat, 1, 360);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Column positions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x0 = left_margin;

contemp_hist_left = x0;
x0 = x0 + contemp_hist_w;

if contemp_zonal_w > 0
    x0 = x0 + gap_hist_to_zonal;
    contemp_zonal_left = x0;
    x0 = x0 + contemp_zonal_w;
else
    contemp_zonal_left = NaN;
end

x0 = x0 + gap_between_pairs;

change_hist_left = x0;
x0 = x0 + change_hist_w;

if change_zonal_w > 0
    x0 = x0 + gap_hist_to_zonal;
    change_zonal_left = x0;
    x0 = x0 + change_zonal_w;
else
    change_zonal_left = NaN;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% KDE grids
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contemp_edges   = linspace(contemp_xmin, contemp_xmax, nbins+1);
contemp_centers = 0.5 * (contemp_edges(1:end-1) + contemp_edges(2:end));
x_contemp       = contemp_centers(:);

change_edges   = linspace(change_xmin, change_xmax, nbins+1);
change_centers = 0.5 * (change_edges(1:end-1) + change_edges(2:end));
x_change       = change_centers(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure('Units','normalized','Position',fig_pos);
set(gcf,'Color','w');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Storage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Contemporary ensemble distribution samples
con_g_all  = [];
con_w_all  = [];
con_ng_all = [];
con_wn_all = [];

% Change ensemble distribution samples
chg_g_all  = [];
chg_w_all  = [];
chg_ng_all = [];
chg_wn_all = [];

% Contemporary zonal ensemble
Z_con_g_all     = nan(numel(latvec), numMods);
Z_con_ng_all    = nan(numel(latvec), numMods);
Z_con_ltot_all  = nan(numel(latvec), numMods);
Z_con_mu_all    = nan(numel(latvec), numMods);

% Change zonal ensemble
Z_chg_g_all     = nan(numel(latvec), numMods);
Z_chg_ng_all    = nan(numel(latvec), numMods);
Z_chg_ltot_all  = nan(numel(latvec), numMods);
Z_chg_mu_all    = nan(numel(latvec), numMods);

% Full spatial fields retained for ensemble-mean spatial correlations
spatialStackSize = [size(area_map), numMods];
M_con_ltot_all = nan(spatialStackSize);
M_con_mu_all   = nan(spatialStackSize);
M_chg_ltot_all = nan(spatialStackSize);
M_chg_mu_all   = nan(spatialStackSize);
M_bio_weights_all = nan(spatialStackSize);

% Axes handles
ax_con_h = gobjects(nRows,1);
ax_con_z = gobjects(nRows,1);
ax_chg_h = gobjects(nRows,1);
ax_chg_z = gobjects(nRows,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Per-model processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii = 1:numMods

    mi  = model_id(ii);
    col = cmap_mod(mi,:);

    S = load(mods(mod_order(mi)).name, ...
        'CoreVar_Inv_adj', ...
        'CoreVar_Inv_rw', ...
        'Bio_Weights');

    A = S.Bio_Weights;
    M_bio_weights_all(:,:,ii) = A;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Contemporary fields
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Vc = S.CoreVar_Inv_adj;

    con_g     = squeeze(mean(Vc(var_g,     depth_range, :, :), 2, 'omitnan'));
    con_ng    = squeeze(mean(Vc(var_ng,    depth_range, :, :), 2, 'omitnan'));
    con_ltot  = squeeze(mean(Vc(var_ltot,  depth_range, :, :), 2, 'omitnan'));
    con_mu    = squeeze(mean(Vc(var_mu,    depth_range, :, :), 2, 'omitnan'));

    con_maskG     = isfinite(con_g)     & isfinite(A) & A > 0 & con_g     ~= 0;
    con_maskNG    = isfinite(con_ng)    & isfinite(A) & A > 0 & con_ng    ~= 0;
    con_maskLtot  = isfinite(con_ltot)  & isfinite(A) & A > 0 & con_ltot  ~= 0;
    con_maskMu    = isfinite(con_mu)    & isfinite(A) & A > 0 & con_mu    ~= 0;

    if any(con_maskG(:))
        con_g_all = [con_g_all; con_g(con_maskG)];
        con_w_all = [con_w_all; A(con_maskG)];
    end

    if any(con_maskNG(:))
        con_ng_all = [con_ng_all; con_ng(con_maskNG)];
        con_wn_all = [con_wn_all; A(con_maskNG)];
    end

    con_g_zonal     = local_zonal_mean(con_g,     area_map, con_maskG);
    con_ng_zonal    = local_zonal_mean(con_ng,    area_map, con_maskNG);
    con_ltot_zonal  = local_zonal_mean(con_ltot,  area_map, con_maskLtot);
    con_mu_zonal    = local_zonal_mean(con_mu,    area_map, con_maskMu);

    Z_con_g_all(:,ii)     = con_g_zonal;
    Z_con_ng_all(:,ii)    = con_ng_zonal;
    Z_con_ltot_all(:,ii)  = con_ltot_zonal;
    Z_con_mu_all(:,ii)    = con_mu_zonal;

    M_con_ltot_all(:,:,ii) = con_ltot;
    M_con_mu_all(:,:,ii)   = con_mu;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Change fields
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Vr = S.CoreVar_Inv_rw;

    chg_g     = squeeze(Vr(var_g,     :, :));
    chg_ng    = squeeze(Vr(var_ng,    :, :));
    chg_ltot  = squeeze(Vr(var_ltot,  :, :));
    chg_mu    = squeeze(Vr(var_mu,    :, :));

    BU = squeeze(Vr(1,:,:));
    TD = squeeze(Vr(6,:,:));
    TDmask = sign(BU) .* sign(TD) == -1;

    chg_maskG = isfinite(chg_g) & isfinite(A) & A > 0;

    if TDmask2
        chg_maskNG = isfinite(chg_ng) & isfinite(A) & A > 0 & TDmask;
    else
        chg_maskNG = isfinite(chg_ng) & isfinite(A) & A > 0;
    end

    chg_maskG_data     = chg_maskG  & chg_g  ~= 0;
    chg_maskNG_data    = chg_maskNG & chg_ng ~= 0;
    chg_maskLtot_data  = isfinite(chg_ltot) & isfinite(A) & A > 0 & chg_ltot ~= 0;
    chg_maskMu_data    = isfinite(chg_mu)   & isfinite(A) & A > 0 & chg_mu   ~= 0;

    if any(chg_maskG_data(:))
        chg_g_all = [chg_g_all; chg_g(chg_maskG_data)];
        chg_w_all = [chg_w_all; A(chg_maskG_data)];
    end

    if any(chg_maskNG_data(:))
        chg_ng_all = [chg_ng_all; chg_ng(chg_maskNG_data)];
        chg_wn_all = [chg_wn_all; A(chg_maskNG_data)];
    end

    chg_g_zonal     = local_zonal_mean(chg_g,     area_map, chg_maskG_data);
    chg_ng_zonal    = local_zonal_mean(chg_ng,    area_map, chg_maskNG_data);
    chg_ltot_zonal  = local_zonal_mean(chg_ltot,  area_map, chg_maskLtot_data);
    chg_mu_zonal    = local_zonal_mean(chg_mu,    area_map, chg_maskMu_data);

    Z_chg_g_all(:,ii)     = chg_g_zonal;
    Z_chg_ng_all(:,ii)    = chg_ng_zonal;
    Z_chg_ltot_all(:,ii)  = chg_ltot_zonal;
    Z_chg_mu_all(:,ii)    = chg_mu_zonal;

    M_chg_ltot_all(:,:,ii) = chg_ltot;
    M_chg_mu_all(:,:,ii)   = chg_mu;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Diagnostics: area-weighted mean absolute grazing vs non-grazing
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Means are area-weighted. Significance is assessed with a paired
    % Wilcoxon signed-rank test across paired grid cells using |g|-|ng|.
    % For CHANGE, this compares |Delta grazing| vs |Delta non-grazing|.

    model_label_print = local_label_to_char(model_names{mi});

    local_print_abs_loss_diagnostic( ...
        model_label_print, ...
        'Contemporary', ...
        con_g, con_ng, A, ...
        con_maskG & con_maskNG);

    local_print_abs_loss_diagnostic( ...
        model_label_print, ...
        'Change', ...
        chg_g, chg_ng, A, ...
        chg_maskG_data & chg_maskNG_data);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Contemporary distribution panel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   
    % CESM-specific contemporary x-axis
    if mi == cesm_model_index
        cxmin = contemp_xmin;
        cxmax = contemp_xmax_cesm;
        cxt   = contemp_xt_cesm;
    
        c_edges   = linspace(cxmin, cxmax, nbins+1);
        c_centers = 0.5 * (c_edges(1:end-1) + c_edges(2:end));
        x_contemp_i = c_centers(:);
    else
        cxmin = contemp_xmin;
        cxmax = contemp_xmax;
        cxt   = contemp_xt;
    
        x_contemp_i = x_contemp;
    end
    
    hpos = row_pos_from_grid(ii, nRows, panel_bottom, panel_height, ...
                             row_gap, contemp_hist_left, contemp_hist_w);
    
    ax_con_h(ii) = axes('Position', hpos);
    hold(ax_con_h(ii),'on');
    
    setup_hist_row_axis(ax_con_h(ii), cxmin, cxmax, cxt, ...
                        ii, nRows, false);
    
    plot_bivariate_distribution_row( ...
        ax_con_h(ii), ...
        x_contemp_i, ...
        con_g, A, con_maskG, ...
        con_ng, A, con_maskNG, ...
        col, rowHalfHeight, cap, lw);

    add_abs_mean_difference_label(ax_con_h(ii), cxmin, cxmax, ...
        con_g(con_maskG), A(con_maskG), ...
        con_ng(con_maskNG), A(con_maskNG));
    
    yline(ax_con_h(ii), 0, '-', 'Color', [0.88 0.88 0.88], 'LineWidth', 0.6);
    set(ax_con_h(ii),'Box','on','LineWidth',1.0,'Layer','top');
    
    if show_contemp_model_labels
        add_model_label(ax_con_h(ii), cxmin, cxmax, model_names{mi}, col);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Contemporary zonal panel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
    if contemp_zonal_w > 0
    
        % CESM-specific contemporary zonal x-axis
        if mi == cesm_model_index
            czxmin = contemp_zxmin;
            czxmax = contemp_zxmax_cesm;
            czxt   = contemp_zxt_cesm;
        else
            czxmin = contemp_zxmin;
            czxmax = contemp_zxmax;
            czxt   = contemp_zxt;
        end
    
        zpos = row_pos_from_grid(ii, nRows, panel_bottom, panel_height, ...
                                 row_gap, contemp_zonal_left, contemp_zonal_w);
    
        ax_con_z(ii) = axes('Position', zpos);
        hold(ax_con_z(ii),'on');
    
        setup_zonal_row_axis(ax_con_z(ii), czxmin, czxmax, czxt, ...
                             ii, nRows, false);
    
        %plot(ax_con_z(ii), con_g_zonal,     latvec, '-',  'Color', col, 'LineWidth', 1.8);
        %plot(ax_con_z(ii), con_ng_zonal,    latvec, ':',  'Color', col, 'LineWidth', 1.5);
        plot(ax_con_z(ii), con_ltot_zonal, latvec, '-', 'Color', col, 'LineWidth', 3);
        plot(ax_con_z(ii), con_mu_zonal,   latvec, '-', 'Color', 'k', 'LineWidth', 1.2);

        r_con_spatial = local_corr_spatial_weighted(con_ltot, con_mu, A, ...
            con_maskLtot & con_maskMu);
        add_zonal_corr_text(ax_con_z(ii), r_con_spatial);
    
        set(ax_con_z(ii),'Box','on','LineWidth',1.0,'Layer','top');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Change distribution panel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    hpos = row_pos_from_grid(ii, nRows, panel_bottom, panel_height, ...
                             row_gap, change_hist_left, change_hist_w);

    ax_chg_h(ii) = axes('Position', hpos);
    hold(ax_chg_h(ii),'on');

    setup_hist_row_axis(ax_chg_h(ii), change_xmin, change_xmax, change_xt, ...
                        ii, nRows, true);

    plot_bivariate_distribution_row( ...
        ax_chg_h(ii), ...
        x_change, ...
        chg_g, A, chg_maskG_data, ...
        chg_ng, A, chg_maskNG_data, ...
        col, rowHalfHeight, cap, lw);

    add_abs_mean_difference_label(ax_chg_h(ii), change_xmin, change_xmax, ...
        chg_g(chg_maskG_data), A(chg_maskG_data), ...
        chg_ng(chg_maskNG_data), A(chg_maskNG_data));

    if show_change_percentages
        add_percent_labels(ax_chg_h(ii), change_xmin, change_xmax, ...
            chg_g(chg_maskG_data), A(chg_maskG_data), ...
            chg_ng(chg_maskNG_data), A(chg_maskNG_data));
    end

    yline(ax_chg_h(ii), 0, '-', 'Color', [0.88 0.88 0.88], 'LineWidth', 0.6);
    set(ax_chg_h(ii),'Box','on','LineWidth',1.0,'Layer','top');

    if show_change_model_labels
        add_model_label(ax_chg_h(ii), change_xmin, change_xmax, model_names{mi}, col);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Change zonal panel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if change_zonal_w > 0

        zpos = row_pos_from_grid(ii, nRows, panel_bottom, panel_height, ...
                                 row_gap, change_zonal_left, change_zonal_w);

        ax_chg_z(ii) = axes('Position', zpos);
        hold(ax_chg_z(ii),'on');

        setup_zonal_row_axis(ax_chg_z(ii), change_zxmin, change_zxmax, change_zxt, ...
                             ii, nRows, true);

        %plot(ax_chg_z(ii), chg_g_zonal,     latvec, '-',  'Color', col, 'LineWidth', 1.8);
        %plot(ax_chg_z(ii), chg_ng_zonal,    latvec, ':',  'Color', col, 'LineWidth', 1.5);
        plot(ax_chg_z(ii), chg_ltot_zonal, latvec, '-', 'Color', col, 'LineWidth', 3);
        plot(ax_chg_z(ii), chg_mu_zonal,   latvec, '-', 'Color', 'k', 'LineWidth', 1.2);

        r_chg_spatial = local_corr_spatial_weighted(chg_ltot, chg_mu, A, ...
            chg_maskLtot_data & chg_maskMu_data);
        add_zonal_corr_text(ax_chg_z(ii), r_chg_spatial);
        
        set(ax_chg_z(ii),'Box','on','LineWidth',1.0,'Layer','top');

       
    end

end % model loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ensemble row
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if include_ensemble

    rowY = numMods + 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Contemporary ensemble distribution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    hpos = row_pos_from_grid(rowY, nRows, panel_bottom, panel_height, ...
                             row_gap, contemp_hist_left, contemp_hist_w);

    ax_con_h(rowY) = axes('Position', hpos);
    hold(ax_con_h(rowY),'on');

    setup_hist_row_axis(ax_con_h(rowY), contemp_xmin, contemp_xmax, contemp_xt, ...
                        rowY, nRows, false);

    plot_bivariate_distribution_row_from_samples( ...
        ax_con_h(rowY), ...
        x_contemp, ...
        con_g_all, con_w_all, ...
        con_ng_all, con_wn_all, ...
        'k', rowHalfHeight, cap, lw);

    add_abs_mean_difference_label(ax_con_h(rowY), contemp_xmin, contemp_xmax, ...
        con_g_all, con_w_all, con_ng_all, con_wn_all);

    yline(ax_con_h(rowY), 0, '-', 'Color', [0.88 0.88 0.88], 'LineWidth', 0.6);
    set(ax_con_h(rowY),'Box','on','LineWidth',1.0,'Layer','top');

    if show_contemp_model_labels
        add_model_label(ax_con_h(rowY), contemp_xmin, contemp_xmax, ...
            '\textbf{Ensemble}', 'k');
    end

    local_print_abs_loss_diagnostic_from_samples( ...
        'Ensemble', 'Contemporary', ...
        con_g_all, con_w_all, ...
        con_ng_all, con_wn_all);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Contemporary ensemble zonal
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if contemp_zonal_w > 0

        zpos = row_pos_from_grid(rowY, nRows, panel_bottom, panel_height, ...
                                 row_gap, contemp_zonal_left, contemp_zonal_w);

        ax_con_z(rowY) = axes('Position', zpos);
        hold(ax_con_z(rowY),'on');

        setup_zonal_row_axis(ax_con_z(rowY), contemp_zxmin, contemp_zxmax, contemp_zxt, ...
                             rowY, nRows, false);

        ens_con_g     = mean(Z_con_g_all,     2, 'omitnan');
        ens_con_ng    = mean(Z_con_ng_all,    2, 'omitnan');
        ens_con_ltot  = mean(Z_con_ltot_all,  2, 'omitnan');
        ens_con_mu    = mean(Z_con_mu_all,    2, 'omitnan');

        %plot(ax_con_z(rowY), ens_con_g,    latvec, '-',  'Color', 'k', 'LineWidth', 1.9);
        %plot(ax_con_z(rowY), ens_con_ng,   latvec, ':',  'Color', 'k', 'LineWidth', 1.5);
        plot(ax_con_z(rowY), ens_con_ltot, latvec, '-', 'Color', 'k', 'LineWidth', 3);
        plot(ax_con_z(rowY), ens_con_mu,   latvec, '-', 'Color', 'w', 'LineWidth', 1.5);

        ens_con_ltot_map = mean(M_con_ltot_all, 3, 'omitnan');
        ens_con_mu_map   = mean(M_con_mu_all,   3, 'omitnan');
        ens_bio_weights_map = mean(M_bio_weights_all, 3, 'omitnan');
        r_ens_con_spatial = local_corr_spatial_weighted( ...
            ens_con_ltot_map, ens_con_mu_map, ens_bio_weights_map, ...
            isfinite(ens_con_ltot_map) & isfinite(ens_con_mu_map) & ...
            isfinite(ens_bio_weights_map) & ens_bio_weights_map > 0);
        add_zonal_corr_text(ax_con_z(rowY), r_ens_con_spatial);

        set(ax_con_z(rowY),'Box','on','LineWidth',1.0,'Layer','top');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Change ensemble distribution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    hpos = row_pos_from_grid(rowY, nRows, panel_bottom, panel_height, ...
                             row_gap, change_hist_left, change_hist_w);

    ax_chg_h(rowY) = axes('Position', hpos);
    hold(ax_chg_h(rowY),'on');

    setup_hist_row_axis(ax_chg_h(rowY), change_xmin, change_xmax, change_xt, ...
                        rowY, nRows, true);

    plot_bivariate_distribution_row_from_samples( ...
        ax_chg_h(rowY), ...
        x_change, ...
        chg_g_all, chg_w_all, ...
        chg_ng_all, chg_wn_all, ...
        'k', rowHalfHeight, cap, lw);

    add_abs_mean_difference_label(ax_chg_h(rowY), change_xmin, change_xmax, ...
        chg_g_all, chg_w_all, chg_ng_all, chg_wn_all);

    if show_change_percentages
        add_percent_labels(ax_chg_h(rowY), change_xmin, change_xmax, ...
            chg_g_all, chg_w_all, ...
            chg_ng_all, chg_wn_all);
    end

    yline(ax_chg_h(rowY), 0, '-', 'Color', [0.88 0.88 0.88], 'LineWidth', 0.6);
    set(ax_chg_h(rowY),'Box','on','LineWidth',1.0,'Layer','top');

    if show_change_model_labels
        add_model_label(ax_chg_h(rowY), change_xmin, change_xmax, ...
            '\textbf{Ensemble}', 'k');
    end

    local_print_abs_loss_diagnostic_from_samples( ...
        'Ensemble', 'Change', ...
        chg_g_all, chg_w_all, ...
        chg_ng_all, chg_wn_all);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Change ensemble zonal
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if change_zonal_w > 0

        zpos = row_pos_from_grid(rowY, nRows, panel_bottom, panel_height, ...
                                 row_gap, change_zonal_left, change_zonal_w);

        ax_chg_z(rowY) = axes('Position', zpos);
        hold(ax_chg_z(rowY),'on');

        setup_zonal_row_axis(ax_chg_z(rowY), change_zxmin, change_zxmax, change_zxt, ...
                             rowY, nRows, true);

        ens_chg_g     = mean(Z_chg_g_all,     2, 'omitnan');
        ens_chg_ng    = mean(Z_chg_ng_all,    2, 'omitnan');
        ens_chg_ltot  = mean(Z_chg_ltot_all,  2, 'omitnan');
        ens_chg_mu    = mean(Z_chg_mu_all,    2, 'omitnan');

        
        %plot(ax_chg_z(rowY), ens_chg_g,    latvec, '-',  'Color', 'k', 'LineWidth', 1.9);
        %plot(ax_chg_z(rowY), ens_chg_ng,   latvec, ':',  'Color', 'k', 'LineWidth', 1.5);
        plot(ax_chg_z(rowY), ens_chg_ltot, latvec, '-', 'Color', 'k', 'LineWidth', 3);
        plot(ax_chg_z(rowY), ens_chg_mu,   latvec, '-', 'Color', 'w', 'LineWidth', 1.5);

        ens_chg_ltot_map = mean(M_chg_ltot_all, 3, 'omitnan');
        ens_chg_mu_map   = mean(M_chg_mu_all,   3, 'omitnan');
        ens_bio_weights_map = mean(M_bio_weights_all, 3, 'omitnan');
        r_ens_chg_spatial = local_corr_spatial_weighted( ...
            ens_chg_ltot_map, ens_chg_mu_map, ens_bio_weights_map, ...
            isfinite(ens_chg_ltot_map) & isfinite(ens_chg_mu_map) & ...
            isfinite(ens_bio_weights_map) & ens_bio_weights_map > 0);
        add_zonal_corr_text(ax_chg_z(rowY), r_ens_chg_spatial);

        set(ax_chg_z(rowY),'Box','on','LineWidth',1.0,'Layer','top');
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Group titles above contemporary and change panel pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Y position just above the row-column titles
group_title_y = panel_bottom + panel_height + 0.005;
group_title_h = 0.030;

% Contemporary group span
contemp_group_left  = contemp_hist_left;
contemp_group_right = contemp_hist_left + contemp_hist_w;

if contemp_zonal_w > 0
    contemp_group_right = contemp_zonal_left + contemp_zonal_w;
end

contemp_group_w = contemp_group_right - contemp_group_left;

annotation('textbox', ...
    [contemp_group_left, group_title_y, contemp_group_w, group_title_h], ...
    'String', {'Distribution of $\overline{2015-2022}$';'Phytoplankton Loss Rates (d$^{-1}$)'}, ...
    'Interpreter', 'latex', ...
    'HorizontalAlignment', 'center', ...
    'VerticalAlignment', 'middle', ...
    'EdgeColor', 'none', ...
    'FontSize', 13);

% Change group span
change_group_left  = change_hist_left;
change_group_right = change_hist_left + change_hist_w;

if change_zonal_w > 0
    change_group_right = change_zonal_left + change_zonal_w;
end

change_group_w = change_group_right - change_group_left;

annotation('textbox', ...
    [change_group_left, group_title_y, change_group_w, group_title_h], ...
    'String', {'Distribution of Change in';'Phytoplankton Loss Rates ($\Delta$ d$^{-1}$)'}, ...
    'Interpreter', 'latex', ...
    'HorizontalAlignment', 'center', ...
    'VerticalAlignment', 'middle', ...
    'EdgeColor', 'none', ...
    'FontSize', 13);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Print
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(gcf,'Renderer','painters')
set(gcf,'Color','none')
set(gcf,'InvertHardcopy','off')

copygraphics(gcf, ...
    'ContentType','vector', ...
    'BackgroundColor','none')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pos = row_pos_from_grid(rowIndex, nRows, panelBottom, panelHeight, rowGap, left, width)
% Return normalized axes position for rowIndex, counting from top to bottom.

    rowH = panelHeight / nRows;

    yBottom = panelBottom + panelHeight - rowIndex * rowH + rowGap/2;
    height  = max(rowH - rowGap, 1e-6);

    pos = [left, yBottom, width, height];
end

function setup_hist_row_axis(ax, xmin, xmax, xt, rowIndex, nRows, signed_background)
% Format one boxed bivariate histogram row.

    hold(ax,'on');

    xlim(ax,[xmin xmax]);
    ylim(ax,[-0.5 0.5]);

    if signed_background
        % Positive changes: right side, blue
        patch(ax, ...
            [0 xmax xmax 0], ...
            [-0.5 -0.5 0.5 0.5], ...
            [0.3 0.5 0.85], ...
            'EdgeColor','none', ...
            'FaceAlpha',0.03);

        % Negative changes: left side, red
        patch(ax, ...
            [xmin 0 0 xmin], ...
            [-0.5 -0.5 0.5 0.5], ...
            [0.85 0.3 0.3], ...
            'EdgeColor','none', ...
            'FaceAlpha',0.03);

        xline(ax,0,'k-','LineWidth',1.0);
    end

    yline(ax,0,'-','Color',[0.88 0.88 0.88],'LineWidth',0.6);

    set(ax, ...
        'YTick',[], ...
        'YTickLabel',[], ...
        'XTick',xt, ...
        'TickLabelInterpreter','latex', ...
        'FontSize',10, ...
        'Box','on', ...
        'Layer','top', ...
        'TickDir','in', ...
        'TickLength',[0.010 0.010], ...
        'LineWidth',1.0, ...
        'XGrid','on', ...
        'GridColor',[0.85 0.85 0.85], ...
        'GridAlpha',0.5);

    if rowIndex < nRows
        set(ax,'XTickLabel',[]);
    end
end

function setup_zonal_row_axis(ax, zxmin, zxmax, zxt, rowIndex, nRows, signed_background)
% Format one boxed zonal-average row.

    hold(ax,'on');

    xlim(ax,[zxmin zxmax]);
    ylim(ax,[-75 75]);

    if signed_background
        patch(ax, ...
            [0 zxmax zxmax 0], [-90 -90 90 90], ...
            [0.3 0.5 0.85], ...
            'EdgeColor','none', ...
            'FaceAlpha',0.03);

        patch(ax, ...
            [zxmin 0 0 zxmin], [-90 -90 90 90], ...
            [0.85 0.3 0.3], ...
            'EdgeColor','none', ...
            'FaceAlpha',0.03);

        xline(ax,0,'k-','LineWidth',0.8);
    end

    set(ax, ...
        'XTick',zxt, ...
        'YTick',[-60 0 60], ...
        'YTickLabel',{'-60','0','60'}, ...
        'TickLabelInterpreter','latex', ...
        'FontSize',10, ...
        'Box','on', ...
        'Layer','top', ...
        'TickDir','in', ...
        'XGrid','on', ...
        'YGrid','on', ...
        'GridColor',[0.86 0.86 0.86], ...
        'GridAlpha',0.5, ...
        'LineWidth',1.0);

    if rowIndex < nRows
        set(ax,'XTickLabel',[]);
    end
end

function plot_bivariate_distribution_row(ax, xgrid, g, wg, maskG, ng, wng, maskNG, col, rowHalfHeight, cap, lw)
% Plot one boxed bivariate KDE row from gridded fields.

    if any(maskG(:))
        gvals = g(maskG);
        gw    = wg(maskG);
    else
        gvals = [];
        gw    = [];
    end

    if any(maskNG(:))
        ngvals = ng(maskNG);
        ngw    = wng(maskNG);
    else
        ngvals = [];
        ngw    = [];
    end

    plot_bivariate_distribution_row_from_samples(ax, xgrid, ...
        gvals, gw, ngvals, ngw, col, rowHalfHeight, cap, lw);
end

function plot_bivariate_distribution_row_from_samples(ax, xgrid, gvals, gw, ngvals, ngw, col, rowHalfHeight, cap, lw)
% Plot one boxed bivariate KDE row from sample vectors.

    rowY = 0;

    fG  = zeros(size(xgrid));
    fNG = zeros(size(xgrid));

    if ~isempty(gvals)
        fG = ksdensity(gvals, xgrid, ...
            'Weights', gw, ...
            'Support','unbounded', ...
            'Function','pdf');
    end

    if ~isempty(ngvals)
        fNG = ksdensity(ngvals, xgrid, ...
            'Weights', ngw, ...
            'Support','unbounded', ...
            'Function','pdf');
    end

    mG = max(fG);
    if ~(isfinite(mG) && mG > 0)
        mG = 1;
    end

    mNG = max(fNG);
    if ~(isfinite(mNG) && mNG > 0)
        mNG = 1;
    end

    hG  =  (fG  ./ mG)  * rowHalfHeight;
    hNG = -(fNG ./ mNG) * rowHalfHeight;

    % Grazing pressure: top, filled
    if ~isempty(gvals) && any(isfinite(hG))
        patch(ax, ...
            [xgrid; flipud(xgrid)], ...
            [rowY + hG; rowY * ones(size(xgrid))], ...
            col, ...
            'FaceAlpha',0.70, ...
            'EdgeColor','none');
    end

    % Non-grazing pressure: bottom, translucent
    if ~isempty(ngvals) && any(isfinite(hNG))
        patch(ax, ...
            [xgrid; flipud(xgrid)], ...
            [rowY * ones(size(xgrid)); rowY + flipud(hNG)], ...
            col, ...
            'FaceAlpha',0.10, ...
            'EdgeColor',col, ...
            'LineWidth',0.5);
    end

    % Grazing median/IQR
    if ~isempty(gvals)

        [g_med, g_25, g_75] = local_wquant3(gvals, gw);

        if isfinite(g_med)
            hGm = interp1(xgrid, hG, g_med, 'linear','extrap');

            if isfinite(hGm)
                yIQR = rowY + 0.5*hGm;

                h1 = plot(ax, [g_25 g_75], [yIQR yIQR], ...
                    'Color','k', ...
                    'LineWidth',lw);

                h2 = plot(ax, [g_25 g_25], [yIQR-cap/2 yIQR+cap/2], ...
                    'Color','k', ...
                    'LineWidth',0.7);

                h3 = plot(ax, [g_75 g_75], [yIQR-cap/2 yIQR+cap/2], ...
                    'Color','k', ...
                    'LineWidth',0.7);

                h4 = scatter(ax, g_med, yIQR, 45, 'w', 'filled', ...
                    'Marker','d', ...
                    'MarkerEdgeColor','k', ...
                    'LineWidth',1);

                uistack([h1 h2 h3 h4],'top');
            end
        end
    end

    % Non-grazing median/IQR
    if ~isempty(ngvals)

        [ng_med, ng_25, ng_75] = local_wquant3(ngvals, ngw);

        if isfinite(ng_med)
            hNGm = interp1(xgrid, hNG, ng_med, 'linear','extrap');

            if isfinite(hNGm)
                yIQR = rowY + 0.5*hNGm;

                h1 = plot(ax, [ng_25 ng_75], [yIQR yIQR], ...
                    'Color','k', ...
                    'LineWidth',lw);

                h2 = plot(ax, [ng_25 ng_25], [yIQR-cap/2 yIQR+cap/2], ...
                    'Color','k', ...
                    'LineWidth',0.7);

                h3 = plot(ax, [ng_75 ng_75], [yIQR-cap/2 yIQR+cap/2], ...
                    'Color','k', ...
                    'LineWidth',0.7);

                h4 = scatter(ax, ng_med, yIQR, 45, 'k', 'filled', ...
                    'Marker','d', ...
                    'MarkerEdgeColor','k', ...
                    'LineWidth',1);

                uistack([h1 h2 h3 h4],'top');
            end
        end
    end
end

function add_percent_labels(ax, xmin, xmax, gvals, gw, ngvals, ngw)
% Add positive/negative percentage labels for top and bottom halves.

    fsn = 8;

    x_neg_text = xmin + 0.03 * range([xmin xmax]);
    x_pos_text = xmax - 0.03 * range([xmin xmax]);

    yTop =  0.30;
    yBot = -0.30;

    if ~isempty(gvals)
        f_pos_g = weighted_frac(gvals, gw, gvals > 0);
        f_neg_g = weighted_frac(gvals, gw, gvals < 0);

        if round(100*f_neg_g) > 0
            text(ax, x_neg_text, yTop, sprintf('%d\\%%', round(100*f_neg_g)), ...
                'HorizontalAlignment','left', ...
                'VerticalAlignment','middle', ...
                'Interpreter','latex', ...
                'Color','k', ...
                'FontSize',fsn, ...
                'Clipping','on');
        end

        if round(100*f_pos_g) > 0
            text(ax, x_pos_text, yTop, sprintf('%d\\%%', round(100*f_pos_g)), ...
                'HorizontalAlignment','right', ...
                'VerticalAlignment','middle', ...
                'Interpreter','latex', ...
                'Color','k', ...
                'FontSize',fsn, ...
                'Clipping','on');
        end
    end

    if ~isempty(ngvals)
        f_pos_ng = weighted_frac(ngvals, ngw, ngvals > 0);
        f_neg_ng = weighted_frac(ngvals, ngw, ngvals < 0);

        if round(100*f_neg_ng) > 0
            text(ax, x_neg_text, yBot, sprintf('%d\\%%', round(100*f_neg_ng)), ...
                'HorizontalAlignment','left', ...
                'VerticalAlignment','middle', ...
                'Interpreter','latex', ...
                'Color','k', ...
                'FontSize',fsn, ...
                'Clipping','on');
        end

        if round(100*f_pos_ng) > 0
            text(ax, x_pos_text, yBot, sprintf('%d\\%%', round(100*f_pos_ng)), ...
                'HorizontalAlignment','right', ...
                'VerticalAlignment','middle', ...
                'Interpreter','latex', ...
                'Color','k', ...
                'FontSize',fsn, ...
                'Clipping','on');
        end
    end
end

function add_model_label(ax, xmin, xmax, labelText, col)
% Add rotated color-coded model label to left of histogram panel.

    label_x = xmin - 0.09 * range([xmin xmax]);

    text(ax, label_x, 0, labelText, ...
        'Color', col, ...
        'FontSize',11, ...
        'Interpreter','latex', ...
        'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', ...
        'Rotation',90, ...
        'Clipping','off');
end


function add_abs_mean_difference_label(ax, xmin, xmax, gvals, gw, ngvals, ngw)
% Annotate the relative difference between the area-weighted absolute means.
% The percentage states how much the larger magnitude exceeds the smaller:
% 100 * (larger/smaller - 1). The label always uses > and is placed on
% the grazing (upper) or non-grazing (lower) row according to the larger mean.

    validG  = isfinite(gvals)  & isfinite(gw)  & gw  > 0;
    validNG = isfinite(ngvals) & isfinite(ngw) & ngw > 0;

    if ~any(validG) || ~any(validNG)
        return
    end

    gabs  = abs(gvals(validG));
    gw    = gw(validG);
    ngabs = abs(ngvals(validNG));
    ngw   = ngw(validNG);

    meanAbsG  = sum(gabs(:)  .* gw(:),  'omitnan') / sum(gw(:),  'omitnan');
    meanAbsNG = sum(ngabs(:) .* ngw(:), 'omitnan') / sum(ngw(:), 'omitnan');

    smaller = min(meanAbsG, meanAbsNG);
    larger  = max(meanAbsG, meanAbsNG);

    if ~(isfinite(smaller) && isfinite(larger)) || smaller <= 0
        return
    end

    relativeDifference = 100 * (larger / smaller - 1);

    labelX = xmin - 0.075 * range([xmin xmax]);
    if meanAbsG >= meanAbsNG
        labelY = 0.30;
    else
        labelY = -0.30;
    end
    labelText = sprintf('%.0f\\%% $>$', relativeDifference);

    text(ax, labelX, labelY, labelText, ...
        'Color','k', ...
        'FontSize',8, ...
        'Interpreter','latex', ...
        'HorizontalAlignment','center', ...
        'VerticalAlignment','middle', ...
        'Rotation',90, ...
        'Clipping','off');
end


function label = local_label_to_char(labelIn)
% Convert model label that may be char/string/cellstr into printable char.

    if iscell(labelIn)
        label = strjoin(labelIn(:)', ' ');
    elseif isstring(labelIn)
        label = char(labelIn);
    else
        label = char(labelIn);
    end
end

function local_print_abs_loss_diagnostic(modelLabel, periodLabel, g, ng, weights, mask)
% Print area-weighted mean absolute grazing/non-grazing values and a
% paired significance test across valid paired grid cells.
%
% Test: paired Wilcoxon signed-rank test on |g|-|ng|.
% For CHANGE, pass Delta fields; the function uses abs(Delta).

    valid = isfinite(g) & isfinite(ng) & isfinite(weights) & weights > 0 & mask;

    gabs  = abs(g(valid));
    ngabs = abs(ng(valid));
    w     = weights(valid);

    local_print_abs_loss_diagnostic_core(modelLabel, periodLabel, gabs, w, ngabs, w);
end

function local_print_abs_loss_diagnostic_from_samples(modelLabel, periodLabel, gvals, gw, ngvals, ngw)
% Print diagnostic from two sample vectors. This is used for the ensemble
% row where paired grid-cell structure is not retained.

    gvalid  = isfinite(gvals)  & isfinite(gw)  & gw  > 0;
    ngvalid = isfinite(ngvals) & isfinite(ngw) & ngw > 0;

    gabs  = abs(gvals(gvalid));
    gw    = gw(gvalid);
    ngabs = abs(ngvals(ngvalid));
    ngw   = ngw(ngvalid);

    local_print_abs_loss_diagnostic_core(modelLabel, periodLabel, gabs, gw, ngabs, ngw);
end

function local_print_abs_loss_diagnostic_core(modelLabel, periodLabel, gabs, gw, ngabs, ngw)
% Shared printer. Means are weighted. Significance test is paired Wilcoxon
% where paired vectors have equal lengths, otherwise a two-sample Wilcoxon
% rank-sum test is used as an ensemble-only fallback.

    if isempty(gabs) || isempty(ngabs)
        fprintf('[%s | %s] Insufficient data for grazing/non-grazing diagnostic.\n', ...
            modelLabel, periodLabel);
        return
    end

    mean_abs_g  = sum(gabs(:)  .* gw(:),  'omitnan') ./ sum(gw(:),  'omitnan');
    mean_abs_ng = sum(ngabs(:) .* ngw(:), 'omitnan') ./ sum(ngw(:), 'omitnan');

    if mean_abs_g > mean_abs_ng
        larger = 'grazing';
    elseif mean_abs_ng > mean_abs_g
        larger = 'non-grazing';
    else
        larger = 'equal';
    end

    p = NaN;
    testName = 'not tested';
    nTest = NaN;

    if numel(gabs) == numel(ngabs)
        d = gabs(:) - ngabs(:);
        d = d(isfinite(d) & d ~= 0);
        nTest = numel(d);

        if nTest >= 3 && exist('signrank','file') == 2
            if mean_abs_g > mean_abs_ng
                p = signrank(d, 0, 'tail', 'right');
            elseif mean_abs_ng > mean_abs_g
                p = signrank(d, 0, 'tail', 'left');
            else
                p = signrank(d, 0);
            end

            testName = 'paired Wilcoxon signed-rank on |g|-|ng|';
        end
    else
        % Ensemble fallback: samples were pooled separately, so there is no
        % paired grid-cell structure.
        nTest = min(numel(gabs), numel(ngabs));

        if nTest >= 3 && exist('ranksum','file') == 2
            if mean_abs_g > mean_abs_ng
                p = ranksum(gabs, ngabs, 'tail', 'right');
            elseif mean_abs_ng > mean_abs_g
                p = ranksum(gabs, ngabs, 'tail', 'left');
            else
                p = ranksum(gabs, ngabs);
            end

            testName = 'two-sample Wilcoxon rank-sum on absolute values';
        end
    end

    if isfinite(p)
        sigText = local_sig_text(p);
        fprintf(['[%s | %s] mean(|grazing|)=%.4g, mean(|non-grazing|)=%.4g; ', ...
                 'larger=%s; %s, n=%d, p=%.3g (%s).\n'], ...
            modelLabel, periodLabel, mean_abs_g, mean_abs_ng, larger, ...
            testName, nTest, p, sigText);
    else
        fprintf(['[%s | %s] mean(|grazing|)=%.4g, mean(|non-grazing|)=%.4g; ', ...
                 'larger=%s; significance not tested or Statistics Toolbox function unavailable.\n'], ...
            modelLabel, periodLabel, mean_abs_g, mean_abs_ng, larger);
    end
end

function sigText = local_sig_text(p)
% Short text for p-value threshold.

    if p < 0.001
        sigText = 'p<0.001';
    elseif p < 0.01
        sigText = 'p<0.01';
    elseif p < 0.05
        sigText = 'p<0.05';
    else
        sigText = 'not significant at p<0.05';
    end
end

function r = local_corr_spatial_weighted(ltot, mu, weights, mask)
% Bio_Weights-weighted Pearson correlation across paired spatial cells.

    if nargin < 4 || isempty(mask)
        mask = true(size(ltot));
    end

    ok = isfinite(ltot) & isfinite(mu) & isfinite(weights) & ...
        weights > 0 & mask;

    if nnz(ok) < 3
        r = NaN;
    else
        x = ltot(ok);
        y = mu(ok);
        w = weights(ok);
        sumW = sum(w);
        meanX = sum(w .* x) / sumW;
        meanY = sum(w .* y) / sumW;
        dx = x - meanX;
        dy = y - meanY;
        varX = sum(w .* dx.^2) / sumW;
        varY = sum(w .* dy.^2) / sumW;
        if varX <= 0 || varY <= 0
            r = NaN;
        else
            covXY = sum(w .* dx .* dy) / sumW;
            r = covXY / sqrt(varX * varY);
        end
    end
end

function add_zonal_corr_text(ax, r)
% Display full-spatial Pearson correlation on the zonal panel.

    if isfinite(r)
        str = sprintf('$r=%.2f$', r);
    else
        str = '$r=NA$';
    end

    text(ax, 0.97, 0.06, str, ...
        'Units','normalized', ...
        'HorizontalAlignment','right', ...
        'VerticalAlignment','bottom', ...
        'Interpreter','latex', ...
        'FontSize',8, ...
        'Color',[0.15 0.15 0.15], ...
        'Clipping','on');
end

function z = local_zonal_mean(field, weights, mask)
% Compute weighted zonal mean for [lat x lon] field.

    valid = isfinite(field) & isfinite(weights) & weights > 0 & mask;

    f = field;
    w = weights;

    f(~valid) = NaN;
    w(~valid) = NaN;

    num = sum(f .* w, 2, 'omitnan');
    den = sum(w,     2, 'omitnan');

    z = num ./ den;
    z(den <= 0) = NaN;
end

function [q50,q25,q75] = local_wquant3(x,w)
% Robust weighted 50/25/75 percentiles.

    ok = isfinite(x) & isfinite(w) & (w > 0);
    x  = x(ok);
    w  = w(ok);

    if isempty(x) || sum(w) <= 0
        q50 = NaN;
        q25 = NaN;
        q75 = NaN;
        return
    end

    [xs,ix] = sort(x(:));
    ws = w(:);
    ws = ws(ix);

    cw = cumsum(ws);
    cw = cw ./ cw(end);

    [cu,iu] = unique(cw,'stable');
    xsu = xs(iu);

    if numel(cu) == 1
        q50 = xsu;
        q25 = xsu;
        q75 = xsu;
        return
    end

    q25 = interp1(cu, xsu, 0.25, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.75, 'linear','extrap');
end

function f = weighted_frac(x, w, cond)
% Weighted fraction of samples satisfying cond.

    valid = isfinite(x) & isfinite(w) & (w > 0);

    x    = x(valid);
    w    = w(valid);
    cond = cond(valid);

    den = sum(w);

    if den == 0
        f = NaN;
    else
        f = sum(w(cond)) / den;
    end
end
