Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 359)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 360)
@@ -21,5 +21,5 @@
 	
 	%Mesh
-	md.type='2d';
+	md.type='';
 	md.numberofelements=0;
 	md.numberofgrids=0;
@@ -61,13 +61,13 @@
 
 	%Stokes mesh
-	md.stokesreconditioning=10^13;
+	md.stokesreconditioning=0;
 	
 	%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.penalty_offset=0;
+	md.penalty_melting=0;
+	md.penalty_lock=0;
+	md.segments=NaN;
+	md.segmentmarkers=NaN;
 	md.rifts=NaN;
 	md.numrifts=0;
@@ -84,5 +84,5 @@
 	md.gridonbed=NaN;
 	md.gridonsurface=NaN;
-	md.minh=1;
+	md.minh=0;
 	md.firn_layer=NaN;
 	md.dof=NaN;
@@ -93,19 +93,19 @@
 	
 	%Materials parameters
-	md.rho_ice=917;
-	md.rho_water=1023;
-	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
+	md.rho_ice=0;
+	md.rho_water=0;
+	md.heatcapacity=0;
+	md.latentheat=0;
+	md.thermalconductivity=0;
+	md.meltingpoint=0;
+	md.beta=0;
+	md.mixed_layer_capacity=0;
+	md.thermal_exchange_velocity=0;
 	
 	%Physical parameters
-	md.g=9.81;
-	md.yts=365*24*3600;
+	md.g=0;
+	md.yts=0;
 	md.drag=NaN;
-	md.drag_type=2;
+	md.drag_type=0;
 	md.p=NaN;
 	md.q=NaN;
@@ -132,5 +132,5 @@
 
 	%Thermal
-	md.min_thermal_constraints=0; %minimum number of unstable thermal constraints
+	md.min_thermal_constraints=0;
 	md.gridondirichlet_thermal=NaN;
 	md.dirichletvalues_thermal=NaN;
@@ -153,10 +153,10 @@
 
 	%Statics parameters
-	md.eps_rel=0.01;
-	md.eps_abs=10;
-	md.sparsity=0.001;
-	md.connectivity=16;
+	md.eps_rel=0;
+	md.eps_abs=0;
+	md.sparsity=0;
+	md.connectivity=0;
 	md.lowmem=0;
-	md.viscosity_overshoot=0.3;
+	md.viscosity_overshoot=0;
 	
 	%Transient parameters
@@ -166,15 +166,14 @@
 	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;
-
+	md.timestepping=0;
+	md.deltaH=0;
+	md.DeltaH=0;
+	md.deltaT=0;
+	md.DeltaT=0;
+	md.timeacc=0;
+	md.timedec=0;
 
 	%Control
-	md.control_type='drag';
+	md.control_type='';
 	md.cont_vx=NaN;
 	md.cont_vy=NaN;
@@ -185,11 +184,11 @@
 	md.nsteps=0;
 	md.maxiter=[];
-	md.tolx=10^-4;
+	md.tolx=0;
 	md.optscal=[];
 	md.mincontrolconstraint=0;
-	md.maxcontrolconstraint=2000;
+	md.maxcontrolconstraint=0;
 	md.fit=[];
-	md.epsvel=eps;
-	md.meanvel=1000/(365*24*3600); %1000 m/an
+	md.epsvel=0;
+	md.meanvel=0;
 
 	%Output parameters
@@ -202,8 +201,8 @@
 	md.deviatoricstress=NaN;
 	md.strainrate=NaN;
-	md.plot=1;
+	md.plot=0;
 
 	%debugging
-	md.debug=1;
+	md.debug=0;
 	md.element_debug=0;
 	md.element_debugid=NaN;
@@ -219,26 +218,26 @@
 
 	%Parallelisation parameters
-	md.cluster='none';
+	md.cluster='';
 	md.np=0;
 	md.time=0;
 	md.queue='';
 	md.exclusive=0;
-	md.scheduler_configuration='local';
-	md.outlog=' ';
-	md.errlog=' ';
+	md.scheduler_configuration='';
+	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_parameters={};
+	md.mesh2grid_interpolation={};
+	md.mesh2grid_filter={};
 	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_xposting=0;
+	md.mesh2grid_yposting=0;
+	md.mesh2grid_nlines=0;
+	md.mesh2grid_ncolumns=0;
+	md.mesh2grid_windowsize=0;
 	md.mesh2grid_results=NaN;
 	md.dummy=NaN;
@@ -249,8 +248,6 @@
 	md.mesh_resolution=NaN;
 
-	%Cielo parameters 
-	md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40 ';
-
-	%Trash, still to organize
+	%Cielo solver string
+	md.solverstring='';
 	md.analysis_type='';
 
@@ -264,21 +261,24 @@
 
 	%qmu
-	md.variables =struct();
-	md.responses =struct();
+	md.variables=struct();
+	md.responses=struct();
 	md.qmu_method=struct();
 	md.qmu_params=struct();
 	md.dakotaresults=struct();
-	md.dakotain     ='';
-	md.dakotaout    ='';
-	md.dakotadat    ='';
-
-	md.npart   =1;
-
-	%Solver options
-	%ice
-	md.solver_type='general';
+	md.dakotain='';
+	md.dakotaout='';
+	md.dakotadat='';
+
+	md.npart=0;
+
+	%Ice solver string
+	md.solver_type='';
+
 	%output
 	md=class(md,'model');
-				
+
+	%set default parameters
+	md=setdefaultparameters(md);
+
 case 1
 	%If single argument of class model, we have a copy constructor. 
Index: /issm/trunk/src/m/classes/@model/setdefaultparameters.m
===================================================================
--- /issm/trunk/src/m/classes/@model/setdefaultparameters.m	(revision 360)
+++ /issm/trunk/src/m/classes/@model/setdefaultparameters.m	(revision 360)
@@ -0,0 +1,262 @@
+function md=setdefaultparameters(md);
+%SETDEFAULTPARAMETERS - plug default parameters onto model
+%
+%   Although the model parameterization should be done in
+%   the parameter file, some parameters are initialized her
+%   with a default value.
+%   These default values can be changed if necessary.
+%
+%   Usage:
+%      md=setdefaultparameters(md);
+
+
+%Materials parameters
+
+%ice density (kg/m^3)
+md.rho_ice=917;
+
+%water density (kg/m^3)
+md.rho_water=1023;
+
+%ice heat capacity cp (J/kg/K)
+md.heatcapacity=2093;
+
+%ice latent heat of fusion L (J/kg)
+md.latentheat=3.34*10^5;
+
+%ice thermal conductivity lamda (W/m/K)
+md.thermalconductivity=2.4;
+
+%the melting point of ice at 1 atmosphere of pressure in K
+md.meltingpoint=273.15;
+
+%rate of change of melting point with pressure (K/Pa)
+md.beta=9.8*10^-8;
+
+%mixed layer (ice-water interface) heat capacity (J/kg/K)
+md.mixed_layer_capacity=3974;
+
+%thermal exchange velocity (ice-water interface) (m/s)
+md.thermal_exchange_velocity=1.00*10^-4;
+
+%Physical parameters
+
+%acceleration due to gravity (m/s^2)
+md.g=9.81;
+
+%converstion from year to seconds
+md.yts=365*24*3600;
+
+%drag type (2=elastic)
+md.drag_type=2;
+
+%Solver parameters
+
+%relative convergence criterion ((vel(n-1)-vel(n))/vel(n))
+md.eps_rel=0.01;
+
+%absolute convergence criterion (max(vel(n-1)-vel(n)) (m/year)
+md.eps_abs=10;
+
+%sparsity
+md.sparsity=0.001;
+
+%the connectivity is the avergaded number of nodes linked to a
+%given node through an edge. This connectivity is used to initially
+%allocate memory to the stiffness matrix. A value of 16 seems to
+%give a good memory/time ration. This value can be checked in
+%trunk/test/Miscellaneous/runme.m
+md.connectivity=16;
+
+%lowmem??
+md.lowmem=0;
+
+%coefficient to update the viscosity between each iteration of
+%a diagnostic according to the following formula
+%viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
+md.viscosity_overshoot=0.3;
+
+%parameter used to print temporary results (convergence criterion,
+%current step,...)
+md.debug=1;
+
+%Stokes mesh
+md.stokesreconditioning=10^13;
+
+%Penalties
+
+%the penalty offset is the power used to compute the value
+%of the penalty as follows
+%kappa=max(stiffness matrix)*10^penalty_offset
+md.penalty_offset=3;
+
+%penalty_melting is the penalty used to keep T<Tpmp
+%kappa=penalty_melting
+md.penalty_melting=10^7;
+
+%in some solutions, it might be needed to stop a run when only
+%a few constraints remain unstable. For thermal computation, this
+%parameter is often used.
+md.penalty_lock=10;
+
+%minimum thickness to avoid stiffness singularity, used in
+%the parameter file
+md.minh=1;
+
+%in some cases, it might be needed to stop a run when only
+%a few constraints remain unstable. For thermal computation, this
+%parameter is often used.
+md.min_thermal_constraints=0;
+
+%Transient parameters
+
+%time between 2 time steps
+md.dt=md.yts/2;
+
+%number of time steps = md.ndt / md.dt
+md.ndt=10*md.dt;
+
+%the artificial diffusivity is used in prognostic to avoid
+%numerical wiggles of the solution.
+md.artificial_diffusivity=1;
+
+%do we use an adaptative time stepping?
+md.timestepping=0;
+
+%minimum thickness change between 2 time steps (m)
+md.deltaH=10^-2;
+
+%maximum thickness change between 2 time steps (m)
+md.DeltaH=50;
+
+%minimum temperature change between 2 time steps (K)
+md.deltaT=10^-2;
+
+%maximum temperature change between 2 time steps (K)
+md.DeltaT=1;
+
+%When the previous criteria are fulfilled, if the time steps
+%can be accelerated as follows
+%dt = timeacc * dt
+md.timeacc=2;
+
+%When the previous criteria are fulfilled, if the time steps
+%can be reduced as follows
+%dt = timedec * dt
+md.timedec=.5;
+
+%Control
+
+%parameter to be inferred by control methods (only
+%drag and B are supported yet)
+md.control_type='drag';
+
+%number of steps in the control methods
+md.nsteps=20;
+
+%maximum number of iteration in the optimization algorithm for
+%each step
+md.maxiter=20;
+
+%tolerance used by the optimization algorithm
+md.tolx=10^-4;
+
+%the inversed parameter is updated as follows:
+%new_par=old_par + optscal(n)*C*gradient with C in [0 1];
+%usually the optscal must be of the order of magnitude of the 
+%inversed parameter (10^8 for B, 50 for drag) and can be decreased
+%after the first iterations
+md.optscal=50*ones(md.nsteps,1);
+
+%mincontrolconstraint is the minimum acceptable value of the
+%inversed parameter (B>0 for example)
+md.mincontrolconstraint=0;
+
+%maxcontrolconstraint is the maximum acceptable value of the
+%inversed parameter (drag<200 for example)
+md.maxcontrolconstraint=200;
+
+%several misfits can be used:
+%0 = absolute misfit
+%1 = relative misfit
+%2 = logarithmic misfit
+md.fit=0*ones(md.nsteps,1);
+
+%minimum velocity to avoid the misfit to be singular
+md.epsvel=eps;
+
+%averaged velocity used to scale the logarithmic Misfit (000 m/an)
+md.meanvel=1000/(365*24*3600);
+
+%the value of plot can be set to 1 if the user wants to see
+%the direction and the updated parameter after each control
+%method iteration (only for cielo serial)
+md.plot=1;
+
+%Parallelisation parameters
+
+%cluster set as 'none' for serial
+md.cluster='none';
+
+%number of processors to be used (if parallel)
+md.np=0;
+
+%?????? (s)
+md.time=0;
+
+%???
+md.queue='';
+
+%?????
+md.exclusive=0;
+
+%?????
+md.scheduler_configuration='local';
+
+%?????
+md.alloc_cleanup=1;
+
+%this option can be activated (1) to load automatically the results
+%onto the model after a parallel run by waiting for the lock file
+%tht is generated once the solution has converged
+md.waitonlock=0;
+
+%mesh2grid parameter
+md.mesh2grid_parameters={'vel'};
+
+%interpolation type
+md.mesh2grid_interpolation={'node'};
+
+%filter
+md.mesh2grid_filter={'average'};
+
+%coordinate of the corner east (actually WEST...)
+md.mesh2grid_cornereast=0;
+
+%coordinate of the corner north
+md.mesh2grid_cornernorth=0;
+
+%x posting
+md.mesh2grid_xposting=1000;
+
+%y posting
+md.mesh2grid_yposting=1000;
+
+%number of lines
+md.mesh2grid_nlines=1000;
+
+%number of columns
+md.mesh2grid_ncolumns=1000;
+
+%window size
+md.mesh2grid_windowsize=1;
+
+%the string of solverstring is used directly by PETSc to solve finite element
+%systems KU=F. By default, we use MUMPS solver
+md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40 ';
+
+%Ice solver: 'general' for Matlab's default solver (or 'lu' or 'sholesky')
+md.solver_type='general';
+
+%?????
+md.npart   =1;