Index: /issm/trunk/src/m/classes/balancethickness.m
===================================================================
--- /issm/trunk/src/m/classes/balancethickness.m	(revision 9859)
+++ /issm/trunk/src/m/classes/balancethickness.m	(revision 9860)
@@ -43,8 +43,7 @@
 			if solution~=BalancethicknessSolutionEnum, return; end
 
-			if ~ismember(md.balancethickness.stabilization,[0 1 3]),
-				checkmessage('model not consistent: balancethickness.stabilization should be a scalar (0 or 1 or 3)');
-			end
-
+			checkfield(md,'balancethickness.spcthickness','forcing',1);
+			checkfield(md,'balancethickness.thickening_rate','size',[md.mesh.numberofvertices 1],'NaN',1);
+			checkfield(md,'balancethickness.stabilization','size',[1 1],'values',[0 1 3]);
 		end % }}}
 		function disp(obj) % {{{
Index: /issm/trunk/src/m/classes/basalforcings.m
===================================================================
--- /issm/trunk/src/m/classes/basalforcings.m	(revision 9859)
+++ /issm/trunk/src/m/classes/basalforcings.m	(revision 9860)
@@ -37,19 +37,13 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			fields={'geothermalflux'};
-			%checksize(md,'basalforcings',fields,[md.mesh.numberofvertices 1]);
 			if ismember(PrognosticAnalysisEnum,analyses),
-				fields={'melting_rate'};
-				checkforcing(md,'basalforcings',fields);
+				checkfield(md,'basalforcings.melting_rate','NaN',1,'forcing',1);
 			end
 			if ismember(BalancethicknessAnalysisEnum,analyses),
-				fields={'melting_rate'};
-				checksize(md,'basalforcings',fields,[md.mesh.numberofvertices 1]);
-				checknan(md, 'basalforcings',fields);
+				checkfield(md,'basalforcings.melting_rate','NaN',1,'size',[md.mesh.numberofvertices 1]);
 			end
 			if ismember(ThermalAnalysisEnum,analyses),
-				fields={'geothermalflux'};
-				checksize(md,'basalforcings',fields,[md.mesh.numberofvertices 1]);
-				checknan(md, 'basalforcings',fields);
+				checkfield(md,'basalforcings.melting_rate','NaN',1,'forcing',1);
+				checkfield(md,'basalforcings.geothermalflux','NaN',1,'forcing',1,'>=',0);
 			end
 		end % }}}
Index: /issm/trunk/src/m/classes/constants.m
===================================================================
--- /issm/trunk/src/m/classes/constants.m	(revision 9859)
+++ /issm/trunk/src/m/classes/constants.m	(revision 9860)
@@ -47,7 +47,7 @@
 		function flag = checkconsistency(obj,md,solution,analyses) % {{{
 
-			fields={'g'};
-			checknan(md,'constants',fields);
-			checkgreaterstrict(md,'constants',fields,0);
+			checkfield(md,'constants.g','>',0,'size',[1 1]);
+			checkfield(md,'constants.yts','>',0,'size',[1 1]);
+			checkfield(md,'constants.referencetemperature','size',[1 1]);
 
 		end % }}}
Index: /issm/trunk/src/m/classes/diagnostic.m
===================================================================
--- /issm/trunk/src/m/classes/diagnostic.m	(revision 9859)
+++ /issm/trunk/src/m/classes/diagnostic.m	(revision 9860)
@@ -79,57 +79,50 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			fields={'restol','maxiter'};
-			checknan(md,'diagnostic',fields);
+			%Early return
+			if ~ismember(DiagnosticHorizAnalysisEnum,analyses), return; end
 
-			fields={'abstol','restol','reltol'};
-			checkgreaterstrict(md,'diagnostic',fields,0);
+			checkfield(md,'diagnostic.spcvx','forcing',1);
+			checkfield(md,'diagnostic.spcvy','forcing',1);
+			if md.mesh.dimension==3, checkfield(md,'diagnostic.spcvz','forcing',1); end
+			checkfield(md,'diagnostic.restol','size',[1 1],'>',0);
+			checkfield(md,'diagnostic.reltol','size',[1 1]);
+			checkfield(md,'diagnostic.abstol','size',[1 1]);
+			checkfield(md,'diagnostic.stokesreconditioning','size',[1 1],'NaN',1);
+			checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1);
+			if md.mesh.dimension==2,
+				checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1);
+			else
+				checkfield(md,'diagnostic.icefront','size',[NaN 6],'NaN',1);
+			end
+			checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]);
+			checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1);
+			checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]);
+			if ~isempty(md.diagnostic.requested_outputs),
+				checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]);
+			end
 
-			fields={'icefront'};
-			if (md.mesh.dimension==2),
-				checksize(md,'diagnostic',fields,[NaN 4]);
-			elseif md.mesh.dimension==3,
-				checksize(md,'diagnostic',fields,[NaN 6]);
+			%singular solution
+			if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2),
+				checkmessage(['model ' md.miscellaneous.name ' is not well posed (singular). You need at least one node with fixed velocity!'])
 			end
-			fields={'icefront(:,end)'};
-			checkvalues(md,'diagnostic',fields,[0 1 2]);
-
-			fields={'referential'};
-			checksize(md,'diagnostic',fields,[md.mesh.numberofvertices 6]);
-
-			if ~isempty(md.diagnostic.requested_outputs),
-				if(size(md.diagnostic.requested_outputs,2)~=1),
-					checkmessage(['model ' md.miscellaneous.name ' requested outputs should be a column vector']);
+			%CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
+			if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6),
+				checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. Each line of diagnostic.referential should contain either only NaN values or no NaN values']);
+			end
+			%CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
+			if any(sum(isnan(md.diagnostic.referential),2)==0),
+				pos=find(sum(isnan(md.diagnostic.referential),2)==0);
+				if any(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))),
+					dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))
+					checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. Vectors in diagnostic.referential (colums 1 to 3 and 4 to 6) must be orthogonal']);
 				end
 			end
-
-			if ismember(DiagnosticHorizAnalysisEnum,analyses),
-				%singular solution
-				if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2),
-					checkmessage(['model ' md.miscellaneous.name ' is not well posed (singular). You need at least one node with fixed velocity!'])
+			%CHECK THAT ROTATION IS IN THE (X,Y) PLANE FOR 2D MODELS
+			if md.mesh.dimension==2,
+				pos=find(sum(isnan(md.diagnostic.referential),2)==0  & md.flowequation.vertex_equation==2);
+				if any(md.diagnostic.referential(pos,3:5)~=0);
+					checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. The rotation should be in the (x,y) plane for 2D diagnostic models (nodeonmacayeal)']);
 				end
-				%ROTATED SPC
-				%CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
-				if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6),
-					checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. Each line of diagnostic.referential should contain either only NaN values or no NaN values']);
-				end
-				%CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
-				if any(sum(isnan(md.diagnostic.referential),2)==0),
-					pos=find(sum(isnan(md.diagnostic.referential),2)==0);
-					if any(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))),
-						dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))
-						checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. Vectors in diagnostic.referential (colums 1 to 3 and 4 to 6) must be orthogonal']);
-					end
-
-				end
-				%CHECK THAT ROTATION IS IN THE (X,Y) PLANE FOR 2D MODELS
-				if md.mesh.dimension==2,
-					pos=find(sum(isnan(md.diagnostic.referential),2)==0  & md.flowequation.vertex_equation==2);
-					if any(md.diagnostic.referential(pos,3:5)~=0);
-						checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. The rotation should be in the (x,y) plane for 2D diagnostic models (nodeonmacayeal)']);
-					end
-				end
-
 			end
-
 		end % }}}
 		function disp(obj) % {{{
Index: /issm/trunk/src/m/classes/flaim.m
===================================================================
--- /issm/trunk/src/m/classes/flaim.m	(revision 9859)
+++ /issm/trunk/src/m/classes/flaim.m	(revision 9860)
@@ -43,17 +43,9 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			if solution~=FlaimSolutionEnum,
-				return;
-			end
+			%Early return
+			if solution~=FlaimSolutionEnum, return; end
 
-			if ~exist(obj.tracks,'file')
-				checkmessage(['flaim.tracks file ''' md.flaim.tracks ''' does not exist.']);
-			end
-			%   probably going to need some checks on flaim.flightreqs here
-			if (numel(obj.criterion)~=md.mesh.numberofvertices) & (numel(obj.criterion)~=md.mesh.numberofelements),
-				checkmessage(['flaim.criterion vector must have number of nodes (' ...
-					int2str(md.mesh.numberofvertices) ') or elements (' int2str(md.mesh.numberofelements)...
-					') values, not ' int2str(numel(md.flaim.criterion)) ' values.']);
-			end
+			checkfield(md,'flaim.tracks','file',1);
+			checkfield(md,'flaim.criterion','numel',[md.mesh.numberofvertices md.mesh.numberofelements]);
 
 		end % }}}
Index: /issm/trunk/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk/src/m/classes/flowequation.m	(revision 9859)
+++ /issm/trunk/src/m/classes/flowequation.m	(revision 9860)
@@ -44,24 +44,14 @@
 
 			if ismember(DiagnosticHorizAnalysisEnum,analyses),
-				%Check the size of flowequation.element_equation
-				fields={'element_equation'};
-				checksize(md,'flowequation',fields,[md.mesh.numberofelements 1]);
-				%Check the values of flowequation.element_equation
-				checkvalues(md,'flowequation',fields,[0:7]);
+
+				checkfield(md,'flowequation.ismacayealpattyn','numel',1,'values',[0 1]);
+				checkfield(md,'flowequation.ishutter','numel',1,'values',[0 1]);
+				checkfield(md,'flowequation.isstokes','numel',1,'values',[0 1]);
 				if (md.mesh.dimension==2),
-					checkvalues(md,'flowequation',fields,[1 2]);
-				end
-				if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0),
-					checkmessage(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
-				end
-
-				%VERTICESTYPE
-				%Check the size of verticess_type
-				fields={'vertex_equation'};
-				checksize(md,'flowequation',fields,[md.mesh.numberofvertices 1]);
-				%Check the values of flowequation.vertex_equation
-				checkvalues(md,'flowequation',fields,[0:7]);
-				if (md.mesh.dimension==2),
-					checkvalues(md,'flowequation',fields,[1 2]);
+					checkfield(md,'flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]);
+					checkfield(md,'flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1:2]);
+				else
+					checkfield(md,'flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[0:7]);
+					checkfield(md,'flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[0:7]);
 				end
 				if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0),
Index: /issm/trunk/src/m/classes/friction.m
===================================================================
--- /issm/trunk/src/m/classes/friction.m	(revision 9859)
+++ /issm/trunk/src/m/classes/friction.m	(revision 9860)
@@ -38,13 +38,10 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			fields={'coefficient','p','q'};
-			checknan(md,'friction',fields);
-			checkgreater(md,'friction',fields,0);
+			%Early return
+			if ~ismember(DiagnosticHorizAnalysisEnum,analyses) & ~ismember(ThermalAnalysisEnum,analyses), return; end
 
-			fields={'p','q'};
-			checksize(md,'friction',fields,[md.mesh.numberofelements 1]);
-
-			fields={'coefficient'};
-			checksize(md,'friction',fields,[md.mesh.numberofvertices 1]);
+			checkfield(md,'friction.coefficient','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'friction.q','NaN',1,'size',[md.mesh.numberofelements 1]);
+			checkfield(md,'friction.p','NaN',1,'size',[md.mesh.numberofelements 1]);
 		end % }}}
 		function disp(obj) % {{{
Index: /issm/trunk/src/m/classes/geometry.m
===================================================================
--- /issm/trunk/src/m/classes/geometry.m	(revision 9859)
+++ /issm/trunk/src/m/classes/geometry.m	(revision 9860)
@@ -39,19 +39,11 @@
 		end % }}}
 		function checkconsistency(obj,md,solution,analyses) % {{{
-			fields={'surface','thickness','bed'};
-			checknan(md,'geometry',fields);
-			checksize(md,'geometry',fields,[md.mesh.numberofvertices 1]);
 
-			fields={'thickness'};
-			checkgreater(md,'geometry',fields,0);
-
+			checkfield(md,'geometry.surface'  ,'NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'geometry.bed'      ,'NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'geometry.thickness','NaN',1,'size',[md.mesh.numberofvertices 1],'>',0);
 			if any((obj.thickness-obj.surface+obj.bed)>10^-9),
 				checkmessage(['equality thickness=surface-bed violated']);
 			end 
-
-			if ismember(PrognosticAnalysisEnum,analyses),
-				fields={'thickness'};
-				checkgreaterstrict(md,'geometry',fields,0);
-			end
 		end % }}}
 		function disp(obj) % {{{
Index: /issm/trunk/src/m/classes/groundingline.m
===================================================================
--- /issm/trunk/src/m/classes/groundingline.m	(revision 9859)
+++ /issm/trunk/src/m/classes/groundingline.m	(revision 9860)
@@ -44,7 +44,6 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			if ~ismember({obj.migration},{'None' 'AgressiveMigration' 'SoftMigration'}),
-				checkmessage(['groundingline.migration field should be ''None'' ''AgressiveMigration'' or ''SoftMigration''']);
-			end
+			checkfield(md,'groundingline.migration','values',{'None' 'AgressiveMigration' 'SoftMigration'});
+
 			if ~strcmp(obj.migration,'None'),
 				if (md.mesh.dimension==3 | strcmpi(md.cluster.name,'none')),
Index: /issm/trunk/src/m/classes/mesh.m
===================================================================
--- /issm/trunk/src/m/classes/mesh.m	(revision 9859)
+++ /issm/trunk/src/m/classes/mesh.m	(revision 9860)
@@ -80,39 +80,29 @@
 		function checkconsistency(obj,md,solution,analyses) % {{{
 
-			fields={'x','y','z','elements','elementonbed','elementonsurface','vertexonbed','vertexonsurface','elementconnectivity'};
-			checknan(md,'mesh',fields);
-
-			fields={'elements','elementonbed','elementonsurface','vertexonbed','vertexonsurface','elementconnectivity'};
-			checkgreater(md,'mesh',fields,0);
-
-			fields={'dimension','elements','numberofelements','numberofvertices'};
-			checkgreaterstrict(md,'mesh',fields,0);
-			
-			fields={'vertexonbed','vertexonsurface','x','y','z'};
-			checksize(md,'mesh',fields,[md.mesh.numberofvertices 1]);
-
-			fields={'elementonbed','elementonsurface'};
-			checksize(md,'mesh',fields,[md.mesh.numberofelements 1]);
-
-			fields={'elements'};
+			checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'mesh.z','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			checkfield(md,'mesh.elements','NaN',1,'>',0,'values',1:md.mesh.numberofvertices);
 			if(md.mesh.dimension==2),
-				checksize(md,'mesh',fields,[md.mesh.numberofelements 3]);
+				checkfield(md,'mesh.elements','size',[md.mesh.numberofelements 3]);
 			else
-				checksize(md,'mesh',fields,[md.mesh.numberofelements 6]);
+				checkfield(md,'mesh.elements','size',[md.mesh.numberofelements 6]);
 			end
 			if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.mesh.elements(:)))));
-				checkmessage('orphan nodes have been found. Check the mesh');
+				checkmessage('orphan nodes have been found. Check the mesh outline');
 			end
-
-
+			checkfield(md,'mesh.dimension','values',[2 3]);
+			checkfield(md,'mesh.numberoflayers','>=',0);
+			checkfield(md,'mesh.numberofelements','>',0);
+			checkfield(md,'mesh.numberofvertices','>',0);
+			%no checks for numberofedges lat long and hemisphere
+			checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]);
+			checkfield(md,'mesh.elementonsurface','size',[md.mesh.numberofelements 1],'values',[0 1]);
+			checkfield(md,'mesh.vertexonbed','size',[md.mesh.numberofvertices 1],'values',[0 1]);
+			checkfield(md,'mesh.vertexonsurface','size',[md.mesh.numberofvertices 1],'values',[0 1]);
 			if (md.mesh.dimension==2),
-				if md.mesh.average_vertex_connectivity<9, 
-					checkmessage('model not consistent: connectivity should be at least 9 for 2d models');
-				end
-			end
-			if md.mesh.dimension==3,
-				if md.mesh.average_vertex_connectivity<24, 
-					checkmessage('model not consistent: connectivity should be at least 24 for 3d models');
-				end
+				checkfield(md,'mesh.average_vertex_connectivity','>=',9,'message','''mesh.average_vertex_connectivity'' should be at least 9 in 2d');
+			else
+				checkfield(md,'mesh.average_vertex_connectivity','>=',24,'message','''mesh.average_vertex_connectivity'' should be at least 24 in 3d');
 			end
 
@@ -121,31 +111,22 @@
 				case PrognosticSolutionEnum,
 					if md.prognostic.stabilization==3,
-						if md.mesh.dimension~=2, checkmessage('Discontinuous Galerkin only supported for 2d meshes'); end
-						fields={'edges'};
-						checksize(md,'mesh',fields,[NaN 4]);
-						fields={'edges(:,1:3)'};
-						checknan(md,'mesh',fields);
-						checkgreaterstrict(md,'mesh',fields,0);
+						checkfield(md,'mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
+						checkfield(md,'mesh.edges','size',[NaN 4]);
+						checkfield(md,'mesh.edges(:,1:3)','>',0);
 					end
 				case BalancethicknessSolutionEnum,
 					if md.balancethickness.stabilization==3,
-						if md.mesh.dimension~=2, checkmessage('Discontinuous Galerkin only supported for 2d meshes'); end
-						fields={'edges'};
-						checksize(md,'mesh',fields,[NaN 4]);
-						fields={'edges(:,1:3)'};
-						checknan(md,'mesh',fields);
-						checkgreaterstrict(md,'mesh',fields,0);
+						checkfield(md,'mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
+						checkfield(md,'mesh.edges','size',[NaN 4]);
+						checkfield(md,'mesh.edges(:,1:3)','>',0);
 					end
 				case TransientSolutionEnum,
 					if md.transient.isprognostic & md.prognostic.stabilization==3,
-						if md.mesh.dimension~=2, checkmessage('Discontinuous Galerkin only supported for 2d meshes'); end
-						fields={'edges'};
-						checksize(md,'mesh',fields,[NaN 4]);
-						fields={'edges(:,1:3)'};
-						checknan(md,'mesh',fields);
-						checkgreaterstrict(md,'mesh',fields,0);
+						checkfield(md,'mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
+						checkfield(md,'mesh.edges','size',[NaN 4]);
+						checkfield(md,'mesh.edges(:,1:3)','>',0);
 					end
 				case ThermalSolutionEnum,
-					if md.mesh.dimension~=3, checkmessage('thermal solution only supported on 3d meshes'); end
+					checkfield(md,'mesh.dimension','values',3,'message','thermal solution only supported on 3d meshes');
 			end
 		end % }}}
Index: /issm/trunk/src/m/utils/consistency/checkfield.m
===================================================================
--- /issm/trunk/src/m/utils/consistency/checkfield.m	(revision 9860)
+++ /issm/trunk/src/m/utils/consistency/checkfield.m	(revision 9860)
@@ -0,0 +1,166 @@
+function checkfield(md,fieldname,varargin)
+%CHECKFIELD - check field consistency
+%
+%   Used to check model consistency.
+%   Available options:
+%      - NaN: 1 if check that there is no NaN
+%      - size: [lines cols], NaN for non checked dimensions
+%      - >:  greater than provided value
+%      - >=: greater or equal to provided value
+%      - <:  smallerthan provided value
+%      - <=: smaller or equal to provided value
+%      - forcing: 1 if check forcing consistency (size and time)
+%      - values: cell of strings or vector of acceptable values
+%      - numel: list of acceptable number of elements
+%      - message: overloaded error message
+%
+%   Usage:
+%      checkfield(md,fieldname,options);
+%
+%   Example:
+%      checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]);
+%      checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1);
+%      checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]);
+
+%get options
+options=pairoptions(varargin{:});
+
+%get field from model
+eval(['field=md.' fieldname ';']);
+
+%Check size
+if exist(options,'size')
+	fieldsize=getfieldvalue(options,'size');
+	if isnan(fieldsize(1)),
+		if (size(field,2)~=fieldsize(2)),
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' should have ' num2str(fieldsize(2)) ' columns']));
+		end
+	elseif isnan(fieldsize(2)),
+		if (size(field,1)~=fieldsize(1)),
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' should have ' num2str(fieldsize(1)) ' lines']));
+		end
+	else
+		if ((size(field)~=fieldsize(1)) |  (size(field,2)~=fieldsize(2)))
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]));
+		end
+	end
+end
+
+%Check numel
+if exist(options,'numel')
+	fieldnumel=getfieldvalue(options,'numel');
+	if ~ismember(numel(field),fieldnumel),
+		if length(fieldnumel)==1
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' size should be ' sprintf('%g ',fieldnumel) ]));
+		elseif length(fieldnumel)==2
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' size should be ' num2str(fieldnumel(1)) ' or ' num2str(fieldnumel(2)) ]));
+		else
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' size should be ' sprintf('%g, ',fieldnumel(1:end-1)) ' or ' num2str(fieldnumel(end)) ]));
+		end
+	end
+end
+
+%check NaN
+if getfieldvalue(options,'NaN',0);
+	field2=reshape(field,prod(size(field)),1);
+	if any(isnan(field2)),
+		checkmessage(getfieldvalue(options,'message',...
+			['NaN values found in field ''' fieldname '''']));
+	end
+end
+
+%check values
+if exist(options,'values')
+	fieldvalues=getfieldvalue(options,'values');
+	if iscell(fieldvalues), %strings
+		if any(~ismember(field,fieldvalues)),
+			if length(fieldvalues)==1
+				checkmessage(getfieldvalue(options,'message',...
+					['field ''' fieldname ''' value should be ' fieldvalues{1} ]));
+			elseif length(fieldvalues)==2
+				checkmessage(getfieldvalue(options,'message',...
+					['field ''' fieldname ''' values should be ' fieldvalues{1} ' or ' fieldvalues{2} ]));
+			else
+				checkmessage(getfieldvalue(options,'message',...
+					['field ''' fieldname ''' should have values in ' sprintf('''%s'', ',fieldvalues{1:end-1}) 'or ''' fieldvalues{end} '''']));
+			end
+		end
+	else
+		field2=reshape(field,prod(size(field)),1);
+		if any(~ismember(field2,fieldvalues)),
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' should have values in [' num2str(values) ']']));
+		end
+	end
+end
+
+%check greater
+if exist(options,'>=')
+	lowerbound=getfieldvalue(options,'>=');
+	field2=reshape(field,prod(size(field)),1);
+	if any(field2<lowerbound),
+		checkmessage(getfieldvalue(options,'message',...
+			['field ''' fieldname ''' should have values above ' num2str(lowerbound)]));
+	end
+end
+if exist(options,'>')
+	lowerbound=getfieldvalue(options,'>');
+	field2=reshape(field,prod(size(field)),1);
+	if any(field2<=lowerbound),
+		checkmessage(getfieldvalue(options,'message',...
+			['field ''' fieldname ''' should have values above ' num2str(lowerbound)]));
+	end
+end
+
+%check smaller
+if exist(options,'<=')
+	upperbound=getfieldvalue(options,'<=');
+	field2=reshape(field,prod(size(field)),1);
+	if any(field2>upperbound),
+		checkmessage(getfieldvalue(options,'message',...
+			['field ''' fieldname ''' should have values above ' num2str(upperbound)]));
+	end
+end
+if exist(options,'<')
+	upperbound=getfieldvalue(options,'<');
+	field2=reshape(field,prod(size(field)),1);
+	if any(field2>=upperbound),
+		checkmessage(getfieldvalue(options,'message',...
+			['field ''' fieldname ''' should have values above ' num2str(upperbound)]));
+	end
+end
+
+%check file
+if getfieldvalue(options,'file',0),
+	if ~exist(field,'file')
+		checkmessage(['file profided in ''' fieldname ''': ''' field ''' does not exist']);
+	end
+end
+
+%Check forcings (size and times)
+if getfieldvalue(options,'forcing',0),
+	if size(field,1)==md.mesh.numberofvertices,
+		if ~size(field,2)==1,
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' should have only one column as there are md.mesh.numberofvertices lines']));
+		end
+	elseif size(field,1)==md.mesh.numberofvertices+1
+		if any(field(end,:)~=sort(field(end,:))),
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' columns should be chronological']));
+		end
+		if any(field(end,1:end-1)==field(end,2:end)),
+			checkmessage(getfieldvalue(options,'message',...
+				['field ''' fieldname ''' columns must not contain duplicate timesteps']));
+		end
+	else
+		checkmessage(getfieldvalue(options,'message',...
+			['field ''' fieldname ''' should have md.mesh.numberofvertices or md.mesh.numberofvertices+1 lines']));
+	end
+end