function Sediment_Koch_m(Para_trans,time_set,mod_comp)
% parameter for sediment transport
v_s     = Para_trans(1); % fall velocity [m/s]
rho_b   = Para_trans(2); % bank material density [kg/m3]
tau_e	= Para_trans(3); % particle erosion threshold [N/m2]
tau_m	= Para_trans(4); % mass erosion theshold [N/m2]
m_pe	= Para_trans(5); % Particle erosion rate [kg/m/d]
m_me	= Para_trans(6); % Particle erosion rate [kg/m/d]
kappa   = Para_trans(7); % erodibility coefficient [m3/N/d]
tau_b   = Para_trans(8); % shear stress threshold for bank erosion [N/m2]

% set time
dt_Hec      = time_set(1); % time interval in Hec-Ras [d]
t0          = time_set(2); % starting time
t_span      = time_set(3); % [day]
te          = time_set(4); % end time
dt_out      = time_set(5); % time step of output [d]
i_tset      = time_set(6); % number for time sepration
n_out_h     = time_set(7); % output steps per hour
t_Hec       = t0:dt_Hec:te; % time from Hec-Ras [d]
no_data_Hec = t_span/dt_Hec+1; % number of data set

% number of compounds and coresponging age
n_mob      = mod_comp(1);	% mobile *** in this module, n_im must be equal to n_mob. IN THIS CASE: TSS
n_im       = mod_comp(2);	% immobile    IN THIS CASE: SEDIMENT ON THE BED 
n_mob_age  = mod_comp(3);	% AGE of TSS, not compute in this case
n_im_age   = mod_comp(4);   % AGE of SEDIMENT, not compute in this case
n_dep      = mod_comp(5);   % deposition
n_bde      = mod_comp(6);   % bed erosion
n_bke      = mod_comp(7);   % bank erosion
n_all      = mod_comp(8);	% all components 

% data for all of XS from Hec-Ras
hdfdata     = 'Amm_final.p07.hdf';
hdfgeo      = 'Amm_final.g04.hdf';
t_HR        = 735235+double(h5read(hdfdata,'/Results/Sediment/Output Blocks/Sediment/Sediment Time Series/Time')); % Time of Hec-Ras [d]
Q           = double(h5read(hdfdata,'/Results/Sediment/Output Blocks/Sediment/Sediment Time Series/Cross Sections/Flow')); % Flow [m3/s]
u           = double(h5read(hdfdata,'/Results/Sediment/Output Blocks/Sediment/Sediment Time Series/Cross Sections/Velocity')); % Velocity [m/s]
dep_w       = double(h5read(hdfdata,'/Results/Sediment/Output Blocks/Sediment/Sediment Time Series/Cross Sections/Effective Depth'));% Effective Depth [m]
tau         = double(h5read(hdfdata,'/Results/Sediment/Output Blocks/Sediment/Sediment Time Series/Cross Sections/Shear Stress'));% Shear stress [Pa]
L_XS_all	= double(h5read(hdfgeo,'/Geometry/Cross Sections/Lengths'));% cell length [m]
L_XS        = double(L_XS_all(2,:)');
A           = double(Q ./u); % XS area [m2]

% find first time point to extract data from Hec-Ras
[m_HR,~]=find(t_HR==t0);

% input data for river Kochartgraben
XS_koch_sta = 1; % start number of koch from all of XS list
XS_koch_end = 109; % end number of koch from all of XS list

Q_koch      = Q(XS_koch_sta:XS_koch_end,m_HR:no_data_Hec+m_HR-1);
u_koch      = u(XS_koch_sta:XS_koch_end,m_HR:no_data_Hec+m_HR-1);
dep_w_koch	= dep_w(XS_koch_sta:XS_koch_end,m_HR:no_data_Hec+m_HR-1);
tau_koch	= tau(XS_koch_sta:XS_koch_end,m_HR:no_data_Hec+m_HR-1);
L_XS_koch	= L_XS(XS_koch_sta:XS_koch_end,1);% length after each XS
A_koch      = A(XS_koch_sta:XS_koch_end,m_HR:no_data_Hec+m_HR-1);
n_XS_koch	= length(A_koch(:,1)); % number of XS
n_cell_koch	= n_XS_koch-1; % number of cell
L_cell_koch = L_XS_koch(1:end-1,1);% cell length
V_koch      = zeros(n_cell_koch,length(A_koch(1,:))); % cell volume between 2 XS [m3]

% use average XS area of 2 XS * cell length to calculate Volume for a cell
for i_A_koch=1:length(A_koch(1,:))
    V_koch(:,i_A_koch) = 0.5*(A_koch(1:end-1,i_A_koch)+A_koch(2:end,i_A_koch)).* L_XS_koch(1:end-1);
end

V_0=V_koch(:,1); % initial volume [m3]

% lateral flow and conc. from catchment generating model, hourly data
lat_koch        = [20,34,40,53,62]; % locations of lateral flow of koch
no_lat_koch     = length(lat_koch); %
lat_koch_data	= 'lat_koch.txt'; % boundary condition file
lat_koch_input	= importdata(lat_koch_data,',');
t_lat_koch      = lat_koch_input(:,1); % time series of lateral flow [d] 
[m_koch,~]      = find(t_lat_koch==t0);
Q_lat_koch      = lat_koch_input(m_koch:m_koch+no_data_Hec-1,2:no_lat_koch+2)'; % lateral flow
c_lat_koch      = lat_koch_input(m_koch:m_koch+no_data_Hec-1,no_lat_koch+3:(n_mob+1)*(no_lat_koch+1)+1)'; % lateral pollutant conc.
% recalculate the Q at inflow face of cell with lateral flow
Q_koch(lat_koch,:)=Q_koch(lat_koch-1,:); % reconstruct the flow at inflow face

Mat=[Q_lat_koch;c_lat_koch;Q_koch;u_koch;A_koch;dep_w_koch;tau_koch;V_koch];
dx_midpoints_koch=[0.5*L_XS_koch(1);0.5*(L_XS_koch(1:end-1)+L_XS_koch(2:end))];

% initial condition
ini_txt	= 'c0_koch.txt'; % boundary condition file
c0_koch	= importdata(ini_txt,',');
ini_ncomp = 10*ones(n_cell_koch,n_all);  % SS mass in water[g] of ncomp, mass[kg/m] of n_im, age[s] of n_mob_age, and age of n_im_age
                                       % depostion flux, bed erosion flux and bank erosion flux
ini_ncomp(:,1) = c0_koch(:,1).* V_0; % initial mass [g]
ini_ncomp(:,2) = c0_koch(:,2);%[kg/m]
ini_ncomp(:,3) = 0.5; %[d]
ini_ncomp(:,4) = 365*3; %[d]
ini_ncomp(:,5) = 0; %[kg/m]
ini_ncomp(:,6) = 0; %[kg/m]
ini_ncomp(:,7) = 0; %[kg/m]
c0=zeros(n_cell_koch*n_all,1);
% initialization
if i_tset==1   
    c0 = ini_ncomp; % use initial ncomp1 to initialize first column of c
    c0 = reshape(c0',n_cell_koch*n_all,1);
else
    res_koch =['result_koch',num2str(i_tset-1),'.mat'];
    load(res_koch);
    c0 = c_all(end,:)'; % first column is SS conc. [mg/L]
    c0(1:n_all:end,1)=c0(1:n_all:end,1).*V_0; % change SS conc. to SS mass [g]
end

% set options of ODE solver
opt_ode=odeset('reltol',1e-3,'abstol',1e-3,'jpattern',transsparse,'maxstep',dt_out);
    
disp([datestr(now) ': Begin time loop of Koch'])

[tvec,Call]=ode23s(@ode_trans,t0:dt_out:te,c0(:),opt_ode);
t_all=tvec;
c_all=Call;

% get volume and Q at each ode output step
n_all_sim=24*(te-t0)*n_out_h+1; % number of simulation steps

V_koch_out=[];
Q_koch_out=[];

for ii_v=1:1:length(V_koch(1,:))
    V_koch_ode=repmat(V_koch(:,ii_v),[1,n_out_h]);
    V_koch_out=[V_koch_out V_koch_ode];
    Q_koch_ode=repmat(Q_koch(:,ii_v),[1,n_out_h]);
    Q_koch_out=[Q_koch_out Q_koch_ode];
end

% calculate SS conc. for each cell
c_all(:,1:n_all:end)=c_all(:,1:n_all:end)./V_koch_out(:,1:n_all_sim)';% make water SS mass to SS conc. [mg/L]
Q_koch_amm=Q_koch_out(:,1:n_all_sim)'; % Q for lateral flow from kas to amm [m3/s]

% generate lateral flow file for Koch to Amm using the data of last cell
lat_KochToAmm=[t_all Q_koch_amm(:,end) c_all(:,end-n_all+1:end-n_all+n_mob)]';
fnameKTA=['lat_kochamm',num2str(i_tset),'.txt'];
fileID = fopen(fnameKTA,'w');% not consider age
fprintf(fileID,'%11.4f,	%5.3f,	%6.2f\r\n',lat_KochToAmm);
fclose(fileID);

f_koch=['result_koch',num2str(i_tset),'.mat'];
save(f_koch,'t_all','c_all')

disp([datestr(now) ': finish Sediment_Koch'])

% =========================================================================
% ODE Function
% =========================================================================
function dcdt = ode_trans(t,c)  
    % get q,u,A.... at time t
    if t==t0
        Mat_t=Mat(:,1);
    elseif t==te
        Mat_t=Mat(:,end);     
    else
        ind_Hec = t_Hec > t;
        i_ind = find(ind_Hec,1,'first')-1;
        Mat_t = Mat(:,i_ind);
    end

%     Q_ups           =Mat_t(1);
    Q_t_lat_koch    =Mat_t(2:no_lat_koch+1);
    c_ups           =Mat_t(no_lat_koch+2:n_mob+no_lat_koch+1);
    c_t_lat_koch    =Mat_t(n_mob+no_lat_koch+2:(n_mob+1)*(no_lat_koch+1));
    Q_t_koch        =Mat_t((n_mob+1)*(no_lat_koch+1)+1:(n_mob+1)*(no_lat_koch+1)+n_XS_koch*1);
%     u_t_koch        =Mat_t(n_XS_koch*1+(n_mob+1)*(no_lat_koch+1)+1:n_XS_koch*2+(n_mob+1)*(no_lat_koch+1));
    A_t_koch        =Mat_t(n_XS_koch*2+(n_mob+1)*(no_lat_koch+1)+1:n_XS_koch*3+(n_mob+1)*(no_lat_koch+1));
    dep_w_XS_koch	=Mat_t(n_XS_koch*3+(n_mob+1)*(no_lat_koch+1)+1:n_XS_koch*4+(n_mob+1)*(no_lat_koch+1));
    tau_XS_koch     =Mat_t(n_XS_koch*4+(n_mob+1)*(no_lat_koch+1)+1:n_XS_koch*5+(n_mob+1)*(no_lat_koch+1));
    % water volume in each cell
    V_t_koch        =Mat_t(n_XS_koch*5+(n_mob+1)*(no_lat_koch+1)+1:end);

    % use the average value of two cross sections as the value for the cell
    dep_w_t_koch	= 0.5*(dep_w_XS_koch(1:end-1)+dep_w_XS_koch(2:end));
    tau_t_koch      = 0.5*(tau_XS_koch(1:end-1)+tau_XS_koch(2:end));
    
    % reshape the concentration to a matrix
    c    = reshape(c,n_all,n_cell_koch)';
    % calculate current conc.
    c_ss = c(:,1:n_mob)./V_t_koch; % SS conc. [mg/L]
    % extend ss conc. to include ss from first XS
    c_up = [c_ups;c_ss]; % SS conc. [mg/L] 
       
    % =====================================================================
    % ADVECTION
    % =====================================================================
    % cell structure: n_mob=n_mob_age=n_im=n_im_age
    % 1:n_mob is mobile component; 
    % n_mob+1:n_mob+n_im is immobile component; 
    % n_mob+n_im+1:n_mob+n_im+n_mob_age is the age of corresponding mobile component
    % n_mob+n_im+n_mob_age+1:n_mob+n_im+n_mob_age+n_im_age is the age of corresponding immobile component
    % n_mob+1:n_mob+n_im is immobile component; 
    % n_mob+n_im+1:n_mob+n_im+n_mob_age is the age of corresponding mobile component
    % n_mob+n_im+n_mob_age+1:n_mob+n_im+n_mob_age+n_im_age is the age of corresponding immobile component
    
    dcdt_adv        = zeros(n_cell_koch,n_all);
    flux_cell       = 86400 * Q_t_koch .* c_up; % ss flux at each XS [g/d]
    dcdt_adv(:,1:n_mob) = flux_cell(1:end-1,:)-flux_cell(2:end,:);% advection [g/d]     
%     dcdt_adv(:,n_mob+n_im+1:n_mob+n_im+n_mob_age) = 0;% do not cumpute age change due to advection
    
    % =====================================================================
    % DISPERSION
    % =====================================================================
    D = 0.1; % set constant dispersion coefficient [m2/s]
    c_disp = [c_up;c_up(end,:)]; % extend concentration for dispersion of last cell
    J_disp = -A_t_koch*ones(1,n_mob)*D.*(c_disp(2:end,:)-c_disp(1:end-1,:))./(dx_midpoints_koch*ones(1,n_mob)); % dispersion flux [g/s]
    
    dcdt_disp = zeros(n_cell_koch,n_all);
    dcdt_disp(:,1:n_mob) = 86400*(J_disp(1:end-1,:)-J_disp(2:end,:)); % TSS dispersion [g/d]
%     dcdt_disp(:,n_mob+n_im+1:n_mob+n_im+n_mob_age) = 0; % do not cumpute age change due to dispersion 
                  
    % =====================================================================
    % REACTION
    % =====================================================================
    % deposition
    dcdt_dep = zeros(n_cell_koch,n_all); % deposition TSS [mg/L/d], Sediment [kg/m/d], age [d/d]
    r_dep = ((max(1-abs(tau_t_koch)/tau_e,0)./dep_w_t_koch)*ones(1,n_mob))*v_s.*c_ss; %[mg/L/s]
    dcdt_dep(:,1:n_mob) = - 86400 * r_dep .* V_t_koch; %[g/d]
    dcdt_dep(:,n_mob+1:n_mob+n_im) = 86.4* r_dep .* (V_t_koch*ones(1,n_im))./ (L_cell_koch*ones(1,n_im)); % [kg/m/d]
    dcdt_dep(:,n_mob+n_im+n_mob_age+n_im_age+1:n_mob+n_im+n_mob_age+n_im_age+n_dep) = ...
        dcdt_dep(:,n_mob+1:n_mob+n_im);%deposition flux [kg/m/d]
    
    % bed erosion
    dcdt_e_bed = zeros(n_cell_koch,n_all); % bed_erosion
    r_e_bed = max(max(m_pe*(abs(tau_t_koch)/tau_e-1),0),m_me*(abs(tau_t_koch)/tau_m-1+m_pe/m_me*(tau_m/tau_e-1)))*ones(1,n_im);%[kg/m/d]
    r_e_bed = r_e_bed .* (c(:,n_mob+1:n_mob+n_im)>0.1);% if sediment mass per meter is smaller than 10 gram, no bed erosion
    dcdt_e_bed(:,n_mob+1:n_mob+n_im) = -r_e_bed; %[kg/m/d]
    dcdt_e_bed(:,1:n_mob) = 1000 * r_e_bed .* (L_cell_koch*ones(1,n_im)); % [g/d]
    dcdt_e_bed(:,n_mob+n_im+n_mob_age+n_im_age+n_dep+1:n_mob+n_im+n_mob_age+n_im_age+n_dep+n_bde) = ...
        r_e_bed; % bed erosion flux [kg/m/d]

    % bank erosion
    dcdt_e_bank = zeros(n_cell_koch,n_all); % bank erosion
    r_e_bank = 1000*max(0,rho_b*(abs(tau_t_koch)-tau_b).*L_cell_koch .* dep_w_t_koch)*ones(1,n_mob)*kappa; %[g/d]
    dcdt_e_bank(:,1:n_mob) = r_e_bank;
    dcdt_e_bank(:,n_mob+n_im+n_mob_age+n_im_age+n_dep+n_bde+1:n_mob+n_im+n_mob_age+n_im_age+n_dep+n_bde+n_bke) = ...
        0.001*r_e_bank./(L_cell_koch*ones(1,n_im));%bank erosion flux [kg/m/d]
    
    % =====================================================================
    % lateral source and sink
    % =====================================================================
    dcdt_lat = zeros(n_cell_koch,n_all);
    flux_lat = 86400*Q_t_lat_koch .* c_t_lat_koch; %[g/d]
    dcdt_lat(lat_koch,1:n_mob) = flux_lat;
%     dcdt_lat(lat_koch,n_mob+n_im+1:n_mob+n_im+n_mob_age) = 0;% do not cumpute age change due to lateral flow
    
    % =====================================================================
    % AGE GROWTH
    % =====================================================================
    dcdt_age = zeros(n_cell_koch,n_all);
%     dcdt_age(:,n_mob+n_im+1:n_mob+n_im+n_mob_age+n_im_age) = 1; % age growth of TSS and Sed with time
                
    % =====================================================================
    % ADD IT UP
    % ===================================================================== 
    dcdt = reshape((dcdt_adv+dcdt_disp+dcdt_dep+dcdt_e_bed+dcdt_e_bank...
        +dcdt_lat+dcdt_age)', n_cell_koch*n_all,1);
end

% =========================================================================
% Sparsity pattern
% =========================================================================
function S = transsparse
    % computes sparsity pattern for transport
    % order of entries
    % (ii-1)*(n_mob+n_im) + 1,2...n_mob : 
    % aqueous-phase concentrations in cell ii
    % (ii-1)*(n_mob+n_im) + n_mob + 1,2...n_im :
    % immobile concentrations in cell ii
    % vectors of the sparse matrix
    ivec=zeros(n_cell_koch*n_all^2+(n_cell_koch-1)*4,1);
    jvec=zeros(size(ivec));
    avec= ones(size(ivec));
    counter=1;
    
    for ii=1:n_cell_koch
        % connection between all components within a cell
        for jj=1:n_all
            for kk=1:n_all
                ivec(counter)=(ii-1)* n_all+jj;
                jvec(counter)=(ii-1)* n_all+kk;
                counter=counter+1;
            end
        end
       
        % transport connection to upstream node
        if (ii>1)
            for jj=1:n_mob
                for kk=1:n_mob
                    ivec(counter)=(ii-1)* n_all+jj;
                    jvec(counter)=(ii-2)* n_all+kk;
                    counter=counter+1;
                end
            end
            
            for jj=1:n_mob_age
                for kk=1:n_mob
                    ivec(counter)=(ii-1)* n_all+n_mob+n_im+jj;
                    jvec(counter)=(ii-2)* n_all+kk;
                    counter=counter+1;
                end
                
                for kk=1:n_mob_age
                    ivec(counter)=(ii-1)* n_all+n_mob+n_im+jj;
                    jvec(counter)=(ii-2)* n_all+n_mob+n_im+kk;
                    counter=counter+1;
                end                
            end            
        end
        
        % transport connection to downstream node
        if (ii<n_cell_koch)
            for jj=1:n_mob
                for kk=1:n_mob
                    ivec(counter)=(ii-1)* n_all+jj;
                    jvec(counter)=ii* n_all+kk;
                    counter=counter+1;
                end
            end
        end
    end
    S=sparse(ivec,jvec,avec);
end

% =========================================================================
% Plotting function
% =========================================================================
function status=plot_profile(t,c,flag)
switch flag
    case 'init'
    case 'done'
    otherwise
    figure(1)
    disp(['t = ' num2str(t)])
    for ii=1:size(t,1)
        cplot = reshape(c(:,ii),n_all,n_cell_koch);
        subplot(n_all,1,1);
        plot(cplot(1,:)');
        ylabel('Tss [mg/L]');
        title(sprintf('t = %10.0f \n TSS',(t(ii))));

        subplot(n_all,1,2);
        plot(cplot(2,:)');
        ylabel('Sed Mass [kg/m]');
        title('Sediment Mass');
        
        subplot(n_all,1,3);
        plot(cplot(3,:)');
        ylabel('Tss age [d]');
        title('Tss age');
        
        subplot(n_all,1,4);
        plot(cplot(4,:)'/365);
        xlabel('Cell [-]');
        ylabel('Sed age [y]');
        title('Sediment age');
    end 
    drawnow;
end
status=0;
end
end
