function analyze_H298_runs(info, dirs)

%% boundary conditions and measurements:

Trep = [5.4 5.4 5.4];
Tout = {[0 0.6 1.3 3.5 9.5];  %T1
        [0 1.3 3.5 9.5 16];   %T2
        [0 1.3 3.5 9.5 16];}; %T3

zrunup = [0.0  0.25 0.5  0.75 1.0  1.25 1.5  1.75 2.0  2.25 2.5  2.75 3.00;
          0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4;
          0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4;
          0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4;];

runupm = load([dirs.data filesep 'runup.txt']);

zrunup2 = 0:0.05:2.4;

% wave height measurements:

x_Hm0 =    {[15 150 15 172 15 15 150 15 15 187 15 172;                                                   %T1
            1.70 1.25 1.72 1.18 1.70 1.72 1.34 1.70 1.75 0.88 1.70 1.09;];
            [15   187  15   15   150  15   172  15   187  15   150  15   40;                             %T2
             1.70 0.94 1.64 1.65 1.18 1.65 1.07 1.63 0.79 1.64 1.16 1.67 1.48;];
            [15   15   15   15   15   15   150  15   172  15   187  15   150  15   40;                   %T3
             1.63 1.64 1.65 1.68 1.64 1.65 0.92 1.67 0.96 1.67 0.70 1.67 1.09 1.69 1.27;];};

%% Plot settings

% Deltares colors
Dcols = [0 139 191;
         128 169 189;
         235 0 115;
         0 196 214;
         180 211 0;
         255 114 0;
         160 158 141];         
         
first = 1;

%% Analysis loop

id = str2double(strrep(info.run,'T',''));

xbi     = xb_read_params(pwd);
tint    = xs_get(xbi,'tint');

eps     = 0.005;
hswitch = 0.1;

indp    = round(Tout{id}*3600/tint);
indp(1) = 1;

% plot colors
if first == 1
    colors  = colormap(jet(length(indp)));
    close;
    first   = 0;

    xp      = 1-0.0001;
    rxp     = 2.5;
    b       = sqrt(-rxp.^2/(2*log(1-xp)));
end

% load model output
xbo         = xb_read_output('xboutput.nc','vars',{'zs','zb','dzav','H','xrunup','Hrunup','runup','u','DIMS','H_mean','Susg_mean'});
xb_output   = xs_peel(xbo);
xb_dims     = xs_peel(xb_output.DIMS);

%squeeze
xb_output.zs    = squeeze(xb_output.zs);
xb_output.zb    = squeeze(xb_output.zb);
xb_output.H    = squeeze(xb_output.H);
xb_output.u    = squeeze(xb_output.u);
xb_output.dzav    = squeeze(xb_output.dzav);
xb_output.H_mean    = squeeze(xb_output.H_mean);

% get profile without short wave runup.
xbo2 = load([dirs.data filesep 'old_results_' info.run '.mat']);
zb_old = xs_get(xbo2.xbo,'zb');

hh = squeeze(xb_output.zs-xb_output.zb);

wetz = zeros(size(hh));
wetz(hh>=hswitch) = 1;
[ind,dum] = find(diff(wetz')==-1);

% load hard layer
hardlayer = load('nebed.dep');
indmin = min(find(hardlayer<=0.05));
indmax = max(find(hardlayer<=0.05));
zbh = squeeze(xb_output.zb(1,:))-hardlayer(1:xb_dims.nx+1);

% initialize
istruct     = NaN(size(xb_output.zb,1),1);
istruct2    = NaN(size(xb_output.zb,1),1);
zsscour     = NaN(size(xb_output.zb,1),1);

% scour depth
for j = 1:size(xb_output.zb,1)
    dz = xb_output.zb(j,:)-zbh;
    idxNoEro = find(dz<0.01);
    istruct(j) = min(idxNoEro(idxNoEro>10))-1;
    dz = diff(xb_output.zb(j,:));
    istruct2(j) = min(find(dz == min(abs(dz))));
    zsscour(j) = min(0, min(xb_output.zb(j,istruct2(j)-2:1:min(istruct2(j)+10,xb_dims.nx)))-xb_output.zb(1,istruct(1)));
end

% check on availability data
indpp = indp(find(indp<=size(xb_output.zb,1)));
time = [0:1:length(istruct)-1]*tint/3600;

% erosion above storm surge level:
temp = xb_dims.globalx_DATA(2:end)-xb_dims.globalx_DATA(1:end-1);
dsu = xb_dims.globalx_DATA*0;
dsu(2:end-1) = 0.5*(temp(1:end-1)+temp(2:end));
dsu(1) = temp(1); dsu(2) = temp(end);
hav = zeros(size(squeeze(xb_output.zb)));
zsr3 = squeeze(xb_dims.point_zs_DATA);
for  j = 1:size(xb_output.zb,1)
    A(j) = sum((max(xb_output.zb(1,:),0)-max(xb_output.zb(j,:),0)).*dsu');
    ind_wl(j) = min(find(squeeze(wetz(j,:,:))==0))-1;
    x_wl(j) = xb_dims.globalx_DATA(ind_wl(j));
    zsr2(j) =  xb_output.zs(j,ind_wl(j));
    ur2(j) = xb_output.u(j,ind_wl(j));
    hav(j,:) = squeeze(wetz(j,:).*(hh(j,:) + xb_output.runup(j,1))) + squeeze(1.d0-wetz(j,:)).*max(eps,xb_output.runup(j,1)+squeeze(xb_output.zs(j,ind_wl(j))-xb_output.zb(j,:)));
end
hav = max(hav,0);
zsr = squeeze(xb_output.zb)+hav;

% get runup measurements
ind = find(runupm(:,1)==id);

%% runup analyis
for j = 1:length(ind)

    % make rayleigh scale
    y = rayl_cdf(zrunup(id,:),b);
    yaxis = [100:-10:50 40 30 20 10 5 4 3 2 1 0.1 0.01];
    yscale = (100-yaxis)/100;
    xtemp = rayl_inv(yscale,b);

    %%% remark: first need to find zero crossings in signal and define
    %%% runup per wave

    % get runup simulations
    time_ind = ceil([runupm(ind(j),2) runupm(ind(j),3)]*3600/tint);
    time_ind = max(1,time_ind);

    % short + long wave runup
    temp1 = zsr2(time_ind(1):time_ind(2))'+xb_output.runup(time_ind(1):time_ind(2)); % select relevant part of timeseries
    % select relevant part of timeseries
    idx1 = signal_findpeaks(temp1); % define local maxima as runup heights
    rHs_tot = temp1(idx1);
    locs1 = find(idx1);
    nr_rtot = histc(rHs_tot,zrunup2);
    nr_rtot = cumsum(nr_rtot);
    nr_rtot = 100*(1-nr_rtot/nr_rtot(end));
    % short wave runup
    temp2 = xb_output.runup(time_ind(1):time_ind(2));
    idx2 = signal_findpeaks(temp2);
    rHs_sw = temp2(idx2);
    locs2 = find(idx2);
    nr_rsw = histc(rHs_sw,zrunup2);
    nr_rsw = cumsum(nr_rsw);
    nr_rsw = 100*(1-nr_rsw/nr_rsw(end));
    % long wave runup
    temp3 = zsr2(time_ind(1):time_ind(2));      % select relevant part of timeseries
    % define local maxima as runup heights
    idx3 = signal_findpeaks(temp3);
    rHs_lw = temp3(idx3);
    locs3 = find(idx3);
    nr_rlw = histc(rHs_lw,zrunup2);
    nr_rlw = cumsum(nr_rlw);
    nr_rlw = 100*(1-nr_rlw/nr_rlw(end));
    % short wave height distribution at toe structure:
    temp4 = xb_output.Hrunup(time_ind(1):time_ind(2));
    nr_rHs = histc(temp4,zrunup2);
    nr_rHs = cumsum(nr_rHs);
    nr_rHs = 100*(1-nr_rHs/nr_rHs(end));

    % plot results
    figure; set(gca,'FontSize',16);
    if id~=5
        xv = rayl_inv((100-runupm(ind(j),6:end))./100,b);
        plot(xv,zrunup(id,:),'Color',Dcols(7,:)/255,'LineWidth',1.5); hold on;
    end
    xv = rayl_inv((100-nr_rtot)./100,b);
    plot(xv,zrunup2,'-o','Color',Dcols(1,:)/255,'LineWidth',1.5,'MarkerSize',8); hold on;
    xv = rayl_inv((100-nr_rsw)./100,b);
    plot(xv,zrunup2,'-^','Color',Dcols(5,:)/255,'LineWidth',1.5,'MarkerSize',8);
    xv = rayl_inv((100-nr_rlw)./100,b);
    plot(xv,zrunup2,'-v','Color',Dcols(3,:)/255,'LineWidth',1.5,'MarkerSize',8);
    set(gca,'XTick',xtemp);
    set(gca,'XTicklabel',num2str(yaxis'));
    xlabel('P [%]'); ylabel('R [m]');
    if id~=5
        legend('R_{tot,measured}','R_{tot}','R_{hf}','R_{lf}')
    else
        legend('R_{tot}','R_{hf}','R_{lf}')
    end
    title(['runup statistcs between t = ',num2str(runupm(ind(j),2)),' till ',num2str(runupm(ind(j),3)),' [hours]']);
    ylim([0 2]);
    grid on;
    print([dirs.output filesep 'runup_statistics_',info.run,'_at_t',num2str(time_ind(1)*10/3600),'hours.png'],'-dpng');
end

% runup statistics
temp1 = zsr2'+xb_output.runup;
idx1 = signal_findpeaks(temp1);
rHs_tot = temp1(idx1);
locs1 = find(idx1);
temp2 = xb_output.runup;
idx2 = signal_findpeaks(temp2);
rHs_sw = temp2(idx2);
locs2 = find(idx2);
temp3 = zsr2;
idx3 = signal_findpeaks(temp3);
rHs_lw = temp3(idx3);
locs3 = find(idx3);

% runup indexes
for j = 1:length(temp1)
    ind_wl2(j) = min(find(temp1(j)-squeeze(xb_output.zb(j,:))<=0.02))-1;
end

% more detailed runup analysis
if id==3
    % fig1: show that long wave runup does not exceed height revetment
    figure;
    subplot(211); set(gca,'FontSize',16);
    plot(xb_dims.globaltime_DATA/3600,zsr2,'Color',Dcols(1,:)/255,'LineWidth',1.5); hold on;
    plot([0 17],[0.6 0.6],'k--','LineWidth',3);
    xlabel('time [hours]'); ylabel('R [m]');
    axis([6.2 6.4 -0.5 1.1]);
    legend('R_{lf}','top revetment');
    pname = [dirs.output filesep,info.run,'_long_wave_runup','.jpg'];
    eval(['print -djpeg ' pname]);

    % fig2: decay runup height a test progresses
    windowSize = 300;
    rHs_tot_sm = filter(ones(1,windowSize)/windowSize,1,rHs_tot);
    rHs_tot_sm(1:300)=rHs_tot_sm(300);

    figure; set(gca,'FontSize',16);
    plot(xb_dims.globaltime_DATA/3600,temp1','Color',Dcols(1,:)/255); hold on;
    plot(xb_dims.globaltime_DATA(locs1)/3600,rHs_tot','Color',Dcols(7,:)/255,'LineWidth',1.5);
    plot(xb_dims.globaltime_DATA(locs1)/3600,rHs_tot_sm','Color',Dcols(3,:)/255,'LineWidth',3);
    axis([1000/3600 16+1000/3600 -1 1.5]); grid on;
    xlabel('time [hours]'); ylabel('R [m]');
    legend('instantaneous runup elevation','runup elevation','trend in runup elevation');
    pname = [dirs.output filesep,info.run,'_runup_trend','.jpg'];
    eval(['print -djpeg ' pname]);

    % fig3: show short wave runup
    runupsw = 0.8*xb_output.Hrunup*2.3;

    figure;
    subplot(211); set(gca,'FontSize',16);
    plot(xb_dims.globaltime_DATA/3600,xb_output.Hrunup','Color',Dcols(1,:)/255,'LineWidth',1.5); hold on;
    plot(xb_dims.globaltime_DATA/3600,runupsw','Color',Dcols(5,:)/255,'LineWidth',1.5);
    axis([6.30 6.375 0 1]);
    xlabel('time [hours]'); ylabel('H_{rms} & R [m]');
    legend('H_{rms} at revetments base','maximum R_{hf}');
    subplot(212); set(gca,'FontSize',16);
    plot(xb_dims.globaltime_DATA/3600,runupsw','Color',Dcols(5,:)/255,'LineWidth',1.5); hold on;
    plot(xb_dims.globaltime_DATA/3600,xb_output.runup','Color',Dcols(1,:)/255,'LineWidth',1.5);
    axis([6.30 6.375 -1 1]);
    xlabel('time [hours]'); ylabel('R_{hf} [m]');
    legend('maximum R_{hf}','instantaneous R_{hf}');
    pname = [dirs.output filesep,info.run,'_short_wave_runup','.jpg'];
    eval(['print -djpeg ' pname]);

    % fig4: show total runup
    figure;
    subplot(211)
    plot(xb_dims.globaltime_DATA/3600,xb_output.runup','Color',Dcols(5,:)/255,'LineWidth',1.5); hold on;
    plot(xb_dims.globaltime_DATA/3600,zsr2,'Color',Dcols(1,:)/255,'LineWidth',1.5);
    plot([0 17],[0.6 0.6],'k--','LineWidth',3);
    axis([6.30 6.375 -1 1.1]);
    ylabel('R [m]');
    legend('R_{hf}','R_{lf}','top revetment');
    subplot(212)
    plot(xb_dims.globaltime_DATA/3600,zsr2+xb_output.runup','Color',Dcols(1,:)/255,'LineWidth',1.5); hold on;
    plot([0 17],[0.6 0.6],'k--','LineWidth',3);
    axis([6.30 6.375 -1 1.1]);
    xlabel('time [hours]'); ylabel('R [m]');
    legend('R_{tot}','top revetment');
    pname = [dirs.output filesep,info.run,'_total_runup','.jpg'];
    eval(['print -djpeg ' pname]);
end

%%  morphology

% measurements
load([dirs.data filesep 'data_' info.run '.mat']);
zt = z; 
clear z;
z0 = xb_output.zb(1,1);
z = zt;

% intrepolate measured profiles to xbgrid
for j = 1:min(length(z),length(indpp))
    zm = interp1(x,z{j}+z0-z{j}(1),xb_dims.globalx_DATA)';
    zmscour(j) = min(0, min(zm(istruct2(indpp(j))-5:1:istruct2(indpp(j))+1))-xb_output.zb(1,istruct(1)));
end

%% visualize profile evolution

% old approach
figure; set(gca,'FontSize',16);
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:)),'-','Color',Dcols(7,:)/255,'LineWidth',1.5); hold on;for j = 2:length(z)
    plot(x,z{j}+(z0-z{j}(1)),'-','LineWidth',1.5,'Color',Dcols(7,:)/255); hold on;
end

for j = 2:length(z)
    plot(xb_dims.globalx_DATA,squeeze(zb_old(indp(j),1,:)),'k--','LineWidth',1.5); hold on;
end
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:))-hardlayer(1:xb_dims.nx+1)','Color',Dcols(7,:)/255,'LineWidth',5,'MarkerSize',4); hold on;
% wave heights
plot(x_Hm0{id}(1,:),x_Hm0{id}(2,:)/sqrt(2),'o','MarkerEdgeColor',Dcols(4,:)/255,'LineWidth',1.5); hold on;
plot(xb_dims.globalx_DATA,squeeze(xb_output.H_mean),'-','Color',Dcols(4,:)/255,'LineWidth',1.5);
axis([max(x)-70 max(x) -2 3]); xlabel('x [m]'); ylabel('z_b [m]');
pname = [dirs.output filesep,info.run,'_profile_evolution_old','.jpg'];
eval(['print -djpeg ' pname]);

% new approach
figure; set(gca,'FontSize',16);
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:)),'-','Color',Dcols(7,:)/255,'LineWidth',1.5); hold on;
for j = 2:length(z)
    plot(x,z{j}+(z0-z{j}(1)),'-','LineWidth',1.5,'Color',Dcols(7,:)/255); hold on;
end

for j = 2:length(z)
    plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(indp(j),1,:)),'k--','LineWidth',1.5); hold on;
end

plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:))-hardlayer(1:xb_dims.nx+1)','-','Color',Dcols(7,:)/255,'LineWidth',5,'MarkerSize',4); hold on;
% wave heights
plot(x_Hm0{id}(1,:),x_Hm0{id}(2,:)/sqrt(2),'o','MarkerEdgeColor',Dcols(4,:)/255,'LineWidth',1.5); hold on;
plot(xb_dims.globalx_DATA,squeeze(xb_output.H_mean),'-','Color',Dcols(4,:)/255,'LineWidth',1.5);
axis([max(x)-70 max(x) -2 3]); xlabel('x [m]'); ylabel('z_b [m]');
pname = [dirs.output filesep,info.run,'_profile_evolution','.jpg'];
eval(['print -djpeg ' pname]);

% remove all but the first and last outputs
for iv = 2:length(xbo.data)
    if length(size(xbo.data(iv).value)) == 3
        xbo.data(iv).value = xbo.data(iv).value([1 end],:,:);
    end
end

zbh = max(zbh,-2);
xb_plot_profile(xbo,'measured',[x' z{end}'+z0], 'nonerodible', [xb_dims.globalx_DATA zbh'])
title(['profiles ' info.test ' ' info.run])
xb_write_plot(gcf, dirs.output, 'final_profiles')

%% visualize scour depth evolution

figure; set(gca,'FontSize',16);
plot(Tout{id},zmscour,'-','Color',Dcols(7,:)/255,'LineWidth',1.5); hold on;
plot(time,zsscour,'-','Color',Dcols(1,:)/255,'LineWidth',1.5);
xlabel('Time [hours]'); ylabel('Scour Depth [m]')
title('Scour depth w.r.t. initial bed level at revetments toe');
legend('observed','simulated');
pname = [dirs.output filesep,info.run,'_scour_evolution','.jpg'];
eval(['print -djpeg ' pname]);

%% erosion above SSL

% measured erosion
dxm = diff(x); dxm(1) = 0.5*dxm(1); dxm(end+1) = 0.5*dxm(end);

Am = zeros(size(z));
for j = 2:length(z)
    Am(j) = -sum((max(z{j},-z0)-max(z{1},-z0)).*dxm);
end

% store all erosion volumes in on earray
Amall = Am;
Asall = A';
tAall = xb_dims.globaltime_DATA/3600;

figure; set(gca,'FontSize',16);
plot(xb_dims.globaltime_DATA/3600-1000/3600,A','Color',Dcols(1,:)/255,'LineWidth',2); hold on;
plot(Tout{id},Am,'*','MarkerEdgeColor',Dcols(1,:)/255,'LineWidth',2);
xlabel('Time [hours]')
ylabel('Erosion Volume [m^3/m]')
axis([0 Tout{id}(end) 0 1.2*max(Am)]);
pname = [dirs.output filesep,info.run,'_erosion_evolution','.jpg'];
eval(['print -djpeg ' pname]);

%% visualize sediment transports

% Sediment transport from profile evolution
dx  = diff(xb_dims.globalx_DATA); dx(1) = 0.5*dx(1); dx(end+1)=0.5*(dx(end));
dz  = -squeeze(xb_output.zb(end,1,:)-xb_output.zb(1,1,:));
dz2 = -squeeze(xb_output.dzav(end,1,:));
dzm = z{end}-z{1};

Ss  = fliplr(cumsum(fliplr(dz.*dx))*(1-0.4));
Ss2  = fliplr(cumsum(fliplr(dz2.*dx))*(1-0.4));

Ss = Ss/(xb_dims.globaltime_DATA(end)-xb_dims.globaltime_DATA(1));
Ss2 = Ss2/(xb_dims.globaltime_DATA(end)-xb_dims.globaltime_DATA(1));

figure;
subplot(2,1,1); set(gca,'FontSize',16);
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:)),'k:','LineWidth',2); hold on;
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(end,1,:)),'k-','LineWidth',2);
plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:))-hardlayer','k-o','LineWidth',2,'MarkerSize',4);
axis([150 200 -2.5 3]); ylabel('z_b [m]');
legend('initial profile','final profile','revetment');
subplot(2,1,2); set(gca,'FontSize',16);
%transports from bed level changes
plot(xb_dims.globalx_DATA,Ss,'k--','LineWidth',2); hold on;
plot(xb_dims.globalx_DATA,Ss2,'Color',Dcols(1,:)/255,'LineWidth',2);
% hydraulic transpports
plot(xb_dims.globalx_DATA,squeeze(xb_output.Susg_mean(end,1,:)),'Color',Dcols(3,:)/255,'LineWidth',2); hold on;
xlabel('x [m]'); ylabel('S [m^2/s]');
legend('S_{\deltaz_b}','S_{avalanching}','S_{suspended}');
xlim([150 200]);

pname = [dirs.output filesep,info.run,'_sediment_transport','.jpg'];
eval(['print -djpeg ' pname]);

close all;

clearvars -except Trep Tout zrunup runupm zrunup2 x_Hm0 first id i b colors Dcols Amall Asall tAall

% movie = 0;
% 
% zhl = squeeze(zb(1,1,:))-hardlayer(1,:)';
% 
% figure;
% for i =1:1:3600
%     plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(1,1,:)),'--','Color',Dcols(7,:)/255); hold on;
%     plot(xb_dims.globalx_DATA,zhl,'Color',Dcols(7,:)/255,'LineWidth',5);
%     plot(xb_dims.globalx_DATA,squeeze(xb_output.zs(i,1,:)),'Color',Dcols(1,:)/255);
%     plot(xb_dims.globalx_DATA,squeeze(zb_old(i,1,:)),'-','Color',Dcols(6,:)/255,'LineWidth',1.5);
%     plot(xb_dims.globalx_DATA,squeeze(xb_output.zb(i,1,:)),'-o','Color',Dcols(3,:)/255,'LineWidth',1.5); 
%     plot(xb_dims.globalx_DATA,squeeze(xb_output.H(i,1,:)),'Color',Dcols(5,:)/255);
%     plot([squeeze(xb_dims.xHrunup_DATA(i)) squeeze(xb_dims.xHrunup_DATA(i))],[-4.2 2.5],'--','Color',Dcols(6,:)/255,'LineWidth',2);
%     plot(xb_dims.globalx_DATA(ind_wl(i)),xb_output.zb(i,ind_wl(i)+1),'o','MarkerEdgeColor',Dcols(3,:)/255,'LineWidth',3,'MarkerSize',12);
%     axis([170 210 -2 2.5]); xlabel('x [m]'); ylabel('z_b [m]');
%     title(['t = ',num2str(xb_dims.globaltime_DATA(i)/3600),' [hours]']);
%     if movie==1
%         name = ['results\avi\hard_',num2str(i,'%04.0f'),'.png'];
%         print('-dpng','-r300',name);
%         pause(0.1); hold off;
%     else
%         pause(0.1); hold off;
%     end
% end
% print endprofile.png -dpng
% pname = ['..\..\report\',testid '_' runid '_fig1' '.jpg'];
% eval(['print -djpeg ' pname]); 