EARTH=6371000;
lat=ncread('volcSynth_mem1-60_2006-2099.nc','lat');
for n=1:length(lat)
 if n==1
   lats=-90;
   latn=0.5*(lat(n)+lat(n+1));
 elseif n==length(lat)
   lats=0.5*(lat(n)+lat(n-1));
   latn=90;
 else
   lats=0.5*(lat(n)+lat(n-1));
   latn=0.5*(lat(n)+lat(n+1));
 end
 area(n)=2*pi*EARTH^2*(sin(latn/180*pi)-sin(lats/180*pi));
end
tdm=1142;
colmass=zeros(tdm,NMEM);
data=ncread('volcSynth_mem1-60_2006-2099.nc','colmass');
for m=1:size(data,2)
  colmass(:,m)=area*squeeze(data(:,m,:))*1e-9; % kg -> Tg ; 
end

