clear; close all; 

variable_ids={'pp'};
vs=[1,1];
var_all{1,1}='pp';     var_all{1,2}='m/s';        var_all{1,3}=['eastward sea surface wind'];
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 including trend+decadal']};
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';
%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=[139 286 -16 16]; 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);

for i=26 %1:length(fall); %first two are ./ ../ and are removed
    fin=[fall(i).folder,'/',fall(i).name];
    %fn_ua=dir([fin,'/ua.*.nc']); if length(fn_ua)>=2; variable_ids(1)={'ua'}; end
    %fn_va=dir([fin,'/va.*.nc']); if length(fn_va)>=2; variable_ids(2)={'va'}; end
    %if length(fn_dpco2)<=2 | length(fn_alk)<=2 | length(fn_dic)<=2; continue; end
    disp([num2str(i),'/',num2str(length(fall)),' ',fall(i).name]);

    fn2_out=[fin,'/all_pp_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn2_ano_out=[fin,'/all_pp_ano_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn2_mean_out=[fin,'/all_pp_mean_',fall(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_out=[fin,'/all_pp_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn5_ano_out=[fin,'/all_pp_ano_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn5_mean_out=[fin,'/all_pp_mean_',fall(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    disp(['output file ...', fn2_ano_out]);
    disp(['output file ...', fn5_ano_out]);
    disp(['output file ...', fn2_mean_out]);
    disp(['output file ...', fn5_mean_out]);

    %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']);
        fin_histos =dir([fin,'/tos.*',char(experiment_ids(1)),'*.nc']);
        fin_ssp5tos =dir([fin,'/tos.*',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];
	   fname_histos=[fin,'/',fin_histos.name];
	   fname_ssp5tos=[fin,'/',fin_ssp5tos.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,data_all_mean]=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);
%	      create_mom6_1var3d(fn2_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(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);
%	      create_mom6_1var3d_append(fn2_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

           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

           if j==length(variable_ids)
              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, decdadal variability from 1901-2100');
           end
        end
        toc;
     end
end


