Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 4755)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 4756)
@@ -7,14 +7,12 @@
 %tolerance we use in litmus checks for the consistency of the model
 tolerance=10^-12;
-
-%check usage
+%check usage {{{1
 if nargin~=1,
 	help ismodelselfconsistent
 	error('ismodelselfconsistent error message: wrong usage');
 end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   TRANSIENT   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%check analysis
+%}}}
+
+%recursive call for TRANSIENT {{{1
 if (md.analysis_type==Transient2DSolutionEnum | md.analysis_type==Transient3DSolutionEnum),
 	if md.dt<=0,
@@ -34,19 +32,26 @@
 	end
 end
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMMON CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%COUNTER
+%}}}
+
+% Common checks
+%COUNTER {{{1
 if md.counter<3,
 	error(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize,setelementstype)!']);
 end
-
-%NAME
+%}}}
+%NAME{{{1
 if isempty(md.name),
 	error(['model is not correctly configured: missing name!']);
 end
-
-%ELEMENTSTYPE
+%}}}
+%ELEMENTS{{{1
+fields={'elements'};
+if (md.dim==2),
+	checksize(md,fields,[md.numberofelements 3]);
+else
+	checksize(md,fields,[md.numberofelements 6]);
+end
+%}}}
+%ELEMENTSTYPE{{{1
 if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
 	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
@@ -72,6 +77,14 @@
 	error(['model not consistent: Stokes transient not implemented yet']);
 end
-
-%ICEFRONT
+%}}}
+%DG {{{1
+if md.prognostic_DG=1;
+	%CHECK THE COMPATIBILITY OF THE EDGES
+	fields={'edges'};
+	checksize(md,fields,[NaN 4]);
+	checknan(md,fields);
+end
+%}}}
+%PRESSURELOAD{{{1
 if (md.dim==2),
 	fields={'pressureload'};
@@ -83,49 +96,41 @@
 	error('type should be either ''2d'' or ''3d''');
 end
-
-%ELEMENTS
-fields={'elements'};
-if (md.dim==2),
-	checksize(md,fields,[md.numberofelements 3]);
-else
-	checksize(md,fields,[md.numberofelements 6]);
-end
-
-%NO NAN
+%}}}
+%NO NAN {{{1
 fields={'numberofelements','numberofgrids','x','y','z','drag_coefficient','drag_type','drag_p','drag_q',...
 	'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
 	'tolx','np','eps_res','max_nonlinear_iterations','exclusive','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 checknan(md,fields);
-
-%FIELDS >= 0 
+%}}}}
+%FIELDS >= 0 {{{1
 fields={'numberofelements','numberofgrids','elements','drag_coefficient','drag_type','drag_p','drag_q',...
 	'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
 	'sparsity','lowmem','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 checkgreater(md,fields,0);
-
-%FIELDS > 0
+%}}}
+%FIELDS > 0 {{{1
 fields={'numberofelements','numberofgrids','elements','drag_type','drag_p',...
 	'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter','tolx',...
 	'sparsity','deltaH','DeltaH','timeacc','timedec'};
 checkgreaterstrict(md,fields,0);
-
-%SIZE NUMBEROFELEMENTS
+%}}}
+%SIZE NUMBEROFELEMENTS {{{1
 fields={'drag_p','drag_q','elementoniceshelf','rheology_n','elementonbed'};
 checksize(md,fields,[md.numberofelements 1]);
-
-%SIZE NUMBEROFGRIDS
+%}}}
+%SIZE NUMBEROFGRIDS {{{1
 fields={'x','y','z','rheology_B','drag_coefficient','melting_rate','accumulation_rate','surface','thickness','bed','gridonbed','gridonsurface'};
 checksize(md,fields,[md.numberofgrids 1]);
-
-%SIZE 6
+%}}}
+%OTHER SIZES {{{1
 fields={'spcvelocity'};
 checksize(md,fields,[md.numberofgrids 6]);
-
-%THICKNESS = SURFACE - BED
+%}}}
+%THICKNESS = SURFACE - BED {{{1
 if any((md.thickness-md.surface+md.bed)>tolerance),
 	error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
 end
-
-%RIFTS
+%}}}
+%RIFTS{{{1
 if md.numrifts,
 	if ~(md.dim==2),
@@ -146,16 +151,22 @@
 	end
 end
-
-%ARTIFICIAL DIFFUSIVITY
+%}}}
+%FLAGS (0 or 1){{{1
 if ~ismember(md.artificial_diffusivity,[0 1]),
 	error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
 end
-
-%PARAMETEROUTPUT
+if ~ismember(md.prognostic_DG,[0 1]),
+	error('model not consistent: prognostic_DG should be a scalar (1 or 0)');
+end
+if ~ismember(md.lowmem,[0 1]),
+	error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
+end
+%}}}
+%PARAMETEROUTPUT {{{1
 if md.numoutput~=length(md.parameteroutput),
 	error('model not consistent: numoutput should be the same size as parameteroutput');
 end
-
-%CONNECTIVITY
+%}}}
+%CONNECTIVITY {{{1
 if (md.dim==2),
 	if md.connectivity<9, 
@@ -168,11 +179,6 @@
 	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
+%}}}
+%PARALLEL{{{1
 if ~strcmpi(md.cluster,'none'),
 
@@ -186,41 +192,8 @@
 
 end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SOLUTION CHECKS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%QMU
-if md.qmu_analysis,
-	if md.qmu_params.evaluation_concurrency~=1,
-		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,
-			error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
-		end
-		if find(md.part)>=md.numberofgrids,
-			error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
-		end
-		if min(md.part)~=0,
-			error(['model not consistent: partition vector not indexed from 0 on']);
-		end
-		if max(md.part)>=md.numberofgrids,
-			error(['model not consistent: partition vector cannot have maximum index larger than number of grids']);
-		end
-		if ~isempty(find(md.part<0)),
-			error(['model not consistent: partition vector cannot have values less than 0']);
-		end
-		if ~isempty(find(md.part>=md.npart)),
-			error(['model not consistent: partition vector cannot have values more than md.npart-1']);
-		end
-		if max(md.part)>=md.npart,
-			error(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']);
-		end
-	end
-	if md.eps_rel>1.1*10^-5,
-		error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
-	end
-end
-
-%DIAGNOSTIC
+%}}}
+
+% Solution checks
+%DIAGNOSTIC{{{1
 if md.analysis_type==DiagnosticSolutionEnum,
 
@@ -259,6 +232,6 @@
 	end
 end
-
-%PROGNOSTIC
+%}}}
+%PROGNOSTIC{{{1
 if md.analysis_type==PrognosticSolutionEnum,
 
@@ -276,6 +249,6 @@
 
 end
-
-%STEADYSTATE
+%}}}
+%STEADYSTATE{{{1
 if md.analysis_type==SteadystateSolutionEnum,
 
@@ -300,5 +273,6 @@
 	end
 end
-
+%}}}
+%THERMAL {{{1
 %THERMAL STEADY AND THERMAL TRANSIENT
 if md.analysis_type==ThermalSolutionEnum,
@@ -338,10 +312,10 @@
 
 end
-
-%BALANCEDTHICKNESS
+%}}}
+%BALANCEDTHICKNESS{{{1
 if md.analysis_type==BalancedthicknessSolutionEnum
 
 	%VELOCITIES MELTING AND ACCUMULATION
-	fields={'vx','vy','accumulation_rate','melting_rate'};
+	fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
 	checksize(md,fields,[md.numberofgrids 1]);
 	checknan(md,fields);
@@ -352,20 +326,6 @@
 	end 
 end
-
-%BALANCEDTHICKNESS2
-if md.analysis_type==Balancedthickness2SolutionEnum
-
-	%VELOCITIES MELTING AND ACCUMULATION
-	fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
-	checksize(md,fields,[md.numberofgrids 1]);
-	checknan(md,fields);
-
-	 %CHECK THE COMPATIBILITY OF THE EDGES
-	 fields={'edges'};
-	 checksize(md,fields,[NaN 4]);
-	 checknan(md,fields);
-end
-
-%BALANCEDVELOCITIES
+%}}}
+%BALANCEDVELOCITIES{{{1
 if md.analysis_type==BalancedvelocitiesSolutionEnum
 
@@ -380,6 +340,6 @@
 	end
 end
-
-%CONTROL
+%}}}
+%CONTROL{{{1
 if md.control_analysis,
 
@@ -418,6 +378,38 @@
 	checknan(md,fields);
 end
-
-%QMU
+%}}}
+%QMU {{{1
+if md.qmu_analysis,
+	if md.qmu_params.evaluation_concurrency~=1,
+		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,
+			error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
+		end
+		if find(md.part)>=md.numberofgrids,
+			error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
+		end
+		if min(md.part)~=0,
+			error(['model not consistent: partition vector not indexed from 0 on']);
+		end
+		if max(md.part)>=md.numberofgrids,
+			error(['model not consistent: partition vector cannot have maximum index larger than number of grids']);
+		end
+		if ~isempty(find(md.part<0)),
+			error(['model not consistent: partition vector cannot have values less than 0']);
+		end
+		if ~isempty(find(md.part>=md.npart)),
+			error(['model not consistent: partition vector cannot have values more than md.npart-1']);
+		end
+		if max(md.part)>=md.npart,
+			error(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']);
+		end
+	end
+	if md.eps_rel>1.1*10^-5,
+		error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
+	end
+end
+
 if strcmpi(md.analysis_type,'qmu'),
 	if ~strcmpi(md.cluster,'none'),
@@ -427,9 +419,10 @@
 	end
 end
-
-end %end function
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+%}}}
+
+end
+
+%checks additional functions
+%checklength {{{1
 function checklength(md,fields,fieldlength)
 	%CHECKSIZE - check length of a field
@@ -440,5 +433,6 @@
 	end
 end
-
+%}}}
+%checksize {{{1
 function checksize(md,fields,fieldsize)
 	%CHECKSIZE - check size of a field
@@ -475,5 +469,6 @@
 	end
 end
-
+%}}}
+%checknan {{{1
 function checknan(md,fields)
 	%CHECKNAN - check nan values of a field
@@ -484,5 +479,6 @@
 	end
 end
-
+%}}}
+%checkreal{{{1
 function checkreal(md,fields)
 	%CHECKREAL - check real values of a field
@@ -493,5 +489,6 @@
 	end
 end
-
+%}}}
+%checkgreaterstrict{{{1
 function checkgreaterstrict(md,fields,lowerbound)
 	%CHECKGREATERSTRICT - check values of a field
@@ -502,5 +499,6 @@
 	end
 end
-
+%}}}
+%checkgreater{{{1
 function checkgreater(md,fields,lowerbound)
 	%CHECKGREATER - check values of a field
@@ -511,5 +509,6 @@
 	end
 end
-
+%}}}
+%checklessstrict{{{1
 function checklessstrict(md,fields,upperbound)
 	%CHECKLESSSTRICT - check values of a field
@@ -520,5 +519,6 @@
 	end
 end
-
+%}}}
+%checkless{{{1
 function checkless(md,fields,upperbound)
 	%CHECKLESS - check values of a field
@@ -529,5 +529,6 @@
 	end
 end
-
+%}}}
+%checkvalues {{{1
 function checkvalues(md,fields,values)
 	%CHECKVALUE - check that a field has a certain value
@@ -538,2 +539,3 @@
 	end
 end
+%}}}
