%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. 4 Plot Loss Rates
%%%%%%%%%%%%%%%%%%%%%%%%%
% Need to change var id to plot panel d of Fig 4. - and adjust axis limits
% Modified to plot any vairbale. 
% Modifued to include stats and ensemble
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_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 -- 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];

    % 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}$'; 'First-order Effect of' ; 'Temperture'}, {'$\Delta l_{VD}$'; 'First-order Effect of' ; 'Vertical Distribution'},...
                     {'$\Delta l_{CC}$';'First-order Effect of' ; 'Community Composition'}, {'$\Delta l_{\Sigma P}$';'First-order Effect of' ; 'Phytoplankton Biomass'},...
                     {'$\Delta l_{\Sigma Z}$';'First-order Effect of';'Zooplankton Biomass'}, {'$\Delta l_{I}$';'First-order Effect of';'Interaction Effects'}};
     decomp_symbols = {{'$\Delta l_{T}$ (d$^{-1}$)'}, {'$\Delta l_{VD}$ (d$^{-1}$)'},...
                     {'$\Delta l_{CC}$ (d$^{-1}$)'}, {'$\Delta l_{\Sigma P}$ (d$^{-1}$)'},...
                     {'$\Delta l_{\Sigma Z}$ (d$^{-1}$)'},{'$\Delta l_{Int}$ (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 P}}{\Delta \; \Sigma\,{\mathrm{P}}}$' ,...
                     '$\frac{\Delta \; l_{\Sigma Z} }{\Delta\;  \Sigma \mathrm{P} }$' };

     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}$)','dummy'};


end



%% === Plot One Variable + Right Zonal Panel ===
 %%%% Pick plotting vairables
         %   1      2     3    4    5    6    7   8      9     10   11   12   13    14    15    16    17     18      19      20     21      22      23       24      25      26      
         % P Bio  Z bio  NPP   G    EP   u    g  GSP   Zoo_M   NGL  gp   fsP  fsZ   SST   MLD    l    ngl   e-rat  eco-rat  [P]S  [P]DBM  Z_DBM  [phy_zoo]  fdaiz  fdiat  fnondiat 
      
         %   27    28    29     30    31    32-36   37-41     42         43        44            45        46       70         71       72      100    102       119   
         %  dgdT  dgdVD dgdCC  dgdP  dgdZ   dngldx  dldx    Dil Dep  co-loc in  Clear Rate   Clear (VD) [Phy]MLD Mechanism  ngl & gp  P and Z  circle  Stbail   Zoo_m     
 
% varibale for map and zonal 
mvar_id  = 16;              
% variables for violin
    % g and ng
    % var_id   = mvar_id -10;        % Grazing - 10         
    % var_id2  = mvar_id -5;      % non-graizng - 5
    % TDmask2  = true;

    % % other 
    var_id   =11;         % 11 for Fig 4d,3; 16 For Sup Fig 5.       
    var_id2  =17;          % 17 for Fig4d;     6 For Sup Fig 5.
    TDmask2  = false;

model_id = 1:11;                % which models to include

% bound
vi_pr = [7.5 92.5];   % Percentile bounds for Violin plots
mp_pc = [5,95]; 
zn_pr = [.3 99.8];

cd /Users/twrohr/Desktop/data.nosync/SSP585
load('All_Stats_Norm.mat')       % for l 
%load('All_Stats_Raw.mat')         % 

figure('Units','normalized','Position',[0 0 1 1.5]); 
set(gcf,'color','w');

% Resolve names from your stats struct that holds zonal results
NM      = resolve_names_from_stats(zon_stats);   % same helper as before
vidx    = mvar_id;
vField  = NM.varFields{vidx};
%vLabel  = decomp_symbols{mvar_id-36};     % latex-ready variable name
vLabel  = CoreVar_Names{var_id};     % latex-ready variable name

% Try to grab a generic variable name from CoreVar_Names if present
baseName = vLabel;

    
map_left  = 0.06;      % map left
map_width = 0.48;      % map width
map_height= 0.40;      % map height
map_top   = 0.95;      % top edge of the map (so y = top - height)
gap_h     = 0.035;     % gap between map and zonal panel
zonal_w   = 0.22;      % Zonal panel width


%%%%%%%%%%%%%%%%%%%%%
%%% --- MAP ----
%%%%%%%%%%%%%%%%%%%%%
for Plot_Map = [1]
    map_pos0  = [map_left, map_top-map_height, map_width, map_height];
    axm = axes('Units','normalized','Position',map_pos0, ...
               'PositionConstraint','innerposition');  % colorbar won't resize the map
    axes(axm); axesm robinson; axis off; tightmap
    axesm('robinson','MapLonLimit',[0 360],'MapLatLimit',[-75 75],'Origin',[0 180 0]);

    % Create Map and Raster 
    plotvar = squeeze(InterMod_mean_maps(vidx,:,:));
    Ra = georefcells([min(Lat_t(:,1)) max(Lat_t(:,1))], ...
                     [min(Lon_t(1,:)) max(Lon_t(1,:))], size(plotvar));

    geoshow(plotvar, Ra, 'DisplayType','texturemap');
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7]);

    % --- Dynamic colors (reverse the colormap, symmetric around 0) ---
    pv = plotvar(:);
    pv = pv(isfinite(pv));
    if isempty(pv)
        clim([-1 1]);   % fallback
    else
        pr = prctile(pv,mp_pc);
        cmins = [-30 30];
        clim(cmins);
    end
    colormap(axm, flipud(cmocean('balance')));

    % -------- Colorbar BELOW the map on a separate strip --------
    drawnow;                       % let tightmap settle the final position
    map_pos = get(axm,'Position'); % [x y w h] normalized

    cb_h  = 0.025;   % colorbar height
    cb_gap= 0.018;   % gap between map and colorbar

    
    cb = colorbar(axm,'southoutside');   
    cb.Units = 'normalized'; 
    cb.FontSize = 18; 
    cb.TickLabelInterpreter = 'latex';
    cb.Position = [map_pos(1)+.11, map_pos(2)-cb_gap-cb_h, map_pos(3)-.22, cb_h];

    pos = cb.Position;  % [x y w h], normalized
    % Dynamic label based on variable name
    annotation('textbox', [pos(1) pos(2)+pos(4) pos(3) 0.02], ...
               'String', '$\Delta_{\%} l (\%)$', ...
               'Interpreter','latex', ...
               'HorizontalAlignment','center', ...
               'VerticalAlignment','bottom', ...
               'EdgeColor','none', ...
               'FontSize', 18);
end


%%%%%%%%%%%%%%%%%%%
%%% Zonal Averages 
%%%%%%%%%%%%%%%%%%%
for Zonal_plot = [1]

% --- Frame → exact globe vertical span (map data units) ---
setm(axm,'Frame','on');
hFrame = handlem('Frame');  if isempty(hFrame), hFrame = framem; end
if numel(hFrame) > 1, hFrame = hFrame(1); end
Yf = get(hFrame,'YData'); if iscell(Yf), Yf = Yf{1}; end
ybot = min(Yf);  ytop = max(Yf);

% Map axes position -> zonal axes position (normalized units)
map_pos = get(axm,'Position');     % [x y w h]
yl      = get(axm,'YLim');         % map-data Y
toNorm  = @(y) map_pos(2) + map_pos(4) * (y - yl(1)) / (yl(2)-yl(1));
ybot_n  = toNorm(ybot);  ytop_n = toNorm(ytop);

gap_h   = 0.015;                                  % horizontal gap to the right panel
zonal_x = map_pos(1) + map_pos(3) + gap_h;
zon_pos = [zonal_x, ybot_n, zonal_w, max(ytop_n - ybot_n, 1e-9)];

axz = axes('Units','normalized','Position',zon_pos,'Clipping','on'); 
hold(axz,'on');

% --- Project latitude centers to same Y as the map (for vertical alignment) ---
mstruct  = gcm(axm);
lat_cent = latvec(:);
[~, y_cent] = mfwdtran(mstruct, lat_cent, zeros(size(lat_cent)));

% --- Collect model curves (Nlat x Nmod) from zonal stats (medians) ---
Nlat = numel(lat_cent);
mods = numel(model_id);
Z    = nan(Nlat, mods);
Zstd    = nan(Nlat, mods);
Z2   = nan(Nlat, mods);   % companion field (e.g., division)

for ii = 1:mods
    mi = model_id(ii); 
    mField = NM.modelFields{mi};
    if isfield(zon_stats.(vField), mField) && ...
       isfield(zon_stats.(vField).(mField),'mean')
        Z(:,ii) = zon_stats.(vField).(mField).mean(:);
        Zstd(:,ii) = zon_stats.(vField).(mField).std(:);
    end
    % Division (or other companion) — guard field existence
    if isfield(zon_stats,'Phytotoplankton_Specific_Division_Rate') && ...
       isfield(zon_stats.Phytotoplankton_Specific_Division_Rate, mField) && ...
       isfield(zon_stats.Phytotoplankton_Specific_Division_Rate.(mField),'median')
        Z2(:,ii) = zon_stats.Phytotoplankton_Specific_Division_Rate.(mField).median(:); 
    end
end

if all(isnan(Z),'all')
    set(axz,'YLim',[ybot ytop]);
    ylabel(axz,'Latitude','Interpreter','latex');
    title(axz,'No zonal data','Interpreter','latex');
    continue
end

% --- Ensemble mean ± std across models (per latitude) ---
ens_mean  = mean(Z, 2, 'omitnan');
ens_med   = median(Z, 2, 'omitnan');
%ens_std   = std(Z, 0, 2, 'omitnan');
%ens_std = sqrt( sum(Zstd.^2, 2, 'omitnan') ) / 11;
ens_std = sqrt( std(Z, 0, 2, 'omitnan').^2 + (sum(Zstd.^2,2,'omitnan')/11^2) ); % Ensemble total uncertinaty (inter and intra model std)

ens_mean2 = mean(Z2, 2, 'omitnan');

% --- Dynamic X limits (based on actual data) ---
xlim(axz, [-50, 50]);
set(axz, 'YLim', [ybot ytop]);      % perfect vertical alignment with map

% ===== Background shading =====
xl = xlim(axz);  
patch('XData', [xl(1) 0 0 xl(1)], 'YData', [ybot ybot ytop ytop], ...
      'FaceColor',[0.85 0.30 0.30],'FaceAlpha',0.04,'EdgeColor','none','Parent',axz);
patch('XData', [0 xl(2) xl(2) 0], 'YData', [ybot ybot ytop ytop], ...
      'FaceColor',[0.30 0.50 0.85],'FaceAlpha',0.04,'EdgeColor','none','Parent',axz);


% ===== Individual models =====
for ii = 1:mods
    mi = model_id(ii);
    if ~all(isnan(Z(:,ii)))
        plot(axz, Z(:,ii), y_cent, 'LineWidth', 1, 'Color', cmap_mod(mi,:));
    end
end

% % ===== ±1&2 σ band + mean =====
 good = isfinite(ens_mean) & isfinite(ens_std) & isfinite(y_cent);
 if any(good)
    % 1 STD
    X1 = [ens_mean(good)-ens_std(good); flipud(ens_mean(good)+ens_std(good))];
    Y1 = [y_cent(good);                  flipud(y_cent(good))];
    patch('XData',X1,'YData',Y1,'FaceColor',[0 0 0], ...
          'FaceAlpha',0.2,'EdgeColor','none','Parent',axz);
   
    % Means
     plot(axz, ens_mean,  y_cent, 'k-', 'LineWidth', 6);
     plot(axz, ens_mean2, y_cent, 'b-', 'LineWidth', 2);
 end

% ===== Axis cosmetics =====
% find nearest "nice" round step: 1, 2, 5 × 10^k
tickStep = 25;

% compute rounded tick limits
tickMin = -50;
tickMax = 50 ;

% final ticks
xt = tickMin : tickStep : tickMax;
set(axz,'XTick',xt);

ticklats = (-80:20:80)';
[~, ytick_map] = mfwdtran(mstruct, ticklats, zeros(size(ticklats)));
set(axz,'YTick', ytick_map, ...
        'YTickLabel', arrayfun(@(v)sprintf('%g',v), ticklats, 'uni',0), ...
        'Layer','top','FontSize',16,'TickLabelInterpreter','latex', ...
        'Box','on','YGrid','on','XGrid','on', ...
        'XAxisLocation','top');
xline(axz, 0, '--', 'Color', 'k', 'LineWidth', 1.2);

% =======================================================================
% (2) Small axis directly BELOW the zonal plot for global medians “ticks”
% =======================================================================
glob_med = nan(mods,1);
for ii = 1:mods
    mi    = model_id(ii);
    mFld  = NM.modelFields{mi};
    if isfield(mod_stats.(vField), mFld) && ...
       isfield(mod_stats.(vField).(mFld),'mean')
        glob_med(ii) = mod_stats.(vField).(mFld).mean;
    end
end

% Narrow tick axis (touching, same width)
posZ    = get(axz,'Position');
hTickAx = 0.022;         % height of tick strip (small)
gapTick = 0.00;          % tiny gap between zonal and tick strip
axTick  = axes('Units','normalized', ...
               'Position',[posZ(1), posZ(2)-gapTick-hTickAx, posZ(3), hTickAx], ...
               'Color','none','Box','off');
set(axTick,'XLim',xlim(axz),'YLim',[0 1], ...
           'XColor','none','YColor','none','HitTest','off');  % invisible frame
hold(axTick,'on');

% Draw the short colored ticks
for ii = 1:mods
    mi = model_id(ii);
    if ~isfinite(glob_med(ii)), continue; end
    plot(axTick, [glob_med(ii) glob_med(ii)], [0 0.8], ...
         '-', 'Color', cmap_mod(mi,:), 'LineWidth', 1.5);
end

% =======================================================================
% (3) Label strip below ticks: names evenly spaced, sorted by median
% =======================================================================
valid   = isfinite(glob_med);
idxKeep = find(valid);
[~, ord] = sort(glob_med(valid), 'ascend');
idxSrt  = idxKeep(ord);                    % model indices in ascending median
nLab    = numel(idxSrt);
if nLab > 0
    hLabAx   = 0.034;     % label strip height (slightly taller)
    gapLab   = 0.002;     % tiny gap between tick strip and label strip
    extraDrop = 0.030;    % additional drop

    axLab  = axes('Units','normalized', ...
                  'Position',[posZ(1), ...
                              posZ(2) - gapTick - hTickAx - gapLab - hLabAx - extraDrop, ...
                              posZ(3), hLabAx], ...
                  'Color','none','Box','off');
    set(axLab,'XLim',xlim(axz),'YLim',[0 1], 'XColor','none','YColor','none', 'HitTest','off');
    hold(axLab,'on');

    % label x-positions: evenly spaced across the panel width
    xl   = xlim(axz);
    pad  = 0.03*range(xl);
    xlab = linspace(xl(1)+pad, xl(2)-pad, nLab);

    for j = 1:nLab
        mi   = model_id(idxSrt(j));
        name = model_names{mi};
        col  = cmap_mod(mi,:);

        text(axLab, xlab(j), 0.01, name, ...
             'Color', col, 'FontSize', 14, 'Interpreter','latex', ...
             'HorizontalAlignment','center', 'VerticalAlignment','bottom', ...
             'Rotation', 90);

        [x1n, y1n] = axData2norm(axTick, glob_med(idxSrt(j)), 0.0);
        [x2n, y2n] = axData2norm(axLab,  xlab(j), 0.98);
        annotation('line',[x1n x2n],[y1n y2n+.007], 'Color', col, 'LineWidth', 0.8);
    end
end

end % Zonal_plot (for)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 3. Add Violin Distribution (area/biomass-weighted)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/data.nosync/SSP585

%load('All_Stats_Norm.mat')       % for l 
load('All_Stats_Raw.mat')         % 


for Violin_dist = 1
    % ============================================
    % First pass: determine if any data exist and
    % compute dynamic y-limits from pooled values.
    % ============================================
    cd /Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2
    mods    = dir('*.mat');
    numMods = numel(model_id);

    global_min = inf;
    global_max = -inf;
    hasData    = false;

    for ii = 1:numMods
        S = load(mods(mod_order(model_id(ii))).name, ...
                 'CoreVar_Inv_rw', 'Bio_Weights');
        V = S.CoreVar_Inv_rw;
        A = S.Bio_Weights;

        g  = squeeze(V(var_id ,:,:));      % full-field
        ng = squeeze(V(var_id2 ,:,:));      % TD-only subset on right
    
        BU = squeeze(V(1,:,:));
        TD = squeeze(V(6,:,:));
        TDmask = sign(BU) .* sign(TD) == -1;
    
        % basic validity (finite + positive weight)
        maskG  = isfinite(g)  & isfinite(A) & (A > 0);
        if TDmask2
        maskNG = isfinite(ng) & isfinite(A) & (A > 0) & TDmask;
        else
        maskNG = isfinite(ng) & isfinite(A) & (A > 0);
        end
    
        % “real data” masks: ignore pure zeros
        maskG_data  = maskG  & (g  ~= 0);
        maskNG_data = maskNG & (ng ~= 0);
    
        vals = [g(maskG_data); ng(maskNG_data)];
        vals = vals(isfinite(vals));
        if ~isempty(vals)
            hasData = true;
            prY = prctile(vals,vi_pr);
            global_min = min(global_min, prY(1));
            global_max = max(global_max, prY(2));
        end
     end

    % If no valid data anywhere, skip making the violin axis
    if ~hasData || ~isfinite(global_min) || ~isfinite(global_max)
        continue;
    end

    % Dynamic y-limits with small padding
    ymin  = -.13;
    ymax  = .25;

    % --- Now create the axis and plot violins (second pass) ---
    dist_pos = [map_pos(1), ...
                ybot_n - max(ytop_n - ybot_n, 1e-9) - 0.05, ...
                zonal_w + map_width + gap_h, ...
                max(ytop_n - ybot_n, 1e-9)-.1];
    axd = axes('Position', dist_pos); 
    hold(axd,'on');
    set(axd, 'TickDir','in', ...          % ticks go inward (cleaner for violins)
         'TickLength',[0.007 0.007], ...
         'LineWidth',1.0);

    nbins    = 250;
    centers  = 1:(numMods+1);      % +1 for Ensemble
    maxBarW  = 0.50;
    y_edges   = linspace(ymin, ymax, nbins+1);
    y_centers = 0.5*(y_edges(1:end-1) + y_edges(2:end));
    y = y_centers(:);

    % background halves (blue up / red down)
    xlim(axd, [.25, numMods+1.75]);
    ylim(axd, [ymin, ymax]);
    patch(axd,[0 numMods+1.75 numMods+1.75 0], [0 0 ymax ymax], [0.3 0.5 0.85], ...
          'EdgeColor','none','FaceAlpha',0.06);
    patch(axd,[0 numMods+1.75 numMods+1.75 0], [ymin ymin 0 0], [0.85 0.3 0.3], ...
          'EdgeColor','none','FaceAlpha',0.06);
    yline(axd,0,'k-');

    set(axd,'Layer','top','FontSize',15,'TickLabelInterpreter','latex');
    set(axd,'XTick',centers,'XTickLabel',[]);
    ylabel(axd, baseName,'FontSize', 20, 'Interpreter','latex');


    for ii = 1:numMods
        text(axd, centers(ii), ymin - 0.04*(ymax-ymin), model_names{ii}, ...
            'HorizontalAlignment','center', 'VerticalAlignment','top', ...
            'FontSize',16, 'Color', cmap_mod(ii,:), 'Interpreter','latex');
    end
    % Ensemble label (black)
    text(axd, centers(numMods+1), ymin - 0.04*(ymax-ymin), '\textbf{Ensemble}', ...
        'HorizontalAlignment','center', 'VerticalAlignment','top', ...
        'FontSize',16, 'Color','k', 'Interpreter','latex');

    xline(axd, centers, ':', 'Color', [0 0 0 0.25]);

    % -------- per-model biomass-weighted KDE violins (second pass) --------
    g_all  = [];  w_all  = [];
    ng_all = [];  wn_all = [];

    for ii = 1:numMods
        S = load(mods(mod_order(model_id(ii))).name, ...
             'CoreVar_Inv_rw', 'Bio_Weights');
        V = S.CoreVar_Inv_rw;
        A = S.Bio_Weights;
    
        g  = squeeze(V(var_id ,:,:));      % full-field distribution
        ng = squeeze(V(var_id2 ,:,:));      % TD-only subset
        
        BU = squeeze(V(1,:,:));            
        TD = squeeze(V(6,:,:));           
        TDmask = sign(BU) .* sign(TD) == -1;    % opposite signs = top-down
        
        % basic validity
        
        maskG  = isfinite(g)  & isfinite(A) & (A > 0);   
        if TDmask2
        maskNG = isfinite(ng) & isfinite(A) & (A > 0) & TDmask;   
        else
        maskNG = isfinite(ng) & isfinite(A) & (A > 0) ;  
        end

        % “real data” = nonzero + valid
        maskG_data  = maskG  & (g  ~= 0);
        maskNG_data = maskNG & (ng ~= 0);
                
        % accumulate for ENSEMBLE (only nonzero data)
        if any(maskG_data(:))
            g_all = [g_all; g(maskG_data)];   %#ok<AGROW>
            w_all = [w_all; A(maskG_data)];   %#ok<AGROW>
        end
        if any(maskNG_data(:))
            ng_all = [ng_all; ng(maskNG_data)]; %#ok<AGROW>
            wn_all = [wn_all; A(maskNG_data)];  %#ok<AGROW>
        end
    
        % nothing to draw for this model?
        if ~any(maskG_data(:)) && ~any(maskNG_data(:))
            continue;
        end


        % ---- biomass-weighted KDEs ----
            fG  = zeros(size(y));
            fNG = zeros(size(y));
            if any(maskG_data(:))
                fG = ksdensity(g(maskG_data),  y, 'Weights', A(maskG_data), ...
                               'Support','unbounded', 'Function','pdf');
            end
            if any(maskNG_data(:))
                fNG = ksdensity(ng(maskNG_data), y, 'Weights', A(maskNG_data), ...
                               'Support','unbounded', 'Function','pdf');
            end


        col = cmap_mod(model_id(ii),:);
        mG   = max(fG);  if ~(isfinite(mG)  && mG>0),  mG  = 1; end
        mNG  = max(fNG); if ~(isfinite(mNG) && mNG>0), mNG = 1; end

        wG  = -(fG  / mG)  * maxBarW;   % left (≤0)
        wNG =  (fNG / mNG) * maxBarW;   % right (≥0)

        % LEFT (g)
        if any(maskG_data(:)) && any(isfinite(wG))
            xL = centers(ii) + wG;                 
            xC = centers(ii) + zeros(size(y));     
            patch('Parent',axd, ...
                  'XData',[xL; flipud(xC)], 'YData',[y;  flipud(y)], ...
                  'FaceColor',col, 'FaceAlpha',0.7, 'EdgeColor','none');
        end

        % RIGHT (ng)
        if any(maskNG_data(:)) && any(isfinite(wNG))
            xC = centers(ii) + zeros(size(y));   
            xR = centers(ii) + wNG;              
            patch('Parent',axd, ...
                  'XData',[xC; flipud(xR)], 'YData',[y;  flipud(y)], ...
                  'FaceColor',col, 'FaceAlpha',0.1, ...
                  'EdgeColor',cmap_mod(ii,:),'linewidth',.5);
        end

        % ---- biomass-weighted medians & IQR ----
        g_med = NaN; g_25 = NaN; g_75 = NaN;
        ng_med = NaN; ng_25 = NaN; ng_75 = NaN;
        if any(maskG_data(:))
            [g_med, g_25, g_75] = local_wquant3(g(maskG), A(maskG));
        end
        if any(maskNG_data(:))
            [ng_med, ng_25, ng_75] = local_wquant3(ng(maskNG), A(maskNG));
        end

        cap  = 0.08;   lw = 1.2;

        % LEFT caps/median (g)
        if isfinite(g_med)
            wGm = interp1(y, wG, g_med, 'linear','extrap');  
            if isfinite(wGm)
                xIQR = centers(ii) + 0.5*wGm;
                h1 = plot(axd, [xIQR xIQR], [g_25 g_75], 'Color','k','LineWidth', lw);
                h2 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [g_25 g_25], 'Color','k','LineWidth',0.7);
                h3 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [g_75 g_75], 'Color','k','LineWidth',0.7);
                h4 = scatter(axd, xIQR, g_med, 50, 'w', 'filled', ...
                             'Marker','d', 'MarkerEdgeColor','k', 'LineWidth',1);
                uistack([h1 h2 h3 h4],'top');
            end
        end

        % RIGHT caps/median (ng)
        if isfinite(ng_med)
            wNGm = interp1(y, wNG, ng_med, 'linear','extrap');
            if isfinite(wNGm)
                xIQR = centers(ii) + 0.5*wNGm;
                h1 = plot(axd, [xIQR xIQR], [ng_25 ng_75], 'Color','k','LineWidth', lw);
                h2 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [ng_25 ng_25], 'Color','k','LineWidth',0.7);
                h3 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [ng_75 ng_75], 'Color','k','LineWidth',0.7);
                h4 = scatter(axd, xIQR, ng_med, 50, 'k', 'filled', ...
                             'Marker','d', 'MarkerEdgeColor','k', 'LineWidth',1);
                uistack([h1 h2 h3 h4],'top');
            end
        end

        % ---- % positive / negative annotations for each model ----
        fsn = 16; 
        dy  = .005*(ymax - ymin); 
        x_off = 0.45;

        if any(maskG_data(:))
            xg = g(maskG); wg = A(maskG);
            f_pos_g = weighted_frac(xg, wg, xg > 0);
            f_neg_g = weighted_frac(xg, wg, xg < 0);
            xL = centers(ii) - x_off*maxBarW;
            if round(100*f_pos_g) > 0
                text(axd, xL, ymax - dy, sprintf('%d\\%%', round(100*f_pos_g)), ...
                    'HorizontalAlignment','center','VerticalAlignment','top', ...
                    'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
            end
            if round(100*f_neg_g) > 0
                text(axd, xL, ymin + dy, sprintf('%d\\%%', round(100*f_neg_g)), ...
                    'HorizontalAlignment','center','VerticalAlignment','bottom', ...
                    'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
            end
        end

        if any(maskNG_data(:))
            xng = ng(maskNG); wng = A(maskNG);
            f_pos_ng = weighted_frac(xng, wng, xng > 0);
            f_neg_ng = weighted_frac(xng, wng, xng < 0);
            xR = centers(ii) + x_off*maxBarW;
            if round(100*f_pos_ng) > 0
                text(axd, xR, ymax - dy, sprintf('%d\\%%', round(100*f_pos_ng)), ...
                    'HorizontalAlignment','center','VerticalAlignment','top', ...
                    'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
            end
            if round(100*f_neg_ng) > 0
                text(axd, xR, ymin + dy, sprintf('%d\\%%', round(100*f_neg_ng)), ...
                    'HorizontalAlignment','center','VerticalAlignment','bottom', ...
                    'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
            end
        end

    end % loop over models

    % ---------- ENSEMBLE VIOLIN (at centers(numMods+1)) ----------
    cx          = centers(numMods+1);
    maxBarW_ens = maxBarW;

    % Left (g_all) in solid black
    if ~isempty(g_all) && any(w_all>0)
        [fg,~] = ksdensity(g_all, y, 'Weights', w_all, 'NumPoints', numel(y), ...
                           'Support','unbounded','Function','pdf');
        if any(fg>0), fg = fg / max(fg); else, fg = zeros(size(fg)); end
        wG = -maxBarW_ens * fg;
        patch(axd, [cx + wG; cx*ones(size(y))], [y; flipud(y)], ...
              'k', 'EdgeColor','none', 'FaceAlpha', 0.70);
        plot(axd, cx + wG, y, 'Color',[0 0 0 0.30], 'LineWidth',0.25);

        [g_medE, g_q25E, g_q75E] = local_wquant3(g_all, w_all);
        cap = .08; lw = 1.1; colk = [0 0 0]; xIQR = cx - .5*maxBarW_ens;
        plot(axd, [xIQR-cap/2, xIQR+cap/2],[g_q25E g_q25E],'Color',colk,'LineWidth',lw);
        if ymax>g_q75E
            plot(axd, [xIQR-cap/2, xIQR+cap/2],[g_q75E g_q75E],'Color',colk,'LineWidth',lw);
        end
        plot(axd, [xIQR, xIQR],[g_q25E min(ymax,g_q75E)],'Color',colk,'LineWidth',lw);
        scatter(axd, xIQR, g_medE, 55, 'w', 'filled', ...
                'Marker','d', 'MarkerEdgeColor','k', 'LineWidth',1.1);
    end

    % Right (ng_all) in translucent black
    if ~isempty(ng_all) && any(wn_all>0)
        [fn,~] = ksdensity(ng_all, y, 'Weights', wn_all, 'NumPoints', numel(y), ...
                           'Support','unbounded','Function','pdf');
        if any(fn>0), fn = fn / max(fn); else, fn = zeros(size(fn)); end
        wN =  maxBarW_ens * fn;
        patch(axd, [cx + wN; cx*ones(size(y))], [y; flipud(y)], ...
              'k', 'EdgeColor','k', 'FaceAlpha', 0.10, 'LineWidth', 0.5);
        plot(axd, cx + wN, y, 'Color',[0 0 0 0.25], 'LineWidth',0.4);

        [ng_medE, ng_q25E, ng_q75E] = local_wquant3(ng_all, wn_all);
        cap = .08; lw = 1.1; colk = [0 0 0]; xIQR = cx + .5*maxBarW_ens;
        plot(axd, [xIQR-cap/2, xIQR+cap/2],[ng_q25E ng_q25E],'Color',colk,'LineWidth',lw);
        if ymax>ng_q75E
            plot(axd, [xIQR-cap/2, xIQR+cap/2],[ng_q75E ng_q75E],'Color',colk,'LineWidth',lw);
        end
        plot(axd, [xIQR, xIQR],[ng_q25E min(ymax,ng_q75E)],'Color',colk,'LineWidth',lw);
        scatter(axd, xIQR, ng_medE, 55, 'k', 'filled', ...
                'Marker','d', 'MarkerEdgeColor','k', 'LineWidth',1.1);
    end

    % Ensemble % positive / negative annotations
    dy = .005*(ymax - ymin); x_off = 0.45;
    if ~isempty(g_all)
        f_pos_gE = weighted_frac(g_all, w_all, g_all > 0);
        f_neg_gE = weighted_frac(g_all, w_all, g_all < 0);
        xL = cx - x_off*maxBarW_ens;
        if round(100*f_pos_gE) > 0
            text(axd, xL, ymax - dy, sprintf('%d\\%%', round(100*f_pos_gE)), ...
                'HorizontalAlignment','center','VerticalAlignment','top', ...
                'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
        end
        if round(100*f_neg_gE) > 0
            text(axd, xL, ymin + dy, sprintf('%d\\%%', round(100*f_neg_gE)), ...
                'HorizontalAlignment','center','VerticalAlignment','bottom', ...
                'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
        end
    end
    if ~isempty(ng_all)
        f_pos_ngE = weighted_frac(ng_all, wn_all, ng_all > 0);
        f_neg_ngE = weighted_frac(ng_all, wn_all, ng_all < 0);
        xR = cx + x_off*maxBarW_ens;
        if round(100*f_pos_ngE) > 0
            text(axd, xR, ymax - dy, sprintf('%d\\%%', round(100*f_pos_ngE)), ...
                'HorizontalAlignment','center','VerticalAlignment','top', ...
                'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
        end
        if round(100*f_neg_ngE) > 0
            text(axd, xR, ymin + dy, sprintf('%d\\%%', round(100*f_neg_ngE)), ...
                'HorizontalAlignment','center','VerticalAlignment','bottom', ...
                'Interpreter','latex','Color','k','FontSize',fsn,'Clipping','on');
        end
    end


      % Make a full box and inward ticks
        set(axd, 'Box','on', ...
                 'TickDir','in', ...
                 'TickLength',[0.010 0.010], ...
                 'LineWidth',1.1);
        
        % Choose y-ticks explicitly (nice, avoids weird auto ticks)
        yt = -.1:0.05:.2;   % <-- adjust if you want (must sit within [ymin ymax])
        set(axd,'YTick',yt);
        
        % LEFT axis (grazing)
        yyaxis(axd,'left');
        set(axd,'YColor','k');   % make sure left spine stays black
        ylabel(axd,'Grazing Pressure, $\Delta l_g$ (d$^{-1}$)','Interpreter','latex','fontsize',20);   % Figure 4d,e
        %ylabel(axd,'Loss Rates, $\Delta l$ (d$^{-1}$)','Interpreter','latex','fontsize',20);           % Supplementary Figure 5b

        % RIGHT axis (non-grazing): same ticks/limits, ticks inward
        yyaxis(axd,'right');
        set(axd,'YColor','k', 'YTick',yt, 'YLim',[ymin ymax]);   % match left exactly
        ylabel(axd,'Non-grazing Pressure, $\Delta l_{ng}$ (d$^{-1}$)','Interpreter','latex','fontsize',20);  % Figure 4de,e
        %ylabel(axd,'Division Rate, $\Delta \mu$ (d$^{-1}$)','Interpreter','latex','fontsize',20);           % Supplementary Figure 5b
       
        
        % Go back to left so any later plotting uses left axis
        yyaxis(axd,'left');


end % Violin_dist


    




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  Helper Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xn, yn] = axData2norm(ax, x, y)
% Convert data coords (x,y) in axes AX to figure normalized coords
    pos = get(ax,'Position');     % [x y w h] in normalized figure units
    xl  = get(ax,'XLim');  yl = get(ax,'YLim');
    xn  = pos(1) + (x - xl(1)) / (xl(2)-xl(1)) * pos(3);
    yn  = pos(2) + (y - yl(1)) / (yl(2)-yl(1)) * pos(4);
end

function counts = histcounts_weighted(x, edges, w)
    % x: values (vector), edges: bin edges, w: weights (same size as x)
    % returns 1 x (numBins) row vector of weighted sums per bin
    idx   = discretize(x, edges);                 % bin index per sample
    valid = ~isnan(idx) & isfinite(w);
    counts_col = accumarray(idx(valid), w(valid), [numel(edges)-1, 1], @sum, 0);
    counts = counts_col.';                        % make it a row vector
end

function q = local_wquant(x,w,p)
    % returns NaN if not enough information
    ok = isfinite(x) & isfinite(w) & (w > 0);
    x = x(ok); w = w(ok);
    if numel(x) == 0 || sum(w) <= 0
        q = NaN; return
    end
    % sort by x
    [xs, ix] = sort(x(:));
    ws = w(:); ws = ws(ix);

    % cumulative weights in [0,1]
    cw = cumsum(ws);
    cw = cw ./ cw(end);

    % collapse duplicate cumulative weights before interp1
    [cu, iu] = unique(cw, 'stable');
    xsu = xs(iu);

    % if still only one point, return it
    if numel(cu) == 1
        q = xsu(1); return
    end

    % interpolate the inverse CDF
    q = interp1(cu, xsu, p, 'linear', 'extrap');
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
    % IQR
    q25 = interp1(cu, xsu, 0.25, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.75, 'linear','extrap');
    % pm 10%
    % q25 = interp1(cu, xsu, 0.40, 'linear','extrap');
    % q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    % q75 = interp1(cu, xsu, 0.60, 'linear','extrap');
end


function f = weighted_frac(x, w, cond)
    % biomass/area-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