function     [new_mesh new_velocity_g new_temperature_g new_pressure_g]=UpdateMesh(md,grids,old_mesh,old_velocity_g,old_temperature_g,old_pressure_g,bed,thickness);
%UPDATEMESH - update the mesh after a time step
%
%   This routine interpolates results onto new mesh
%
%   Usage:
%      [new_mesh new_velocity_g new_temperature_g new_pressure_g]=UpdateMesh(md,grids,old_mesh,old_velocity_g,old_temperature_g,old_pressure_g,bed,thickness)

%grids old vertical coordinate
old_z=old_mesh.grids(:,3);

%initialization
new_z             =old_z;
new_velocity_g    =old_velocity_g;
new_temperature_g =old_temperature_g;
new_pressure_g    =old_pressure_g;

for i=1:md.numberofgrids2d
        for j=1:md.numlayers
                count=(j-1)*md.numberofgrids2d+i;
                new_z(count)=bed(count)+(j-1)/(md.numlayers-1)*(thickness(count));     

		%Linearly interpolate results
		if new_z(count)>old_z(count), %The grid is upper
			if ~isnan(md.uppergrids(count)) %check that this grid is not already at the top (otherwise: same value)
				old_upper_z=old_z(md.uppergrids(count));
				coeff=(new_z(count)-old_z(count))/(old_upper_z-old_z(count));
				new_velocity_g(grids(count).grid.doflist)=old_velocity_g(grids(count).grid.doflist)+coeff*(old_velocity_g(grids(count).grid.doflist)-old_velocity_g(grids(md.uppergrids(count)).grid.doflist));
				new_temperature_g(grids(count).grid.doflist)=old_temperature_g(grids(count).grid.doflist)+coeff*(old_temperature_g(grids(count).grid.doflist)-old_temperature_g(grids(md.uppergrids(count)).grid.doflist));
				new_pressure_g(grids(count).grid.doflist)=old_pressure_g(grids(count).grid.doflist)+coeff*(old_pressure_g(grids(count).grid.doflist)-old_pressure_g(grids(md.uppergrids(count)).grid.doflist));
			end
		else                          %The grid is lower
			if ~isnan(md.lowergrids(count)) %check that this grid is not already at the bottom (otherwise: same value)
				old_lower_z=old_z(md.lowergrids(count));
				coeff=(new_z(count)-old_z(count))/(old_lower_z-old_z(count));
				new_velocity_g(grids(count).grid.doflist)=old_velocity_g(grids(count).grid.doflist)+coeff*(old_velocity_g(grids(count).grid.doflist)-old_velocity_g(grids(md.lowergrids(count)).grid.doflist));
				new_temperature_g(grids(count).grid.doflist)=old_temperature_g(grids(count).grid.doflist)+coeff*(old_temperature_g(grids(count).grid.doflist)-old_temperature_g(grids(md.lowergrids(count)).grid.doflist));
				new_pressure_g(grids(count).grid.doflist)=old_pressure_g(grids(count).grid.doflist)+coeff*(old_pressure_g(grids(count).grid.doflist)-old_pressure_g(grids(md.lowergrids(count)).grid.doflist));
			end
		end
        end
end

%Build new mesh
new_mesh.elements=old_mesh.elements;
new_mesh.grids=old_mesh.grids;
new_mesh.grids(:,3)=new_z;
