%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. S11 Delta MLD vs Dilution Index %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))

%% Scatter: 
 %%%% 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    101  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  coex    Stbail   Zoo_m     
  
clear; clc;
bnds = [-50,50] ;

% --- Paths
basedir    = '/Users/twrohr/Desktop/data.nosync/SSP585';
regrid_dir = fullfile(basedir,'Abridged/regridded2');


% --- Model names & colors (same order you’ve been using)
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];  % same as your scripts
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];
cmap_mod = cmap_mod(1:11,:);  % keep first 11 model colors

% --- Gather files (expects one .mat per model)
cd(regrid_dir);
mods = dir('*.mat');


%% =======================
%  Panel 1: per-model scatter (as-is)
%  =======================
mld_th_sh = 30;
mld_th_dp = 50;

% Accumulator
Xall_dp = []; Yall_dp = [];
Xall_sh = []; Yall_sh = [];

% --- Plot each model
for ii = 1:11
    ii
    mi  = mod_order(ii);                % file index for this model
    S   = load(mods(mi).name, 'CoreVar_Inv_rw');
    load (mods(mi).name,'CoreVar_Inv_adj')
    V   = S.CoreVar_Inv_rw;             % (var, lat, lon)
    MLD_base = squeeze(mean(CoreVar_Inv_adj(15,:,:,:),2,'omitnan'));
    
    dMLD_dp = squeeze(V(22,:,:));      % ΔPhy biomass
    dDL_dp = squeeze(V(21,:,:));      % ΔZoo biomass
    dMLD_sh = squeeze(V(22,:,:));      % ΔPhy biomass
    dDL_sh = squeeze(V(21,:,:));      % ΔZoo biomass

    % Sahllow
    dMLD_dp(MLD_base>mld_th_sh) = nan;
    dDL_dp(MLD_base>mld_th_sh) = nan;

    % Flatten
    x_dp = dMLD_dp(:);  
    y_dp = dDL_dp(:);

    % Save for panel 2
    Xall_dp = [Xall_dp; x_dp];
    Yall_dp = [Yall_dp; y_dp];

    % Deep
    dMLD_sh(MLD_base<mld_th_dp) = nan;
    dDL_sh(MLD_base<mld_th_dp) = nan;

    % Flatten
    x_sh = dMLD_sh(:);  
    y_sh = dDL_sh(:);

    % Save for panel 2
    Xall_sh = [Xall_sh; x_sh];
    Yall_sh = [Yall_sh; y_sh];


end






%% 
% --- Figure with two panels
figure('Units','normalized','Position',[0 0 0.9 0.7],'Color','w');
tlo = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
xbnds = [-50,50];
ybnds = [-5,5];


% =======================
%  Panel 1: all-models bin-scatter heat map
% =======================
ax1 = nexttile; hold(ax1,'on');
set(ax1,'FontSize',14,'Box','on','Layer','top','TickLabelInterpreter','latex');
xlabel(ax1,'$\Delta$ MLD (m)','FontSize',16,'interpreter','latex');
ylabel(ax1,'$\Delta$ Dilution Length Scale (m)','FontSize',16,'interpreter','latex');
xline(ax1,0,'k-','LineWidth',0.75); 
yline(ax1,0,'k-','LineWidth',0.75);
axis(ax1,'square'); grid(ax1,'on');

% Remove NaNs/Infs before binning
valid1 = isfinite(Xall_sh) & isfinite(Yall_sh);
Xv1    = Xall_sh(valid1);
Yv1    = Yall_sh(valid1);

% Bin-scatter
hbs1 = binscatter(ax1, Xv1, Yv1, [120 120], ...
    'XLimits',xbnds, 'YLimits',ybnds, ...
    'ShowEmptyBins','off');           % avoid zero-count bins on log scale

set(ax1,'ColorScale','log');          % log color scaling
cl1 = caxis(ax1);                     % push lower limit to 1 (or >0)
caxis(ax1, [max(1,cl1(1))  cl1(2)]);

cb1 = colorbar(ax1);
ylabel(cb1,'Grid-cell count (log scale)','interpreter','latex', 'FontSize',16);
colormap(ax1, cmocean('thermal'));
set(ax1,'TickDir','out','LineWidth',1.0);

xlim(ax1, xbnds); 
ylim(ax1, ybnds);

% Fit (robust linear)
mdl1 = fitlm(Xv1, Yv1, 'linear', 'RobustOpts','on');
R2_1 = mdl1.Rsquared.Ordinary; 

% Prediction grid over your panel limits
xg = linspace(-100, 100, 400)';
[yhat1, yCI1] = predict(mdl1, xg, 'Prediction','curve');   % 95% CI

% Thicker regression line & CI
plot(ax1, xg, yhat1,    'k-',  'LineWidth', 3);   % main fit line
plot(ax1, xg, yCI1(:,1),'k--', 'LineWidth', 2);   % lower CI
plot(ax1, xg, yCI1(:,2),'k--', 'LineWidth', 2);   % upper CI

% Correlation coefficient (Pearson r)
r1 = corr(Xv1, Yv1, 'rows','complete');

% Title with r in parentheses
title(ax1, sprintf('Deep MLDs ($R^2 = %.2f$)', R2_1), ...
    'interpreter','latex', 'FontSize',20);


% =======================
% Quadrant percentages for Panel 1
% =======================
N1   = numel(Xv1);
q1_1 = (Xv1 > 0) & (Yv1 > 0);   % Q1: x>0, y>0 (top-right)
q2_1 = (Xv1 < 0) & (Yv1 > 0);   % Q2: x<0, y>0 (top-left)
q3_1 = (Xv1 < 0) & (Yv1 < 0);   % Q3: x<0, y<0 (bottom-left)
q4_1 = (Xv1 > 0) & (Yv1 < 0);   % Q4: x>0, y<0 (bottom-right)

pct1_1 = 100 * sum(q1_1) / N1;
pct2_1 = 100 * sum(q2_1) / N1;
pct3_1 = 100 * sum(q3_1) / N1;
pct4_1 = 100 * sum(q4_1) / N1;

% Positions very close to actual corners (small padding)
xL1 = xlim(ax1); 
yL1 = ylim(ax1);
xr1 = diff(xL1);
yr1 = diff(yL1);

padX = 0.02 * xr1;
padY = 0.02 * yr1;

x_left_1   = xL1(1) + padX;
x_right_1  = xL1(2) - padX;
y_bottom_1 = yL1(1) + padY;
y_top_1    = yL1(2) - padY;

textArgs = {'FontSize',18,'FontWeight','bold','Interpreter','latex', ...
            'Color','k','BackgroundColor','w','Margin',3, ...
            'EdgeColor','k','LineWidth',1.2};
% Q2 (top-left)
text(ax1, x_left_1,  y_top_1,    sprintf('%.1f\\%%', pct2_1), ...
     textArgs{:}, 'HorizontalAlignment','left',  'VerticalAlignment','top');

% Q1 (top-right)
text(ax1, x_right_1, y_top_1,    sprintf('%.1f\\%%', pct1_1), ...
     textArgs{:}, 'HorizontalAlignment','right', 'VerticalAlignment','top');

% Q3 (bottom-left)
text(ax1, x_left_1,  y_bottom_1, sprintf('%.1f\\%%', pct3_1), ...
     textArgs{:}, 'HorizontalAlignment','left',  'VerticalAlignment','bottom');

% Q4 (bottom-right)
text(ax1, x_right_1, y_bottom_1, sprintf('%.1f\\%%', pct4_1), ...
     textArgs{:}, 'HorizontalAlignment','right', 'VerticalAlignment','bottom');



% =======================
%  Panel 2: all-models bin-scatter heat map
% =======================
ax2  = nexttile; hold(ax2,'on');
set(ax2,'FontSize',14,'Box','on','Layer','top','TickLabelInterpreter','latex');
xlabel(ax2,'$\Delta$ MLD (m)','FontSize',16,'interpreter','latex');
ylabel(ax2,'$\Delta$ Dilution Length Scale (m)','FontSize',16,'interpreter','latex');
xline(ax2,0,'k-','LineWidth',0.75); 
yline(ax2,0,'k-','LineWidth',0.75);
axis(ax2,'square'); grid(ax2,'on');

% Remove NaNs/Infs before binning
valid2 = isfinite(Xall_dp) & isfinite(Yall_dp);
Xv2    = Xall_dp(valid2);
Yv2    = Yall_dp(valid2);

% Bin-scatter
hbs2 = binscatter(ax2, Xv2, Yv2, [120 120], ...
    'XLimits',xbnds, 'YLimits',ybnds, ...
    'ShowEmptyBins','off');           % avoid zero-count bins on log scale

set(ax2,'ColorScale','log');          % log color scaling
cl2 = caxis(ax2);                     % push lower limit to 1 (or >0)
caxis(ax2, [max(1,cl2(1))  cl2(2)]);

cb2 = colorbar(ax2);
ylabel(cb2,'Grid-cell count (log scale)','interpreter','latex', 'FontSize',16);
colormap(ax2, cmocean('thermal'));
set(ax2,'TickDir','out','LineWidth',1.0);

xlim(ax2, xbnds); 
ylim(ax2, ybnds);

% Fit (robust linear)
mdl2 = fitlm(Xv2, Yv2, 'linear', 'RobustOpts','on');
R2_2 = mdl2.Rsquared.Ordinary; 

% Prediction grid over your panel limits (reuse xg)
[yhat2, yCI2] = predict(mdl2, xg, 'Prediction','curve');   % 95% CI

% Thicker regression line & CI
plot(ax2, xg, yhat2,    'k-',  'LineWidth', 3);   % main fit line
plot(ax2, xg, yCI2(:,1),'k--', 'LineWidth', 2);   % lower CI
plot(ax2, xg, yCI2(:,2),'k--', 'LineWidth', 2);   % upper CI

% Correlation coefficient (Pearson r)
r2 = corr(Xv2, Yv2, 'rows','complete');

% Title with r in parentheses
title(ax2, sprintf('Shallow MLDs ($R^2 = %.2f$)', R2_2), ...
    'interpreter','latex', 'FontSize',20);


% =======================
% Quadrant percentages for Panel 2
% =======================
N2   = numel(Xv2);
q1_2 = (Xv2 > 0) & (Yv2 > 0);   % Q1: x>0, y>0 (top-right)
q2_2 = (Xv2 < 0) & (Yv2 > 0);   % Q2: x<0, y>0 (top-left)
q3_2 = (Xv2 < 0) & (Yv2 < 0);   % Q3: x<0, y<0 (bottom-left)
q4_2 = (Xv2 > 0) & (Yv2 < 0);   % Q4: x>0, y<0 (bottom-right)

pct1_2 = 100 * sum(q1_2) / N2;
pct2_2 = 100 * sum(q2_2) / N2;
pct3_2 = 100 * sum(q3_2) / N2;
pct4_2 = 100 * sum(q4_2) / N2;

% Positions close to corners
xL2 = xlim(ax2); 
yL2 = ylim(ax2);
xr2 = diff(xL2);
yr2 = diff(yL2);

padX2 = 0.02 * xr2;
padY2 = 0.02 * yr2;

x_left_2   = xL2(1) + padX2;
x_right_2  = xL2(2) - padX2;
y_bottom_2 = yL2(1) + padY2;
y_top_2    = yL2(2) - padY2;

% reuse textArgs

% Q2 (top-left)
text(ax2, x_left_2,  y_top_2,    sprintf('%.1f\\%%', pct2_2), ...
     textArgs{:}, 'HorizontalAlignment','left',  'VerticalAlignment','top');

% Q1 (top-right)
text(ax2, x_right_2, y_top_2,    sprintf('%.1f\\%%', pct1_2), ...
     textArgs{:}, 'HorizontalAlignment','right', 'VerticalAlignment','top');

% Q3 (bottom-left)
text(ax2, x_left_2,  y_bottom_2, sprintf('%.1f\\%%', pct3_2), ...
     textArgs{:}, 'HorizontalAlignment','left',  'VerticalAlignment','bottom');

% Q4 (bottom-right)
text(ax2, x_right_2, y_bottom_2, sprintf('%.1f\\%%', pct4_2), ...
     textArgs{:}, 'HorizontalAlignment','right', 'VerticalAlignment','bottom');


%% 
% --- Figure with two panels
figure('Units','normalized','Position',[0 0 0.9 0.7],'Color','w');
tlo = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');


% =======================
%  Panel 1: all-models bin-scatter heat map
%  =======================
ax1 = nexttile; hold(ax1,'on');
set(ax1,'FontSize',14,'Box','on','Layer','top','TickLabelInterpreter','latex');
xlabel(ax1,'$\Delta$ MLD (m)','FontSize',16,'interpreter','latex');
ylabel(ax1,'$\Delta$ Dilution Length Scale (m)','FontSize',16,'interpreter','latex');
xline(ax1,0,'k-','LineWidth',0.75); 
yline(ax1,0,'k-','LineWidth',0.75);
axis(ax1,'square'); grid(ax1,'on');

% Remove NaNs/Infs before binning
valid = isfinite(Xall_sh) & isfinite(Yall_sh);
Xv = (Xall_sh(valid));  Yv = (Yall_sh(valid));

% Bin-scatter
% after creating ax2 and Xv,Yv ...
hbs = binscatter(ax1, Xv, Yv, [120 120], ...
    'XLimits',bnds, 'YLimits',bnds, ...
    'ShowEmptyBins','off');           % avoid zero-count bins on log scale

 set(ax1,'ColorScale','log');          % log color scaling
 cl = caxis(ax1);                      % push lower limit to 1 (or >0)
 caxis(ax1, [max(1,cl(1))  cl(2)]);
% caxis(ax1, [0,300]);

cb = colorbar(ax1);
ylabel(cb,'Grid-cell count (log scale)','interpreter','latex', 'FontSize',16);
colormap(ax1, cmocean('thermal'));
set(ax1,'TickDir','out','LineWidth',1.0);

xlim(bnds); ylim(bnds)
%clim([1,200])

% Fit
mdl = fitlm(Xv, Yv, 'linear', 'RobustOpts','on');

% Prediction grid over your panel limits
xg = linspace(-100, 100, 400)';
[yhat, yCI] = predict(mdl, xg, 'Prediction','curve');   % 95% CI
h_fit = plot(ax1, xg, yhat, 'k-',  'LineWidth', 2);

% Optional: annotate slope/intercept and add legend
b = mdl.Coefficients.Estimate;  % [b0; b1]
title(ax1, 'Deep MLDs','interpreter','latex', 'FontSize',20);

hold(ax1,'on');
plot(ax1, xg, yhat,   'k-',  'LineWidth', 2);
plot(ax1, xg, yCI(:,1),'k--', 'LineWidth', 1.2);
plot(ax1, xg, yCI(:,2),'k--', 'LineWidth', 1.2);


% =======================
%  Panel 2: all-models bin-scatter heat map
%  =======================
ax2  = nexttile; hold(ax2,'on');
set(ax2,'FontSize',14,'Box','on','Layer','top','TickLabelInterpreter','latex');
xlabel(ax2,'$\Delta$ MLD (m)','FontSize',16,'interpreter','latex');
ylabel(ax2,'$\Delta$ Diltuion Length Scale (m)','FontSize',16,'interpreter','latex');
xline(ax2,0,'k-','LineWidth',0.75); 
yline(ax2,0,'k-','LineWidth',0.75);
axis(ax2,'square'); grid(ax2,'on');

% Remove NaNs/Infs before binning
valid = isfinite(Xall_dp) & isfinite(Yall_dp);
Xv = (Xall_dp(valid));  Yv = (Yall_dp(valid));

% Bin-scatter
% after creating ax2 and Xv,Yv ...
hbs = binscatter(ax2, Xv, Yv, [120 120], ...
    'XLimits',bnds, 'YLimits',bnds, ...
    'ShowEmptyBins','off');           % avoid zero-count bins on log scale

 set(ax2,'ColorScale','log');          % log color scaling
 cl = caxis(ax2);                      % push lower limit to 1 (or >0)
 caxis(ax2, [max(1,cl(1))  cl(2)]);
 %caxis(ax2, [0,200]);

cb = colorbar(ax2);
ylabel(cb,'Grid-cell count (log scale)','interpreter','latex', 'FontSize',16);
colormap(ax2, cmocean('thermal'));
set(ax2,'TickDir','out','LineWidth',1.0);

xlim(bnds); ylim(bnds)
%clim([1,200])


% Fit
mdl = fitlm(Xv, Yv, 'linear', 'RobustOpts','on')

% Prediction grid over your panel limits
xg = linspace(-100, 100, 400)';
[yhat, yCI] = predict(mdl, xg, 'Prediction','curve');   % 95% CI
h_fit = plot(ax2, xg, yhat, 'k-',  'LineWidth', 2);

% Optional: annotate slope/intercept and add legend
b = mdl.Coefficients.Estimate;  % [b0; b1]
title(ax2, 'Deep MLDs','interpreter','latex', 'FontSize',20);

hold(ax2,'on');
plot(ax2, xg, yhat,   'k-',  'LineWidth', 2);
plot(ax2, xg, yCI(:,1),'k--', 'LineWidth', 1.2);
plot(ax2, xg, yCI(:,2),'k--', 'LineWidth', 1.2);
