%MODEL class definition
%
%   Usage:
%      md = model(varargin)

classdef model
    properties (SetAccess=public) %Model fields
		 % {{{
		 %Careful here: no other class should be used as default value this is a bug of matlab
		 mesh             = 0;
		 mask             = 0;

		 geometry         = 0;
		 constants        = 0;
		 surfaceforcings  = 0;
		 basalforcings    = 0;
		 materials        = 0;
		 friction         = 0;
		 flowequation     = 0;
		 timestepping     = 0;
		 initialization   = 0;
		 rifts            = 0;

		 debug            = 0;
		 verbose          = 0;
		 settings         = 0;
		 solver           = 0;
		 cluster          = 0;

		 balancethickness = 0;
		 diagnostic       = 0;
		 groundingline    = 0;
		 hydrology        = 0;
		 prognostic       = 0;
		 thermal          = 0;
		 steadystate      = 0;
		 transient        = 0;

		 autodiff         = 0;
		 flaim            = 0;
		 inversion        = 0;
		 qmu              = 0;

		 results          = 0;
		 radaroverlay     = 0;
		 miscellaneous    = 0;
		 private          = 0;

		 %}}}
	 end
	 methods (Static)
		 function md = loadobj(md) % {{{
			 % This function is directly called by matlab when a model object is
			 % loaded. If the input is a struct it is an old version of model and
			 % old fields must be recovered (make sure they are in the deprecated
			 % model properties)

			 if verLessThan('matlab','7.9'),
				 disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
				 disp('         if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');

				 % This is a Matlab bug: all the fields of md have their default value
				 % Example of error message:
				 % Warning: Error loading an object of class 'model':
				 % Undefined function or method 'exist' for input arguments of type 'cell'
				 %
				 % This has been fixed in MATLAB 7.9 (R2009b) and later versions
			 end

			 if isstruct(md)
				 disp('Recovering model object from a previous version');
				 md = structtomodel(model,md);
			 end
		 end% }}}
	 end
	 methods
		 function md = model(varargin) % {{{

			 switch nargin
				 case 0
					 md=setdefaultparameters(md);
				 otherwise
					 error('model constructor error message: 0 of 1 argument only in input.');
				 end
		 end
		 %}}}
		 function md = checkmessage(md,string) % {{{
			 if(nargout~=1) error('wrong usage, model must be an output'); end
			 disp(['model not consistent: ' string]);
			 md.private.isconsistent=false;
		 end
		 %}}}
		 function md = structtomodel(md,structmd) % {{{

			 if ~isstruct(structmd) error('input model is not a structure'); end

			 %loaded model is a struct, initialize output and recover all fields
			 md = structtoobj(model,structmd);

			 %Old field now classes
			 if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end
			 if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end

			 %Field name change
			 if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
			 if isfield(structmd,'p'), md.friction.p=structmd.p; end
			 if isfield(structmd,'q'), md.friction.q=structmd.p; end
			 if isfield(structmd,'melting'), md.basalforcings.melting_rate=structmd.melting; end
			 if isfield(structmd,'melting_rate'), md.basalforcings.melting_rate=structmd.melting_rate; end
			 if isfield(structmd,'accumulation'), md.surfaceforcings.mass_balance=structmd.accumulation; end
			 if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end
			 if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end
			 if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end
			 if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end
			 if isfield(structmd,'gridonbed'), md.mesh.vertexonbed=structmd.gridonbed; end
			 if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end
			 if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end
			 if isfield(structmd,'gridoniceshelf'), md.mask.vertexonfloatingice=structmd.gridoniceshelf; end
			 if isfield(structmd,'gridonicesheet'), md.mask.vertexongroundedice=structmd.gridonicesheet; end
			 if isfield(structmd,'gridonwater'), md.mask.vertexonwater=structmd.gridonwater; end
			 if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end
			 if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.solver=structmd.petscoptions; end
			 if isfield(structmd,'g'), md.constants.g=structmd.g; end
			 if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end
			 if isfield(structmd,'surface_mass_balance'), md.surfaceforcings.mass_balance=structmd.surface_mass_balance; end
			 if isfield(structmd,'basal_melting_rate'), md.basalforcings.melting_rate=structmd.basal_melting_rate; end
			 if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end
			 if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end
			 if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
			 if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end
			 if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end
			 if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end
			 if isfield(structmd,'riftproperties'), %old implementation
				 md.rifts=rifts();
				 md.rifts.riftproperties=structmd.riftproperties; 
				 md.rifts.riftstruct=structmd.rifts;
				 md.rifts.riftproperties=structmd.riftinfo;
			 end
			 if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end
			 if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end
			 if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end
			 if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end
			 if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end
			 if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end
			 if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end
			 if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end
			 if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end
			 if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end
			 if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end
			 if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end
			 if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end
			 if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end
			 if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end
			 if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end
			 if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end
			 if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end
			 if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end
			 if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end
			 if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end
			 if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end
			 if isfield(structmd,'spcthickness'), md.prognostic.spcthickness=structmd.spcthickness; end
			 if isfield(structmd,'artificial_diffusivity'), md.prognostic.stabilization=structmd.artificial_diffusivity; end
			 if isfield(structmd,'hydrostatic_adjustment'), md.prognostic.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end
			 if isfield(structmd,'penalties'), md.prognostic.vertex_pairing=structmd.penalties; end
			 if isfield(structmd,'penalty_offset'), md.prognostic.penalty_factor=structmd.penalty_offset; end
			 if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
			 if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
			 if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
			 if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
			 if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end
			 if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end
			 if isfield(structmd,'elementonwater'), md.mask.elementonwater=structmd.elementonwater; end
			 if isfield(structmd,'nodeoniceshelf'), md.mask.vertexonfloatingice=structmd.nodeoniceshelf; end
			 if isfield(structmd,'nodeonicesheet'), md.mask.vertexongroundedice=structmd.nodeonicesheet; end
			 if isfield(structmd,'nodeonwater'), md.mask.vertexonwater=structmd.nodeonwater; end
			 if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end
			 if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end
			 if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end
			 if isfield(structmd,'ismacayealpattyn'), md.flowequation.ismacayealpattyn=structmd.ismacayealpattyn; end
			 if isfield(structmd,'ishutter'), md.flowequation.ishutter=structmd.ishutter; end
			 if isfield(structmd,'isstokes'), md.flowequation.isstokes=structmd.isstokes; end
			 if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end
			 if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end
			 if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end
			 if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end
			 if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end
			 if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end
			 if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end
			 if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end
			 if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end
			 if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end
			 if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end
			 if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end
			 if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end
			 if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end
			 if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end
			 if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end
			 if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end
			 if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end
			 if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end
			 if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end
			 if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end
			 if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end
			 if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end
			 if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end
			 if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end
			 if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end
			 if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end
			 if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end
			 if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end
			 if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end
			 if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end
			 if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end
			 if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end
			 if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end
			 if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end
			 if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end
			 if isfield(structmd,'long'), md.mesh.long=structmd.long; end
			 if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end
			 if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end
			 if isfield(structmd,'dim'), md.mesh.dimension=structmd.dim; end
			 if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end
			 if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end
			 if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end
			 if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end
			 if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end
			 if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end
			 if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end
			 if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end
			 if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end
			 if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end
			 if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end
			 if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end
			 if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end
			 if isfield(structmd,'elementonbed'), md.mesh.elementonbed=structmd.elementonbed; end
			 if isfield(structmd,'elementonsurface'), md.mesh.elementonsurface=structmd.elementonsurface; end
			 if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end
			 if isfield(structmd,'nodeonbed'), md.mesh.vertexonbed=structmd.nodeonbed; end
			 if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end
			 if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end
			 if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end
			 if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end
			 if isfield(structmd,'edges'), md.mesh.edges=structmd.edges; end
			 if isfield(structmd,'y'), md.mesh.y=structmd.y; end
			 if isfield(structmd,'x'), md.mesh.x=structmd.x; end
			 if isfield(structmd,'z'), md.mesh.z=structmd.z; end
			 if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end
			 if isfield(structmd,'pressureload'), md.diagnostic.icefront=structmd.pressureload; end
			 if isfield(structmd,'diagnostic_ref'), md.diagnostic.referential=structmd.diagnostic_ref; end
			 if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end
			 if isfield(structmd,'part'); md.qmu.partition=structmd.part; end
				 
			 %Field changes
			 if (isfield(structmd,'type') & ischar(structmd.type)), 
				 if strcmpi(structmd.type,'2d'), md.mesh.dimension=2; end
				 if strcmpi(structmd.type,'3d'), md.mesh.dimension=3; end
			 end
			 if isnumeric(md.verbose),
				 md.verbose=verbose;
			 end
			 if size(md.diagnostic.icefront,2)==3 || size(md.diagnostic.icefront,2)==5,
				 front=md.diagnostic.icefront;
				 md.diagnostic.icefront=[front 1*md.mask.elementonfloatingice(front(:,end))];
			 end
			 if isfield(structmd,'spcvelocity'), 
				 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
				 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
				 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
				 pos=find(structmd.spcvelocity(:,1)); md.diagnostic.spcvx(pos)=structmd.spcvelocity(pos,4); 
				 pos=find(structmd.spcvelocity(:,2)); md.diagnostic.spcvy(pos)=structmd.spcvelocity(pos,5); 
				 pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); 
			 end
			 if isfield(structmd,'spcvx'), 
				 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
				 pos=find(~isnan(structmd.spcvx)); md.diagnostic.spcvx(pos)=structmd.spcvx(pos); 
			 end
			 if isfield(structmd,'spcvy'),
				 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
				 pos=find(~isnan(structmd.spcvy)); md.diagnostic.spcvy(pos)=structmd.spcvy(pos);     
			 end
			 if isfield(structmd,'spcvz'),
				 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
				 pos=find(~isnan(structmd.spcvz)); md.diagnostic.spcvz(pos)=structmd.spcvz(pos);     
			 end
			 if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]),
				 pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,end)=0;
				 pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,end)=1;
				 pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,end)=2;
			 end
			 if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50,
				 pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0;
				 pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1;
				 pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2;
				 pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3;
				 pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4;
				 pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5;
				 pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6;
				 pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7;
			 end
			 if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50,
				 pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0;
				 pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1;
				 pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2;
				 pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3;
				 pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4;
				 pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5;
				 pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6;
				 pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7;
			 end
			 if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law),
				 if (structmd.rheology_law==272), md.materials.rheology_law='None';      end
				 if (structmd.rheology_law==368), md.materials.rheology_law='Paterson';  end
				 if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end
			 end
			 if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration),
				 if (structmd.groundingline_migration==272), md.groundingline.migration='None';      end
				 if (structmd.groundingline_migration==273), md.groundingline.migration='AgressiveMigration';  end
				 if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end
			 end
			 if isfield(structmd,'control_type') & isnumeric(structmd.control_type),
				 if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end
				 if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end
				 if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end
			 end
			 if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]),
				 pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101;
				 pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102;
				 pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103;
				 pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104;
				 pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105;
				 pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201;
				 pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501;
				 pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502;
				 pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503;
			 end

			 if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2,
					 md.thermal.stabilization=2;
					 md.prognostic.stabilization=1;
					 md.balancethickness.stabilization=1;
			 end
			 if isnumeric(md.prognostic.hydrostatic_adjustment)
				 if md.prognostic.hydrostatic_adjustment==269,
					 md.prognostic.hydrostatic_adjustment='Incremental';
				 else
					 md.prognostic.hydrostatic_adjustment='Absolute';
				 end
			 end

			 %New fields
			 if ~isfield(structmd,'upperelements');
				 md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d;
				 md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN;
			 end
			 if ~isfield(structmd,'lowerelements');
				 md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d;
				 md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN;
			 end
				 
			 if ~isfield(structmd,'diagnostic_ref');
				 md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
			 end

		 end% }}}
		 function md = setdefaultparameters(md) % {{{

			 %initialize subclasses
			 md.mesh             = mesh();
			 md.mask             = mask();
			 md.constants        = constants();
			 md.geometry         = geometry();
			 md.initialization   = initialization();
			 md.surfaceforcings  = surfaceforcings();
			 md.basalforcings    = basalforcings();
			 md.friction         = friction();
			 md.rifts            = rifts();
			 md.timestepping     = timestepping();
			 md.groundingline    = groundingline();
			 md.materials        = materials();
			 md.flowequation     = flowequation();
			 md.debug            = debug();
			 md.verbose          = verbose('solution',true,'qmu',true,'control',true);
			 md.settings         = settings();
			 md.solver           = solver();
			 if ismumps,
				 md.solver           = addoptions(md.solver,DiagnosticVertAnalysisEnum,mumpsoptions);
			 else
				 md.solver           = addoptions(md.solver,DiagnosticVertAnalysisEnum,iluasmoptions);
			 end
			 md.cluster          = generic();
			 md.balancethickness = balancethickness();
			 md.diagnostic       = diagnostic();
			 md.hydrology        = hydrology();
			 md.prognostic       = prognostic();
			 md.thermal          = thermal();
			 md.steadystate      = steadystate();
			 md.transient        = transient();
			 md.autodiff         = autodiff();
			 md.flaim            = flaim();
			 md.inversion        = inversion();
			 md.qmu              = qmu();
			 md.radaroverlay     = radaroverlay();
			 md.results          = struct();
			 md.miscellaneous    = miscellaneous();
			 md.private          = private();
		 end
		 %}}}
		 function disp(obj) % {{{
			 disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(obj.mesh) ']'],'mesh properties'));
			 disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(obj.mask) ']'],'defines grounded and floating elements'));
			 disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(obj.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
			 disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(obj.constants) ']'],'physical constants'));
			 disp(sprintf('%19s: %-22s -- %s','surfaceforcings' ,['[1x1 ' class(obj.surfaceforcings) ']'],'surface forcings'));
			 disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings'));
			 disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(obj.materials) ']'],'material properties'));
			 disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties'));
			 disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(obj.flowequation) ']'],'flow equations'));
			 disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(obj.timestepping) ']'],'time stepping for transient models'));
			 disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(obj.initialization) ']'],'initial guess/state'));
			 disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(obj.rifts) ']'],'rifts properties'));
			 disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(obj.debug) ']'],'debugging tools (valgrind, gprof)'));
			 disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve'));
			 disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(obj.settings) ']'],'settings properties'));
			 disp(sprintf('%19s: %-22s -- %s','solver'          ,['[1x1 ' class(obj.solver) ']'],'PETSc options for each solution'));
			 disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)'));
			 disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution'));
			 disp(sprintf('%19s: %-22s -- %s','diagnostic'      ,['[1x1 ' class(obj.diagnostic) ']'],'parameters for diagnostic solution'));
			 disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution'));
			 disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution'));
			 disp(sprintf('%19s: %-22s -- %s','prognostic'      ,['[1x1 ' class(obj.prognostic) ']'],'parameters for prognostic solution'));
			 disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution'));
			 disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution'));
			 disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution'));
			 disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters'));
			 disp(sprintf('%19s: %-22s -- %s','flaim'           ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters'));
			 disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(obj.inversion) ']'],'parameters for inverse methods'));
			 disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(obj.qmu) ']'],'dakota properties'));
			 disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(obj.results) ']'],'model results'));
			 disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay'));
			 disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields'));
		 end % }}}
	 end
 end
