%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Case Studies - Figure 6b-g and supporting Figures S7-S9 %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To plot a case study, the use must select the model, region, and varibales.
% This can be starting in line 100 in the section 'Select use inputs'
% Please be sure to update all file directories and data paths as required (see 'Read Me' available online at: https://data.imas.utas.edu.au/attachments/ffbba75d-56dd-49b5-8f69-458b731ce91e/

clc; close all
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs

%% Startup
 clearvars last_load
for see_startup =[1]
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/SOLACE'))
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))
clearvars plot sum

% Model Info and Colormaps
mods = dir;
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]];

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


CoreVar_Names = {'$\Sigma$ Phytoplankton Biomass',     '$\Sigma$ Zooplankton Biomass',         'NPP',                  'Grazing Losses',            'EP',                 ...
                 {'Phytotoplankton Specific';'Division Rate'},   {'Zooplankton Specific';'Grazing Rate'},          'GSP',                  '{Zooplankton Loss Rate}',               'Non-Grazing Losses', ...
                 'Specific Grazing Loss Rate',   {'Small Phytotoplankton';'Compositon'},      {'Small Zooplankton';'Compositon'},     'SST' ,                      'MLD' ,               ...
                 {'Phytotoplankton Specific';'Loss Rate'},           {'Phytoplankton Specific'; 'Non-grazing Loss Rate'}  'Export Efficeincy'     'Ecosystem Efficeincy'       '[Phyto]$_{\mathrm{Surf}}$' ,  ...
                 'DBM [Phyto]',               'DBM Depth',                   '$\overline{\textrm{[Phy]}}$',       'Diaz Frac',                 'Diatom Frac',    'Non- Diatom Frac'...
                 'dgdT',                      'dgdVD',                       'dgdCC',                'dgdP',                      'dgdZ',...
                 'dngdT',                     'dngdVD',                      'dngdCC',               'dngdP',                     'dngdZ',...
                 'dldT',                     'dldVD',                      'dldCC',               'dldP',                     'dldZ',...
                 'Dilution Index'};


CoreVar_Sym    = {'$\Sigma$P',                '$\Sigma$Z',                   'NPP',                  'G',                         'EP',   ...
                  '$\mu$',                    'g',                           'GSP',                  'Z$_M$',                     'NGL', ...
                  'gp',                       'f$_{SP}$',                    'f_{SZ}',               'SST'                        'MLD', ...
                  'l',                        'ngl'                          'e-rat'                 'eco-rat'};


CoreVar_Inv_Units = {'(mmol m$^{-2}$)',            '(mmol m$^{-2}$)',     '(mmol m$^{-2}$ d$^{-1}$)',    '(mmol m$^{-2}$ d$^{-1}$)',     '(mmol m$^{-2}$ d$^{-1}$)'...
                     '(d$^{-1}$)',                 '(d$^{-1}$)',          '(mmol m$^{-2}$ d$^{-1}$)',    '(mmol m$^{-2}$ d$^{-1}$)',     '(mmol m$^{-2}$ d$^{-1}$)' ...
                      '(d$^{-1}$)',                '($\%$)',               '($\%$)'  ,                  '($^{\circ}$C)'                     '(m)'...
                      '(d$^{-1}$)',             '(d$^{-1}$)',              '(frac)'                        '(frac)'   ,                  '(mmol m$^{-3}$)',...
                      '(mmol m$^{-3}$)',           '(m)',                  '(mmol m$^{-3}$)'            ,  '(frac)',                       '(frac)',               '(frac)',...
                      '(d$^{-1}$)',              '(d$^{-1}$)',             '(d$^{-1}$)',               '(d$^{-1}$)',                      '(d$^{-1}$)',... 
                      '(d$^{-1}$)',              '(d$^{-1}$)',             '(d$^{-1}$)',               '(d$^{-1}$)',                      '(d$^{-1}$)'...
                      '(d$^{-1}$)',              '(d$^{-1}$)',             '(d$^{-1}$)',               '(d$^{-1}$)',                      '(d$^{-1}$)' ...
                      'm'}; 
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
];
lighter_colors = wedge_colors + 0.6 .* (1 - wedge_colors);  % shift toward white, preserve hue


end


%% Plot Case Studies
clc; 
delete(findall(groot,'Type','figure'));  % closes visible + hidden figures
drawnow;                                  

for i=[1]
% Specify User Inputs
    
    %%%% Pick Model Index: 1. iHAMMOCC 2. CMOC 3. WOMBAT 4. OECOv2 5. MARBL 6. CanOE 7. MEDUSA2.1 8. PISCESv2-IPSL 9. PISCESv2-CNRM 10. BFM5.2 11. COBALTv2'
         model = 10; 
    
    %%%% 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  12. Tropics
         region_index = 6; 
    
    %%%% Plotting Specs:
         % Plot Map? 
         map = true;
         % Plot Profiles?
         profiles = false; mon = 2; %l 1. Jan 2. July

    %%%% Pick Filters
         % Phytoplankton Biomass increase or decrease 
         %inc =true;   % Filter for Increasing P Biomass
         inc =false;   % Filter for Decreasing P Biomass

         % Top-down or Bottom-up Control
         TD = true;  % Filter for Top down control
         %TD = false;  % Filter for Bottom up control
   
    %%%% Pick 8 plotting vairables
         %    1      2     3    4    5    6    7   8      9     10   11   12   13    14    15    
         %  P Bio  Z bio  NPP   G    EP   u    g  GSP   Zoo_M   NGL  gp   fsP  fsZ   SST   MLD    
      
         %   16    17     18      19      20     21      22      23       24      25      26     
         %    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            
         %  dgdT  dgdVD dgdCC  dgdP  dgdZ   dngldx  dldx    Dil Dep  co-loc in  Clear Rate  

         %   45          46        70         71       72      100    101       102       119   
         % Clear (VD) [Phy]MLD  Mechanism  ngl & gp  P and Z  circle  coex    Stbail     Zoo_m 
         vlist = [15, 6,  12, 11, 70 ...  
                 14, 71, 13, 17, 100];  
         
    %%%% Pick 8 plotting vairables 
         %  1/6. Temp | 2/7 VD |  3. CC | 4/9 Phy | 5/10. Zoo
         bor   = [nan, nan, nan, 3, nan, ...
                  nan, nan, nan, 3, nan];

% Specify Label Names
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];
monL = {'J','F','M','A','M','J','J','A','S','O','N','D'};

% Step 0. Load Model if Required
    if exist('last_load','var') && last_load == model 
          cd /Users/twrohr/Desktop/Data.nosync/SSP585/Abridged/regridded2; 
          mods = dir; mods = mods(~[mods.isdir]);
         fprintf(['\n Already Dowloaded ' mods(mod_order(model)).name])
         cd /Users/twrohr/Desktop/data.nosync/SSP585
          load('Longhurst360.mat');
          Lon =  Lon_t; Lat = Lat_t;
    else
        cd /Users/twrohr/Desktop/data.nosync/SSP585
          load('Longhurst360.mat');
          Lon =  Lon_t; Lat = Lat_t;
        cd /Users/twrohr/Desktop/Data.nosync/SSP585/Abridged/regridded2; 
          mods = dir; mods = mods(~[mods.isdir]);
          load(mods(mod_order(model)).name, 'CoreVar_Inv_adj','CoreVar_Inv_pc', 'CoreVar_Names', 'CoreVar_Inv_Units','Sens'); 
          fprintf(['\n Dowloaded ' mods(mod_order(model)).name])
          last_load=model;
    end
    CoreVar_Names{23}='$\overline{\textrm{[Phy]}}$';


% Step 1: Apply Longhurst mask: 
    if region_index ==10
        region_mask = double(Long_Prov == 8 | Long_Prov == 2);
        %region_mask(  Lat > -46 ) = logical(0);
        %region_mask(  Lat < -68    | Lat>-40  ) = logical(0); %For Pisces VD panel
    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);
    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); % Remove SH% Remove Pacific !!Latmin bound 40 for CESM; 20 normalls
    elseif region_index == 12
        region_mask = double(Long_Prov == 3 | Long_Prov == 4 | Long_Prov == 5);
        region_mask = logical(region_mask);
    else
    region_mask = double(Long_Prov == region_index);
    region_mask = logical(region_mask);
    end
    %region_mask((Lon>95 & Lon<150)   |   Lat<-26 ) =logical(0);

    %region_mask(Lat<0) =logical(0); % Remove NH
    
% Step 3: Apply other logical filers 
     % 3a. Biomass Change Masks 
       incPmask = squeeze(CoreVar_Inv_pc(1,:,:));  incPmask(incPmask<0)=nan; % Increasing Biomass (removes Decreasing Biomass) % 42 for MARBL 
       decPmask = squeeze(CoreVar_Inv_pc(1,:,:));  decPmask(decPmask>0)=nan; % Decreasing Biomass (removes Increasing Biomass) % 42 for MARBL 
       % Extra Filters
        % Piscesv2
        LMmask1 = squeeze(CoreVar_Inv_pc(28,:,:)); 
        LMmask1(LMmask1<15) =nan;
        %dldz=squeeze(Sens(4,:,:));
        %decPmask = dldz;  decPmask(decPmask>0)=nan; % Decreasing Biomass (removes Increasing Biomass) % 42 for MARBL 
       
     % 3b. Bottom-Up v. Top down Control Masks 
       TDvar   =  sign(squeeze(CoreVar_Inv_pc(1,:,:))...
                  .*squeeze(CoreVar_Inv_pc(6,:,:)));
       TDmask = TDvar;  TDmask(TDmask>0)=nan; % Top Down Control 
       BUmask = TDvar;  BUmask(BUmask<0)=nan; % Bottom-Up Control

     % 3c. Dominate Loss Driver Mask 
       %LMmask1 = Dom_Mech_all;  LMmask(LMmask~=1)=nan; % T on Grazing
       %LMmask2 = Dom_Mech_all;  LMmask(LMmask~=2)=nan; % Vd on Grazing
       %LMmask3 = Dom_Mech_all;  LMmask(LMmask~=3)=nan; % CC on Grazing

     % Select Other Filters
       BinDum  = ones(size(Lat));
       % Pick Filters 
       if inc 
        BinDum(isnan(incPmask))  = 0; % Increasing Biomass (Removes Decreasing) 
       else
        BinDum(isnan(decPmask))  = 0;  cb_map = cmapTD;% Decreasing Biomass (Removes Increasing) 
       end
       if TD 
        BinDum(isnan(TDmask))    = 0; cb_map = cmapTD;   % Top-Down Control (removes Bottom-up) 
       else
        BinDum(isnan(BUmask))    = 0; cb_map = cmapBU;   % Bottom-Up Control (removes Top-down) 
       end
        % BinDum(isnan(LMmask1))   = 0; % Specify extra filter
     
    % Combine with Spatial Mask
      FullMask = region_mask & logical(BinDum);
      %FullMask = region_mask;%  & logical(BinDum);

      % Step 4. Create Climatologies (biomass-weighted spatial MEDIANS)
clearvars TSrw_2016 TSrw_2096 TSrw_reg

% Dimensions
[nVar, nT, nLat, nLon] = size(CoreVar_Inv_adj);

% Spatial mask as linear index into (lat,lon)
mask_lin = FullMask(:);                    % logical nLat*nLon
nMask    = nnz(mask_lin);

% Preallocate: var x time
TSrw_reg = nan(nVar, nT);

% ---- Spatial median (weighted by phy biomass = var #1) for each time ----
% Step 4. Create Climatologies (biomass-weighted spatial MEANS)
clearvars TSrw_2016 TSrw_2096 TSrw_reg

% Dimensions
[nVar, nT, nLat, nLon] = size(CoreVar_Inv_adj);

% Spatial mask as linear index into (lat,lon)
mask_lin = FullMask(:);                    % logical nLat*nLon

% Preallocate: var x time
TSrw_reg = nan(nVar, nT);


% Step 4. Create Climatologies
  clearvars TSrw_2016 TSrw_2096 TSrw_reg

    % 4a. Average over Spatial Mask
      TSrw_dum = CoreVar_Inv_adj(:,:,FullMask);
      TSrw_dum(isinf(TSrw_dum)) = nan; 
      tmp = TSrw_dum(42,:,:); tmp(tmp >200) = NaN; TSrw_dum(42,:,:) = tmp; 
      tmp = TSrw_dum(44,:,:); tmp(tmp >.2) = NaN; TSrw_dum(44,:,:) = tmp; 
      TSrw_reg = squeeze(nanmean(TSrw_dum,3));
      
    % 4b.  Make Climatoloies
      for m =1:12
           TSrw_2016(:,m)  =  mean(TSrw_reg(:,m:12:12*5),2);  TSrw_2096(:,m)  = mean(TSrw_reg(:,60+m:12:end),2);
      end

    % 4c. Shift to Center on Summer if South of Equator 
      if mean(Lat(FullMask),'omitnan')<0
          pTSrw_2096 = circshift(TSrw_2096,6,2); pTSrw_2016 = circshift(TSrw_2016,6,2); pmonL = circshift(monL,6);
      else
          pTSrw_2096 = TSrw_2096; pTSrw_2016 = TSrw_2016; pmonL =  monL;
      end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%%% Step 5. Plot Timeseries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
for plotreg = [1]
    % tsize=35; ysize =30; ticsize = 15;
     tsize=20; ysize =15; ticsize = 10;

  % Set Figure Axis
   figure('Units','normalized','Position',[2 0 1.75 1]); %clc 8 panel (2 0 1.75 1) for full thing (then shrink); (2 0 1.75 .55) for top row 
   set(gcf,'color','w');
   set(gca, 'FontSize', 30,'TickLabelInterpreter','latex') % Set font size for ticks
     
  % Loop Through Panels
  t = tiledlayout(2, 5); % 5 %,'Padding', 'compact'
  for s =1:10
      nexttile

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%% Circle Plot %%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if vlist(s) == 100
            
                %--- set up the tile axes ---
             ax = gca; cla(ax); hold(ax,'on');

            axis(ax,'off');
            ax.XColor = 'none'; ax.YColor = 'none';
            
            % Key bits:
            ax.Units = 'normalized';
            ax.PositionConstraint = 'outerposition';   % important for tiledlayout
            ax.LooseInset = [0 0 0 0];                 % remove padding
            ax.Clipping = 'off';
            
            axis(ax,'equal');                          % keep circle circular
            
            
            %ax = gca; cla(ax); hold(ax,'on');
           % axis(ax,'equal'); axis(ax,'off'); set(ax,'XColor','none','YColor','none');
           
            % Compute Res
            phy_W = pTSrw_2016(1,:);  
            % phy_W = 1+0*pTSrw_2016(1,:); 
            res=pTSrw_2096(16,:)-pTSrw_2016(16,:);
            for pp=1:11
            if pp < 11
            dum          =  pTSrw_2096(26+pp,:)-pTSrw_2016(26+pp,:);
            dmeans(pp)  =  wmean(dum,phy_W);
            res = res - dum;
            else 
            dmeans(11)  =  wmean(res,phy_W);
            end
            end
            % Driver Magntiude
            tot_l = sum(abs(dmeans) ); 
            %tot_l = sum((dmeans) ); 
            rel_mag = 100.*dmeans./tot_l;

            %--- colors (make sure you have 5) ---
            if ~exist('wedge_colors','var') || size(wedge_colors,1) < 5
                wedge_colors  = lines(5);
            end
            if ~exist('lighter_colors','var') || size(lighter_colors,1) < 5
                lighter_colors = min(1, wedge_colors + 0.45);
            end
        
            % ---------- Rename ----------
            vals_g  = rel_mag(1:5);    % 1x5
            vals_ng = rel_mag(6:10);   % 1x5
            val_res = rel_mag(11);         % scalar
            
            % ---------- geometry (exactly 6 wedges) ----------
            n_wedges = 6;                        % 5 stacked + 1 residual
            the0     = deg2rad(90);            % start angle
            th       = linspace(the0, the0 - 2*pi, n_wedges+1);  % 7 edges -> 6 sectors
            th_bg    = linspace(0, 2*pi, 300);
            
            % reference radius from magnitudes (include residual)
            %max_neg = max(abs([vals_g, vals_ng, val_res]), [], 'omitnan');
            max_neg = abs(min(([vals_g, vals_ng, val_res])));
            ref     = max(10, ceil(max_neg/10)*10);
            max_val = max(abs([vals_g, vals_ng, val_res]), [], 'omitnan');
            max_r   = ref + 1.1*max_val;
            
            % faint inner disk
            for k = 1:numel(th_bg)-1
                t  = linspace(th_bg(k), th_bg(k+1), 3);
                fill([0*cos(t), ref*cos(fliplr(t))], ...
                     [0*sin(t), ref*sin(fliplr(t))], ...
                     0.5*[1 1 1], 'EdgeColor','none','FaceAlpha',0.25, 'Parent', ax);
            end
            
           
            for drawwedge=[1]
            % --- helper: draw a single stacked wedge from r0→r1 over angle t0→t1 ---
            draw_wedge = @(r0, r1, t0, t1, col) fill( ...
                [r0*cos(linspace(t0,t1,50)), ...
                 r1*cos(linspace(t1,t0,50))], ...
                [r0*sin(linspace(t0,t1,50)), ...
                 r1*sin(linspace(t1,t0,50))], ...
                col, 'EdgeColor','k', 'Parent', ax);

            % ----- desired order of variables around the circle -----
            order = [4 5 1 3 2];  % maps sector positions 1..5 to var indices
    
            % ---------- draw 5 stacked sectors (positions 1..5) ----------
            for pos = 1:5
                idx = order(pos);        % which variable goes into this sector
                t0  = th(pos);
                t1  = th(pos+1);
            
                g   = vals_g(idx);
                ng  = vals_ng(idx);
            
                % colors follow variable index
                col_g  = wedge_colors(idx,:);
                col_ng = lighter_colors(idx,:);
            
                % --- initialize radii for this sector ---
                r_inner = ref;   % baseline (center reference)
                r_outer = ref;   % running outer boundary
            
                % --- handle stacking logic ---
                vals_to_draw = [];  cols_to_draw = [];
            
                % separate positive and negative stacking directions
                if g >= 0
                    vals_to_draw = [vals_to_draw; g];
                    cols_to_draw = [cols_to_draw; col_g];
                else
                    vals_to_draw = [vals_to_draw; g];
                    cols_to_draw = [cols_to_draw; col_g];
                end
            
                if ng >= 0
                    vals_to_draw = [vals_to_draw; ng];
                    cols_to_draw = [cols_to_draw; col_ng];
                else
                    vals_to_draw = [vals_to_draw; ng];
                    cols_to_draw = [cols_to_draw; col_ng];
                end
            
                % --- draw positive stack outward, negative stack inward ---
                % positives
                r_pos = ref;
                for ii = find(vals_to_draw > 0)'
                    r1 = r_pos;
                    r2 = r_pos + vals_to_draw(ii);
                    draw_wedge(r1, r2, t0, t1, cols_to_draw(ii,:));
                    r_pos = r2;
                end
            
                % negatives (draw inward)
                r_neg = ref;
                for ii = find(vals_to_draw < 0)'
                    r1 = r_neg;
                    r2 = r_neg + vals_to_draw(ii);
                    draw_wedge(r2, r1, t0, t1, cols_to_draw(ii,:));  % note reversed order
                    r_neg = r2;
                end
            end
            
            
                        % ---------- residual stays in sector 6 ----------
            res_idx = 6;
            t0 = th(res_idx); 
            t1 = th(res_idx+1);
            
            if isfinite(val_res)
                col_res = [0.3 0.3 0.3];   % dark grey/black
            
                if val_res >= 0
                    % stack outward from ref
                    r1 = ref;
                    r2 = ref + val_res;
                    draw_wedge(r1, r2, t0, t1, col_res);
                else
                    % stack inward from ref
                    r1 = ref + val_res;
                    r2 = ref;
                    draw_wedge(r1, r2, t0, t1, col_res);
                end
            end
  end

        
            %--- radial tick guides centered on ref (−, 0, +) ---
            r_extent = ref + max(max_val,[],'omitnan');
            pref     = [5 25 50 100 150];
            [~,k]    = min(abs(pref - (r_extent/3)));
            tick_spacing = pref(k);
            r_ticks  = -tick_spacing:tick_spacing:tick_spacing;
            r_ticks(1) = max(-ref, r_ticks(1));
            % for rr = r_ticks
            %     rt = ref + rr; if rt <= 0, continue; end
            %     plot(ax, rt*cos(th_bg), rt*sin(th_bg), 'k:', 'LineWidth', 0.35);
            %     text(ax, rt, 0, sprintf('%+d', rr), 'HorizontalAlignment','center', ...
            %         'VerticalAlignment','bottom', 'Interpreter','latex','FontSize',20);
            % end
            tol = max(1e-9, 1e-6*ref);     % numerical tolerance
            didZero = false;               % ensure we only draw the center once
            
            for rr = r_ticks
                rt = ref + rr;             % absolute radius for this tick
                if rt <= tol && ~didZero
                    % dot at center
                    plot(ax, 0, 0, 'k.', 'MarkerSize', 16);
                
                    % center label "0%" right at the origin
                    tstr = sprintf('$%d\\%%$', rr);   % rr = 0 here
                    text(ax, 0, 0, tstr, ...
                        'Interpreter','latex', ...
                        'HorizontalAlignment','center', ...
                        'VerticalAlignment','bottom', ...
                        'Clipping','off', ...
                        'FontSize', 20, 'Color','k');
                
                    didZero = true;
                elseif rt > tol
                    % normal circular tick + label for ± values
                    plot(ax, rt*cos(th_bg), rt*sin(th_bg), 'k:', 'LineWidth', 0.35);
                    if rr==0; 
                        tstr = sprintf('$%d\\%%$', rr);
                        text(ax, 1.2*rt, 0, tstr, ...
                        'Interpreter','latex', ...
                        'HorizontalAlignment','left', ...
                        'VerticalAlignment','bottom', ...
                        'FontSize', 20);
                    else
                    tstr = sprintf('$%+d\\%%$', rr);
                    text(ax, 1.2*rt, 0, tstr, ...
                        'Interpreter','latex', ...
                        'HorizontalAlignment','center', ...
                        'VerticalAlignment','bottom', ...
                        'FontSize', 20);
                    end
                    
                end
            end

           
            title({'Mean Annual';'Magnitude Share'}, 'interpreter', 'latex','fontsize',tsize)
     

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%% Composition of P & Z Together %%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  elseif vlist(s) == 72
  % Add patches between the two lines
    x = 1:12; x_interp = 1:.1:12;
    for ii =[12,13]
    y1 = interp1(x, pTSrw_2016(ii,:), x_interp, 'linear');  % Linear interpolation for y1
    y2 = interp1(x, pTSrw_2096(ii,:), x_interp, 'linear');  % Linear interpolation for y2

    plot(x_interp,y1,'k:','linewidth',2); hold on
    plot(x_interp,y2,'k','linewidth',3); hold on

    
    % Store figure axes limits for arrow coordinate conversion
    ax = gca;
    xlim_ = xlim(ax); 
    ylim_ = ylim(ax); 
    axpos = get(ax, 'Position');  % [x0 y0 width height]

    % Add shaded patch
    for i = 1:length(x_interp)-1
        x_patch = [x_interp(i), x_interp(i+1), x_interp(i+1), x_interp(i)];
    
        if y1(i) > y2(i)
            y_patch = [y2(i), y2(i+1), y1(i+1), y1(i)];
            patch(x_patch, y_patch, 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
        elseif y2(i) > y1(i)
            y_patch = [y1(i), y1(i+1), y2(i+1), y2(i)];
            patch(x_patch, y_patch, 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
        end
    end

        % Add 12 custom arrows between lines with triangle heads
        for add_arrows=1
        x_arrow_pos = linspace(x_interp(1), x_interp(end), 12);
        
        arrow_head_height = range(ylim_)./20;  % in y-axis units
        arrow_head_width = .3; % in x-axis units
        
        ylims = ylim(gca);               % get current y-axis limits
        y_threshold = 0.08 * range(ylim_);  % 15% of y-range
        
        for j = 2:length(x_arrow_pos)
            [~, idx] = min(abs(x_interp - x_arrow_pos(j)));
        
            x0 = x_interp(idx);
            y1_val = y1(idx);
            y2_val = y2(idx);
            dy = y2_val - y1_val;
        
            if abs(dy) < y_threshold
                continue  % Skip if vertical difference is too small
            end
        
            % Arrow shaft: always from y1 to y2, minus arrowhead height
            shaft_top = y2_val - sign(dy)*arrow_head_height;
            
            % Arrowhead triangle
            x_head = [x0 - arrow_head_width/2, x0 + arrow_head_width/2, x0];
            if dy > 0  % Upward
                line([x0, x0], [y1_val, shaft_top], 'Color', 'b', 'LineWidth', 1);
        
                y_head = [y2_val - arrow_head_height, ...
                          y2_val - arrow_head_height, ...
                          y2_val];
                 fill(x_head, y_head, 'b', 'EdgeColor', 'none');
        
            else       % Downward
                line([x0, x0], [y1_val, shaft_top], 'Color', 'r', 'LineWidth', 1);
        
                y_head = [y2_val + arrow_head_height, ...
                          y2_val + arrow_head_height, ...
                          y2_val];
        
                 fill(x_head, y_head, 'r', 'EdgeColor', 'none');
            end
           
        end
        end

    end
      
    % Bring the lines back on top of the patches
      uistack(findobj(gca, 'Type', 'line'), 'top');
      % Adjust labels
      xlim([1,12]); %ylim([0,.33])
      xticklabels(pmonL)
      ax = gca; ax.TickLabelInterpreter = 'latex'; ax.FontSize = ticsize;
      title(strcat(num2str(round(pc)), 'Ccommunity Composition  '), 'interpreter', 'latex','fontsize',16)
      title('Community Composition', 'interpreter', 'latex','fontsize',tsize)
      ylabel(CoreVar_Inv_Units(12), 'interpreter', 'latex','fontsize',ysize)
      



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%% Patrials %%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  elseif vlist(s) == 70
  res = (pTSrw_2096(16,:)-pTSrw_2016(16,:))    ;
  for pp=1:5
  plot(pTSrw_2096(26+pp,:)-pTSrw_2016(26+pp,:),'color',wedge_colors(pp,:),'linewidth',5); hold on   % normally 'linewidth',7 
  plot(pTSrw_2096(31+pp,:)-pTSrw_2016(31+pp,:),'color',wedge_colors(pp,:),'linewidth',3,'LineStyle',':'); hold on %cmap(2*pp,:) % normally 'linewidth',5
  % compute Residual    Grazing                                     non grazing
  res = res - (pTSrw_2096(26+pp,:)-pTSrw_2016(26+pp,:)) - (pTSrw_2096(31+pp,:)-pTSrw_2016(31+pp,:));
  end
  plot(res,'color','k','linewidth',4); hold on

 %legend('T','VD','CC','P','Z')
  yline(0,'k:','linewidth',2)
  xlim([1,12]); 
  xticklabels(pmonL)
  ax = gca; ax.TickLabelInterpreter = 'latex'; ax.FontSize = ticsize;
  %title({'Partial Derivatives';'of Loss Terms'}, 'interpreter', 'latex','fontsize',35)
 
   % Add Background
  % Get the current axis limits
    ylim_values = ylim;  xlim_values = xlim;  % Get the current x-limits of the plot

  % Create a red patch for the area above 0
     patch([xlim_values(1), xlim_values(2), xlim_values(2), xlim_values(1)], ...  % X-coordinates
      [0, 0, ylim_values(2), ylim_values(2)], ...  % Y-coordinates for red area
      'b', 'FaceAlpha', 0.05, 'EdgeColor', 'none');  % Red color with transparency

  % Create a blue patch for the area below 0
     patch([xlim_values(1), xlim_values(2), xlim_values(2), xlim_values(1)], ...  % X-coordinates
      [ylim_values(1), ylim_values(1), 0, 0], ...  % Y-coordinates for blue area
      'r', 'FaceAlpha', 0.05, 'EdgeColor', 'none');  % Blue color with transparency

  % Bring the plot back on top of the patches
     uistack(findobj(gca, 'Type', 'line'), 'top');
  % adjust axis 
    ylim(ylim_values )

    title({'Loss Rate Decompositon:';'Seasonal Contribution'}, 'interpreter', 'latex','fontsize',tsize)
    ylabel('(d$^{-1}$)', 'interpreter', 'latex','fontsize',ysize )
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%% Loss RATES together %%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  elseif vlist(s) == 71
 

  % Add patches between the two lines
    x = 1:12; x_interp = 1:.1:12;
    for ii =[11,17]
    y1 = interp1(x, pTSrw_2016(ii,:), x_interp, 'linear');  % Linear interpolation for y1
    y2 = interp1(x, pTSrw_2096(ii,:), x_interp, 'linear');  % Linear interpolation for y2

    plot(x_interp,y1,'k:','linewidth',2); hold on
    plot(x_interp,y2,'k','linewidth',3); hold on

    
    % Store figure axes limits for arrow coordinate conversion
    ax = gca;
    xlim_ = xlim(ax); 
    ylim_ = ylim(ax); 
    axpos = get(ax, 'Position');  % [x0 y0 width height]

    % Add shaded patch
    for i = 1:length(x_interp)-1
        x_patch = [x_interp(i), x_interp(i+1), x_interp(i+1), x_interp(i)];
    
        if y1(i) > y2(i)
            y_patch = [y2(i), y2(i+1), y1(i+1), y1(i)];
            patch(x_patch, y_patch, 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
        elseif y2(i) > y1(i)
            y_patch = [y1(i), y1(i+1), y2(i+1), y2(i)];
            patch(x_patch, y_patch, 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
        end
    end

        % Add 12 custom arrows between lines with triangle heads
        for add_arrows=1
        x_arrow_pos = linspace(x_interp(1), x_interp(end), 12);
        
        arrow_head_height = range(ylim_)./20;  % in y-axis units
        arrow_head_width = .3; % in x-axis units
        
        ylims = ylim(gca);               % get current y-axis limits
        y_threshold = 0.08 * range(ylim_);  % 15% of y-range
        
        for j = 2:length(x_arrow_pos)
            [~, idx] = min(abs(x_interp - x_arrow_pos(j)));
        
            x0 = x_interp(idx);
            y1_val = y1(idx);
            y2_val = y2(idx);
            dy = y2_val - y1_val;
        
            if abs(dy) < y_threshold
                continue  % Skip if vertical difference is too small
            end
        
            % Arrow shaft: always from y1 to y2, minus arrowhead height
            shaft_top = y2_val - sign(dy)*arrow_head_height;
            
            % Arrowhead triangle
            x_head = [x0 - arrow_head_width/2, x0 + arrow_head_width/2, x0];
            if dy > 0  % Upward
                line([x0, x0], [y1_val, shaft_top], 'Color', 'b', 'LineWidth', 1);
        
                y_head = [y2_val - arrow_head_height, ...
                          y2_val - arrow_head_height, ...
                          y2_val];
                 fill(x_head, y_head, 'b', 'EdgeColor', 'none');
        
            else       % Downward
                line([x0, x0], [y1_val, shaft_top], 'Color', 'r', 'LineWidth', 1);
        
                y_head = [y2_val + arrow_head_height, ...
                          y2_val + arrow_head_height, ...
                          y2_val];
        
                 fill(x_head, y_head, 'r', 'EdgeColor', 'none');
            end
           
        end
        end

    end
      
    % Bring the lines back on top of the patches
      uistack(findobj(gca, 'Type', 'line'), 'top');
      % Adjust labels
      xlim([1,12]); %ylim([0,.33])
      xticklabels(pmonL)
      ax = gca; ax.TickLabelInterpreter = 'latex'; ax.FontSize = ticsize;
      pc = 100.*(mean(TSrw_2096(16,:))-mean(TSrw_2016(16,:)))./mean(TSrw_2016(16,:));
        text(ax,  .99,.99, [num2str(round(pc)), '$\%$'], 'units','normalized',...
                'HorizontalAlignment','right','VerticalAlignment','top', ...
                'FontSize',20,'Interpreter','latex','Color','k');
      title({'Phytotoplankton Specific';'Loss Rate'}, 'interpreter', 'latex','fontsize',tsize)
      ylabel(CoreVar_Inv_Units(6), 'interpreter', 'latex','fontsize',ysize )
     

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%% ALL OTHER VAIRABLES %%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  else

    x = 1:12; 
    x_interp = 1:.1:12;
    % Coexistance
    if vlist(s) ==101 
    y1_dum = interp1(x, pTSrw_2016(12,:), x_interp, 'linear');  
    y2_dum = interp1(x, pTSrw_2096(12,:), x_interp, 'linear'); 
    y1 = abs(2*y1_dum -100); y2 = abs(2*y2_dum -100);
    CoreVar_Inv_Units{vlist(s)} = '($\%$)'; CoreVar_Names{vlist(s)} = {'Coexistence';'(abs($f_{\mathrm{sp}}$ - $f_{\mathrm{diat}}$))'};
    % First Order Stability 
    elseif  vlist(s) == 102 
    % MLd mean
    %y1_dum = interp1(x, pTSrw_2016(1,:)./pTSrw_2016(15,:), x_interp, 'linear'); 
    %y2_dum = interp1(x, pTSrw_2096(1,:)./pTSrw_2096(15,:), x_interp, 'linear');
    y1_dum = interp1(x, pTSrw_2016(46,:), x_interp, 'linear'); 
    y2_dum = interp1(x, pTSrw_2096(46,:), x_interp, 'linear');
    
    % Dilution indec
    %y1_dum = interp1(x, pTSrw_2016(23,:), x_interp, 'linear'); 
    %y2_dum = interp1(x, pTSrw_2096(23,:), x_interp, 'linear');
    do=.005; % differential operator 
    y1_do  = [y1_dum-do/2; y1_dum+do/2 ];  y2_do  = [y2_dum-do/2; y2_dum+do/2 ];
    Tf    = 1.7.^((18 - 30)/10); Z=10;
      % Clearance rates 
      if model == 5;  % Marble
      clear1 = (3.3.*Tf.*(y1_do./(1.2+y1_do)))./y1_do;
      clear2 = (3.3.*Tf.*(y2_do./(1.2+y2_do)))./y2_do;
      elseif model == 1;  % iHammoc
      clear1 = (1.2.*(y1_do./(9.7+y1_do)))./y1_do;
      clear2 = (1.2.*(y2_do./(9.7+y2_do)))./y2_do;
      end
    y1 = ( clear1(2,:) - clear1(1,:) )./do; % First Order Stability
    y2 = ( clear2(2,:) - clear2(1,:) )./do; 
    CoreVar_Inv_Units{vlist(s)} = '(m$^6$ mmol$^{-2}$ d$^{-1}$) '; CoreVar_Names{vlist(s)} = {'Stabilizing Influence of';'Grazing Formulation'};
    pc = 100.*(mean(y2)-mean(y1))./abs(mean(y1));
   
    % [Phyto]MLD 
    elseif  vlist(s) == 103
    y1 = interp1(x, pTSrw_2016(1,:)./pTSrw_2016(15,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(1,:)./pTSrw_2096(15,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = 'log$_{10}$(mmol m$^{-3}$)'; CoreVar_Names{vlist(s)} = {'[Phytoplankton]$_{\mathrm{MLD}}$'};
    CoreVar_Inv_Units{vlist(s)} = '(mmol m$^{-3}$)'; CoreVar_Names{vlist(s)} = {'[Phytoplankton]$_{\mathrm{MLD}}$'};
    elseif  vlist(s) == 104 % Clearance Rate 
    y1 = interp1(x, pTSrw_2016(28,:)./pTSrw_2016(2,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(28,:)./pTSrw_2096(2,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = 'm$^3$ mmol$^{-1}$ d$^{-1}$'; CoreVar_Names{vlist(s)} = {'Clearance Rate';'(Grazing Efficiency)'};
    elseif  vlist(s) == 105 % VD - Clearance Rate 
    y1 = interp1(x, pTSrw_2016(4,:)./pTSrw_2016(2,:)./pTSrw_2016(1,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(4,:)./pTSrw_2096(2,:)./pTSrw_2096(1,:), x_interp, 'linear');
    %y1 = interp1(x, pTSrw_2016(7,:)./pTSrw_2016(23,:), x_interp, 'linear'); 
    %y2= interp1(x, pTSrw_2096(7,:)./pTSrw_2096(23,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = 'm$^2$ mmol$^{-1}$ d$^{-1}$'; CoreVar_Names{vlist(s)} = {'Clearance Rate'};
    elseif  vlist(s) == 119
    y1 = interp1(x, pTSrw_2016(9,:)./pTSrw_2016(2,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(9,:)./pTSrw_2096(2,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = '(d$^{-1}$)'; CoreVar_Names{vlist(s)} = {'Zooplankton Specific';'Mortality Rate'};
    pc = 100.*(mean(y1)-mean(y2))./mean(y1);
       
    elseif  vlist(s) == 120
    y1 = interp1(x, pTSrw_2016(28,:)./pTSrw_2016(2,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(28,:)./pTSrw_2016(2,:), x_interp, 'linear');
    %y1 = interp1(x, pTSrw_2016(7,:)./pTSrw_2016(23,:), x_interp, 'linear'); 
    %y2= interp1(x, pTSrw_2096(7,:)./pTSrw_2096(23,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = '(m$^{-3}$ mmol$^{-1}$ d$^{-1}$)'; CoreVar_Names{vlist(s)} = {'Clearance Index'}; % higher means more dilution  (effective thickness)
    
    elseif  vlist(s) == 49
    y1 = interp1(x, pTSrw_2016(29,:)./pTSrw_2016(2,:), x_interp, 'linear'); 
    y2= interp1(x, pTSrw_2096(29,:)./pTSrw_2016(2,:), x_interp, 'linear');
    CoreVar_Inv_Units{vlist(s)} = '(m$^{-3}$ mmol$^{-1}$ d$^{-1}$)'; CoreVar_Names{vlist(s)} = {'Clearance Rates (CC)'}; % higher means more dilution  (effective thickness)
   
    else
    y1 = interp1(x, pTSrw_2016(vlist(s),:), x_interp, 'linear');  
    y2 = interp1(x, pTSrw_2096(vlist(s),:), x_interp, 'linear');  
    pc = 100.*(mean(y2)-mean(y1))./mean(y1);
    % pc = 100.*(mean(y2(27:65))-mean(y1(27:65)))./abs(mean(y1(27:65)));
    end

    plot(x_interp, y1, 'k:', 'linewidth', 2); hold on
    plot(x_interp, y2, 'k', 'linewidth', 3); hold on
    
    % Add Arrows and Patches
        ax = gca;
        isLog = false; %(vlist(s) == 103);
        if isLog, set(ax, 'YScale', 'log'); end
        
        % Helpers to move between data-space and axis-space
        if isLog
            y2t  = @(y) log10(y);      % data -> transformed (log axis)
            t2y  = @(t) 10.^t;         % transformed -> data
        else
            y2t  = @(y) y;
            t2y  = @(t) t;
        end
        
        % Axis ranges in transformed space
        yl      = ylim(ax);
        yl_t    = y2t(yl);
        yrangeT = range(yl_t);
        
        % Shaded patch between lines (skip non-positive values on log axes)
        for i = 1:length(x_interp)-1
            x_patch = [x_interp(i), x_interp(i+1), x_interp(i+1), x_interp(i)];
        
            y1p = [y1(i), y1(i+1)];
            y2p = [y2(i), y2(i+1)];
        
            if isLog
                if any(y1p <= 0) || any(y2p <= 0)
                    continue;  % can't draw in log space if <= 0
                end
            end
        
            if y1(i) > y2(i)
                y_patch = [y2(i), y2(i+1), y1(i+1), y1(i)];
                patch(x_patch, y_patch, 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
            elseif y2(i) > y1(i)
                y_patch = [y1(i), y1(i+1), y2(i+1), y2(i)];
                patch(x_patch, y_patch, 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
            end
        end
        
        % --- Custom arrows between lines (consistent size on linear or log axes) ---
        x_arrow_pos       = linspace(x_interp(1), x_interp(end), 12);
        arrow_head_heightT = 0.05 * yrangeT;   % head height in transformed units (5% of axis)
        arrow_head_width   = 0.3;              % in x-axis units
        y_thresholdT       = 0.08 * yrangeT;   % threshold in transformed units (8% of axis)
        
        for j = 2:length(x_arrow_pos)
            % locate nearest sample on the interpolated grid
            [~, idx] = min(abs(x_interp - x_arrow_pos(j)));
        
            x0     = x_interp(idx);
            y1_val = y1(idx);
            y2_val = y2(idx);
        
            % skip invalid values in log space
            if isLog && (y1_val <= 0 || y2_val <= 0), continue; end
        
            % compute distance in transformed coordinates
            y1T  = y2t(y1_val);
            y2T  = y2t(y2_val);
            dyT  = y2T - y1T;
        
            if abs(dyT) < y_thresholdT
                continue  % difference too small to draw an arrow
            end
        
            % Shaft end (just under the head) in transformed coords
            shaftTopT = y2T - sign(dyT)*arrow_head_heightT;
            shaftTop  = t2y(shaftTopT);
        
            % Draw shaft (in data space)
            if dyT > 0
                line([x0, x0], [y1_val, shaftTop], 'Color', 'b', 'LineWidth', 1);
            else
                line([x0, x0], [y1_val, shaftTop], 'Color', 'r', 'LineWidth', 1);
            end
        
            % Arrowhead triangle (convert the head base from transformed -> data space)
            x_head = [x0 - arrow_head_width/2, x0 + arrow_head_width/2, x0];
            if dyT > 0  % Upward arrow
                baseY = t2y(y2T - arrow_head_heightT);
                tipY  = y2_val;
                y_head = [baseY, baseY, tipY];
                fill(x_head, y_head, 'b', 'EdgeColor', 'none');
            else        % Downward arrow
                baseY = t2y(y2T + arrow_head_heightT);
                tipY  = y2_val;
                y_head = [baseY, baseY, tipY];
                fill(x_head, y_head, 'r', 'EdgeColor', 'none');
            end
        end

      
      % Bring the lines back on top
      uistack(findobj(gca, 'Type', 'line'), 'top');
      % Change lable
        text(ax,  .99,.99, [num2str(round(pc)), '$\%$'], 'units','normalized',...
                'HorizontalAlignment','right','VerticalAlignment','top', ...
                'FontSize',20,'Interpreter','latex','Color','k');
      xlim([1,12])
      
      xticklabels(pmonL)
      ax = gca; ax.TickLabelInterpreter = 'latex'; ax.FontSize = ticsize;
      title(CoreVar_Names{vlist(s)}, 'interpreter', 'latex','fontsize',tsize)
      ylabel(CoreVar_Inv_Units(vlist(s)), 'interpreter', 'latex','fontsize',ysize )
      
      % Subplot label
      labs = {'i','ii','iii','iv','v','vi','vii','viii','ix','x'};
      

    
  end
       
        % Subplot label
        text(ax,  .02,.99, labs(s), 'units','normalized',...
                'HorizontalAlignment','left','VerticalAlignment','top', ...
                'FontSize',20,'Interpreter','latex','Color','k');


        % DEfine Axis for boarders 
        yl = ylim(ax); 
          % Special limits
          if vlist(s)==72 ; yl=[22,63];end %yl =[30,90]; end
           %if s==9 ; yl=[65,130]; end 
        xl = xlim(ax);

        % Draw a frame that matches the current axes limits
        if ~isnan(bor(s)) && bor(s) >5
        ax = gca; box(ax,'off'); ax.Color = 'none';
        ax.Layer = 'top';                 % keep data/grid beneath the frame
        rectangle('Position',[xl(1) yl(1) diff(xl) diff(yl)], ...
                  'EdgeColor',wedge_colors(bor(s)-5,:), 'LineWidth',3,'linestyle','--');
        ylim(yl)
        elseif ~isnan(bor(s)) && bor(s) <6
        ax = gca; box(ax,'off'); ax.Color = 'none';
        ax.Layer = 'top';                 % keep data/grid beneath the frame
        rectangle('Position',[xl(1) yl(1) diff(xl) diff(yl)], ...
                  'EdgeColor',wedge_colors(bor(s),:), 'LineWidth',3,'linestyle','-');
        ylim(yl)
        end

      if vlist(s)==15 ||vlist(s)==22 ||vlist(s)==50 ||vlist(s)==101 || vlist(s)==42 ; ax.YDir = 'reverse'; end 
      if vlist(s)==101 ; ylim([0,100]) ;
      elseif vlist(s)==72 ; ylim([22,63]); %ylim([30,90]) 
      end 
      % if s==3 || s ==4 ; ylim([0,100]); end 
      % if s==5 || s ==6 ; ylim([0,.3]); end    
      %if s==9 ; ylim([65,130]); end 
 end 
end

set(gcf,'HandleVisibility','callback'); 






%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Plot Map Region
%%%%%%%%%%%%%%%%%%%%%%%%%%

if map
for plotreg = 1
h=figure('Units','normalized','Position',[.5 .5 .6 .6]); clf
set(gcf,'color','w'); set(gca,'Visible','off')
map_proj = 'robinson'; clims = [-60,60];
%map_proj = 'Mollweide'
    % === Step 1. Data & masking (unchanged) ===
    plotvar = squeeze(CoreVar_Inv_pc(1,:,:));     % or *_rw
    % Dil_ind_2015  = mean( squeeze(CoreVar_Inv_adj(1,1:60,:,:)) ./ squeeze(CoreVar_Inv_adj(23,1:60,:,:))  ); 
    %      Dil_ind_2100  = mean( squeeze(CoreVar_Inv_adj(1,61:end,:,:)) ./ squeeze(CoreVar_Inv_adj(23,61:end,:,:))  ); 
    %      plotvar       = squeeze( Dil_ind_2100 - Dil_ind_2015);
    % 
    Lon_p   = wrapTo360(Lon);                     % make sure grid is [0,360)
    plotvar_masked = plotvar;
    plotvar_masked(~FullMask) = NaN;
    %plotvar_masked(~region_mask) = NaN;  % Only mask Region

    
    % --- 1D lat/lon vectors for the reference object (now 0–360) ---
    latvec = linspace(-89.5,  89.5, 180);
    lonvec = linspace(  0.5, 359.5, 360);

    % === Step 2. Region longitudes on the circle ===
    Lon_region = wrapTo360(Lon_p(region_mask));
    Lat_region =      Lat(region_mask);

    % circular center and minimal-span window
    [map_center, lonmin, lonmax, span, lonlims] = regionLonWindow(Lon_region);

    % === Step 3. Geo ref (full global raster; you’ll mask to show region) ===
    R = georefcells([min(latvec), max(latvec)], [0 360], size(plotvar));

    % === Step 4. Bounding box in latitude  ===  
         % latmin latmax lomin lonmax
    offsets = [ 0,0,0,0; ...
                0,2,0,0;...
              -1,-6,15,1; ... % Trades Atlantic
                0,0,0,0;...
                -5,5,5,0; ...
               -2,2,-2,2;... % NA
                0,0,0,0; ...
                0,0,0,0;...
                0,0,0,0; ...
                0,-5,0,0;... % Southern Ocean
                +10,5,4,-4]; % North Pacific
    
    latmin = min(Lat_region) +offsets(region_index,1); 
    latmax = max(Lat_region) +offsets(region_index,2); 
    lonmin = lonmin          +offsets(region_index,3); 
    lonmax = lonmax          +offsets(region_index,4); 
    map_center = mod((lonmin + lonmax) / 2, 360);

    
    % “circumpolar” heuristic; 
    is_circumpolar = (abs(latmin) > 55 && span > 300) || region_index == 1;


    % ---- bottom - up layer axis ----
    h1 = axes('Position',[0.1 0.2 0.75 0.8]); axis(h1,'off')
    if is_circumpolar
        if mean(Lat_region) < 0
            axesm('stereo','Origin',[-90 0 0], ...
                  'MapLatLimit',[-90 latmax], 'MapLonLimit', lonlims);
        else
            axesm('stereo','Origin',[ 90 0 0], ...
                  'MapLatLimit',[ 40  90],   'MapLonLimit', lonlims);
        end
    else
        % Global projection; center the world at the circular center
        axesm(map_proj,'Origin',[0 map_center 0], ...
              'MapLatLimit',[latmin latmax],'MapLonLimit',[lonmin lonmax]);
    end
    axis off; tightmap

    gi1 = geoshow(plotvar, R, 'DisplayType','texturemap');
    gi1.AlphaDataMapping = 'none'; gi1.FaceAlpha = 'texturemap';
    combined_mask = ~isnan(BUmask) & ~isnan(plotvar); 
    combined_mask =  ~isnan(plotvar_masked); % To Plot Bottom up layer: Comment Out This + and Comemnt Out TD section 
    gi1.AlphaData = double(combined_mask);
    colormap(h1, cmapBU); caxis(h1, [clims]);
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7]); plotm(coastlat,coastlon,'k');

    % ---- top - down layer axis  ----
    h2 = axes('Position',get(h1,'Position'),'Color','none'); axis(h2,'off')
    if is_circumpolar
        if mean(Lat_region) < 0
            axesm('stereo','Origin',[-90 0 0], ...
                  'MapLatLimit',[-90 latmax], 'MapLonLimit', lonlims);
        else
            axesm('stereo','Origin',[ 90 0 0], ...
                  'MapLatLimit',[ 40  90],   'MapLonLimit', lonlims);
        end
    else
        axesm(map_proj,'Origin',[0 map_center 0], ...
              'MapLatLimit',[latmin latmax],'MapLonLimit',[lonmin lonmax]);
    end
    axis off; tightmap; setm(gca,'Frame','on','Grid','off')

    if TD    
    gi2 = geoshow(plotvar, R, 'DisplayType','texturemap');
    gi2.AlphaDataMapping = 'none'; gi2.FaceAlpha = 'texturemap';
    combined_mask = ~isnan(TDmask) & ~isnan(plotvar); 
    combined_mask =  ~isnan(plotvar_masked); 
    %alpha(gi2, double(combined_mask));
    %alpha(gi2, double(~isnan(plotvar_masked)));
    gi2.AlphaData = double(combined_mask);  
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7]); plotm(coastlat,coastlon,'k');
    colormap(h2, cmapTD); caxis(h2, [clims(1) clims(2)]);
    end

    % % ---- contours (make longitudes consistent with the chosen center) ----
    % % If you see contours “break” at the seam, wrap relative to map_center:
    Lon_for_contour = wrapTo180(Lon_p - map_center) + map_center; % optional
    %contourm(Lat_t, Lon_for_contour, double(region_mask), [1 1], 'k','LineWidth',1);
    contourm(Lat_t, Lon_t, double(FullMask),   [1 1], 'k','LineWidth',3);

    % ---- thick colorbar strip (unchanged) ----
    if is_circumpolar
        h_cb = axes('Position',[0.2 0.06 0.55 0.08]);
    else
        h_cb = axes('Position',[0.1 0.06 0.75 0.08]);
    end
    set(h_cb,'Visible','off')
    imagesc(linspace(clims(1), clims(2),256)); colormap(h_cb, cb_map); caxis(h_cb,[clims(1), clims(2)]);
    set(h_cb,'YTick',[],'XTick',linspace(1,256,9),'XTickLabel',linspace(clims(1), clims(2),9), ...
             'TickLabelInterpreter','latex','FontSize',30);
    %title('Change in biomass (\%)','Interpreter','latex','FontSize',40)
end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Include Profiles %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if profiles
for make_profiles=[1]    

% --- Interpolate FullMask (on Lat_t/Lon_t grid) to native grid (Lat_n/Lon_n)
% Inputs available:
%   FullMask : nLat_t x nLon_t logical/double (0/1)
%   Lat_t, Lon_t : nLat_t x nLon_t 2-D grids for the mask
%   Lat, Lon from the profile file (saved as Lat_n, Lon_n below)

cd /Users/twrohr/Desktop/data.nosync/SSP585/Diag_Profs
mods = dir; mods = mods(~[mods.isdir]);
load(mods(mod_order(model)).name, 'ppt','zpt','ppt_fr','zpt_fr','Lat','Lon','depthlev');

% Keep a copy of the native grid (don’t overwrite these later)
Lat_n = Lat; 
Lon_n = Lon;

% ---- Build coordinate vectors from the mask grid (Lat_t/Lon_t are 2-D) ----
% Convert longitudes to 0..360 on BOTH grids
Lon_t360 = mod(Lon_t, 360);
Lon_n360 = mod(Lon_n, 360);

% Extract 1-D coordinate vectors from the 2-D mask grid
latv = Lat_t(:,1);
lonv = Lon_t360(1,:);

% Ensure latitude is ascending for griddedInterpolant
flip_lat = false;
if any(diff(latv) < 0)
    latv = flipud(latv);
    FullMask = flipud(FullMask);
    flip_lat = true;
end

% Ensure longitude is strictly ascending; if not, sort and permute mask
if any(diff(lonv) <= 0)
    [lonv, j] = sort(lonv);
    FullMask  = FullMask(:, j);
end

% Build interpolant on the *gridded* domain (lat, lon); nearest keeps mask crisp
F_mask = griddedInterpolant({latv, lonv}, double(FullMask), 'nearest', 'nearest');

% Make sure the native grid we evaluate on is 2-D (lat x lon)
if isvector(Lat_n) && isvector(Lon_n) && numel(Lat_n) > 1 && numel(Lon_n) > 1
    % Common case: 1-D vectors for coords
    [LonN2, LatN2] = meshgrid(Lon_n360, Lat_n);
else
    % Already 2-D arrays
    LatN2 = Lat_n;
    LonN2 = Lon_n360;
end

% Evaluate the mask and make it logical
InterpMask = F_mask(LatN2, LonN2);
InterpMask = isfinite(InterpMask) & (InterpMask >= 0.5);

% -------------------------------------------------------------------------
% Depth-wise regional means over the interpolated mask
% (kept your original model-dependent slicing)
% -------------------------------------------------------------------------
clearvars P_2015 P_2100 Z_2015 Z_2100

if model < 5 % Ignores Month (mon) option and plots 2015-2020 January Clim along with 2100 Interannua January 
    % ppt:    (syn_Prof,            mon,          depth, lat, lon)
    %          1: 2015:2020 clim    1: 1st Jan 
    % ppt_fr: (mon,                              depth, lat, lon)
    %          1:5 first 5 jan, 6: last jan    
    [~, ndepth, ~, ~] = size(zpt_fr);
    P_2015_n = squeeze(ppt(1, 1, :, :, :));      % Climatology Jan of 2015-2020
    P_2100_n = squeeze(ppt_fr(6, :, :, :));      % 2100 (last) Free running Jan
    P_2100_n_vd = squeeze(ppt(2, 2, :, :, :));   % VD syn Jan of 2100
    Z_2015_n = squeeze(zpt(1, 1, :, :, :));      % Climatology Jan of 2015-2020
    Z_2100_n = squeeze(zpt_fr(6,:, :, :));       % 2100 (last) Free running Jan
    Z_2100_n_vd = squeeze(zpt(2, 2, :, :, :));  % VD syn Jan of 2100
    month_tit = 'January';

elseif model == 5
    % ppt:    (Pgroup, syn_prof, mon, depth, lat, lon), zpt_fr: (mon, depth, lat, lon); mon:1 first january, July ; last januaray july
    % ppt_fr: (Pgroup, mon, depth, lat, lon), zpt_fr: (mon, depth, lat, lon)
    [~, ndepth, ~, ~] = size(zpt_fr);
    P_2015_n = squeeze(sum(ppt(:, 1, mon,  :, :, :), 1));       % Sum P of Climatology Jan/July of 2015-2020
    P_2100_n = squeeze(sum(ppt_fr(:, mon+2, :, :, :), 1));      % Sum P 2100 (last) Free running Jan or July
    P_2100_n_vd = squeeze(sum(ppt(:, 2, mon+2,  :, :, :), 1));  % Sum P of VD Syn Jan/July 2100
    Z_2015_n = squeeze(zpt(1, mon,   :, :, :));                 % Climatology Jan/July of 2015-2020
    Z_2100_n = squeeze(zpt_fr(mon+2, :, :, :));                 % 2100 (last) Free running Jan or July
    Z_2100_n_vd = squeeze(zpt(2, mon+2,   :, :, :));            % Z of VD Syn Jan/July 2100
    if mon==1; month_tit = 'January';  
    elseif mon==2; month_tit = 'July';
    end

else
    % ppt_fr: (Pgroup, mon, depth, lat, lon), zpt_fr: (Zgroup?, mon, depth, lat, lon) or similar
    [~, ~, ndepth, ~, ~] = size(ppt_fr);
    P_2015_n = squeeze(sum(ppt(:, 1, mon,   :, :, :), 1));      % Sum P of Climatology Jan/July of 2015-2020
    P_2100_n = squeeze(sum(ppt_fr(:, mon+2, :, :, :), 1));      % Sum P 2100 (last) Free running Jan or July
    P_2100_n_vd = squeeze(sum(ppt(:, 2, mon+2,  :, :, :), 1));  % Sum P of VD Syn Jan/July 2100
    Z_2015_n = squeeze(sum(zpt(:, 1,  mon,   :, :, :), 1));     % Sum P of Climatology Jan/July of 2015-2020
    Z_2100_n = squeeze(sum(zpt_fr(:, mon+2, :, :, :), 1));      % Sum Z 2100 (last) Free running Jan or July
    Z_2100_n_vd = squeeze(sum(zpt(:, 2, mon+2,  :, :, :), 1));  % Sum Z of VD Syn Jan/July 2100
    if mon==1; month_tit = 'January';
    elseif mon==2; month_tit = 'July';
    end
end

% Sanity check: mask size must match spatial slice size
spatial_sz = size(P_2015_n); spatial_sz = spatial_sz(end-1:end);
if ~isequal(size(InterpMask), spatial_sz)
    error('InterpMask size %s does not match field size %s.', mat2str(size(InterpMask)), mat2str(spatial_sz));
end

P_2015 = nan(ndepth,1); P_2100 = nan(ndepth,1); P_2100_vd = nan(ndepth,1);
Z_2015 = nan(ndepth,1); Z_2100 = nan(ndepth,1); Z_2100_vd = nan(ndepth,1);

for k = 1:ndepth
    % Phyto 2015
    layer = squeeze(P_2015_n(k,:,:));
    layer(~InterpMask) = NaN;
    P_2015(k) = mean(layer(:), 'omitnan');

    % Phyto 2100
    layer = squeeze(P_2100_n(k,:,:));
    layer(~InterpMask) = NaN;
    P_2100(k) = mean(layer(:), 'omitnan');

    % Phyto 2100 VD syn
    layer = squeeze(P_2100_n_vd(k,:,:));
    layer(~InterpMask) = NaN;
    P_2100_vd(k) = mean(layer(:), 'omitnan');

    % Zoo 2015
    layer = squeeze(Z_2015_n(k,:,:));
    layer(~InterpMask) = NaN;
    Z_2015(k) = mean(layer(:), 'omitnan');

    % Zoo 2100
    layer = squeeze(Z_2100_n(k,:,:));
    layer(~InterpMask) = NaN;
    Z_2100(k) = mean(layer(:), 'omitnan');

    % Zoo 2100 VD syn
    layer = squeeze(Z_2100_n_vd(k,:,:));
    layer(~InterpMask) = NaN;
    Z_2100_vd(k) = mean(layer(:), 'omitnan');

end


%%%%%%%%%%%%%%%%%%%%%%%
%%%% Figure
figure('Units','normalized','Position',[2 0 1.75 1]); %clc 8 panel (.5,2
   set(gcf,'color','w');
   
  t = tiledlayout(2, 5); %,'Padding', 'compact'
  for s =1:10
      nexttile

        if s ==1
        
        plot(P_2015,depthlev(1:ndepth),'k:','linewidth',4); hold on; % P
        %plot(P_2100,depthlev(1:ndepth),'k','linewidth',4); hold on; % P
        plot(P_2100_vd,depthlev(1:ndepth),'k','linewidth',2); hold on; % P
        set(gca, 'YDir', 'reverse');
        %set(gca, 'XScale', 'log'); 
        set(gca, 'FontSize', ticsize,'TickLabelInterpreter','latex') % Set font size for ticks
        ylabel('Depth (m)','fontsize',ysize,'interpreter','latex'); ylim([0,300]);  
        %xlabel(' (mmolC m$^{-3}$)','fontsize',14,'interpreter','latex')
        title('[Phytoplankton]','fontsize',tsize,'interpreter','latex');
        
        ylim([0,150])
        
     % --- Shade differences between profiles & draw horizontal arrows ---
        ax = gca;
        
        % Depth vector (make it a column)
        y_prof = depthlev(1:ndepth);
        y_prof = y_prof(:);
        
        % Profiles as column vectors too
        x2015  = P_2015(:);
        x2100  = P_2100_vd(:);
        
        % A fine depth grid for smooth patches
        yq = linspace(min(y_prof,[],'omitnan'), max(y_prof,[],'omitnan'), 400);
        x2015q = interp1(y_prof, x2015, yq, 'linear', 'extrap');
        x2100q = interp1(y_prof, x2100, yq, 'linear', 'extrap');
        
        % Clip to current y-limits so we don't draw outside the panel
        yl = ylim(ax);
        ymin = min(yl); ymax = max(yl);
        % With 'YDir' reversed this is fine—these are still data coordinates
        maskY = (yq >= min(ymin,ymax)) & (yq <= max(ymin,ymax));
        yq     = yq(maskY);
        x2015q = x2015q(maskY);
        x2100q = x2100q(maskY);
        
        % Draw short horizontal quads between adjacent depth samples
        for i = 1:numel(yq)-1
            x1a = x2015q(i);   x1b = x2015q(i+1);
            x2a = x2100q(i);   x2b = x2100q(i+1);
            ya  = yq(i);       yb  = yq(i+1);
        
            if any(~isfinite([x1a x1b x2a x2b ya yb]))
                continue
            end
        
            % Decide color by midpoint difference at this slab
            dx_mid = 0.5*((x2a+x2b) - (x1a+x1b));
            if dx_mid > 0
                col = 'b';   % blue (increase)
            elseif dx_mid < 0
                col = 'r';   % red (decrease)
            else
                continue  % exactly equal → nothing to shade
            end
        
            patch([x1a x1b x2b x2a], [ya yb yb ya], col, ...
                  'FaceAlpha', 0.18, 'EdgeColor', 'none', 'Parent', ax);
        end
        
        % Bring the two profile lines back on top
        uistack(findobj(ax,'Type','line'), 'top');
        
        % ---- Horizontal arrows at selected depths ----
        % Choose 10 evenly spaced depths within the visible range
        n_arrows = 10;
        y_ar = linspace(min(yq), max(yq), n_arrows);
        
        % Sizes (in data units)
        xl = xlim(ax);
        xr = range(xl);
        yr = range(ylim(ax));
        arrow_head_len   = 0.05 * xr;      % head length in x-units
        arrow_head_halfh = 0.015 * yr;     % head half-height in y-units
        dx_threshold     = 0.06 * xr;      % minimum change to draw an arrow
        
        for k = 1:numel(y_ar)
            y0 = y_ar(k);
            x1 = interp1(y_prof, x2015, y0, 'linear', 'extrap');
            x2 = interp1(y_prof, x2100, y0, 'linear', 'extrap');
        
            if ~isfinite(x1) || ~isfinite(x2)
                continue
            end
        
            dx = x2 - x1;
            if abs(dx) < dx_threshold
                continue
            end
        
            if dx > 0
                % Rightward (increase) – blue
                col = 'b';
        
                x_start = x1;
                x_end   = x2;
                x_shaft_end = x_end - arrow_head_len;
        
                % Shaft
                line([x_start x_shaft_end], [y0 y0], 'Color', col, 'LineWidth', 1.1, 'Parent', ax);
        
                % Head (triangle)
                fill([x_end-arrow_head_len, x_end-arrow_head_len, x_end], ...
                     [y0-arrow_head_halfh,  y0+arrow_head_halfh,  y0], ...
                     col, 'EdgeColor', 'none', 'Parent', ax);
        
            else
                % Leftward (decrease) – red
                col = 'r';
        
                x_start = x1;
                x_end   = x2;
                x_shaft_end = x_end + arrow_head_len;  % moving left, head points left
        
                % Shaft
                line([x_start x_shaft_end], [y0 y0], 'Color', col, 'LineWidth', 1.1, 'Parent', ax);
        
                % Head (triangle)
                fill([x_end+arrow_head_len, x_end+arrow_head_len, x_end], ...
                     [y0-arrow_head_halfh,  y0+arrow_head_halfh,  y0], ...
                     col, 'EdgeColor', 'none', 'Parent', ax);
            end
        end
        %xline(1.33,'k:','linewidth',3)


    

        elseif s==2
        plot(Z_2015,depthlev(1:ndepth),'k:','linewidth',4); hold on; % Zoo
        %plot(Z_2100,depthlev(1:ndepth),'b','linewidth',4); hold on; 
        plot(Z_2100_vd,depthlev(1:ndepth),'k','linewidth',2); hold on; % Zoo
        set(gca, 'YDir', 'reverse'); 
        %set(gca, 'XScale', 'log'); 
        %ylabel('Depth (m)','fontsize',25,'interpreter','latex'); ylim([0,300]);  xlabel('Zooplankton Biomass (mmolC m$^{-3}$)','fontsize',25,'interpreter','latex')
        title('[Zooplankton]','fontsize',tsize,'interpreter','latex');
        
        %title(month_tit,'fontsize',25,'interpreter','latex');
        ylim([0,150])
        

         % --- Shade differences between profiles & draw horizontal arrows ---
        ax = gca;
        
        % Depth vector (make it a column)
        y_prof = depthlev(1:ndepth);
        y_prof = y_prof(:);
        
        % Profiles as column vectors too
        x2015  = Z_2015(:);
        x2100  = Z_2100_vd(:);
        
        % A fine depth grid for smooth patches
        yq = linspace(min(y_prof,[],'omitnan'), max(y_prof,[],'omitnan'), 400);
        x2015q = interp1(y_prof, x2015, yq, 'linear', 'extrap');
        x2100q = interp1(y_prof, x2100, yq, 'linear', 'extrap');
        
        % Clip to current y-limits so we don't draw outside the panel
        yl = ylim(ax);
        ymin = min(yl); ymax = max(yl);
        % With 'YDir' reversed this is fine—these are still data coordinates
        maskY = (yq >= min(ymin,ymax)) & (yq <= max(ymin,ymax));
        yq     = yq(maskY);
        x2015q = x2015q(maskY);
        x2100q = x2100q(maskY);
        
        % Draw short horizontal quads between adjacent depth samples
        for i = 1:numel(yq)-1
            x1a = x2015q(i);   x1b = x2015q(i+1);
            x2a = x2100q(i);   x2b = x2100q(i+1);
            ya  = yq(i);       yb  = yq(i+1);
        
            if any(~isfinite([x1a x1b x2a x2b ya yb]))
                continue
            end
        
            % Decide color by midpoint difference at this slab
            dx_mid = 0.5*((x2a+x2b) - (x1a+x1b));
            if dx_mid > 0
                col = 'b';   % blue (increase)
            elseif dx_mid < 0
                col = 'r';   % red (decrease)
            else
                continue  % exactly equal → nothing to shade
            end
        
            patch([x1a x1b x2b x2a], [ya yb yb ya], col, ...
                  'FaceAlpha', 0.18, 'EdgeColor', 'none', 'Parent', ax);
        end
        
        % Bring the two profile lines back on top
        uistack(findobj(ax,'Type','line'), 'top');
        
        % ---- Horizontal arrows at selected depths ----
        % Choose 10 evenly spaced depths within the visible range
        n_arrows = 10;
        y_ar = linspace(min(yq), max(yq), n_arrows);
        
        % Sizes (in data units)
        xl = xlim(ax);
        xr = range(xl);
        yr = range(ylim(ax));
        arrow_head_len   = 0.05 * xr;      % head length in x-units
        arrow_head_halfh = 0.015 * yr;     % head half-height in y-units
        dx_threshold     = 0.06 * xr;      % minimum change to draw an arrow
        
        for k = 1:numel(y_ar)
            y0 = y_ar(k);
            x1 = interp1(y_prof, x2015, y0, 'linear', 'extrap');
            x2 = interp1(y_prof, x2100, y0, 'linear', 'extrap');
        
            if ~isfinite(x1) || ~isfinite(x2)
                continue
            end
        
            dx = x2 - x1;
            if abs(dx) < dx_threshold
                continue
            end
        
            if dx > 0
                % Rightward (increase) – blue
                col = 'b';
        
                x_start = x1;
                x_end   = x2;
                x_shaft_end = x_end - arrow_head_len;
        
                % Shaft
                line([x_start x_shaft_end], [y0 y0], 'Color', col, 'LineWidth', 1.1, 'Parent', ax);
        
                % Head (triangle)
                fill([x_end-arrow_head_len, x_end-arrow_head_len, x_end], ...
                     [y0-arrow_head_halfh,  y0+arrow_head_halfh,  y0], ...
                     col, 'EdgeColor', 'none', 'Parent', ax);
        
            else
                % Leftward (decrease) – red
                col = 'r';
        
                x_start = x1;
                x_end   = x2;
                x_shaft_end = x_end + arrow_head_len;  % moving left, head points left
        
                % Shaft
                line([x_start x_shaft_end], [y0 y0], 'Color', col, 'LineWidth', 1.1, 'Parent', ax);
        
                % Head (triangle)
                fill([x_end+arrow_head_len, x_end+arrow_head_len, x_end], ...
                     [y0-arrow_head_halfh,  y0+arrow_head_halfh,  y0], ...
                     col, 'EdgeColor', 'none', 'Parent', ax);
            end
        end
        end
 end


end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end











%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Functions %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------- helper: circular window for longitudes ----------
function [center, lonmin, lonmax, span, lonlims] = regionLonWindow(lon_deg)
    lon_deg = mod(lon_deg, 360);             % [0,360)
    ang = deg2rad(lon_deg);
    c = mean(exp(1i*ang));                   % complex mean on the circle
    center = mod(rad2deg(angle(c)), 360);    % circular center in [0,360)

    rel = mod(lon_deg - center + 180, 360) - 180;  % relative in [-180,180]
    lo = min(rel); hi = max(rel);
    span = hi - lo;                           % minimal arc length (≤360)
    lonmin = mod(center + lo, 360);
    lonmax = mod(center + hi, 360);

    % For MapLonLimit, allow crossing by adding 360 to the upper bound
    if lonmax < lonmin
        lonlims = [lonmin, lonmax + 360];
    else
        lonlims = [lonmin, lonmax];
    end
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.40, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.60, 'linear','extrap');
end



%% compute profiles
clc

Tf    = 1.7.^((18 - 30)/10); Z=1;
g1 = (3.3.*Tf.*(P_2015./(1.2+P_2015)));
g2 = (3.3.*Tf.*(P_2100_vd./(1.2+P_2100_vd)));

GSP1 = sum(g1.*Z_2015*10)
GSP2 = sum(g2.*Z_2100_vd*10)

Pzoo1 = sum(P_2015.*Z_2015 *10 )./sum(Z_2015.*10)
Pzoo2 = sum(P_2100_vd.*Z_2100_vd *10 )./sum(Z_2100_vd.*10)

DilZ1 = sum(P_2015.*10)./Pzoo1
DilZ2 = sum(P_2100_vd.*10)./Pzoo2

l_vd1 = GSP1./sum(P_2015.*10)
l_vd2 = GSP2./sum(P_2100_vd.*10)

cl_vd1 = l_vd1/sum(Z_2015.*10)
cl_vd2 = l_vd2/sum(Z_2100_vd.*10)


clexp_vd1 = g1./P_2015; 
clexp_vd2 = g2./P_2100_vd; 

CL_vd1 = sum(clexp_vd1.*Z_2015 *10)./sum(Z_2015.*10)
CL_vd2 = sum(clexp_vd2.*Z_2100_vd*10 )./sum(Z_2100_vd.*10)

% CLb_1 = (3.3.*Tf.*(Pzoo1./(1.2+Pzoo1)));
% CLb_2 = (3.3.*Tf.*(Pzoo2./(1.2+Pzoo2)));

%%
dz = 10;
Tf = 1.7.^((18 - 30)/10);

% Local predator-specific ingestion (1/d)
g1 = 3.3 .* Tf .* (P_2015    ./ (1.2 + P_2015));
g2 = 3.3 .* Tf .* (P_2100_vd ./ (1.2 + P_2100_vd));
g3 = 3.3 .* Tf .* (P_2100 ./ (1.2 + P_2100));

% Depth-integrated grazing loss (mmol m^-2 d^-1)
GSP1 = sum(g1 .* Z_2015    .* dz);
GSP2 = sum(g2 .* Z_2100_vd .* dz);
GSP3 = sum(g3 .* Z_2100 .* dz);

% Column totals (mmol m^-2)
SigmaP1 = sum(P_2015    .* dz);
SigmaP2 = sum(P_2100_vd .* dz);
SigmaP3 = sum(P_2100 .* dz);
SigmaZ1 = sum(Z_2015    .* dz);
SigmaZ2 = sum(Z_2100_vd .* dz);
SigmaZ3 = sum(Z_2100 .* dz);

% Z-weighted mean P (mmol m^-3) = Phy_CZ
PhyCZ1 = sum(P_2015    .* Z_2015    .* dz) ./ SigmaZ1;
PhyCZ2 = sum(P_2100_vd .* Z_2100_vd .* dz) ./ SigmaZ2;
PhyCZ3 = sum(P_2100 .* Z_2100 .* dz) ./ SigmaZ3;

% Dilution Index (m) = SigmaP / Phy_CZ
DilZ1 = SigmaP1 / PhyCZ1;
DilZ2 = SigmaP2 / PhyCZ2;
DilZ3 = SigmaP3 / PhyCZ3;

% Per-prey specific loss (1/d): l_vd = gp = GSP / SigmaP
l_vd1 = GSP1 / SigmaP1;
l_vd2 = GSP2 / SigmaP2;

% Mean predator-specific ingestion (1/d): g_bar = GSP / SigmaZ
gbar1 = GSP1 / SigmaZ1;
gbar2 = GSP2 / SigmaZ2;
gbar3 = GSP3 / SigmaZ3;


% Clearance metrics
%  A) Z-weighted mean clearance (per predator): <g/P>_Z
epsP = 1e-12;
clexp1 = g1 ./ max(P_2015,    epsP);
clexp2 = g2 ./ max(P_2100_vd, epsP);
CL_Z1 = sum(clexp1 .* Z_2015    .* dz) ./ SigmaZ1   % m^3 mmol^-1 d^-1
CL_Z2 = sum(clexp2 .* Z_2100_vd .* dz) ./ SigmaZ2

%  B) "CI" = gbar / Phy_CZ = <g/P> weighted by ZP (ratio-of-sums form)
CI1 = gbar1 / PhyCZ1  % = (sum gZ dz)/(sum PZ dz)
CI2 = gbar2 / PhyCZ2
CI3 = gbar3 / PhyCZ3



%% Stability & overlap diagnostics for two profile sets + c'(Phy_CZ) comparison
% Needs in workspace: P_2015, Z_2015 (or IZ_2015), P_2100, Z_2100
% Type-II: g(P)=r*P/(K+P); clearance c(P)=r/(K+P); c'(P)=-r/(K+P)^2

% ---- Tunables / constants ----
dz_scalar = 10;                        % m (set to your layer thickness)
Tf = 1.7.^((18 - 30)/10);              % your temperature factor
r  = 3.3 .* Tf;                        % d^-1
K  = 1.2;                               % mmol m^-3
tol = 1e-9;

% Use IZ_2015 if present
if exist('IZ_2015','var'), Z_2015_eff = IZ_2015; else, Z_2015_eff = Z_2015; end

% ---- Compute metrics for each (P,Z) pair ----
M2015 = metrics_for(P_2015, Z_2015_eff, "2015", r, K, dz_scalar, tol);
M2100 = metrics_for(P_2100, Z_2100,      "2100", r, K, dz_scalar, tol);

% ---- Pretty print ----
fprintf('\n=== Stability & overlap metrics (Type-II) ===\n');
fprintf('Units: CI, CL_Z [m^3 mmol^-1 d^-1]; S_* [m^3 mmol^-1 d^-2]; Elasticities [-]; Phy_CZ [mmol m^-3]; DilZ [m]\n\n');

print_metrics(M2015);
print_metrics(M2100);

% ---- Relative changes 2100 vs 2015 ----
dCI    = rel_change(M2100.CI,     M2015.CI);
dPhyCZ = rel_change(M2100.Phy_CZ, M2015.Phy_CZ);
dDilZ  = rel_change(M2100.DilZ,   M2015.DilZ);

fprintf('--- 2100 vs 2015 relative changes ---\n');
fprintf('  CI (encounter-weighted clearance):   %+6.2f %%\n', 100*dCI);
fprintf('  Phy_CZ (Z-weighted prey conc.):      %+6.2f %%\n', 100*dPhyCZ);
fprintf('  Dilution depth (DilZ = 1/NER):       %+6.2f %%\n\n', 100*dDilZ);

% ---- Interpretation hints (elasticity bands) ----
labels = ["strongly saturated","moderately saturated","near-linear"];
edges  = [-Inf, -0.66, -0.33, Inf];
band = @(x) labels(discretize(x, edges));

fprintf('Interpretation:\n');
fprintf('  [%s] E_enc = %.2f -> %s regime at encounter depths.\n', char(M2015.label), M2015.E_enc, char(band(M2015.E_enc)));
fprintf('  [%s] E_enc = %.2f -> %s regime at encounter depths.\n\n', char(M2100.label), M2100.E_enc, char(band(M2100.E_enc)));

%% ----------------- Local functions -----------------
function M = metrics_for(P,Z,label,r,K,dz_scalar,tol)
  % Mask & dz handling
  mask = isfinite(P) & isfinite(Z);
  P = P(mask); Z = Z(mask);
  if isscalar(dz_scalar)
    dz = dz_scalar * ones(size(P));
  else
    dz = dz_scalar(mask);
  end
  epsP = 1e-12;
  P = max(P, epsP);

  % Weights & totals
  SigmaP = sum(P .* dz);       % mmol m^-2
  SigmaZ = sum(Z .* dz);       % mmol m^-2
  wZP    = Z .* P .* dz;  WZP  = sum(wZP);
  wZ     = Z       .* dz;  WZ  = sum(wZ);

  % Local clearance & ingestion
  c   = r ./ (K + P);                    % m^3 mmol^-1 d^-1
  % g = r .* P ./ (K + P);               % (not needed below)

  % Core means
  CI     = r * (sum( wZP ./ (K + P) ) / WZP);   % ZP-weighted <c>
  CL_Z   = sum( c .* wZ ) / WZ;                 % Z-weighted  <c>
  Phy_CZ = sum( P .* Z .* dz ) / SigmaZ;        % Z-weighted prey conc.
  NER    = Phy_CZ / SigmaP;                     % m^-1
  DilZ   = 1 / NER;                             % m

  % Analytic slopes (encounter- and predator-weighted)
  S_enc  = -r * (sum( wZP ./ (K + P).^2 ) / WZP);
  S_Z    = -r * (sum( wZ  ./ (K + P).^2 ) / WZ );
  % Elasticities
  E_enc  = - (sum( wZP .* (P./(K+P)) ) / WZP);
  E_Z    = - (sum( wZ  .* (P./(K+P)) ) / WZ );

  % ---- NEW: slope/elasticity evaluated at Phy_CZ ----
  cprime_at_PhyCZ = -r / (K + Phy_CZ)^2;           % m^3 mmol^-1 d^-2
  E_at_PhyCZ      = -Phy_CZ / (K + Phy_CZ);        % dimensionless
  % Jensen check: <c'(P)>_{ZP} >= c'(Phy_CZ) (both negative; ">= " means less negative)
  jensen_ok = (S_enc - cprime_at_PhyCZ) >= -1e-12;

  % Finite differences
  % (A) Frozen-weights central difference for S_enc
  epsA   = 1e-4; % mmol m^-3 absolute step
  CI0_fw = r * sum( wZP ./ (K + P) ) / WZP;
  CIp_fw = r * sum( wZP ./ (K + (P + epsA)) ) / WZP;
  CIm_fw = r * sum( wZP ./ (K + (P - epsA)) ) / WZP;
  S_fd_frozen = (CIp_fw - CIm_fw) / (2*epsA);

  % (B) Full recompute multiplicative (% step) for elasticity
  epsM   = 1e-4; % 0.01% step
  CI0_full = CI;
  Pm       = P * (1 + epsM);
  wZPp     = Z .* Pm .* dz;  WZPp = sum(wZPp);
  CIp_full = r * sum( wZPp ./ (K + Pm) ) / WZPp;
  E_fd_full = (log(CIp_full) - log(CI0_full)) / log(1 + epsM);

  % Optional closure to catch mistakes
  CI_check = (sum((r.*P./(K+P)) .* Z .* dz) / sum(P.*Z.*dz));
  closure_err_CI = abs(CI - CI_check);
  if closure_err_CI > max(tol, 10*eps)
    warning('[%s] CI closure off by %.3g (check dz / NaNs).', label, closure_err_CI);
  end

  % Pack results
  M = struct('label',label,'SigmaP',SigmaP,'SigmaZ',SigmaZ,...
             'CI',CI,'CL_Z',CL_Z,'Phy_CZ',Phy_CZ,'NER',NER,'DilZ',DilZ,...
             'S_enc',S_enc,'S_Z',S_Z,'E_enc',E_enc,'E_Z',E_Z,...
             'S_fd_frozen',S_fd_frozen,'E_fd_full',E_fd_full,...
             'cprime_at_PhyCZ',cprime_at_PhyCZ,'E_at_PhyCZ',E_at_PhyCZ,...
             'jensen_ok',jensen_ok);
end

function print_metrics(M)
  fprintf('[%s]\n', char(M.label));
  fprintf('  Phy_CZ      = %.6g  | Dilution Depth = %.6g m  | NER = %.6g m^-1\n', M.Phy_CZ, M.DilZ, M.NER);
  fprintf('  CI (ZP-avg) = %.6g  | CL_Z (Z-avg)   = %.6g\n', M.CI, M.CL_Z);
  fprintf('  S_enc (ana) = %.6g  | S_fd_frozen    = %.6g  | Δ=%.3g\n', M.S_enc, M.S_fd_frozen, M.S_fd_frozen - M.S_enc);
  fprintf('  E_enc (ana) = %.6g  | E_fd_full      = %.6g  | Δ=%.3g\n', M.E_enc, M.E_fd_full,   M.E_fd_full   - M.E_enc);

  % ---- NEW: c'(Phy_CZ) comparison & Jensen check ----
  fprintf('  c''(Phy_CZ)  = %.6g  | E_at_PhyCZ     = %.6g\n', M.cprime_at_PhyCZ, M.E_at_PhyCZ);
  sign_note = "OK";
  if ~M.jensen_ok, sign_note = "(!)"; end
  % Remember: numbers are negative; S_enc >= c'(Phy_CZ) means S_enc is less negative (closer to zero)
  fprintf('  Jensen check <c''>_{ZP} >= c''(Phy_CZ): %s  (S_enc - c''(Phy_CZ) = %.3g)\n', char(sign_note), M.S_enc - M.cprime_at_PhyCZ);

  fprintf('  S_Z (pred)  = %.6g  | E_Z (pred)     = %.6g\n\n', M.S_Z, M.E_Z);
end

function r = rel_change(x1,x0)
  r = (x1 - x0) / x0;
end
