function fout = network_timestep_v5(v_in,parameters)
%Computes functions whose zeros define backward Euler step for drainage
%network problem.
%Input variables are:
%v_in:          Concatenated column vectors for S_K, S_R, N and
%               (optionally) V_K, see below
%parameters:    Parameter structure with the following fields
%   grid:       substructure with fields
%               n_nodes: number of nodes in network
%               n_edges: number of network edges
%               up_node: list of upstream nodes for each edge
%               down_node: list of downstream nodes for each edge
%                   (up_node and down_node together give the data expected
%                   in a connectivity array)
%               bdy_nodes: list of nodes on the domain boundary, where
%               Dirichlet conditions are to be applied
%   dt:         Size of time step
%   L:          List (n_edges-by-one vector) of edge lengths
%   n_c:        List (n_edges-by-one vector) of conduits per edge
%   Phi_0:      Reduced hydraulic potential at each node (n_nodes-by-one
%               vector), given by rho_i g s + (rho_w-rho_i) g b, s being
%               local ice surface elevation and b local bed elevation at
%               the node. rho_i and rho_w are ice and water density and g
%               is acceleration due to gravity
%   N_bdy_nodes:Effective pressure at boundary nodes (column vector with
%               same dimensions as bdy_nodes)
%   c_1:        Melting term is c_1*Q*Psi, where Q is discharge and Psi
%               hydraulic gradient along each network edge
%   c_2:        Creep closure is c_2*S*N^n_Glen, where S is conduit
%               cross-section or storage cavity size, N is effective
%               pressure
%   c_3:        Coefficient in Manning or Darcy-Weisbach law, discharge is
%               Q = c_3*S^alpha*Psi^(beta-1)
%   alpha:      Exponent relating conduit cross-section to discharge
%   beta:       Exponent relating Hydraulic gradient to discharge
%   n_Glen:     Glen's law exponent in creep closure relation
%   epsilon:    Regualrizer for the square-root dependence of discharge on
%               hydraulic gradient
%   uh:         List (n_edges-by-one vector) of 'cavity' opening rates for
%               each edge
%   S_0_R:      cut-off size for cavity opening term for R conduits
%   S_0_K:      cut-off size for cavity opening term for pure cavity K
%               conduits
%   V_0_K:      cut-off size for cavity opening term for storage cavity K
%               terms
%   S_P_R:      percolation cut-off size for R conduits
%   S_P_K:      percolation cut-off size for K conduits
%   epsilon_P_R: regularizer for percolation cut-off for R conduits
%   epsilon_P_K: regularizer for precolation cut-off for K conduits
%   T:          List (n_edges-by-one vector) of 'cavity' tortuosities along
%               each edge 
%   nu:         List (n_edges-by-one vector) of reduction factors in
%               opening rates from 'cavity' to  'channel' conduits for each
%               edge
%   gamma_store:List (n_nodes-by-one vector) of 'compressibiities' of
%               storage system at each node
%   gamma_S:    List(n_nodes-y-one-vector) of parameters defining contribution of storage in conduits to mass balance at each node
%   r:          ice to water density ratio
%   k_leak:     leakage permeability
%   q_in:       List (n_nodes-by-one vector) of water supply rates at each
%               node
%   Kamb_storage:Flag that sets whether a Kamb-type storage cavity volume
%               needs to be solved for at each node
%   uh_K:       List (n_nodes-by-one vector) of storage cavity opening
%               rates at each node
%   v_in_prev:  vector v_in at last time step
%
%The vectors S_R and S_K in v_in have dimensions n_edges-by-one, N and
%(optionally) V_K are n_nodes-by-one, and v_in = [S_R; S_K; N; V_K] if
%Kamb_storage is true, v_in = [S_R; S_K,N] otherwise. S_R and S_K are
%defined on network edges, N and V_K are defined at nodes.
%
%Physically, S_R is a generic drainage element as in Schoof (Ice sheet
%acceleration driven by melt supply variability, 2010), evolving as
%(constraint #1)
% dS_R / dt - ( c_1*Q_R*Psi + nu*uh*(1-S_R/S_0_R -
%           c_2*S_R*|N_edge|^(n_Glen-1)*N_edge ) = 0
%where N_edge is the mean of N at the neighbouring nodes,
% Q_R = c_3*heaviside(S_R-S_P_R)*((S_R-S_P_R).^2+epsilon_P_R^2)^(1/2)-epsilon_P_R)^alpha*|Psi_abs|^(beta-2)*Psi
%and
% Psi = (N(up_node) - N(down_node) + Phi_0(down_node) - Phi_0(up_node))/L,
% Psi_abs = sqrt(Psi*2+epsilon^2)
%where L is node spacing between 'up_node' and 'down_node'
%S_K is a cavity element with no melt enlargement effect, so (constraint
%#2)
% dS_K / dt - ( uh*(1-S_K/S_0_K) - c_2*S_K*|N_edge|^(n_Glen-1)*N_edge ) = 0
% Discharge in the S_K conduits is given by
% Q_K = c_3*heaviside(S_K-S_P_K)*((S_K-S_P_K).^2+epsilon_P_K^2)^(1/2)-epsilon_P_K)^alpha*|Psi_abs_K|^(beta-2)*Psi_K
%where Psi_K = Psi/T and Psi_abs_K = sqrt(Psi_K^2+epsilon^2)
%S_R and S_K are defined on the same network edges. In addition, we permit a leakage flux as
%Q_leak = k_leak*Psi
%The total discharge
%along an edge is
%Q = Q_R + (n_c-1)*Q_K+Q_leak
%
%Conservation of mass at each node gives (constraint #3)
% d(storage at node) / dt + 1/2* gamma_S * sum (L*(dS_R/dt+(n_c-1)*T*dS_K/dt)) -  ( q_in +
% sum Q_in - sum Q_out - sum L*melt) = 0
%where Q_in is  discharges along an edge into the node, and Q_out is
%discharge along an edge out of the node, and the sum over dS_R/dt nd
%S_K/dt is taken over all edges into or out of the node. "melt" is given by
% 1/2*gamma_S*r*L*Q*Psi
% and the sum is taken over all edges into or out of a node; r is the
% ice-to-water density ratio.
%Storage at a node has an 'elastic' componennt gamma_store*(N-N_0), where
%N_0 is a local constant
%
%If Kamb_storage is true, there is also storage in a cavity of size V_K(t).
%This then evolves according to (constraint #4)
%dV_K / dt - ( uh_K*(1-V_K/V_0_K) - c_3*V_K*|N|^(n_Glen-1)*N ) = 0
%
%The output of this function is a concatenated list of of the left-hand
%sides of the constraints given
%above, evaluated on each edge (first two constraints) and on each node
%(mass conservation and optionally storage cavity evolution constraints),
%concatenated in the order given above.
%
%Works as network_timestep_v3 but adds functionality for storage in conduits,
%percolation cut-off and a leakage flux

%unpack parameters
n_nodes = parameters.grid.n_nodes;   %number of nodes
n_edges = parameters.grid.n_edges;   %number of edges
up_node = parameters.grid.up_node;   %list (n_edges-by-1 vector) of 'upstream' node for each edge
down_node = parameters.grid.down_node;   %list (n_edges-by-1 vector) of 'downstream' nodes for each edge
bdy_nodes = parameters.grid.bdy_nodes;   %list (n_bdy_nodes-by-1 vector) of nodes at domain boundary

dt = parameters.dt;         %time step

L = parameters.L;           %list (n_edges-by-1 vector) of edge lengths
n_c = parameters.n_c;       %list (n_edges-by-1 vector) of number of conduits along edge
Phi_0 = parameters.Phi_0;   %reduced hydraulic potential rho_i s + (rho_w-rho_i) b at each node (n_nodes-by-1 vector) 
N_bdy_nodes = parameters.N_bdy_nodes;   %Dirichlet conditions (n_bdy_nodes-by-1 vector) for effective pressure at domain boundary

c_1 = parameters.c_1;       %relates opening rate to Q*Psi
c_2 = parameters.c_2;       %relates closure rate to s*N*n
c_3 = parameters.c_3;       %relates discharge Q to S^alpha*Psi^beta
alpha = parameters.alpha;   %exponent in dependence of discharge Q on cross-section S
beta = parameters.beta;     %exponent in dependence of discharge Q on hydraulic gradient Psi, same convention as in Schoof et al 2012 / Hewitt et al 2012 so 'frozen-time' problenm for N only becomes beta-Laplacian (i.e. p-Laplacian with p=beta)
epsilon = parameters.epsilon;   %regularization parameter for hydraulic gradient
n_Glen = parameters.n_Glen; %Glen's law exponent
uh = parameters.uh;         %list (n_edges-by-1 vector) of cavity opening rates for each network edge
S_0_R = parameters.S_0_R;       %cut-off size for cavity opening
S_0_K  = parameters.S_0_K;  %cut-off size for Kamb cavity opening
S_P_R = parameters.S_P_R;   %percolation cut-off size for R conduit
S_P_K = parameters.S_P_K;   %percolation cut-off size for K conduit
epsilon_P_R = parameters.epsilon_P_R;   %percolation cut-off regularizer for R-conduits;
epsilon_P_K = parameters.epsilon_P_K;   %percolation cut-off regularizer for K-conduits;
T = parameters.T;           %list (n_edges-by-1 vector) of tortuosities for 'cavities' along each network edge
nu = parameters.nu;         %list (n_edges-by-1 vector) of step size ratios for 'channel' along each network edge
gamma_S = parameters.gamma_S;   %contribution of conduit evolution to mass balance
gamma_store = parameters.gamma_store;   %list (n_nodes-by-1 vector) of 'ealstic' water storage capacities at each node (i.e. storage is affine in N with coefficient gamma_store)
r = parameters.r;           %scalar density ratio value
q_in = parameters.q_in;     %list (n_nodes-by-1 vector) of water supply rates to each node
k_leak = parameters.k_leak; %leakage diffusivity

Kamb_storage = parameters.Kamb_storage;     %Boolean that determines whether Kamb-style storage cavities are included at each node
if Kamb_storage
    V_0_K = parameters.V_0_K;   %cut-off size for storage cavity opening
    uh_K = parameters.uh_K; %list (n_nodes-bv-1 vector) of cavity opening rates for storage cavities at each node
end

%get state variable v_in from previous time step
v_in_prev = parameters.v_in_prev;

%unpack input variable v_in
if (Kamb_storage && length(v_in) ~= 2*(n_edges+n_nodes)) || (~Kamb_storage && length(v_in) ~= 2*n_edges+n_nodes), error('input size incompatible with network geometry'), end
                                                        %check input size
S_R = v_in(1:n_edges);
S_K = v_in(n_edges+1:2*n_edges);
N = v_in(2*n_edges+1:2*n_edges+n_nodes);
N(bdy_nodes) = N_bdy_nodes;
if Kamb_storage
    V_K = v_in(2*n_edges+n_nodes+1:2*n_edges+2*n_nodes);  %optional Kamb-style storage cavity at node
end

%unpack old input variable
S_R_prev = v_in_prev(1:n_edges);
S_K_prev = v_in_prev(n_edges+1:2*n_edges);
N_prev = v_in_prev(2*n_edges+1:2*n_edges+n_nodes);
if Kamb_storage
    V_K_prev = v_in_prev(2*n_edges+n_nodes+1:2*n_edges+2*n_nodes);  %optional Kamb-style storage cavity at node
end

%initialize fout
fout = zeros(size(v_in));

%compute hydraulic gradient along each network edge
Psi = (Phi_0(up_node) - N(up_node) - Phi_0(down_node) + N(down_node))./L;
Psi_abs_R = (Psi.^2+epsilon.^2).^(1/2);
Psi_abs_K = (Psi.^2./T.^2+epsilon.^2).^(1/2);

%compute effective pressure along each network edge
N_edge = (N(up_node)+N(down_node))/2;

%compute 'channel' flux along each edge
Q_R = c_3*(((S_R-S_P_R).^2+epsilon_P_R.^2).^(1/2)-epsilon_P_R).^alpha.*Psi_abs_R.^(beta-2).*Psi;
Q_R(S_R<S_P_R)=0;
%compute 'cavity' flux along each edge
Q_K = c_3*(((S_K-S_P_K).^2+epsilon_P_K.^2).^(1/2)-epsilon_P_K).^alpha.*Psi_abs_K.^(beta-2).*Psi./T;
Q_K(S_K<S_P_K)=0;
%compute leakage flux along each edge
Q_leak = k_leak*Psi;

%compute channel opening rate
fout(1:n_edges) = (S_R-S_R_prev)/dt - c_1*Q_R.*Psi - nu.*uh.*(1-S_R/S_0_R) + c_2*S_R.*abs(N_edge).^(n_Glen-1).*N_edge;

%compute cavity opening rate
fout(n_edges+1:2*n_edges) = (S_K-S_K_prev)/dt - uh.*(1-S_K/S_0_K) + c_2*S_K.*abs(N_edge).^(n_Glen-1).*N_edge;

%compute net fluxes out of each node
netflux = accumarray(up_node,Q_R+(n_c-1).*Q_K+Q_leak,[n_nodes 1]) - accumarray(down_node,Q_R+(n_c-1).*Q_K+Q_leak,[n_nodes 1]);

%compute net rate of volume change associated with node
netstorage = 1/(2*dt)*gamma_S.*accumarray(up_node,L.*((S_R-S_R_prev)+T.*(n_c-1).*(S_K-S_K_prev)),[n_nodes 1]) ...
    + 1/(2*dt)*gamma_S.*accumarray(down_node,L.*((S_R-S_R_prev)+T.*(n_c-1).*(S_K-S_K_prev)),[n_nodes 1]);

meltrate = r/2*gamma_S.*accumarray(up_node,L.*Q_R.*Psi,[n_nodes 1]) +  r/2*gamma_S.*accumarray(down_node,L.*Q_R.*Psi,[n_nodes 1]);

%impose mass conservation at each node
fout(2*n_edges+1:2*n_edges+n_nodes) = -gamma_store.*(N-N_prev)/dt + netstorage + netflux - q_in - meltrate;

%correct for Kamb storage, compute Kamb storage change if necessary
if Kamb_storage
    fout(2*n_edges+1:2*n_edges+n_nodes) = fout(2*n_edges+1:2*n_edges+n_nodes) + (V_K-V_K_prev)/dt;
    fout(2*n_edges+n_nodes+1:2*n_edges+2*n_nodes) = (V_K-V_K_prev)/dt - uh_K.*(1-V_K./V_0_K) + c_3*V_K.*abs(N).^(n_Glen-1).*N;
end

%impose Dirichlet conditions at boundary nodes instead of mass conservation
fout(2*n_edges+bdy_nodes) = N(bdy_nodes)-N_bdy_nodes;

fout = real(fout);

