Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 20406)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 20407)
@@ -141,57 +141,57 @@
 
 		end % }}}
-function md = checkconsistency(self,md,solution,analyses) % {{{
-
-
-		md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isthermal','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isaccumulation','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.ismelt','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
-		md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
-		
-		md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
-		md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h
-		md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
-		md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0);
-		md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100);
-		md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1);
-		
-		md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
-		md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); 
-		md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 
-		md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 
-		
-		md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]);
-		md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
-		md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);
-		
-		md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'>=',0);
-		md = checkfield(md,'fieldname','smb.dzTop','NaN',1,'Inf',1,'>',0);
-		md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0);
-		md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'>=',1);
-		md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365); %10 years max 
-		md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
-		
-		switch self.aIdx,
-		case {1 2}
-			md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
-			md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
-		case 3
-			md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'>=',0,'<=',1);
-		case 4
-			md = checkfield(md,'fieldname','smb.t0wet','NaN',1,'Inf',1,'>=',15,'<=',21.9);
-			md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'>=',30,'<=',30);
-			md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'>=',7,'<=',7);
-		end
-		
-		%check zTop is < local thickness:
-		he=sum(md.geometry.thickness(md.mesh.elements),2)/size(md.mesh.elements,2);
-		if any(he<self.zTop),
-			error('SMBgemb consistency check error: zTop should be smaller than local ice thickness');
-		end
-		md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+
+
+			md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isthermal','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isaccumulation','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.ismelt','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
+			md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
+
+			md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
+			md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h
+			md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
+			md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0);
+			md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100);
+			md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1);
+
+			md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
+			md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); 
+			md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 
+			md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 
+
+			md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]);
+			md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
+			md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);
+
+			md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'>=',0);
+			md = checkfield(md,'fieldname','smb.dzTop','NaN',1,'Inf',1,'>',0);
+			md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0);
+			md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'>=',1);
+			md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365); %10 years max 
+			md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
+
+			switch self.aIdx,
+				case {1 2}
+					md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
+					md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
+				case 3
+					md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'>=',0,'<=',1);
+				case 4
+					md = checkfield(md,'fieldname','smb.t0wet','NaN',1,'Inf',1,'>=',15,'<=',21.9);
+					md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'>=',30,'<=',30);
+					md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'>=',7,'<=',7);
+			end
+
+			%check zTop is < local thickness:
+			he=sum(md.geometry.thickness(md.mesh.elements),2)/size(md.mesh.elements,2);
+			if any(he<self.zTop),
+				error('SMBgemb consistency check error: zTop should be smaller than local ice thickness');
+			end
+			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 
 		end % }}}
@@ -273,35 +273,35 @@
 			WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean');
 			
-			WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
-			
-			WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2,'scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2,'scale',1);
-		
-			WriteData(fid,'object',self,'class','smb','fieldname','aIdx','format','Integer','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','swIdx','format','Integer','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','denIdx','format','Integer','scale',1);
+			WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1);
+			
+			WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2);
+		
+			WriteData(fid,'object',self,'class','smb','fieldname','aIdx','format','Integer');
+			WriteData(fid,'object',self,'class','smb','fieldname','swIdx','format','Integer');
+			WriteData(fid,'object',self,'class','smb','fieldname','denIdx','format','Integer');
 			WriteData(fid,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double');
 
-			WriteData(fid,'object',self,'class','smb','fieldname','outputFreq','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','aSnow','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','aIce','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','cldFrac','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','t0wet','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','t0dry','format','Double','scale',1);
-			WriteData(fid,'object',self,'class','smb','fieldname','K','format','Double','scale',1);
+			WriteData(fid,'object',self,'class','smb','fieldname','outputFreq','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','aSnow','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','aIce','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','cldFrac','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','t0wet','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','t0dry','format','Double');
+			WriteData(fid,'object',self,'class','smb','fieldname','K','format','Double');
 
 			%figure out dt from forcings: 
@@ -321,3 +321,3 @@
 		end % }}}
 	end
-end		
+end
Index: /issm/trunk-jpl/src/m/classes/stressbalance.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 20406)
+++ /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 20407)
@@ -107,5 +107,4 @@
 			%singular solution
 			if ((~(any(~isnan(md.stressbalance.spcvx)) | any(~isnan(md.stressbalance.spcvy)))) & ~any(md.mask.groundedice_levelset>0)),
-				md = checkmessage(md,['model is not well posed (singular). You need at least one node with fixed velocity!']);
 				disp(sprintf('\n !!! Warning: no spc applied, model might not be well posed if no basal friction is applied, check for solution crash\n'));
 			end
Index: /issm/trunk-jpl/src/m/classes/thermal.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/thermal.m	(revision 20406)
+++ /issm/trunk-jpl/src/m/classes/thermal.m	(revision 20407)
@@ -75,9 +75,15 @@
 
 			md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]);
-			md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1);
+			md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0);
 			if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy & dimension(md.mesh)==3),
-				pos=find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN);
+
+				%Make sure the spc are less than melting point
+				TEMP=md.thermal.spctemperature(1:md.mesh.numberofvertices,:);
 				replicate=repmat(md.geometry.surface-md.mesh.z,1,size(md.thermal.spctemperature,2));
-				md = checkfield(md,'fieldname','thermal.spctemperature(find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN))','<=',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos)+10^-5,'message','spctemperature should be less or equal than the adjusted melting point');
+				pos=find(~isnan(TEMP));
+				md = checkfield(md,'fieldname','thermal.spctemperature','field',TEMP(pos)',...
+					'<=',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos)+10^-5,...
+					'message','spctemperature should be less or equal than the adjusted melting point');
+
 				md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0 1]);
 				md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[0 1]);