function [modeltrain_number storms] = modeltrainHS(varargin)

modeltraindir = cd;
%% THE MODELTRAIN INPUT PARAMETERS

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',11,...
    'threshold_factor',3,...
    'ngd',1,...
    'nd',20,...
    'dzg',0.1,...
    'sedcal1',1,...
    'sedcal2',0.1);


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

OPT.inputdir = fullfile(fileparts(which(mfilename)),'..','InputData','Grids');
OPT.stormdir = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')],['Storm_',num2str(1,'%02i')]);
OPT.xbeachdir = fullfile(fileparts(which(mfilename)),'..','XBeach');
OPT.xbeachoutputdir = fullfile(fileparts(which(mfilename)),'..','Output');
OPT.xbeachrundir = fullfile(fileparts(which(mfilename)),'..','Runs');
OPT.recoverydir = fullfile(fileparts(which(mfilename)),'..','Barrier_Island_Recovery');

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);

%% PARAM.txt

% dep-file
i=1;
copyfile(fullfile(OPT.inputdir,'z.dep'),fullfile(OPT.stormdir,'z_run_00.dep'));
crsection='z_run_00.dep';

% make the first param.txt
boundaryConditions = fullfile(['jonswaplist_storm_',num2str(i,'%02i'),'.txt']);
surgeConditions=fullfile(['Scaled_Surge_Storm_',num2str(i,'%02i'),'.dat']);
OPT.ngd,OPT.nd,OPT.dzg,OPT.sedcal1,OPT.sedcal2

[params_run]=CreateParam(i,crsection,ngridp,boundaryConditions,surgeConditions,OPT.stormdir,OPT.ngd,OPT.nd,OPT.dzg,OPT.sedcal1,OPT.sedcal2);
% [params_run]=CreateParam(1,crsection,ngridp,boundaryConditions,surgeCondi
% tions,OPT.stormdir);

%% SELECT THE EVENTS
[storms recovery]=selectStorms('modeltrain_number',OPT.modeltrain_number,'threshold_factor',OPT.threshold_factor);
recoveryTimes=[(recovery.time(:)'/30)];
i=1;
while i <= size(storms,2) & D_low>OPT.Tide
    
    OPT.xbeachrunstormdir = fullfile(OPT.xbeachrundir,['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
    OPT.xbeachoutputstormdir = fullfile(OPT.xbeachoutputdir,['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
    mkdir(OPT.xbeachrunstormdir);
    mkdir(OPT.xbeachoutputstormdir);
    
    copyfile(fullfile(OPT.stormdir,['params_run_',num2str(i,'%02i'),'.txt']),fullfile(OPT.xbeachdir,'params.txt'));

    
    % derive the threshold values between the regimes and note the
    % event regime
    Max_WVHT_storm = max(storms{1,i}.Katrinas.waveheight);
    Slope = 0.0035;
    Wavesteepness = 0.016;
    Xi = Slope/sqrt(Wavesteepness);
    
    Hswash = ((D_low - max(storms{i}.Katrinas.surge.one))/(0.83*Xi+0.2));
%     Hswash=(4*(D_low-0.45))/((4*0.83*(0.0035/sqrt(0.016))+0.8)+0.5);
    Hinundation = (D_low -max(storms{i}.Katrinas.surge.one))/(-0.02*Xi+0.14);
%     Hinundation=(0.56*(D_low-0.45)+3.72)/(0.56*(-0.02*(0.0035/sqrt(0.016))+0.14)+1);
    if Max_WVHT_storm<=Hswash
        Hurricane.regime='Swash';
    else if Max_WVHT_storm<=Hinundation
            Hurricane.regime='Overwash';
        else
            Hurricane.regime='Inundation';
        end
    end
    
      
    %% RUN XBEACH
    copyfile(OPT.stormdir,OPT.xbeachdir,'f');
    cd(OPT.xbeachdir);
    
%   Here I tried to work with both cores, which isn't working so far    
%     ! mpiexec -n 2 MPIXBeach.exe Kammi
%     ! xbeach.exe
    ! mpiexec -n 8 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=ceil(recoveryTimes(i));
    
    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;
    
    cd(OPT.xbeachoutputstormdir);
    RunProperties=Analyses{1,i};
    save('RunProperties','RunProperties');
    cd(modeltraindir);
    
    %% PREPARE THE NEXT RUN
    crsection=['z_run_',num2str(i,'%02i'),'.dep'];
    D_low=max(recoveryprofilet);

    
    i=i+1;
    OPT.stormdir = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
    if i <= size(storms,2)
        % make the param_run_i.txt
        boundaryConditions = fullfile(['jonswaplist_storm_',num2str(i,'%02i'),'.txt']);
        surgeConditions=fullfile(['Scaled_Surge_Storm_',num2str(i,'%02i'),'.dat']);
%         [params_run]=CreateParam(i,crsection,ngridp,boundaryConditions,su
%         rgeConditions,OPT.stormdir);
        [params_run]=CreateParam(i,crsection,ngridp,boundaryConditions,surgeConditions,OPT.stormdir,OPT.ngd,OPT.nd,OPT.dzg,OPT.sedcal1,OPT.sedcal2);
        % write the next pre-storm profile
        dlmwrite(fullfile(OPT.stormdir,crsection),[recoveryprofilet;recoveryprofilet;recoveryprofilet],'delimiter',' ','precision','%.10f');
        [value kolom] = max(recoveryprofilet);
        make_gdist_input(ngridp,3,OPT.nd,OPT.ngd,kolom,OPT.stormdir);
        

    end
    
    %% CLEANING
    cd(OPT.xbeachdir);
    try
        delete('Storm*.inp');
        movefile('*.dat',OPT.xbeachrunstormdir);
        movefile('*.bcf',OPT.xbeachrunstormdir);
        movefile('*Storm*.inp',OPT.xbeachrunstormdir);
        movefile('*.txt',OPT.xbeachrunstormdir);
        movefile('*.dep',OPT.xbeachrunstormdir);
    catch
        delete('Storm*.inp');
        movefile('*.txt',OPT.xbeachrunstormdir);
        movefile('*.dep',OPT.xbeachrunstormdir);
    end
    cd(modeltraindir);
end
OPT.xbeachoutputtraindir = fullfile(OPT.xbeachoutputdir,['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')]);
cd(OPT.xbeachoutputtraindir)
save('Analyses','Analyses')
cd(modeltraindir);
