%%%
%
% To get the mass of sediment on the bed from a number of ROMS output
% files.
%
%       This file is outdated, but uses the older netcdf toolbox notation
%       and can be used as an example.
%          Aaron Bever
%
river_sed = 1;
total_mass=1;
start = [2006 1 1 0 0 0];    % must be in gregrorian dates
stop = [2006 9 1 0 0 0];

diry='/data/users-tmp/abever/MODELS/WAIPAOA2/OUTPUT/FEB2009_B/'
files=dir([diry,'*his*.nc']);

if(total_mass)
grid=netcdf('/data/users-tmp/abever/MODELS/WAIPAOA2/FORCING/INP_3AUG/modern_test_grid.nc');
end

nb = 180;  % northern boundary
sb = 1;   % southern boundary
wb = 1;
eb = 115;   % eastern boundary

TIME=[];
BM=[];
for jj=1:length(files)
    disp(jj)
    Fname = files(jj).name;
    nc = netcdf(Fname);
    time = nc{'ocean_time'}(:);
    time = roms_timeA(time,1);
    ind(jj).t = tind_L(time,start,stop);  % getting time indicators
    TIME = [TIME;time(ind(jj).t)];
    
    if(river_sed)
        bm=nc{'mudmass_04'}(ind(jj).t,:,sb:nb,wb:eb);
        bm=bm+nc{'mudmass_05'}(ind(jj).t,:,sb:nb,wb:eb);
    else
        bm=nc{'mudmass_01'}(ind(jj).t,:,sb:nb,wb:eb);
        bm=bm+nc{'mudmass_02'}(ind(jj).t,:,sb:nb,wb:eb);
        bm=bm+nc{'mudmass_03'}(ind(jj).t,:,sb:nb,wb:eb);
    end
    
    if(length(ind(jj).t)==1)
        bm=squeeze(sum(bm,1));
        bm2(1,:,:)=bm(:,:);
    else
        bm=squeeze(sum(bm,2));
        bm2=bm;
    end
    
    BM=[BM;bm2];
    close(nc)
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

    BM_2D=BM;
    clear BM
    save bed_massI.mat BM_2D
    for i=1:length(BM_2D(:,1,1))
        % out of memory here
        BM_2D(i,:,:)=squeeze(BM_2D(i,:,:)).*area2;
    end

    BMT_1D=squeeze(sum(sum(BM_2D,2),3));
    BM_TIME=TIME;
    load bed_massI;
    if(river_sed)
        save bed_mass_river.mat BMT_1D BM_2D BM_TIME nb sb eb wb
    else
        save bed_mass_bed.mat BMT_1D BM_2D BM_TIME nb sb eb wb
    end

