Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1360)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1361)
@@ -328,4 +328,25 @@
 		end
 	end
+
+	%SINGULAR
+	if ~any(md.gridondirichlet_diag),
+		disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
+		bool=0;return;
+	end
+
+	%DIRICHLET VALUES
+	if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
+		disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
+		bool=0;return;
+	end
+
+	%DIRICHLET IF THICKNESS <= 0
+	if any(md.thickness<=0),
+		pos=find(md.thickness<=0);
+		if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
+			disp(['model ' md.name ' has some grids with 0 thickness']);
+			bool=0; return;
+		end
+	end
 end
 
Index: /issm/trunk/src/m/classes/public/modelextract2.m
===================================================================
--- /issm/trunk/src/m/classes/public/modelextract2.m	(revision 1360)
+++ /issm/trunk/src/m/classes/public/modelextract2.m	(revision 1361)
@@ -70,5 +70,13 @@
 end
 
-%element and grid position
+%kick out all elements with 3 dirichlets
+spc_elem=find(~flag_elem);
+spc_grid=sort(unique(md1.elements(spc_elem,:)));
+flag=ones(md1.numberofgrids,1);
+flag(spc_grid)=0;
+pos=find(sum(flag(md1.elements),2)==0);
+flag_elem(pos)=0;
+
+%extracted elements and nodes lists
 pos_elem=find(flag_elem);
 pos_grid=sort(unique(md1.elements(pos_elem,:)));
@@ -103,4 +111,24 @@
 	md2=md1;
 
+	%automatically modify fields
+
+	%loop over model fields
+	model_fields=fields(md1);
+	for i=1:length(model_fields),
+
+		%get field
+		field=md1.(model_fields(i));
+		fieldsize=size(field);
+
+		%size = number of grids * n
+		if fieldsize(1)==numberofgrids1
+			md2.(model_fields(i))=field(pos_grid,:);
+
+		%size = number of elements * n
+		elseif fieldsize(1)==numberofelements1
+			md2.(model_fields(i))=field(pos_elem,:);
+		end
+	end
+
 	%modify some specific fields
 
@@ -144,5 +172,5 @@
 	if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
 
-	%Dirichlets Diag
+	%Boundary conditions: Dirichlets on new boundary
 	%Catch the elements that have not been extracted
 	orphans_elem=find(~flag_elem);
@@ -151,12 +179,26 @@
 	gridstoflag1=intersect(orphans_grid,pos_grid);
 	gridstoflag2=Pgrid(gridstoflag1);
-	if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
-		md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
-		md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
-	else
-		md2.dirichletvalues_diag=zeros(numberofgrids2,2);
-		disp(' ')
-		disp('!! modelextract warning: dirichlet values should be checked !!')
-		disp(' ')
+	if ~isnan(md1.gridondirichlet_diag),
+		md2.gridondirichlet_diag(gridstoflag2)=1;
+		if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
+			md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
+			md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
+		else
+			md2.dirichletvalues_diag=zeros(numberofgrids2,2);
+			disp(' ')
+			disp('!! modelextract warning: dirichlet values should be checked !!')
+			disp(' ')
+		end
+	end
+	if ~isnan(md1.gridondirichlet_thermal),
+		md2.gridondirichlet_thermal(gridstoflag2)=1;
+		if ~isnan(md1.observed_temperature)
+			md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2); 
+		else
+			md2.dirichletvalues_thermal=zeros(numberofgrids2,2);
+			disp(' ')
+			disp('!! modelextract warning: dirichlet values should be checked !!')
+			disp(' ')
+		end
 	end
 
@@ -215,28 +257,4 @@
 	md2.strainrate=NaN;
 
-	%other fields
-
-	%loop over model fields
-	model_fields=fields(md1);
-	for i=1:length(model_fields),
-
-		%get field
-		field=md1.(model_fields(i));
-
-		%size = number of grids * 1
-		if (size(field,1)==numberofgrids1 & size(field,2)==1),
-			md2.(model_fields(i))=field(pos_grid);
-
-		%size = number of grids * 2
-		elseif (size(field,1)==numberofgrids1 & size(field,2)==2), 
-			md2.(model_fields(i))=[field(pos_grid,1) field(pos_grid,2)];
-
-		%size = number of elements * 1
-		elseif (size(field,1)==numberofelements1 & size(field,2)==1),
-			md2.(model_fields(i))=field(pos_elem);
-
-		end
-	end
-
 %Keep track of pos_grid and pos_elem
 md2.extractedgrids=pos_grid;