Index: /issm/trunk/src/m/classes/model/model.m
===================================================================
--- /issm/trunk/src/m/classes/model/model.m	(revision 9980)
+++ /issm/trunk/src/m/classes/model/model.m	(revision 9981)
@@ -232,4 +232,5 @@
 			 if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end
 			 if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end
+			 if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end
 			 if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end
 			 if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end
@@ -270,4 +271,6 @@
 			 if isfield(structmd,'z'), md.mesh.z=structmd.z; end
 			 if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end
+          if isfield(structmd,'pressureload'), md.diagnostic.icefront=structmd.pressureload; end
+			 if isfield(structmd,'diagnostic_ref'), md.diagnostic.referential=structmd.diagnostic_ref; end
 
 			 %Field changes
@@ -287,26 +290,40 @@
 				 pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); 
 			 end
-			 if isfield(structmd,'spcvelocity'), 
+			 if isfield(structmd,'spcvx'), 
 				 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+				 pos=find(~isnan(structmd.spcvx)); md.diagnostic.spcvx(pos)=structmd.spcvx(pos); 
+			 end
+			 if isfield(structmd,'spcvy'),
 				 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+				 pos=find(~isnan(structmd.spcvy)); md.diagnostic.spcvy(pos)=structmd.spcvy(pos);     
+			 end
+			 if isfield(structmd,'spcvz'),
 				 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-				 pos=find(structmd.spcvelocity(:,1)); md.diagnostic.spcvx(pos)=structmd.spcvelocity(pos,4); 
-				 pos=find(structmd.spcvelocity(:,2)); md.diagnostic.spcvy(pos)=structmd.spcvelocity(pos,5); 
-				 pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); 
+				 pos=find(~isnan(structmd.spcvz)); md.diagnostic.spcvz(pos)=structmd.spcvz(pos);     
 			 end
 			 if ismember(structmd.pressureload(end,end),[118 119 120]),
-				 pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,end)=0;
-				 pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,end)=1;
-				 pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,end)=2;
-			 end
-			 if (structmd.elements_type(end,end)>100),
-				 pos=find(structmd.elements_type==59); md.mesh.elements_type(pos,end)=0;
-				 pos=find(structmd.elements_type==55); md.mesh.elements_type(pos,end)=1;
-				 pos=find(structmd.elements_type==56); md.mesh.elements_type(pos,end)=2;
-				 pos=find(structmd.elements_type==60); md.mesh.elements_type(pos,end)=3;
-				 pos=find(structmd.elements_type==62); md.mesh.elements_type(pos,end)=4;
-				 pos=find(structmd.elements_type==57); md.mesh.elements_type(pos,end)=5;
-				 pos=find(structmd.elements_type==58); md.mesh.elements_type(pos,end)=6;
-				 pos=find(structmd.elements_type==61); md.mesh.elements_type(pos,end)=7;
+				 pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,:)=0;
+				 pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,:)=1;
+				 pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,:)=2;
+			 end
+			 if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50,
+				 pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0;
+				 pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1;
+				 pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2;
+				 pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3;
+				 pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4;
+				 pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5;
+				 pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6;
+				 pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7;
+			 end
+			 if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50,
+				 pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0;
+				 pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1;
+				 pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2;
+				 pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3;
+				 pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4;
+				 pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5;
+				 pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6;
+				 pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7;
 			 end
 			 if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law),
@@ -335,4 +352,9 @@
 				 pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502;
 				 pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503;
+			 end
+			 if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2,
+					 md.thermal.stabilization=2;
+					 md.prognostic.stabilization=1;
+					 md.balancethickness.stabilization=1;
 			 end