Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m	(revision 1040)
@@ -90,4 +90,5 @@
 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','s',(md.z(pos)-md.bed(pos))./md.thickness(pos));
+grids(pos)=SetStructureField(grids(pos),'grid','surface',md.surface(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
Index: /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m	(revision 1040)
@@ -87,4 +87,5 @@
 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','s',(md.z(pos)-md.bed(pos))./md.thickness(pos));
+grids(pos)=SetStructureField(grids(pos),'grid','surface',md.surface(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
Index: /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m	(revision 1040)
@@ -80,4 +80,5 @@
 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','s',(md.z(pos)-md.bed(pos))./md.thickness(pos));
+grids(pos)=SetStructureField(grids(pos),'grid','surface',md.surface(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
Index: /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m	(revision 1040)
@@ -80,4 +80,5 @@
 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','s',(md.z(pos)-md.bed(pos))./md.thickness(pos));
+grids(pos)=SetStructureField(grids(pos),'grid','surface',md.surface(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
Index: /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m	(revision 1040)
@@ -104,4 +104,5 @@
 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','s',(md.z(pos)-md.bed(pos))./md.thickness(pos));
+grids(pos)=SetStructureField(grids(pos),'grid','surface',md.surface(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
Index: /issm/trunk/src/m/solutions/ice/UpdateGridPosition.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/UpdateGridPosition.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/UpdateGridPosition.m	(revision 1040)
@@ -14,4 +14,5 @@
 	for j=1:size(count,1),
 		grids(count(j)).grid.z=bed(count(j))+grids(count(j)).grid.s*thickness(count(j));
+		grids(count(j)).grid.surface=bed(count(j))+thickness(count(j));
 	end
 end
Index: /issm/trunk/src/m/solutions/ice/diagnostic3d.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/diagnostic3d.m	(revision 1039)
+++ /issm/trunk/src/m/solutions/ice/diagnostic3d.m	(revision 1040)
@@ -80,5 +80,9 @@
 
 %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
-u_g(4:6:m_dv.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z)/md.stokesreconditioning;
+pressure=zeros(length(fem.m_dv.grids),1);
+for i=1:length(fem.m_dv.grids),
+	pressure(i)=md.rho_ice*md.g*(fem.m_dv.grids(i).grid.surface-fem.m_dv.grids(i).grid.z)/md.stokesreconditioning;
+end
+	u_g(4:6:m_dv.gridset.gsize)=pressure;
 
 if fem.isstokes,