function [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)
%UPDATEGEOMETRY - update the geometry during after a time step
%
%   Update geometry once given a new thickness, in the transient solutions
%
%   Usage:
%      [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)

%initialize bed
new_bed=zeros(md.numberofgrids,1);
new_surface=zeros(md.numberofgrids,1);

%First, check that thickness is >0
pos=find(new_thickness<0);
new_thickness(pos)=10; %minimum thickness

%For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. 
pos=find(md.gridonicesheet);
new_bed(pos)=bed(pos);
new_surface(pos)=bed(pos)+new_thickness(pos);

%For grids on ice shelt, we have hydrostatic equilibrium (for now)
pos=find(md.gridoniceshelf);
new_bed(pos)=-md.rho_ice/md.rho_water*new_thickness(pos);
new_surface(pos)=(1-md.rho_ice/md.rho_water)*new_thickness(pos);

%Deal with grounding line migration: later on.
