clear; close all; tic;

fgco2_region=[120 300 -30 30]; %read_1nccmipnd could remove atlantic flux;[160 280 -15 15];
fgco2_regions=[140 285 -15 15];
nino3_region=[210 270 -5 5];
nino34_region=[190 240 -5 5];
trange=[datenum(1901,1,1) datenum(2101,1,1)];

fin_obs={'obs_ls.nc','obs_rb.nc'}; title_obs={'pCO2-based SOMFFN','pCO2-based MLS'};
fin_dir='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member5/data_all_bandpass_smallregion_dicalk';
fall=dir([fin_dir,'/*']);
dirFlags = [fall.isdir]; fall = fall(dirFlags); % Extract only those that are directories.
fall=fall(3:end); num=0;

%source_ids_rm={'MPI-ESM1-2-HR','MPI-ESM1-2-LR','BCC-CSM2-MR','INM-CM4-8','INM-CM5-0'};
source_ids_rm={'MPI-ESM1-2-HR','MPI-ESM1-2-LR','BCC-CSM2-MR','INM-CM4-8','INM-CM5-0','CanESM5','MIROC-ES2L','NorESM2-LM','UKESM1-0-LL'};
for i=1:length(fall); %first two are ./ ../ and are removed
    fin=[fall(i).folder,'/',fall(i).name];
    [source_member_id0]=strsplit(fall(i).name,'.');
    var_name='dpco2'; %'dpco2';%'talkos';
    fn_var=dir([fin,'/',var_name(1:end-2),'*.nc']);;fn_varos=dir([fin,'/',var_name,'*.nc']);
    if length(fn_var)<=2 & length(fn_varos)<=2; continue; end
    if contains(source_member_id0(1),source_ids_rm); continue; end
 
    num=num+1;
    [source_member_ids(num,:)]=strsplit(fall(i).name,'.');

    fin=[fall(i).folder,'/',fall(i).name];
    toc;disp([num2str(i),'/',num2str(length(fall)),'   ',fin]);tic;
    fn2_ano=[fin,'/all_ano_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_ano=[fin,'/all_ano_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn2 =[fin,'/all_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5 =[fin,'/all_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn2_dep=[fin,'/all_decomposition_fgco2_spco2_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_dep=[fin,'/all_decomposition_fgco2_spco2_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];

    [sst5(num,:),lon,lat,time_sst]=read_1nccmipnd(fn5_ano,'tos',nino3_region,trange,1,1);
    [rains5(num,:),lon,lat,time_rain]=read_1nccmipnd(fn5,'pr',nino3_region,trange,1,1);

    [spco25_temp(num,:),lon,lat,time_rain]=read_1nccmipnd(fn5_dep,'spco2_temp',fgco2_regions,trange,1,1);
    [spco25_nont(num,:),lon,lat,time_rain]=read_1nccmipnd(fn5_dep,'spco2_nont',fgco2_regions,trange,1,1);
    [alk5(num,:),lon,lat,time_rain]=read_1nccmipnd(fn5_ano,'talkos',fgco2_regions,trange,1,1);
    [dic5(num,:),lon,lat,time_rain]=read_1nccmipnd(fn5_ano,'dissicos',fgco2_regions,trange,1,1);

    sst_mons=[12 14];
    fgco2_mons=[7 18];
    [sst5_annall(num,:),time_ann]=compute_annmean(sst5(num,:),time_sst,sst_mons);
    [rains5_annall(num,:),time_ann]=compute_annmean(rains5(num,:),time_rain,sst_mons);
    [spco25_temp_annall(num,:),time_ann]=compute_annmean(spco25_temp(num,:),time_rain,fgco2_mons);
    [spco25_nont_annall(num,:),time_ann]=compute_annmean(spco25_nont(num,:),time_rain,fgco2_mons);
    [alk5_annall(num,:),time_ann]=compute_annmean(alk5(num,:),time_rain,fgco2_mons);
    [dic5_annall(num,:),time_ann]=compute_annmean(dic5(num,:),time_rain,fgco2_mons);

    rain5_skewness(num)=skewness(rains5(num,:));
    nino_strength5_all(num,:)=compute_strength(sst5_annall(num,:),rains5_annall(num,:));

end

[source_ids_all,ia,ic]=unique(source_member_ids(:,1));
perall=[1900 1901 2051]; perdur=[190 40 40];
strength_criteria_nino=1;
strength_criteria_nina=-1;
num=0;
for i=1:length(source_ids_all)
    source=char(source_ids_all(i));
    if contains(source,source_ids_rm); continue; end
    num=num+1;
    source_ids{num}=source;
    indm=find(strcmp(source_member_ids(:,1),char(source_ids(num))));
    if length(indm)>=3; indm=indm(1:3); end

    [data,spco25_temp_nino_all(num,:,:)]=compute_nino_mean(spco25_temp_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nino);
    [data,spco25_nont_nino_all(num,:,:)]=compute_nino_mean(spco25_nont_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nino);

    [data,spco25_temp_nina_all(num,:,:)]=compute_nino_mean(spco25_temp_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nina);
    [data,spco25_nont_nina_all(num,:,:)]=compute_nino_mean(spco25_nont_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nina);
end
mm=num;

save('matres_spco2tnt_pdf_nino_nina.mat','spco25_nont_nino_all','spco25_temp_nino_all','spco25_nont_nina_all','spco25_temp_nina_all');

mk= ['+oox*sshdvv^^>'];mks=1;mks2=50; mkf=['nfnnnfnnnfnfnn'];
C = linspecer(4);
col=[0 0.4470 0.7410; 0.8500 0.3250 0.0980];

%refined figure
fid=figure(1);clf
    set(gcf,'paperOrientation','landscape');
    set_figure_size(gcf, '', [0.1 0.1 0.1 0.1],1,0.85);
    %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:dx1:0];  xlimf1=[min(edges1) max(edges1)]; xtickf1=[min(edges1):0.5:max(edges1)]; 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];

    subaxis(2,2,1,'sv',0.08,'sh',0.09,'MR',0.035,'ML',0.09,'MT',0.07,'MB',0.17);hold on
    h1=histogram(data1,edges1,'normalization','probability'); h1.FaceColor=col(1,:);
    h2=histogram(data2,edges1,'normalization','probability'); h2.FaceColor=col(2,:);
    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);
    ylabel('Probablity','fontsize',12); xlabel('ratio of pCO2-nont vs pCO2-thermal','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1)),' El Ninos)'],['2051-2090 (',num2str(length(data2)),' El Ninos)'],'fontsize',10,'location','northwest');
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf);
    set(gca,'fontsize',14); grid off; % set(l1,'fontsize',10);

    subaxis(2,2,2);hold on
    h6=cdfplot([-3; data1; 0]); h6.LineWidth=3;
    h7=cdfplot([-3; data2; 0]); h7.LineWidth=3;
    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);
    ylabel('Probablity','fontsize',12); xlabel('ratio of pCO2-nont vs pCO2-thermal','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h6 h7],['1901-1940 (',num2str(length(data1)),' El Ninos)'],['2051-2090 (',num2str(length(data2)),' El Ninos)'],'fontsize',10,'location','northwest');
    xlabel('Ratio of pCO2-nont vs pCO2-thermal','fontsize',12); ylabel('Cumulative Probablity','fontsize',12);
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1); yticks([0:0.25:1]); ylim([0 1]); %yticks(ytickf);
    set(gca,'fontsize',14); title(''); grid off

 %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:dx1:0];  xlimf1=[min(edges1) max(edges1)]; xtickf1=[min(edges1):0.5:max(edges1)]; 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];

    subaxis(2,2,3);hold on
    h1=histogram(data1,edges1,'normalization','probability'); h1.FaceColor=col(1,:);
    h2=histogram(data2,edges1,'normalization','probability'); h2.FaceColor=col(2,:);
    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);
    ylabel('Probablity','fontsize',12); xlabel('ratio of pCO2-nont vs pCO2-thermal','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h1 h2],['1901-1940 (',num2str(length(data1)),' La Ninas)'],['2051-2090 (',num2str(length(data2)),' La Ninas)'],'fontsize',10,'location','northwest');
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1);ylim(ylimf);yticks(ytickf);
    set(gca,'fontsize',14); grid off; % set(l1,'fontsize',10);

    subaxis(2,2,4);hold on
    h6=cdfplot([-3; data1; 0]); h6.LineWidth=3;
    h7=cdfplot([-3; data2; 0]); h7.LineWidth=3;
    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);
    ylabel('Probablity','fontsize',12); xlabel('ratio of pCO2-nont vs pCO2-thermal','fontsize',12); %set(gca,'xticklabels','')
    l1=legend([h6 h7],['1901-1940 (',num2str(length(data1)),' La Ninas)'],['2051-2090 (',num2str(length(data2)),' La Ninas)'],'fontsize',10,'location','northwest');
    xlabel('Ratio of pCO2-nont vs pCO2-thermal','fontsize',12); ylabel('Cumulative Probablity','fontsize',12);
    set(gca,'xtick',xtickf1,'tickdir','out');xlim(xlimf1); yticks([0:0.25:1]); ylim([0 1]); %yticks(ytickf);
    set(gca,'fontsize',14); title(''); grid off









function [corr]=compute_corr(fco22,sst2,fco25,sst5,t);

for l=1:3;
    if l==1; ind=find(t>=datenum(1901,1,1) & t<=datenum(1941,1,1)); end
    if l==2; ind=find(t>=datenum(1980,1,1) & t<=datenum(2020,1,1)); end
    if l==3; ind=find(t>=datenum(2061,1,1) & t<=datenum(2101,1,1)); end

    sst0=sst2(ind);
    f0=fco22(ind);
    r = corrcoef(sst0,f0);
    corr_2(l)=r(2,1);

    sst0=sst5(ind);
    f0=fco25(ind);
    r = corrcoef(sst0,f0);
    corr_5(l)=r(2,1);

end
corr(1)=corr_5(1);
corr(2:3)=corr_2(2:3);
corr(4)=corr_5(3);
end




