clear; close all; tic; 

variable_ids={'tos','fgco2','dpco2','spco2','talkos','dissicos','sos'};
experiment_ids={'historical','ssp245','ssp585'};

dir_home='/tigress/enhuil/analysis/co2_enso_projection/ana_fgco2_by_member2/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=9%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']);
    if length(fn_dpco2)==0; 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_out=[dir_subn,'/all_decomposition_fgco2_spco2_',dir_suball(i).name,'.his_ssp245_1901_2100_tropical_pacific.nc'];
    fn5_out=[dir_subn,'/all_decomposition_fgco2_spco2_',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_out=fn2_out; end 
        if j==2; fn_allin=fn5_allin; fn_anoin=fn5_anoin; fn_out=fn5_out; 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);
        fgco2=ncread(fn_allin,'fgco2'); %gC/m2/year
        dpco2=ncread(fn_allin,'dpco2');
        spco2=ncread(fn_allin,'spco2');
        tos=ncread(fn_allin,'tos');
        kw=fgco2./dpco2; %gC*m/m3/uatm/yr
        %at some point in some time, the dpco2 is very small (~0.001), this could lead to a huge kw
        kw(kw>1)=1;
        kw(kw<0.01)=0.01;

        kw_ano=compute_anomaly_bandpass(kw,time,12*15,13);
        fgco2_ano=ncread(fn_anoin,'fgco2'); %gC/m2/year
        dpco2_ano=ncread(fn_anoin,'dpco2');
        spco2_ano=ncread(fn_anoin,'spco2');
        
        [jm,im,lm]=size(dpco2);
        dpco2_m=nan(jm,im,lm);
        kw_m=nan(jm,im,lm);
        spco2_m=nan(jm,im,lm);
        for j=1:jm; for i=1:im
            d0=squeeze(dpco2(j,i,:));
            if length(d0(~isnan(d0)))==lm
               dpco2_m(j,i,:)=smooth(squeeze(dpco2(j,i,:)),30*12);%dpco2-dpco2_ano;
               kw_m(j,i,:)=smooth(squeeze(kw(j,i,:)),30*12);%dpco2-dpco2_ano;
               spco2_m(j,i,:)=smooth(squeeze(spco2(j,i,:)),30*12);%dpco2-dpco2_ano;
            end
        end;end
        dpco2_m(:,:,end-23:end)=nan;
        kw_m(:,:,end-23:end)=nan;
        spco2_m(:,:,end-23:end)=nan;

        %fgco2 decomposition 
        %f'=k'p_m+k_m*p'+k'p'-(k'p')_m
        fgco2_kw=dpco2_m.*kw_ano;
        fgco2_dpco2=dpco2_ano.*kw_m; 
        a=dpco2_ano.*kw_ano;
        [im,jm,lm]=size(a);
        b=repmat(nanmean(a,3),[1,1,lm]);
        fgco2_ct=a-b; clear a b;
        fgco2_resi=fgco2_ano-fgco2_kw-fgco2_dpco2-fgco2_ct;
        
        %spco2 decomposition 
        %spco2_temp=0.0423*spco2_m*delta_temp
        tos_ano=ncread(fn_anoin,'tos');
        spco2_temp=0.0423*(spco2_m.*tos_ano);
        spco2_nont=spco2_ano-spco2_temp;
        fgco2_temp=kw_m.*spco2_temp;
        fgco2_nont=kw_m.*spco2_nont;

        [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),[],fgco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(2,:),lon,lat,time-datenum(1900,1,1),[],fgco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(3,:),lon,lat,time-datenum(1900,1,1),[],dpco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(4,:),lon,lat,time-datenum(1900,1,1),[],dpco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(5,:),lon,lat,time-datenum(1900,1,1),[],spco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(6,:),lon,lat,time-datenum(1900,1,1),[],spco2_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(7,:),lon,lat,time-datenum(1900,1,1),[],kw);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(8,:),lon,lat,time-datenum(1900,1,1),[],kw_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(9,:),lon,lat,time-datenum(1900,1,1),[],fgco2_kw);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(10,:),lon,lat,time-datenum(1900,1,1),[],fgco2_dpco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(11,:),lon,lat,time-datenum(1900,1,1),[],fgco2_ct);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(12,:),lon,lat,time-datenum(1900,1,1),[],fgco2_resi);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(13,:),lon,lat,time-datenum(1900,1,1),[],spco2_temp);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(14,:),lon,lat,time-datenum(1900,1,1),[],spco2_nont);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(15,:),lon,lat,time-datenum(1900,1,1),[],dpco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(16,:),lon,lat,time-datenum(1900,1,1),[],spco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(17,:),lon,lat,time-datenum(1900,1,1),[],kw_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(18,:),lon,lat,time-datenum(1900,1,1),[],tos);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(19,:),lon,lat,time-datenum(1900,1,1),[],tos_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(20,:),lon,lat,time-datenum(1900,1,1),[],fgco2_temp);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(21,:),lon,lat,time-datenum(1900,1,1),[],fgco2_nont);
        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_nont; spco2_temp=0.0423*spco2_m.*tos_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}='fgco2';        var_all{1,2}='gC/m2/yr';         var_all{1,3}=['CO2 flux (positive to air)'];
      var_all{2,1}='fgco2_ano';    var_all{2,2}='gC/m2/yr';         var_all{2,3}=['CO2 flux (positive to air) anomaly against 30 yrs mean'];
      var_all{3,1}='dpco2';        var_all{3,2}='uatm';             var_all{3,3}=['Delta CO2 (ocean-atm) Partial Pressure'];
      var_all{4,1}='dpco2_ano';    var_all{4,2}='uatm';             var_all{4,3}=['Delta CO2 (ocean-atm) Partial Pressure anomaly against 30 yrs mean'];
      var_all{5,1}='spco2';        var_all{5,2}='uatm';             var_all{5,3}=['sea surface CO2 Partial Pressure'];
      var_all{6,1}='spco2_ano';    var_all{6,2}='uatm';             var_all{6,3}=['sea surface CO2 Partial Pressure anomaly against 30 yrs mean'];
      var_all{7,1}='kw';           var_all{7,2}='gC*m/m3/uatm/yr';  var_all{7,3}=['kw (fgco2/dpco2)'];
      var_all{8,1}='kw_ano';       var_all{8,2}='gC*m/m3/uatm/yr';  var_all{8,3}=['kw (fgco2/dpco2) anomaly against 30 yrs mean'];

      var_all{9,1}='fgco2_kw';     var_all{9,2}='gC/m2/yr';         var_all{9,3}=['fgco2 anomaly by kw anomaly'];
      var_all{10,1}='fgco2_dpco2'; var_all{10,2}='gC/m2/yr';        var_all{10,3}=['fgco2 anomaly by dpco2 anomaly'];
      var_all{11,1}='fgco2_ct';    var_all{11,2}='gC/m2/yr';        var_all{11,3}=['fgco2 anomaly by cross term anomaly'];
      var_all{12,1}='fgco2_resi';  var_all{12,2}='gC/m2/yr';        var_all{12,3}=['fgco2 anomaly residual'];

      var_all{13,1}='spco2_temp';  var_all{13,2}='uatm';            var_all{13,3}=['spco2 anomaly by temperature'];
      var_all{14,1}='spco2_nont';  var_all{14,2}='uatm';            var_all{14,3}=['spco2 anomaly by non-temperature'];

      var_all{15,1}='dpco2m';      var_all{15,2}='uatm';            var_all{15,3}=['Delta CO2 (ocean-atm) Partial Pressure 30 yrs mean'];
      var_all{16,1}='spco2m';      var_all{16,2}='uatm';            var_all{16,3}=['sea surface CO2 Partial Pressure 30 yrs mean'];
      var_all{17,1}='kwm';         var_all{17,2}='gC*m/m3/uatm/yr'; var_all{17,3}=['kw (fgco2/dpco2) 30 yrs mean'];

      var_all{18,1}='tos';         var_all{18,2}='degree';          var_all{18,3}=['sea surface temperature'];
      var_all{19,1}='tos_ano';     var_all{19,2}='degree';          var_all{19,3}=['sea surface temperature anomaly against 30 yrs mean'];

      var_all{20,1}='fgco2_temp';  var_all{20,2}='gC/m2/yr';        var_all{20,3}=['fgco2 anomaly by dpco2-temp anomaly'];
      var_all{21,1}='fgco2_nont';  var_all{21,2}='gC/m2/yr';        var_all{21,3}=['fgco2 anomaly by dpco2-nont anomaly'];
   end

end

