function [sealevelprofile sealevelprofiletrain seatransl]=SeaLevelRise(seaprof_single,MSLevel,wavebase,sealevelchange,xas,recoveryTime)
% function [sealevelprofile sealevelprofiletrain seatransl stormprofile equilibriumprofile]=SeaLevelRise(seaprof_single,MSLevel,wavebase,sealevelchange,xas,recoveryTime,equilibriumprofile,stormprofile)
% This function takes into account the sealevel rise and creates a new
% profile in such a way that the amount of sediment 'eroded' on the
% foreshore is equal to the amount of accretion of the back side of the
% island. The input parameters needed for this function are listed below;
%   seaprof_single      = pre- sealevelrise profile
%   MSLevel             = Mean Sea Level
%   wavebase            = wavebase/closure depth
%   sealevelchange      = monthly sealevel change rate
%   xas                 = x-points (needed to calculate the area)
%   recoveryTime        = The recovery time after each hurricane
% recoveryTime=1;
% The initial pre-sealevel profile
endvalue=seaprof_single;

% The initial pre-sealevel profile is raised with the same amount as the
% sealevelchange = sealevelrate(century) / 1200 (months) *
% recoveryTime(months)
sealevelprof=seaprof_single+sealevelchange*recoveryTime;

% Find the wavebase intersection, below this intersection the profile
% doesn't change.
[IntersectA IntersectB intersections]=CalcintersectionsArea(seaprof_single,sealevelprof,MSLevel,wavebase,xas);
[rd cd]=find(xas>intersections.A.wavebase(1,1),1);
sealevelprof(1:(cd-1))=sealevelprof(1:(cd-1))-sealevelchange*recoveryTime;

x1=intersections.A.wavebase(1,1);
x2=xas(end);

SeaVol=-100;
j=1;
while SeaVol<0
    
    xtrans=j*0.1;
    xas1=xas(1:cd-1);
    xas2=xas(cd:end)+xtrans;
    xascombi=[xas1(1:cd-1),xas2];
    
    sealevelprof=interp1(xascombi,sealevelprof,xascombi);
    sealevelprofx=interp1(xascombi,sealevelprof,xas);
    [rowf columnf]=find(abs(sealevelprofx-sealevelprofx(end))>0.01);
    sealevelprofx(columnf(end)+1:end)=endvalue(columnf(end)+1:end);
    
    [IntersectA IntersectB intersections]=CalcintersectionsArea(seaprof_single,sealevelprofx,MSLevel,wavebase,xas);
    
    x1=intersections.A.wavebase(1,1)-xtrans;
%         x2=intersections.B.MSL(1,2);
%     x2=intersections.A.MSL(1,2)+xtrans;
    x2=xas(end);

    [rowsealastpart columnsealastpart]=find(xas>=x2,1);
    sealevelprofx(columnsealastpart:end)=seaprof_single(columnsealastpart:end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    SeaVol=Volume.total; %zou 0 moeten zijn.
    

    j=j+1;
end
sealevelprofile=sealevelprofx;
sealevelprofiletrain=sealevelprofx-sealevelchange*recoveryTime;
sealevelprofiletrain(columnf(end)+1:end)=endvalue(columnf(end)+1:end);
seatransl=xtrans;

