Index: issm/trunk/src/m/classes/model.m
===================================================================
--- issm/trunk/src/m/classes/model.m	(revision 9465)
+++ issm/trunk/src/m/classes/model.m	(revision 9505)
@@ -12,280 +12,280 @@
 
 		 %Model general information
-		 notes           = {'',false};
-		 name            = {'',true,'String'};
-		 runtimename     = {'',false}; %name used when running a parallel job
-		 counter         = {0,false};
+		 notes           = modelfield('default','','marshall',false);
+		 name            = modelfield('default','','marshall',true,'format','String');
+		 runtimename     = modelfield('default','','marshall',false); %name used when running a parallel job
+		 counter         = modelfield('default',0,'marshall',false);
 
 		 %Mesh
-		 bamg                = {struct(),false};
-		 dim                 = {0,true,'Integer'};
-		 numberofelements    = {0,true,'Integer'};
-		 numberofnodes       = {0,true,'Integer'};
-		 elements            = {NaN,true,'DoubleMat',2};
-		 elements_type       = {NaN,true,'DoubleMat',2};
-		 vertices_type       = {NaN,true,'DoubleMat',1};
-		 x                   = {NaN,true,'DoubleMat',1};
-		 y                   = {NaN,true,'DoubleMat',1};
-		 z                   = {NaN,true,'DoubleMat',1};
-		 nodeconnectivity    = {NaN,false};
-		 elementconnectivity = {NaN,true,'DoubleMat',3};
-		 edges               = {NaN,true,'DoubleMat',3};
+		 bamg                = modelfield('default',struct(),'marshall',false);
+		 dim                 = modelfield('default',0,'marshall',true,'format','Integer');
+		 numberofelements    = modelfield('default',0,'marshall',true,'format','Integer');
+		 numberofnodes       = modelfield('default',0,'marshall',true,'format','Integer');
+		 elements            = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 elements_type       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 vertices_type       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 x                   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 y                   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 z                   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeconnectivity    = modelfield('default',NaN,'marshall',false);
+		 elementconnectivity = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 edges               = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
 
 		 %I/O
-		 io_gather           = {NaN,true,'Boolean'};
+		 io_gather           = modelfield('default',0,'marshall',true,'format','Boolean');
 
 		 %Initial 2d mesh 
-		 numberofelements2d = {0,true,'Integer'};
-		 numberofnodes2d    = {0,true,'Integer'};
-		 elements2d         = {NaN,true,'DoubleMat',3};
-		 elements_type2d    = {NaN,false};
-		 vertices_type2d    = {NaN,false};
-		 x2d                = {NaN,false};
-		 y2d                = {NaN,false};
+		 numberofelements2d = modelfield('default',0,'marshall',true,'format','Integer');
+		 numberofnodes2d    = modelfield('default',0,'marshall',true,'format','Integer');
+		 elements2d         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 elements_type2d    = modelfield('default',NaN,'marshall',false);
+		 vertices_type2d    = modelfield('default',NaN,'marshall',false);
+		 x2d                = modelfield('default',NaN,'marshall',false);
+		 y2d                = modelfield('default',NaN,'marshall',false);
 
 		 %latlon of the coorinates
-		 lat        = {NaN,false};
-		 long       = {NaN,false};
-		 hemisphere = {NaN,false};
+		 lat        = modelfield('default',NaN,'marshall',false);
+		 long       = modelfield('default',NaN,'marshall',false);
+		 hemisphere = modelfield('default',NaN,'marshall',false);
 
 		 %Elements type
-		 ishutter             = {0,true,'Boolean'};
-		 ismacayealpattyn     = {0,true,'Boolean'};
-		 isstokes             = {0,true,'Boolean'};
+		 ishutter             = modelfield('default',0,'marshall',true,'format','Boolean');
+		 ismacayealpattyn     = modelfield('default',0,'marshall',true,'format','Boolean');
+		 isstokes             = modelfield('default',0,'marshall',true,'format','Boolean');
 
 		 %Elements
-		 elementonhutter      = {NaN,false};
-		 elementonmacayeal    = {NaN,false};
-		 elementonpattyn      = {NaN,false};
-		 elementonstokes      = {NaN,false};
+		 elementonhutter      = modelfield('default',NaN,'marshall',false);
+		 elementonmacayeal    = modelfield('default',NaN,'marshall',false);
+		 elementonpattyn      = modelfield('default',NaN,'marshall',false);
+		 elementonstokes      = modelfield('default',NaN,'marshall',false);
 
 		 %Nodes
-		 nodeonhutter         = {NaN,true,'DoubleMat',1};
-		 nodeonmacayeal       = {NaN,true,'DoubleMat',1};
-		 nodeonpattyn         = {NaN,true,'DoubleMat',1};
-		 nodeonstokes         = {NaN,true,'DoubleMat',1};
-		 borderstokes         = {NaN,true,'DoubleMat',3};
+		 nodeonhutter         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeonmacayeal       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeonpattyn         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeonstokes         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 borderstokes         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
 
 		 %Stokes
-		 stokesreconditioning = {0,true,'Double'};
-		 shelf_dampening      = {0,true,'Integer'};
+		 stokesreconditioning = modelfield('default',0,'marshall',true,'format','Double');
+		 shelf_dampening      = modelfield('default',0,'marshall',true,'format','Integer');
 
 		 %Penalties
-		 penalties            = {NaN,true,'DoubleMat',3};
-		 penalty_offset       = {0,true,'Double'};
-		 penalty_lock         = {0,true,'Integer'};
-		 segments             = {NaN,false};
-		 segmentmarkers       = {NaN,false};
-		 rifts                = {NaN,false};
-		 riftinfo             = {NaN,true,'DoubleMat',3};
-		 riftproperties       = {NaN,false};
-		 numrifts             = {0,true,'Integer'};
+		 penalties            = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 penalty_offset       = modelfield('default',0,'marshall',true,'format','Double');
+		 penalty_lock         = modelfield('default',0,'marshall',true,'format','Integer');
+		 segments             = modelfield('default',NaN,'marshall',false);
+		 segmentmarkers       = modelfield('default',NaN,'marshall',false);
+		 rifts                = modelfield('default',NaN,'marshall',false);
+		 riftinfo             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 riftproperties       = modelfield('default',NaN,'marshall',false);
+		 numrifts             = modelfield('default',0,'marshall',true,'format','Integer');
 
 		 %Projections
-		 uppernodes           = {NaN,false};
-		 upperelements        = {NaN,true,'DoubleMat',2};
-		 lowerelements        = {NaN,true,'DoubleMat',2};
-		 lowernodes           = {NaN,false};
+		 uppernodes           = modelfield('default',NaN,'marshall',false);
+		 upperelements        = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 lowerelements        = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 lowernodes           = modelfield('default',NaN,'marshall',false);
 
 		 %Extrusion
-		 numlayers         = {0,true,'Integer'};
-		 elementonbed      = {NaN,true,'BooleanMat',2};
-		 elementonsurface  = {NaN,true,'BooleanMat',2};
-		 nodeonbed         = {NaN,true,'BooleanMat',1};
-		 nodeonsurface     = {NaN,true,'BooleanMat',1};
-		 minh              = {0,false};
+		 numlayers         = modelfield('default',0,'marshall',true,'format','Integer');
+		 elementonbed      = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',2);
+		 elementonsurface  = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',2);
+		 nodeonbed         = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',1);
+		 nodeonsurface     = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',1);
+		 minh              = modelfield('default',0,'marshall',false);
 
 		 %Extraction
-		 extractednodes    = {NaN,false};
-		 extractedelements = {NaN,false};
+		 extractednodes    = modelfield('default',NaN,'marshall',false);
+		 extractedelements = modelfield('default',NaN,'marshall',false);
 
 		 %Materials parameters
-		 rho_ice                    = {0,true,'Double'};
-		 rho_water                  = {0,true,'Double'};
-		 heatcapacity               = {0,true,'Double'};
-		 latentheat                 = {0,true,'Double'};
-		 thermalconductivity        = {0,true,'Double'};
-		 meltingpoint               = {0,true,'Double'};
-		 referencetemperature       = {0,true,'Double'}; %for enthalpy
-		 beta                       = {0,true,'Double'};
-		 mixed_layer_capacity       = {0,true,'Double'};
-		 thermal_exchange_velocity  = {0,true,'Double'};
-		 min_thermal_constraints    = {0,true,'Integer'};
-		 min_mechanical_constraints = {0,true,'Integer'};
-		 stabilize_constraints      = {0,true,'Integer'};
+		 rho_ice                    = modelfield('default',0,'marshall',true,'format','Double');
+		 rho_water                  = modelfield('default',0,'marshall',true,'format','Double');
+		 heatcapacity               = modelfield('default',0,'marshall',true,'format','Double');
+		 latentheat                 = modelfield('default',0,'marshall',true,'format','Double');
+		 thermalconductivity        = modelfield('default',0,'marshall',true,'format','Double');
+		 meltingpoint               = modelfield('default',0,'marshall',true,'format','Double');
+		 referencetemperature       = modelfield('default',0,'marshall',true,'format','Double'); %for enthalpy
+		 beta                       = modelfield('default',0,'marshall',true,'format','Double');
+		 mixed_layer_capacity       = modelfield('default',0,'marshall',true,'format','Double');
+		 thermal_exchange_velocity  = modelfield('default',0,'marshall',true,'format','Double');
+		 min_thermal_constraints    = modelfield('default',0,'marshall',true,'format','Integer');
+		 min_mechanical_constraints = modelfield('default',0,'marshall',true,'format','Integer');
+		 stabilize_constraints      = modelfield('default',0,'marshall',true,'format','Integer');
 
 		 %Physical parameters
-		 g                = {0,true,'Double'};
-		 yts              = {0,true,'Double'};
-		 drag_coefficient = {NaN,true,'DoubleMat',1};
-		 drag_p           = {NaN,true,'DoubleMat',2};
-		 drag_q           = {NaN,true,'DoubleMat',2};
-		 rheology_B       = {NaN,true,'DoubleMat',1};
-		 rheology_n       = {NaN,true,'DoubleMat',2};
-		 rheology_law     = {'',true,'String'};
+		 g                = modelfield('default',0,'marshall',true,'format','Double');
+		 yts              = modelfield('default',0,'marshall',true,'format','Double');
+		 drag_coefficient = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 drag_p           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 drag_q           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 rheology_B       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 rheology_n       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 rheology_law     = modelfield('default','','marshall',true,'format','String');
 
 		 %Geometrical parameters
-		 elementoniceshelf = {NaN,true,'BooleanMat',2};
-		 elementonicesheet = {NaN,true,'BooleanMat',2};
-		 elementonwater    = {NaN,true,'BooleanMat',2};
-		 nodeoniceshelf    = {NaN,true,'DoubleMat',1};
-		 nodeonicesheet    = {NaN,true,'DoubleMat',1};
-		 nodeonwater       = {NaN,true,'DoubleMat',1};
-		 surface           = {NaN,true,'DoubleMat',1};
-		 thickness         = {NaN,true,'DoubleMat',1};
-		 thickness_coeff   = {NaN,true,'DoubleMat',1};
-		 bed               = {NaN,true,'DoubleMat',1};
-		 bathymetry        = {NaN,true,'DoubleMat',1};
+		 elementoniceshelf = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',2);
+		 elementonicesheet = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',2);
+		 elementonwater    = modelfield('default',NaN,'marshall',true,'format','BooleanMat','sizetype',2);
+		 nodeoniceshelf    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeonicesheet    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nodeonwater       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 surface           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 thickness         = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 thickness_coeff   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 bed               = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 bathymetry        = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
 
 		 %Boundary conditions
-		 nodeonboundary = {NaN,false};
-		 pressureload   = {NaN,true,'DoubleMat',3};
-		 spcvx          = {NaN,true,'DoubleMat',1};
-		 spcvy          = {NaN,true,'DoubleMat',1};
-		 spcvz          = {NaN,true,'DoubleMat',1};
-		 spctemperature = {NaN,true,'DoubleMat',1};
-		 spcthickness   = {NaN,true,'DoubleMat',1};
-		 spcwatercolumn = {NaN,true,'DoubleMat',1};
-		 diagnostic_ref = {NaN,true,'DoubleMat',1};
+		 nodeonboundary = modelfield('default',NaN,'marshall',false);
+		 pressureload   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 spcvx          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 spcvy          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 spcvz          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 spctemperature = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 spcthickness   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 spcwatercolumn = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 diagnostic_ref = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
 
 		 %Observations 
-		 vx_obs                    = {NaN,true,'DoubleMat',1};
-		 vy_obs                    = {NaN,true,'DoubleMat',1};
-		 vel_obs                   = {NaN,false};
-		 surface_accumulation_rate = {NaN,true,'DoubleMat',1};
-		 surface_ablation_rate     = {NaN,true,'DoubleMat',1};
-		 surface_mass_balance      = {NaN,true,'DoubleMat',1};
-		 dhdt                      = {NaN,true,'DoubleMat',1};
-		 geothermalflux            = {NaN,true,'DoubleMat',1};
-		 thickness_obs             = {NaN,true,'DoubleMat',1};
+		 vx_obs                    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 vy_obs                    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 vel_obs                   = modelfield('default',NaN,'marshall',false);
+		 surface_accumulation_rate = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 surface_ablation_rate     = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 surface_mass_balance      = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 dhdt                      = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 geothermalflux            = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 thickness_obs             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
 
 		 %Statics parameters
-		 eps_res                  = {0,true,'Double'};
-		 eps_rel                  = {0,true,'Double'};
-		 eps_abs                  = {0,true,'Double'};
-		 max_nonlinear_iterations = {0,true,'Integer'};
-		 max_steadystate_iterations = {0,true,'Integer'};
-		 connectivity             = {0,true,'Integer'};
-		 lowmem                   = {0,true,'Boolean'};
-		 viscosity_overshoot      = {0,true,'Double'};
+		 eps_res                  = modelfield('default',0,'marshall',true,'format','Double');
+		 eps_rel                  = modelfield('default',0,'marshall',true,'format','Double');
+		 eps_abs                  = modelfield('default',0,'marshall',true,'format','Double');
+		 max_nonlinear_iterations = modelfield('default',0,'marshall',true,'format','Integer');
+		 max_steadystate_iterations = modelfield('default',0,'marshall',true,'format','Integer');
+		 connectivity             = modelfield('default',0,'marshall',true,'format','Integer');
+		 lowmem                   = modelfield('default',0,'marshall',true,'format','Boolean');
+		 viscosity_overshoot      = modelfield('default',0,'marshall',true,'format','Double');
 
 		 %Transient 
-		 dt                     = {0,true,'Double'};
-		 ndt                    = {0,true,'Double'};
-		 time_adapt             = {0,true,'Boolean'};
-		 cfl_coefficient        = {0,true,'Double'};
-		 artificial_diffusivity = {0,true,'Integer'};
-		 prognostic_DG          = {0,true,'Integer'};
-		 hydrostatic_adjustment = {0,true,'Integer'};
-		 isprognostic = {0,true,'Boolean'};
-		 isdiagnostic = {0,true,'Boolean'};
-		 isthermal    = {0,true,'Boolean'};
+		 dt                     = modelfield('default',0,'marshall',true,'format','Double');
+		 ndt                    = modelfield('default',0,'marshall',true,'format','Double');
+		 time_adapt             = modelfield('default',0,'marshall',true,'format','Boolean');
+		 cfl_coefficient        = modelfield('default',0,'marshall',true,'format','Double');
+		 artificial_diffusivity = modelfield('default',0,'marshall',true,'format','Integer');
+		 prognostic_DG          = modelfield('default',0,'marshall',true,'format','Integer');
+		 hydrostatic_adjustment = modelfield('default',0,'marshall',true,'format','Integer');
+		 isprognostic = modelfield('default',0,'marshall',true,'format','Boolean');
+		 isdiagnostic = modelfield('default',0,'marshall',true,'format','Boolean');
+		 isthermal    = modelfield('default',0,'marshall',true,'format','Boolean');
 		 %Control
-		 control_analysis = {0,true,'Boolean'};
-		 control_type     = {0,true,'DoubleMat',3};
-		 weights          = {[],true,'DoubleMat',1};
-		 nsteps           = {0,true,'Integer'};
-		 maxiter          = {[],true,'DoubleMat',3};
-		 cm_responses     = {[],true,'DoubleMat',3};
-		 tolx             = {0,true,'Double'};
-		 optscal          = {[],true,'DoubleMat',3};
-		 eps_cm           = {0,true,'Double'};
-		 cm_min           = {NaN,true,'DoubleMat',3};
-		 cm_max           = {NaN,true,'DoubleMat',3};
-		 cm_jump          = {[],true,'DoubleMat',3};
-		 cm_gradient      = {0,true,'Boolean'};
-		 epsvel                              = {0,true,'Double'};
-		 meanvel                             = {0,true,'Double'};
-		 num_control_type                    = {0,true,'Integer'};
-		 num_cm_responses                    = {0,true,'Integer'};
+		 control_analysis = modelfield('default',0,'marshall',true,'format','Boolean');
+		 control_type     = modelfield('default',0,'marshall',true,'format','DoubleMat','sizetype',3);
+		 weights          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 nsteps           = modelfield('default',0,'marshall',true,'format','Integer');
+		 maxiter          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 cm_responses     = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 tolx             = modelfield('default',0,'marshall',true,'format','Double');
+		 optscal          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 eps_cm           = modelfield('default',0,'marshall',true,'format','Double');
+		 cm_min           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 cm_max           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 cm_jump          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
+		 cm_gradient      = modelfield('default',0,'marshall',true,'format','Boolean');
+		 epsvel                         = modelfield('default',0,'marshall',true,'format','Double');
+		 meanvel                        = modelfield('default',0,'marshall',true,'format','Double');
+		 num_control_type               = modelfield('default',0,'marshall',true,'format','Integer');
+		 num_cm_responses               = modelfield('default',0,'marshall',true,'format','Integer');
 		 %Output
-		 requested_outputs                    = {[],true,'DoubleMat',3};
+		 requested_outputs               = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',3);
 		 %Debugging
-		 verbose                             = {0,true,'Integer'};
-		 mem_debug                           = {0,false};
-		 gprof                               = {0,false};
+		 verbose                        = modelfield('default',0,'marshall',true,'format','Integer');
+		 mem_debug                      = modelfield('default',0,'marshall',false);
+		 gprof                          = modelfield('default',0,'marshall',false);
 		 %Results fields
-		 output_frequency                    = {0,true,'Integer'};
-		 results_on_vertices                 = {0,true,'Boolean'};
-		 inputfilename                       = {'',true,'String'};
-		 outputfilename                      = {'',true,'String'};
-		 results                             = {struct(),false};
-		 vx                                  = {NaN,true,'DoubleMat',1};
-		 vy                                  = {NaN,true,'DoubleMat',1};
-		 vz                                  = {NaN,true,'DoubleMat',1};
-		 vel                                 = {NaN,false};
-		 temperature                         = {NaN,true,'DoubleMat',1}; %temperature solution vector
-		 waterfraction                       = {NaN,true,'DoubleMat',1};
-		 groundingline_melting_rate                     = {NaN,true,'Double'};
-		 basal_melting_rate                  = {NaN,true,'DoubleMat',1};
-		 basal_melting_rate_correction       = {NaN,true,'DoubleMat',1};
-		 pressure                            = {NaN,true,'DoubleMat',1};
+		 output_frequency               = modelfield('default',0,'marshall',true,'format','Integer');
+		 results_on_vertices            = modelfield('default',0,'marshall',true,'format','Boolean');
+		 inputfilename                  = modelfield('default','','marshall',true,'format','String');
+		 outputfilename                 = modelfield('default','','marshall',true,'format','String');
+		 results                        = modelfield('default',struct(),'marshall',false);
+		 vx                             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 vy                             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 vz                             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 vel                            = modelfield('default',NaN,'marshall',false);
+		 temperature                    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1); %temperature solution vector
+		 waterfraction                  = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 groundingline_melting_rate                = modelfield('default',NaN,'marshall',true,'format','Double');
+		 basal_melting_rate             = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 basal_melting_rate_correction  = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
+		 pressure                       = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
 		 %Hydrology
-		 watercolumn                         = {NaN,true,'DoubleMat',1};
+		 watercolumn                    = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',1);
 
 		 %Hydrology
-		 hydro_n                             = {0,true,'Double'}; 
-		 hydro_CR                            = {0,true,'Double'};
-		 hydro_p                             = {0,true,'Double'};
-		 hydro_q                             = {0,true,'Double'};
-		 hydro_kn                            = {0,true,'Double'};
+		 hydro_n                        = modelfield('default',0,'marshall',true,'format','Double'); 
+		 hydro_CR                       = modelfield('default',0,'marshall',true,'format','Double');
+		 hydro_p                        = modelfield('default',0,'marshall',true,'format','Double');
+		 hydro_q                        = modelfield('default',0,'marshall',true,'format','Double');
+		 hydro_kn                       = modelfield('default',0,'marshall',true,'format','Double');
 		 
 		 %Parallelisation
-		 cluster       = {NaN,false};
-		 outlog        = {'',false};
-		 errlog        = {'',false};
-		 waitonlock    = {0,true,'Boolean'};
+		 cluster       = modelfield('default',NaN,'marshall',false);
+		 outlog        = modelfield('default','','marshall',false);
+		 errlog        = modelfield('default','','marshall',false);
+		 waitonlock    = modelfield('default',0,'marshall',true,'format','Boolean');
 
 		 %dummy
-		 dummy = {NaN,false};
+		 dummy = modelfield('default',NaN,'marshall',false);
 
 		 %PETSc and MATLAB solver string
-		 petscoptions  = {NaN,false};
+		 petscoptions  = modelfield('default',NaN,'marshall',false);
 
 		 %Analysis
-		 solution_type = {'',false};
+		 solution_type = modelfield('default','','marshall',false);
 
 		 %radar power images
-		 sarpwr = {NaN,false};
-		 sarxm  = {NaN,false};
-		 sarym  = {NaN,false};
+		 sarpwr = modelfield('default',NaN,'marshall',false);
+		 sarxm  = modelfield('default',NaN,'marshall',false);
+		 sarym  = modelfield('default',NaN,'marshall',false);
 
 		 %qmu
-		 variables                       = {struct(),false};
-		 responses                       = {struct(),false};
-		 qmu_method                      = {struct(),false};
-		 qmu_params                      = {struct(),false};
-		 dakotaresults                   = {struct(),false};
-		 dakotain                        = {'',false};
-		 dakotaout                       = {'',false};
-		 dakotadat                       = {'',false};
-		 qmu_analysis                    = {0,true,'Boolean'};
-		 part                            = {[],true,'DoubleMat',2};
-		 npart                           = {0,true,'Integer'};
-		 numberofvariables               = {0,true,'Integer'};
-		 numberofresponses               = {0,true,'Integer'};
-		 variabledescriptors             = {{},true,'StringArray'};
-		 responsedescriptors             = {{},true,'StringArray'};
-		 qmu_mass_flux_profile_directory = {NaN,false};
-		 qmu_mass_flux_profiles          = {NaN,false};
-		 qmu_mass_flux_segments          = {{},true,'MatArray'};
-		 qmu_relax                       = {0,false};
+		 variables                       = modelfield('default',struct(),'marshall',false);
+		 responses                       = modelfield('default',struct(),'marshall',false);
+		 qmu_method                      = modelfield('default',struct(),'marshall',false);
+		 qmu_params                      = modelfield('default',struct(),'marshall',false);
+		 dakotaresults                   = modelfield('default',struct(),'marshall',false);
+		 dakotain                        = modelfield('default','','marshall',false);
+		 dakotaout                       = modelfield('default','','marshall',false);
+		 dakotadat                       = modelfield('default','','marshall',false);
+		 qmu_analysis                    = modelfield('default',0,'marshall',true,'format','Boolean');
+		 part                            = modelfield('default',NaN,'marshall',true,'format','DoubleMat','sizetype',2);
+		 npart                           = modelfield('default',0,'marshall',true,'format','Integer');
+		 numberofvariables               = modelfield('default',0,'marshall',true,'format','Integer');
+		 numberofresponses               = modelfield('default',0,'marshall',true,'format','Integer');
+		 variabledescriptors             = modelfield('default',{},'marshall',true,'format','StringArray');
+		 responsedescriptors             = modelfield('default',{},'marshall',true,'format','StringArray');
+		 qmu_mass_flux_profile_directory = modelfield('default',NaN,'marshall',false);
+		 qmu_mass_flux_profiles          = modelfield('default',NaN,'marshall',false);
+		 qmu_mass_flux_segments          = modelfield('default',{},'marshall',true,'format','MatArray');
+		 qmu_relax                       = modelfield('default',0,'marshall',false);
 
 		 %flaim
-		 fm_tracks             = {'',false};
-		 fm_flightreqs         = {struct(),false};
-		 fm_criterion          = {[],false};
-		 fm_gridsatequator     = {200000,false};
-		 fm_usevalueordering   = {true,false};
-		 fm_split_antimeridian = {true,false};
-		 fm_solution           = {'',false};
-		 fm_quality            = {0,false};
+		 fm_tracks             = modelfield('default','','marshall',false);
+		 fm_flightreqs         = modelfield('default',struct(),'marshall',false);
+		 fm_criterion          = modelfield('default',NaN,'marshall',false);
+		 fm_gridsatequator     = modelfield('default',200000,'marshall',false);
+		 fm_usevalueordering   = modelfield('default',true,'marshall',false);
+		 fm_split_antimeridian = modelfield('default',true,'marshall',false);
+		 fm_solution           = modelfield('default','','marshall',false);
+		 fm_quality            = modelfield('default',0,'marshall',false);
 
 		 %grounding line migration: 
-		 groundingline_migration = {0,true,'Integer'};
+		 groundingline_migration = modelfield('default',0,'marshall',true,'format','Integer');
 
 		 %partitioner:
-		 adjacency = {[],false};
-		 vwgt      = {[],false};
+		 adjacency = modelfield('default',NaN,'marshall',false);
+		 vwgt      = modelfield('default',NaN,'marshall',false);
 		 %}}}
 	 end
@@ -511,5 +511,5 @@
 			 for i=1:length(fieldnames),
 				 fieldname=fieldnames{i};
-				 md.(fieldname)=md.(fieldname){1};
+				 md.(fieldname)=md.(fieldname).default;
 			 end
 
