function [dVero, x_all, z1_all, z2_all, x1_ero, z1_ero, x2_ero, z2_ero] = compute_Vero2(x1, z1, x2, z2, zs0, flag) 
    % INPUT
    % x1 = distances along initial profile 
    % z1 = elevation along initial profile
    % x2 = distances along final profile
    % z2 = elevation along final profile
    % zs0 = max. storm surge level (above this level, the erosion volume will be determined)
    % flag (optional): 
    %   flag = 1: only first erosion zone above zs0 (in the first dune row) will be used to calculate erosion volume
    %   flag = 0 (default): no landward limit to calculate erosion volume
    %   NOTE: flag = 1 does not work in all cases, requires check on output

    % OUTPUT
    % dVero = erosion volume 
    % @@_all = original profiles, but including cross-points of initial and final profile
    % x1_ero,z1_ero,x2_ero,z2_ero = sections of the profiles used for calculating the erosion volume 
    %                               (can be used for plotting the erosion area with fill([x1; flipud(x2)], [z1; flipud(z2)])
    
    %% CALCULATIONS
    % Zorg voor goede format van profielen
    if size(x1,1) > 1;      x1 = x1';   z1 = z1';   end         % make row vectors
    if size(x2,1) > 1;      x2 = x2';   z2 = z2';   end
    if mean(z1(1:6)) > mean(z1(end-6:end));     z1 = fliplr(z1);     x1 = fliplr(x1);  end  % from sea to land (low to high)
    if mean(z2(1:6)) > mean(z2(end-6:end));     z2 = fliplr(z2);     x2 = fliplr(x2);  end
    if x1(1)>x1(end);       x1 = -x1;   end                     % distance positive towards land
    if x2(1)>x2(end);       x2 = -x2;   end                     % distance positive towards land
    
    % Projecteer profielen op dezelfde punten en voeg kruispunt toe als punt (nodig voor gemeten profielen met verschillend grid)
    [xcross, zcross, x1, z1, ~, z2, dircross] = findCrossings(x1, z1, x2, z2, 'synchronizegrids');
    
    % Verwijder nans 
    ind      = ~isnan(z1)&~isnan(z2);
    z1_all   = z1(ind);        z2_all = z2(ind);     x_all = x1(ind);
    
    if ~exist('flag', 'var'); flag = 0; end
    if flag
        % Bepaal kruispunt van einde erosiezone in 1e duinrij (1e crossing boven zs0 waarbij profiel 2 
        % weer gelijk (crossdir = 1) of groter aan profiel 1 (crossdir = 2) is).
        ind1 = find((zcross>zs0 & dircross' > 0));
        ind = find(x_all == xcross(ind1)); 
        clear ind1
        % Gebruik deze index om punten landwaarts van dit punt niet mee te nemen in het erosievolume.
    else
        ind = length(z1_all);   % Index van meest landwaartse punt; alle punten meenemen in erosievolume.
    end
    z1       = z1_all(1:ind);  z2 = z2_all(1:ind);   x = x_all(1:ind);
    
    % Verwijder punten waar sedimentatie plaatsvindt
    ind      = z2 <= z1;  
    x        = x(ind);         z1 = z1(ind);         z2 = z2(ind);
    
    % Verwijder punten onder stormvloedpeil (svp)
    [x_svp1, ~, x1, z1, ~, ~, ~] = findCrossings(x, z1, [min(x) max(x)], [zs0 zs0], 'keeporiginalgrid');    % Voeg punt toe op kruising met zs0
    [x_svp2, ~, x2, z2, ~, ~, ~] = findCrossings(x, z2, [min(x) max(x)], [zs0 zs0], 'keeporiginalgrid');    % Voeg punt toe op kruising met zs0
    
    try
    z1_ero = z1(x1 >= x_svp1(1));     x1_ero = x1(x1 >= x_svp1(1));
    catch
        z1_ero = z1;     x1_ero = x1;
    end
    try 
        z1_ero = z1_ero(x1_ero <= x_svp1(2));     x1_ero = x1_ero(x1_ero <= x_svp1(2));
    end
    try 
        z2_ero = z2(x2 >= x_svp2(1)); 
        x2_ero = x2(x2 >= x_svp2(1));
        try 
            z2_ero = z2_ero(x2_ero <= x_svp2(2));     x2_ero = x2_ero(x2_ero <= x_svp2(2));
        end
    catch
        if sum(z2>zs0) > 1
            z2_ero = z2;
            x2_ero = x1_ero;
        else
            z2_ero = ones(length(x1_ero),1)*zs0;
            x2_ero = x1_ero;
        end
    end
    
        % Bereken erosievolume
    dVero    = polyarea([x1_ero; flipud(x2_ero)],[z1_ero; flipud(z2_ero)]); % m3/m
 
end