function AnalyseXBeachOutput_Hydro(info,dirs)
    test_name = info.run;
    
    %% Settings
    refdate = datetime(2017,2,1);
    refdate0 = datetime(2017,2,3,3,0,0);
    tstop = 255000;
    endref = refdate + seconds(tstop);
    
    %% Load measurement data
    load([dirs.data filesep 'obs_hydro.mat']);
%     pos_meas = load([basepath '\scripts\bathy.mat']);
    for i_obs=1:numel(obs_points)
        obs_time.(obs_points{i_obs}) = datenum(data.(obs_points{i_obs})(:,1:6));
    end


    %% Load xb data
    fname = [dirs.output filesep 'xboutput.nc'];
    pointtime   = nc_varget(fname,'pointtime');
    point_zs    = nc_varget(fname,'point_zs');
    point_zb    = nc_varget(fname,'point_zb');
    point_H     = nc_varget(fname,'point_H');
    pointx      = nc_varget(fname,'pointx');
    pointy      = nc_varget(fname,'pointy');

    meantime    = nc_varget(fname,'meantime');
    globalx     = nc_varget(fname,'globalx');
    H_mean      = nc_varget(fname,'H_mean'); % Hrms short waves
    H_min      = nc_varget(fname,'H_min'); % Hrms short waves
    H_max      = nc_varget(fname,'H_max'); % Hrms short waves
    zs1_var     = nc_varget(fname,'zs1_var'); % water level without tide
    zs_mean     = nc_varget(fname,'zs_mean'); % mean water level 
    zb_mean     = nc_varget(fname,'zb_mean'); % mean zb

    id_peak = 110; % of 110 | 3 february 7:50 | refdate+seconds(meantime(id_peak))+seconds(3600)
    Hm0_HF_peak = squeeze(H_mean(id_peak,1,:) * sqrt(2));
    Hm0_min_HF_peak = squeeze(H_min(id_peak,1,:) * sqrt(2));
    Hm0_max_HF_peak = squeeze(H_max(id_peak,1,:) * sqrt(2));
    Hm0_LF_peak = squeeze(4*sqrt(zs1_var(id_peak,1,:)));
    zs_mean_peak = squeeze(zs_mean(id_peak,1,:));
    zb_mean_peak = squeeze(zb_mean(1,1,:));
    
    %% Process model data
    ts = 513000:18000:tstop*10; %51300 corresponds to first measurement
    %ts2 = tstop*10-(20*36000):18000:tstop*10; %51300 corresponds to first measurement
    shift = 11990;
    for ii=1:size(point_zs,2)
        for jj = 1:length(ts)
            try
                zstmp = point_zs(ts(jj)-shift:ts(jj+1)-shift,ii);
                Hrms_hf(ii,jj) = rms(point_H(ts(jj)-shift:ts(jj+1)-shift,ii));
                mean_wl(ii,jj) = mean(point_zs(ts(jj)-shift:ts(jj+1)-shift,ii));
                mean_zb(ii,jj) = mean(point_zb(ts(jj)-shift:ts(jj+1)-shift,ii));

                y_IG            = lfhf2(detrend(point_zs(ts(jj)-shift:ts(jj+1)-shift,ii)),0,0.05,0.1);
                [freqS,PwrS] = spectrumsimple(1:length(y_IG),y_IG,'figures',0);
                Hm0_IG(ii,jj)      = spec2hs(freqS,PwrS);
            catch
                Hrms_hf(ii,jj) = NaN;
                mean_wl(ii,jj) = NaN;
                mean_zb(ii,jj) = NaN;
                Hm0_IG(ii,jj)= NaN;
            end
        end
    end

    %% Compare model results with measurement data
    % 1 hour shift between observed model data en measurement data, hence seconds(3600)
    % added to model data
    ts=ts/10;
    for j=1:numel(obs_points)                     
        if j==1 %only plot entire time series for offshore ADCP        
            idStart = 1;
            idEnd = length(obs_time.(obs_points{2})); 
            plottime_start = datetime(obs_time.(obs_points{2})(idStart), 'ConvertFrom', 'datenum');
            plottime_end = datetime(obs_time.(obs_points{2})(idEnd), 'ConvertFrom', 'datenum');
            id = find(refdate+seconds(ts)+seconds(3600) >= plottime_end, 1, 'first');

            figure(j);
            set(gcf, 'Position', [200 200 900 500])
            subplot(4,1,1)
            hold on
            a = plot(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'), data.(obs_points{j})(:,7), 'ko', 'markersize', 3);
            b = plot(refdate+seconds(ts(1:id))+seconds(3600), mean_wl(j,1:id)-mean_zb(j,1:id), 'linewidth', 1.5);
            plot([plottime_start plottime_start], [7 14], 'k--')
            plot([plottime_end plottime_end], [7 14], 'k--')
            ylabel('Water depth [m]')
            set(gca,'fontsize',11,'XTicklabel',[]);
            legend([a,b], {'data', 'model'}, 'location', 'northeast')
            xlim([min(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum')), max(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'))])
            ylim([7 14]);grid on;

            subplot(4,1,2)
            hold on
            plot(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'), data.(obs_points{j})(:,8), 'ko', 'markersize', 3)
            plot(refdate+seconds(ts(1:id))+seconds(3600), Hrms_hf(j,1:id)*sqrt(2), 'linewidth', 1.5)
            plot([plottime_start plottime_start], [0 7], 'k--')
            plot([plottime_end plottime_end], [0 7], 'k--')
            ylabel('H_{m0,HF} [m]')
            set(gca,'fontsize',11,'YTick',[0:2:8],'YTicklabel',[0:2:8],'XTicklabel',[]);
            xlim([min(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum')), max(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'))])
            ylim([0 7]); grid on;

            subplot(4,1,3)
            hold on
            plot(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'), data.(obs_points{j})(:,14), 'ko', 'markersize', 3)
            plot(refdate+seconds(ts(1:id))+seconds(3600), Hm0_IG(j,1:id), 'linewidth', 1.5)
            plot([plottime_start plottime_start], [0 1.5], 'k--')
            plot([plottime_end plottime_end], [0 1.5], 'k--')
            ylabel('H_{m0,LF} [m]')
            set(gca,'fontsize',11,'XTicklabel',[]);     
            xlim([min(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum')), max(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'))])
            ylim([0 1.5]); grid on;      

            subplot(4,1,4)
            hold on
            plot(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'), data.(obs_points{j})(:,12), 'ko', 'markersize', 3)
            plot([plottime_start plottime_start], [0 25], 'k--')
            plot([plottime_end plottime_end], [0 25], 'k--')
            ylabel('Wave periods [s]')
            xlim([min(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum')), max(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum'))])
            set(gca,'fontsize',11,'ylim',[0 25]);
            xlabel('Date')
            grid on;
            
            rms_skill(j,1:3) = NaN; sci(j,1:3) = NaN; bias(j,1:3) = NaN; relbias(j,1:3) = NaN; r2(j,1:3) = NaN; bss(j,1:3) = NaN;

        else
            idStart = find(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum') <= refdate, 1, 'last');
            if isempty(idStart); idStart = 1; end
            idEnd = find(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum') >= endref, 1, 'first');
            if isempty(idEnd); idEnd = length(obs_time.(obs_points{j})); end

            id = ~isnan(data.(obs_points{j})(idStart:idEnd-1,8));
            id(1:2)=0;

            timeshift = 3599;
            if (j == 2) || (j == 3) || (j == 5) || (j == 6)|| (j == 8)
                id2 = 3;
             elseif j == 4 
                id2 = 5;
            else
                id2 = 1;
            end

            figure(j)
            set(gcf, 'Position', [200 200 900 500])
            subplot(3,1,1)
            hold on
            plot(datetime(obs_time.(obs_points{j})(idStart:idEnd-1), 'ConvertFrom', 'datenum'), data.(obs_points{j})(idStart:idEnd-1,7), 'ko', 'markersize', 5)
            plot(refdate+seconds(ts)+seconds(timeshift), mean_wl(j,:)-mean_zb(j,:), 'linewidth', 1.5)

            d_time = obs_time.(obs_points{j})(idStart+id2:idEnd); 
            d_data = data.(obs_points{j})(idStart+id2:idEnd,7);
            m_time = datenum(refdate+seconds(ts(id(id2:end)))+seconds(timeshift))';
            m_data = (mean_wl(j,(id(id2:end)))-mean_zb(j,(id(id2:end))))';
%             plot(d_time,d_data,'r.');hold on; plot(m_time,m_data,'.');

            try
            [r2(j,1) sci(j,1) relbias(j,1) bss(j,1) rms_skill(j,1) bias(j,1)] = xb_skill_BOI([d_time,d_data], [m_time,m_data]);
            text(0.15,0.9,['RMSE = ' num2str(round(rms_skill(j,1),2))], 'Units', 'normalized');           
            text(0.15,0.8,['sci = ' num2str(round(sci(j,1),2))], 'Units', 'normalized'); 
            text(0.15,0.7,['bias = ' num2str(round(bias(j,1),2))], 'Units', 'normalized'); 
            text(0.15,0.6,['relbias = ' num2str(round(relbias(j,1),2))], 'Units', 'normalized'); 
%             clear r2 sci relbias bss rms bias
            catch
                rms_skill(j,1) = NaN; sci(j,1) = NaN; bias(j,1) = NaN; relbias(j,1) = NaN; r2(j,1) = NaN; bss(j,1) = NaN;
            end

            xlim([(datetime(obs_time.(obs_points{j})(idStart), 'ConvertFrom', 'datenum')), (datetime(obs_time.(obs_points{j})(idEnd), 'ConvertFrom', 'datenum'))])
            ylabel('Water depth [m]')
            set(gca,'fontsize',11,'ylim',[0 3],'XTicklabel',[]);
            grid on;

            subplot(3,1,2)
            hold on
            plot(datetime(obs_time.(obs_points{j})(idStart:idEnd-1), 'ConvertFrom', 'datenum'), data.(obs_points{j})(idStart:idEnd-1,8), 'ko', 'markersize', 5)
            plot(refdate+seconds(ts)+seconds(timeshift), Hrms_hf(j,:)*sqrt(2), 'linewidth', 1.5)  

            % skill scores
            d_time = obs_time.(obs_points{j})(idStart+id2:idEnd); 
            d_data = data.(obs_points{j})(idStart+id2:idEnd,8);
            m_time = datenum(refdate+seconds(ts(id(id2:end)))+seconds(timeshift))';
            m_data = (Hrms_hf(j,id(id2:end))*sqrt(2))';
%             plot(d_time,d_data,'r.');hold on; plot(m_time,m_data,'.');

            try
            [r2(j,2) sci(j,2) relbias(j,2) bss(j,2) rms_skill(j,2) bias(j,2)] = xb_skill_BOI([d_time,d_data], [m_time,m_data]);
            text(0.15,0.9,['RMSE = ' num2str(round(rms_skill(j,2),2))], 'Units', 'normalized');           
            text(0.15,0.8,['sci = ' num2str(round(sci(j,2),2))], 'Units', 'normalized'); 
            text(0.15,0.7,['bias = ' num2str(round(bias(j,2),2))], 'Units', 'normalized'); 
            text(0.15,0.6,['relbias = ' num2str(round(relbias(j,2),2))], 'Units', 'normalized'); 
%             clear r2 sci relbias bss rms bias
            catch
                rms_skill(j,2) = NaN; sci(j,2) = NaN; bias(j,2) = NaN; relbias(j,2) = NaN; r2(j,2) = NaN; bss(j,2) = NaN;
            end

            xlim([(datetime(obs_time.(obs_points{j})(idStart), 'ConvertFrom', 'datenum')), (datetime(obs_time.(obs_points{j})(idEnd), 'ConvertFrom', 'datenum'))])
            ylabel('H_{m0,HF} [m]')
            set(gca,'fontsize',11,'ylim',[0 2],'XTicklabel',[]);
            grid on;

            subplot(3,1,3)
            hold on
            a = plot(datetime(obs_time.(obs_points{j})(idStart:idEnd), 'ConvertFrom', 'datenum'), data.(obs_points{j})(idStart:idEnd,14), 'ko', 'markersize', 5);
            b = plot(refdate+seconds(ts)+seconds(timeshift), Hm0_IG(j,:), 'linewidth', 1.5);
            ylabel('H_{m0,LF} [m]')
            xlabel('Date')
            legend([a,b], {'data', 'model'}, 'location', 'northeast')           
            xlim([(datetime(obs_time.(obs_points{j})(idStart), 'ConvertFrom', 'datenum')), (datetime(obs_time.(obs_points{j})(idEnd), 'ConvertFrom', 'datenum'))])
            set(gca,'fontsize',11,'ylim',[0 2]);

            d_time = obs_time.(obs_points{j})(idStart+id2:idEnd); 
            d_data = data.(obs_points{j})(idStart+id2:idEnd,14);
            m_time = datenum(refdate+seconds(ts(id(id2:end)))+seconds(timeshift))';
            m_data = (Hm0_IG(j,id(id2:end)))';
%             plot(d_time,d_data,'r.');hold on; plot(m_time,m_data,'.');

            try
            [r2(j,3) sci(j,3) relbias(j,3) bss(j,3) rms_skill(j,3) bias(j,3)] = xb_skill_BOI([d_time,d_data], [m_time,m_data]);
            text(0.15,0.9,['RMSE = ' num2str(round(rms_skill(j,3),2))], 'Units', 'normalized');           
            text(0.15,0.8,['sci = ' num2str(round(sci(j,3),2))], 'Units', 'normalized'); 
            text(0.15,0.7,['bias = ' num2str(round(bias(j,3),2))], 'Units', 'normalized'); 
            text(0.15,0.6,['relbias = ' num2str(round(relbias(j,3),2))], 'Units', 'normalized'); 
%             clear r2 sci relbias bss rms bias
            catch
                rms_skill(j,3) = NaN; sci(j,3) = NaN; bias(j,3) = NaN; relbias(j,3) = NaN; r2(j,3) = NaN; bss(j,3) = NaN;
            end
            grid on;
        end
        xb_write_plot(gcf, dirs.output, ['Hydro_at_' obs_points{j}])
    end
    save([dirs.output filesep 'skill.mat'], 'obs_points','rms_skill', 'sci','bias', 'relbias')
    
    dlm = {'$WL$\newline$_{RMSE}$','$WL$\newline$_{SCI}$','$WL$\newline$_{bias}$','$WL$\newline$_{relbias}$','$HF$\newline$_{RMSE}$','$HF$\newline$_{SCI}$',...
    '$HF$\newline$_{bias}$','$HF$\newline$_{relbias}$','$IG$\newline$_{RMSE}$','$IG$\newline$_{SCI}$','$IG$\newline$_{bias}$','$IG$\newline$_{relbias}$'}; 
    matrix2latex([rms_skill(2:end-1,1),sci(2:end-1,1),bias(2:end-1,1),relbias(2:end-1,1),rms_skill(2:end-1,2),sci(2:end-1,2),bias(2:end-1,2),relbias(2:end-1,2),rms_skill(2:end-1,3),sci(2:end-1,3),bias(2:end-1,3),relbias(2:end-1,3)], 'filename', 'skill.tex', ...
    'rowjustification','p{0.75cm}',...
    'rowlabel', obs_points(2:end-1), ...
    'collabel', dlm, ...
    'format',	{'%s', '%4.2f', '%4.2f','%4.2f','%4.2f','%4.2f', '%4.2f','%4.2f','%4.2f','%4.2f', '%4.2f','%4.2f','%4.2f'});

    
    close all    
    %% Scatter new

    WL_total_ddata = [];
    WL_total_mdata = [];
    HF_total_ddata = [];
    HF_total_mdata = [];
    IG_total_ddata = [];
    IG_total_mdata = [];

    CM = jet(9);
    figure(11);
    set(gcf, 'Position', [200 200 1000 300])

    for j=2:numel(obs_points)-1 
            idStart = find(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum') <= refdate, 1, 'last');
            if isempty(idStart); idStart = 1; end
            idEnd = find(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum') >= endref, 1, 'first');
            if isempty(idEnd); idEnd = length(obs_time.(obs_points{j})); end

            id = ~isnan(data.(obs_points{j})(idStart:idEnd-1,8));
            id(1:2)=0;
            timeshift = 3599;
            if (j == 2) || (j == 3) || (j == 5) || (j == 6)|| (j == 8)
                id2 = 3;
             elseif j == 4 
                id2 = 5;
            else
                id2 = 1;
            end

            d_time = obs_time.(obs_points{j})(id);

            subplot(1,4,1) 
            hold on 
            d_data = data.(obs_points{j})(id,7);
            m_data = (mean_wl(j,(id(id2:end)))-mean_zb(j,(id(id2:end))))';
            WL_total_ddata = [WL_total_ddata; d_data];
            WL_total_mdata = [WL_total_mdata; m_data];
            d = scatter(d_data,m_data, 5, CM(j,:), 'filled');


            subplot(1,4,2) 
            hold on
            d_data = data.(obs_points{j})(id,8);
            m_data = (Hrms_hf(j,id(id2:end))*sqrt(2))';
            HF_total_ddata = [HF_total_ddata; d_data];
            HF_total_mdata = [HF_total_mdata; m_data];
            d = scatter(d_data,m_data, 5, CM(j,:), 'filled');  
            plot([0 2],[0 2],'k:','linewidth',1.5);


            subplot(1,4,3); 
            hold on 
            d_data = data.(obs_points{j})(id,14);
            m_data = (Hm0_IG(j,id(id2:end)))';
            IG_total_ddata = [IG_total_ddata; d_data];
            IG_total_mdata = [IG_total_mdata; m_data];
            d = scatter(d_data,m_data, 5, CM(j,:), 'filled');


            subplot(1,4,4); hold on; 
            d_data = data.(obs_points{j})(id,14);
            m_data = (Hm0_IG(j,id(id2:end)))';
            d = scatter(d_data,m_data, 5, CM(j,:), 'filled');
            xlim([100 101]); ylim([100 101])

    end

    % best fits
    WL_p = polyfit(WL_total_mdata,WL_total_ddata,1); 
    WL_f = polyval(WL_p,WL_total_ddata); 
    HF_p = polyfit(HF_total_mdata,HF_total_ddata,1); 
    HF_f = polyval(HF_p,HF_total_ddata); 
    IG_p = polyfit(IG_total_mdata,IG_total_ddata,1); 
    IG_f = polyval(IG_p,IG_total_ddata); 

%     try
%     % skillscores
%     [WL_r2 WL_sci WL_relbias WL_bss WL_rms WL_bias] = xb_skill_BOI([[1:1:length(WL_total_ddata)]; WL_total_ddata'],[[1:1:length(WL_total_mdata)]; WL_total_mdata']); 
%     [HF_r2 HF_sci HF_relbias HF_bss HF_rms HF_bias] = xb_skill_BOI([[1:1:length(WL_total_ddata)]; HF_total_ddata'],[[1:1:length(WL_total_mdata)]; HF_total_mdata']); 
%     [IG_r2 IG_sci IG_relbias IG_bss IG_rms IG_bias] = xb_skill_BOI([[1:1:length(WL_total_ddata)]; IG_total_ddata'],[[1:1:length(WL_total_mdata)]; IG_total_mdata']); 
%     catch
%     end
    
    figure(11); hold on; 
    subplot(1,4,1);hold on;  axis equal; box on; grid on;
    % plot(WL_total_ddata,total_mdata,'.');hold on;
    k2 = plot(WL_f,WL_total_ddata,'k','linewidth',1.5);
    k1 = plot([0 3],[0 3],'k:','linewidth',1.5);
    set(gca,'xtick',0:0.5:3.0); set(gca,'ytick',0:0.5:3.0); 
    set(gca,'fontsize',9)
    xlim([0 3]); ylim([0 3]);
    xlabel('Measured water depth [m]', 'interpreter', 'tex')
    ylabel('Modeled water depth [m]', 'interpreter', 'none')
    % legend([d(2:9)';k1;k2],'PT1','PT2','PT3','PT4','PT5','PT6','PT7','PT8','best fit','linear fit','Position',[0.75 0.47 0.1 0.1]);
    subplot(1,4,2);hold on;  axis equal; box on; grid on;
    % plot(HF_total_ddata,total_mdata,'.');hold on;
    plot(HF_f,HF_total_ddata,'k','linewidth',1.5);
    xlim([0 2]); ylim([0 2]);
    xlabel('Measured H_{m0,G} [m]', 'interpreter', 'tex')
    ylabel('Modeled H_{m0,G} [m]', 'interpreter', 'tex')
    set(gca,'xtick',0:0.5:2.0); set(gca,'ytick',0:0.5:2.0); 
    set(gca,'fontsize',9)
    subplot(1,4,3);hold on;  axis equal; box on; grid on;
    % plot(IG_total_ddata,total_mdata,'.');hold on;
    plot([0 2],[0 2],'k:','linewidth',1.5);
    plot(IG_f,IG_total_ddata,'k','linewidth',1.5);
    xlabel('Measured H_{m0,IG} [m]', 'interpreter', 'tex')
    ylabel('Modeled H_{m0,IG} [m]', 'interpreter', 'tex')
    xlim([0 2]); ylim([0 2]);
    set(gca,'xtick',0:0.5:2.0); set(gca,'ytick',0:0.5:2.0); 
    set(gca,'fontsize',9)    
    subplot(1,4,4);hold on; axis off
    % plot(IG_total_ddata,total_mdata,'.');hold on;
    plot([0 2],[0 2],'k:','linewidth',1.5);
    plot(IG_f,IG_total_ddata,'k','linewidth',1.5);
    legend([obs_points(2:9),'best fit','linear fit'], 'location', 'west')
    
    xb_write_plot(gcf, dirs.output, ['scatter_hydro'])

    %%  Figure of cross-shore waterlevel, HF and LF wave height at peak of the storm
    load([dirs.data filesep 'distgr.mat'])
    % observed
    for j=1:numel(obs_points)-1  
        id_peak_obs(j) = find(datetime(obs_time.(obs_points{j}), 'ConvertFrom', 'datenum') < (refdate+seconds(meantime(id_peak))+seconds(3600)),1,'last');
        datetime(obs_time.(obs_points{j})(id_peak_obs(j),1), 'ConvertFrom', 'datenum')
        data_obs_wl(j) = data.(obs_points{j})(id_peak_obs(j),7)-obs_position(j,3);
        data_obs_HF(j) = data.(obs_points{j})(id_peak_obs(j),8);
        data_obs_LF(j) = data.(obs_points{j})(id_peak_obs(j),14);
    end

    WL = zs_mean_peak;
    id_zero = find((WL-zb_mean_peak) <= 0.4, 1, 'first');

    figure;
    subplot(3,1,1);hold on;
    fill([distgr(1:id_zero)'; distgr(id_zero);distgr(1)],[zb_mean_peak(1:id_zero); min(zb_mean_peak(1:id_zero)); min(zb_mean_peak(1:id_zero))] , [255 200 0]./255, 'FaceAlpha', 0.4, 'EdgeColor', 'none', 'HandleVisibility','off')
    hold on;
    fill([distgr(1:(id_zero))';flipud(distgr(1:(id_zero))');],[zb_mean_peak(1:id_zero);flipud(WL(1:id_zero));], [59 171 227]./255, 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'HandleVisibility','off')
    plot(distgr(index_points(1:9)),data_obs_wl(1:9),'ko','markersize', 5,'MarkerFaceColor','k'); hold on;
    plot(distgr(1:id_zero),zb_mean_peak(1:id_zero),'r','linewidth',1);
    set(gca,'fontsize',10,'xlim',[-200 120],'xtick',[-200:50:120],'xticklabel',[-200:50:120],'ylim',[-2 3]);
    grid on;
    ylabel('Elevation [m+MSL]'); 
    subplot(3,1,2)
    plot(distgr(index_points(1:9)),data_obs_HF(1:9),'ko', 'markersize', 5); hold on;
    plot(distgr(1:id_zero),Hm0_HF_peak(1:id_zero), 'linewidth', 1.5); 
    fill([distgr(1:id_zero)';flipud(distgr(1:id_zero)');],[Hm0_min_HF_peak(1:id_zero);flipud(Hm0_max_HF_peak(1:id_zero));], [0.7 0.7 0.7], 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'HandleVisibility','off')
    set(gca,'fontsize',10,'xlim',[-200 120],'xtick',[-200:50:120],'xticklabel',[-200:50:120]);
    set(gca,'ylim',[0 3]);grid on;
    ylabel('H_{m0,HF} [m]'); 
    subplot(3,1,3)
    plot(distgr(index_points(1:9)),data_obs_LF(1:9),'ko', 'markersize', 5); hold on;
    plot(distgr(1:id_zero),Hm0_LF_peak(1:id_zero), 'linewidth', 1.5);
    set(gca,'xlim',[-200 120],'xtick',[-200:50:120],'xticklabel',[-200:50:120]);
    set(gca,'fontsize',10,'ylim',[0 3]);grid on;
    ylabel('H_{m0,LF} [m]'); xlabel('Distance along transect (m)')
    
    xb_write_plot(gcf, dirs.output, ['cross-shore-at-peak'])
    
end