clear

load('matres_fgco2_decomposition_trend_nino_nina.mat'); load('matres_fgco2_trend_nino_nina.mat'); 
load('matres_spco2tnt_pdf_nino_nina.mat'); load('matres_rainfall_pdf_nino_nina.mat')

inds=find(lon>=140 & lon<=285 & lat>=-15 & lat<=15);
inds_nan=find(lon>=276 & lon<=285 & lat>=9 & lat<=15);
area=grid_area(lon,lat);
trend_unit=100; %per century
idx=[2:8]; %[2 3 4 6 7 8 10 11 12]; %model with decomposition terms

fgco2_trend_model_nino=fgco2_trend_model_nino*trend_unit;
fgco2_trend_ste_model_nino=fgco2_trend_ste_model_nino*trend_unit;%*2;
%fgco2_trend_ste_model_nino(end-1:end)=fgco2_trend_ste_model_nino(end-1:end)/2;

fgco2_trend_model_nina=fgco2_trend_model_nina*trend_unit;
fgco2_trend_ste_model_nina=fgco2_trend_ste_model_nina*trend_unit;%*2;
%fgco2_trend_ste_model_nina(end-1:end)=fgco2_trend_ste_model_nina(end-1:end)/2;

data_varall={'fgco2_dpco2_trend_model_kw_m_nino','fgco2_kw_trend_model_kw_ano_nino','fgco2_kw_trend_model_dpco2_m_nino','fgco2_dpco2_trend_model_dpco2_ano_nino',...
             'fgco2_dpco2_trend_model_kw_m_nina','fgco2_kw_trend_model_kw_ano_nina','fgco2_kw_trend_model_dpco2_m_nina','fgco2_dpco2_trend_model_dpco2_ano_nina',};
for v=1:length(data_varall);
    eval(['data=',char(data_varall(v)),';']);
    [mm,jm,im]=size(data);
    for i=1:mm;
        datax=squeeze(data(i,:,:))*trend_unit;
        datax=datax.*area; datax(inds_nan)=nan; 
        fgco2_trend_comp(v,i)=nansum(datax(inds))/1e15;
    end
end

data=fgco2_trend_comp; clear fgco2_trend_comp; fgco2_trend_comp(:,1)=data(:,end); 
data=fgco2_trend_model_nino;     clear fgco2_trend_model_nino;     fgco2_trend_model_nino(1)=data(end-1); % model mean is 8 model mean including access
data=fgco2_trend_ste_model_nino; clear fgco2_trend_ste_model_nino; fgco2_trend_ste_model_nino(1)=data(end-1);
data=fgco2_trend_model_nina;     clear fgco2_trend_model_nina;     fgco2_trend_model_nina(1)=data(end-1); % model mean is 8 model mean including access
data=fgco2_trend_ste_model_nina; clear fgco2_trend_ste_model_nina; fgco2_trend_ste_model_nina(1)=data(end-1);

mk= ['+ooxsshdvv^>'];mks=1;mks2=120;
mkf=['nfnnnfnnnfnf'];
col=[0.9922    0.6667    0.2824];
C = linspecer(7);
trend_unit=100;
leg_text={'F''_{\deltat}','F''_{\deltatK^o}','F''_{\deltatK''}','F''_{\deltat\Deltap^o}','F''_{\deltat\Deltap^''}'};

bsize=0.45; fsize=12; lfsize=10;
fid=figure(1);clf
set(gcf,'paperOrientation','landscape');
set_figure_size(gcf, '', [0.1 0.1 0.1 0.1],1.4,0.8);
subaxis(1,3,1,'ml',0.1,'mr',0.04,'mb',0.12,'sv',0.04,'sh',0.07,'mt',0.08); hold on
    datax=[1:4]; x_tic=[datax];
    b1=bar(x_tic(1),fgco2_trend_model_nino,0.02,'FaceAlpha',0.1); b1.BarWidth=bsize;
    b2=errorbar(x_tic(1),fgco2_trend_model_nino,fgco2_trend_ste_model_nino,'o','MarkerSize',1,'Color','k','linewidth',0.5);
    datax2=x_tic(3); %datax2=[1 idx+1]+0.15;
    b3=bar(datax2,fgco2_trend_comp(1:4,:),0.02,'stacked','FaceAlpha',0.8); b3(1).BarWidth=bsize;
    b3(1).FaceColor=[85 164 242]/255; %light blue
    b3(2).FaceColor=[12 2 96]/255; %dark blue
    b3(3).FaceColor=[247 136 1]/255; %[0    0.0118    0.3569]; %dark blue
    b3(4).FaceColor=[221 45 45]/255; %[0.5843    0.8157    0.9882]; %light blue
    leg_d=legend([b1 b3],leg_text,'Orientation','horizontal','location','northeast');
    leg_d.FontSize=lfsize+1; leg_d.NumColumns=3; leg_d.ItemTokenSize=[10 10];

    ba1=bar(x_tic(2),fgco2_trend_model_nina,0.02,'FaceAlpha',0.1); ba1.BarWidth=bsize; ba1.FaceColor=b1.FaceColor;
    ba2=errorbar(x_tic(2),fgco2_trend_model_nina,fgco2_trend_ste_model_nina,'o','MarkerSize',1,'Color','k','linewidth',0.5);
    datax2=x_tic(4); %datax2=[1 idx+1]+0.15;
    ba3=bar(datax2,fgco2_trend_comp(5:8,:),0.02,'stacked','FaceAlpha',0.8); ba3(1).BarWidth=bsize;
    ba3(1).FaceColor=[85 164 242]/255; %light blue
    ba3(2).FaceColor=[12 2 96]/255; %dark blue
    ba3(3).FaceColor=[247 136 1]/255; %[0    0.0118    0.3569]; %dark blue
    ba3(4).FaceColor=[221 45 45]/255; %[0.5843    0.8157    0.9882]; %light blue

    xlim([0.5 4.5]); xticks(x_tic);set(gca,'xticklabels',''); yticks([-0.08:0.02:0.08]); ylim([-0.03 0.04]);
    labels={{'F''_{\deltat} ';'El Nino '} {'F''_{\deltat}';'La Nina'} {'F''_{\deltat} ';'El Nino'} {'F''_{\deltat}';'La Nina'}};
    h=my_xticklabels(gca,x_tic,labels);
    set(h,'fontsize',fsize);
    ylabel({'CO_2 flux anomaly trend F''_{\deltat} (PgC/yr/century)'},'fontsize',13)
    text(0.66,-0.038,'Multi-model mean'); text(2.84,-0.038,'Decomposition')
    set(gca,'fontsize',fsize); grid off;
    text(0.6,0.0385,'a','fontsize',22);

subaxis(2,3,2);hold on
    %1901-1940
    x0=spco25_nont_nino_all(:,:,2); x0=x0(~isnan(x0)); y0=spco25_temp_nino_all(:,:,2); y0=y0(~isnan(y0));
    data1=x0./y0; data1(data1<=-3)=nan; data1(data1>0)=nan;
    %2051-2090
    x0=spco25_temp_nino_all(:,:,3); x0=x0(~isnan(x0)); y0=spco25_nont_nino_all(:,:,3); y0=y0(~isnan(y0));
    data2=x0./y0; data2(data2<=-3)=nan; data2(data2>0)=nan;
    [h3,p3,ci3,zval3] = ztest(data2,mean(data1),std(data1));
    dx1=0.2;edges1=[-3.1:dx1:0];  xlimf1=[min(edges1) max(edges1)]; xtickf1=[-3:0.5:0]; xlimf1=[-3 0];
    x1=[min(edges1)+dx1/2:dx1:max(edges1)+dx1/2];
    data1_line=pdf_line(x1,data1); data2_line=pdf_line(x1,data2);
    ylimf=[0 0.5]; ytickf=[0:0.1:1];

    h1=histogram(data1,edges1,'normalization','probability');
    h2=histogram(data2,edges1,'normalization','probability');
    h3=plot(x1,data1_line*dx1,'b','linewidth',1.5); h4=plot(x1,data2_line*dx1,'r','linewidth',1.5);
    %ha=patch([-1.2 -0.8 -0.8 -1.2],[0 0 1 1],'y'); ha.FaceAlpha = 0.2;ha.FaceColor=[0.9290 0.6940 0.1250];ha.EdgeAlpha=0.2;
    plot([-1 -1],[-1 1],'k--','linewidth',1);
    ylabel('Probablity (El Nino)','fontsize',fsize); %xlabel('pCO_2 anomaly non-therm:thermal ratio','fontsize',fsize); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1)),' El Ninos)'],['2051-2090 (',num2str(length(data2)),' El Ninos)'],'fontsize',lfsize,'location','northwest');
    l1.ItemTokenSize=[10 8]; l1.Position(2)=l1.Position(2)-0.035;
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf); set(gca,'xticklabels','');
    set(gca,'fontsize',fsize); grid off; % set(l1,'fontsize',10);
    text(-2.97,0.47,'b','fontsize',22);

subaxis(2,3,5);hold on
    %1901-1940
    x0=spco25_nont_nina_all(:,:,2); x0=x0(~isnan(x0)); y0=spco25_temp_nina_all(:,:,2); y0=y0(~isnan(y0));
    data1=x0./y0; data1(data1<=-3)=nan; data1(data1>0)=nan;
    %2051-2090
    x0=spco25_temp_nina_all(:,:,3); x0=x0(~isnan(x0)); y0=spco25_nont_nina_all(:,:,3); y0=y0(~isnan(y0));
    data2=x0./y0; data2(data2<=-3)=nan; data2(data2>0)=nan;
    [h3,p3,ci3,zval3] = ztest(data2,mean(data1),std(data1));
    dx1=0.2;edges1=[-3.1:dx1:0];  xlimf1=[min(edges1) max(edges1)]; xtickf1=[-3:0.5:0]; xlimf1=[-3 0];
    x1=[min(edges1)+dx1/2:dx1:max(edges1)+dx1/2];
    data1_line=pdf_line(x1,data1); data2_line=pdf_line(x1,data2);
    ylimf=[0 0.5]; ytickf=[0:0.1:1];

    h1=histogram(data1,edges1,'normalization','probability');
    h2=histogram(data2,edges1,'normalization','probability');
    h3=plot(x1,data1_line*dx1,'b','linewidth',1.5); h4=plot(x1,data2_line*dx1,'r','linewidth',1.5);
    %ha=patch([-1.2 -0.8 -0.8 -1.2],[0 0 1 1],'y'); ha.FaceAlpha = 0.2;ha.FaceColor=[0.9290 0.6940 0.1250];ha.EdgeAlpha=0.2;
    %hatchfill2(ha,'single','HatchAngle',45,'HatchDensity',15,'HatchColor','k','HatchLineWidth',0.5);
    plot([-1 -1],[-1 1],'k--','linewidth',1);
    ylabel('Probablity (La Nina)','fontsize',fsize); xlabel('pCO_2 anomaly non-therm:thermal ratio (\gamma)','fontsize',fsize); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1)),' La Ninas)'],['2051-2090 (',num2str(length(data2)),' La Ninas)'],'fontsize',lfsize,'location','northwest');
    l1.ItemTokenSize=[10 8]; l1.Position(2)=l1.Position(2)-0.035;
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf);
    set(gca,'fontsize',fsize); grid off; % set(l1,'fontsize',10);
    text(-2.97,0.48,'c','fontsize',22);

%rainfall pdf
    data1=raina_ann_nino_all(:,:,2); data2=raina_ann_nino_all(:,:,3); %%2:1901-1940; 3:2050-2090
    data1=data1(~isnan(data1)); data2=data2(~isnan(data2));
    [h3,p3,ci3,zval3] = ztest(data2,mean(data1),std(data1));
    data1_nino=data1; data2_nino=data2;

    data1=raina_ann_nina_all(:,:,2); data2=raina_ann_nina_all(:,:,3); %%2:1901-1940; 3:2050-2090
    data1=data1(~isnan(data1)); data2=data2(~isnan(data2));
    [h3,p3,ci3,zval3] = ztest(data2,mean(data1),std(data1));
    data1_nina=data1; data2_nina=data2;

    dx1=0.5;edges1=[-6:dx1:6];  xlimf1=[min(edges1) max(edges1)]; xtickf1=[min(edges1):2:max(edges1)]; xlimf1=[-6 6];
    x1=[min(edges1)+dx1/2:dx1:max(edges1)+dx1/2];
    data1_line_nino=pdf_line(x1,data1_nino); data2_line_nino=pdf_line(x1,data2_nino);
    data1_line_nina=pdf_line(x1,data1_nina); data2_line_nina=pdf_line(x1,data2_nina);
    ylimf=[0 0.5]; ytickf=[0:0.1:1];

    subaxis(2,3,3);hold on
    h1=histogram(data1_nino,edges1,'normalization','probability');
    h2=histogram(data2_nino,edges1,'normalization','probability');
    h3=plot(x1,data1_line_nino*dx1,'b','linewidth',1.5); h4=plot(x1,data2_line_nino*dx1,'r','linewidth',1.5);
    plot([0 0],[-1 1],'k--','linewidth',1);
    ylabel({'Probablity (El Nino)'},'fontsize',fsize); %xlabel('Nino3 rainfall anomaly (mm/d)','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1_nino)),' El Ninos)'],['2051-2090 (',num2str(length(data2_nino)),' El Ninos)'],'fontsize',lfsize,'location','northeast');
    l1.ItemTokenSize=[10 8];
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf);
    set(gca,'fontsize',fsize); set(gca,'xticklabels','');grid off;
    text(-5.8,0.47,'d','fontsize',22);

    subaxis(2,3,6);hold on
    h1=histogram(data1_nina,edges1,'normalization','probability');
    h2=histogram(data2_nina,edges1,'normalization','probability');
    h3=plot(x1,data1_line_nina*dx1,'b','linewidth',1.5); h4=plot(x1,data2_line_nina*dx1,'r','linewidth',1.5);
    plot([0 0],[-1 1],'k--','linewidth',1);
    ylabel({'Probablity (La Nina)'},'fontsize',12); xlabel('Nino3 rainfall anomaly (mm/d)','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1_nina)),' La Ninas)'],['2051-2090 (',num2str(length(data2_nina)),' La Ninas)'],'fontsize',lfsize,'location','northeast');
    l1.ItemTokenSize=[10 8];
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf);
    set(gca,'fontsize',fsize); grid off;
    text(-5.8,0.48,'e','fontsize',22);

%print_pdf_liao(fid,['res_figs/fig03_bar_fCO2_trend_pdf_pCO2nont_vs_pCO2t.pdf'])
%print('-dpng','-r500', ['res_figs/fig03_bar_fCO2_trend_pdf_pCO2nont_vs_pCO2t.png']);









