function md = model(varargin)
%MODEL - constructor for a model
%
%   Usage:
%      md = model(varargin)

switch nargin

case 0
	% if no input arguments, create a default object

	%Model general information
	md.notes='';
	md.name='';
	md.counter=0; %1 mesh, 2 geography, 3 parameterize
	md.domainoutline=NaN;
	md.riftoutline=NaN;
	md.iceshelfoutline=NaN;
	md.icesheetoutline=NaN;
	md.parameterfile=NaN;
	
	%Mesh
	md.type='2d';
	md.numberofelements=0;
	md.numberofgrids=0;
	md.elements=NaN;
	md.elements_type=NaN;
	md.x=NaN;
	md.y=NaN;
	md.z=NaN;
	md.bed_slopex=NaN;
	md.bed_slopey=NaN;
	md.surface_slopex=NaN;
	md.surface_slopey=NaN;

	%Initial 2d mesh 
	md.numberofelements2d=0;
	md.numberofgrids2d=0;
	md.elements2d=NaN;
	md.elements_type2d=NaN;
	md.x2d=NaN;
	md.y2d=NaN;
	md.z2d=NaN;
		
	%Elements
	md.elementonhutter=NaN;
	md.elementonmacayeal=NaN;
	md.elementonpattyn=NaN;
	md.elementonstokes=NaN;

	%Nodes
	md.gridonhutter=NaN;
	md.gridonmacayeal=NaN;
	md.gridonpattyn=NaN;
	md.gridonstokes=NaN;

	%Stokes mesh
	md.stokesreconditioning=10^13;
	
	%Penalties
	md.penalties=NaN;
	md.penalty_offset=3;   %penalty=Kmax*10^penalty_offset
	md.penalty_melting=10^7;%penalty used to keep T<Tpmp
	md.penalty_lock=10;
	md.segments=NaN;       %private member, not displayed
	md.segmentmarkers=NaN; %private member, not displayed
	md.rifts=NaN;
	md.numrifts=0;

	%Projections
	md.uppergrids=NaN;
	md.lowergrids=NaN;
	md.deadgrids=NaN;
	
	%Extrusion
	md.numlayers=0;
	md.elementonbed=NaN;
	md.elementonsurface=NaN;
	md.gridonbed=NaN;
	md.gridonsurface=NaN;
	md.minh=1;
	md.firn_layer=NaN;
	md.dof=NaN;

	%Extraction
	md.extractedgrids=NaN;
	md.extractedelements=NaN;
	
	%Materials parameters
	md.rho_ice=0;
	md.rho_water=0;
	md.heatcapacity=2093;
	md.latentheat=3.34*10^5; %(J/kg);
	md.thermalconductivity=2.4; %(W/m/K)
	md.meltingpoint=273.15; %The melting point of ice at 1 atmosphere of pressure in K
	md.beta=9.8*10^-8; %dependence of melting on pressure.
	md.mixed_layer_capacity=3974; %J/kg/K
	md.thermal_exchange_velocity=1.00*10^-4; %m/s
	
	%Physical parameters
	md.g=0;
	md.yts=365*24*3600;
	md.drag=NaN;
	md.drag_type=2;
	md.p=NaN;
	md.q=NaN;
	md.B=NaN;
	md.n=NaN;

	%Geometrical parameters
	md.elementoniceshelf=NaN;
	md.elementonicesheet=NaN;
	md.gridoniceshelf=NaN;
	md.gridonicesheet=NaN;
	md.surface=NaN;
	md.thickness=NaN;
	md.new_thickness=NaN;
	md.bed=NaN;

	%Boundary conditions
	md.gridonboundary=NaN;
	%Diagnostic
	md.segmentonneumann_diag=NaN;
	md.neumannvalues_diag=NaN;
	md.gridondirichlet_diag=NaN;
	md.dirichletvalues_diag=NaN;

	%Thermal
	md.min_thermal_constraints=0; %minimum number of unstable thermal constraints
	md.gridondirichlet_thermal=NaN;
	md.dirichletvalues_thermal=NaN;

	%Transient
	md.segmentonneumann_prog=NaN;
	md.neumannvalues_prog=NaN;
	md.segmentonneumann_prog2=NaN;
	md.neumannvalues_prog2=NaN;
	md.gridondirichlet_prog=NaN;
	md.dirichletvalues_prog=NaN;

	%Observations
	md.vx_obs=NaN;
	md.vy_obs=NaN;
	md.vel_obs=NaN;
	md.accumulation=NaN;
	md.geothermalflux=NaN;
	md.observed_temperature=NaN;

	%Statics parameters
	md.eps_rel=0;
	md.eps_abs=0;
	md.acceleration=0;
	md.sparsity=0;
	md.connectivity=10;
	md.lowmem=0;
	md.viscosity_overshoot=0.3;
	
	%Transient parameters
	md.dt=0;
	md.ndt=0;
	md.artificial_diffusivity=0;
	md.transient_results=NaN;
	md.thermaltransient_results=NaN;
	md.timestepping=0; %no adaptative time stepping by default.
	md.deltaH=10^-2; %meters
	md.DeltaH=50; %meters
	md.deltaT=10^-2; %K
	md.DeltaT=1; %K
	md.timeacc=2;
	md.timedec=.5;


	%Control
	md.control_type='drag';
	md.cont_vx=NaN;
	md.cont_vy=NaN;
	md.cont_vz=NaN;
	md.cont_vel=NaN;
	md.cont_parameters=NaN;
	md.cont_J=[];
	md.nsteps=0;
	md.maxiter=[];
	md.tolx=0;
	md.optscal=[];
	md.mincontrolconstraint=0;
	md.maxcontrolconstraint=2000;
	md.fit=[];
	md.epsvel=eps;
	md.meanvel=0;
	
	%Output parameters
	md.parameteroutput={};
	md.viscousheating=NaN;
	md.pressure_elem=NaN;
	md.stress=NaN;
	md.stress_surface=NaN;
	md.stress_bed=NaN;
	md.deviatoricstress=NaN;
	md.strainrate=NaN;
	md.plot=1;
	
	%debugging
	md.debug=1;
	md.element_debug=0;
	md.element_debugid=NaN;

	%Results fields
	md.vx=NaN;
	md.vy=NaN;
	md.vz=NaN;
	md.vel=NaN;
	md.temperature=NaN; %temperature solution vector
	md.melting=NaN;
	md.pressure=NaN;

	%Parallelisation parameters
	md.cluster='none';
	md.np=0;
	md.time=0;
	md.queue='';
	md.exclusive=0;
	md.scheduler_configuration='local';
	md.outlog=' ';
	md.errlog=' ';
	md.alloc_cleanup=1;
	md.waitonlock=0;

	%mesh2grid
	md.mesh2grid_parameters={'vel'};
	md.mesh2grid_interpolation={'node'};
	md.mesh2grid_filter={'average'};
	md.mesh2grid_cornereast=0;
	md.mesh2grid_cornernorth=0;
	md.mesh2grid_xposting=1000;
	md.mesh2grid_yposting=1000;
	md.mesh2grid_nlines=1000;
	md.mesh2grid_ncolumns=1000;
	md.mesh2grid_windowsize=1;
	md.mesh2grid_results=NaN;
	md.dummy=NaN;

	%mesh 
	md.mesh_domainoutline='';
	md.mesh_riftoutline='';
	md.mesh_resolution=NaN;

	%Cielo parameters 
	md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40 ';
	
	%Needed for running in parallel
	md.analysis_type='';

	%management of large models
	md.repository='';

	%radar power images
	md.sarpwr=NaN;
	md.sarxm=NaN;
	md.sarym=NaN;

	%qmu
	md.variables=struct([]);
	md.responses=struct([]);
	md.dakotaresults=NaN;
	md.dakotain =NaN;
	md.dakotaout=NaN;
	md.dakotadat=NaN;
%	md.method               ='sampling';
	md.seed                 =1234;
	md.samples              =10;
	md.sample_type          ='lhs';
	md.method               ='local reliability';
	md.interval_type        ='forward';
	md.fd_gradient_step_size= 0.01;
	md.evaluation_concurrency= 2;
	md.npart                 = 1;

	md.analysis_driver    ='';
	md.analysis_components='';

	%Solver options
	%ice
	md.solver_type='general';
	%output
	md=class(md,'model');
				
case 1
	%If single argument of class model, we have a copy constructor. 
	if (isa(varargin{1},'model'))
		md = varargin{1};
	elseif (isa(varargin{1},'char'))
		%user gave us the name of the repository, this model is going to be large! 
		md=model;
		md.repository=varargin{1};
		md=activaterepository(md);
	else
		error('model constructor error message: copy constructor called on a non ''model'' class object');
	end 
otherwise
	error('model constructor error message: 0 of 1 argument only in input.');
end
