function floodwave()
%
% This m-file consists of a rudimentary "flood routing" algorithm.  It
% divides a river length into node-centered elements, then computes the
% flow into and out of each element during successive time steps, changing
% the stage in each element after each time step.  Flow is computed using a
% Manning-like formula.  Input requires the mean channel slope, the
% roughness coefficient, and the flood period and stages at the upstream
% inflow node.  The input hydrograph is "triangular."  The channel is
% assumed assumed to have uniform width with no tributaries.
%
xmin = 0;
xmax = 20000; % downstream length (m)
imax = 101;
dx = (xmax - xmin)/(imax - 1); % space step (m)
meanslope = 0.001; % average downstream slope (rise/run)
N = 0.07; % roughness coefficient

stagemin = 0.3; % minimum stage (base flow stage) (m)
stagemax = 0.5; % peak flood stage (m)
flood_period = 10000; % even number (seconds)
tmax = 100000; % on the order of 10 times flood_period
dt = 2.0; % time step (seconds); not too large
kmax = tmax/dt;
tinit = tmax/10;
kinit = tinit/dt;
krise = kinit + (flood_period/2)/dt;
kfall = krise + (flood_period/2)/dt;

for i=1:imax % initial downstream distances and stages
  x(i) = (i-1)*dx;
  stageold(i) = stagemin;
  stagenew(i) = stagemin;
end

rise_rate = (stagemax - stagemin)/(flood_period/2); % rate of flood rise (m/s)
fall_rate = -rise_rate; % rate of flood fall (m/s)

for k=1:kinit % set upstream stage hydrograph in these loops
  t(k) = k*dt;
  stage0(k) = stagemin;
end
for k=kinit+1:krise
  t(k) = k*dt;
  stage0(k) = stage0(k-1) + rise_rate*dt;
end
for k=krise+1:kfall
  t(k) = k*dt;
  stage0(k) = stage0(k-1) + fall_rate*dt;
end
for k=kfall+1:kmax
  t(k) = k*dt;
  stage0(k) = stagemin;
end

%plot(t,stage0)
%axis([0 tmax 0 stagemax])

for k=1:10000 % step through flood
  stageold(1) = stage0(k);
  for i=2:imax-1
    Hin = (stageold(i-1) + stageold(i))/2; % upstream element stage
    Hout = (stageold(i) + stageold(i+1))/2; % downstream element stage
    Sin = meanslope + (stageold(i-1) - stageold(i))/dx; % upstream slope
    Sout = meanslope + (stageold(i) - stageold(i+1))/dx; % downstream slope
    stagenew(i) = stageold(i) ...
                + dt*(1/N)*(Hin^0.67)*(Sin^0.5)*Hin/dx ... % IN per unit width
                - dt*(1/N)*(Hout^0.67)*(Sout^0.5)*Hout/dx; % OUT per unit width
  end
  for i=2:imax-1
    stageold(i) = stagenew(i);
  end
end

plot(x,stageold)
axis([0 xmax 0 stagemax])