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_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_alk=dir([dir_subn,'/talk*.nc']);
    fn_dic=dir([dir_subn,'/dissic*.nc']);
    %if length(fn_dpco2)<=2 | length(fn_alk)<=2 | length(fn_dic)<=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_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_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);
        fgco2=ncread(fn_allin,'fgco2'); %gC/m2/year

        var_exist=check_var_nc(fn_allin,'dpco2');
        if var_exist; dpco2=ncread(fn_allin,'dpco2');
        else; dpco2=nan(size(fgco2)); end

        spco2=ncread(fn_allin,'spco2');
        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, kw_m]=compute_anomaly_bandpass(kw,time,12*15,13);
        fgco2_ano=ncread(fn_anoin,'fgco2'); %gC/m2/year

        var_exist=check_var_nc(fn_anoin,'dpco2');
        if var_exist; dpco2_ano=ncread(fn_anoin,'dpco2');
        else; dpco2_ano=nan(size(fgco2_ano)); end

        spco2_ano=ncread(fn_anoin,'spco2');
        tos         =ncread(fn_allin,'tos');
        sos         =ncread(fn_allin,'sos');
        dissicos    =ncread(fn_allin,'dissicos');
        talkos      =ncread(fn_allin,'talkos');

        [jm,im,lm]=size(dpco2);

        %fgco2_m=ncread(fn_meanin,'fgco2'); %gC/m2/year
        spco2_m=ncread(fn_meanin,'spco2');

        var_exist=check_var_nc(fn_meanin,'dpco2');
        if var_exist; dpco2_m=ncread(fn_meanin,'dpco2');
        else; dpco2_m=nan(size(spco2_m)); end

        %tos_m     =ncread(fn_meanin,'tos');
        sos_m     =ncread(fn_meanin,'sos');
        dissicos_m=ncread(fn_meanin,'dissicos');
        talkos_m  =ncread(fn_meanin,'talkos');
       
        %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; 

        fgco2_kwa_dpco2a=dpco2_ano.*kw_ano;
        [fgco2_kwa_dpco2a_ano, fgco2_kwa_dpco2a_m]=compute_anomaly_bandpass(fgco2_kwa_dpco2a,time,12*15,13);
        fgco2_ct=fgco2_kwa_dpco2a-fgco2_kwa_dpco2a_m;
        clear fgco2_kwa_dpco2a fgco2_kwa_dpco2a_ano fgco2_kwa_dpco2a_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');
        sos_ano     =ncread(fn_anoin,'sos');
        dissicos_ano=ncread(fn_anoin,'dissicos');
        talkos_ano  =ncread(fn_anoin,'talkos');

        %dissicos_m=movmean(dissicos,30*12,3);%spco2-spco2_ano;
        %talkos_m  =movmean(talkos,30*12,3);%spco2-spco2_ano;
        %sos_m     =movmean(sos,30*12,3);%spco2-spco2_ano;
        %spco2_m   =movmean(spco2,30*12,3);%spco2-spco2_ano;

        alpha_dic=(3*(talkos_m.*dissicos_m)-2*(dissicos_m.*dissicos_m))./((2*dissicos_m-talkos_m).*(talkos_m-dissicos_m));
        alpha_alk=(talkos_m.*talkos_m)./((2*dissicos_m-talkos_m).*(talkos_m-dissicos_m));
        dpddic=(spco2_m./dissicos_m).*(3*(talkos_m.*dissicos_m)-2*(dissicos_m.*dissicos_m))./((2*dissicos_m-talkos_m).*(talkos_m-dissicos_m));
        dpdalk=-1*(spco2_m./talkos_m).*(talkos_m.*talkos_m)./((2*dissicos_m-talkos_m).*(talkos_m-dissicos_m));
        dpdtemp=spco2_m*0.0423;
        dpdsalt=spco2_m./sos_m;

        spco2_temp=dpdtemp.*tos_ano;     
        spco2_salt=dpdsalt.*sos_ano;     
        spco2_dic=dpddic.*dissicos_ano; 
        spco2_alk=dpdalk.*talkos_ano;   
        spco2_resi=spco2_ano-spco2_temp-spco2_salt-spco2_dic-spco2_alk;
        spco2_nont=spco2_ano-spco2_temp;

        fgco2_temp=kw_m.*spco2_temp;
        fgco2_salt=kw_m.*spco2_salt;
        fgco2_dic=kw_m.*spco2_dic;
        fgco2_alk=kw_m.*spco2_alk;
        fgco2_nont=kw_m.*spco2_nont;


        %clear spco2_m dpdtemp dpddic dpdalk dpdsalt dissicos_m talkos_m sos_m

        [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),[],tos);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(10,:),lon,lat,time-datenum(1900,1,1),[],tos_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(11,:),lon,lat,time-datenum(1900,1,1),[],sos);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(12,:),lon,lat,time-datenum(1900,1,1),[],sos_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(13,:),lon,lat,time-datenum(1900,1,1),[],dissicos);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(14,:),lon,lat,time-datenum(1900,1,1),[],dissicos_ano);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(15,:),lon,lat,time-datenum(1900,1,1),[],talkos);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(16,:),lon,lat,time-datenum(1900,1,1),[],talkos_ano);

        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(17,:),lon,lat,time-datenum(1900,1,1),[],fgco2_kw);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(18,:),lon,lat,time-datenum(1900,1,1),[],fgco2_dpco2);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(19,:),lon,lat,time-datenum(1900,1,1),[],fgco2_ct);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(20,:),lon,lat,time-datenum(1900,1,1),[],fgco2_resi);

        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(21,:),lon,lat,time-datenum(1900,1,1),[],spco2_temp);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(22,:),lon,lat,time-datenum(1900,1,1),[],spco2_salt);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(23,:),lon,lat,time-datenum(1900,1,1),[],spco2_dic);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(24,:),lon,lat,time-datenum(1900,1,1),[],spco2_alk);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(25,:),lon,lat,time-datenum(1900,1,1),[],spco2_resi);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(26,:),lon,lat,time-datenum(1900,1,1),[],alpha_dic);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(27,:),lon,lat,time-datenum(1900,1,1),[],alpha_alk);

        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(28,:),lon,lat,time-datenum(1900,1,1),[],kw_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(29,:),lon,lat,time-datenum(1900,1,1),[],spco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(30,:),lon,lat,time-datenum(1900,1,1),[],dpco2_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(31,:),lon,lat,time-datenum(1900,1,1),[],dissicos_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(32,:),lon,lat,time-datenum(1900,1,1),[],talkos_m);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(33,:),lon,lat,time-datenum(1900,1,1),[],sos_m);

        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(34,:),lon,lat,time-datenum(1900,1,1),[],fgco2_temp);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(35,:),lon,lat,time-datenum(1900,1,1),[],fgco2_salt);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(36,:),lon,lat,time-datenum(1900,1,1),[],fgco2_dic);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(37,:),lon,lat,time-datenum(1900,1,1),[],fgco2_alk);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(38,:),lon,lat,time-datenum(1900,1,1),[],fgco2_nont);
        create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(39,:),lon,lat,time-datenum(1900,1,1),[],spco2_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_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}='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}='tos';          var_all{9,2}='degree';           var_all{9,3}=['sea surface temperature'];
      var_all{10,1}='tos_ano';     var_all{10,2}='degree';          var_all{10,3}=['sea surface temperature anomaly against 30 yrs mean'];
      var_all{11,1}='sos';         var_all{11,2}='psu';             var_all{11,3}=['sea surface salinity'];
      var_all{12,1}='sos_ano';     var_all{12,2}='psu';             var_all{12,3}=['sea surface salinity anomaly against 30 yrs mean'];
      var_all{13,1}='dissicos';    var_all{13,2}='umol/kg';         var_all{13,3}=['sea surface dissolved inorganic carbon'];
      var_all{14,1}='dissicos_ano';var_all{14,2}='umol/kg';         var_all{14,3}=['sea surface dissolved inorganic carbon anomaly against 30 yrs mean'];
      var_all{15,1}='talkos';      var_all{15,2}='umol/kg';         var_all{15,3}=['sea surface alkalinity'];
      var_all{16,1}='talkos_ano';  var_all{16,2}='umol/kg';         var_all{16,3}=['sea surface alkalinity anomaly against 30 yrs mean'];

      var_all{17,1}='fgco2_kw';    var_all{17,2}='gC/m2/yr';        var_all{17,3}=['fgco2 anomaly by kw anomaly'];
      var_all{18,1}='fgco2_dpco2'; var_all{18,2}='gC/m2/yr';        var_all{18,3}=['fgco2 anomaly by dpco2 anomaly'];
      var_all{19,1}='fgco2_ct';    var_all{19,2}='gC/m2/yr';        var_all{19,3}=['fgco2 anomaly by cross term anomaly'];
      var_all{20,1}='fgco2_resi';  var_all{20,2}='gC/m2/yr';        var_all{20,3}=['fgco2 anomaly residual'];

      var_all{21,1}='spco2_temp';  var_all{21,2}='uatm';            var_all{21,3}=['spco2 anomaly by temperature'];
      var_all{22,1}='spco2_salt';  var_all{22,2}='uatm';            var_all{22,3}=['spco2 anomaly by salinity'];
      var_all{23,1}='spco2_dic';   var_all{23,2}='uatm';            var_all{23,3}=['spco2 anomaly by dic'];
      var_all{24,1}='spco2_alk';   var_all{24,2}='uatm';            var_all{24,3}=['spco2 anomaly by alk'];
      var_all{25,1}='spco2_resi';  var_all{25,2}='uatm';            var_all{25,3}=['spco2 anomaly by resi'];

      var_all{26,1}='alpha_dic';   var_all{26,2}='';                var_all{26,3}=['Buffer factor (dic coefficient)'];
      var_all{27,1}='alpha_alk';   var_all{27,2}='';                var_all{27,3}=['alk coefficient'];
      var_all{28,1}='kw_m';        var_all{28,2}='';                var_all{28,3}=['kw_m 30 yrs moving mean'];
      var_all{29,1}='spco2_m';     var_all{29,2}='';                var_all{29,3}=['spco2_m 30 yrs moving mean'];
      var_all{30,1}='dpco2_m';     var_all{30,2}='';                var_all{30,3}=['dpco2_m 30 yrs moving mean'];
      var_all{31,1}='dissicos_m';  var_all{31,2}='';                var_all{31,3}=['DIC_m 30 yrs moving mean'];
      var_all{32,1}='talkos_m';    var_all{32,2}='';                var_all{32,3}=['Alk_m 30 yrs moving mean'];
      var_all{33,1}='sos_m';       var_all{33,2}='';                var_all{33,3}=['SSS 30 yrs moving mean'];

      var_all{34,1}='fgco2_temp';  var_all{34,2}='gC/m2/yr';        var_all{34,3}=['fgco2 anomaly by temperature'];
      var_all{35,1}='fgco2_salt';  var_all{35,2}='gC/m2/yr';        var_all{35,3}=['fgco2 anomaly by salinity'];
      var_all{36,1}='fgco2_dic';   var_all{36,2}='gC/m2/yr';        var_all{36,3}=['fgco2 anomaly by dic'];
      var_all{37,1}='fgco2_alk';   var_all{37,2}='gC/m2/yr';        var_all{37,3}=['fgco2 anomaly by alk'];
      var_all{38,1}='fgco2_nont';  var_all{38,2}='gC/m2/yr';        var_all{38,3}=['fgco2 anomaly by nont'];
      var_all{39,1}='spco2_nont';  var_all{39,2}='uatm';            var_all{39,3}=['spco2 anomaly by nont'];
   end

end

