
area=ncread('areacella_fx_NorESM1-M_piControl_r0i0p0.nc','areacella');
land=ncread('sftlf_fx_NorESM1-M_piControl_r0i0p0.nc','sftlf')/100;

% CMIP5 
PREFIX='/tos-project1/NS2345K/noresm/inputdata/atm/cam/camoslo/emissions/IPCC_emiss_';
SUFFIX='_1850-2005.nc';
VARS_CMIP5={'bcair','bcbb','bcff','pombb','pomff','so2bb','so2ff'};
figure(1)
clf
for n=1:numel(VARS_CMIP5)
  VNAME=char(VARS_CMIP5{n});
  switch VNAME 
  case {'bcair','bcff'} 
    VNAME2='emiss_bc';
    FAC=1;
  case 'bcbb' 
    VNAME2='emiss_bc_bb';
    FAC=1;
  case 'pombb'
    VNAME2='emiss_pom_bb';
    FAC=2.6;
  case 'pomff' 
    VNAME2='emiss_pom';
    FAC=1.4;
  case 'so2bb'
    VNAME2='emiss_so2_bb';
    FAC=0.5;
  case 'so2ff' 
    VNAME2='emiss_so2';
    FAC=0.5;
  end 
  data=squeeze(sum(ncread([PREFIX VNAME SUFFIX],VNAME2,[1 1 1 1 17],[inf inf inf 1 1]),3))*24*3600*365*1e-9;  % Tg/s
  gtot2(n)=sum(sum(data.*area));
  ind=find(land>0.01); 
  glnd2(n)=sum(sum(data(ind).*area(ind)));
  % 
  subplot(4,4,n) 
  bar([gtot2(n) glnd2(n) gtot2(n)-glnd2(n)]); 
  xlim=get(gca,'xlim');ylim=get(gca,'ylim');
  text(xlim(2),ylim(2),VNAME,'horiz','right','vert','top','interpreter','none'); 
  disp([VNAME ' = ' num2str(gtot2(n))])
  disp([VNAME '_scaled = ' num2str(gtot2(n)*FAC)])
end 

% CMIP5 new 
disp([])
disp('--- CMIP5 NorESM1.2 ---')
PREFIX='/cluster/shared/noresm/inputdata/noresm-only/atm/cam/chem/trop_mozart_aero/emis/ipcc_ar5/ar5_mam3_';
SUFFIX='_1850-2005_c090804.nc';
VARS_CMIP5={'BCFFX','BCFFN','bc','so2','so4_pr','oc','POM'};
kgs_to_Tgyr=24*3600*365*1e-9;
for n=1:numel(VARS_CMIP5)
  VNAME=char(VARS_CMIP5{n});
  switch VNAME
  case {'BCFFX','BCFFN'}
    VARS={'emiss_dom','emiss_ene','emiss_ind','emiss_shp','emiss_tra','emiss_wst','emiss_awb'};
    VARV=[];
    FAC=1;
    MASS=12.011;
  case {'bc'}
    VARV={'forestfire','grassfire'};
    VARS=[];
    FAC=1;
    MASS=12.011;
  case {'so2'}
    VARS={'emiss_dom','emiss_awb','emiss_shp','emiss_tra','emiss_wst'};
    VARV={'emiss_ene','emiss_ind','contvolc','forestfire','grassfire'};
    %VARV={'contvolc'};
    FAC=1;
    MASS=64.0648;
  case {'so4_pr'}
    VARV={'emiss_ene','emiss_ind','contvolc','forestfire','grassfire'};
    VARS={'emiss_dom','emiss_awb','emiss_shp','emiss_tra','emiss_wst'};
    FAC=1;
    MASS=98.0784;
  case {'oc'} 
    VARV={'forestfire','grassfire'};
    VARS=[];
    FAC=2.6;
    MASS=12.011;
  case {'POM'}
    VARS={'emiss_dom','emiss_ene','emiss_ind','emiss_shp','emiss_tra','emiss_wst'};
    VARV=[];
    FAC=1.4;
    MASS=12.011;
    SUFFIX='_1850-2005_c100429.nc';
  end
  data_surf=zeros(size(area));
  data_vert=zeros(size(area));
  if ~isempty(VARV) 
    altitude_int=ncread([PREFIX VNAME '_elev' SUFFIX],'altitude_int');
  end
  moleccm2_to_kgm2=1.0e4/6.02e23*MASS*1.e-3;
  for m=1:numel(VARS)
    disp(VARS{m})
    %data_surf=data_surf+FAC*mean(ncread([PREFIX VNAME '_surf' SUFFIX],VARS{m},[1 1 193],[inf inf 12]),3)*...
    data_surf=data_surf+FAC*mean(ncread([PREFIX VNAME '_surf' SUFFIX],VARS{m},[1 1 193],[inf inf 1]),3)*...
              moleccm2_to_kgm2*kgs_to_Tgyr;
  end
  for m=1:numel(VARV)
    for k=1:numel(altitude_int)-1
      %data_vert=data_vert+FACS*squeeze(mean(ncread([PREFIX VNAME '_elev' SUFFIX],VARV{m},[1 1 k 193],[inf inf 1 12]),4))...
      data_vert=data_vert+FAC*squeeze(mean(ncread([PREFIX VNAME '_elev' SUFFIX],VARV{m},[1 1 k 193],[inf inf 1 1]),4))...
              *(altitude_int(k+1)-altitude_int(k))*1e5*moleccm2_to_kgm2*kgs_to_Tgyr;
    end
  end
  tot_surf(n)=sum(data_surf(:).*area(:));
  tot_vert(n)=sum(data_vert(:).*area(:));
  disp([])
  disp([VNAME ' surf = ' num2str(tot_surf(n))])
  disp([VNAME ' vert = ' num2str(tot_vert(n))])
  disp([VNAME ' both = ' num2str(tot_surf(n)+tot_vert(n))])
end

% CMIP6 
disp([])
disp('--- CMIP6 ---')
PREFIX='/cluster/shared/noresm/inputdata/atm/cam/chem/emis/CMIP6_emissions_1750_2015_240717_CAM5.3_Oslo_version20171122/emissions-cmip6_';
SUFFIX='_1849-2015_1.9x2.5_version20171122.nc';
VARS_CMIP6={'BC_AX','BC_N','BC_NI','OM_NI','SO2','SO4_PR'}; 
kgs_to_Tgyr=24*3600*365*1e-9;
for n=1:numel(VARS_CMIP6)
  VNAME=char(VARS_CMIP6{n});
  switch VNAME
  case {'BC_AX','BC_N'}
    VARS={'emiss_anthro'};
    VARV={'emiss_aircraft'};
    FAC=1;
    MASS=12.011;
  case {'BC_NI'}
    VARV={'emiss_bb'};
    VARS=[];
    FAC=1;
    MASS=12.011;
  case {'SO2'}
    VARV={'emiss_aircraft','emiss_volcanoes','emiss_bb','emiss_ene_ind'};
    %VARV={'emiss_volcanoes'};
    VARS={'emiss_ag_sol_was','emiss_ship','emiss_res_tran'};
    FAC=1;
    MASS=64.0648;
  case {'SO4_PR'}
    VARV={'emiss_volcanoes','emiss_bb','emiss_ene_ind'};
    %VARV={'emiss_volcanoes'};
    VARS={'emiss_ag_sol_was','emiss_shipping','emiss_res_tran'};
    FAC=1;
    MASS=98.0784;
  case {'OM_NI'}
    VARV={'emiss_bb'};
    VARS={'emiss_anthro'};
    FAC=1;
    MASS=12.011; 
  end 
  data_surf=zeros(size(area));
  data_vert=zeros(size(area));
  if ~isempty(VARS)   
    altitude_int=ncread([PREFIX VNAME '_all_vertical' SUFFIX],'altitude_int');
  end 
  for m=1:numel(VARS)
    moleccm2_to_kgm2=1.0e4/6.02e23*MASS(min(m,numel(MASS)))*1.e-3;
    data_surf=data_surf+FAC*mean(ncread([PREFIX VNAME '_all_surface' SUFFIX],VARS{m},[1 1 (2004-119)],[inf inf 120]),3)*...
              moleccm2_to_kgm2*kgs_to_Tgyr;
  end
  for m=1:numel(VARV)
    moleccm2_to_kgm2=1.0e4/6.02e23*MASS(min(m,numel(MASS)))*1.e-3;
    for k=1:numel(altitude_int)-1 
      data_vert=data_vert+FAC*squeeze(mean(ncread([PREFIX VNAME '_all_vertical' SUFFIX],VARV{m},[1 1 k (2004-119)],[inf inf 1 120]),4))...
                *(altitude_int(k+1)-altitude_int(k))*1e5*moleccm2_to_kgm2*kgs_to_Tgyr;
    end 
  end
  tot_surf(n)=sum(data_surf(:).*area(:)); 
  tot_vert(n)=sum(data_vert(:).*area(:)); 
  disp([])
  disp([VNAME ' surf = ' num2str(tot_surf(n))])
  disp([VNAME ' vert = ' num2str(tot_vert(n))])
  disp([VNAME ' both = ' num2str(tot_surf(n)+tot_vert(n))])
end      



