clear; close all; tic;

%fgco2_region=[120 300 -30 30]; %read_1nccmipnd could remove atlantic flux;[160 280 -15 15];
fgco2_region=[140 285 -15 15];
nino3_region=[210 270 -5 5];%[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';
%fin_dir='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member2/data_all_bandpass_bigregion';
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','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='fgco2'; %'nbp'; %'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'];
    fn5_ano_global=[fin,'/all_ano_',fall(i).name,'.his_ssp585_1901_2100_global.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'];
    fn5_global =[fin,'/all_',fall(i).name,'.his_ssp585_1901_2100_global.nc'];
    fn2_decom=[fin,'/all_decomposition_fgco2_spco2_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_decom=[fin,'/all_decomposition_fgco2_spco2_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    sst_mons=[12 14];
    fgco2_mons=[7 18];
    a=5;
    if a==2;
       ssplab='SSP 245';
       fn_decom=fn2_decom; fn_ano=fn2_ano; fn=fn2;
    elseif a==5;
       ssplab='SSP 585';
       fn_decom=fn5_decom; fn_ano=fn5_ano; fn=fn5; fn_ano_global=fn5_ano_global; fn_global=fn5_global;
    end

    [sst(num,:),lon,lat,time_sst]=read_1nccmipnd(fn_ano,'tos',nino3_region,trange,1,1);
    [rain(num,:),lon,lat,time_rain]=read_1nccmipnd(fn,'pr',nino3_region,trange,1,1);
    [raina(num,:),lon,lat,time_rain]=read_1nccmipnd(fn_ano,'pr',nino3_region,trange,1,1);
    [fgco2(num,:),lon,lat,time_fgco2]=read_1nccmipnd(fn_decom,'fgco2_ano',fgco2_region,trange,1,1);

    [sst_ann(num,:),time_ann]=compute_annmean(sst(num,:),time_sst,sst_mons);
    [rain_ann(num,:),time_ann]=compute_annmean(rain(num,:),time_rain,sst_mons);
    [raina_ann(num,:),time_ann]=compute_annmean(raina(num,:),time_rain,sst_mons);
    [fgco2_ann(num,:),time_ann]=compute_annmean(fgco2(num,:),time_fgco2,fgco2_mons);

    rain_skewness(num)=skewness(rain(num,:));
    nino_strength(num,:)=compute_strength(sst_ann(num,:),rain_ann(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,raina_ann_nino_all(num,:,:)]=compute_nino_mean(raina_ann,nino_strength,time_ann,indm,perall,perdur,strength_criteria_nino);
    [data,fgco2_ann_nino_all(num,:,:)]=compute_nino_mean(fgco2_ann,nino_strength,time_ann,indm,perall,perdur,strength_criteria_nino);
    [data,raina_ann_nina_all(num,:,:)]=compute_nino_mean(raina_ann,nino_strength,time_ann,indm,perall,perdur,strength_criteria_nina);
    [data,fgco2_ann_nina_all(num,:,:)]=compute_nino_mean(fgco2_ann,nino_strength,time_ann,indm,perall,perdur,strength_criteria_nina);

end
mm=num;
disp(['load obs Landschutzer and Rodenbeck data ...']);
%yr_obs=[1982 1986 1991 1994 1997 2002 2006 2009 2015];
yr_obs_nino=[1982 1986 1991      1997           2009 2015]; %obs nino year compute by plot_obs_sst_rain_nino_strength.m
yr_obs_nina=[1984 1988 1996 1998 1999 2005 2007 2010]; %obs nina year computed by plot_obs_sst_rain_nino_strength.m

for i=1:2
    fn_in=char(fin_obs(i));
    [fgco2o,lon_obs,lat_obs,time_obs]=read_1nccmipnd(fn_in,'fco2_ano',fgco2_region,trange,1,1);
    [fgco2o,time_obs_ann]=compute_annmean(fgco2o,time_obs,fgco2_mons);

    [yrall,mo,da]=datevec(time_obs_ann); clear mo da
    [tf,idx_nino] = ismember(yr_obs_nino,yrall);
    fgco2o_nino_all(i,:)=fgco2o(idx_nino);
    [tf,idx_nina] = ismember(yr_obs_nina,yrall);
    fgco2o_nina_all(i,:)=fgco2o(idx_nina);
end

save('matres_rainfall_pdf_nino_nina.mat','raina_ann_nino_all','raina_ann_nina_all');

 fid=figure(1);clf
    set(gcf,'paperOrientation','landscape');
    set_figure_size(gcf, '', [0.1 0.1 0.1 0.1],1.1,0.5);

    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.35]; ytickf=[0:0.1:1]; fsize=12; lfsize=9;

    subaxis(1,2,1,'sv',0.08,'sh',0.06,'MR',0.035,'ML',0.09,'MT',0.07,'MB',0.17);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({'El Nino';'Probablity'},'fontsize',12); %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.317,'a','fontsize',22);

    subaxis(1,2,2);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({'La Nina';'Probablity'},'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.317,'b','fontsize',22);













