clear; close all; tic; 

dir_home='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member5/data_all_bandpass_smallregion_dicalk';
dir_suball=dir([dir_home,'/*']);
dirFlags = [dir_suball.isdir]; dir_suball = dir_suball(dirFlags); % Extract only those that are directories.
dir_suball=dir_suball(3:end);

for i=1:length(dir_suball); %first two are ./ ../ and are removed
    [source_member_ids(i,:)]=strsplit(dir_suball(i).name,'.');

    dir_subn=[dir_suball(i).folder,'/',dir_suball(i).name];
    fn_dpco2=dir([dir_subn,'/*dpco2*.nc']);
    fn_spco2=dir([dir_subn,'/*spco2*.nc']);
    if length(fn_dpco2)<=2 | length(fn_spco2)<=2 ; continue; end

    toc;disp([num2str(i),'/',num2str(length(dir_suball)),' ',dir_suball(i).name]);tic;

    fn2_allin=[dir_subn,'/all_',dir_suball(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_allin=[dir_subn,'/all_',dir_suball(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn2_anoin=[dir_subn,'/all_ano_',dir_suball(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_anoin=[dir_subn,'/all_ano_',dir_suball(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn2_meanin=[dir_subn,'/all_mean_',dir_suball(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_meanin=[dir_subn,'/all_mean_',dir_suball(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    fn2_out=[dir_subn,'/all_apco2_',dir_suball(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_out=[dir_subn,'/all_apco2_',dir_suball(i).name,'.his_ssp585_1901_2100_tropical_pacific.nc'];
    %if isfile(fn5_out) & isfile(fn2_out); continue; end

    for j=1:2
        if j==1; fn_allin=fn2_allin; fn_anoin=fn2_anoin; fn_meanin=fn2_meanin; fn_out=fn2_out; end 
        if j==2; fn_allin=fn5_allin; fn_anoin=fn5_anoin; fn_meanin=fn5_meanin; fn_out=fn5_out; end 
        %if ~isfile(fn_out); continue; end
        %fstat=dir(fn_out);fage=(datenum(clock)-datenum(fstat.date))*24*3600;
        %if fage<3600*24; continue; end
        %if j==1 & isfile(fn2_out); continue; end; if j==2 & isfile(fn5_out); continue; end

        disp(['creating file ... ',fn_out]);
        lon=ncread(fn_allin,'lon'); lat=ncread(fn_allin,'lat'); time=ncread(fn_allin,'time')+datenum(1900,1,1);

        spco2=ncread(fn_allin,'spco2');
        dpco2=ncread(fn_allin,'dpco2');
        apco2=spco2-dpco2;
        [apco2_ano, apco2_m]=compute_anomaly_bandpass(apco2,time,12*15,13);

        spco2_m=ncread(fn_meanin,'spco2');
        dpco2_m=ncread(fn_meanin,'dpco2');

        spco2_ano=ncread(fn_anoin,'spco2');
        dpco2_ano=ncread(fn_anoin,'dpco2');

        [var_all]=load_var_all(1);
        %write fgco2, fgco2_ano, dpco2, dpco2_ano, spco2, spco2_ano, kw, kw_ano, fgco2_kw, fgco2_dpco2, fgco2_ct, fgco2_resi
        create_mom6_1var3d(fn_out,'lon','lat','time','time_bnds',var_all(1,:),lon,lat,time-datenum(1900,1,1),[],dpco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(2,:),lon,lat,time-datenum(1900,1,1),[],dpco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(3,:),lon,lat,time-datenum(1900,1,1),[],dpco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(4,:),lon,lat,time-datenum(1900,1,1),[],spco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(5,:),lon,lat,time-datenum(1900,1,1),[],spco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(6,:),lon,lat,time-datenum(1900,1,1),[],spco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(7,:),lon,lat,time-datenum(1900,1,1),[],apco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(8,:),lon,lat,time-datenum(1900,1,1),[],apco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(9,:),lon,lat,time-datenum(1900,1,1),[],apco2_m);

        ncwriteatt(fn_out,'/','data description','this is CMIP6 results from 1901-2100');
        ncwriteatt(fn_out,'/','budget1',['kw=fgco2/dpco2']);
        ncwriteatt(fn_out,'/','budget2',['fgco2_ano=fgco2_kw+fgco2_dpco2+fgco2_ct+fgco2_resi']);
        ncwriteatt(fn_out,'/','budget3',['spco2_ano=spco2_temp+spco2_salt+spco2_dic+dpco2_alk; spco2_temp=0.0423*spco2_m.*tos_ano; spco2_salt=(spco2_m/salt_m)*salt_ano']);

    end
end

function [var_all]=load_var_all(load_ind)
   if load_ind==1;
      %write fgco2, fgco2_ano, dpco2, dpco2_ano, spco2, spco2_ano, kw, kw_ano, fgco2_kw, fgco2_dpco2, fgco2_ct, fgco2_resi
      var_all{1,1}='dpco2';        var_all{1,2}='uatm';             var_all{1,3}=['Delta CO2 (ocean-atm) Partial Pressure'];
      var_all{2,1}='dpco2_ano';    var_all{2,2}='uatm';             var_all{2,3}=['Delta CO2 (ocean-atm) Partial Pressure anomaly'];
      var_all{3,1}='dpco2_m';      var_all{3,2}='uatm';             var_all{3,3}=['Delta CO2 (ocean-atm) Partial Pressure mean'];
      var_all{4,1}='spco2';        var_all{4,2}='uatm';             var_all{4,3}=['ocean pco2'];
      var_all{5,1}='spco2_ano';    var_all{5,2}='uatm';             var_all{5,3}=['ocean pco2 anomaly'];
      var_all{6,1}='spco2_m';      var_all{6,2}='uatm';             var_all{6,3}=['ocean pco2 mean'];
      var_all{7,1}='apco2';        var_all{7,2}='uatm';             var_all{7,3}=['air pco2'];
      var_all{8,1}='apco2_ano';    var_all{8,2}='uatm';             var_all{8,3}=['air pco2 anomaly'];
      var_all{9,1}='apco2_m';      var_all{9,2}='uatm';             var_all{9,3}=['air pco2 mean'];
   end

end

