Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 2325)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 2326)
@@ -1,17 +1,33 @@
-function bool=ismodelselfconsistent(md),
+function ismodelselfconsistent(md),
 %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
 %
 %   Usage:
-%      bool=ismodelselfconsistent(md),
+%      ismodelselfconsistent(md),
 
 %tolerance we use in litmus tests for the consistency of the model
 tolerance=10^-12;
-if (nargin~=1  )
+
+%check usage
+if nargin~=1,
 	help ismodelselfconsistent
 	error('ismodelselfconsistent error message: wrong usage');
 end
 
-%initialize output
-bool=1;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   TRANSIENT   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%check analysis
+if md.analysis_type==TransientAnalysisEnum,
+	if md.dt<=0,
+		error('model not consistent: field dt must be positive for a transient run')
+	end
+
+	%recursive call to ismodelselfconsistent
+	analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum ThermalAnalysisEnum];
+	for i=1:length(analysis),
+		md.analysis_type=analysis;
+		ismodelselfconsistent(md);
+	end
+end
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMMON CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -19,67 +35,44 @@
 %COUNTER
 if md.counter<3,
-	disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
-	bool=0;return;
+	error(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize,setelementstype)!']);
 end
 
 %NAME
 if isempty(md.name),
-	disp(['model is not correctly configured: missing name!']);
-	bool=0;return;
-end
-
-%MESH
-if md.numberofelements<=0,
-	disp(['model ' md.name ' does not have any elements!']);
-	bool=0; return;
-end
-if md.numberofgrids<=0,
-	disp(['model ' md.name ' does not have any grids!']);
-	bool=0; return;
+	error(['model is not correctly configured: missing name!']);
 end
 
 %ELEMENTSTYPE
 if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
-	disp(['Types of elements have not been set properly, run setelementstype first'])
-	bool=0;return;
+	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
 end
 if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum)  + (md.elements_type(:,1)==PattynFormulationEnum)))
-	disp(['Types of elements have not been set properly, run setelementstype first'])
-	bool=0;return;
+	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
 end
 if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
-	disp(['Types of elements have not been set properly, run setelementstype first'])
-	bool=0;return;
+	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
 end
 if strcmpi(md.type,'2d'),
 	if (ismember(PattynFormulationEnum,md.elements_type(:,1)) |  ismember(StokesFormulationEnum,md.elements_type(:,2))),
-		disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
-		bool=0;return;
+		error(['model not consistent: for a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
 	end
 end
 if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
-	disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
-	bool=0;return;
+	error(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
 end
 if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
-	disp(['coupling MacAyeal/Pattyn not implemented yet']);
-	bool=0;return;
+	error(['model not consistent: coupling MacAyeal/Pattyn not implemented yet']);
 end
 if (md.isstokes & md.analysis_type==TransientAnalysisEnum());
-	disp(['Stokes transient not implemented yet']);
-	bool=0;return;
+	error(['model not consistent: Stokes transient not implemented yet']);
 end
 
 %ICEFRONT
 if strcmpi(md.type,'2d'),
-	if size(md.pressureload,2)~=3 ,
-		disp(['pressureload in ' md.name ' shoud have three columns']);
-		bool=0;return;
-	end
+	fields={'pressureload'};
+	testsize(md,fields,[NaN 3]);
 elseif strcmpi(md.type,'3d'),
-	if size(md.pressureload,2)~=5 ,
-		disp(['pressureload in ' md.name ' shoud have five columns']);
-		bool=0;return;
-	end
+	fields={'pressureload'};
+	testsize(md,fields,[NaN 5]);
 end
 
@@ -88,12 +81,5 @@
 	'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
 	'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(isnan(md.(fields{i}))),
-			disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
+testnan(md,fields);
 
 %FIELDS >= 0 
@@ -101,52 +87,23 @@
 	'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
 	'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(md.(fields{i})<0),
-			disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
-if any(md.p<=0),
-	disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
-	bool=0; return;
-end
-
-%FIELDS ~=0
-fields={'numberofelements','numberofgrids','elements','drag_type',...
+testgreater(md,fields,0);
+
+%FIELDS > 0
+fields={'numberofelements','numberofgrids','elements','drag_type','p',...
 	'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
 	'sparsity','deltaH','DeltaH','timeacc','timedec'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(md.(fields{i})==0),
-			disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
+testgreaterstrict(md,fields,0);
 
 %SIZE NUMBEROFELEMENTS
 fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
-for i=1:size(fields,2),
-	if (size(md.(fields{i}),1)~=md.numberofelements),
-		disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
-		bool=0; return;
-	end
-end
+testlength(md,fields,md.numberofelements);
 
 %SIZE NUMBEROFGRIDS
 fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
-for i=1:length(fields),
-	if length(md.(fields{i}))~=md.numberofgrids,
-		disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
-		bool=0; return;
-	end
-end
+testlength(md,fields,md.numberofgrids);
 
 %THICKNESS = SURFACE - BED
 if any((md.thickness-md.surface+md.bed)>tolerance),
-	disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
-	bool=0; return;
+	error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
 end
 
@@ -154,6 +111,5 @@
 if md.numrifts,
 	if ~strcmpi(md.type,'2d'),
-		disp(['Models with rifts are only supported in 2d for now!']);
-		bool=0;return;
+		error(['model not consistent: models with rifts are only supported in 2d for now!']);
 	end
 end
@@ -162,6 +118,5 @@
 		if ~isempty(find(md.segmentmarkers>=2)),
 			%We have segments with rift markers, but no rift structure!
-			disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
-			bool=0; return;
+			error(['model not consistent: model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
 		end
 	end
@@ -169,40 +124,63 @@
 
 %ARTIFICIAL DIFFUSIVITY
-if ~isscalar(md.artificial_diffusivity),
-	disp('artificial_diffusivity should be a scalar (1 or 0)');
-	bool=0;return;
+if ~ismember(md.artificial_diffusivity,[0 1]),
+	error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
 end
 
 %PARAMETEROUTPUT
 if md.numoutput~=length(md.parameteroutput),
-	disp('numoutput should be the same size as parameteroutput');
-	bool=0;return;
+	error('model not consistent: numoutput should be the same size as parameteroutput');
+end
+
+%CONNECTIVITY
+if strcmpi(md.type,'2d'),
+	if md.connectivity<9, 
+		error('model not consistent: connectivity should be at least 9 for 2d models');
+	end
+end
+if strcmpi(md.type,'3d'),
+	if md.connectivity<24, 
+		error('model not consistent: connectivity should be at least 24 for 3d models');
+	end
+end
+
+%LOWMEM = 0 or 1
+if ((md.lowmem ~= 1) & (md.lowmem~=0)),
+	error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
+end
+
+%PARALLEL
+if ~strcmpi(md.cluster,'none'),
+
+	%NAN VALUES
+	fields={'time','np'};
+	testnan(md,fields);
+
+	%FIELD > 0
+	fields={'time','np'};
+	testgreaterstrict(md,fields,0);
+
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SOLUTION CHECKS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 %QMU
 if md.qmu_analysis,
 	if md.qmu_params.evaluation_concurrency~=1,
-		disp(['concurrency should be set to 1 when running dakota in library mode']);
-		bool=0;return;
+		error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
 	end
 	if ~isempty(md.part),
 		if numel(md.part)~=md.numberofgrids,
-			disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
-			bool=0;return;
+			error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
 		end
 		if find(md.part)>=md.numberofgrids,
-			disp(['user supplied partition should be indexed from 0 (c-convention)']);
-			bool=0;return;
+			error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
 		end
 		if md.npart~=md.numberofgrids,
-			disp(['user supplied partition should have same size as md.npart']);
-			bool=0;return;
-		end
-
+			error(['model not consistent: user supplied partition should have same size as md.npart']);
+		end
 	end
 	if md.eps_rel>10^-3,
-		disp(['for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
-		bool=0;return;
+		error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
 	end
 end
@@ -214,6 +192,5 @@
 	if strcmpi(md.type,'2d'),
 		if isempty(find(~md.spcvelocity(:,1:2))),
-			disp(['model ' md.name ' is totally constrained, no need to solve!']);
-			bool=0;return;
+			error(['model not consistent: model ' md.name ' is totally constrained, no need to solve!']);
 		end
 	end
@@ -226,6 +203,5 @@
 	%SINGULAR
 	if ~any(sum(md.spcvelocity(:,1:2),2)==2),
-		disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
-		bool=0;return;
+		error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
 	end
 
@@ -234,6 +210,5 @@
 		pos=find(md.thickness<=0);
 		if any(find(md.spcthickness(pos,1)==0)),
-			disp(['model ' md.name ' has some grids with 0 thickness']);
-			bool=0; return;
+			error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
 		end
 	end
@@ -245,6 +220,5 @@
 	%VELOCITIES
 	if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
-		disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
-		bool=0; return;
+		error(['model not consistent: a 3d velocity is required. Run ''diagnostic'' solution first!'])
 	end
 end
@@ -252,25 +226,23 @@
 %STEADYSTATE
 if md.analysis_type==SteadystateAnalysisEnum,
+
 	%NDT
 	if md.dt~=0,
-		disp(['For a steadystate computation, dt must be zero.']);
-		bool=0;return;
-	end
+		error(['model not consistent: for a steadystate computation, dt must be zero.']);
+	end
+
 	%PRESSURE
 	if isnans(md.pressure),
-		disp(['For a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
-		bool=0;return;
+		error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
 	end
 
 	%eps: 
 	if isnan(md.eps_rel),
-		disp(['For a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
-		bool=0;return;
+		error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
 	end
 
 	%dim: 
 	if strcmpi(md.type,'2d'),
-		disp(['For a steadystate computation, model needs to be 3d']);
-		bool=0;return;
+		error(['model not consistent: for a steadystate computation, model needs to be 3d']);
 	end
 end
@@ -281,16 +253,13 @@
 	%EXTRUSION
 	if strcmp(md.type,'2d'),
-		disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
-		bool=0;return;
+		error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
 	end
 
 	%VELOCITIES AND PRESSURE
 	if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
-		disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
-		bool=0;return;
+		error(['model not consistent: a 3d velocity is required. Run ''diagnostic'' solution first!'])
 	end
 	if (length(md.pressure)~=md.numberofgrids),
-		disp(['pressure is required. Run ''diagnostic'' solution first!'])
-		bool=0;return;
+		error(['model not consistent: pressure is required. Run ''diagnostic'' solution first!'])
 	end
 end
@@ -301,48 +270,10 @@
 	%DT and NDT
 	fields={'dt','ndt'};
-	for i=1:length(fields),
-		if any(md.(fields{i})<0),
-			disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
+	testgreaterstrict(md,fields,0);
 
 	%INITIAL TEMPERATURE, MELTING AND ACCUMULATION
-	if (length(md.temperature)~=md.numberofgrids),
-		disp(['An  initial temperature is needed for a transient thermal computation'])
-		bool=0;return;
-	end
-	if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
-		disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
-		bool=0;return;
-	end
-end
-
-%PARAMETERS
-if md.analysis_type==ParametersAnalysisEnum,
-
-	%OUTPUT
-	if ~iscell(md.parameteroutput)
-		disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
-		bool=0; return;
-	end
-	for i=1:length(md.parameteroutput)
-		if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress')  & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
-				& ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed')  & ~strcmpi(md.parameteroutput(i),'stress_surface')
-			disp(['one of the parameteroutput is not supported yet']);
-			bool=0; return;
-		end
-	end
-	%VELOCITY
-	if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
-		disp(['velocities are required!']);
-		bool=0;return;
-	end
-
-	%HUTTER
-	if any(md.elements_type(:,1)==HutterFormulationEnum); 
-		disp(['The model has Hutter''s elements. Impossible to compute parameters']);
-		bool=0;return;
-	end
+	fields={'temperature','accumulation','melting'};
+	testsize(md,fields,[md.numberofgrids 1]);
+
 end
 
@@ -352,34 +283,19 @@
 	%CONTROL TYPE
 	if ~ischar(md.control_type),
-		disp('control_type should be a string');
-		bool=0;return;
+		error('model not consistent: control_type should be a string');
 	end
 
 	%LENGTH CONTROL FIELDS
-	if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps | length(md.cm_jump)~=md.nsteps)
-		disp('maxiter, optscal, fit and cm_jump must have the length specified by nsteps')
-		bool=0;return;
-	end
+	fields={'maxiter','optscal','fit','cm_jump'};
+	testlength(md,fields,md.nsteps);
 
 	%FIT
 	if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
-		disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
-		bool=0;return;
+		error('model not consistent: wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
 	end
 
 	%OBSERVED VELOCITIES
 	fields={'vx_obs','vy_obs'};
-	for i=1:length(fields),
-		if any(length(md.(fields{i}))~=md.numberofgrids),
-			disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
-			bool=0; return;
-		end
-	end
-
-	%SINGULAR
-	if ~any(sum(md.spcvelocity(:,1:2),2)==2),
-		disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
-		bool=0;return;
-	end
+	testsize(md,fields,[md.numberofgrids 1]);
 
 	%DIRICHLET IF THICKNESS <= 0
@@ -387,6 +303,5 @@
 		pos=find(md.thickness<=0);
 		if any(find(md.spcthickness(pos,1)==0)),
-			disp(['model ' md.name ' has some grids with 0 thickness']);
-			bool=0; return;
+			error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
 		end
 	end
@@ -397,107 +312,92 @@
 	if ~strcmpi(md.cluster,'none'),
 		if md.waitonlock==0,
-			disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
-			bool=0;return;
-		end
-	end
-end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PACKAGE CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%NAN VALUES
-fields={'sparsity'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(isnan(md.(fields{i}))),
-			disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
-
-%FIELD > 0
-fields={'sparsity'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(md.(fields{i})<0),
-			disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
-
-%FIELD ~= 0
-fields={'sparsity'};
-for i=1:length(fields),
-	if ~isempty(md.(fields{i})),
-		if any(md.(fields{i})==0),
-			disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
-			bool=0; return;
-		end
-	end
-end
-
-%SPARSITY BETWEEN 0 AND 1
-if ( (md.sparsity<=0) | (md.sparsity>1)),
-	disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
-	bool=0; return;
-end
-
-%CONNECTIVITY
-if strcmpi(md.type,'2d'),
-	if md.connectivity<9, 
-		disp('connectivity should be at least 9 for 2d models');
-		bool=0;return;
-	end
-end
-if strcmpi(md.type,'3d'),
-	if md.connectivity<24, 
-		disp('connectivity should be at least 24 for 3d models');
-		bool=0;return;
-	end
-end
-
-%LOWMEM = 0 or 1
-if ((md.lowmem ~= 1) & (md.lowmem~=0)),
-	disp(['model ' md.name ' lowmem field should be 0 or 1']);
-	bool=0; return;
-end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PARALLEL CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-if ~strcmpi(md.cluster,'none'),
-
-	%NAN VALUES
-	fields={'time','np'};
-	for i=1:length(fields),
-		if ~isempty(md.(fields{i})),
-			if any(isnan(md.(fields{i}))),
-				disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
-				bool=0; return;
+			error(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
+		end
+	end
+end
+
+end %end function
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TEST FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function testlength(md,fields,fieldlength)
+	%TESTSIZE - test length of a field
+	for i=1:length(fields),
+		if length(eval(['md.' fields{i}]))~=fieldlength,
+			error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
+		end
+	end
+end
+
+function testsize(md,fields,fieldsize)
+	%TESTSIZE - test size of a field
+	for i=1:length(fields),
+		if isnan(fieldsize(1)),
+			if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
+				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
 			end
-		end
-	end
-
-	%FIELD > 0
-	fields={'time','np'};
-	for i=1:length(fields),
-		if ~isempty(md.(fields{i})),
-			if any(md.(fields{i})<0),
-				disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
-				bool=0; return;
+		elseif isnan(fieldsize(2)),
+			if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
+				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
 			end
-		end
-	end
-
-	%FIELD ~= 0
-	fields={'time','np'};
-	for i=1:length(fields),
-		if ~isempty(md.(fields{i})),
-			if any(md.(fields{i})==0),
-				disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
-				bool=0; return;
+		else
+			if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) |  (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
+				error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
 			end
 		end
 	end
-
-end
+end
+
+function testnan(md,fields)
+	%TESTNAN - test nan values of a field
+	for i=1:length(fields),
+		if any(isnan(eval(['md.' fields{i}]))),
+			error(['model not consistent: NaN values in field ' fields{i}]);
+		end
+	end
+end
+
+function testreal(md,fields)
+	%TESTREAL - test real values of a field
+	for i=1:length(fields),
+		if any(eval(['~isreal(md.' fields{i} ')'])),
+			error(['model not consistent: complex values in field ' fields{i}]);
+		end
+	end
+end
+
+function testgreaterstrict(md,fields,lowerbound)
+	%TESTGREATERSTRICT - test values of a field
+	for i=1:length(fields),
+		if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
+			error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
+		end
+	end
+end
+
+function testgreater(md,fields,lowerbound)
+	%TESTGREATER - test values of a field
+	for i=1:length(fields),
+		if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
+			error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
+		end
+	end
+end
+
+function testlessstrict(md,fields,upperbound)
+	%TESTLESSSTRICT - test values of a field
+	for i=1:length(fields),
+		if any(eval(['md.' fields{i} '>=' num2str(upperbound) ])),
+			error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
+		end
+	end
+end
+
+function testless(md,fields,upperbound)
+	%TESTLESS - test values of a field
+	for i=1:length(fields),
+		if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
+			error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
+		end
+	end
+end
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 2325)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 2326)
@@ -27,7 +27,7 @@
 md.sub_analysis_type=options.sub_analysis_type;
 
-if ~ismodelselfconsistent(md),
-	error(' '); %previous error messages should explain what is going on.
-end
+%check model consistency
+displaystring(md.debug,'\n%s\n','checking model consistency');
+ismodelselfconsistent(md),
 
 %preprocesses model before solving