clear; close all; 

variable_ids={'fgco2','spco2','nbp'};
vs=[-3600*24*365.25*1e3,1e6/101325,-3600*24*365.25*1e3];
var_all{1,1}='fgco2';   var_all{1,2}='gC/m2/yr';   var_all{1,3}=['CO2 flux (positive to air)'];
var_all{2,1}='spco2';   var_all{2,2}='uatm';       var_all{2,3}=['sea surface CO2 Partial Pressure'];
var_all{3,1}='nbp';     var_all{3,2}='gC/m2/yr';   var_all{3,3}=['air-land CO2 flux (positive to air)'];
var_all{4,1}='tos';     var_all{4,2}='degree';     var_all{4,3}=['sea surface temperature'];
var_all{5,1}='pr';      var_all{5,2}='mm/dd';      var_all{5,3}=['precipitation'];
var_ano=var_all;
var_mean=var_all;
for i=1:length(variable_ids);
    var_ano(i,3)={[char(var_all(i,3)),' anomaly against moving mean']};
    var_mean(i,3)={[char(var_all(i,3)),' moving mean']};
end

experiment_ids={'historical','ssp245','ssp585'};

%extract_region=[119 301 -31 31]; fin_dir='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member5/data_all_bandpass_bigregion_dicalk';
extract_region=[0 360 -90 90]; fin_dir='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member5/data_all_bandpass_bigregion';

fall=dir([fin_dir,'/*']);
dirFlags = [fall.isdir]; fall = fall(dirFlags); % Extract only those that are directories.
fall=fall(3:end);

for i=1:6%length(fall); %first two are ./ ../ and are removed
    fin=[fall(i).folder,'/',fall(i).name];

    disp([num2str(i),'/',num2str(length(fall)),' ',fall(i).name]);

    fn5_out=[fin,'/all_',fall(i).name,'.his_ssp585_1901_2100_global.nc'];
    fn5_ano_out=[fin,'/all_ano_',fall(i).name,'.his_ssp585_1901_2100_global.nc'];
    fn5_mean_out=[fin,'/all_mean_',fall(i).name,'.his_ssp585_1901_2100_global.nc'];
    disp(['output file ...', fn5_ano_out]);
    disp(['output file ...', fn5_mean_out]);

    variable_id='nbp';
    fin_his =dir([fin,'/',variable_id,'.*',char(experiment_ids(1)),'*.nc']);
    fin_ssp5=dir([fin,'/',variable_id,'.*',char(experiment_ids(3)),'*.nc']);
    if isfile(fn5_out) & isfile(fn5_ano_out);
       fsize=dir(fn5_ano_out); fsize=fsize.bytes/1024/1024/1024;
       %if fsize>4; continue; end
    end
    if length(fin_his)<1 | length(fin_ssp5)<1; continue; end
    %if isfile(fn5_out) & isfile(fn5_ano_out); continue; end
    
    for j=1:length(variable_ids);
        variable_id=char(variable_ids(j));
        disp(['   ',num2str(j),'/',num2str(length(variable_ids)),' ',variable_id]);tic;
        fin_his =dir([fin,'/',variable_id,'.*',char(experiment_ids(1)),'*.nc']);
        fin_ssp2=dir([fin,'/',variable_id,'.*',char(experiment_ids(2)),'*.nc']);
        fin_ssp5=dir([fin,'/',variable_id,'.*',char(experiment_ids(3)),'*.nc']);
        if length(fin_his)==1 & length(fin_ssp5)==1;
           fname_his=[fin,'/',fin_his.name];
           fname_ssp2=[fin,'/',fin_ssp2.name];
	   fname_ssp5=[fin,'/',fin_ssp5.name];

%           disp(['       SSP245']);
%	   [data_all,lon,lat,time]=read_nccmipnd({fname_his,fname_ssp2},variable_id,'',extract_region,vs(j),3);
%           data_all=permute(data_all,[2,1,3]);
%           data_all_ano=compute_anomaly_bandpass(data_all,time,12*15,13);
%           if j==1;
%              create_mom6_1var3d(fn2_out,'lon','lat','time','time_bnds',var_all(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all);
%	      create_mom6_1var3d(fn2_ano_out,'lon','lat','time','time_bnds',var_ano(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_ano);
%           else
%              create_mom6_1var3d_append(fn2_out,'lon','lat','time','time_bnds',var_all(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all);
%	      create_mom6_1var3d_append(fn2_ano_out,'lon','lat','time','time_bnds',var_ano(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_ano);
%           end
%           clear data_all data_all_ano

           disp(['       SSP585']);
           [data_all,lon,lat,time]=read_nccmipnd({fname_his,fname_ssp5},variable_id,'',extract_region,vs(j),3);
           data_all=permute(data_all,[2,1,3]);
           [data_all_ano, data_all_mean]=compute_anomaly_bandpass(data_all,time,12*15,13);
           if j==1;
              create_mom6_1var3d(fn5_out,'lon','lat','time','time_bnds',var_all(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all);
              create_mom6_1var3d(fn5_ano_out,'lon','lat','time','time_bnds',var_ano(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_ano);
              create_mom6_1var3d(fn5_mean_out,'lon','lat','time','time_bnds',var_mean(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_mean);
           else
              create_mom6_1var3d_append(fn5_out,'lon','lat','time','time_bnds',var_all(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all);
              create_mom6_1var3d_append(fn5_ano_out,'lon','lat','time','time_bnds',var_ano(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_ano);
              create_mom6_1var3d_append(fn5_mean_out,'lon','lat','time','time_bnds',var_mean(j,:),lon(1,:),lat(:,1),time-datenum(1900,1,1),[],data_all_mean);
           end
           clear data_all data_all_ano data_all_mean

        end
        toc;
     end
     ncwriteatt(fn5_out,'/','data description','this is CMIP6 results from 1901-2100');
     ncwriteatt(fn5_ano_out,'/','data description','this is CMIP6 results anomaly against moving 30 years mean from 1901-2100');
     ncwriteatt(fn5_mean_out,'/','data description','this is CMIP6 results moving mean including trend and decadal variability from 1901-2100');
end


