function sed_transport()
%% coefficients and settings
% [runid,testid,datadir]=testinfo
par.morfac = 1;
par.np = 0.4;
par.tint = 1;

addpath P:\1002266-xbeach\code

dirs = {'run000'};
eval(['cd ',dirs{1}]); 
XB(1) = getdimensions;
cd ..

%% simulation results
i=1;
xw  = XB(i).xc;
eval(['cd ',dirs{i}]); 
s.zb = readvar('zb.dat',XB(i),'2D');    
s.zb = s.zb(:,2,:); 
s.dzav = readvar('dzav.dat',XB(i),'2D');  
s.dzav = s.dzav(:,2,:); 
cd ..


% load measurements
xm = load('x.asc'); 
zm = load('z.asc'); 
measured = [1 2 3 4];
extrapoints = unique([min(xm(:,measured)) max(xm(:,measured))]);
xwm = sort([xw(:,2); extrapoints']);
temp=zm(~isnan(zm(:,1)));
temp2=round(length(xwm(:,1))/2);temp3=length(xwm(:,1))-temp2;
s.Sdzm = [repmat(0,temp2,3); repmat(temp(end),temp3,3)];
for i=1:4;
    s.zbm(:,i) = interp1(xm(~isnan(xm(:,measured(i))),measured(i)),zm(~isnan(xm(:,measured(i))),measured(i)),xwm(:,1));
    s.Sdzm(~isnan(s.zbm(:,i)),i) = s.zbm(~isnan(s.zbm(:,i)),i);
end

temp=0.5.*diff(xw(:,2));temp2=[temp(2:end);0];
dx = ([temp(1);temp+temp2]);
temp=0.5.*diff(xwm(:,1));temp2=[temp(2:end);0];
dxm = ([temp(1);temp+temp2]);
% Get interval averaged values
Tout = round([0 0.2 0.6 6]*3600/par.morfac/par.tint)+1;
for i = 1:length(Tout)-1
    dz = s.zb(:,:,Tout(i+1))-s.zb(:,:,Tout(i));
    dzav = s.dzav(:,:,Tout(i+1))-s.dzav(:,:,Tout(i));
    dzm = s.Sdzm(:,i+1)-s.Sdzm(:,i);
    % sediment transport from bed level changes
    s.Sdzb(:,i) = (1-par.np)*flipud(cumsum(flipud(dz.*([0.5*diff(xw(1:2,2));diff(xw(:,2))]))))/((Tout(i+1)-Tout(i))*par.morfac*par.tint);
    % sediment transport due to avalanching
    s.Sdzav(:,i) = (1-par.np)*flipud(cumsum(flipud(dzav.*(dx))))/((Tout(i+1)-Tout(i))*par.morfac*par.tint);
    % sediment transport from bed level changes in lab
    s.Sdzbm(:,i) = (1-par.np)*flipud(cumsum(flipud(dzm.*(dxm))))/((Tout(i+1)-Tout(i))*par.morfac*par.tint);
end

%% figures;
% sediment transports in XBeach
figure();
for i = 1:length(Tout)-1
    subplot(2,2,i)
    % transports from bed level changes
    plot(xw(:,2),s.Sdzb(:,i),'k','LineWidth',1.5), hold on 
    plot(xwm(:,1),s.Sdzbm(:,i),'r','LineWidth',2)
    plot(xw(:,2),s.Sdzav(:,i),'k--','LineWidth',2); 
    plot(xw(:,2),(s.Sdzb(:,i)-s.Sdzav(:,i)),':b','LineWidth',2)
    xlabel('x [m]'); ylabel('S [m^3/m/s]');
    axis([85 max(max(xw)) 1.25*min(s.Sdzb(:,i)) 1.25*max(s.Sdzb(:,i)-s.Sdzav(:,i))+1e-6]);
    legend('S_{\Deltazb,XB}','S_{\Deltazb,m}','S_{avalanch}','S_{tot,XBeach}',3);
end
% pname = ['sediment_transports_intervals.jpg'];
% eval(['print -djpeg ' pname]);
end