function AnalyseXBeachOutput_Morpho(info, dirs)
    test_name = info.run;

    %% Notes
    % Morfologische analyse output

    %% ------------------------------------------------------------------ %%
    %  ---------------------------- Settings ---------------------------- %
    %  ------------------------------------------------------------------ %
    dxmin = 1;        % minimale gridresolutie[m]
    
    % Hoogste offshore waterstand tijdens storm
    zs      = load([dirs.output filesep 'tide.txt']); 
    zs0     = max(zs(:,2)); 
    
    %% CALCULATE DUNE EROSION STATS IN MEASURED PROFILES
    
    VeroMeas = [];
    retreatMeas = [];  
    BermSlopeMeas = [];
    
    %% READ AND PROCESS XBEACH OUTPUT
    ncfile  = [dirs.output filesep 'xboutput.nc']; 
    errfile = [dirs.output filesep 'XBerror.txt'];
    
    if isfile(ncfile) %&& ~isfile(errfile)
        x           = squeeze(ncread(ncfile,'globalx'));
        zb          = squeeze(ncread(ncfile,'zb'));

        %% Calculate dune erosion statistics       
        [VeroXB, ~, ~, ~, d1_XB, z1_XB, d2_XB, z2_XB] = compute_Vero2(x, zb(:,1), x, zb(:,end), zs0, 0); 
        BermSlopeXB = compute_bedslope_berm_based_on_accretion_zone(x,zb(:,1),zb(:,end),zs0);
        retreatXB   = compute_retreat2(x, zb(:,1), zb(:,end), dxmin, zs0+1);

        height = zs0:0.1:zs0+2;
        for ii = 1:length(height)
        retreatXB2(ii,1) = compute_retreat2(x, zb(:,1), zb(:,end), dxmin, height(ii));
        end
        slretreat = retreatXB2(1); % storm level retreat
        maxretreat = max(retreatXB2); clear retreatXB2
        
        %% Plot figure cross-section
        
        % Let op: zelf aanpassen: xlims, referentieniveau in label y-as en tekst met retreat distance, hoogte mean water level
        
        % voor xlims
        idx_init    = find(zb(:,1)>zs0, 1, 'first');
        idx_end     = find(zb(:,end)>zs0, 1, 'first');        
        
        figure; 
        set(gcf, 'Position', [300 200 900 600])
        subplot('Position', [0.06 0.45 0.9 0.5]); hold on; grid on; box on;          
        title([test_name], 'interpreter', 'none')
        xlim([x(idx_init)-120 x(idx_end)+30])  % Oorspronkelijk  % Range zelf instellen
            
        xlims = xlim;
        % Erosievolume 
        h1 = fill([d1_XB; flipud(d2_XB)],  [z1_XB; flipud(z2_XB)], 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'displayname', 'Dune erosion XBeach'); % [255 200 0]./255 = oranjegeel
%        h2 = fill([d1;    flipud(d2)],     [z1; flipud(z2)],       'k', 'FaceAlpha', 0.1, 'EdgeColor', 'none', 'displayname', 'Dune erosion measured');
        % Max. waterstand
        h3 = plot(xlims,  [zs0, zs0],  '-b', 'linewidth', 1.2, 'DisplayName', 'Max. storm surge level');
%             h7 = plot(xlims,  [mean(zs(:,2)), mean(zs(:,2))],':b', 'linewidth', 1.2, 'DisplayName', 'Mean water level'); % Let op: case-specifiek. Dit is gemiddelde over enkele maanden rondom de storm in TAW.
        % Profielen
        h5 = plot(x,                zb(:,end),   '-r', 'linewidth', 1.2, 'DisplayName', 'XBeach post-storm bed level');
%        h6 = plot(-Meas_PostStorm.x(:,1), Meas_PostStorm.z(:,1), '-k', 'linewidth', 1.2, 'DisplayName', 'Measured post-storm bed level');
        h4 = plot(x,                zb(:,1),     '--k', 'linewidth', 1.2, 'DisplayName', 'Initial bed level');
        % Tekstblok met duinerosie-eigenschappen
        ylims = ylim;   
        txty  = ylims(1)+0.90.*(ylims(2)-ylims(1));     
        txtdy = 0.05.*(ylims(2)-ylims(1));   
        txtx  = xlims(1)+0.35.*(xlims(2)-xlims(1));
        text(txtx,txty,        ['Dune erosion volume - XBeach            = '    num2str(VeroXB, '%.0f')         ' m^3/m']);           
        text(txtx,txty-txtdy,  ['Dune erosion volume - measured         = 64 - 116 m^3/m']); 
        try text(txtx,txty-2*txtdy,['Dune retreat at' num2str(zs0+1,'%1.1f') 'm +NAP - XBeach    = '       num2str(retreatXB, '%.0f')      ' m']  ); end
        text(txtx,txty-3*txtdy,['Dune retreat at' num2str(zs0+1,'%1.1f') 'm +NAP - measured = 15 - 18 m']  );
                      
        leg = legend([h4 h5 h1 h3], 'location', 'northwest');
        set(gca, 'FontSize', 11)

        xlabel('Distance along transect [m]');
        ylabel('Elevation [m +NAP]');

        subplot('Position', [0.06 0.1 0.9 0.25]); hold on; grid on; box on;
        % Max. waterstand
        h3 = plot([min(x) max(x)],  [zs0, zs0],  '-b', 'linewidth', 1.2, 'DisplayName', 'Max. storm surge level');
%       h7 = plot([min(x) max(x)],  [mean(zs(:,2)), mean(zs(:,2))],'--b','linewidth', 1.2, 'DisplayName', 'Mean water level'); % Let op: case-specifiek. Dit is gemiddelde over enkele maanden rondom de storm in NAP.
        % Profielen
        h5 = plot(x,                zb(:,end),   '-r', 'linewidth', 1.2, 'DisplayName', 'XBeach post-storm bed level');
%        h6 = plot(Meas_PostStorm.x(:,1), Meas_PostStorm.z(:,1), '-k', 'linewidth', 1.2, 'DisplayName', 'Measured post-storm bed level');
        h4 = plot(x,                zb(:,1),     '--k', 'linewidth', 1.2, 'DisplayName', 'Initial bed level');
        % Zoombox for upper subplot 
        rectangle('position',[xlims(1) ylims(1) xlims(2)-xlims(1) ylims(2)-ylims(1)+0.5], 'FaceColor', 'none', 'EdgeColor', 'k', 'Linestyle', '-.', 'LineWidth', 0.5)
                
        set(gca, 'FontSize', 11)
        xlim([min(x) max(x)])
        xlabel('Distance along transect [m]');
        ylabel('Elevation [m +NAP]');
        
        xb_write_plot(gcf, dirs.output, ['BL_BeforeAndAfterStorm_' test_name])
        save([dirs.output 'values.mat'], 'VeroXB','VeroMeas', 'retreatXB','retreatMeas'); 
        
%         dlm = {'VeroXB','VeroMeas', 'retreatXB','retreatMeas'};
%         matrix2latex([VeroXB, VeroMeas, retreatXB, retreatMeas], 'filename', 'skill.tex', ...
%             'rowjustification','p{1.5cm}',...
%             'caption',  'Erosion volume and dune retreat at 3m+NAP for the Egmond case', ...
%             'rowlabel', test_name, ...
%             'collabel', dlm, ...
%             'format',	{'%s', '%4.2f', '%4.2f','%4.2f','%4.2f'});
    end  
end    

   

