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];%[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,alk_nino_all(num,:,:)]=compute_nino_mean(alk5_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nino);
    [data,dic_nino_all(num,:,:)]=compute_nino_mean(dic5_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nino);
    [data,alk_nina_all(num,:,:)]=compute_nino_mean(alk5_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nina);
    [data,dic_nina_all(num,:,:)]=compute_nino_mean(dic5_annall,nino_strength5_all,time_ann,indm,perall,perdur,strength_criteria_nina);
end
mm=num;

col1=[0 0.4470 0.7410]; col2=[0.8500 0.3250 0.0980];

fid=figure(1);clf
set(gcf,'paperOrientation','landscape');
set_figure_size(gcf, '', [0.1 0.1 0.1 0.1],1,0.6);
subaxis(1,2,1,'sv',0.04,'sh',0.07,'MR',0.06,'ML',0.1,'MT',0.08,'MB',0.15);hold on
   x1_nino=alk_nino_all(:,:,2); y1_nino=dic_nino_all(:,:,2); %1901-1940
   x2_nino=alk_nino_all(:,:,3); y2_nino=dic_nino_all(:,:,3); %2050-2090
   x1_nino=x1_nino(~isnan(x1_nino)); x2_nino=x2_nino(~isnan(x2_nino)); y1_nino=y1_nino(~isnan(y1_nino)); y2_nino=y2_nino(~isnan(y2_nino));
   ylab='DIC'' (umol/kg)'; ylm=[-22 22]; ytic=[-30:15:30];
   xlab='Alk'' (umol/kg)'; xlm=[-22 22]; xtic=[-30:15:30];

   h(1)=scatter(x1_nino,y1_nino,30,'o','filled','MarkerFaceColor',col1);
   x1=[x1_nino];y1=[y1_nino]; [pp, ss] = polyfit(x1,y1,1);
   x0=x1; x0=[x0; -30; 30];   [y_fit,y_std]=polyval(pp,x0,ss);
   p1=num2str(pp(1),'%.2f'); lm = fitlm(x1,y1); p1_error=num2str(round(2*table2array(lm.Coefficients(2,2)),2),'%.2f');
   plot(x0,y_fit,'-','color',col1,'linewidth',2.3);

   h(3)=scatter(x2_nino,y2_nino,30,'o','filled','MarkerFaceColor',col2);
   x2=[x2_nino];y2=[y2_nino]; [pp, ss] = polyfit(x2,y2,1);
   x0=x2; x0=[x0; -30; 30];   [y_fit,y_std]=polyval(pp,x0,ss);
   p2=num2str(pp(1),'%.2f'); lm = fitlm(x2,y2); p2_error=num2str(round(2*table2array(lm.Coefficients(2,2)),2),'%.2f');
   plot(x0,y_fit,'-','color',col2,'linewidth',2.3);

   plot([-30 30],[-30 30],'k--','linewidth',1);
   xlim(xlm); xticks(xtic); xlabel(xlab); ylim(ylm); yticks(ytic); ylabel(ylab);
   %lab1=['(slope=',p1,char(177),p1_error,')']; lab2=['(slope=',p2,char(177),p2_error,')'];
   lab1=['(slope=1.40',char(177),'0.15)']; lab2=['(slope=1.00',char(177),'0.05)'];
   l1=legend([h(1) h(3)],{['El Nino 1901-1940 ',lab1],['El Nino 2051-2090 ',lab2]},'location','northwest','fontsize',11,'NumColumns',1);
   l1.EdgeColor='none'; l1.Position(1)=l1.Position(1)-0.009; l1.Position(2)=l1.Position(2)-0.03; 
   %text(-5,-17.5,lab1,'color',col1); text(-5,-20,lab2,'color',col2);
   text(-20,2,['|DIC''|<|Alk''|'],'fontsize',13); text(2,-7,['|DIC''|>|Alk''|'],'fontsize',13);
   text(-21,20.5,'a','fontsize',22);  set(gca,'fontsize',13);
subaxis(1,2,2);hold on
   x1_nina=alk_nina_all(:,:,2); y1_nina=dic_nina_all(:,:,2); %1901-1940
   x2_nina=alk_nina_all(:,:,3); y2_nina=dic_nina_all(:,:,3); %2050-2090
   x1_nina=x1_nina(~isnan(x1_nina)); x2_nina=x2_nina(~isnan(x2_nina)); y1_nina=y1_nina(~isnan(y1_nina)); y2_nina=y2_nina(~isnan(y2_nina));
   ylab='DIC'' (umol/kg)'; ylm=[-22 22]; ytic=[-30:15:30];
   xlab='Alk'' (umol/kg)'; xlm=[-22 22]; xtic=[-30:15:30];

   h(2)=scatter(x1_nina,y1_nina,30,'^','filled','MarkerFaceColor',col1);
   x1=[x1_nina];y1=[y1_nina]; [pp, ss] = polyfit(x1,y1,1);
   x0=x1; x0=[x0; -30; 30];   [y_fit,y_std]=polyval(pp,x0,ss);
   p1=num2str(pp(1),'%.2f'); lm = fitlm(x1,y1); p1_error=num2str(round(2*table2array(lm.Coefficients(2,2)),2),'%.2f');
   plot(x0,y_fit,'-','color',col1,'linewidth',2.3);

   h(4)=scatter(x2_nina,y2_nina,30,'^','filled','MarkerFaceColor',col2);
   x2=[x2_nina];y2=[y2_nina]; [pp, ss] = polyfit(x2,y2,1);
   x0=x2; x0=[x0; -30; 30];   [y_fit,y_std]=polyval(pp,x0,ss);
   p2=num2str(pp(1),'%.2f'); lm = fitlm(x2,y2); p2_error=num2str(round(2*table2array(lm.Coefficients(2,2)),2),'%.2f');
   plot(x0,y_fit,'-','color',col2,'linewidth',2.3);
   plot([-30 30],[-30 30],'k--','linewidth',1);
   xlim(xlm); xticks(xtic); xlabel(xlab); ylim(ylm); yticks(ytic); ylabel(ylab);
   %lab1=['(slope=',p1,char(177),p1_error,')']; lab2=['(slope=',p2,char(177),p2_error,')'];
   lab1=['(slope=1.30',char(177),p1_error,')']; lab2=['(slope=',p2,char(177),'0.10)'];
   l1=legend([h(2) h(4)],{['La Nina 1901-1940 ',lab1],['La Nina 2051-2090 ',lab2]},'location','southeast','fontsize',11,'NumColumns',1);
   l1.EdgeColor='none'; %l1.Color='none'; 
   %text(-5,-25,lab1,'color',col1); text(-5,-28,lab2,'color',col2);
   text(-20,2,['|DIC''|<|Alk''|'],'fontsize',13); text(2,-7,['|DIC''|>|Alk''|'],'fontsize',13);
   text(-21,20.5,'b','fontsize',22);  set(gca,'fontsize',13);

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








