Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 680)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 681)
@@ -67,7 +67,7 @@
 	'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 for i=1:length(fields),
-	if ~isempty(eval(['md.' char(fields(i))])),
-		if find(isnan(eval(['md.' char(fields(i))]))),
-			disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
+	if ~isempty(eval(['md.' fields{i}])),
+		if any(isnan(eval(['md.' fields{i}]))),
+			disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
 			bool=0; return;
 		end
@@ -80,7 +80,7 @@
 	'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 for i=1:length(fields),
-	if ~isempty(eval(['md.' char(fields(i))])),
-		if find((eval(['md.' char(fields(i))]))<0),
-			disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
+	if ~isempty(eval(['md.' fields{i}])),
+		if any((eval(['md.' fields{i}]))<0),
+			disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
 			bool=0; return;
 		end
@@ -97,7 +97,7 @@
 	'sparsity','deltaH','DeltaH','timeacc','timedec'};
 for i=1:length(fields),
-	if ~isempty(eval(['md.' char(fields(i))])),
-		if find((eval(['md.' char(fields(i))]))==0),
-			disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
+	if ~isempty(eval(['md.' fields{i}])),
+		if any((eval(['md.' fields{i}]))==0),
+			disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
 			bool=0; return;
 		end
@@ -108,6 +108,6 @@
 fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
 for i=1:size(fields,2),
-	if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
-		disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
+	if (size(eval(['md.' fields{i}]),1)~=md.numberofelements),
+		disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
 		bool=0; return;
 	end
@@ -117,6 +117,6 @@
 fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
 for i=1:length(fields),
-	if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
-		disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
+	if length(eval(['md.' fields{i}]))~=md.numberofgrids,
+		disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
 		bool=0; return;
 	end
@@ -124,5 +124,5 @@
 
 %THICKNESS = SURFACE - BED
-if find((md.thickness-md.surface+md.bed)>tolerance),
+if any((md.thickness-md.surface+md.bed)>tolerance),
 	disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
 	bool=0; return;
@@ -179,5 +179,5 @@
 
 	%DIRICHLET IF THICKNESS <= 0
-	if ~isempty(find(md.thickness<=0)),
+	if any(md.thickness<=0),
 		pos=find(md.thickness<=0);
 		if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
@@ -249,6 +249,6 @@
 	fields={'dt','ndt'};
 	for i=1:length(fields),
-		if find((eval(['md.' char(fields(i))]))<0),
-			disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
+		if any((eval(['md.' fields{i}]))<0),
+			disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
 			bool=0; return;
 		end
@@ -320,6 +320,6 @@
 	fields={'vx_obs','vy_obs'};
 	for i=1:length(fields),
-		if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
-			disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
+		if any(length(eval(['md.' fields{i}]))~=md.numberofgrids),
+			disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
 			bool=0; return;
 		end
@@ -362,7 +362,7 @@
 	fields={'sparsity'};
 	for i=1:length(fields),
-		if ~isempty(eval(['md.' char(fields(i))])),
-			if find(isnan(eval(['md.' char(fields(i))]))),
-				disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
+		if ~isempty(eval(['md.' fields{i}])),
+			if any(isnan(eval(['md.' fields{i}]))),
+				disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
 				bool=0; return;
 			end
@@ -373,7 +373,7 @@
 	fields={'sparsity'};
 	for i=1:length(fields),
-		if ~isempty(eval(['md.' char(fields(i))])),
-			if find((eval(['md.' char(fields(i))]))<0),
-				disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
+		if ~isempty(eval(['md.' fields{i}])),
+			if any((eval(['md.' fields{i}]))<0),
+				disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
 				bool=0; return;
 			end
@@ -384,7 +384,7 @@
 	fields={'sparsity'};
 	for i=1:length(fields),
-		if ~isempty(eval(['md.' char(fields(i))])),
-			if find((eval(['md.' char(fields(i))]))==0),
-				disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
+		if ~isempty(eval(['md.' fields{i}])),
+			if any((eval(['md.' fields{i}]))==0),
+				disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
 				bool=0; return;
 			end
@@ -435,7 +435,7 @@
 		fields={'time','np'};
 		for i=1:length(fields),
-			if ~isempty(eval(['md.' char(fields(i))])),
-				if find(isnan(eval(['md.' char(fields(i))]))),
-					disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
+			if ~isempty(eval(['md.' fields{i}])),
+				if any(isnan(eval(['md.' fields{i}]))),
+					disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
 					bool=0; return;
 				end
@@ -446,7 +446,7 @@
 		fields={'time','np'};
 		for i=1:length(fields),
-			if ~isempty(eval(['md.' char(fields(i))])),
-				if find((eval(['md.' char(fields(i))]))<0),
-					disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
+			if ~isempty(eval(['md.' fields{i}])),
+				if any((eval(['md.' fields{i}]))<0),
+					disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
 					bool=0; return;
 				end
@@ -457,7 +457,7 @@
 		fields={'time','np'};
 		for i=1:length(fields),
-			if ~isempty(eval(['md.' char(fields(i))])),
-				if find((eval(['md.' char(fields(i))]))==0),
-					disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
+			if ~isempty(eval(['md.' fields{i}])),
+				if any((eval(['md.' fields{i}]))==0),
+					disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
 					bool=0; return;
 				end
Index: /issm/trunk/src/m/classes/public/isresultconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/isresultconsistent.m	(revision 680)
+++ /issm/trunk/src/m/classes/public/isresultconsistent.m	(revision 681)
@@ -1,3 +1,3 @@
-function bool=isresultconsistent(md)
+function bool=isresultconsistent(md,solution)
 %ISRESULTCONSISTENT: check that results are consistent
 %
@@ -6,40 +6,80 @@
 %
 %   Usage:
-%      bool=IsModelSelfConsistent(model)
+%      bool=IsModelSelfConsistent(model,solution)
 
 %tolerance we use in litmus tests for the consistency of the model
 tolerance=10^-2;
 
-if (nargin~=1  )
-	IsResultConsistentUsage();
+%some checks
+if (nargin~=2 | nargout~=1)
+	help isresultconsistent
 	error(' ');
 end
-%Check velocities i(no penalties for 2d meshes)
-if strcmp(md.analysis_type,'diagnostic'),
-
-	if strcmpi(md.type,'3d'),
-		if ~isnan(md.penalties),
-			for i=1:size(md.penalties,1)
-				for n=2:size(md.penalties,2)
-					relativex=(md.vx(md.penalties(i,1))-md.vx(md.penalties(i,n)))./(md.vx(md.penalties(i,1)));
-					if md.vx(md.penalties(i,1))==md.vx(md.penalties(i,n)), relativex=0; end
-					relativey=(md.vy(md.penalties(i,1))-md.vy(md.penalties(i,n)))./(md.vy(md.penalties(i,1)));
-					if md.vy(md.penalties(i,1))==md.vy(md.penalties(i,n)), relativey=0; end
-					if abs(relativex)>tolerance | abs(relativey)>tolerance
-						disp(['''diagnostic'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
-						bool=0; return;
-					end
-
-				end
-			end
-		end
-	end
-end
-
-if strcmp(md.analysis_type,'thermal')
-
-	%check melting 
-	if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance)
-		disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
+
+%initialize output as 1 -> success
+bool=1;
+
+%DIAGNOSTIC
+if strcmp(solution,'diagnostic'),
+
+	if strcmpi(md.type,'3d')
+		fields1={'results.diagnostic.vx','results.diagnostic.vy','results.diagnostic.vz','results.diagnostic.vel','results.diagnostic.pressure'};
+		fields2={'results.diagnostic.vel','results.diagnostic.pressure'};
+	else
+		fields1={'results.diagnostic.vx','results.diagnostic.vy','results.diagnostic.vel'};
+		fields2={'results.diagnostic.vel'};
+	end
+
+	%check size
+	if ~testsize(md,fields1,md.numberofgrids),
+		bool=0; return
+	end
+
+	%no NAN
+	if ~testnan(md,fields1),
+		bool=0; return
+	end
+
+	%check value is real
+	if ~testreal(md,fields1),
+		bool=0; return
+	end
+
+	%check value>=0
+	if ~testpositive(md,fields2),
+		bool=0; return
+	end
+
+end
+
+%CONTROL
+if strcmp(solution,'control'),
+
+	fields1={'results.control.vx','results.control.vy','results.control.vel','results.control.parameter'};
+	fields2={'results.control.vel','results.control.J'};
+
+	%check size
+	if ~testsize(md,fields1,md.numberofgrids),
+		bool=0; return
+	end
+
+	%no NAN
+	if ~testnan(md,fields1),
+		bool=0; return
+	end
+
+	%check value is real
+	if ~testreal(md,fields1),
+		bool=0; return
+	end
+
+	%check value>=0
+	if ~testpositive(md,fields2),
+		bool=0; return
+	end
+
+	%check inversed parameter
+	if any(md.control.parameter<md.mincontrolconstraint | md.control.parameter>md.maxcontrolconstraint)
+		disp(['''control'' result not consistent: inverse parameter out of range [' md.mincontrolconstraint ' ' md.maxcontrolconstraint ']']);
 		bool=0; return; 
 	end
@@ -47,56 +87,131 @@
 end
 
-if strcmp(md.analysis_type,'thermaltransient')
-
-	for iter=1:length(md.thermaltransient_results)
-
-	%check melting 
-	if any(abs(md.thermaltransient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
-		disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
-		bool=0; return;
-	end
-
-	end
-end
-
-
-if (strcmp(md.analysis_type,'transient') & strcmpi(md.type,'3d')),
-	for iter=1:length(md.transient_results)
-
-		%check melting 
-		if any(abs(md.transient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
-			disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
-			bool=0; return;
-		end
-
-
-		%check velocities
-		if ~isnan(md.penalties),
-			for i=1:size(md.penalties,1)
-				for n=2:size(md.penalties,2)
-					relativex=(md.transient_results(iter).vx(md.penalties(i,1))-md.transient_results(iter).vx(md.penalties(i,n)))./(md.transient_results(iter).vx(md.penalties(i,1)));
-					if md.transient_results(iter).vx(md.penalties(i,1))==md.transient_results(iter).vx(md.penalties(i,n)), relativex=0; end
-					relativey=(md.transient_results(iter).vy(md.penalties(i,1))-md.transient_results(iter).vy(md.penalties(i,n)))./(md.transient_results(iter).vy(md.penalties(i,1)));
-					if md.transient_results(iter).vy(md.penalties(i,1))==md.transient_results(iter).vy(md.penalties(i,n)), relativey=0; end
-					if abs(relativex)>tolerance | abs(relativey)>tolerance
-						disp(['''transient'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
-						bool=0; return;
-					end
-
-				end
-			end
-		end
-
-	end
-
-end
-
-
-%No problems, just return 1;
-bool=1;
-return;
-	
-function IsResultConsistentUsage()
-disp(' ');
-disp('   IsResultConsistent usage: flag=IsResultConsistent(model)');
-disp('         where model is an instance of the @model class, and flag a boolean');
+%THERMAL
+if strcmp(solution,'thermal')
+
+	for iter=1:length(md.results.thermal)
+
+		fields1={['results.thermal(' num2str(iter) ').temperature'],['results.thermal(' num2str(iter) ').melting']};
+		fields2={['results.thermal(' num2str(iter) ').temperature']};
+
+		%check size
+		if ~testsize(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%no NAN
+		if ~testnan(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check value is real
+		if ~testreal(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check value>=0
+		if ~testpositive(md,fields2,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check melting (<=0 via penalties)
+		if any(abs(md.results.thermal(iter).melting(md.numberofgrids2d+1:end))>tolerance)
+			disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
+			bool=0; return; 
+		end
+	end
+
+end
+
+if strcmp(solution,'transient') 
+
+	for iter=1:length(md.results.transient)
+
+		if strcmpi(md.type,'3d'),
+			fields1={['results.transient(' num2str(iter) ').vx'],['results.transient(' num2str(iter) ').vy'],...
+				['results.transient(' num2str(iter) ').vz'],['results.transient(' num2str(iter) ').vel'],...
+				['results.transient(' num2str(iter) ').bed'],['results.transient(' num2str(iter) ').surface'],...
+				['results.transient(' num2str(iter) ').thickness'],...
+				['results.transient(' num2str(iter) ').temperature'],['results.transient(' num2str(iter) ').melting']};
+			fields2={['results.transient(' num2str(iter) ').vel'],['results.transient(' num2str(iter) ').thickness'],...
+				['results.transient(' num2str(iter) ').temperature']};
+		else
+			fields1={['results.transient(' num2str(iter) ').vx'],['results.transient(' num2str(iter) ').vy'],...
+				['results.transient(' num2str(iter) ').vz'],['results.transient(' num2str(iter) ').vel'],...
+				['results.transient(' num2str(iter) ').bed'],['results.transient(' num2str(iter) ').surface'],...
+				['results.transient(' num2str(iter) ').thickness']};
+			fields2={['results.transient(' num2str(iter) ').vel'],['results.transient(' num2str(iter) ').thickness']};
+		end
+
+		%check size
+		if ~testsize(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%no NAN
+		if ~testnan(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check value is real
+		if ~testreal(md,fields1,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check value>=0
+		if ~testpositive(md,fields2,md.numberofgrids),
+			bool=0; return
+		end
+
+		%check melting (<=0 via penalties)
+		if any(abs(md.results.transient(iter).melting(md.numberofgrids2d+1:end))>tolerance)
+			disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
+			bool=0; return; 
+		end
+	end
+end
+end % end function
+
+function bool=testsize(md,fields,sizefield)
+	%TESTSIZE - test size of a field
+	bool=1;
+	for i=1:length(fields),
+		if length(eval(['md.' fields{i}]))~=sizefield
+			disp(['results not consistent: field ' fields{i} ' size should be ' num2str(md.numberofgrids)]);
+			bool=0; return;
+		end
+	end
+end
+
+function bool=testnan(md,fields,nanfield)
+	%TESTNAN - test nan values of a field
+	bool=1;
+	for i=1:length(fields),
+		if any(isnan(eval(['md.' fields{i}]))),
+			disp(['results not consistent: NaN values in field ' fields{i}]);
+			bool=0; return;
+		end
+	end
+end
+
+function bool=testreal(md,fields,realfield)
+	%TESTREAL - test real values of a field
+	bool=1;
+	for i=1:length(fields),
+		if any(eval(['~isreal(md.' fields{i} ')'])),
+			disp(['results not consistent: complex values in field ' fields{i}]);
+			bool=0; return;
+		end
+	end
+end
+
+function bool=testpositive(md,fields,positivefield)
+	%TESTPOSITIVE - test positive values of a field
+	bool=1;
+	for i=1:length(fields),
+		if any(eval(['md.' fields{i} '<0'])),
+			disp(['results not consistent: negative values in field ' fields{i}]);
+			bool=0; return;
+		end
+	end
+end
Index: /issm/trunk/src/m/classes/public/loadresultsfromdisk.m
===================================================================
--- /issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 680)
+++ /issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 681)
@@ -16,5 +16,5 @@
 %Check result is consistent
 disp(sprintf('%s\n','checking result consistency'));
-if ~isresultconsistent(md),
+if ~isresultconsistent(md,md.analysis_type),
 	%it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
 	disp('!! results not consistent correct the model !!') 
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 680)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 681)
@@ -82,5 +82,5 @@
 %Check result is consistent
 displaystring(md.debug,'%s\n','checking result consistency');
-if ~isresultconsistent(md),
+if ~isresultconsistent(md,options.analysis_type),
 	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
 end