function Main_script_RivSed
% =====================================================================
% Model setting
% =====================================================================
% parameter setting
v_s     = 4E-6; % fall velocity [m/s]
rho_b   = 2650; % bank material density [kg/m3]
tau_e	= 2.5; % particle erosion threshold [N/m2]
tau_m	= 3.5; % mass erosion theshold [N/m2]
m_pe	= 30; % Particle erosion rate [kg/m/d]
m_me	= 40; % Particle erosion rate [kg/m/d]
kappa   = 0.0018; % erodibility coefficient [m3/N/d]
tau_b   = 5; % shear stress threshold for bank erosion [N/m2]
% parameters for transferring to function Sediment_Koch,Sediment_Kas,Sediment_Amm
Para_trans = [v_s,rho_b,tau_e,tau_m,m_pe,m_me,kappa,tau_b];

% set time
dt_Hec      = 1/24; % time interval in Hec-Ras [d]
t0          = datenum(2013,1,1,0,0,0); % starting time
t_span      = 365*4; % [day]
% te          = t0+t_span; % end time
dt_out      = 5/60/24;    % time step of output [d]
n_out_h     = 12; % how many output steps per hour
% no_data_Hec = t_span/dt_Hec+1; % number of data set
n_sep_d     = 120; % seperations for a model run [d]

% number of compounds and coresponging age
n_mob      = 1;	  % mobile *** in this module, n_im must be equal to n_mob. IN THIS CASE: TSS
n_im       = 1;	  % immobile    IN THIS CASE: SEDIMENT ON THE BED 
n_mob_age  = 1;	  % AGE of TSS, not compute in this case
n_im_age   = 1;   % AGE of SEDIMENT, not compute in this case
n_dep      = 1;   % deposition
n_bde      = 1;   % bed erosion
n_bke      = 1;   % bank erosion
n_all      = 7;	  % all components 
% model components for transferring to function Sediment_Koch,Sediment_Kas,Sediment_Amm
mod_comp   =[n_mob,n_im,n_mob_age,n_im_age,n_dep,n_bde,n_bke,n_all];

% =====================================================================
% Model computation
% =====================================================================

for i_tset=1:ceil(t_span/n_sep_d)
    % time setting for transferring to function Sediment_Koch,Sediment_Kas,Sediment_Amm
    % every 120 days as a simulation period to avoid matrix exceeding limit
    if i_tset==1
        t_begin = t0;
    else
        t_begin = t0+(i_tset-1)*n_sep_d;
    end
    
    if i_tset < ceil(t_span/n_sep_d) 
        t_end = t_begin+n_sep_d;
    else
        t_end = t0+t_span;
    end
    
    t_sim   = t_end-t_begin;
    time_set=[dt_Hec,t_begin,t_sim,t_end,dt_out,i_tset,n_out_h];
       
    % =====================================================================
    % parallel computation for Kochart and Kasbach
    % =====================================================================
    disp([datestr(now) ': start parallel computing for Koch and Kas']);
    
    funList = {@Sediment_Koch_m,@Sediment_Kas_m};
    parfor i=1:length(funList)
        % call the function
        funList{i}(Para_trans,time_set,mod_comp);
    end
    disp([datestr(now) ': finish Kas and Koch, start Amm']);
    
    % run sediment for Ammer stem
    Sediment_Amm_m(Para_trans,time_set,mod_comp);
    disp([datestr(now) ': finish sediment simulation'])   
end
end