function [State]=ice(State,Var);

% Solves flow equation with a semi-implicit
% finite differences method.
% In particular,
% h_t = (D(t) h_x)_x
% is approximated by Crank-Nicolson type method
% but D(t*) is approximated by
% 3/2 D(t_l) - 1/2D(t_l-1).

%Constants
Rice=Var.Rice; 
g=Var.g; 
rho_i=Var.rho_i;
Nx=Var.Nx;
dx=Var.dx;
Gam=Var.Gam;
n=Var.n;
Do=Var.Do;
ns=Var.ns;
bs=Var.bs;
uo=Var.uo;
hs=Var.hs;
a=Var.a;
Tb=Var.Tb;
Heq=Var.Heq(:);
rho_b=Var.rho_b;

%State variables
hh=State.hh(:); 
H=State.H(:); 
M=State.M(:)*Var.dt_ice;

%enforce boundary conditions at start
D=zeros(Nx+1,1); 

%first step of lower order
delh=(H(3:Nx+1)+hh(3:Nx+1)-H(1:Nx-1)-hh(1:Nx-1))/(2*dx);
Dold(2:Nx,1)=(Gam/(n+2))*hh(2:Nx).^(n+2).*abs(delh).^(n-1);
Dold([1 Nx+1],1)=0;
Dmid=(Dold(2:Nx+1)+Dold(1:Nx))/2;
Ao=[-Rice*Dmid(2:Nx) (ones(Nx-1,1)+Rice*(Dmid(1:Nx-1)+Dmid(2:Nx))) -Rice*Dmid(1:Nx-1)];
MM=spdiags(Ao, -1:0, Nx-1, Nx-1);
bb=M(2:Nx)+hh(2:Nx);
hh(2:Nx)=MM\bb; %note backslash: tridiagonal solve

%time-stepping loop
for year=0:Var.dt_ice:Var.ice_step;

  delh=(H(3:Nx+1)+hh(3:Nx+1)-H(1:Nx-1)-hh(1:Nx-1))/(2*dx);

  %Sediment deformation
  %Write S as a diffusion term (assuming that sediment thickness is very
  %large a/b always > h_s) 
  as=abs(rho_i*g*hh(2:Nx).*delh);
  S=    2*Do*rho_i*g*hh(2:Nx)/((ns+1)*bs) .* (as/(2*Do*uo)).^ns;
  pl=find(hs < as/bs);
  S(pl)=2*Do*rho_i*g*hh(pl+1)/((ns+1)*bs) .* (as(pl)/(2*Do*uo)).^ns.*(1-(1-bs*hs./as(pl)).^(n+1));

  
  %Ice deformation
  D(2:Nx)=S+(Gam/(n+2))*hh(2:Nx).^(n+2).*abs(delh).^(n-1);
  Dmid=.75*(D(2:Nx+1)+D(1:Nx))-.25*(Dold(2:Nx+1)+Dold(1:Nx));
  Dm2=Dmid(1:Nx-1)+Dmid(2:Nx);
  MM=spdiags([-Rice*Dmid(2:Nx) (ones(Nx-1,1)+Rice*Dm2) -Rice*Dmid(1:Nx-1)], -1:1, Nx-1, Nx-1);
  bb=M(2:Nx)+(ones(Nx-1,1)-Rice*Dm2).*hh(2:Nx) + Rice*Dmid(2:Nx).*hh(3:Nx+1)+Rice*Dmid(1:Nx-1).*hh(1:Nx-1);
  hh(2:Nx)=MM\bb; %note backslash: tridiagonal solve
  hh(hh<0)=0; 
  hh([1 end])=0;
  Dold=D;

  %Bed depression, eq8.7.2 p245 Van Der Veen
  dH=1/Tb*(Heq-H-(rho_i/rho_b)*hh);
  H=H+dH*Var.dt_ice; 
end;



%Calculate flow
hh1=hh; 
hh2=hh;

%Ice deformation
for iter=1:5,
  D(2:Nx)=S+(Gam/(n+2))*hh2(2:Nx).^(n+2).*abs(delh).^(n-1);
  Dmid=.75*(D(2:Nx+1)+D(1:Nx))-.25*(Dold(2:Nx+1)+Dold(1:Nx));
  Dm2=Dmid(1:Nx-1)+Dmid(2:Nx);
  MM=spdiags([-Rice*Dmid(2:Nx) (ones(Nx-1,1)+Rice*Dm2) -Rice*Dmid(1:Nx-1)], -1:1, Nx-1, Nx-1);
  bb=(ones(Nx-1,1)-Rice*Dm2).*hh2(2:Nx) + Rice*Dmid(2:Nx).*hh2(3:Nx+1)+Rice*Dmid(1:Nx-1).*hh2(1:Nx-1);
  hh2(2:Nx)=MM\bb; %note backslash: tridiagonal solve
  hh2(hh2<0)=0; 
end;
State.Flow=(hh2-hh1)'/iter;

State.year=State.year+year;
State.hh=hh';
State.H=H';