clear
close all


fname_fco2_ano='/tigress/GEOCLIM/LRGROUP/shared_data/pco2_flux_all_mom6kw/model_spinup82yrs_ls_rb_co2_data_noarea_detrend_deseason.nc';
fname_fco2='/tigress/GEOCLIM/LRGROUP/shared_data/pco2_flux_all_mom6kw/model_spinup82yrs_ls_rb_co2_data.nc';
fname_sst_ano='/tigress/GEOCLIM/LRGROUP/shared_data/sst_OISSTmonthly/allsst.avhrr-only-v2.1981_2018_monthly_detrend_deseason.nc';
fname_sst='/tigress/GEOCLIM/LRGROUP/shared_data/sst_OISSTmonthly/allsst.avhrr-only-v2.1981_2018_monthly.nc';

fco2_region=[0 360 -90 90]; 
trange=[datenum(1982,1,1) datenum(2018,1,1)];
nobs_all={'ls','rb'}; 

for n=1:2
    nobs=char(nobs_all(n));
    xname=['xh_',nobs]; yname=['yh_',nobs];
    fn_out=['obs_',nobs,'.nc'];
    
    %landschutzer data
    tic;disp(['load fCO2 data ...'])
    [fco2,lon,lat,time]=read_1ncobsnd(fname_fco2,['fco2_',nobs],xname,yname,fco2_region,trange,12,3); %mol/m2/yr-->gC/m2/yr
    [dpco2,lon,lat,time]=read_1ncobsnd(fname_fco2,['deltap_',nobs],xname,yname,fco2_region,trange,1,3); %uatm
    [kw,lon,lat,time]=read_1ncobsnd(fname_fco2,['kw_sol_',nobs],xname,yname,fco2_region,trange,12,3); %mol*m/m3/uatm/yr --> gC*m/m3/uatm/yr
    [pco2,lon,lat,time]=read_1ncobsnd(fname_fco2,['pco2w_',nobs],xname,yname,fco2_region,trange,1,3); %mol*m/m3/uatm/yr --> gC*m/m3/uatm/yr
    
    [fco2_ano,lon,lat,time]=read_1ncobsnd(fname_fco2_ano,['fco2_',nobs],xname,yname,fco2_region,trange,12,3); %mol/m2/yr-->gC/m2/yr
    [dpco2_ano,lon,lat,time]=read_1ncobsnd(fname_fco2_ano,['deltap_',nobs],xname,yname,fco2_region,trange,1,3); %uatm
    [kw_ano,lon,lat,time]=read_1ncobsnd(fname_fco2_ano,['kw_sol_',nobs],xname,yname,fco2_region,trange,12,3); %mol*m/m3/uatm/yr --> gC*m/m3/uatm/yr
    [pco2_ano,lon,lat,time]=read_1ncobsnd(fname_fco2_ano,['pco2w_',nobs],xname,yname,fco2_region,trange,1,3); %mol*m/m3/uatm/yr --> gC*m/m3/uatm/yr
    
    toc;disp(['compute fCO2 decomposition ...']); tic
    kwstd=nanstd(kw(:)); kwmean=nanmeankw(:)); 
    stop
    kw(kw>(kwmean+kwstd))=; kw(kw<(kwmean-kwstd)0.01)=0.01; %kw in rodenbeck are calculated from fco2/dpco2, for small dpco2 point, the kw could be extremely large.
    %kw(kw>(kwmean+kwstd))=1; kw(kw<(kwmean-kwstd)0.01)=0.01; %kw in rodenbeck are calculated from fco2/dpco2, for small dpco2 point, the kw could be extremely large.
    [jm,im,lm]=size(fco2);
    dpco2_m=repmat(nanmean(dpco2,3),[1,1,lm]);
    kw_m=repmat(nanmean(kw,3),[1,1,lm]);
    fco2_kw=kw_ano.*dpco2_m;
    fco2_dpco2=kw_m.*dpco2_ano;
    a=dpco2_ano.*kw_ano;
    b=repmat(nanmean(a,3),[1,1,lm]);
    fco2_ct=a-b; 
    clear a b 
    fco2_resi=fco2_ano-fco2_kw-fco2_dpco2-fco2_ct;
    
    toc;disp(['load sst data ...']); tic
    %sst data
    [sst,lon_sst,lat_sst,time_sst]=read_1ncobsnd(fname_sst,'sst','lon','lat',fco2_region,trange,1,3); %uatm
    [sst_ano,lon_sst,lat_sst,time_sst]=read_1ncobsnd(fname_sst_ano,'sst','lon','lat',fco2_region,trange,1,3); %uatm
    [jm_sst,im_sst,lm_sst]=size(sst_ano);
    sst_c=nan(jm,im,lm); 
    sst_c_ano=nan(jm,im,lm); 
    toc;disp(['regrid sst to fCO2 point ...']); tic
    for l=1:lm_sst
        if mod(l,100)==0; toc; disp([num2str(l),'/',num2str(lm_sst)]); tic; end
        sst_c(:,:,l)=interp2(lon_sst,lat_sst,squeeze(sst(:,:,l)),lon,lat);
        sst_c_ano(:,:,l)=interp2(lon_sst,lat_sst,squeeze(sst_ano(:,:,l)),lon,lat);
    end
    
    toc;disp(['compute pCO2 decomposition ...']); tic
    pco2_m=repmat(nanmean(pco2,3),[1,1,lm]);
    pco2_temp=0.0423*(pco2_m.*sst_c_ano); 
    pco2_nont=pco2_ano-pco2_temp;

    fco2_temp=kw_m.*pco2_temp;
    fco2_nont=kw_m.*pco2_nont;
 
    toc;disp(['write nc...']); tic
    [var_all]=load_var_all(1); lon=lon(1,:); lat=lat(:,1);
    %write fco2, fco2_ano, dpco2, dpco2_ano, pco2, pco2_ano, kw, kw_ano, fco2_kw, fco2_dpco2, fco2_ct, fco2_resi
    create_mom6_1var3d(fn_out,'lon','lat','time','time_bnds',var_all(1,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(2,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_ano,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(3,:),lon,lat,time-datenum(1900,1,1),[],permute(dpco2,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(4,:),lon,lat,time-datenum(1900,1,1),[],permute(dpco2_ano,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(5,:),lon,lat,time-datenum(1900,1,1),[],permute(pco2,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(6,:),lon,lat,time-datenum(1900,1,1),[],permute(pco2_ano,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(7,:),lon,lat,time-datenum(1900,1,1),[],permute(kw,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(8,:),lon,lat,time-datenum(1900,1,1),[],permute(kw_ano,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(9,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_kw,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(10,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_dpco2,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(11,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_ct,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(12,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_resi,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(13,:),lon,lat,time-datenum(1900,1,1),[],permute(pco2_temp,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(14,:),lon,lat,time-datenum(1900,1,1),[],permute(pco2_nont,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(15,:),lon,lat,time-datenum(1900,1,1),[],permute(sst_c,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(16,:),lon,lat,time-datenum(1900,1,1),[],permute(sst_c_ano,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(17,:),lon,lat,time-datenum(1900,1,1),[],permute(kw_m,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(18,:),lon,lat,time-datenum(1900,1,1),[],permute(dpco2_m,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(19,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_temp,[2,1,3]));
    create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(20,:),lon,lat,time-datenum(1900,1,1),[],permute(fco2_nont,[2,1,3]));
    ncwriteatt(fn_out,'/','data description',['this is ',nobs,' observational result as a reference for CMIP6 analysis']);
    ncwriteatt(fn_out,'/','budget1',['kw=fco2/dpco2']);
    ncwriteatt(fn_out,'/','budget2',['fco2_ano=fco2_kw+fco2_dpco2+fco2_ct+fco2_resi']);
    ncwriteatt(fn_out,'/','budget3',['pco2_ano=pco2_temp+pco2_nont; pco2_temp=0.0423*pco2_m.*tos_ano']);
end


function [var_all]=load_var_all(load_ind)

   if load_ind==1;
      %write fco2, fco2_ano, dpco2, dpco2_ano, pco2, pco2_ano, kw, kw_ano, fco2_kw, fco2_dpco2, fco2_ct, fco2_resi
      var_all{1,1}='fco2';        var_all{1,2}='gC/m2/yr';         var_all{1,3}=['CO2 flux (positive to air)'];
      var_all{2,1}='fco2_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}='pco2';        var_all{5,2}='uatm';             var_all{5,3}=['sea surface CO2 Partial Pressure'];
      var_all{6,1}='pco2_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 (fco2/dpco2)'];
      var_all{8,1}='kw_ano';       var_all{8,2}='gC*m/m3/uatm/yr'; var_all{8,3}=['kw (fco2/dpco2) anomaly against 30 yrs mean'];
      var_all{9,1}='fco2_kw';     var_all{9,2}='gC/m2/yr';         var_all{9,3}=['fco2 anomaly by kw anomaly'];
      var_all{10,1}='fco2_dpco2'; var_all{10,2}='gC/m2/yr';        var_all{10,3}=['fco2 anomaly by dpco2 anomaly'];
      var_all{11,1}='fco2_ct';    var_all{11,2}='gC/m2/yr';        var_all{11,3}=['fco2 anomaly by cross term anomaly'];
      var_all{12,1}='fco2_resi';  var_all{12,2}='gC/m2/yr';        var_all{12,3}=['fco2 anomaly residual'];
      var_all{13,1}='pco2_temp';  var_all{13,2}='uatm';            var_all{13,3}=['pco2 anomaly by temperature'];
      var_all{14,1}='pco2_nont';  var_all{14,2}='uatm';            var_all{14,3}=['pco2 anomaly by non-temperature'];
      var_all{15,1}='sst';        var_all{15,2}='degree';          var_all{15,3}=['sea surface temperature from OISST'];
      var_all{16,1}='sst_ano';    var_all{16,2}='degree';          var_all{16,3}=['sea surface temperature anomaly from OISST'];
      var_all{17,1}='kwm';        var_all{17,2}='gC*m/m3/uatm/yr'; var_all{17,3}=['kw (fco2/dpco2) mean'];
      var_all{18,1}='dpco2m';     var_all{18,2}='uatm';            var_all{18,3}=['Delta CO2 (ocean-atm) Partial Pressure mean'];
      var_all{19,1}='fco2_temp';  var_all{19,2}='gC/m2/yr';        var_all{19,3}=['fco2 anomaly by pco2 temp anomaly'];
      var_all{20,1}='fco2_nont';  var_all{20,2}='gC/m2/yr';        var_all{20,3}=['fco2 anomaly by pco2 nont anomaly'];
   end

end


function [data,lon,lat,time]=read_1ncobsnd(fname_all,vname,xname,yname,nino_region,trange,vs,nd);
   fname=fname_all;vname=char(vname);
   fstat=dir(fname);
   lonmin=nino_region(1); lonmax=nino_region(2);
   latmin=nino_region(3); latmax=nino_region(4);
   lon0=double(ncread(fname,xname));
   [indi]=find(lon0>=lonmin & lon0<=lonmax);
   ib=min(indi);ie=max(indi);in=ie-ib+1;

   lat0=double(ncread(fname,yname));
   [indj]=find(lat0>=latmin & lat0<=latmax);
   jb=min(indj);je=max(indj);jn=je-jb+1;
   [lon,lat]=meshgrid(lon0(ib:ie),lat0(jb:je));

   tb=datenum(1900,1,1);
   if contains(vname,'sst'); tb=datenum(1981,9,1); end
   time0=double(ncread(fname,'time'))+tb;
   [indl]=find(time0>=trange(1) & time0<=trange(2));
   lb=min(indl);le=max(indl);ln=le-lb+1;
   time=time0(lb:le);
   
   data=ncread(fname,vname,[ib,jb,lb],[in,jn,ln]);
   data=permute(data,[2,1,3])*vs;
   [im,jm,mm,lm]=size(data);
   if nd==1; %get 1d time series
      if contains(vname,'fco2'); %if vname include string "fco2"
         area=grid_area(lon,lat);
         area3d=repmat(area,[1,1,lm]);
         data0=data;
         data=squeeze(nansum(nansum(data.*area3d,1),2))/1e15;

         %remove flux in the atlantic gulf of mexico
         if lonmin==160 & lonmax==285 & latmin==-15 & latmax==15
            lonmin2=276; lonmax2=285; latmin2=9; latmax2=15;
            lon0=lon(1,:);lat0=lat(:,1);
            [indi]=find(lon0>=lonmin2 & lon0<=lonmax2);
            ib=min(indi);ie=max(indi);
            [indj]=find(lat0>=latmin2 & lat0<=latmax2);
            jb=min(indj);je=max(indj);

            data_rm=data0(jb:je,ib:ie,:);
            area3d_rm=area3d(jb:je,ib:ie,:);
            data_rm_sum=squeeze(nansum(nansum(data_rm.*area3d_rm,1),2))/1e15;
            data=data-data_rm_sum;
         end
         clear area3d data0
      else
         data=squeeze(nanmean(nanmean(data,1),2));
      end
   elseif nd==3
      data=data;
   else
      error('nd>1 does not work')
   end
   %toc;
end


