cd /project2/tas1/tanzh/Codes/GCM_analysis/from_Orli/Output_2018_01a/
addpath('/project2/tas1/tanzh/Codes/Output/');

%% Load datasets
load('aqua_data_20200316.mat');
load('cmip_data_20200316_JJA.mat')
load('cmip_data_20200316_JJA_raw.mat')

%% Figure 1: Line plots
do_eqn = 1;
if do_eqn
    ncdata_synth = aqua_eqn;
    with_p6k = 1;
else
    ncdata_synth = aqua_ann;
    with_p6k = 0;
end

lat = ncdata_synth.ctrl.lat;
latb = [-90; (lat(1:end-1)+lat(2:end))/2; 90];
areafac = diff(sind(latb))./sum(diff(sind(latb)));  

fieldnm_all = {'t_surf', 'mse_atm', 'lhf', 'acre_lw'};
fieldnm_out = {'\DeltaSST (K)', '\DeltaMSE (kJ/kg)', ...
               '\DeltaLHF (W/m^2)', '\Delta LW ACRE (W/m^2)'};
           
figure;
for ifield =1:4
    fieldnm = fieldnm_all{ifield};

    val_co2  = ncdata_synth.co2.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    val_cld  = ncdata_synth.cld.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    val_qsfc = ncdata_synth.qsfc.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    val_wnd  = ncdata_synth.wnd.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    val_tot  = ncdata_synth.warm.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    if with_p6k
        val_p6k  = ncdata_synth.p6k.(fieldnm) - ncdata_synth.ctrl.(fieldnm);
    end
    
    plt(ifield) = subplot(2,2,ifield); hold all;
    lne(1) =  plot(lat, val_tot,  'k-',  'LineWidth', 1.5);
    lne(2) =  plot(lat, val_co2,  'r-',  'LineWidth', 1.5);
    lne(3) =  plot(lat, val_cld,  'b-',  'LineWidth', 1.5);
    lne(4) =  plot(lat, val_wnd, '-', 'Color', [0 .7 0],  'LineWidth', 1.5);
    lne(5) =  plot(lat, val_qsfc, 'm-',  'LineWidth', 1.5);
    if with_p6k
        lne(6) =  plot(lat, val_p6k,  'm--',  'LineWidth', 1.5);
    end
    plot(lat, 0*lat, '-', 'Color', [.5 .5 .5],  'LineWidth', 1.5)
    
    xlim([0 90]);
    set(gca, 'TickDir', 'both', 'XTick', [0 20 40 60 90]);
    title(fieldnm_out{ifield});
    if ifield == 1
        if with_p6k
            ylim([-10 20]);
        else
            ylim([-10 10]);
        end
        axis_tmp1 = gca;
        axis_tmp2 = axes('position',get(axis_tmp1,'position'),'visible','off');
        if with_p6k
        	lgd_hdl1 = legend(axis_tmp1, lne(1:3), ...
                {['Total: ', num2str(areafac'*val_tot, '%+4.1f')], ...
                 ['CO2: ', num2str(areafac'*val_co2, '%+4.1f')], ...
                 ['Cloud: ', num2str(areafac'*val_cld, '%+4.1f')]}, ...
                 'Position', [0.09 0.81 0.13 0.11]);
            lgd_hdl2 = legend(axis_tmp2, lne(4:6), ...
                {['SFC wind: ', num2str(areafac'*val_wnd, '%+4.1f')], ...
                 ['SFC q: ', num2str(areafac'*val_qsfc, '%+4.1f')], ...
                 ['SFC q +6K: ', num2str(areafac'*val_p6k, '%+4.1f')]}, ...
                 'Position', [0.255 0.81 0.13 0.11]);
            lgd_hdl1.ItemTokenSize = [18 18]; lgd_hdl1.Box = 'off'; lgd_hdl1.FontSize = 10;
            lgd_hdl2.ItemTokenSize = [18 18]; lgd_hdl2.Box = 'off'; lgd_hdl2.FontSize = 10;
        else
        	lgd_hdl1 = legend(axis_tmp1, lne(1:3), ...
                {['Total: ', num2str(areafac'*val_tot, '%+4.1f')], ...
                 ['CO2: ', num2str(areafac'*val_co2, '%+4.1f')], ...
                 ['Cloud: ', num2str(areafac'*val_cld, '%+4.1f')]}, ...
                 'Position', [0.16 0.56 0.13 0.11]);
            lgd_hdl2 = legend(axis_tmp2, lne(4:5), ...
                {['SFC wind: ', num2str(areafac'*val_wnd, '%+4.1f')], ...
                 ['SFC q: ', num2str(areafac'*val_qsfc, '%+4.1f')]}, ...
                 'Position', [0.325 0.56+0.11/3 0.13 0.11/3*2]);
            lgd_hdl1.ItemTokenSize = [18 18]; lgd_hdl1.Box = 'off'; lgd_hdl1.FontSize = 10;
            lgd_hdl2.ItemTokenSize = [18 18]; lgd_hdl2.Box = 'off'; lgd_hdl2.FontSize = 10;

        end
        set(axis_tmp1, 'YTick', -10:5:30, 'XTickLabel', {}); 

    elseif ifield == 2
        ylim([-5 30]);
        set(gca, 'YTick', -10:5:30, 'XTickLabel', {});
    elseif ifield == 3
        if with_p6k
            ylim([-40 80]);
        else
            ylim([-30 60]);
        end
        xlb_hdl = xlabel('Latitude (^o)');
    elseif ifield == 4
        if with_p6k
            ylim([-25 20])
        else
            ylim([-15 10])
        end
    end
end
set(gcf, 'Position', [200 100 800 480])
set(plt(1), 'Position', [.08 .55 .4 .38])
set(plt(2), 'Position', [.55 .55 .4 .38])
set(plt(3), 'Position', [.08 .1 .4 .38])
set(plt(4), 'Position', [.55 .1 .4 .38])
if with_p6k
    set(xlb_hdl, 'Position', [100 -55 -1])
else
    set(xlb_hdl, 'Position', [100 -40 -1])
end
ann_hdl(1) = annotation('textbox', 'String', '(a)', 'linestyle', 'none', ...
                        'position', [0.08 0.94 0.05 0.05], 'FontSize', 14);
ann_hdl(2) = annotation('textbox', 'String', '(b)', 'linestyle', 'none', ...
                        'position', [0.55 0.94 0.05 0.05], 'FontSize', 14);
ann_hdl(3) = annotation('textbox', 'String', '(c)', 'linestyle', 'none', ...
                        'position', [0.08 0.49 0.05 0.05], 'FontSize', 14);
ann_hdl(4) = annotation('textbox', 'String', '(d)', 'linestyle', 'none', ...
                        'position', [0.55 0.49 0.05 0.05], 'FontSize', 14);


%% Figure 2: Contour plot of responses

do_eqn = 1;
if do_eqn
    ncdata_synth = aqua_eqn;
    with_p6k = 1; ncol = 6;
    case_all_synth = {'warm', 'co2', 'cld', 'wnd', 'qsfc', 'p6k'};
    title_all = {'Total', 'CO2', 'Cloud', 'SFC wind', 'SFC q', 'SFC q +6K'};
else
    ncdata_synth = aqua_ann;
    with_p6k = 0; ncol = 5;
    case_all_synth = {'warm', 'co2', 'cld', 'wnd', 'qsfc'};
    title_all = {'Total', 'CO2', 'Cloud', 'SFC wind', 'SFC q'};
end

lat = ncdata_synth.ctrl.lat;
pfull = ncdata_synth.ctrl.pfull;  

figure; 

% Plot Temperature
fieldnm = 'temp'; 
for i = 1:ncol
    plt(i) = subplot(3,6,i); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull(5:end)/1000, val_resp(:,5:end)', -20:20, 'w-','LineWidth',0.5)
    contour(lat, pfull/1000, val_ref', 180:20:320, 'k-')
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    set(gca, 'XTickLabel', {});
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel('T (K)');
    end
    colormap(cmap_tmp);
    title(title_all{i});
end
for i = 1:ncol
    set(plt(i), 'Position', [-0.08+i*0.15, 0.68, 0.13, 0.27]);
end

% Plot Zonal Wind
fieldnm = 'ucomp'; 
for i = 1:ncol
    plt(i+6) = subplot(3,6,i+6); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull(5:end)/1000, val_resp(:,5:end)', -20:20, 'w-','LineWidth',0.5)
    contour(lat, pfull/1000, val_ref', -50:10:-10, 'k--')
    contour(lat, pfull/1000, val_ref', 10:10:100, 'k-')
    contour(lat, pfull/1000, val_ref', [0 0], 'k-', 'LineWidth', 1.5)
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    set(gca, 'XTickLabel', {});
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel('U (m/s)');
    end
    colormap(cmap_tmp);
end
for i = 1:ncol
    set(plt(i+6), 'Position', [-0.08+i*0.15, 0.38, 0.13, 0.27]);
end


% Plot Streamfunction
fieldnm = 'psi'; 
for i = 1:ncol
    plt(i+12) = subplot(3,6,i+12); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull/1000, val_ref', -150:5:-5, 'k--')
    contour(lat, pfull/1000, val_ref', 5:5:150, 'k-')
    contour(lat, pfull/1000, val_ref', [0 0], 'k-', 'LineWidth', 1.5)
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel([char(936), ' (10^{10} kg/s)']);
    end
    colormap(cmap_tmp);
    if i == 3
        xlb_hdl = xlabel('Latitude (^o)');
    end
end
for i = 1:ncol
    set(plt(i+12), 'Position', [-0.08+i*0.15, 0.08, 0.13, 0.27]);
end
cbar_hdl = colorbar();
set(cbar_hdl, 'Position', [0.06+ncol*0.15 0.08 0.015 0.87])
set(xlb_hdl, 'Position', [105.0    1.15   -1.0000])
set(gcf, 'Position', [100 100 800 480])


%% Figure 3: Quantitative shifts

figure; hold all; 

do_eqn = 1;
if do_eqn
    ncdata_synth = aqua_eqn;
    with_p6k = 1; 
    case_all_synth = {'warm', 'co2', 'cld', 'wnd', 'qsfc', 'p6k'};
    color_val = [0 0 0; 1 0 0; 0 0 1; 0 .5 0; 1 0 1; 1 0.5 1];
    sym_val = {'x', 'x', 'x', 'x', 'x', '+'};
else
    ncdata_synth = aqua_ann;
    with_p6k = 0; 
    case_all_synth = {'warm', 'co2', 'cld', 'wnd', 'qsfc'};
    color_val = [0 0 0; 1 0 0; 0 0 1; 0 .5 0; 1 0 1];
    sym_val = {'x', 'x', 'x', 'x', 'x'};
end

rectangle('Position',[0.05,-0.99,4.95,1.99],'FaceColor', [.8 .8 .8], 'EdgeColor', 'none');
for icase = 1:length(case_all_synth)
    casenm = case_all_synth{icase};
    plt(icase) = plot(1, ...
         ncdata_synth.(casenm).us0 - ncdata_synth.ctrl.us0, sym_val{icase}, ...
         'Color', color_val(icase,:), 'LineWidth', 2);
    plt(icase+10) = plot(2, ...
         ncdata_synth.(casenm).pme0 - ncdata_synth.ctrl.pme0, sym_val{icase}, ...
         'Color', color_val(icase,:), 'LineWidth', 2);
    plt(icase+20) = plot(3, ...
         ncdata_synth.(casenm).psi0 - ncdata_synth.ctrl.psi0, sym_val{icase}, ...
         'Color', color_val(icase,:), 'LineWidth', 2);
    plt(icase+30) = plot(4, ...
         ncdata_synth.(casenm).umax0 - ncdata_synth.ctrl.umax0, sym_val{icase}, ...
         'Color', color_val(icase,:), 'LineWidth', 2);
end
plot([0 5], [0 0], 'k-');
if with_p6k
    lgd_hdl = legend(plt(1:6), {'Total', 'CO2', 'Cloud', 'SFC wind', 'SFC q', 'SFC q +6K'});
else
    lgd_hdl = legend(plt(1:5), {'Total', 'CO2', 'Cloud', 'SFC wind', 'SFC q'});
end
xlim([0.5 4.5]); ylim([-1 3]);
ylabel('Shift (\circ poleward)')
set(gca, 'Position', [0.12 0.12 0.58 0.8], 'TickDir', 'both', 'XTick', 1:4, 'XTickLabel', ...
    {[char(981),'_{u=0}'], [char(981),'_{P-E=0}'], [char(981),'_{', char(936), '=0}'], [char(981),'_{u=max}']});
set(gcf, 'Position', [100 100 400 240])
% set(lgd_hdl, 'Position', [0.77 0.52 0.2 0.4])
set(lgd_hdl, 'Position', [0.75 0.12 0.2 0.4])


%% Supp. Figure 1: Contour: total, sum, residual, Interactive minus Locked

do_eqn = 1;
if do_eqn
    ncdata_synth = aqua_eqn;
else
    ncdata_synth = aqua_ann;
end

case_all_synth = {'warm', 'sum', 'res', 'iml'};
title_all = {'Total Locked', 'Sum of 4 Components', 'Residual', 'Interactive - Locked'};

lat = ncdata_synth.ctrl.lat;
pfull = ncdata_synth.ctrl.pfull; 

figure; 
% Plot Temperature
fieldnm = 'temp'; 
for i = 1:4
    plt(i) = subplot(3,4,i); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull(5:end)/1000, val_resp(:,5:end)', -20:20, 'w-','LineWidth',0.5)
    contour(lat, pfull/1000, val_ref', 180:20:320, 'k-')
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    set(gca, 'XTickLabel', {});
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel('T (K)');
    end
    colormap(cmap_tmp);
    title(title_all{i});
end
for i = 1:4
    set(plt(i), 'Position', [-0.15+i*0.22, 0.68, 0.2, 0.27]);
end

% Plot Zonal Wind
fieldnm = 'ucomp'; 
for i = 1:4
    plt(i+4) = subplot(3,4,i+4); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull(5:end)/1000, val_resp(:,5:end)', -20:20, 'w-','LineWidth',0.5)
    contour(lat, pfull/1000, val_ref', -50:10:-10, 'k--')
    contour(lat, pfull/1000, val_ref', 10:10:100, 'k-')
    contour(lat, pfull/1000, val_ref', [0 0], 'k-', 'LineWidth', 1.5)
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    set(gca, 'XTickLabel', {});
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel('U (m/s)');
    end
    colormap(cmap_tmp);
end
for i = 1:4
    set(plt(i+4), 'Position', [-0.15+i*0.22, 0.38, 0.2, 0.27]);
end


% Plot Streamfunction
fieldnm = 'psi'; 
for i = 1:4
    plt(i+8) = subplot(3,4,i+8); hold all;
    casenm = case_all_synth{i};
    val_ref  = ncdata_synth.ctrl.(fieldnm);
    val_resp = ncdata_synth.(casenm).(fieldnm) - val_ref;
    contourf(lat, pfull/1000, val_resp', -20:20, 'LineStyle', 'none')
    contour(lat, pfull/1000, val_ref', -150:5:-5, 'k--')
    contour(lat, pfull/1000, val_ref', 5:5:150, 'k-')
    contour(lat, pfull/1000, val_ref', [0 0], 'k-', 'LineWidth', 1.5)
    set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
    cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
    caxis([-8 8]); xlim([0 90]); ylim([0 1]); 
    if i > 1
        set(gca, 'YTickLabel', {});
    else
        ylabel([char(936), ' (10^{10} kg/s)']);
    end
    colormap(cmap_tmp);
    if i == 2
        xlb_hdl = xlabel('Latitude (^o)');
    end
end
for i = 1:4
    set(plt(i+8), 'Position', [-0.15+i*0.22, 0.08, 0.2, 0.27]);
end
cbar_hdl = colorbar();
set(cbar_hdl, 'Position', [0.95 0.08 0.015 0.87])
set(xlb_hdl, 'Position', [105.0    1.15   -1.0000])
set(gcf, 'Position', [100 100 800 480])


%% Supp. Figure 2: Cloud changes in interactive runs

do_eqn = 1;
if do_eqn
    ncdata_synth = aqua_eqn;
else
    ncdata_synth = aqua_ann;
end

lat = ncdata_synth.ctrl.lat;
latb = [-90; (lat(1:end-1)+lat(2:end))/2; 90];
areafac = diff(sind(latb))./sum(diff(sind(latb)));
areatot = 4*pi*(6371e3)^2;
pfull = ncdata_synth.ctrl.pfull; 

figure;
plt(1) = subplot(1,3,1); hold all;
val_ref  = ncdata_synth.ctrl_fr.cld_amt*100;
val_resp = ncdata_synth.warm_fr.cld_amt*100 - val_ref;
contourf(lat, pfull/1000, val_resp', -40:40, 'LineStyle', 'none')
contour(lat, pfull/1000, val_ref', 10:10:300, 'k-')
set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
caxis([-16 16]); xlim([0 90]); ylim([0 1]); 
ylabel('Sigma level');
colormap(cmap_tmp);
title('Cloud Fraction (%)');

plt(2) = subplot(1,3,2); hold all;
val_ref  = ncdata_synth.ctrl_fr.liq_wat*1e6;
val_resp = ncdata_synth.warm_fr.liq_wat*1e6 - val_ref;
contourf(lat, pfull/1000, val_resp', -400:400, 'LineStyle', 'none')
contour(lat, pfull/1000, val_ref', 10:10:300, 'k-')
set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
set(gca, 'YTickLabel', {}); xlb_hdl = xlabel('Latitude (^o)');
cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
caxis([-16 16]); xlim([0 90]); ylim([0 1]); 
title('Cloud Liquid (mg/kg)');
colormap(cmap_tmp);

plt(3) = subplot(1,3,3); hold all;
val_ref  = ncdata_synth.ctrl_fr.ice_wat*1e6;
val_resp = ncdata_synth.warm_fr.ice_wat*1e6 - val_ref;
contourf(lat, pfull/1000, val_resp', -400:400, 'LineStyle', 'none')
contour(lat, pfull/1000, val_ref', 10:10:300, 'k-')
set(gca, 'YDir', 'reverse', 'TickDir', 'both', 'XTick', [0 20 60 90], 'YTick', 0.2:0.2:1);
set(gca, 'YTickLabel', {});
cmap_tmp = redblue(20); cmap_tmp = cmap_tmp([2 4:17 19],:);
caxis([-16 16]); xlim([0 90]); ylim([0 1]); 
title('Cloud Ice (mg/kg)');
colormap(cmap_tmp);
for i = 1:3
    set(plt(i), 'Position', [-0.2+i*0.285, 0.155, 0.25, 0.75]);
end
set(xlb_hdl, 'Position', [45 1.1 -1])
cbar_hdl = colorbar();
set(cbar_hdl, 'Position', [0.92 0.155 0.03 0.75])
set(gcf, 'Position', [100 100 600 240])

 
%% Figure 4: CMIP5 scatter plots

band_width_all = diff(sind(cmip_analysis.band_boundary));
model_nm_all = cmip_analysis.model_nm_all;

model_idx = 1:length(model_nm_all); 
band_width = zeros(size(band_width_all));
band_width(2:3) = -band_width_all(2:3)/sum(band_width_all(2:3));
band_width(7:8) =  band_width_all(7:8)/sum(band_width_all(7:8));
data_lwacre = band_width*squeeze(cmip_analysis.LWACREa(:,model_idx,2) - ...
    cmip_analysis.LWACREa(:,model_idx,1));  % 1: Historical; 2: RCP 8.5
data_lhf    = band_width*squeeze(cmip_analysis.LHFa(:,model_idx,2) - ...
    cmip_analysis.LHFa(:,model_idx,1));
data_mse925 = band_width*squeeze(cmip_analysis.T925a(:,model_idx,2) - ...
    cmip_analysis.T925a(:,model_idx,1)) * 1004.7/1000 + ...
    band_width*squeeze(cmip_analysis.q925a(:,model_idx,2) - ...
    cmip_analysis.q925a(:,model_idx,1)) * 2.5e6/1000;
data_hce_clim = squeeze(cmip_analysis.hce_lat(model_idx,1));
data_hce_resp = squeeze(cmip_analysis.hce_lat(model_idx,2) - ...
    cmip_analysis.hce_lat(model_idx,1));
data_edj_clim = squeeze(cmip_analysis.edj_lat(model_idx,1));
data_edj_resp = squeeze(cmip_analysis.edj_lat(model_idx,2) - ...
    cmip_analysis.edj_lat(model_idx,1));

plt_edj = 1; % 1 for Fig 4 (EDJ); 0 for Fig S9 (HCE)
if plt_edj
    title_all = {'JJA SH EDJ Climatology (^o)', ...
                 '\DeltaLHF contrast (W/m^2)', '\Deltam_{925} contrast (kJ/kg)', ...
                 '\DeltaACRE_{LW} contrast (W/m^2)', '\Delta(ACRE_{LW}+1.09LHF) contrast', ...
                 '\Delta(ACRE_{LW}+3.72m_{925}) contrast'};
else
    title_all = {'JJA SH HC Edge Climatology (^o)', ...
                 '\DeltaLHF contrast (W/m^2)', '\Deltam_{925} contrast (kJ/kg)', ...
                 '\DeltaACRE_{LW} contrast (W/m^2)', '\Delta(ACRE_{LW}+1.06LHF) contrast', ...
                 '\Delta(ACRE_{LW}+3.20m_{925}) contrast'};
end

figure; 
for i = 1:6
    plt(i) = subplot(2,3,i); hold all;
    switch i
        case 1
            if plt_edj
                plot_data_x = data_edj_clim; xrange = [-55 -35];
            else
                plot_data_x = data_hce_clim; xrange = [-30 -22];
            end
        case 2
            plot_data_x = data_lhf';      xrange = [-2.5 17.5];
        case 3
            plot_data_x = data_mse925';   xrange = [-2.5 7.5];
        case 4
            plot_data_x = data_lwacre';   xrange = [-4 6];
        case 5
            if plt_edj
                plot_data_x = data_lwacre' + 1.0851*data_lhf';    xrange = [-2.5 22.5];
            else
                plot_data_x = data_lwacre' + 1.0649*data_lhf';    xrange = [-2.5 22.5];
            end
        case 6
            if plt_edj
                plot_data_x = data_lwacre' + 3.7167*data_mse925'; xrange = [-5 25];
            else
                plot_data_x = data_lwacre' + 3.1989*data_mse925'; xrange = [-5 25];
            end
    end
    
    if plt_edj
        plot_data_y = data_edj_resp; 
    else
        plot_data_y = data_hce_resp; 
    end
    
    % Plot the regression line
    regre_param = [ones(size(plot_data_x)) plot_data_x]\plot_data_y;
    hold all; plot(xrange, regre_param(1) + regre_param(2)*xrange, '-', 'Color', [.5 .5 .5])
    
    plot(plot_data_x, plot_data_y, 'kx')
    plot(xrange, [0 0], 'k:')
    [rval, pval] = corrcoef(plot_data_x, plot_data_y); rval=rval(1,2); pval=pval(1,2);
    title(title_all{i});
    if plt_edj
        ylim([-6 2]); 
    else
        
    end
    xlim(xrange);
    set(gca, 'TickDir', 'both');
    if i == 1
        if plt_edj
            ylb_hdl = ylabel('JJA Southern Hemishphere EDJ shift (^o)');
            set(ylb_hdl, 'Position', [-58 -8 -1])
        else
            ylb_hdl = ylabel('JJA Southern Hemishphere HC Edge shift (^o)');
            set(ylb_hdl, 'Position', [-31.5 -3 -1])
        end
    end
    if mod(i,3)~=1
        set(gca, 'YTickLabel', {});
    end
    ann_hdl(i) = annotation('textbox', 'String', ['r = ', num2str(rval, '%4.2f')], ...
        'linestyle', 'none');
end
for i = 1:3
    set(plt(i),       'Position', [-0.23+0.31*i, 0.56, 0.29, 0.36])
    set(ann_hdl(i),   'Position', [-0.05+0.31*i, 0.79, 0.10, 0.10])
    set(plt(i+3),     'Position', [-0.23+0.31*i, 0.07, 0.29, 0.36])
    set(ann_hdl(i+3), 'Position', [-0.05+0.31*i, 0.30, 0.10, 0.10])
end
ann_hdl(1) = annotation('textbox', 'String', '(a)', 'linestyle', 'none', ...
                        'position', [0.31*1-0.22 0.58 0.05 0.05], 'FontSize', 14);
ann_hdl(2) = annotation('textbox', 'String', '(b)', 'linestyle', 'none', ...
                        'position', [0.31*2-0.22 0.58 0.05 0.05], 'FontSize', 14);
ann_hdl(3) = annotation('textbox', 'String', '(c)', 'linestyle', 'none', ...
                        'position', [0.31*3-0.22 0.58 0.05 0.05], 'FontSize', 14);
ann_hdl(4) = annotation('textbox', 'String', '(d)', 'linestyle', 'none', ...
                        'position', [0.31*1-0.22 0.09 0.05 0.05], 'FontSize', 14);
ann_hdl(5) = annotation('textbox', 'String', '(e)', 'linestyle', 'none', ...
                        'position', [0.31*2-0.22 0.09 0.05 0.05], 'FontSize', 14);
ann_hdl(6) = annotation('textbox', 'String', '(f)', 'linestyle', 'none', ...
                        'position', [0.31*3-0.22 0.09 0.05 0.05], 'FontSize', 14);
set(gcf, 'Position', [100 100 800 480])


%% Figure S8: CMIP5 line plots for LWACRE, LHF, MSE925
title_all = {'JJA \DeltaACRE_{LW} (W/m^2)', 'JJA \DeltaLHF (W/m^2)', 'JJA \Deltam_{925} (kJ/kg)'};

model_nm_all = fieldnames(cmip_raw);
figure; 
for imodel = 1:length(model_nm_all)
    cmip_data = cmip_raw.(model_nm_all{imodel}); 

    subplot(1,3,1); hold all;
    plot(cmip_data.lat2d, diff(cmip_data.acre_lw,1,2), 'b-');
    subplot(1,3,2); hold all;
    plot(cmip_data.lat2d, diff(cmip_data.lhf,1,2), 'b-');
    subplot(1,3,3); hold all;
    plot(cmip_data.lat3d, diff(cmip_data.mse925,1,2), 'b-');
end
%
set(gcf, 'Position', [100 100 800 240])
for i = 1:3
    subplot(1,3,i); hold all;
    set(gca, 'Position', [0.33*i-0.28, 0.18, 0.28, 0.7]);
    if i == 2
        xlabel('Latitude (^o)'); 
    end
    xlim([-90 10]); title(title_all{i});
    set(gca, 'XTick', [-80 -60 -30 -10 10]);
    yrange = get(gca, 'YLim');
    plot([-80 -80], yrange, '--', 'Color', [.3 .3 .3])
    plot([-60 -60], yrange, '--', 'Color', [.3 .3 .3])
    plot([-30 -30], yrange, '--', 'Color', [.3 .3 .3])
    plot([-10 -10], yrange, '--', 'Color', [.3 .3 .3])
    plot([-90 10], [0 0], '--', 'Color', [.3 .3 .3])
end

ann_hdl(1) = annotation('textbox', 'String', '(a)', 'linestyle', 'none', ...
                        'position', [0.33*1-0.05 0.26 0.05 0.05], 'FontSize', 14);
ann_hdl(2) = annotation('textbox', 'String', '(b)', 'linestyle', 'none', ...
                        'position', [0.33*2-0.05 0.26 0.05 0.05], 'FontSize', 14);
ann_hdl(3) = annotation('textbox', 'String', '(c)', 'linestyle', 'none', ...
                        'position', [0.33*3-0.05 0.26 0.05 0.05], 'FontSize', 14);