%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. S11 Coexistence as driver of CC in MARBL %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))
for startup =1  
% Load regridded lat lon
cd /Users/twrohr/Desktop/data.nosync/SSP585
         load('Longhurst.mat');
          Lon =  LonGrid; Lat = LatGrid;

% Constants
Ra = 6.371e6;                  % Earth radius [m]
deg2rad = pi / 180;
dphi = deg2rad;               % 1 degree in radians
dlambda = deg2rad;

% Latitude vector (center of cells)
lat = linspace(-89.5, 89.5, 180);  % 180 latitudes from -89.5 to +89.5
lat_rad = lat * deg2rad;

% Area per 1x1 degree cell as a function of latitude
area_per_lat = Ra^2 * dphi * dlambda * cos(lat_rad);  % [m^2]

% Repeat across 360 longitudes to make full 180x360 grid
area_map = repmat(area_per_lat', 1, 360);  % [180x360]

% Convert lon for warppign --- Shift longitudes to [0, 360]
        lonvec = linspace(-179.5, 179.5, 360);
        latvec = linspace(-89.5, 89.5, 180);
        lonvec_wrapped = lonvec;
        lonvec_wrapped(lonvec_wrapped < 0) = lonvec_wrapped(lonvec_wrapped < 0) + 360;
        [lonvec_sorted, lon_order] = sort(lonvec_wrapped);

% Colormaps
    % Get full bam colormap (e.g., 256 colors)
    n_colors = 256;

    % --- Define BU Map  ---
    blue_rgb   = [0.35, 0.55, 0.85];  % muted slate blue -- 0.45, 0.60, 0.80
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green -- 0.40, 0.60, 0.40
        % --- Interpolation ---
        n_half = floor(n_colors/2);
        blue_to_white = [linspace(blue_rgb(1), white_rgb(1), n_half)', ...
                         linspace(blue_rgb(2), white_rgb(2), n_half)', ...
                         linspace(blue_rgb(3), white_rgb(3), n_half)'];
        white_to_green = [linspace(white_rgb(1), green_rgb(1), n_colors - n_half)', ...
                          linspace(white_rgb(2), green_rgb(2), n_colors - n_half)', ...
                          linspace(white_rgb(3), green_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapBU = [blue_to_white; white_to_green];
    
    % --- Define TD map  ---
    red_rgb     = [0.78, 0.25, 0.25];  % soft classy red -- 0.75, 0.30, 0.30
    mustard_rgb = [0.90, 0.7, 0.25];  % strong end member yellow -- 0.80, 0.65, 0.20
        % --- Interpolation ---
        red_to_white = [linspace(red_rgb(1), white_rgb(1), n_half)', ...
                        linspace(red_rgb(2), white_rgb(2), n_half)', ...
                        linspace(red_rgb(3), white_rgb(3), n_half)'];
        white_to_mustard = [linspace(white_rgb(1), mustard_rgb(1), n_colors - n_half)', ...
                            linspace(white_rgb(2), mustard_rgb(2), n_colors - n_half)', ...
                            linspace(white_rgb(3), mustard_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapTD = [red_to_white; white_to_mustard];


   
    %%% Discete
    % Interpolate 6 steps from strong to weak BU
    green_start = [0.40, 0.60, 0.40];
    green_end   = [0.85, 0.95, 0.85];
    greens = [linspace(green_start(1), green_end(1), 6)', ...
              linspace(green_start(2), green_end(2), 6)', ...
              linspace(green_start(3), green_end(3), 6)'];
    % Interpolate 6 steps from weak to strong TD
    red_start = [0.90, 0.70, 0.70];
    red_end   = [0.70, 0.00, 0.00];
    reds = [linspace(red_start(1), red_end(1), 6)', ...
            linspace(red_start(2), red_end(2), 6)', ...
            linspace(red_start(3), red_end(3), 6)'];
    % Combine
    custom_cmap = [greens; reds];

    % Model Colors
    model_names = {'iHAMMOCC','CMOC','WOMBAT','OECOv2','MARBL','CanOE','MEDUSA2.1',{'PISCESv2';'IPSL'},{'PISCESv2';'CNRM'},'BFM5.2','COBALTv2'};
    mod_order=[10,5,1,9,2,6,11,8,4,3,7,12];
    cmap1 = brewermap(13,'Oranges');
    cmap2 = brewermap(15,'Greens');
    cmap3 = brewermap(11,'Purples');
    cmap_mod = [cmap2([5,7,11,13,15],:); cmap1([5,7,11,13],:); cmap3([7,11],:); 0 0 0];

    CoreVar_Names = {'Phytoplankton Biomass',     'Zooplankton Biomass',         'NPP',                        'Grazing Losses',            'EP',                 ...
                     'Spec. Division Rate',       'Spec. Grazing Rate',          'GSP',                        'Z Loss rate',               'Non-Grazing Losses', ...
                     'Spec. Grazing Loss Rate',   'Fraction Small P',            'Fraction Small Z',           'SST' ,                      'MLD' ,               ...
                     'Spec. Loss Rate',           'Spec. Non-grazing Loss Rate'  'Export Efficeincy'           'Ecosystem Efficeincy'       'Surface [Phyto]' ,  ...
                     'DBM [Phyto]',               'DBM Depth',                   'Dilution Index',             'Diaz Frac',                 'Diatom Frac',           'Non- Diatom Frac'...
                     'Temperature (g)',           'Vertical Dist (g) ',          'Community Composition (g)',  'Phyto. Biomass (g)',        'Zooplankton Biomass (g)',...
                     'Temperature (n)',           'Vertical Dist (n) ',          'Community Composition (n)',  'Phyto. Biomass (n)',        'Zooplankton Biomass (n)',...
                     'Temperature (T)',           'Vertical Dist (T) ',          'Community Composition (T)',  'Phyto. Biomass (T)',        'Zooplankton Biomass (T)',};

     decomp_names = {{'$\Delta$ Loss Rate via' ; '$\Delta$ Temperture'}, {'$\Delta$ Loss Rate via' ; '$\Delta$ Vert. Dist.'},...
                     {'$\Delta$ Loss Rate via' ; '$\Delta$ Comm. Comp.'}, {'$\Delta$ Loss Rate via' ; '$\Delta$ Phyto Biomass'},...
                     {'$\Delta$ Loss Rate via' ; '$\Delta$ Zoo Biomass'}}
     % decomp_names =  { '$\frac{\partial loss rate}{partial Temeprature}$' , '$\frac{\partial loss rate}{partial Vert. Dist}$' ,...
     %                '$\frac{\partial loss rate}{partial Comm. Comp.}$' , '$\frac{\partial loss rate}{partial Phyto Biomass}$' ,...
     %                 '$\frac{\partial loss rate}{partial Zoo Biomass}$' };


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'}; 
end


% =====================================================
%% Plot
% =====================================================
clc; close all
%%%% 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 = 9; % DBM Formaiton in CNRM
model = 5; 

% Load Model
for Load_mod =[1]

  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_Inv_rw'); 
          fprintf(['\n Dowloaded ' mods(mod_order(model)).name])
          last_load=model;
  end
end

  %%%% Pick 6 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 

% Variables
vs = [101, 29, 12, 29+5] ;  % DBM Formation in CNRM 
clims_low = [-50, -.05, -70, -.005];

% Create Figure
figure('Units','normalized','Position',[0 0 .6 1], 'Color','w');
tl = tiledlayout(3,2,'TileSpacing','compact','Padding','compact');
plotvar       =  squeeze(CoreVar_Inv_adj(1,1,:,:));
Ra = georefcells([min(Lat_t(:,1)) max(Lat_t(:,1))], ...
                 [min(Lon_t(1,:)) max(Lon_t(1,:))], size(plotvar));
for i = 1:4
    
    % Assign Variables
    for Assign = 1
        if vs(i) == 120
         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);
         vunits = '(m)'; vname = {'Dilution Index'};
         climvd  = prctile(plotvar(:),[2,98]);
         climv   = [-max(abs(climvd)), max(abs(climvd))];
        elseif vs(i) == 105
         ClearR_2015  = mean( squeeze(CoreVar_Inv_adj(28,1:60,:,:)) ./ squeeze(CoreVar_Inv_adj(2,1:60,:,:))  ); 
         ClearR_2100  = mean( squeeze(CoreVar_Inv_adj(28,61:end,:,:)) ./ squeeze(CoreVar_Inv_adj(2,61:end,:,:))  ); 
         plotvar       = squeeze( ClearR_2100 - ClearR_2015);
         vunits = 'm$^2$ mmol$^{-1}$ d$^{-1}$'; vname = {'Clearance Rate'};
         climvd  = prctile(plotvar(:),[2,98]);
         climv   = [-max(abs(climvd)), max(abs(climvd))];
         elseif vs(i) == 101
         % Change in co-existence
         plotvar = squeeze(mean(CoreVar_Inv_adj(12,1:60,:,:),2));
         plotvar1 = abs( 2*plotvar-100);
         plotvar = squeeze(mean(CoreVar_Inv_adj(12,61:end,:,:),2));
         plotvar2 = abs( 2*plotvar-100);
         plotvar = -(plotvar2-plotvar1);
         climvd  = prctile(plotvar(:),[2,98]);
         climv   = [-max(abs(climvd)), max(abs(climvd))];
         vunits = '$ \% $'; vname = {'Coexistence'};
         else
         plotvar = squeeze(CoreVar_Inv_rw(vs(i),:,:));
         climvd  = prctile(plotvar(:),[2,98]);
         climv   = [-max(abs(climvd)), max(abs(climvd))];
         vunits  = CoreVar_Inv_Units{vs(i)}; 
         vname = CoreVar_Names{vs(i)};
        end
    end

    % Plot
    ax = nexttile(tl); axes(ax); axis off
    axesm('robinson','MapLonLimit',[0 360],'MapLatLimit',[-80 80],'Origin',[0 180 0]);
    tightmap

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

    % Styling
    colormap(flipud(cmocean('balance')));
    caxis([clims_low(i),-clims_low(i)]);
    title(vname, 'Interpreter','latex', 'FontSize', 20);

    % Colorbar 
    cb = colorbar(ax,'southoutside');
    cb.TickLabelInterpreter = 'latex';
    cb.Label.Interpreter   = 'latex';
    cb.Label.String        = vunits;
    cb.FontSize = 14;
end



%% Response Curve plot
for respons=[1]
phy    = [0:.01:2];
Tf       = 1.7.^((18 - 30)/10);  % @ 18 degrees
g        = (3.3.*Tf.*phy) ./(phy+1.2);
agg      =  .01.*phy.^(1.75);

 figure('Units','normalized','Position',[0 0 .3 .3], 'Color','w');
 
 subplot(2,1,1)
 plot(phy,g,'k','linewidth',3); hold on 
 plot(phy,agg,'r','linewidth',3)
 yline(.02,'k:'); yline(.01,'k:')
 ylabel('Loss Rate (mmol m$^{-3}$)','Interpreter','latex','FontSize',12)
 xlabel('[Phyto] (mmol m$^{-3}$)','Interpreter','latex','FontSize',12)
 subplot(2,1,2)
 %plot(phy,agg,'k','linewidth',3); hold on 
 plot(phy,agg,'r','linewidth',3)
 yline(.02,'k:','linewidth',2); yline(.01,'k:','linewidth',2)
 ylabel('Loss Rate (mmol m$^{-3}$)','Interpreter','latex','FontSize',12)
 xlabel('[Phyto] (mmol m$^{-3}$)','Interpreter','latex','FontSize',12)
end