function bool=ismodelselfconsistent(md,solutiontype,package),
%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
%
%   Usage:
%      bool=ismodelselfconsistent(md,solutiontype,package),

if strcmpi(solutiontype,'mesh'),
	%this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!

	bool=1;
	return;
end

%tolerance we use in litmus tests for the consistency of the model
tolerance=10^-12;
if (nargin~=3  )
	IsModelSelfConsistentUsage();
	error(' ');
end

if strcmpi(solutiontype,'qmu'),
	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

if isempty(md.name),
	disp(['model is not correctly configured: missing name!']);
	bool=0;return;
end

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;
end

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;
end

%check 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;
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum)  + (md.elements_type(:,1)==pattynenum)))
	disp(['Types of elements have not been set properly, run setelementstype first'])
	bool=0;return;
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum)))
	disp(['Types of elements have not been set properly, run setelementstype first'])
	bool=0;return;
end
if strcmpi(md.type,'2d'),
	if (ismember(pattynenum,md.elements_type(:,1)) |  ismember(stokesenum,md.elements_type(:,2))),
		disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
		bool=0;return;
	end
end
%check that if there are rifts, model is 2d
if md.numrifts,
	if ~strcmpi(md.type,'2d'),
		disp(['Models with rifts are only supported in 2d for now!']);
		bool=0;return;
	end
end

%If there are rifts, check that the solver is set to mumps in parallel
if md.numrifts, 
	if ~strcmpi(md.cluster,'none')
		if isempty(findstr('aijmumps',md.solverstring)),
			disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
			bool=0;return;
		end
	end
end

%check elementstype in diagnostic
if strcmpi(solutiontype,'diagnostic')
	if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
		disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
	end
end

%Check that mesh is 3d
if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),
	if strcmp(md.type,'2d'),
		disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
		bool=0;return;
	end
end

%check artificial_diffusivity is an integer
if ~isscalar(md.artificial_diffusivity),
	disp('artificial_diffusivity should be a scalar (1 or 0)');
	bool=0;return;
end

%check on 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





%Check previous computation
if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
	if (isempty(md.vx) | isempty(md.vy) | isempty(md.vz)),
		disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
		bool=0;return;
	end
end

if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
	if (isempty(md.pressure) | isnan(md.pressure)),
		disp(['pressure is required. Run ''diagnostic'' solution first!'])
		bool=0;return;
	end
	
	if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy)  | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
		disp(['velocities are required!']);
		bool=0;return;
	end
end

if strcmp(solutiontype,'thermaltransient')
	if isempty(md.temperature),
		disp(['An  initial temperature is needed for a transient thermal computation'])
		bool=0;return;
	end
	if isstruct(md.temperature) |  isstruct(md.melting),
		disp(['The initial temperature or melting should be a list and not a structure'])
		bool=0;return;
	end

end

if strcmp(solutiontype,'transient')
	if isstruct(md.temperature),
		disp(['The initial temperature should be empty or a list but not a structure'])
		bool=0;return;
	end

end

if strcmp(solutiontype,'parameters')
	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
	if strcmpi(md.type,'2d') & md.acceleration
		disp(['solution requires acceleration=0']);
		bool=0;return;
	end
	if any(md.elements_type(:,1)==hutterenum); 
		disp(['The model has Hutter''s elements. Impossible to compute parameters']);
		bool=0;return;
	end
end


%Check that problem is not singular
if strcmp(solutiontype,'diagnostic'),
	if isempty(find(md.gridondirichlet_diag)),
		disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
		bool=0;return;
	end
end

%Check we have values for dirichletvalues_diag
if strcmp(solutiontype,'diagnostic'),
	if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
		disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
		bool=0;return;
	end
end



if strcmp(solutiontype,'prognostic'),
	%old testing
%	if isempty(find(md.gridondirichlet_diag)),
%		disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
%		bool=0;return;
%	end
%	if isempty(find(md.gridondirichlet_prog)),
%		disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
%		bool=0;return;
%	end
	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;
	end
end
%some checks for control methods
if strcmp(solutiontype,'control')
	%check fields used by control
	if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
		disp('maxiter, optscal and fit must have the length specified by nsteps')
		bool=0;return;
	end
	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;
	end
end


%Check that no NaN values popped up: 
fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
        'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
		'tolx','np','time','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)) '!']);
			bool=0; return;
		end
	end
end

%Check that these fields are > 0 
fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
         'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','np','time','exclusive',...
		 '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)) '!']);
			bool=0; return;
		end
	end
end

%Check that these fields are ~=0 
fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
         'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx','np','time',...
		 '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)) '!']);
			bool=0; return;
		end
	end
end

%Check that these fields are of size md.numberofelements
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) '!']);
		bool=0; return;
	end
end

%Check that these fields are of size md.numberofgrids
fields={'melting','x','y','z','B','drag','gridondirichlet_diag','vx_obs','vy_obs','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) '!']);
		bool=0; return;
	end
end

%Make sure lowmem is either 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

%Make sure sparsity is 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

%Make sure thickness=surface-bed;
if find((md.thickness-md.surface+md.bed)>tolerance),
	disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
	bool=0; return;
end

%Check that rifts were processed correctly: 
if ~isstruct(md.rifts),
	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;
	end
end


%Check that thickness is not null. If is it, then the corresponding grids should be spc'd.
if ~isempty(find(md.thickness<=0)),
	pos=find(md.thickness<=0);
	if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
		disp(['model ' md.name ' has some grids with 0 thickness']);
		bool=0; return;
	end
end

%Check that p>0 in the sliding law u = k * sigma ^ p * Neff ^ q :
if ~isempty(find(md.p<=0)),
	disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
	bool=0; return;
end


%Check parameteroutput
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

%check parallel settings
if ~strcmpi(md.cluster,'none'),
	if ~strcmpi(package,'cielo'),
		disp(['parallel runs only allowed for ''cielo'' package']);
		bool=0; return;
	end
end

%Put as many checks as you want here: 


%No problems, just return 1;
bool=1;
return;

function IsModelSelfConsistentUsage()
disp(' ');
disp('   IsModelSelfConsistent usage: flag=IsModelSelfConsistent(model)');
disp('         where model is an instance of the @model class, and flag a boolean');




