function [water_kg,bed_kg,jdmat]=roms_sed_mass(fname,NCS);
% ROMS_SED_MASS calculates total amount of sediment 
%  Usage: [water_kg,bed_kg,jdmat]=roms_sed_mass(fname,NCS);
%    fname= roms history file 
%    NCS= number of sediment classes  (e.g. 3= MUD_01, MUD_02, MUD_03)
%  Example: roms_sed_mass('ocean_his.nc',3); 
%          
jdmat=roms_time(fname,'ocean_time');
nt=length(jdmat);

grd=roms_get_grid(fname,fname,1,1);
[Mp,Lp]=size(grd.pm);
if isempty(grd.mask_rho),
  mask=ones(size(grd.pm));
else
  mask=grd.mask_rho;
end
mask(1,:)=0;
mask(:,end)=0;
mask(end,:)=0;
mask(:,1)=0;
dx=1./grd.pm;
dy=1./grd.pn;
area2d=dx.*dy;

nc=netcdf(fname);
water_kg=ones(nt,NCS)*NaN;
for i=1:nt,
% Compute 3d field of dz at each time step
 z=roms_depths(fname,fname,5,0,i);
 dz=diff(z,1,3);
 dz=permute(dz,[3 2 1]);
%   For each sediment class, compute total kilograms of sediment 
 for k=1:NCS;
   var1=sprintf('mud_%2.2d',k);       % Water column mass
   c1=nc{var1}(i,:,:,:);
   vertsum1=squeeze(sum(c1.*dz,1));
   water_kg(i,k)=sum(vertsum1(:).*area2d(:).*mask(:));

   var2=sprintf('mudmass_%2.2d',k);   % Bed mass
   c2=nc{var2}(i,:,:,:);
   vertsum2=squeeze(sum(c2,1));
   bed_kg(i,k)=sum(vertsum2(:).*area2d(:).*mask(:));
 end
end
close(nc);