function bed_mass(name,F,diry,rect);%,start,stop);
% function [bm_2d time] = bed_mass(name,F,diry);
%
% This function is used to loop through ROMS netcdf files and pull out the
% mass of sediment on the bed. It currently needs to be hardcoded in to
% change the river/bed sediment variables to match the model formulation.
%
%   name = string to search for files ('*his*')
%   F = flag for using only fluvial (1) or bed (2) sediment.
%   diry = directory to search in
%   rect = the rectangle of points to get, when not getting the whole model
%      domain. rect = [sb nb wb eb], where nb = north boundary, etc. At this
%      point nb/sb are the first dimension, and eb/wb are the second
%      dimension (in matlab).
% IMPORTANT - rect needs to be indexed based on starting at zero!!!!!!!!!
%
%   start = date to start with (not included)
%   stop = date to end with (not included) this feature was in an earlier
%   version of the file but I never used it.
%
%           Aaron Bever May 2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin ==3
    Fname=dir([diry,name]);
else
    diry = pwd;
    diry = [diry,'/'];
    Fname=dir([diry,name]);
end

if(~exist('rect','var'))
     S = get_varsize(Fname(1).name,'mudmass_04');
     rect = [0 S(3) 0 S(4)];
end

files=dir([diry,name]);
  
BM_TIME = get_variable('ocean_time',name,diry);
BM=[];
for jj=1:length(files)
    Fname = files(jj).name;
    disp(Fname);

    S = get_varsize(Fname,'mudmass_04');
    if(F==1)  % fluvial sedimnt
        bm=nc_varget(Fname,'mudmass_04',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        bm=bm+nc_varget(Fname,'mudmass_05',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        bm=nc_varget(Fname,'mudmass_06',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        bm=bm+nc_varget(Fname,'mudmass_07',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        
    else      % bed sediment
        bm=nc_varget(Fname,'mudmass_01',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        bm=bm+nc_varget(Fname,'mudmass_02',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
        bm=bm+nc_varget(Fname,'mudmass_03',[0 0 rect(1) rect(3)],...
            [S(1) S(2) rect(2)-rect(1) rect(4)-rect(3)]);
    end
    
    bm2=squeeze(sum(bm,2));

    
    BM=[BM;bm2];
end

%if(total_mass)
%pm=grid{'pm'}(sb:nb,wb:eb);   %x
%pn=grid{'pn'}(sb:nb,wb:eb);   %y
%dx=1./pm;
%dy=1./pn;
%area2=dx.*dy;
%clear dx dy pm pn
%end

if(F==1)
   save bed_massR.mat BM BM_TIME rect
else
   save bed_massB.mat BM BM_TIME rect
end

