function rt=f_dcdt(t, c, flag, fr, vx, V, alpha, taus)

% function rt=f_dcdt(t, c, flag, fr, vx, V, alpha, taus)
%
% input
% t:        time
% c:        concentrations of shallow & deep [1x2]
% flag:     Matlab ODE solver options flag
% fr:       river source
% vx:       volume exchange shallow/deep
% V:        volumes of shallow & deep reservoirs [1x2]
% alpha:    sedimentation efficiency
% taus:     shallow ocean residence time for c
%
% return value
% rt:       concentrations of shallow & deep [1x2]

rt=zeros(size(c));
phi=c(1)*V(1)/taus;

rt(1)=(fr + vx*c(2) - vx*c(1) - phi)/V(1);
rt(2)=(vx*c(1) - vx*c(2) + (1-alpha)*phi)/V(2);