function [Temp_c,basaltemp,meltrate_s]=TempProfile1D_function(MaterialParameters,W_c,dz,surftemp,gflux,tmelt_bed,istied,meltrate_s)

% TempProfile1D_function

% Mike Wolovick, 10/21/2011

% This script solves the 1D advection-diffusion problem for a steady state
% temperature profile.  It includes basal melting/accretion as necessary.

% The method of solution is to invert a matrix whose terms are given by a
% finite difference approximation to the 1D adv-diff problem in steady
% state.

% There is one unknown for every grid center, plus one more unknown for
% either the basal temperature or the melt rate, and another for the 
% surface temperature, for a total of zsize+2 unknowns.

% This function is based on the script TempProfile1D_matrix.m, which was
% written for FullTempModel_dynamic (the 3D version of FTM_2D). 

% Changes on 12/22/2011:  the function is now compatable with variable
% vertical grid spacing.

% How this function is different from TempProfile1D_column:
% this function has vertical velocity defined on grid centers.  This means
% that the lack of mass conservation in the first guess does not also
% produce a lack of energy conservation.

% Inputs:
% MaterialParameters    structure containing the fields "cond_i", "rho_i",
%                       "specheat_i", and "latentheat"
% W_c                   vertical velocity (m/s) defined on centers.
% dz                    vertical grid size (m) (defined on centers, if it
%                       is variable)
% surftemp              surface temperature (K)
% gflux                 geothermal heat flux (W/m^2)   
% tmelt_bed             melting point at the bed (K)
% istied                controls whether basal temperature is tied to the melting
%                       point.
% meltrate_s            proscribed melt rate (m/s).

% Note:  if istied=1 the given value of meltrate_s is ignored and the
% function calculates the meltrate om its own.  meltrate_s should only be
% specified in the input if the grid cell is a thermal dam (so "melt rate"
% is negative indicating freeze on and the freeze on is limited by water
% supply).

% Outputs:
% Temp_c               temperature evaluated at grid centers (K)
% basaltemp            temperature of the bed (K)
% meltrate_s           melt rate at the bed (m/s)



% Indexing key:
% index=(contributor-1)*(totalsize)+equation
% equations for the bed, the cell just above the bed, and the cell just
% below the surface have different coefficients from the other equations.
% Those three equations are also the only ones to have nonzero constraints.
% When I describe contributors and equations, I am referring to 1:zsize+1,
% ie, the bed plus all the grid cells.


%% Preparation::

% Unpack Inputs:
if isfield(MaterialParameters,'condguess')==0
    cond_i=MaterialParameters.cond_i;
    specheat_i=MaterialParameters.specheat_i;
    rho_i=MaterialParameters.rho_i;
else
    cond_i=MaterialParameters.condguess;
    specheat_i=MaterialParameters.specheatguess;
    rho_i=MaterialParameters.rho_i;
end
latentheat=MaterialParameters.latentheat;
totalsize=length(W_c)+2;

% Pre-allocate memory:
EnergyBalanceMatrix=zeros(totalsize,totalsize);
EnergyConstraint=zeros(totalsize,1);

% Detect variable vertical grid spacing:
if length(dz)>1
    variablevertspace=1;
    dz_ud=[dz(1);.5*(dz(1:end-1)+dz(2:end));dz(end)];
else
    variablevertspace=0;
end

%% Build Components of EnergyBalanceMatrix:

% Change behavior for variable vertical grid spacing:
if variablevertspace==0
    
    % Fully within-ice conduction contributions:
    TheseInds=(linspace(3,totalsize-1,totalsize-3)'-1)*(totalsize)+linspace(2,totalsize-2,totalsize-3)'; % contribution from T(3:end-1) to T(2:end-2)
    EnergyBalanceMatrix(TheseInds)=cond_i/(dz^2);
    TheseInds=(linspace(3,totalsize-2,totalsize-4)-1)'*(totalsize)+linspace(3,totalsize-2,totalsize-4)'; % contribution from T(3:end-2) to (3:end-2)
    EnergyBalanceMatrix(TheseInds)=-2*cond_i/(dz^2);
    TheseInds=(linspace(2,totalsize-2,totalsize-3)-1)'*(totalsize)+linspace(3,totalsize-1,totalsize-3)'; % contribution from T(2:end-2) to T(3:end-1)
    EnergyBalanceMatrix(TheseInds)=cond_i/(dz^2);
    
    % Just below surface conduction:
    thisind=(totalsize-1)*totalsize+totalsize-1;    % contribution from T(end) to T(end-1)
    EnergyBalanceMatrix(thisind)=2*cond_i/(dz^2);
    thisind=(totalsize-2)*totalsize+totalsize-1;    % contribution from T(end-1) to T(end-1)
    EnergyBalanceMatrix(thisind)=-3*cond_i/(dz^2);
    
    % Just above bed conduction:
    thisind=totalsize+2;           % contribution from T(2) to T(2)
    EnergyBalanceMatrix(thisind)=-3*cond_i/(dz^2);
    
    % Surface equation and constraint:
    EnergyBalanceMatrix(end)=1;
    EnergyConstraint(end)=surftemp;
    
    % Fully within-ice advection contributions:
    % Up-moving:
    TheseW=W_c(2:end);
    TheseInds=(linspace(2,totalsize-2,totalsize-3)'-1)*(totalsize)+linspace(3,totalsize-1,totalsize-3)';  % contribution from T(2:end-2) to T(3:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW>=0))=EnergyBalanceMatrix(TheseInds(TheseW>=0))+rho_i*specheat_i*TheseW(TheseW>=0)/dz;
    TheseW=W_c;
    TheseInds=(linspace(2,totalsize-1,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';  % contribution from T(2:end-1) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW>=0))=EnergyBalanceMatrix(TheseInds(TheseW>=0))-rho_i*specheat_i*TheseW(TheseW>=0)/dz;
    % Down-moving:
    TheseW=W_c;
    TheseInds=(linspace(3,totalsize,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';  % contribution from T(3:end) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW<0))=EnergyBalanceMatrix(TheseInds(TheseW<0))-rho_i*specheat_i*TheseW(TheseW<0)/dz;
    TheseInds=(linspace(2,totalsize-1,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';   % contribution from T(2:end-1) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW<0))=EnergyBalanceMatrix(TheseInds(TheseW<0))+rho_i*specheat_i*TheseW(TheseW<0)/dz;
    
    % Basal terms:
    if istied==1
        
        % Just above bed conduction constraint:
        EnergyConstraint(2)=-2*cond_i*tmelt_bed/(dz^2);  % contribution from T(1) to T(2), as a constraint
        
        % Just above bed advection constraint:
        if W_c(1)>=0
            EnergyConstraint(2)=EnergyConstraint(2)-rho_i*specheat_i*W_c(1)*tmelt_bed/dz;  % contribution from T(1) to T(2), as a constraint
        end
        
        % Basal conduction contribution:
        thisind=totalsize+1;    % contribution from T(2) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/dz;
        
        % Basal melting contribution:
        thisind=1;          % contribution from T(1) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)-rho_i*latentheat;
        
        % Basal constraints:
        EnergyConstraint(1)=2*cond_i*tmelt_bed/dz-gflux;
        
    else
        
        % Just above bed conduction contribution:
        thisind=2;    % contribution from T(1) to T(2)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/(dz^2);
        
        % Just above bed advection contribution:
        if W_c(1)>=0
            thisind=2;    % contribution from T(1) to T(2)
            EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+rho_i*specheat_i*W_c(1)/dz;
        end
        
        % Basal conduction contributions:
        thisind=totalsize+1;   % contribution from T(2) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/dz;
        thisind=1;         % contribution from T(1) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)-2*cond_i/dz;
        
        % Basal constraint:
        EnergyConstraint(1)=rho_i*latentheat*meltrate_s-gflux;
        
    end
    
else
    
    % Fully within-ice conduction contributions:
    TheseInds=(linspace(3,totalsize-1,totalsize-3)'-1)*(totalsize)+linspace(2,totalsize-2,totalsize-3)'; % contribution from T(3:end-1) to T(2:end-2)
    EnergyBalanceMatrix(TheseInds)=cond_i./(dz(1:end-1).*dz_ud(2:end-1));
    TheseInds=(linspace(3,totalsize-2,totalsize-4)-1)'*(totalsize)+linspace(3,totalsize-2,totalsize-4)'; % contribution from T(3:end-2) to (3:end-2)
    EnergyBalanceMatrix(TheseInds)=-cond_i.*(1./(dz(2:end-1).*dz_ud(2:end-2))+1./(dz(2:end-1).*dz_ud(3:end-1)));
    TheseInds=(linspace(2,totalsize-2,totalsize-3)-1)'*(totalsize)+linspace(3,totalsize-1,totalsize-3)'; % contribution from T(2:end-2) to T(3:end-1)
    EnergyBalanceMatrix(TheseInds)=cond_i./(dz(2:end).*dz_ud(2:end-1));
    
    % Just below surface conduction:
    thisind=(totalsize-1)*totalsize+totalsize-1;    % contribution from T(end) to T(end-1)
    EnergyBalanceMatrix(thisind)=2*cond_i/(dz(end)*dz_ud(end));
    thisind=(totalsize-2)*totalsize+totalsize-1;    % contribution from T(end-1) to T(end-1)
    EnergyBalanceMatrix(thisind)=-cond_i*(1/(dz(end)*dz_ud(end-1))+2/(dz(end)*dz_ud(end)));
    
    % Just above bed conduction:
    thisind=totalsize+2;           % contribution from T(2) to T(2)
    EnergyBalanceMatrix(thisind)=-cond_i*(2/(dz(1)*dz_ud(1))+1/(dz(1)*dz_ud(2)));
    
    % Surface equation and constraint:
    EnergyBalanceMatrix(end)=1;
    EnergyConstraint(end)=surftemp;
    
    % Fully within-ice advection contributions:
    % Up-moving:
    TheseW=W_c(2:end);
    TheseDZ=dz(2:end);
    TheseInds=(linspace(2,totalsize-2,totalsize-3)'-1)*(totalsize)+linspace(3,totalsize-1,totalsize-3)';  % contribution from T(2:end-2) to T(3:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW>=0))=EnergyBalanceMatrix(TheseInds(TheseW>=0))+rho_i*specheat_i*TheseW(TheseW>=0)./TheseDZ(TheseW>=0);
    TheseW=W_c;
    TheseDZ=dz;
    TheseInds=(linspace(2,totalsize-1,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';  % contribution from T(2:end-1) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW>=0))=EnergyBalanceMatrix(TheseInds(TheseW>=0))-rho_i*specheat_i*TheseW(TheseW>=0)./TheseDZ(TheseW>=0);
    % Down-moving:
    TheseW=W_c;
    TheseDZ=dz;
    TheseInds=(linspace(3,totalsize,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';  % contribution from T(3:end) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW<0))=EnergyBalanceMatrix(TheseInds(TheseW<0))-rho_i*specheat_i*TheseW(TheseW<0)./TheseDZ(TheseW<0);
    TheseInds=(linspace(2,totalsize-1,totalsize-2)'-1)*(totalsize)+linspace(2,totalsize-1,totalsize-2)';   % contribution from T(2:end-1) to T(2:end-1)
    EnergyBalanceMatrix(TheseInds(TheseW<0))=EnergyBalanceMatrix(TheseInds(TheseW<0))+rho_i*specheat_i*TheseW(TheseW<0)./TheseDZ(TheseW<0);
    
    % Basal terms:
    if istied==1
        
        % Just above bed conduction constraint:
        EnergyConstraint(2)=-2*cond_i*tmelt_bed/(dz(1)*dz_ud(1));  % contribution from T(1) to T(2), as a constraint
        
        % Just above bed advection constraint:
        if W_c(1)>=0
            EnergyConstraint(2)=EnergyConstraint(2)-rho_i*specheat_i*W_c(1)*tmelt_bed/dz(1);  % contribution from T(1) to T(2), as a constraint
        end
        
        % Basal conduction contribution:
        thisind=totalsize+1;    % contribution from T(2) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/dz_ud(1);
        
        % Basal melting contribution:
        thisind=1;          % contribution from T(1) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)-rho_i*latentheat;
        
        % Basal constraints:
        EnergyConstraint(1)=2*cond_i*tmelt_bed/dz_ud(1)-gflux;
        
    else
        
        % Just above bed conduction contribution:
        thisind=2;    % contribution from T(1) to T(2)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/(dz(1)*dz_ud(1));
        
        % Just above bed advection contribution:
        if W_c(1)>=0
            thisind=2;    % contribution from T(1) to T(2)
            EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+rho_i*specheat_i*W_c(1)/dz(1);
        end
        
        % Basal conduction contributions:
        thisind=totalsize+1;   % contribution from T(2) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)+2*cond_i/dz_ud(1);
        thisind=1;         % contribution from T(1) to T(1)
        EnergyBalanceMatrix(thisind)=EnergyBalanceMatrix(thisind)-2*cond_i/dz_ud(1);
        
        % Basal constraint:
        EnergyConstraint(1)=rho_i*latentheat*meltrate_s-gflux;
        
    end
    
end

%% Solve Equation and Assign Output:

Temp=EnergyBalanceMatrix\EnergyConstraint;
Temp_c=Temp(2:end-1);
if istied==1
    basaltemp=tmelt_bed;
    meltrate_s=Temp(1);
else
    basaltemp=Temp(1);
end