
area=ncread('areacella_fx_NorESM1-M_piControl_r0i0p0.nc','areacella');
land=ncread('sftlf_fx_NorESM1-M_piControl_r0i0p0.nc','sftlf')/100;
FNAME='/cluster/work/users/ingo/noresm/NF2005AERAMIP_f19_f19_CMIP5aero_01/run/NF2005AERAMIP_f19_f19_CMIP5aero_01.cam2.h1.1979-01-02-00000.nc'

% fluxes read from file 
VARS={'BC_AX','BC_N','BC_NI','OM_NI','SO2','SO4_PR'}; 
for n=1:numel(VARS)
  VNAME=char(VARS{n});
  data=(ncread(FNAME,[VNAME '_CLXF'])+ncread(FNAME,[VNAME '_EMISRF']))*24*3600*365*1e-9;  % Tg/s
  gtot(n)=sum(sum(data.*area));
  disp([VNAME ' = ' num2str(gtot(n))])
end 
disp([])

% total emissions after adding secondary contributions (old implementation)
%VARS={'BC','BC_AIR','BC_N','BC_NI','BC_AX','POM','SO2','SO4'}
VARS={'BC','BC_N','BC_NI','BC_AX','POM','SO2','SO4'}
for n=1:numel(VARS)
  VNAME=char(VARS{n});
  data=ncread(FNAME,['EMI_' VNAME])*24*3600*365*1e-9;  % Tg/s
  gtotold(n)=sum(sum(data.*area));
  disp([VNAME ' = ' num2str(gtotold(n))])
end 
disp([])

% total emissions after adding secondary contributions (new implementation)
%VARS={'BC','POM','SO2','SO4'}
VARS={'BC','BC_N','BC_NI','BC_AX','POM','SO2','SO4'}
for n=1:numel(VARS)
  VNAME=char(VARS{n});
  data=ncread(FNAME,['EMI_' VNAME '_CMIP6'])*24*3600*365*1e-9;  % Tg/s
  gtotnew(n)=sum(sum(data.*area));
  disp([VNAME ' = ' num2str(gtotnew(n)) ' (' int2str(gtotnew(n)/gtotold(n)*100) '%)'])
end
disp([])

