clear
%close all

%model_name='CESM2-WACCM.r1i1p1f1';
model_name='GFDL-ESM4.r1i1p1f1';

fdir=['/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member5/data_all_bandpass_smallregion_dicalk/',model_name,'/'];
fin_ano=[fdir,'all_ano_',model_name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
fin_decomp=[fdir,'all_decomposition_fgco2_spco2_',model_name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
fin_mean=[fdir,'all_mean_',model_name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
fin_allpr=[fdir,'/all_',model_name,'.his_ssp585_1901_2100_tropical_pacific.nc'];

fgco2_region=[140 285 -15 15]; trange=[datenum(1901,1,1) datenum(2101,1,1)];
nino3_region=[210 270 -5 5];
sst_mons=[12 14]; fgco2_mons=[7 18];

[sst,lon,lat,time_sst]=read_1nccmipnd(fin_ano,'tos',nino3_region,trange,1,1);
[rain,lon,lat,time_sst]=read_1nccmipnd(fin_allpr,'pr',nino3_region,trange,1,1);
[sst_ann,time_ann]=compute_annmean(sst,time_sst,sst_mons);
[rain_ann,time_ann]=compute_annmean(rain,time_sst,sst_mons);
ninos=compute_strength(sst_ann,rain_ann);

[fgco2_ano,lon,lat,time_fgco2]=read_1nccmipnd(fin_decomp,'fgco2_ano',fgco2_region,trange,1,3);
[fgco2_mean,lon,lat,time_fgco2]=read_1nccmipnd(fin_mean,'fgco2',fgco2_region,trange,1,3);
[kw_ano,lon,lat,time_fgco2]=read_1nccmipnd(fin_decomp,'kw_ano',fgco2_region,trange,1,3);
[kw_mean,lon,lat,time_fgco2]=read_1nccmipnd(fin_decomp,'kw_m',fgco2_region,trange,1,3);
[dpco2_ano,lon,lat,time_fgco2]=read_1nccmipnd(fin_decomp,'dpco2_ano',fgco2_region,trange,1,3);
[dpco2_mean,lon,lat,time_fgco2]=read_1nccmipnd(fin_decomp,'dpco2_m',fgco2_region,trange,1,3);

inds=find(lon>=140 & lon<=285 & lat>=-15 & lat<=15);
inds_nan=find(lon>=276 & lon<=285 & lat>=9 & lat<=15);
var_all={'fgco2_ano','fgco2_mean','dpco2_ano','dpco2_mean','kw_ano','kw_mean'};
for v=1:length(var_all)
    var=char(var_all(v));
    eval(['data=',var,';']);
    for l=1:length(time_fgco2)
        data2d=squeeze(data(:,:,l));
        data2d(inds_nan)=nan; 
        data_mons(v,l)=nanmean(data2d(inds));
    end
    [data_ann(v,:),time_ann]=compute_annmean(data_mons(v,:),time_fgco2,fgco2_mons);
end

%leg_all(1,:)={'F','F''','F^o','F''-Nino','F^o-Nino'};leg_all(2,:)={'\Deltap','\Deltap''','\Deltap^o','\Deltap''-Nino','\Deltap^o-Nino'};
%leg_all(3,:)={'K','K''','K^o','K''-Nino','K^o-Nino'};
leg_all(1,:)={'F','F''','F^0','El Nino','La Nina'};
leg_all(2,:)={'\Deltap','\Deltap''','\Deltap^0','El Nino','La Nina'};
leg_all(3,:)={'K','K''','K^0','El Nino','La Nina'};
ytext=[19,38,0.55]; laball='abc';
col_blue=[0 0.4470 0.7410]; col_red=[0.8500 0.3250 0.0980]; col_yel=[0.9290 0.6940 0.1250]; col_4=[0.3010 0.7450 0.9330];
fid=figure(1);clf
set(gcf,'paperOrientation','landscape');
set_figure_size(gcf, '', [0.1 0.1 0.1 0.1],1.2,0.8);
for i=1:3;
    subaxis(3,1,i,'ml',0.08,'mr',0.04,'mb',0.1,'sv',0.035,'mt',0.05);hold on; 
    i1=1+(i-1)*2; i2=2+(i-1)*2;
    h1=plot(time_ann,data_ann(i1,:)+data_ann(i2,:),'color','k','Linewidth',0.6); %fgco2_all=fgco2_ano+fgco2_mean
    h2=plot(time_ann,data_ann(i1,:),'color',col_4,'LineWidth',1);
    h3=plot(time_ann,data_ann(i2,:),'color',col_yel,'LineWidth',2);
    s1_nino=plot(time_ann(ninos>=1),  data_ann(i1,ninos>=1),'.','color',col_red,'MarkerSize',20);
    s1_nina=plot(time_ann(ninos<=-1), data_ann(i1,ninos<=-1),'^','MarkerFaceColor',col_blue,'MarkerEdgeColor',col_blue,'MarkerSize',5);
    s2_nino=plot(time_ann(ninos>=1),  data_ann(i2,ninos>=1),'.','color',col_red,'MarkerSize',20);
    s2_nina=plot(time_ann(ninos<=-1), data_ann(i2,ninos<=-1),'^','MarkerFaceColor',col_blue,'MarkerEdgeColor',col_blue,'MarkerSize',5);
    xticks([datenum(1900,1:12*20:12*200,15)]);xlim([datenum(1900,1,1) datenum(2091,1,1)]);
    datetick('x','yyyy','keeplimits','keepticks');
    l1=legend([h1 h2 h3 s1_nino s1_nina],leg_all(i,:),'location','northeast');l1.FontSize=12; l1.NumColumns=2; l1.ItemTokenSize=[10 0.2];
    if i<=2; set(gca,'xticklabels',''); else; xlabel('Years'); end
    if i==1; yticks([-10:10:20]);ylim([-5 20]); ylabel(['CO_2 flux (F, gC/m^{2}/yr)']); end
    if i==2; yticks([-60:20:60]);ylim([-10 40]); ylabel(['\Deltap (uatm)']);  end
    if i==3; yticks([-1:0.25:1]);ylim([-0.05 0.6]); ylabel(['K (F/\Deltap)']); l1.Location='east';end
    set(gca,'fontsize',12,'FontName','Times');
    text(time_ann(1)-680,ytext(i),laball(i),'fontsize',22);
end
print_pdf_liao(fid,['res_figs/figs02_sample_ano_ann_time_series.pdf'])
print('-dpng','-r500', ['res_figs/figs02_sample_ano_ann_time_series.png']);



