function fout = network_percolation_v2(parameters,sol_prev)
%fout = network_percolation_v2(parameters,sol_prev,flags)
%Run backward Euler stepping procedure for network drainage model.
%
%The model is posed on a network with n_nodes nodes and n_edges edges,
%each edge connecting an "up_node" to a "down_node". On each node, an
%effective pressure N is defined. Conduit sizes S_R and S_K are defined on
%each edge. If a flag "Kamb_storage" is set, an additional evolving cavity
%volume V_K is defined at each node.
%
%A subset of nodes is defined as boundary nodes, on which a prescribed
%effective pressure (which can depend on time and position) is defined.
%
%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', Phi_0 is a
%geometrically defined hydraulic potential, and c_1, nu, uh, S_0_R, n_Glen,
%S_P_R, epsilon_P_R, alpha and beta are constants.
%
%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)
% In addition to the parameters already defined, T is also a constant.
%
%In addition to flow along conduits S_R and S_K, 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 a flag Kamb_storage is set true, there is also storage in a cavity of
%size V_K(t), so that
%storage at node = gamma_store*(N-N_0) + V_K
%
%The additional storage volume 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
%
%Input variables for network_percolation are
%   parameters: parameter structure with user-defined fields required by
%       network_timestep_v5 
%       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 to N
%       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
%       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
%       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
%       Note that network_percolation overwrites any entries supplied as
%       part of the parameters structure under fields q_in and N_bdy_nodes
%       using the inputs q_in_fun and N_bdy_fun defined below.
%       Additional parameters used by network_percolation but not used
%       directly innetwork_timestep_v5 are
%       q_in_fun:   function handle for water supply rate function, taking
%                   arguments q_in_fun(parameters,t). Needs to output an
%                   n_nodes-by-1 vector; this is supplied as q_in to
%                   network_timestep_v5
%       N_bdy_fun: function handle for effective pressure boundary value,
%                   taking arguments N_bdy_fun(parameters,t). Outputs a
%                   length(bdy_nodes)-by-one vector that is supplied to
%                   network_timestep_v5 as N_bdy_nodes. Additional
%                   parameters can be passed to N_bdy_fun through the
%                   parameters structure inout to network_percolation
%       t_0:        initial time at which computation is started. 
%       tplot:      time interval for saving
%       stepmax:    number of these time intervals in simulation
%       step_in:    initial step value, only used in if a previous solution
%                   sol_prev is supplied to be concatenated with the
%                   solution computed in current run)
%       saveinterval: number of time intervals elapsed before solution array is
%                       saved again
%       filename:    name of file to save output to
%       init:       substructure containing initial values for N, S_R, S_K
%                   and (used only if Kamb_storage is set) V_K 
%   sol_prev:
%       An optional structure containing previous solutions, to which current output
%       is to be appended. Must contain the same fields as fout (below)
%
%Output: a structure fout, with fields
%       S_R, S_K, N, V_K (optional), t: values of S_R, S_K, N, V_K at time
%                   intervals tplot. These are appended to the equivalently
%                   named fields in sol_prev if sol_prev is supplied.
%       parameters: the parameter structure supplied to network_percolation
%The output structure fout is also saved at regular intervals to a file
%(see saveinterval and filename above)


%Generate default parameters
if nargin == 0 || ~isfield(parameters,'grid')
    [parameters.grid, parameters.L, parameters.Phi_0] = hewitt_grid_v2(201,100,.05271,.02108);
end
if nargin == 0 || ~isfield(parameters,'n_c')
    parameters.n_c = 2*ones(size(parameters.L));
end
if nargin == 0 || ~isfield(parameters,'uh')
    parameters.uh = 0.1844*ones(size(parameters.L));
end  
if nargin == 0 || ~isfield(parameters,'T')
    parameters.T = ones(size(parameters.L));
end    
if nargin == 0 || ~isfield(parameters,'nu')
    parameters.nu = ones(size(parameters.L));
end  
if nargin == 0 || ~isfield(parameters,'gamma_store')
    parameters.gamma_store = .00128*ones(size(parameters.Phi_0));
end
if nargin == 0 || ~isfield(parameters,'q_in_fun')
    q_in_fun = @qfun;
else
    q_in_fun = parameters.q_in_fun;
end
if nargin == 0 || ~isfield(parameters,'Kamb_storage')
    parameters.Kamb_storage = false;
end
if nargin == 0 || (~isfield(parameters,'uh_K') && parameters.Kamb_storage)
    parameters.uh_K = 1*ones(size(parameters.Phi_0));
end
if nargin == 0 || ~isfield(parameters,'alpha')
    parameters.alpha = 5/4;
end       
if nargin == 0 || ~isfield(parameters,'beta')
    parameters.beta = 3/2;
end    
if nargin == 0 || ~isfield(parameters,'epsilon')
    parameters.epsilon = 1e-4;
end  
if nargin == 0 || ~isfield(parameters,'n_Glen')
    parameters.n_Glen = 3;
end      
if nargin == 0 || ~isfield(parameters,'c_1')
    parameters.c_1 = 1;
end   
if nargin == 0 || ~isfield(parameters,'c_2')
    parameters.c_2 = 1;
end  
if nargin == 0 || ~isfield(parameters,'c_3')
    parameters.c_3 = 1;
end
if nargin == 0 || ~isfield(parameters,'S_0_K')
    parameters.S_0_K = 1e6;
end  
if nargin == 0 || ~isfield(parameters,'S_0_R')
    parameters.S_0_R = 1e6;
end
if nargin == 0 || ~isfield(parameters,'V_0_K')
    parameters.V_0_K = 1e6;
end 
if nargin == 0 || ~isfield(parameters,'S_P_R')
    parameters.S_P_R = 0*ones(size(parameters.L));
end
if nargin == 0 || ~isfield(parameters,'S_P_K')
    parameters.S_P_K = 0*ones(size(parameters.L));
end
if nargin == 0 || ~isfield(parameters,'epsilon_P_R')
    parameters.epsilon_P_R = 0*ones(size(parameters.L));
end
if nargin == 0 || ~isfield(parameters,'epsilon_P_K')
    parameters.epsilon_P_K = 0*ones(size(parameters.L));
end
if nargin == 0 || ~isfield(parameters,'gamma_S')
    parameters.gamma_S = 0.0267*ones(size(parameters.Phi_0));
end
if nargin == 0 || ~isfield(parameters,'r')
    parameters.r = 0.9;
end
if nargin == 0 || ~isfield(parameters,'k_leak')
    parameters.k_leak = 0;
end
if nargin == 0 || ~isfield(parameters,'N_bdy_fun')
    N_bdy_fun = @Nfun;
else
    N_bdy_fun = parameters.N_bdy_fun;
end
if nargin == 0 || ~isfield(parameters,'N_bdy_0')
    parameters.N_bdy_0 = zeros(size(parameters.grid.bdy_nodes));
end
if nargin == 0 || ~isfield(parameters,'dt')
    parameters.dt = 1e-2;
end
if nargin == 0 || ~isfield(parameters,'t_0')
    parameters.t_0 = 0;        %Initial time
end
if nargin == 0 || ~isfield(parameters,'tplot')
    parameters.tplot = 0.01132;	%Time step between values of s stored
end
if nargin > 1 && ~isempty(sol_prev)
    step_in = length(sol_prev.t)+1;	                     %Sets counter
else
    step_in = 1;
end
if nargin == 0 || ~isfield(parameters,'stepmax')
    parameters.stepmax = 200;	                        %Limits the number of time steps taken
end

%output saving parameters
if nargin == 0 || ~isfield(parameters,'filename')
    parameters.filename = 'output.mat';
end
if nargin == 0 || ~isfield(parameters,'saveinterval')
    parameters.saveinterval = 10;
end

%auxiliary time stepping variables, adjust time interval dt, pad stepmax
plotgap = round(parameters.tplot/parameters.dt);         %No of steps between storing
parameters.dt = parameters.tplot/plotgap;                %Ensures parameters.tplot is a multiple of dt
step = step_in;
t=parameters.t_0;
parameters.stepmax=parameters.stepmax+step_in-1;

%parameters for Newton iteration
srchparams.toldelta = 50*(parameters.grid.n_edges+parameters.grid.n_nodes)*sqrt(eps);
srchparams.tolgrad = 50*(parameters.grid.n_edges+parameters.grid.n_nodes)*sqrt(eps);
srchparams.tolF = 50*(parameters.grid.n_edges+parameters.grid.n_nodes)*sqrt(eps);
srchparams.itmax = 100;
srchparams.verbose = false;

%Set initial conditions
if nargin == 0 || ~isfield(parameters,'init')
    S_K = ones(size(parameters.L))+1e-3*rand(size(parameters.L));
    S_R = ones(size(parameters.L))+1e-3*rand(size(parameters.L));
    N = zeros(size(parameters.Phi_0))+1e-3*rand(size(parameters.Phi_0));
    if parameters.Kamb_storage
        V_K = ones(size(parameters.Phi_0));
        v_in = [S_K; S_R; N; V_K];
    else
        v_in = [S_K; S_R; N];
    end
else
    S_K = parameters.init.S_K;
    S_R = parameters.init.S_R;
    N = parameters.init.N;
    if parameters.Kamb_storage
        V_K = parameters.init.V_K;
        v_in = [S_R; S_K; N; V_K];
    else
        v_in = [S_R; S_K; N];
    end
end

%set up output
fout.t = zeros(1,parameters.stepmax);
fout.S_R = zeros(parameters.grid.n_edges,parameters.stepmax);
fout.S_K = zeros(parameters.grid.n_edges,parameters.stepmax);
fout.N = zeros(parameters.grid.n_nodes,parameters.stepmax);
if parameters.Kamb_storage
    fout.V_K = zeros(parameters.grid.n_nodes,parameters.stepmax);
end
if nargin > 1 &&  ~isempty(sol_prev)
    fout.t(1:step_in-1)=sol_prev.t;
    fout.S_R(:,1:step_in-1)=sol_prev.S_R;
    fout.S_K(:,1:step_in-1)=sol_prev.S_K;
    fout.N(:,1:step_in-1)=sol_prev.N;
    if parameters.Kamb_storage
        fout.V_K(:,1:step_in-1)=sol_prev.N;
    end
end
clear sol_prev
fout.parameters = parameters;

%save initial conditions
fout.t(1) = t;
fout.S_R(:,1) = S_R;
fout.S_K(:,1) = S_K;
fout.N(:,1) = N;
if parameters.Kamb_storage
    fout.V_K(:,1) = V_K;
end

while step <= parameters.stepmax-1
    t = t+parameters.dt;
    parameters.v_in_prev = v_in;
    parameters.q_in = q_in_fun(parameters,t);
    parameters.N_bdy_nodes = N_bdy_fun(parameters,t);
    [v_in, error_flag] = Newton_v2(@network_timestep_v5,@network_timestep_Jacobian_v5,v_in,parameters,srchparams);
    if error_flag
        warning('Iteration failure: convergence not achieved')
        parameters.dt = parameters.dt/2;
    end
    if t-parameters.t_0 >= parameters.tplot*(step-step_in+1)
        step = step+1   %step counter, displays automatically
        fout.t(step) = t;
        fout.S_R(:,step) = v_in(1:parameters.grid.n_edges);
        fout.S_K(:,step) = v_in(parameters.grid.n_edges+(1:parameters.grid.n_edges));
        fout.N(:,step) = v_in(2*parameters.grid.n_edges+(1:parameters.grid.n_nodes));
        if parameters.Kamb_storage
            fout.V_K(:,step) =  v_in(2*parameters.grid.n_edges+parameters.grid.n_nodes+(1:parameters.grid.n_nodes));
        end
        if floor(step/parameters.saveinterval) == step/parameters.saveinterval, 
            save(parameters.filename,'fout')
        end
    end
end        
            
end
