%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure 4a. Driver Magnitude
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; close all
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))
cd /Users/twrohr/Desktop/data.nosync/SSP585
          load('Longhurst360.mat');

for startup =[1]
% ---------------- Model names / ordering / colors (kept consistent) ------
 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];  % provided order → maps to directory listing
model_id    = 1:11;
b1 = -25; b2 =25;

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];

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

com_col = [wedge_colors(4,:); lighter_colors(4,:);... 
           wedge_colors(5,:); lighter_colors(5,:);... 
           wedge_colors(1,:); lighter_colors(1,:);... 
           wedge_colors(3,:); lighter_colors(3,:);... 
           wedge_colors(2,:); lighter_colors(2,:);
           0.2  0.2  0.2]; 
end  


%% ======================== Build Loss_Driver ========================
clearvars Loss_Drivers
latims = [-60, 60]; % -60 to 60
lat     =  latims(1):latims(2);
latlim  =  90+latims(1) : 90+latims(2);   % index into 180-lat zonal arrays
varid = [30,30+5,31,31+5,27,27+5,29,29+5,28,28+5,47 ];
 decomp_symbols = {'$\Delta l_{g,\Sigma P}$', '$\Delta l_{ng,\Sigma P}$', ...
                    '$\Delta l_{g,\Sigma Z}$','$\Delta l_{ng,\Sigma Z}$', ...
                    '$\Delta l_{g,\ T}$',  '$\Delta l_{ng,T}$', ...
                    '$\Delta l_{g,CC}$', '$\Delta l_{ng, CC}$', ...
                    '$\Delta l_{g, VD}$','$\Delta l_{ng, VD}$', ...
                    '$\Delta l_{g, Int}$'};
  decomp_symbols2 = {'$l_{g,\Sigma P}$', '$l_{ng,\Sigma P}$', ...
                    '$l_{g,\Sigma Z}$','$l_{ng,\Sigma Z}$', ...
                    '$l_{g,T}$',  '$l_{ng,T}$', ...
                    '$l_{g,CC}$', '$l_{ng,CC}$', ...
                    '$l_{g,VD}$','$l_{mg,VD}$', ...
                    '$l_{Int}$'};
 decomp_symbols3 = {'$\Delta l_{\Sigma P}$', '$\Delta l_{\Sigma Z}$', '$\Delta l_{T}$',...
                    '$\Delta  l_{CC}$','$\Delta l_{VD}$', '$\Delta l_{Int}$'};
%%%% Pick Region: 1. Polar North 2. Polar South 3. Trades Atlantic 4. Trades Indian 5. Trades Pacific 6. North Atlantic 7. Westerlies Pacific 8. Westerlies Southern; 10: SO (West S + Polar S) 
                     % 11. North Pacific
 rm = false; pcth =.5; region_index = 1000;
 %rm = true; region_index =6; latlim = 1:180;

for i = 1:length(varid)
    fprintf([num2str(varid(i)), '\n'])
          % Initialize Ensemble
          Emag_all = [];     EFMag_all  = [];    Ew_all = [];
          Emag_ud  = [];     EFMag_ud  = [];     Ew_ud = [];
          Emag_ui = [];      EFMag_ui  = [];     Ew_ui = [];
          Emag_bu_pd  = [];  EFMag_bu_pd  = [];  Ew_bu_pd = [];
          Emag_td_pd  = [];  EFMag_td_pd  = [];  Ew_td_pd = [];
          Emag_bu_pi  = [];  EFMag_bu_pi  = [];  Ew_bu_pi = [];
          Emag_td_pi  = [];  EFMag_td_pi  = [];  Ew_td_pi = [];
          Emag_bu  = [];     EFMag_bu  = [];     Ew_bu = [];
          Emag_td  = [];     EFMag_td  = [];     Ew_td= [];
    
for s = 1:12
    % -------- Single model  ----------
    if s<12
      cd /Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2
      mods = dir('*.mat');
      fname = mods(mod_order(model_id(s))).name;
      cd /Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2
      load(fname, 'CoreVar_Inv_rw', 'Bio_Weights');
          
      % Masks and Weights
          u      = squeeze(CoreVar_Inv_rw(6,latlim,:));                      u = u(:);
          P      = squeeze(CoreVar_Inv_rw(1,latlim,:));                      P = P(:); 
          w      = squeeze(Bio_Weights(latlim,:));                           w =   w(:);
 
      % Loss Decompositon   
          % Metric 1 
          mag    = squeeze(abs(CoreVar_Inv_rw(varid(i),latlim,:)));               mag = mag(:);
          % Metric 2 
          fmag    = squeeze((CoreVar_Inv_rw(varid(i),latlim,:)));                 fmag = fmag(:);
          %FMag = fmag;
           FMag   = squeeze(sum(abs(CoreVar_Inv_rw(varid,latlim,:)),'omitnan'));   FMag = fmag./FMag(:);
          % Apply Longhurst mask: 
                if region_index ==10
                    region_mask = double(Long_Prov == 8 | Long_Prov == 2);
                    %region_mask(  Lat > -53 ) = logical(0);
                region_mask = logical(region_mask);
                elseif region_index ==11
                    region_mask = double(Long_Prov == 7  );
                    region_mask = logical(region_mask);
                    region_mask( Lat < 10 |  Lat > 68 ) = logical(0); % Remove SH% Remove Pacific !!Latmin bound 40 for CESM; 20 normalls
                elseif region_index ==6
                    region_mask = double(Long_Prov == 1 | Long_Prov == 6);
                    region_mask = logical(region_mask); 
                    region_mask((Lon>0 & Lon<250 ) |   Lat>70 |  Lat<40 ) =logical(0); 
                else
                region_mask = double(Long_Prov == region_index);
                region_mask = logical(region_mask);
                end
                region_mask =region_mask(:);
          % Mask 
          if rm; mask   = isfinite(mag) & isfinite(P) & region_mask;
          else;  mask   = isfinite(mag) & isfinite(P) ; end
          
          % All values (no TD or BU filter)
          mag_all = mag(mask);        FMag_all  = FMag(mask);        w_all  = w(mask); 
          mag_ui  = mag(mask & u>0 ); FMag_ui   = FMag(mask & u>0 ); w_ui   = w(mask & u>0 ); 
          mag_ud = mag(mask & u<0 );  FMag_ud   = FMag(mask& u<0 );  w_ud   = w(mask & u<0 ); 

          if  varid(i) == 36 || ... 
             (any(varid(i) == [27]) && any(s==[1,2,4,6,7]) )   || ... % [1,2,3,4,6,7] not on grazing 
             (any(varid(i) == [32]) && any(s==[1,2,4,6,7]) ) || ... % [1,2,4,6,7,8,9] not on non grazing
             (any(varid(i) == [29,34]) && any(s==[1,2,3,4])); 
             Loss_Drivers(s,      :,     [1,6,7] ,       i,      :) = nan;
          else
             % Add To Esemble
              Emag_all =  [Emag_all; mag_all];   EFMag_all =  [EFMag_all; FMag_all]; Ew_all =  [Ew_all; w_all./nansum(w_all)];
              Emag_ui =  [Emag_ui; mag_ui];     EFMag_ui =  [EFMag_ui; FMag_ui]; Ew_ui =  [Ew_ui; w_ui./nansum(w_ui)];
              Emag_ud =  [Emag_ud; mag_ud];     EFMag_ud =  [EFMag_ud; FMag_ud]; Ew_ud =  [Ew_ud; w_ud./nansum(w_ud)];

              % Add To individual Models
                  % All Points
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_all, w_all);
                  mn = wmean(mag_all, w_all);   stdv = std(mag_all, w_all);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        1,       i,     :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_all, w_all);
                  mn = wmean(FMag_all, w_all);   stdv = std(FMag_all, w_all);
                           % model | metric | control | driver | stat
                  Loss_Drivers(s,      2,        1,       i,      :)  = [med, q25, q75, mn, stdv];

                  % u increasing
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_ui, w_ui);
                  mn = wmean(mag_ui, w_ui);   stdv = std(mag_ui, w_ui);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        6,       i,     :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_ui, w_ui);
                  mn = wmean(FMag_ui, w_ui);   stdv = std(FMag_ui, w_ui);
                           % model | metric | control | driver | stat
                  Loss_Drivers(s,      2,        6,       i,      :)  = [med, q25, q75, mn, stdv];

                  % u decreasing
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_ud, w_ud);
                  mn = wmean(mag_ud, w_ud);   stdv = std(mag_ud, w_ud);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        7,       i,     :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_ud, w_ud);
                  mn = wmean(FMag_ud, w_ud);   stdv = std(FMag_ud, w_ud);
                           % model | metric | control | driver | stat
                  Loss_Drivers(s,      2,        7,       i,      :)  = [med, q25, q75, mn, stdv];
          end


          %%%%% ----- u increasing -----
          mag_td = mag(mask & u.*P<0 & P<0);    FMag_td  = FMag(mask & u.*P<0 & P<0);  w_td  = w(mask & u.*P<0 & P<0); 
          mag_bu = mag(mask & u.*P>0 & P>0);    FMag_bu  = FMag(mask & u.*P>0 & P>0);  w_bu  = w(mask & u.*P>0 & P>0);
          if  100.*length(mag_td)/length(mag_all) >= pcth && 100.*length(mag_bu)/length(mag_all) >= pcth; proc = true; else; proc = false; fprintf(['Not enough (u in) mod=',num2str(s),' \n']); end 
          if ~proc; Loss_Drivers(s,      :,        2:3 ,       i,      :) = nan; else

              % TD values & P Dec
              if varid(i) == 36 || ... 
                 (any(varid(i) == [27]) && any(s==[1,2,4,6,7]) )   || ... % [1,2,3,4,6,7] not on grazing 
                 (any(varid(i) == [32]) && any(s==[1,2,4,6,7]) )   || ... % [1,2,4,6,7,8,9] no on non grazing
                 (any(varid(i) == [29,34]) && any(s==[1,2,3,4])); 
                 Loss_Drivers(s,      :,        2:3 ,       i,      :) = nan;
              else

              % Add To ensemble    
              Emag_td_pd =  [Emag_td_pd; mag_td];   EFMag_td_pd =  [EFMag_td_pd; FMag_td]; Ew_td_pd =  [Ew_td_pd; w_td./nansum(w_td)];
              Emag_bu_pi =  [Emag_bu_pi; mag_bu];   EFMag_bu_pi =  [EFMag_bu_pi; FMag_bu]; Ew_bu_pi =  [Ew_bu_pi; w_bu./nansum(w_bu)];
              
              % Add To individual models 
                  % TD values & P Dec    
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_td, w_td);
                  mn = wmean(mag_td, w_td);   stdv = std(mag_td, w_td);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        2,       i,      :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_td, w_td);
                  mn = wmean(FMag_td, w_td);   stdv = std(FMag_td, w_td);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      2,        2,       i,     :)  = [med, q25, q75, mn, stdv];
    
                  % BU values & P Inc
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_bu, w_bu);
                  mn = wmean(mag_bu, w_bu);   stdv = std(mag_bu, w_bu);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        3,       i,     :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_bu, w_bu);
                  mn = wmean(FMag_bu, w_bu);   stdv = std(FMag_bu, w_bu);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      2,        3,       i,      :)  = [med, q25, q75, mn, stdv];
              end

                  

          end


          %%%%% ----- u decreasing -----
          mag_td = mag(mask & u.*P<0 & P>0);  FMag_td  = FMag(mask & u.*P<0 & P>0);   w_td  = w(mask & u.*P<0 & P>0); 
          mag_bu = mag(mask & u.*P>0 & P<0);  FMag_bu  = FMag(mask & u.*P>0 & P<0);   w_bu  = w(mask & u.*P>0 & P<0);
          if  100.*length(mag_td)/length(mag_all) >=pcth & 100.*length(mag_bu)/length(mag_all) >=pcth; proc = true; else; proc = false;  fprintf(['Not enough (u dec) mod=',num2str(s),' \n']); end 
          if ~proc; Loss_Drivers(s,      :,        4:5,       i,      :) = nan; else
    
              % TD values & P Inc
              if varid(i) == 36 || ... 
                 (any(varid(i) == [27]) && any(s==[1,2,4,6,7]) )   || ... % [1,2,3,4,6,7] not on grazing 
                 (any(varid(i) == [32]) && any(s==[1,2,4,6,7]) ) || ... % [1,2,4,6,7,8,9] no on non grazing
                 (any(varid(i) == [29,34]) && any(s==[1,2,3,4])); 
                 Loss_Drivers(s,      :,        4:5,       i,      :) = nan;
              else

              % Add to Ensemble
              Emag_td_pi =  [Emag_td_pi; mag_td];   EFMag_td_pi =  [EFMag_td_pi; FMag_td]; Ew_td_pi =  [Ew_td_pi; w_td./nansum(w_td)];
              Emag_bu_pd =  [Emag_bu_pd; mag_bu];   EFMag_bu_pd =  [EFMag_bu_pd; FMag_bu]; Ew_bu_pd =  [Ew_bu_pd; w_bu./nansum(w_bu)];
              
              % Add to Individual Models
              % TD values & P Inc
                 % Raw
                  [med, q25, q75] = local_wquant3(mag_td, w_td);
                  mn = wmean(mag_td, w_td);   stdv = std(mag_td, w_td);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        4,       i,      :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_td, w_td);
                  mn = wmean(FMag_td, w_td);   stdv = std(FMag_td, w_td);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      2,        4,       i,     :)  = [med, q25, q75, mn, stdv];
                 
              % BU values & P Dec
                  % Raw
                  [med, q25, q75] = local_wquant3(mag_bu, w_bu);
                  mn = wmean(mag_bu, w_bu);   stdv = std(mag_bu, w_bu);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      1,        5,       i,     :)  = [med, q25, q75, mn, stdv];
                  % Magnitude Share
                  [med, q25, q75] = local_wquant3(FMag_bu, w_bu);
                  mn = wmean(FMag_bu, w_bu);   stdv = std(FMag_bu, w_bu);
                           % model | metric | control | driver |stat
                  Loss_Drivers(s,      2,        5,       i,      :)  = [med, q25, q75, mn, stdv];
              end
                 
          end 


    else % Ensemble
         if varid(i) == 36 
                 Loss_Drivers(12,      :,        :,       i,      :) = nan;
         else
         % All values
         % mask   = isfinite(Emag_all) & isfinite(EFMag_all) ;
         


              % Raw
              [med, q25, q75] = local_wquant3(Emag_all, Ew_all);
              mn = wmean(Emag_all, Ew_all);   stdv = std(Emag_all, Ew_all);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        1,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_all, Ew_all);
              mn = wmean(EFMag_all, Ew_all);   stdv = std(EFMag_all, Ew_all);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        1,       i,     :)  = [med, q25, q75, mn, stdv];

          % TD values & P Dec
              % Raw
              [med, q25, q75] = local_wquant3(Emag_td_pd, Ew_td_pd);
              mn = wmean(Emag_td_pd, Ew_td_pd);   stdv = std(Emag_td_pd, Ew_td_pd);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        2,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_td_pd, Ew_td_pd);
              mn = wmean(EFMag_td_pd, Ew_td_pd);   stdv = std(EFMag_td_pd, Ew_td_pd);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        2,       i,     :)  = [med, q25, q75, mn, stdv];
  
          % BU values & P Inc
              % Raw
              [med, q25, q75] = local_wquant3(Emag_bu_pi, Ew_bu_pi);
              mn = wmean(Emag_bu_pi, Ew_bu_pi);   stdv = std(Emag_bu_pi, Ew_bu_pi);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        3,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_bu_pi, Ew_bu_pi);
              mn = wmean(EFMag_bu_pi, Ew_bu_pi);   stdv = std(EFMag_bu_pi, Ew_bu_pi);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        3,       i,     :)  = [med, q25, q75, mn, stdv];

          % TD values & P Inc
              % Raw
              [med, q25, q75] = local_wquant3(Emag_td_pi, Ew_td_pi);
              mn = wmean(Emag_td_pi, Ew_td_pi);   stdv = std(Emag_td_pi, Ew_td_pi);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        4,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_td_pi, Ew_td_pi);
              mn = wmean(EFMag_td_pi, Ew_td_pi);   stdv = std(EFMag_td_pi, Ew_td_pi);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        4,       i,     :)  = [med, q25, q75, mn, stdv];

         % BU values & P Dec
              % Raw
              [med, q25, q75] = local_wquant3(Emag_bu_pd, Ew_bu_pd);
              mn = wmean(Emag_bu_pd, Ew_bu_pd);   stdv = std(Emag_bu_pd, Ew_bu_pd);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        5,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_bu_pd, Ew_bu_pd);
              mn = wmean(EFMag_bu_pd, Ew_bu_pd);   stdv = std(EFMag_bu_pd, Ew_bu_pd);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        5,       i,     :)  = [med, q25, q75, mn, stdv];

         % u Inc
              % Raw
              [med, q25, q75] = local_wquant3(Emag_ui, Ew_ui);
              mn = wmean(Emag_ui, Ew_ui);   stdv = std(Emag_ui, Ew_ui);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        6,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_ui, Ew_ui);
              mn = wmean(EFMag_ui, Ew_ui);   stdv = std(EFMag_ui, Ew_ui);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        6,       i,     :)  = [med, q25, q75, mn, stdv];

         % u Dec
              % Raw
              [med, q25, q75] = local_wquant3(Emag_ud, Ew_ud);
              mn = wmean(Emag_ud, Ew_ud);   stdv = std(Emag_ud, Ew_ud);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      1,        7,       i,     :)  = [med, q25, q75, mn, stdv];
              % Magnitude Share
              [med, q25, q75] = local_wquant3(EFMag_ud, Ew_ud);
              mn = wmean(EFMag_ud, Ew_ud);   stdv = std(EFMag_ud, Ew_ud);
                       % model | metric | control | driver | stat
              Loss_Drivers(12,      2,        7,       i,     :)  = [med, q25, q75, mn, stdv];

         end
       
    end

end
end




%% ====== Absolute Change + TD–BU anomaly — explicit axes (no tiledlayout) ======
for plot = 1
statMean  = 4;                 % mean
statStd   = 5;                 % std
nmod      = 11;                % # individual models (exclude ensemble)
metricIdx = 1;                 % panel 1 metric
lw = 1;  sz = 100;               % line width + dot size
stdsc = 4;                      % (kept; set to 1 = no shrinking)

% Regions (for anomaly panels 2 & 3)
% 1. All 2. TD PD + uI 3. BU PI + uI 4. TD PI +uD 5. BU PD + uD 6. uI 7. uD 
regTD =   [2,4];         % Anomaly Regions [2,4];  Standard - Anomaly  target
regBU =   [3,5];        % Baseline Region [3,5];  Standard  -  Baslein 
x  = [1:10,11.5];   % panel 1
x6 = 1:6;    % paired labels panels 2–3

% Interpreters
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaultTextInterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');

% ---------- Figure & manual axes layout ----------
clf; 
figure('Color','w','Units','normalized','Position',[0.04 0.10 1 1]);

% Three panels left-to-right
W1  = 0.8;          % width per panel
H1 = .5;   % .4 for figure 4
gh = .08;
gw = .035;
X0 = 0.06;          % left margin
Y0 = 0.04;

pos1 = [X0,             Y0+.1+gh,         W1,  H1];                % Panel 1 (Mean -Absolute change)
pos3 = [X0 ,                  Y0, (W1-gw)./2,  H1];                   % Panel 2 (TD–BU improving)
pos4 = [X0 + (W1-gw)./2+gw,   Y0, (W1-gw)./2,  H1];                   % Panel 3 (TD–BU worsening)

% ===== Panel 1: Absolute change (metricIdx = 1) =====
G_mean   = squeeze(Loss_Drivers(1:nmod, metricIdx, 1, :, statMean));  % [nmod x 11]
G_std    = squeeze(Loss_Drivers(1:nmod, metricIdx, 1, :, statStd));
G_mean_E = squeeze(Loss_Drivers(12,     metricIdx, 1, :, statMean));  % [1 x 11]
G_std_E  = squeeze(Loss_Drivers(12,     metricIdx, 1, :, statStd));

Ymod_mean = G_mean;
Ymod_std  = G_std ;
Yens_mean = G_mean_E;
Yens_std  = G_std_E ;

%%%%%%%%%%%%%%
% Panel 1 - Mean Absolute Change
spread  = .32;                                   % half-width of model band around each driver
offsets = linspace(-.5*spread, .5*spread, nmod); % fixed left->right order

ax1 = axes('Position', pos1); hold(ax1,'on'); box(ax1,'off');

% Ensemble bars
b1 = bar(ax1, x, Yens_mean, 0.72, 'FaceColor','flat', 'EdgeColor',[0 0 0], 'LineWidth',1.2);
b1.CData = com_col;

% Per-model dots ± std at fixed offsets (left->right)
for j = 1:nmod
    c  = cmap_mod(j,:);
    cx = x + offsets(j);
    scatter(ax1, cx, Ymod_mean(j,:), sz, c, 'filled', 'MarkerEdgeColor','k', 'LineWidth',0.2);
    errorbar(ax1, cx, Ymod_mean(j,:), Ymod_std(j,:)./stdsc, Ymod_std(j,:)./stdsc, ...
        'Color', c, 'LineStyle','none','LineWidth',lw,'CapSize',6);
end

% X/Y limits and ticks
xlim(ax1, [0.5 - spread, 12 + spread]);
set(ax1,'XTick',x,'XTickLabel',decomp_symbols ,'TickDir','out','FontSize',25);
xtickangle(ax1,0);
% add a Δ just left of the tick labels row
%text(ax1, 0.028, -0.017, '$\Delta$:', 'Units','normalized', ...
%     'Interpreter','latex','HorizontalAlignment','right','VerticalAlignment','top','FontSize',13);

ylim(ax1, [0,.062]);
%ylabel(ax1,'Absolute Change (d$^{-1}$)','FontSize',25);

% Add Annotations
for labels =1
% --- Grazing/Non Grazing Labels ---
yl = ylim(ax1); yTop = yl(2);
yPad = -0.003;               
fs   = 20;                     

for i = 1:numel(x)
    if i==11
        %lbl = '\textbf{all lossess}';
        lbl = '';
    elseif mod(i,2)==1
        lbl = '\textbf{grazing}';
    else
        lbl = '\textbf{non-grazing}';
    end

    text(ax1, x(i), yTop - yPad, lbl, ...
        'Color', com_col(i,:), ...
        'HorizontalAlignment','center', ...
        'VerticalAlignment','top', ...
        'FontSize', fs, ...
        'Clipping','off');    % keep it visible at the top
end


%% ---------- Right-side colored model name list + near-converging leaders ----------

pos1   = get(gca,'Position');         % [x y w h] normalized
gap    = -0.002;                       % gap between panel and list
wList  = min(0.18, 1 - (pos1(1)+pos1(3)+gap) - 0.02);   % keep on-figure
posList = [pos1(1)+pos1(3)+gap, pos1(2), wList, pos1(4)];

axList = axes('Position', posList, 'Color','none', ...
              'XLim',[0 1], 'YLim',[0 1], ...
              'XTick',[], 'YTick',[], 'Box','off', ...
              'XColor','none','YColor','none');   % no axis lines or ticks
axis(axList,'off');

% vertical spread (pad a little at top/bottom)
ys = .02:.98/nmod:.98;
ys(9) = ys(9)-.005;
ys(8) = ys(8)-.02;
ys(7) = ys(7)-.015;
% label placement + leader geometry
xText     = 0.0;          % left margin inside list axis
txtFS     = 18;   % 16
padRight  = 0.035;         % extra gap after text before line
targetX   = 0.74;          % closer right-hand focal x (steeper look)
targetY   = 0.50;          % focal y (middle of stack)
alpha     = 0.92;          % how far toward the focal point (steeper/closer)

for j = 1:nmod
    c = cmap_mod(j,:);

    % draw label
    ht = text(axList, xText, ys(j), model_names{j}, ...
        'Color', c, 'Interpreter','latex', ...
        'HorizontalAlignment','left', 'VerticalAlignment','middle', ...
        'FontSize', txtFS, 'Clipping','off');

    % start line just to the right of rendered text
    drawnow;                            % ensure 'Extent' is valid
    ext = get(ht,'Extent');             % [x y w h] in list-axis units
    x1  = ext(1) + ext(3) + padRight;   % right edge of text + tiny gap
    y1  = ext(2) + 0.5*ext(4);          % vertical center of text
    % bend toward (targetX,targetY) but stop early (alpha)
    x2 = x1 + alpha*(targetX - x1);
    y2 = y1 + alpha*(targetY - y1);

    %line(axList, [x1 x2], [y1 y2], 'Color', c, 'LineWidth', 1.3, 'Clipping','off');
end
end

end



%% Support Functions
% -------------------- Helpers ------------------------
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

% at the very end of the file (after your script code)
function y = sum2(a,b)
    y = sum([a,b], 2, 'omitnan');        % row-wise sum ignoring NaNs
    y(isnan(a) & isnan(b)) = NaN;        % if both NaN -> NaN
end




%%%%% APPENDIX A
%% IGNORE ---- Panel 2: TD - BU (Improving, worsening) u 
% OLD TD−BU anomalies Script =====
%figure('Color','w','Units','normalized','Position',[0.04 0.10 1 1]);
for p =  [];%1:2
metricIdx = 2;
spread  = [.42];                                % half-width of model band around each driver
offsets = linspace(-.5*spread, .5*spread, nmod); % fixed left->right order


    % ----- Per-model means/std -----
    Tm = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regTD(p), :, statMean));  % [nmod x 11]
    Bm = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regBU(p), :, statMean));
    Ts = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regTD(p), :, statStd));
    Bs = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regBU(p), :, statStd));

    % Pair 11->6
    TDpair_mean = [ sum2(Tm(:,1),Tm(:,2)) , sum2(Tm(:,3),Tm(:,4)) , ...
                    sum2(Tm(:,5),Tm(:,6)) , sum2(Tm(:,7),Tm(:,8)) , ...
                    sum2(Tm(:,9),Tm(:,10)), Tm(:,11) ];
    BUpair_mean = [ sum2(Bm(:,1),Bm(:,2)) , sum2(Bm(:,3),Bm(:,4)) , ...
                    sum2(Bm(:,5),Bm(:,6)) , sum2(Bm(:,7),Bm(:,8)) , ...
                    sum2(Bm(:,9),Bm(:,10)), Bm(:,11) ];
   
    % Model anomaly = TDpair − BUpair
    Ymod_mean6 = TDpair_mean - BUpair_mean;
    %Ymod_std6  = prop2std(TDpair_std, BUpair_std);

    % ----- Ensemble paired stacks (odd vs even components) -----
    Te = 100.*squeeze(Loss_Drivers(12, metricIdx, regTD(p), :, statMean));
    Be = 100.*squeeze(Loss_Drivers(12, metricIdx, regBU(p), :, statMean));

    lowSeg  = [ Te(1)-Be(1), Te(3)-Be(3), Te(5)-Be(5), Te(7)-Be(7), Te(9)-Be(9), Te(11)-Be(11) ];
    highSeg = [ Te(2)-Be(2), Te(4)-Be(4), Te(6)-Be(6), Te(8)-Be(8), Te(10)-Be(10), 0 ];
    Yens_stack_low  = lowSeg;
    Yens_stack_high = highSeg;

    % ----- Axes for this panel -----
    posP = pos3; ttl = '\textbf{TD-BU (improving growth)}';
    if p==2, posP = pos4; ttl = '\textbf{TD-BU (worsening growth)}'; end

    % Left axis: ensemble stacks
    axL = axes('Position', posP); hold(axL,'on'); box(axL,'off');
    B = bar(axL, x6, [Yens_stack_low(:) Yens_stack_high(:)], 'stacked', ...
        'BarWidth',0.72, 'EdgeColor',[0 0 0], 'LineWidth',1.0);

    oddMap  = [1 3 5 7 9 11];
    evenMap = [2 4 6 8 10 11];   % last is zero-height but keeps palette
    B(1).FaceColor = 'flat'; B(2).FaceColor = 'flat';
    B(1).CData = com_col(oddMap ,:);
    B(2).CData = com_col(evenMap,:);

    % Axis cosmetics (left / ensemble)
    xlim(axL, [0.5 - spread, 6.5 + spread]);
    set(axL,'XTick',x6,'XTickLabel',decomp_symbols3,'TickDir','out','FontSize',25);
    xtickangle(axL,0);
    yline(axL,0,'-','Color',[0 0 0],'LineWidth',0.9,'HandleVisibility','off');
    laxl = 100.*[-.25 .25]; ylim(axL, laxl);
    ylabel(axL,'Inter-model Magnitude Share Anomaly ($\%$)');
    
    % Right overlay axis: model points (stems from 0 + dots)
    axR = axes('Position', posP, 'Color','none', 'YAxisLocation','right', ...
               'XAxisLocation','bottom', 'TickDir','out', 'Box','off'); 
    hold(axR,'on');

    for j = 1:nmod
        c  = cmap_mod(j,:);
        cx = x6 + offsets(j);
        % stems
        line(axR, [cx; cx], [zeros(1,numel(cx)); Ymod_mean6(j,:)], ...
             'Color', c, 'LineWidth', lw, 'HitTest','off');
        % points
        scatter(axR, cx, Ymod_mean6(j,:), sz, c, 'filled', ...
            'MarkerEdgeColor','k', 'LineWidth',0.25);
    end

    xlim(axR, xlim(axL));
    set(axR,'XTick',x6,'XTickLabel',[],'YTick',100.*[-.6:.2:.6],'FontSize',25);
    ylim(axR, laxl.*2.5);             % give models extra headroom if you like
    ylabel(axR,'Ensemble Magnitude Share Anomaly ($\%$)');

    % share x-axis; show ticks only on left axis unless p==2
    xlim(axR, xlim(axL));
    set(axR,'XTick',x6,'XTickLabel',[]);
    ylim(axR, laxl.*2.5);

    % Tick visibility rules:
    if p == 1
        % Hide RIGHT axis ticks/label on panel 2 (improving)
        set(axR, 'YTick', [], 'YColor', 'none');
        ylabel(axR, '');
    else
        % Hide LEFT axis ticks/label on panel 3 (worsening)
        set(axL, 'YTick', [], 'YColor', 'none');
        ylabel(axL, '');
    end

end