function FlowSedGen

% parameters for sediment generating model
% non-urban area and ineffective urban area
rho = 1000; % [kg/m3];
g = 9.81; %[m/s2]
tau_cri = 0.30 * ones(2,14); %critical shear stress [N/m2]
C_h = 0.001 * ones(2,14); % proportional coefficient [s/m]

% effective urban area
m_max = 23; %[g/m2]
m0=zeros(1,14); %[g/m2]
k=0.33; %[1/d]
k_w=80/1E5;%[d/mm^1.5]
n_w=1.5;%[-]
thr_ur=0;%[mm/d]
r_rw=0.15;%[-] removal ratio of wwtp

% time setting
t0=datenum(2013,1,1,0,0,0); % beginning time
te=datenum(2016,12,31,23,0,0);% end time
dt=1/24; % time step [d]
t=t0:dt:te; % simulation time
n=length(t);% total time steps

% initial condition for year 2013
m_urb(1,1:14) = 2; %initial particle on urban surface[g/m2]

load p_bestfit_sel.mat;
load FlowForSedGen.mat
urb_eff=p_bestfit(end);

% read file from subcatchment_char.txt
TPDf=fopen('subcatchment_char.txt'); % Subc_No.,Total[km2],Urban[km2],Agriculture[km2],Forest[km2],slope[degree]
TPD=textscan(TPDf,'%f %f %f %f %f %f'); 
fclose(TPDf);
A_urb=TPD{3};
A_agr=TPD{4};
A_for=TPD{5};
slope=TPD{6};
A = [A_urb';A_agr';A_for'];
urb_ineff = 1 - urb_eff;
n_sub = 14;

% =========================================================================
%                        Sediment generating routine
% =========================================================================
disp([datestr(now) ': start sediment generating simulation'])
% sediment generating for rural area
slo = tan(slope'/180*pi);% slope coefficient
Slo = [slo;slo]; % slope for every landuse of every sub-catchment
for i_sedn=1:n
    tau_nonurb(:,:,i_sedn)=rho*g * q_of_sg(2:end,:,i_sedn)/1000 .* Slo;% shear stress for agricultural areas [N/m2]
    tau_ineff(:,i_sedn)=rho*g * q_of_urb_ineff_sg(:,i_sedn)/1000 .* slo';% shear stress for ineffective urban areas [N/m2]
    flux_nurb_land(:,:,i_sedn)= C_h .* max(zeros(2,14),tau_nonurb(:,:,i_sedn)-tau_cri) .* A(2:3,:)*10^6;%[kg/d]
    flux_sed_ineff(:,i_sedn)= C_h(1,:) .* max(zeros(1,14),tau_ineff(:,i_sedn)'-tau_cri(1,:)) .* A(1,:)*urb_ineff*10^6;%[kg/d]
end
flux_nurb_subc=permute(sum(flux_nurb_land),[3,2,1]);% rural flux from each subcatchment[kg/d] 

% sediment generating for effective urban area
t_dry=t(1)*ones(1,14);
for i_sed=2:n % time
    for j_sed=1:n_sub % sub-catchment
        if q_of_urb_sg(j_sed,i_sed)<=thr_ur && q_of_urb_sg(j_sed,i_sed-1)>thr_ur % the first dry day
            m0(1,j_sed)=m_urb(i_sed-1,j_sed);%[g/m2]
            t_dry(1,j_sed)=t(i_sed-1);
            m_urb(i_sed,j_sed)= (m_max - m0(1,j_sed))*(1-exp(-k*(t(i_sed)-t_dry(1,j_sed))))+m0(1,j_sed);%[g/m2]
            c_sed_urb_eff(i_sed,j_sed)=0;%[mg/L]
        elseif q_of_urb_sg(j_sed,i_sed)<=thr_ur && q_of_urb_sg(j_sed,i_sed-1)<=thr_ur % particle accumulation during dry period
            m_urb(i_sed,j_sed)= (m_max - m0(1,j_sed))*(1-exp(-k*(t(i_sed)-t_dry(1,j_sed))))+m0(1,j_sed);%[g/m2]
            c_sed_urb_eff(i_sed,j_sed)=0;%[mg/L]
        else % particle washoff during rain periods
            dmdt(i_sed,j_sed) = max(-m_urb(i_sed-1,j_sed)/dt, -k_w*q_of_urb_sg(j_sed,i_sed)^n_w*m_urb(i_sed-1,j_sed)); %[g/m2/d]
            m_urb(i_sed,j_sed) = m_urb(i_sed-1,j_sed)+ dmdt(i_sed,j_sed)*dt;%[g/m2]
            c_sed_urb_eff(i_sed,j_sed) = -1000 * dmdt(i_sed,j_sed)/q_of_urb_sg(j_sed,i_sed);%[mg/L]
        end
    end
end

flux_sed_urb_eff=c_sed_urb_eff .* q_of_urb_sg' .* repmat(A_urb',[length(c_sed_urb_eff(:,1)),1]) * urb_eff; % flux from eff urban areas[kg/d]
flux_urb_sub=flux_sed_urb_eff+flux_sed_ineff'; % flux from urban area [kg/d]
flux_urb_sub=flux_urb_sub*(1-r_rw);% the flux after removal by wwtp [kg/d]
flux_subc=flux_nurb_subc+flux_urb_sub; % [kg/d]

% ss input for low flow condition, which is caculated based on
% interflow and base flow
c_sed_if=3*2.02; % [mg/L]
c_sed_bf=3*2.02; % [mg/L]
flux_subc_if=86.4*c_sed_if*Q_if_sub_sg; % [kg/d]
flux_subc_bf=86.4*c_sed_bf*Q_bf_sub_sg; % [kg/d]
flux_sub_tot=flux_subc'+flux_subc_if+flux_subc_bf; % [kg/d]

flux_sub_out=flux_sub_tot;
c_sed_out=1/86.4*flux_sub_out ./ Q_subc_sg;

% generate lateral flow and SS conc. file for river sediment transport
disp([datestr(now) ': generating files for river sediment transport'])
Q_subc_rs=Q_subc_sg'; % flow of each subcatchment [m3/s]
c_subc_rs=c_sed_out'; % SS conc. of each subcatchment [mg/L]

% Q and c with initial condition
Q_subc_rs_w=[Q_subc_rs(1,:);Q_subc_rs(1:end-1,:)];
c_subc_rs_w=[c_subc_rs(1,:);c_subc_rs(1:end-1,:)];

%Q and ss from wwtp
Q_wwtp_Gau=0.11*ones(length(Q_subc_rs(:,1)),1);% [m3/s]
c_wwtp_Gau=3*2.02*ones(length(Q_subc_rs(:,1)),1); % [mg/L]
Q_wwtp_Hai=0.013*ones(length(Q_subc_rs(:,1)),1);% [m3/s]
c_wwtp_Hai=3*2.02*ones(length(Q_subc_rs(:,1)),1);% [mg/L]

flux_wwtp=86.4*(Q_wwtp_Gau.*c_wwtp_Gau+Q_wwtp_Hai.*c_wwtp_Hai); %kg/d

lat_amm_sub=[t' Q_subc_rs_w(:,1:5) Q_wwtp_Gau Q_subc_rs_w(:,6) Q_subc_rs_w(:,12) c_subc_rs_w(:,1:5) c_wwtp_Gau c_subc_rs_w(:,6) c_subc_rs_w(:,12)]';
fileID = fopen('lat_ammsub.txt','w');
fprintf(fileID,'%11.4f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f\r\n',lat_amm_sub);
fclose(fileID);

lat_koch_sub=[t' Q_subc_rs_w(:,7:10) Q_wwtp_Hai Q_subc_rs_w(:,11) c_subc_rs_w(:,7:10) c_wwtp_Hai c_subc_rs_w(:,11)]';
fileID = fopen('lat_koch.txt','w');
fprintf(fileID,'%11.4f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,	%5.3f,  %6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f,	%6.1f\r\n',lat_koch_sub);
fclose(fileID);

lat_kas_sub=[t' Q_subc_rs_w(:,13:14) c_subc_rs_w(:,13:14)]';
fileID = fopen('lat_kas.txt','w');
fprintf(fileID,'%11.4f,	%5.3f,	%5.3f,	%6.1f,	%6.1f\r\n',lat_kas_sub);
fclose(fileID);

% data for monthly ss input from different processes
flux_u=sum(flux_urb_sub,2);%kg/d
flux_nu=sum(flux_nurb_subc,2);
flux_bf=sum(flux_subc_bf)';
flux_if=sum(flux_subc_if)';
flux_base=flux_bf+flux_if;

flux_urban=flux_u+flux_wwtp; %kg/d
flux_slope=flux_nu+flux_base;%kg/d
save flux_cat_d.mat flux_urban flux_nu flux_base

% select monthly results
t_year=year(t);
t_mon=month(t);
for i_year=2014:2016
    j_year=i_year-2014+1;
    for i_mon=1:12
        sel_y=(t_year==i_year);
        sel_m=(t_mon==i_mon);
        sel=sel_y.*sel_m;
        flux_u_m(i_mon,j_year)=mean(flux_u(sel==1)); % [kg/d]
        flux_nu_m(i_mon,j_year)=mean(flux_nu(sel==1));
        flux_base_m(i_mon,j_year)=mean(flux_base(sel==1));
        flux_wwtp_m(i_mon,j_year)=mean(flux_wwtp(sel==1));
        flux_u_m_sum(i_mon,j_year)=sum(flux_u(sel==1))/1000/24; % [ton/month]
        flux_nu_m_sum(i_mon,j_year)=sum(flux_nu(sel==1))/1000/24;
        flux_base_m_sum(i_mon,j_year)=sum(flux_base(sel==1))/1000/24;
        flux_wwtp_m_sum(i_mon,j_year)=sum(flux_wwtp(sel==1))/1000/24;
    end    
end

save flux_cat.mat flux_u_m flux_nu_m flux_base_m flux_wwtp_m flux_urban flux_slope flux_u_m_sum flux_nu_m_sum flux_base_m_sum flux_wwtp_m_sum

% calculate annual load
for i_year=2014:2016
    j_year=i_year-2014+1;
    sel_y=(t_year==i_year);
    flux_u_y(j_year,1)=sum(flux_u(sel_y==1))/24/1000; % [ton]
    flux_nu_y(j_year,1)=sum(flux_nu(sel_y==1))/24/1000;
    flux_base_y(j_year,1)=sum(flux_base(sel_y==1))/24/1000; 
    flux_wwtp_y(j_year,1)=sum(flux_wwtp(sel_y==1))/24/1000;
end
save flux_cat_y.mat flux_u_y flux_nu_y flux_base_y flux_wwtp_y

disp([datestr(now) ': finish sediment generating simulation'])
end