function [x,H_out,S_out,L_out] = SMsolver2_function(A_in,x_in,user_grid)
%SMsolver2_function
%
%           Function version of the script SMsolver2
%
%			Solves for grounding line position for marine ice sheet with
%			grounding line flux given as per Schoof, 2007, Ice sheet
%			grounding line dynamics: steady states, stability and
%			hysteresis, JGR, 112, F03S28, doi:10.1029/2006JF000664.               
%
%           Input data for this program are
%		A_Glen = value of Glen's law parameter
%		x_in = initial guess for grounding line position
%		user_grid = grid points at which ice thickness, surface and base 
%		elevation are to be computed, listed in ascending order as a row vector. This is optional,
%		default is the MISMIP mode 2 grid
%
%
%           Output from he program is a grounding line position x, grid point positions X_out (= user_grid, if specified),
%	    ice thickness H_out, ice surface elevation S_out and base elevation L_out. For grid
%           points < x, ice is grounded, otherwise floating.
%
%           Ensure that the remaining physical parameters are entered
%           correctly below under `Parameters as defined in Schoof 2007...'; 
%           By specifying the experiment number (1,2,3) as expt_no and the experiment type (a,b)
%	    as expt_type, the relevant parameter values are selected. Specify benchmark as A or B
%	    to select which of the Schoof (2007) boundary layer models you wish to use.
%
%           The depth to bed below sea level is defined in subroutine
%           SMcold_bedheight. DESPITE ITS NAME, this function computes
%           depth BELOW sea level, not height above, i.e. it returns b(x)
%           as in the MISMIP draft. SMcold_bedslope is the first derivative
%           of SMcold_bedheight, and must agree with SMcold_bedheight or
%           the Newton solver below may not converge; you will also get the
%           wrong surface elevation.
%
%           As above, SMcold_bedheight uses the expteriment number expt_no to determine which bed
%	    specification to use


%-----------------------------------------------------------------------
%Definition of global variables
%-----------------------------------------------------------------------


global m n r A C rho_g a theta

global expt_no expt_type benchmark

expt_no = 1;
expt_type = 'a';
benchmark = 'A';

n = 3;              %Glen's law exponent

if expt_type == 'a'
m = 1/3;         %sliding exponent
C = 7.624e6;        %m = 1/3 value
elseif expt_type == 'b'
m = 1;
C = 7.2082e+010;     %m = 1 value 
else error('invalid experiment type')
end

r = 0.9;            %ratio of ice to water density

rho_g = 8820;        %900 kg m^{-3} * 9.8 m s^{-2}

A = A_in;               %Glen's law parameter in SI units
                    % NB Paterson (1994) does not use SI units!

x = x_in;

a = 0.3/(365*24*3600);

if benchmark == 'A'
theta = 1;          %cold (theta = 1)
elseif benchmark == 'B'
theta = 0;	    %warm (theta = 0) boundary layer specifications,
else error('incorrect benchmark model specification')
end                    %models A and B respectively in Schoof 2007



%-----------------------------------------------------------------------
%Newton iteration parameters
%-----------------------------------------------------------------------



delta_x = 10;       %finite difference step size (metres) for gradient calculation

tolf = 1e-4;        %tolerance for finding zeros
normf = tolf + eps;

toldelta = 1e1;     %Newton step size tolerance
dx = toldelta + 1;

%-----------------------------------------------------------------------
%Newton iteration
%-----------------------------------------------------------------------


delta_x;

while normf > tolf | abs(dx) > toldelta
    f = SMcold_function(x);
    normf = abs(f);
    grad = (SMcold_function(x+delta_x)-f)/delta_x;
    dx = -f/grad;
    x = x + dx;
end
                    %computes steady state grounding line position x based
                    %on intial guess in line 35, using equations (20) and
                    %(24) in Schoof 2007


%-----------------------------------------------------------------------
%Calculate ice surface elevation S_Soln for grounded sheet
%-----------------------------------------------------------------------


if nargin < 3, user_grid = linspace(0,1.8e6,1500); end
            %user grid must be a ROW vector in metres

grounded = user_grid(user_grid < x).';

x_grid = [x; grounded(length(grounded):-1:1)];
                        %specifies the grid points on which  ice surface elevations are to be calculated
                        %these are in reverse order, i.e. counting
                        %backwards from the grounding line at x to the ice
                        %divide at zero. FIRST POINT MUST BE GROUNDING LINE
                        %POSITION x COMPUTED ABOVE, EVEN IF THIS IS NOT
                        %INCLUDED IN YOUR COMPUTATIONAL GRID

h_f = SMcold_bedheight(x)/r;
                        %computes ice thickness at the grounding line as
                        %initial condition

options = odeset('AbsTol',1e-6,'RelTol',1e-6);%odeset('AbsTol',f/1e-3);

[X_soln,H_soln] = ode45(@SMsurface,x_grid,[h_f],options);
                        %computes ice THICKNESS H_soln at points with
                        %position X_soln
                        
S_soln = H_soln - SMcold_bedheight(X_soln);
                        %computes ice surface elevation by adding bed
                        %elevation
                        
figure
hold on

plot(X_soln,-SMcold_bedheight(X_soln),'k')
plot(X_soln,S_soln,'b')

%H_soln is now rearranged and the point
%corresponding to the grounding line location is expunged to input these data into your
%computational grid:

S_soln = S_soln(length(S_soln):-1:2);

H_soln = H_soln(length(H_soln):-1:2);



%-----------------------------------------------------------------------
%Calculate ice thickness for shelf from van der Veen (1986)
%-----------------------------------------------------------------------


floating = user_grid(user_grid > x).';

q_0 = a*x;

H_soln2 = h_f*(q_0 + a*(floating-x))./(q_0^(n+1) + h_f^(n+1)*((1-r)*rho_g/4)^n*A*((q_0 + a*(floating-x)).^(n+1) - q_0^(n+1))/a).^(1/(n+1)); 

plot([x; floating],(1-r)*[h_f; H_soln2],'b');

plot([x; floating],-r*[h_f; H_soln2],'b');

H_out = [H_soln;H_soln2];
S_out = [S_soln;(1-r)*H_soln2];
L_out = [-SMcold_bedheight(grounded);-r*H_soln2];


%-----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%   Subfunctions   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------------------------------------------------------


function z = SMcold_function(x)
%Evaluates function whose zeros define x_g in `cold' steady marine sheet
%problem

global m n A C rho_g r a theta

h_f = r^(-1)*SMcold_bedheight(x);
b_x = SMcold_bedslope(x);
s = a*x;

z = theta*a + C*s.^(m+1)./(rho_g*h_f.^(m+2)) - theta*s.*b_x./h_f - A*(rho_g*(1-r)/4)^n*h_f.^(n+1);


function z = SMcold_bedheight(x)
%DEPTH OF BED BELOW SEA LEVEL, i.e. gives b(x) in Schoof 2007.
%NOTE the sign convention, SMcold_bedheight is positive if the bed is below
%sea level, negative if above sea level.

global expt_no

if expt_no == 1 || expt_no == 2
z =  -720 +778.5*(x/7.5e5);
elseif expt_no == 3
z = -(729 - 2184.8*(x/7.5e5).^2 + 1031.72*(x/7.5e5).^4 - 151.72*(x/7.5e5).^6);
else error('Invalid experiment number')
end


function z = SMcold_bedslope(x)
%FIRST DERIVATIVE OF DEPTH OF BED BELOW SEA LEVEL; must agree with
%SMcold_bedheight.

global expt_no

if expt_no == 1 || expt_no == 2
z = 778.5/7.5e5;
elseif expt_no == 3
z = -(-2184.8*2*x/7.5e5^2 + 1031.72*4*x.^3/7.5e5^4 - 151.72*6*x.^5/7.5e5^6);
else error('Invalid experiment number')
end


function z = SMsurface(x,h)

global n m epsilon a C rho_g

b_x = SMcold_bedslope(x);
s = a*x;

z = b_x - (C/rho_g)*s.^m./h.^(m+1);