OPT = struct(...
    'MSLevel',0,...
    'wavebase',-7,...
    'Tide',0.37,...
    'theta',1,...
    'sealevelchangerate',0.64,...
    'alpha',1,...
    'Beta',0.05,...
    'T',100,...
    'add_sediment',0,...
    'param',1,...
    'modeltrain_number',4,...
    'threshold_factor',2);


OPT = struct(...
    'MSLevel',0,...
    'wavebase',-7,...
    'Tide',0.37,...
    'theta',1,...
    'sealevelchangerate',0.64,...
    'alpha',1,...
    'Beta',0.05,...
    'T',100,...
    'add_sediment',0,...
    'param',1,...
    'modeltrain_number',4,...
    'threshold_factor',2,...
    'inputdir',fullfile(fileparts(which(mfilename)),'..','InputData','Grids'),...
    'stormdir',fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')],['Storm_',num2str(1,'%02i')]),...
    'xbeachdir',fullfile(fileparts(which(mfilename)),'..','InputData','SingleStorm'),...
    'xbeachoutputdir',fullfile(fileparts(which(mfilename)),'..','Output'),...
    'xbeachrundir',fullfile(fileparts(which(mfilename)),'..','Runs'),...
    'recoverydir',fullfile(fileparts(which(mfilename)),'..','Barrier_Island_Recovery'));

% OPT = setproperty(OPT, varargin{:});

varnames = {'H','zb','zs'};
sealevelchange=OPT.sealevelchangerate/1200;

%% INITIAL CONDITIONS

% x-grid
x = load(fullfile(OPT.inputdir,'x.grd'));
xas = x(2,:);
ngridp = size(x,2);

% initial profile
InitialProfile = load(fullfile(OPT.inputdir,'z.dep'));
mkdir(fullfile(OPT.stormdir));
dlmwrite(fullfile(OPT.stormdir,'InitialProfile.txt'),InitialProfile);
initialprofile=InitialProfile(2,:);
D_low=max(initialprofile);

% initial roll-over distance
addpath(OPT.recoverydir)
[IntersectA IntersectB intersections]=CalcintersectionsArea(initialprofile,initialprofile,OPT.Tide,OPT.wavebase,xas);
XDistance=intersections.A.MSL(1,1);
i=6
    crsection='z_run_05.dep'
    cd(OPT.xbeachdir);
    copyfile(fullfile(OPT.xbeachdir,['params_run_',num2str(i,'%02i'),'.txt']),fullfile(OPT.xbeachdir,'params.txt'));

%   Here I tried to work with both cores, which isn't working so far    
%     ! mpiexec -n 2 MPIXBeach.exe Kammi
    ! xbeach.exe
    XB = xb_getdimensions;
    
    %% OUTPUT XBEACH
    
    % pre-storm profile
    temp = load(crsection);
    A(i,:) = temp(2,:);
    
    % post-storm profile
    temp = xb_readvar('zb.dat',XB);
    temp = squeeze(temp(:,:,end))';
    B(i,:) = temp(1,:);
    
    % MSL
    MSL=(ones(1,(size(B(i,:),2)))*OPT.MSLevel)';
    
    % create .dat files for each event
    for j = 1:length(varnames)
        varn{j} = [varnames{j},'_run_',num2str(i,'%02i')];
        copyfile([varnames{j},'.dat'],[varn{j},'.dat']);
    end
    
    %% RECOVERY
    [IntersectA IntersectB intersections]=CalcintersectionsArea(A(i,:),B(i,:),OPT.MSLevel,OPT.wavebase,xas);
    
    % calculate overwash volume
    % highest intersection point between pre- and post-storm profiles
    [x y]=max(intersections.A.B(2,:));
    x1=intersections.A.B(1,y);
    x2=xas(end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    if Volume.accretion<0
        OverwashV(i)=0;
    end
    if Volume.accretion>=0
        OverwashV(i)=Volume.accretion;
    end
    x1=0;
    x2=xas(end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    XBeach.erosion=Volume.erosion;
    XBeach.accretion=Volume.accretion;
    XBeach.continuity=Volume.total;
    
    Volumetranslation(i)=sum(OverwashV);
    
    %Create new equilibriumprofile
    [IntersectA IntersectB intersections]=CalcintersectionsArea(A(i,:),B(i,:),OPT.MSLevel,OPT.wavebase,xas);

    [equilibriumprofile stormprofile]=CreateEquilibriumProfile4(i,OverwashV,initialprofile,A,B,OPT.wavebase,OPT.MSLevel,xas);
    allequilibriumprofiles(i,:)=equilibriumprofile;
    
    % recovery
    recoveryTime=100;
    
    if recoveryTime<=0
        recoveryTime=1;
    end
    
    [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl MaxHeight]=recoverprof1910(OPT.alpha,recoveryTime,stormprofile,equilibriumprofile,OPT.MSLevel,OPT.wavebase,OPT.Tide,OPT.theta,OPT.T,OPT.Beta,xas,sealevelchange,OPT.add_sediment);
    
    %% STORE THE RESULTS
    
    zbrun=['zb_run_',num2str(i,'%02i')];
    zb=xb_readvar([zbrun,'.dat'],XB);
    zsrun=['zs_run_',num2str(i,'%02i')];
    zs=xb_readvar([zsrun,'.dat'],XB);
    Hrun=['H_run_',num2str(i,'%02i')];
    H=xb_readvar([Hrun,'.dat'],XB);
        
    Hurricane.pre_storm=A(i,:);
    Hurricane.post_storm=B(i,:);
    Hurricane.XBeach.Volumes=XBeach;
    Hurricane.XBeach.hydro.zb=zb;
    Hurricane.XBeach.hydro.zs=zs;
    Hurricane.XBeach.hydro.H=H;
    Hurricane.recovery.profiles=recoveryprofileall;
    Hurricane.recovery.continuity=Continuity.total;
    Hurricane.recovery.islandheight=HeightIsland;
    Hurricane.recovery.islandwidth=widthIsland;
    Hurricane.recovery.islandwvolume=VolIsland;
    Hurricane.recovery.islandroll_over=distance-XDistance;
    Hurricane.recovery.islandheight_max=MaxHeight;
    Hurricane.equilibriumprofile=equilibriumprofile;
    Analyses{1,i}=Hurricane;
    
    RunProperties=Analyses{1,i};
    save('RunProperties','RunProperties');