Index: sm/trunk/src/m/classes/model.m
===================================================================
--- /issm/trunk/src/m/classes/model.m	(revision 8457)
+++ 	(revision )
@@ -1,943 +1,0 @@
-%MODEL class definition
-%
-%   Usage:
-%      md = model(varargin)
-
-classdef model
-    properties (SetAccess=public) %Model fields
-		 % {{{1
-		 %Careful here: no other class should be used as default value
-		 %this is a bug of matlab
-
-		 %Model general information
-		 notes='';
-		 name=''; 
-		 runtimename=''; %name used when running a parallel job
-		 counter=0; %1 mesh, 2 geography, 3 parameterize
-		 domainoutline=NaN;
-		 riftoutline=NaN;
-		 iceshelfoutline=NaN;
-		 icesheetoutline=NaN;
-		 wateroutline=NaN;
-		 parameterfile=NaN;
-		 runmefile=NaN;
-
-		 %Mesh
-		 bamg=struct();
-		 dim=0;
-		 numberofelements=0;
-		 numberofnodes=0;
-		 elements=NaN;
-		 elements_type=NaN;
-		 vertices_type=NaN;
-		 x=NaN;
-		 y=NaN;
-		 z=NaN;
-		 bed_slopex=NaN;
-		 bed_slopey=NaN;
-		 surface_slopex=NaN;
-		 surface_slopey=NaN;
-		 nodeconnectivity=NaN;
-		 elementconnectivity=NaN;
-		 edges=NaN;
-
-		 %I/O
-		 io_gather=NaN;
-
-		 %Initial 2d mesh 
-		 numberofelements2d=0;
-		 numberofnodes2d=0;
-		 elements2d=NaN;
-		 elements_type2d=NaN;
-		 vertices_type2d=NaN;
-		 x2d=NaN;
-		 y2d=NaN;
-		 z2d=NaN;
-
-		 %latlon of the coorinates
-		 lat=NaN;
-		 long=NaN;
-		 hemisphere=NaN;
-
-		 %Elements type
-		 ishutter=0;
-		 ismacayealpattyn=0;
-		 isstokes=0;
-
-		 %Elements
-		 elementonhutter=NaN;
-		 elementonmacayeal=NaN;
-		 elementonpattyn=NaN;
-		 elementonstokes=NaN;
-
-		 %Nodes
-		 nodeonhutter=NaN;
-		 nodeonmacayeal=NaN;
-		 nodeonpattyn=NaN;
-		 nodeonstokes=NaN;
-		 borderstokes=NaN;
-
-		 %Stokes
-		 stokesreconditioning=0;
-		 shelf_dampening=0;
-
-		 %Penalties
-		 penalties=NaN;
-		 penalty_offset=0;
-		 penalty_melting=0;
-		 penalty_lock=0;
-		 segments=NaN;
-		 segmentmarkers=NaN;
-		 rifts=NaN;
-		 riftinfo=NaN;
-		 riftproperties=NaN;
-		 numrifts=0;
-
-		 %Projections
-		 uppernodes=NaN;
-		 upperelements=NaN;
-		 lowerelements=NaN;
-		 lowernodes=NaN;
-
-		 %Extrusion
-		 numlayers=0;
-		 extrusionexponent=0;
-		 elementonbed=NaN;
-		 elementonsurface=NaN;
-		 nodeonbed=NaN;
-		 nodeonsurface=NaN;
-		 minh=0;
-		 firn_layer=NaN;
-
-		 %Extraction
-		 extractednodes=NaN;
-		 extractedelements=NaN;
-
-		 %Materials parameters
-		 rho_ice=0;
-		 rho_water=0;
-		 rho_firn=0;
-		 heatcapacity=0;
-		 latentheat=0;
-		 thermalconductivity=0;
-		 meltingpoint=0;
-		 beta=0;
-		 
-		 mixed_layer_capacity=0;
-		 thermal_exchange_velocity=0;
-		 min_thermal_constraints=0;
-		 min_mechanical_constraints=0;
-		 stabilize_constraints=0;
-
-		 %Physical parameters
-		 g=0;
-		 yts=0;
-		 drag_type=0;
-		 drag_coefficient=NaN;
-		 drag_p=NaN;
-		 drag_q=NaN;
-		 rheology_B=NaN;
-		 rheology_n=NaN;
-		 rheology_law=0;
-
-		 %Geometrical parameters
-		 elementoniceshelf=NaN;
-		 elementonicesheet=NaN;
-		 elementonwater=NaN;
-		 elementonnuna=NaN;
-		 nodeoniceshelf=NaN;
-		 nodeonicesheet=NaN;
-		 nodeonwater=NaN;
-		 nodeonnuna=NaN;
-		 surface=NaN;
-		 thickness=NaN;
-		 thickness_coeff=NaN;
-		 bed=NaN;
-		 bathymetry=NaN;
-		 mask=NaN;
-
-		 %Boundary conditions
-		 nodeonboundary=NaN;
-		 pressureload=NaN;
-		 spcvelocity=NaN;
-		 spctemperature=NaN;
-		 spcthickness=NaN;
-		 spcwatercolumn=NaN;
-		 diagnostic_ref=NaN;
-
-		 %Observations 
-		 vx_obs=NaN;
-		 vy_obs=NaN;
-		 vel_obs=NaN;
-		 vx_obs_raw=NaN;
-		 vy_obs_raw=NaN;
-		 vx_bal=NaN;
-		 vy_bal=NaN;
-		 vel_bal=NaN;
-		 vel_obs_raw=NaN;
-		 surface_accumulation_rate=NaN;
-		 surface_ablation_rate=NaN;
-		 surface_mass_balance=NaN;
-		 dhdt=NaN;
-		 geothermalflux=NaN;
-		 observed_temperature=NaN;
-		 thickness_obs=NaN;
-
-		 %Forcings
-		 forcings=struct();
-
-		 %Statics parameters
-		 eps_res=0;
-		 eps_rel=0;
-		 eps_abs=0;
-		 max_nonlinear_iterations=0;
-		 sparsity=0;
-		 connectivity=0;
-		 lowmem=0;
-		 viscosity_overshoot=0;
-
-		 %Transient parameters
-		 dt=0;
-		 ndt=0;
-		 time_adapt=0;
-		 cfl_coefficient=1/2;
-		 adaptative_cfl=0;
-		 artificial_diffusivity=0;
-		 prognostic_DG=0;
-		 deltaH=0;
-		 DeltaH=0;
-		 deltaT=0;
-		 DeltaT=0;
-		 timeacc=0;
-		 timedec=0;
-		 hydrostatic_adjustment=0;
-
-		 %Transient flags
-		 isprognostic=0;
-		 isdiagnostic=0;
-		 isthermal=0;
-
-		 %Control
-		 control_analysis=0;
-		 control_type=0;
-		 num_control_type=0;
-		 weights=[];
-		 nsteps=0;
-		 maxiter=[];
-		 cm_responses=[];
-		 tolx=0;
-		 optscal=[];
-		 eps_cm=0;
-		 cm_noisedmp=0;
-		 cm_min=NaN;
-		 cm_max=NaN;
-		 cm_jump=[];
-		 cm_gradient=0;
-		 epsvel=0;
-		 meanvel=0;
-
-		 %Output parameters
-		 numoutput=0;
-		 parameteroutput={};
-		 viscousheating=NaN;
-		 pressure_elem=NaN;
-		 stress=NaN;
-		 stress_surface=NaN;
-		 stress_bed=NaN;
-		 deviatoricstress=NaN;
-		 strainrate=NaN;
-
-		 %debugging
-		 verbose=0;
-		 verbose_binary=0;
-		 element_debug=0;
-		 element_debugid=NaN;
-		 mem_debug=0;
-		 gprof=0;
-		 memoryleaks=NaN;
-
-		 %Results fields
-		 output_frequency=0;
-		 inputfilename='';
-		 outputfilename='';
-		 results=struct();
-		 vx=NaN;
-		 vy=NaN;
-		 vz=NaN;
-		 vel=NaN;
-		 temperature=NaN; %temperature solution vector
-		 gl_melting_rate=NaN;
-		 basal_melting_rate=NaN;
-		 basal_melting_rate_correction=NaN;
-		 basal_melting_rate_correction_apply=0;
-		 pressure=NaN;
-		 
-		 %hydrology
-		 watercolumn=NaN;
-		 hydro_n=NaN;
-		 hydro_CR=NaN;
-		 hydro_p=NaN;
-		 hydro_q=NaN;
-		 hydro_gamma=0;
-		 hydro_kn=0;
-		 
-		 %Parallelisation parameters
-		 cluster=NaN;
-		 outlog='';
-		 errlog='';
-		 alloc_cleanup=1;
-		 waitonlock=0;
-
-		 %mesh2grid
-		 mesh2grid_parameters={};
-		 mesh2grid_data=NaN;
-		 mesh2grid_xm=NaN;
-		 mesh2grid_ym=NaN;
-
-		 %dummy
-		 dummy=NaN;
-
-		 %mesh 
-		 mesh_domainoutline='';
-		 mesh_riftoutline='';
-		 mesh_resolution=NaN;
-
-		 %PETSc and MATLAB solver string
-		 petscoptions=NaN;
-
-		 %Analysis
-		 solution_type='';
-
-		 %management of large models
-		 repository='';
-
-		 %radar power images
-		 sarpwr=NaN;
-		 sarxm=NaN;
-		 sarym=NaN;
-
-		 %qmu
-		 variables=struct();
-		 responses=struct();
-		 qmu_method=struct();
-		 qmu_params=struct();
-		 dakotaresults=struct();
-		 dakotain='';
-		 dakotaout='';
-		 dakotadat='';
-		 qmu_analysis=0;
-		 part=[];
-		 npart=0;
-		 numberofvariables=0;
-		 numberofresponses=0;
-		 variabledescriptors=NaN;
-		 responsedescriptors=NaN;
-		 numvariabledescriptors=0;
-		 numresponsedescriptors=0;
-		 qmu_mass_flux_num_profiles=NaN;
-		 qmu_mass_flux_profile_directory=NaN;
-		 qmu_mass_flux_profiles=NaN;
-		 qmu_mass_flux_segments=NaN;
-		 qmu_relax=0;
-		 qmu_save_femmodel=0;
-
-		 %flaim
-		 fm_tracks='';
-		 fm_flightreqs=struct();
-		 fm_criterion=[];
-		 fm_solution='';
-		 fm_quality=0;
-
-		 %grounding line migration: 
-		 gl_migration=NaN;
-
-		 %solution parameters
-		 kff=0;
-		 loadresults=0;
-
-		 %partitioner:
-		 adjacency=[];
-		 vwgt=[];
-		 %}}}
-	 end
-	 properties(Dependent,Hidden) %Deprecated property
-		 % {{{
-		 % We must keep here ALL deprecated properties
-		 % This ensures that when a model object from a previous version
-		 % all fields are recovered.
-		 % loadobj is then used to assign old properties to new properties
-		 % But if we don't keep the old properties here, they would not be
-		 % recovered when the model is loaded
-		 % KEEP IN CHRONOLOGICAL ORDER (OLDER FIRST)
-		 drag;
-		 p;
-		 q;
-		 B;
-		 n;
-		 melting;
-		 accumulation;
-		 type;
-		 numberofgrids;
-		 numberofgrids2d;
-		 gridonhutter;
-		 gridonmacayeal;
-		 gridonpattyn;
-		 gridonstokes;
-		 uppergrids;
-		 lowergrids;
-		 gridonbed;
-		 gridonsurface;
-		 extractedgrids;
-		 gridoniceshelf;
-		 gridonicesheet;
-		 gridonwater;
-		 gridonnuna;
-		 gridonboundary;
-		 melting_rate;
-		 accumulation_rate;
-		 % }}}
-	 end
-	 methods (Static)
-		 function md = loadobjk(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 isstruct(md)
-				 disp('Recovering model object from a previous version');
-
-				 %Ok the loaded model is a struct, initialize output and recover all fields
-				 structmd=md;
-				 md=model;
-
-				 structfields=fields(structmd);
-				 modelprops  =properties(md);
-
-				 for i=1:length(structfields),
-					 fieldname =structfields{i};
-					 fieldvalue=getfield(structmd,fieldname);
-					 if ismember(fieldname,modelprops),
-						 md=setfield(md,fieldname,fieldvalue);
-					 end
-				 end
-
-				 %Field name change
-				 if isfield(structmd,'drag'), md.drag_coefficient=structmd.drag; end
-				 if isfield(structmd,'p'), md.drag_p=structmd.p; end
-				 if isfield(structmd,'q'), md.drag_q=structmd.p; end
-				 if isfield(structmd,'B'), md.rheology_B=structmd.B; end
-				 if isfield(structmd,'n'), md.rheology_n=structmd.n; end
-				 if isfield(structmd,'melting'), md.basal_melting_rate=structmd.melting; end
-				 if isfield(structmd,'melting_rate'), md.basal_melting_rate=structmd.melting_rate; end
-				 if isfield(structmd,'accumulation'), md.surface_mass_balance=structmd.accumulation; end
-				 if isfield(structmd,'accumulation_rate'), md.surface_mass_balance=structmd.accumulation_rate; end
-				 if isfield(structmd,'numberofgrids'), md.numberofnodes=structmd.numberofgrids; end
-				 if isfield(structmd,'numberofgrids2d'), md.numberofnodes2d=structmd.numberofgrids2d; end
-				 if isfield(structmd,'gridonhutter'), md.nodeonhutter=structmd.gridonhutter; end
-				 if isfield(structmd,'gridonmacayeal'), md.nodeonmacayeal=structmd.gridonmacayeal; end
-				 if isfield(structmd,'gridonpattyn'), md.nodeonpattyn=structmd.gridonpattyn; end
-				 if isfield(structmd,'gridonstokes'), md.nodeonstokes=structmd.gridonstokes; end
-				 if isfield(structmd,'uppergrids'), md.uppernodes=structmd.uppergrids; end
-				 if isfield(structmd,'lowergrids'), md.lowernodes=structmd.lowergrids; end
-				 if isfield(structmd,'gridonbed'), md.nodeonbed=structmd.gridonbed; end
-				 if isfield(structmd,'gridonsurface'), md.nodeonsurface=structmd.gridonsurface; end
-				 if isfield(structmd,'extractedgrids'), md.extractednodes=structmd.extractedgrids; end
-				 if isfield(structmd,'gridoniceshelf'), md.nodeoniceshelf=structmd.gridoniceshelf; end
-				 if isfield(structmd,'gridonicesheet'), md.nodeonicesheet=structmd.gridonicesheet; end
-				 if isfield(structmd,'gridonwater'), md.nodeonwater=structmd.gridonwater; end
-				 if isfield(structmd,'gridonnuna'), md.nodeonnuna=structmd.gridonnuna; end
-				 if isfield(structmd,'gridonboundary'), md.nodeonboundary=structmd.gridonboundary; end
-
-				 %Field class change
-				 if (isfield(structmd,'type') & ischar(structmd.type)), 
-					 if strcmpi(structmd.type,'2d'), md.dim=2; end
-					 if strcmpi(structmd.type,'3d'), md.dim=3; end
-				 end
-				 if isnumeric(md.verbose), md.verbose=verbose; end
-
-				 %New fields
-				 if ~isfield(structmd,'upperelements');
-					 md.upperelements=transpose(1:md.numberofelements)+md.numberofelements2d;
-					 md.upperelements(end-md.numberofelements2d+1:end)=NaN;
-				 end
-				 if ~isfield(structmd,'lowerelements');
-					 md.lowerelements=transpose(1:md.numberofelements)-md.numberofelements2d;
-					 md.lowerelements(1:md.numberofelements2d)=NaN;
-				 end
-			 end
-		 end% }}}
-	 end
-	 methods
-		 function md=model(varargin) % {{{1
-
-			 switch nargin
-				 case 0
-
-					 %set default parameters
-					 md=setdefaultparameters(md);
-
-				 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
-		 end
-		 %}}}
-		 function md=activaterepository(md) % {{{1
-			 %ACTIVATEREPOSITORY - save the model fields separately
-			 %
-			 %   This function is required when a matlab session 
-			 %   is running low on memory, or when a model is very large. 
-			 %   By activating a repository (for example, 'repository.model'), 
-			 %   the model fields will each separately be saved and accessed when 
-			 %   needed. This will limit the memory use at the maximum; 
-			 %
-			 %   Usage:
-			 %      md=activaterepository(md)
-
-
-			 %for each field of the model, we create an equivalent variable, with the same name, which we 
-			 %save in the repository. 
-
-			 modelfields=fields(md);
-
-			 %then save each of the big fields in their respective files
-			 for i=1:length(modelfields),
-				 field=char(modelfields(i));
-				 eval([field '=md.' field ';']);
-
-				 %anything NaN, will be saved 
-
-				 if isnumeric(eval(field)),
-					 if isnan(eval(field)),
-						 if exist(['.' md.repository '.' field],'file');
-							 error(['Repository ' md.repository ' already exists! Delete all files starting with .' md.repository ' in your directory, before activating a new repository']);
-						 end
-						 eval(['save .' md.repository '.' field ' ' field]);
-					 end
-				 end
-			 end
-		 end % }}}
-		 function disp(md) % {{{1
-			 disp(sprintf('\n%s = \n',inputname(1)));
-
-			 %name
-			 disp(sprintf('   Name:'));
-			 if isempty(md.name),
-				 disp(sprintf('      N/A'));
-			 else
-				 disp(sprintf('      %s',md.name));
-			 end
-
-			 %notes: md.notes is a cell array of strings
-			 disp(sprintf('   Notes:'));
-			 if isempty(md.notes),
-				 disp(sprintf('      N/A'));
-			 else
-				 if iscell(md.notes),
-					 strings=md.notes;
-					 for i=1:length(strings),
-						 string=char(strings(i));
-						 if i==1,
-							 disp(sprintf('      %s',string));
-						 else
-							 disp(sprintf('      %s',string));
-						 end
-					 end
-				 elseif ischar(md.notes),
-					 disp(sprintf('      %s',md.notes));
-				 else
-					 disp('display error message: ''notes'' for model is neither a cell array of string, nor a string. Cannot display.');
-				 end
-			 end
-
-			 %repository
-			 if ~isempty(md.repository),
-				 disp(sprintf('   Repository: %s',md.repository));
-			 end
-
-			 %diverse
-			 disp(sprintf('   Mesh: %s',md.repository));
-			 disp(sprintf('%s%s%s','      Mesh: type ''',inputname(1),'.mesh'' to display'));
-			 disp(sprintf('   Parameters: %s',md.repository));
-			 disp(sprintf('%s%s%s','      Boundary conditions: type ''',inputname(1),'.bc'' to display'));
-			 disp(sprintf('%s%s%s','      Observations: type ''',inputname(1),'.obs'' to display'));
-			 disp(sprintf('%s%s%s','      Materials: type ''',inputname(1),'.mat'' to display'));
-			 disp(sprintf('%s%s%s','      Parameters: type ''',inputname(1),'.par'' to display'));
-			 disp(sprintf('%s%s%s','      ExpPar: type ''',inputname(1),'.exppar'' to display'));
-			 disp(sprintf('   Solution parameters: %s',md.repository));
-			 disp(sprintf('%s%s%s','      Qmu: type ''',inputname(1),'.qmu'' to display'));
-			 disp(sprintf('%s%s%s','      FLAIM: type ''',inputname(1),'.flaim'' to display'));
-			 disp(sprintf('%s%s%s','      Diagnostic solution parameters: type ''',inputname(1),'.diagnostic'' to display'));
-			 disp(sprintf('%s%s%s','      Thermal solution parameters: type ''',inputname(1),'.thermal'' to display'));
-			 disp(sprintf('%s%s%s','      Prognostic solution parameters: type ''',inputname(1),'.prognostic'' to display'));
-			 disp(sprintf('%s%s%s','      Transient solution parameters: type ''',inputname(1),'.transient'' to display'));
-			 disp(sprintf('%s%s%s','      Control solution parameters: type ''',inputname(1),'.control'' to display'));
-			 disp(sprintf('   Parallel:'));
-			 disp(sprintf('%s%s%s','      Parallel options: type ''',inputname(1),'.parallel'' to display'));
-			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
-			 disp(sprintf('   Results:'));
-			 disp(sprintf('%s%s%s','      Results: type ''',inputname(1),'.res'' to display'));
-			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
-		 end
-		 %}}}
-		 function mdstruct=struct(md) % {{{1
-			 mdstruct=struct();
-			 modelprops=properties(md);
-
-			 for i=1:length(modelprops),
-				 mdstruct.(modelprops{i})=md.(modelprops{i});
-			 end
-		 end
-		 %}}}
-		 function md=setdefaultparameters(md) % {{{1
-		 %SETDEFAULTPARAMETERS - plug default parameters onto model
-		 %
-		 %   Although the model parameterization should be done in
-		 %   the parameter file, some parameters are initialized here
-		 %   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;
-
-			 %firn density (kg/m^3)
-			 md.rho_firn=830;
-
-			 %ice heat capacity cp (J/kg/K)
-			 md.heatcapacity=2093;
-
-			 %ice latent heat of fusion L (J/kg)
-			 md.latentheat=3.34*10^5;
-
-			 %basal melting rate correction: 
-			 md.gl_melting_rate=0; 
-			 md.basal_melting_rate_correction_apply=0; %do not apply
-
-			 %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
-
-			 %mechanical residue convergence criterion norm(K(uold)uold - F)/norm(F)
-			 md.eps_res=10^-4; %from test233 in nightly runs
-
-			 %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;
-
-			 %maximum of non-linear iterations.
-			 md.max_nonlinear_iterations=100;
-
-			 %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=25;
-
-			 %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;
-
-			 %parameter used to print temporary results (convergence criterion,
-			 %current step,...)
-			 md.verbose=verbose('solution',true,'qmu',true,'control',true);
-
-			 %Stokes
-			 md.stokesreconditioning=10^13;
-			 md.shelf_dampening=0;
-
-			 %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;
-			 md.min_mechanical_constraints=0;
-
-			 %Transient parameters
-
-			 %time between 2 time steps
-			 md.dt=1/2;
-
-			 %number of time steps = md.ndt / md.dt
-			 md.ndt=10*md.dt;
-
-			 %time adaptation? 
-			 md.time_adapt=0;
-			 md.cfl_coefficient=1/2;
-
-			 %Hydrostatic adjustment
-			 md.hydrostatic_adjustment=AbsoluteEnum;
-			 %the artificial diffusivity is used in prognostic to avoid
-			 %numerical wiggles of the solution.
-			 md.artificial_diffusivity=1;
-
-			 %Discontinuous Galerkin is used in prognostic to avoid
-			 %numerical wiggles of the solution and conserve mass
-			 md.prognostic_DG=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;
-
-			 %Solution activated for transient runs. By default we do a
-			 %full analysis: Diagnostic, Prognostic and Thermal
-			 md.isprognostic=1;
-			 md.isdiagnostic=1;
-			 md.isthermal=1;
-
-			 %Control
-
-			 %parameter to be inferred by control methods (only
-			 %drag and B are supported yet)
-			 md.control_type=DragCoefficientEnum();
-
-			 %number of steps in the control methods
-			 md.nsteps=20;
-
-			 %maximum number of iteration in the optimization algorithm for
-			 %each step
-			 md.maxiter=20*ones(md.nsteps,1);
-
-			 %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);
-
-			 %several responses can be used:
-			 md.cm_responses=SurfaceAbsVelMisfitEnum*ones(md.nsteps,1);
-
-			 %cm_jump is used to speed up control method. When
-			 %misfit(1)/misfit(0) < md.cm_jump, we go directly to
-			 %the next step
-			 md.cm_jump=.7*ones(md.nsteps,1); %30 per cent decrement.
-
-			 %stop control solution at the gradient computation and return it? 
-			 md.cm_gradient=0;
-
-			 %eps_cm is a criteria to stop the control methods.
-			 %if J[n]-J[n-1]/J[n] < criteria, the control run stops
-			 %NaN if not applied
-			 md.eps_cm=NaN; %not activated
-
-			 %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);
-
-			 %grounding line migration: 
-			 md.gl_migration=NoneEnum;
-
-			 %How often to save results, default is 1 so save every step
-			 md.output_frequency=1;
-
-			 %Parallelisation parameters
-
-			 %cluster set as none for serial
-			 md.cluster=none;
-
-			 %this option can be activated to load automatically the results
-			 %onto the model after a parallel run by waiting for the lock file
-			 %N minutes that is generated once the solution has converged
-			 %0 to desactivate
-			 md.waitonlock=Inf;
-
-			 %prerun allocation cleanup
-			 md.alloc_cleanup=1;
-
-			 %set petsc options for different analysis
-			 md.petscoptions=petscoptions;
-			 md.petscoptions=addoptions(md.petscoptions,DiagnosticVertAnalysisEnum,mumpsoptions);
-
-			 %solution speed-up
-			 md.kff=1;
-
-			 %hydrology:  from Johnson's 2002 thesis, section 3.5.4
-			 md.hydro_n=.02;
-			 md.hydro_CR=2;
-			 md.hydro_p=2;
-			 md.hydro_q=1;
-			 md.hydro_gamma=1/(md.hydro_n*md.hydro_CR^md.hydro_p*(md.rho_water*md.g)^md.hydro_q);
-			 md.hydro_kn=0;
-
-			 %Rheology law: what is the temperature dependence of B with T
-			 %available: None, Paterson and Arrhenius
-			 md.rheology_law=PatersonEnum;
-
-			 %i/o:
-			 md.io_gather=1;
-		 end
-		 %}}}
-		 function md=subsasgn(md,index,val) % {{{1
-
-			 if ~isempty(md.repository),
-				 field=index(1).subs;
-				 if exist(['.' md.repository '.' field],'file'),
-					 if length(index)==1,
-						 %save val in corresponding repository
-						 eval([field '=val;']);
-						 eval(['save .' md.repository '.' field ' ' field]);
-					 else
-						 %load field from file
-						 load(['.' md.repository '.' field],'-mat');
-						 fieldval=eval(field);
-						 eval([field '=builtin(''subsasgn'',fieldval,index(2),val);']);
-						 eval(['save .' md.repository '.' field ' ' field]);
-					 end
-				 else
-					 md=builtin('subsasgn',md,index,val);
-				 end
-			 else
-				 md=builtin('subsasgn',md,index,val);
-			 end
-	 end %}}}
-		 function result=subsref(md,index) % {{{1
-
-			 if length(index)==1,
-				 index1=index(1);
-				 if(strcmp(index1.subs,'mesh')), displaymesh(md);return; end
-				 if(strcmp(index1.subs,'bc')), displaybc(md);return; end
-				 if(strcmp(index1.subs,'mat')), displaymaterials(md);return; end
-				 if(strcmp(index1.subs,'par')), displayparameters(md);return; end
-				 if(strcmp(index1.subs,'exppar')), displayexppar(md);return; end
-				 if(strcmp(index1.subs,'res')), displayresults(md);return; end
-				 if(strcmp(index1.subs,'obs')), displayobservations(md);return; end
-				 if(strcmp(index1.subs,'qmu')), displayqmu(md);return; end
-				 if(strcmp(index1.subs,'flaim')), displayflaim(md);return; end
-				 if(strcmp(index1.subs,'diagnostic')), displaydiagnostic(md);return; end
-				 if(strcmp(index1.subs,'prognostic')), displayprognostic(md);return; end
-				 if(strcmp(index1.subs,'thermal')), displaythermal(md);return; end
-				 if(strcmp(index1.subs,'transient')), displaytransient(md);return; end
-				 if(strcmp(index1.subs,'control')), displaycontrol(md);return; end
-				 if(strcmp(index1.subs,'mesh2grid')), displaymesh2grid(md);return; end
-				 if(strcmp(index1.subs,'parallel')), displayparallel(md);return; end
-
-				 if ~isempty(md.repository),
-					 %load variable needed, if necessary
-					 if exist(['.' md.repository '.' index1.subs],'file'),
-						 structure=load(['.' md.repository '.' index1.subs],'-mat');
-						 eval(['result=structure.' index1.subs ';']);
-					 else
-						 %the file does not exist, just return built in value
-						 result=builtin('subsref',md,index);
-					 end
-				 else
-					 %return built in value
-					 result=builtin('subsref',md,index);
-				 end
-			 else
-				 if ~isempty(md.repository),
-					 if exist(['.' md.repository '.' index(1).subs],'file'),
-						 %load variable needed
-						 variable=load(['.' md.repository '.' index(1).subs],'-mat');
-
-						 %use temporary structure.
-						 result=builtin('subsref',variable,index);
-					 else
-						 result=builtin('subsref',md,index);
-					 end
-				 else
-					 result=builtin('subsref',md,index);
-				 end
-			 end
-		 end %}}}
-	 end
- end
Index: /issm/trunk/src/m/classes/version/7.6/model.m
===================================================================
--- /issm/trunk/src/m/classes/version/7.6/model.m	(revision 8458)
+++ /issm/trunk/src/m/classes/version/7.6/model.m	(revision 8458)
@@ -0,0 +1,871 @@
+%MODEL class definition
+%
+%   Usage:
+%      md = model(varargin)
+
+classdef model
+    properties (SetAccess=public) %Model fields
+		 % {{{1
+		 %Careful here: no other class should be used as default value
+		 %this is a bug of matlab
+
+		 %Model general information
+		 notes='';
+		 name=''; 
+		 runtimename=''; %name used when running a parallel job
+		 counter=0; %1 mesh, 2 geography, 3 parameterize
+		 domainoutline=NaN;
+		 riftoutline=NaN;
+		 iceshelfoutline=NaN;
+		 icesheetoutline=NaN;
+		 wateroutline=NaN;
+		 parameterfile=NaN;
+		 runmefile=NaN;
+
+		 %Mesh
+		 bamg=struct();
+		 dim=0;
+		 numberofelements=0;
+		 numberofnodes=0;
+		 elements=NaN;
+		 elements_type=NaN;
+		 vertices_type=NaN;
+		 x=NaN;
+		 y=NaN;
+		 z=NaN;
+		 bed_slopex=NaN;
+		 bed_slopey=NaN;
+		 surface_slopex=NaN;
+		 surface_slopey=NaN;
+		 nodeconnectivity=NaN;
+		 elementconnectivity=NaN;
+		 edges=NaN;
+
+		 %I/O
+		 io_gather=NaN;
+
+		 %Initial 2d mesh 
+		 numberofelements2d=0;
+		 numberofnodes2d=0;
+		 elements2d=NaN;
+		 elements_type2d=NaN;
+		 vertices_type2d=NaN;
+		 x2d=NaN;
+		 y2d=NaN;
+		 z2d=NaN;
+
+		 %latlon of the coorinates
+		 lat=NaN;
+		 long=NaN;
+		 hemisphere=NaN;
+
+		 %Elements type
+		 ishutter=0;
+		 ismacayealpattyn=0;
+		 isstokes=0;
+
+		 %Elements
+		 elementonhutter=NaN;
+		 elementonmacayeal=NaN;
+		 elementonpattyn=NaN;
+		 elementonstokes=NaN;
+
+		 %Nodes
+		 nodeonhutter=NaN;
+		 nodeonmacayeal=NaN;
+		 nodeonpattyn=NaN;
+		 nodeonstokes=NaN;
+		 borderstokes=NaN;
+
+		 %Stokes
+		 stokesreconditioning=0;
+		 shelf_dampening=0;
+
+		 %Penalties
+		 penalties=NaN;
+		 penalty_offset=0;
+		 penalty_melting=0;
+		 penalty_lock=0;
+		 segments=NaN;
+		 segmentmarkers=NaN;
+		 rifts=NaN;
+		 riftinfo=NaN;
+		 riftproperties=NaN;
+		 numrifts=0;
+
+		 %Projections
+		 uppernodes=NaN;
+		 upperelements=NaN;
+		 lowerelements=NaN;
+		 lowernodes=NaN;
+
+		 %Extrusion
+		 numlayers=0;
+		 extrusionexponent=0;
+		 elementonbed=NaN;
+		 elementonsurface=NaN;
+		 nodeonbed=NaN;
+		 nodeonsurface=NaN;
+		 minh=0;
+		 firn_layer=NaN;
+
+		 %Extraction
+		 extractednodes=NaN;
+		 extractedelements=NaN;
+
+		 %Materials parameters
+		 rho_ice=0;
+		 rho_water=0;
+		 rho_firn=0;
+		 heatcapacity=0;
+		 latentheat=0;
+		 thermalconductivity=0;
+		 meltingpoint=0;
+		 beta=0;
+		 
+		 mixed_layer_capacity=0;
+		 thermal_exchange_velocity=0;
+		 min_thermal_constraints=0;
+		 min_mechanical_constraints=0;
+		 stabilize_constraints=0;
+
+		 %Physical parameters
+		 g=0;
+		 yts=0;
+		 drag_type=0;
+		 drag_coefficient=NaN;
+		 drag_p=NaN;
+		 drag_q=NaN;
+		 rheology_B=NaN;
+		 rheology_n=NaN;
+		 rheology_law=0;
+
+		 %Geometrical parameters
+		 elementoniceshelf=NaN;
+		 elementonicesheet=NaN;
+		 elementonwater=NaN;
+		 elementonnuna=NaN;
+		 nodeoniceshelf=NaN;
+		 nodeonicesheet=NaN;
+		 nodeonwater=NaN;
+		 nodeonnuna=NaN;
+		 surface=NaN;
+		 thickness=NaN;
+		 thickness_coeff=NaN;
+		 bed=NaN;
+		 bathymetry=NaN;
+		 mask=NaN;
+
+		 %Boundary conditions
+		 nodeonboundary=NaN;
+		 pressureload=NaN;
+		 spcvelocity=NaN;
+		 spctemperature=NaN;
+		 spcthickness=NaN;
+		 spcwatercolumn=NaN;
+		 diagnostic_ref=NaN;
+
+		 %Observations 
+		 vx_obs=NaN;
+		 vy_obs=NaN;
+		 vel_obs=NaN;
+		 vx_obs_raw=NaN;
+		 vy_obs_raw=NaN;
+		 vx_bal=NaN;
+		 vy_bal=NaN;
+		 vel_bal=NaN;
+		 vel_obs_raw=NaN;
+		 surface_accumulation_rate=NaN;
+		 surface_ablation_rate=NaN;
+		 surface_mass_balance=NaN;
+		 dhdt=NaN;
+		 geothermalflux=NaN;
+		 observed_temperature=NaN;
+		 thickness_obs=NaN;
+
+		 %Forcings
+		 forcings=struct();
+
+		 %Statics parameters
+		 eps_res=0;
+		 eps_rel=0;
+		 eps_abs=0;
+		 max_nonlinear_iterations=0;
+		 sparsity=0;
+		 connectivity=0;
+		 lowmem=0;
+		 viscosity_overshoot=0;
+
+		 %Transient parameters
+		 dt=0;
+		 ndt=0;
+		 time_adapt=0;
+		 cfl_coefficient=1/2;
+		 adaptative_cfl=0;
+		 artificial_diffusivity=0;
+		 prognostic_DG=0;
+		 deltaH=0;
+		 DeltaH=0;
+		 deltaT=0;
+		 DeltaT=0;
+		 timeacc=0;
+		 timedec=0;
+		 hydrostatic_adjustment=0;
+
+		 %Transient flags
+		 isprognostic=0;
+		 isdiagnostic=0;
+		 isthermal=0;
+
+		 %Control
+		 control_analysis=0;
+		 control_type=0;
+		 num_control_type=0;
+		 weights=[];
+		 nsteps=0;
+		 maxiter=[];
+		 cm_responses=[];
+		 tolx=0;
+		 optscal=[];
+		 eps_cm=0;
+		 cm_noisedmp=0;
+		 cm_min=NaN;
+		 cm_max=NaN;
+		 cm_jump=[];
+		 cm_gradient=0;
+		 epsvel=0;
+		 meanvel=0;
+
+		 %Output parameters
+		 numoutput=0;
+		 parameteroutput={};
+		 viscousheating=NaN;
+		 pressure_elem=NaN;
+		 stress=NaN;
+		 stress_surface=NaN;
+		 stress_bed=NaN;
+		 deviatoricstress=NaN;
+		 strainrate=NaN;
+
+		 %debugging
+		 verbose=0;
+		 verbose_binary=0;
+		 element_debug=0;
+		 element_debugid=NaN;
+		 mem_debug=0;
+		 gprof=0;
+		 memoryleaks=NaN;
+
+		 %Results fields
+		 output_frequency=0;
+		 inputfilename='';
+		 outputfilename='';
+		 results=struct();
+		 vx=NaN;
+		 vy=NaN;
+		 vz=NaN;
+		 vel=NaN;
+		 temperature=NaN; %temperature solution vector
+		 gl_melting_rate=NaN;
+		 basal_melting_rate=NaN;
+		 basal_melting_rate_correction=NaN;
+		 basal_melting_rate_correction_apply=0;
+		 pressure=NaN;
+		 
+		 %hydrology
+		 watercolumn=NaN;
+		 hydro_n=NaN;
+		 hydro_CR=NaN;
+		 hydro_p=NaN;
+		 hydro_q=NaN;
+		 hydro_gamma=0;
+		 hydro_kn=0;
+		 
+		 %Parallelisation parameters
+		 cluster=NaN;
+		 outlog='';
+		 errlog='';
+		 alloc_cleanup=1;
+		 waitonlock=0;
+
+		 %mesh2grid
+		 mesh2grid_parameters={};
+		 mesh2grid_data=NaN;
+		 mesh2grid_xm=NaN;
+		 mesh2grid_ym=NaN;
+
+		 %dummy
+		 dummy=NaN;
+
+		 %mesh 
+		 mesh_domainoutline='';
+		 mesh_riftoutline='';
+		 mesh_resolution=NaN;
+
+		 %PETSc and MATLAB solver string
+		 petscoptions=NaN;
+
+		 %Analysis
+		 solution_type='';
+
+		 %management of large models
+		 repository='';
+
+		 %radar power images
+		 sarpwr=NaN;
+		 sarxm=NaN;
+		 sarym=NaN;
+
+		 %qmu
+		 variables=struct();
+		 responses=struct();
+		 qmu_method=struct();
+		 qmu_params=struct();
+		 dakotaresults=struct();
+		 dakotain='';
+		 dakotaout='';
+		 dakotadat='';
+		 qmu_analysis=0;
+		 part=[];
+		 npart=0;
+		 numberofvariables=0;
+		 numberofresponses=0;
+		 variabledescriptors=NaN;
+		 responsedescriptors=NaN;
+		 numvariabledescriptors=0;
+		 numresponsedescriptors=0;
+		 qmu_mass_flux_num_profiles=NaN;
+		 qmu_mass_flux_profile_directory=NaN;
+		 qmu_mass_flux_profiles=NaN;
+		 qmu_mass_flux_segments=NaN;
+		 qmu_relax=0;
+		 qmu_save_femmodel=0;
+
+		 %flaim
+		 fm_tracks='';
+		 fm_flightreqs=struct();
+		 fm_criterion=[];
+		 fm_solution='';
+		 fm_quality=0;
+
+		 %grounding line migration: 
+		 gl_migration=NaN;
+
+		 %solution parameters
+		 kff=0;
+		 loadresults=0;
+
+		 %partitioner:
+		 adjacency=[];
+		 vwgt=[];
+		 %}}}
+	 end
+	 properties(Dependent,Hidden) %Deprecated property
+		 % {{{
+		 % We must keep here ALL deprecated properties
+		 % This ensures that when a model object from a previous version
+		 % all fields are recovered.
+		 % loadobj is then used to assign old properties to new properties
+		 % But if we don't keep the old properties here, they would not be
+		 % recovered when the model is loaded
+		 % KEEP IN CHRONOLOGICAL ORDER (OLDER FIRST)
+		 drag;
+		 p;
+		 q;
+		 B;
+		 n;
+		 melting;
+		 accumulation;
+		 type;
+		 numberofgrids;
+		 numberofgrids2d;
+		 gridonhutter;
+		 gridonmacayeal;
+		 gridonpattyn;
+		 gridonstokes;
+		 uppergrids;
+		 lowergrids;
+		 gridonbed;
+		 gridonsurface;
+		 extractedgrids;
+		 gridoniceshelf;
+		 gridonicesheet;
+		 gridonwater;
+		 gridonnuna;
+		 gridonboundary;
+		 melting_rate;
+		 accumulation_rate;
+		 % }}}
+	 end
+	 methods
+		 function md=model(varargin) % {{{1
+
+			 switch nargin
+				 case 0
+
+					 %set default parameters
+					 md=setdefaultparameters(md);
+
+				 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
+		 end
+		 %}}}
+		 function md=activaterepository(md) % {{{1
+			 %ACTIVATEREPOSITORY - save the model fields separately
+			 %
+			 %   This function is required when a matlab session 
+			 %   is running low on memory, or when a model is very large. 
+			 %   By activating a repository (for example, 'repository.model'), 
+			 %   the model fields will each separately be saved and accessed when 
+			 %   needed. This will limit the memory use at the maximum; 
+			 %
+			 %   Usage:
+			 %      md=activaterepository(md)
+
+
+			 %for each field of the model, we create an equivalent variable, with the same name, which we 
+			 %save in the repository. 
+
+			 modelfields=fields(md);
+
+			 %then save each of the big fields in their respective files
+			 for i=1:length(modelfields),
+				 field=char(modelfields(i));
+				 eval([field '=md.' field ';']);
+
+				 %anything NaN, will be saved 
+
+				 if isnumeric(eval(field)),
+					 if isnan(eval(field)),
+						 if exist(['.' md.repository '.' field],'file');
+							 error(['Repository ' md.repository ' already exists! Delete all files starting with .' md.repository ' in your directory, before activating a new repository']);
+						 end
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 end
+				 end
+			 end
+		 end % }}}
+		 function disp(md) % {{{1
+			 disp(sprintf('\n%s = \n',inputname(1)));
+
+			 %name
+			 disp(sprintf('   Name:'));
+			 if isempty(md.name),
+				 disp(sprintf('      N/A'));
+			 else
+				 disp(sprintf('      %s',md.name));
+			 end
+
+			 %notes: md.notes is a cell array of strings
+			 disp(sprintf('   Notes:'));
+			 if isempty(md.notes),
+				 disp(sprintf('      N/A'));
+			 else
+				 if iscell(md.notes),
+					 strings=md.notes;
+					 for i=1:length(strings),
+						 string=char(strings(i));
+						 if i==1,
+							 disp(sprintf('      %s',string));
+						 else
+							 disp(sprintf('      %s',string));
+						 end
+					 end
+				 elseif ischar(md.notes),
+					 disp(sprintf('      %s',md.notes));
+				 else
+					 disp('display error message: ''notes'' for model is neither a cell array of string, nor a string. Cannot display.');
+				 end
+			 end
+
+			 %repository
+			 if ~isempty(md.repository),
+				 disp(sprintf('   Repository: %s',md.repository));
+			 end
+
+			 %diverse
+			 disp(sprintf('   Mesh: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Mesh: type ''',inputname(1),'.mesh'' to display'));
+			 disp(sprintf('   Parameters: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Boundary conditions: type ''',inputname(1),'.bc'' to display'));
+			 disp(sprintf('%s%s%s','      Observations: type ''',inputname(1),'.obs'' to display'));
+			 disp(sprintf('%s%s%s','      Materials: type ''',inputname(1),'.mat'' to display'));
+			 disp(sprintf('%s%s%s','      Parameters: type ''',inputname(1),'.par'' to display'));
+			 disp(sprintf('%s%s%s','      ExpPar: type ''',inputname(1),'.exppar'' to display'));
+			 disp(sprintf('   Solution parameters: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Qmu: type ''',inputname(1),'.qmu'' to display'));
+			 disp(sprintf('%s%s%s','      FLAIM: type ''',inputname(1),'.flaim'' to display'));
+			 disp(sprintf('%s%s%s','      Diagnostic solution parameters: type ''',inputname(1),'.diagnostic'' to display'));
+			 disp(sprintf('%s%s%s','      Thermal solution parameters: type ''',inputname(1),'.thermal'' to display'));
+			 disp(sprintf('%s%s%s','      Prognostic solution parameters: type ''',inputname(1),'.prognostic'' to display'));
+			 disp(sprintf('%s%s%s','      Transient solution parameters: type ''',inputname(1),'.transient'' to display'));
+			 disp(sprintf('%s%s%s','      Control solution parameters: type ''',inputname(1),'.control'' to display'));
+			 disp(sprintf('   Parallel:'));
+			 disp(sprintf('%s%s%s','      Parallel options: type ''',inputname(1),'.parallel'' to display'));
+			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
+			 disp(sprintf('   Results:'));
+			 disp(sprintf('%s%s%s','      Results: type ''',inputname(1),'.res'' to display'));
+			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
+		 end
+		 %}}}
+		 function mdstruct=struct(md) % {{{1
+			 mdstruct=struct();
+			 modelprops=properties(md);
+
+			 for i=1:length(modelprops),
+				 mdstruct.(modelprops{i})=md.(modelprops{i});
+			 end
+		 end
+		 %}}}
+		 function md=setdefaultparameters(md) % {{{1
+		 %SETDEFAULTPARAMETERS - plug default parameters onto model
+		 %
+		 %   Although the model parameterization should be done in
+		 %   the parameter file, some parameters are initialized here
+		 %   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;
+
+			 %firn density (kg/m^3)
+			 md.rho_firn=830;
+
+			 %ice heat capacity cp (J/kg/K)
+			 md.heatcapacity=2093;
+
+			 %ice latent heat of fusion L (J/kg)
+			 md.latentheat=3.34*10^5;
+
+			 %basal melting rate correction: 
+			 md.gl_melting_rate=0; 
+			 md.basal_melting_rate_correction_apply=0; %do not apply
+
+			 %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
+
+			 %mechanical residue convergence criterion norm(K(uold)uold - F)/norm(F)
+			 md.eps_res=10^-4; %from test233 in nightly runs
+
+			 %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;
+
+			 %maximum of non-linear iterations.
+			 md.max_nonlinear_iterations=100;
+
+			 %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=25;
+
+			 %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;
+
+			 %parameter used to print temporary results (convergence criterion,
+			 %current step,...)
+			 md.verbose=verbose('solution',true,'qmu',true,'control',true);
+
+			 %Stokes
+			 md.stokesreconditioning=10^13;
+			 md.shelf_dampening=0;
+
+			 %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;
+			 md.min_mechanical_constraints=0;
+
+			 %Transient parameters
+
+			 %time between 2 time steps
+			 md.dt=1/2;
+
+			 %number of time steps = md.ndt / md.dt
+			 md.ndt=10*md.dt;
+
+			 %time adaptation? 
+			 md.time_adapt=0;
+			 md.cfl_coefficient=1/2;
+
+			 %Hydrostatic adjustment
+			 md.hydrostatic_adjustment=AbsoluteEnum;
+			 %the artificial diffusivity is used in prognostic to avoid
+			 %numerical wiggles of the solution.
+			 md.artificial_diffusivity=1;
+
+			 %Discontinuous Galerkin is used in prognostic to avoid
+			 %numerical wiggles of the solution and conserve mass
+			 md.prognostic_DG=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;
+
+			 %Solution activated for transient runs. By default we do a
+			 %full analysis: Diagnostic, Prognostic and Thermal
+			 md.isprognostic=1;
+			 md.isdiagnostic=1;
+			 md.isthermal=1;
+
+			 %Control
+
+			 %parameter to be inferred by control methods (only
+			 %drag and B are supported yet)
+			 md.control_type=DragCoefficientEnum();
+
+			 %number of steps in the control methods
+			 md.nsteps=20;
+
+			 %maximum number of iteration in the optimization algorithm for
+			 %each step
+			 md.maxiter=20*ones(md.nsteps,1);
+
+			 %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);
+
+			 %several responses can be used:
+			 md.cm_responses=SurfaceAbsVelMisfitEnum*ones(md.nsteps,1);
+
+			 %cm_jump is used to speed up control method. When
+			 %misfit(1)/misfit(0) < md.cm_jump, we go directly to
+			 %the next step
+			 md.cm_jump=.7*ones(md.nsteps,1); %30 per cent decrement.
+
+			 %stop control solution at the gradient computation and return it? 
+			 md.cm_gradient=0;
+
+			 %eps_cm is a criteria to stop the control methods.
+			 %if J[n]-J[n-1]/J[n] < criteria, the control run stops
+			 %NaN if not applied
+			 md.eps_cm=NaN; %not activated
+
+			 %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);
+
+			 %grounding line migration: 
+			 md.gl_migration=NoneEnum;
+
+			 %How often to save results, default is 1 so save every step
+			 md.output_frequency=1;
+
+			 %Parallelisation parameters
+
+			 %cluster set as none for serial
+			 md.cluster=none;
+
+			 %this option can be activated to load automatically the results
+			 %onto the model after a parallel run by waiting for the lock file
+			 %N minutes that is generated once the solution has converged
+			 %0 to desactivate
+			 md.waitonlock=Inf;
+
+			 %prerun allocation cleanup
+			 md.alloc_cleanup=1;
+
+			 %set petsc options for different analysis
+			 md.petscoptions=petscoptions;
+			 md.petscoptions=addoptions(md.petscoptions,DiagnosticVertAnalysisEnum,mumpsoptions);
+
+			 %solution speed-up
+			 md.kff=1;
+
+			 %hydrology:  from Johnson's 2002 thesis, section 3.5.4
+			 md.hydro_n=.02;
+			 md.hydro_CR=2;
+			 md.hydro_p=2;
+			 md.hydro_q=1;
+			 md.hydro_gamma=1/(md.hydro_n*md.hydro_CR^md.hydro_p*(md.rho_water*md.g)^md.hydro_q);
+			 md.hydro_kn=0;
+
+			 %Rheology law: what is the temperature dependence of B with T
+			 %available: None, Paterson and Arrhenius
+			 md.rheology_law=PatersonEnum;
+
+			 %i/o:
+			 md.io_gather=1;
+		 end
+		 %}}}
+		 function md=subsasgn(md,index,val) % {{{1
+
+			 if ~isempty(md.repository),
+				 field=index(1).subs;
+				 if exist(['.' md.repository '.' field],'file'),
+					 if length(index)==1,
+						 %save val in corresponding repository
+						 eval([field '=val;']);
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 else
+						 %load field from file
+						 load(['.' md.repository '.' field],'-mat');
+						 fieldval=eval(field);
+						 eval([field '=builtin(''subsasgn'',fieldval,index(2),val);']);
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 end
+				 else
+					 md=builtin('subsasgn',md,index,val);
+				 end
+			 else
+				 md=builtin('subsasgn',md,index,val);
+			 end
+	 end %}}}
+		 function result=subsref(md,index) % {{{1
+
+			 if length(index)==1,
+				 index1=index(1);
+				 if(strcmp(index1.subs,'mesh')), displaymesh(md);return; end
+				 if(strcmp(index1.subs,'bc')), displaybc(md);return; end
+				 if(strcmp(index1.subs,'mat')), displaymaterials(md);return; end
+				 if(strcmp(index1.subs,'par')), displayparameters(md);return; end
+				 if(strcmp(index1.subs,'exppar')), displayexppar(md);return; end
+				 if(strcmp(index1.subs,'res')), displayresults(md);return; end
+				 if(strcmp(index1.subs,'obs')), displayobservations(md);return; end
+				 if(strcmp(index1.subs,'qmu')), displayqmu(md);return; end
+				 if(strcmp(index1.subs,'flaim')), displayflaim(md);return; end
+				 if(strcmp(index1.subs,'diagnostic')), displaydiagnostic(md);return; end
+				 if(strcmp(index1.subs,'prognostic')), displayprognostic(md);return; end
+				 if(strcmp(index1.subs,'thermal')), displaythermal(md);return; end
+				 if(strcmp(index1.subs,'transient')), displaytransient(md);return; end
+				 if(strcmp(index1.subs,'control')), displaycontrol(md);return; end
+				 if(strcmp(index1.subs,'mesh2grid')), displaymesh2grid(md);return; end
+				 if(strcmp(index1.subs,'parallel')), displayparallel(md);return; end
+
+				 if ~isempty(md.repository),
+					 %load variable needed, if necessary
+					 if exist(['.' md.repository '.' index1.subs],'file'),
+						 structure=load(['.' md.repository '.' index1.subs],'-mat');
+						 eval(['result=structure.' index1.subs ';']);
+					 else
+						 %the file does not exist, just return built in value
+						 result=builtin('subsref',md,index);
+					 end
+				 else
+					 %return built in value
+					 result=builtin('subsref',md,index);
+				 end
+			 else
+				 if ~isempty(md.repository),
+					 if exist(['.' md.repository '.' index(1).subs],'file'),
+						 %load variable needed
+						 variable=load(['.' md.repository '.' index(1).subs],'-mat');
+
+						 %use temporary structure.
+						 result=builtin('subsref',variable,index);
+					 else
+						 result=builtin('subsref',md,index);
+					 end
+				 else
+					 result=builtin('subsref',md,index);
+				 end
+			 end
+		 end %}}}
+	 end
+ end
Index: /issm/trunk/src/m/classes/version/7.7/model.m
===================================================================
--- /issm/trunk/src/m/classes/version/7.7/model.m	(revision 8458)
+++ /issm/trunk/src/m/classes/version/7.7/model.m	(revision 8458)
@@ -0,0 +1,943 @@
+%MODEL class definition
+%
+%   Usage:
+%      md = model(varargin)
+
+classdef model
+    properties (SetAccess=public) %Model fields
+		 % {{{1
+		 %Careful here: no other class should be used as default value
+		 %this is a bug of matlab
+
+		 %Model general information
+		 notes='';
+		 name=''; 
+		 runtimename=''; %name used when running a parallel job
+		 counter=0; %1 mesh, 2 geography, 3 parameterize
+		 domainoutline=NaN;
+		 riftoutline=NaN;
+		 iceshelfoutline=NaN;
+		 icesheetoutline=NaN;
+		 wateroutline=NaN;
+		 parameterfile=NaN;
+		 runmefile=NaN;
+
+		 %Mesh
+		 bamg=struct();
+		 dim=0;
+		 numberofelements=0;
+		 numberofnodes=0;
+		 elements=NaN;
+		 elements_type=NaN;
+		 vertices_type=NaN;
+		 x=NaN;
+		 y=NaN;
+		 z=NaN;
+		 bed_slopex=NaN;
+		 bed_slopey=NaN;
+		 surface_slopex=NaN;
+		 surface_slopey=NaN;
+		 nodeconnectivity=NaN;
+		 elementconnectivity=NaN;
+		 edges=NaN;
+
+		 %I/O
+		 io_gather=NaN;
+
+		 %Initial 2d mesh 
+		 numberofelements2d=0;
+		 numberofnodes2d=0;
+		 elements2d=NaN;
+		 elements_type2d=NaN;
+		 vertices_type2d=NaN;
+		 x2d=NaN;
+		 y2d=NaN;
+		 z2d=NaN;
+
+		 %latlon of the coorinates
+		 lat=NaN;
+		 long=NaN;
+		 hemisphere=NaN;
+
+		 %Elements type
+		 ishutter=0;
+		 ismacayealpattyn=0;
+		 isstokes=0;
+
+		 %Elements
+		 elementonhutter=NaN;
+		 elementonmacayeal=NaN;
+		 elementonpattyn=NaN;
+		 elementonstokes=NaN;
+
+		 %Nodes
+		 nodeonhutter=NaN;
+		 nodeonmacayeal=NaN;
+		 nodeonpattyn=NaN;
+		 nodeonstokes=NaN;
+		 borderstokes=NaN;
+
+		 %Stokes
+		 stokesreconditioning=0;
+		 shelf_dampening=0;
+
+		 %Penalties
+		 penalties=NaN;
+		 penalty_offset=0;
+		 penalty_melting=0;
+		 penalty_lock=0;
+		 segments=NaN;
+		 segmentmarkers=NaN;
+		 rifts=NaN;
+		 riftinfo=NaN;
+		 riftproperties=NaN;
+		 numrifts=0;
+
+		 %Projections
+		 uppernodes=NaN;
+		 upperelements=NaN;
+		 lowerelements=NaN;
+		 lowernodes=NaN;
+
+		 %Extrusion
+		 numlayers=0;
+		 extrusionexponent=0;
+		 elementonbed=NaN;
+		 elementonsurface=NaN;
+		 nodeonbed=NaN;
+		 nodeonsurface=NaN;
+		 minh=0;
+		 firn_layer=NaN;
+
+		 %Extraction
+		 extractednodes=NaN;
+		 extractedelements=NaN;
+
+		 %Materials parameters
+		 rho_ice=0;
+		 rho_water=0;
+		 rho_firn=0;
+		 heatcapacity=0;
+		 latentheat=0;
+		 thermalconductivity=0;
+		 meltingpoint=0;
+		 beta=0;
+		 
+		 mixed_layer_capacity=0;
+		 thermal_exchange_velocity=0;
+		 min_thermal_constraints=0;
+		 min_mechanical_constraints=0;
+		 stabilize_constraints=0;
+
+		 %Physical parameters
+		 g=0;
+		 yts=0;
+		 drag_type=0;
+		 drag_coefficient=NaN;
+		 drag_p=NaN;
+		 drag_q=NaN;
+		 rheology_B=NaN;
+		 rheology_n=NaN;
+		 rheology_law=0;
+
+		 %Geometrical parameters
+		 elementoniceshelf=NaN;
+		 elementonicesheet=NaN;
+		 elementonwater=NaN;
+		 elementonnuna=NaN;
+		 nodeoniceshelf=NaN;
+		 nodeonicesheet=NaN;
+		 nodeonwater=NaN;
+		 nodeonnuna=NaN;
+		 surface=NaN;
+		 thickness=NaN;
+		 thickness_coeff=NaN;
+		 bed=NaN;
+		 bathymetry=NaN;
+		 mask=NaN;
+
+		 %Boundary conditions
+		 nodeonboundary=NaN;
+		 pressureload=NaN;
+		 spcvelocity=NaN;
+		 spctemperature=NaN;
+		 spcthickness=NaN;
+		 spcwatercolumn=NaN;
+		 diagnostic_ref=NaN;
+
+		 %Observations 
+		 vx_obs=NaN;
+		 vy_obs=NaN;
+		 vel_obs=NaN;
+		 vx_obs_raw=NaN;
+		 vy_obs_raw=NaN;
+		 vx_bal=NaN;
+		 vy_bal=NaN;
+		 vel_bal=NaN;
+		 vel_obs_raw=NaN;
+		 surface_accumulation_rate=NaN;
+		 surface_ablation_rate=NaN;
+		 surface_mass_balance=NaN;
+		 dhdt=NaN;
+		 geothermalflux=NaN;
+		 observed_temperature=NaN;
+		 thickness_obs=NaN;
+
+		 %Forcings
+		 forcings=struct();
+
+		 %Statics parameters
+		 eps_res=0;
+		 eps_rel=0;
+		 eps_abs=0;
+		 max_nonlinear_iterations=0;
+		 sparsity=0;
+		 connectivity=0;
+		 lowmem=0;
+		 viscosity_overshoot=0;
+
+		 %Transient parameters
+		 dt=0;
+		 ndt=0;
+		 time_adapt=0;
+		 cfl_coefficient=1/2;
+		 adaptative_cfl=0;
+		 artificial_diffusivity=0;
+		 prognostic_DG=0;
+		 deltaH=0;
+		 DeltaH=0;
+		 deltaT=0;
+		 DeltaT=0;
+		 timeacc=0;
+		 timedec=0;
+		 hydrostatic_adjustment=0;
+
+		 %Transient flags
+		 isprognostic=0;
+		 isdiagnostic=0;
+		 isthermal=0;
+
+		 %Control
+		 control_analysis=0;
+		 control_type=0;
+		 num_control_type=0;
+		 weights=[];
+		 nsteps=0;
+		 maxiter=[];
+		 cm_responses=[];
+		 tolx=0;
+		 optscal=[];
+		 eps_cm=0;
+		 cm_noisedmp=0;
+		 cm_min=NaN;
+		 cm_max=NaN;
+		 cm_jump=[];
+		 cm_gradient=0;
+		 epsvel=0;
+		 meanvel=0;
+
+		 %Output parameters
+		 numoutput=0;
+		 parameteroutput={};
+		 viscousheating=NaN;
+		 pressure_elem=NaN;
+		 stress=NaN;
+		 stress_surface=NaN;
+		 stress_bed=NaN;
+		 deviatoricstress=NaN;
+		 strainrate=NaN;
+
+		 %debugging
+		 verbose=0;
+		 verbose_binary=0;
+		 element_debug=0;
+		 element_debugid=NaN;
+		 mem_debug=0;
+		 gprof=0;
+		 memoryleaks=NaN;
+
+		 %Results fields
+		 output_frequency=0;
+		 inputfilename='';
+		 outputfilename='';
+		 results=struct();
+		 vx=NaN;
+		 vy=NaN;
+		 vz=NaN;
+		 vel=NaN;
+		 temperature=NaN; %temperature solution vector
+		 gl_melting_rate=NaN;
+		 basal_melting_rate=NaN;
+		 basal_melting_rate_correction=NaN;
+		 basal_melting_rate_correction_apply=0;
+		 pressure=NaN;
+		 
+		 %hydrology
+		 watercolumn=NaN;
+		 hydro_n=NaN;
+		 hydro_CR=NaN;
+		 hydro_p=NaN;
+		 hydro_q=NaN;
+		 hydro_gamma=0;
+		 hydro_kn=0;
+		 
+		 %Parallelisation parameters
+		 cluster=NaN;
+		 outlog='';
+		 errlog='';
+		 alloc_cleanup=1;
+		 waitonlock=0;
+
+		 %mesh2grid
+		 mesh2grid_parameters={};
+		 mesh2grid_data=NaN;
+		 mesh2grid_xm=NaN;
+		 mesh2grid_ym=NaN;
+
+		 %dummy
+		 dummy=NaN;
+
+		 %mesh 
+		 mesh_domainoutline='';
+		 mesh_riftoutline='';
+		 mesh_resolution=NaN;
+
+		 %PETSc and MATLAB solver string
+		 petscoptions=NaN;
+
+		 %Analysis
+		 solution_type='';
+
+		 %management of large models
+		 repository='';
+
+		 %radar power images
+		 sarpwr=NaN;
+		 sarxm=NaN;
+		 sarym=NaN;
+
+		 %qmu
+		 variables=struct();
+		 responses=struct();
+		 qmu_method=struct();
+		 qmu_params=struct();
+		 dakotaresults=struct();
+		 dakotain='';
+		 dakotaout='';
+		 dakotadat='';
+		 qmu_analysis=0;
+		 part=[];
+		 npart=0;
+		 numberofvariables=0;
+		 numberofresponses=0;
+		 variabledescriptors=NaN;
+		 responsedescriptors=NaN;
+		 numvariabledescriptors=0;
+		 numresponsedescriptors=0;
+		 qmu_mass_flux_num_profiles=NaN;
+		 qmu_mass_flux_profile_directory=NaN;
+		 qmu_mass_flux_profiles=NaN;
+		 qmu_mass_flux_segments=NaN;
+		 qmu_relax=0;
+		 qmu_save_femmodel=0;
+
+		 %flaim
+		 fm_tracks='';
+		 fm_flightreqs=struct();
+		 fm_criterion=[];
+		 fm_solution='';
+		 fm_quality=0;
+
+		 %grounding line migration: 
+		 gl_migration=NaN;
+
+		 %solution parameters
+		 kff=0;
+		 loadresults=0;
+
+		 %partitioner:
+		 adjacency=[];
+		 vwgt=[];
+		 %}}}
+	 end
+	 properties(Dependent,Hidden) %Deprecated property
+		 % {{{
+		 % We must keep here ALL deprecated properties
+		 % This ensures that when a model object from a previous version
+		 % all fields are recovered.
+		 % loadobj is then used to assign old properties to new properties
+		 % But if we don't keep the old properties here, they would not be
+		 % recovered when the model is loaded
+		 % KEEP IN CHRONOLOGICAL ORDER (OLDER FIRST)
+		 drag;
+		 p;
+		 q;
+		 B;
+		 n;
+		 melting;
+		 accumulation;
+		 type;
+		 numberofgrids;
+		 numberofgrids2d;
+		 gridonhutter;
+		 gridonmacayeal;
+		 gridonpattyn;
+		 gridonstokes;
+		 uppergrids;
+		 lowergrids;
+		 gridonbed;
+		 gridonsurface;
+		 extractedgrids;
+		 gridoniceshelf;
+		 gridonicesheet;
+		 gridonwater;
+		 gridonnuna;
+		 gridonboundary;
+		 melting_rate;
+		 accumulation_rate;
+		 % }}}
+	 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 isstruct(md)
+				 disp('Recovering model object from a previous version');
+
+				 %Ok the loaded model is a struct, initialize output and recover all fields
+				 structmd=md;
+				 md=model;
+
+				 structfields=fields(structmd);
+				 modelprops  =properties(md);
+
+				 for i=1:length(structfields),
+					 fieldname =structfields{i};
+					 fieldvalue=getfield(structmd,fieldname);
+					 if ismember(fieldname,modelprops),
+						 md=setfield(md,fieldname,fieldvalue);
+					 end
+				 end
+
+				 %Field name change
+				 if isfield(structmd,'drag'), md.drag_coefficient=structmd.drag; end
+				 if isfield(structmd,'p'), md.drag_p=structmd.p; end
+				 if isfield(structmd,'q'), md.drag_q=structmd.p; end
+				 if isfield(structmd,'B'), md.rheology_B=structmd.B; end
+				 if isfield(structmd,'n'), md.rheology_n=structmd.n; end
+				 if isfield(structmd,'melting'), md.basal_melting_rate=structmd.melting; end
+				 if isfield(structmd,'melting_rate'), md.basal_melting_rate=structmd.melting_rate; end
+				 if isfield(structmd,'accumulation'), md.surface_mass_balance=structmd.accumulation; end
+				 if isfield(structmd,'accumulation_rate'), md.surface_mass_balance=structmd.accumulation_rate; end
+				 if isfield(structmd,'numberofgrids'), md.numberofnodes=structmd.numberofgrids; end
+				 if isfield(structmd,'numberofgrids2d'), md.numberofnodes2d=structmd.numberofgrids2d; end
+				 if isfield(structmd,'gridonhutter'), md.nodeonhutter=structmd.gridonhutter; end
+				 if isfield(structmd,'gridonmacayeal'), md.nodeonmacayeal=structmd.gridonmacayeal; end
+				 if isfield(structmd,'gridonpattyn'), md.nodeonpattyn=structmd.gridonpattyn; end
+				 if isfield(structmd,'gridonstokes'), md.nodeonstokes=structmd.gridonstokes; end
+				 if isfield(structmd,'uppergrids'), md.uppernodes=structmd.uppergrids; end
+				 if isfield(structmd,'lowergrids'), md.lowernodes=structmd.lowergrids; end
+				 if isfield(structmd,'gridonbed'), md.nodeonbed=structmd.gridonbed; end
+				 if isfield(structmd,'gridonsurface'), md.nodeonsurface=structmd.gridonsurface; end
+				 if isfield(structmd,'extractedgrids'), md.extractednodes=structmd.extractedgrids; end
+				 if isfield(structmd,'gridoniceshelf'), md.nodeoniceshelf=structmd.gridoniceshelf; end
+				 if isfield(structmd,'gridonicesheet'), md.nodeonicesheet=structmd.gridonicesheet; end
+				 if isfield(structmd,'gridonwater'), md.nodeonwater=structmd.gridonwater; end
+				 if isfield(structmd,'gridonnuna'), md.nodeonnuna=structmd.gridonnuna; end
+				 if isfield(structmd,'gridonboundary'), md.nodeonboundary=structmd.gridonboundary; end
+
+				 %Field class change
+				 if (isfield(structmd,'type') & ischar(structmd.type)), 
+					 if strcmpi(structmd.type,'2d'), md.dim=2; end
+					 if strcmpi(structmd.type,'3d'), md.dim=3; end
+				 end
+				 if isnumeric(md.verbose), md.verbose=verbose; end
+
+				 %New fields
+				 if ~isfield(structmd,'upperelements');
+					 md.upperelements=transpose(1:md.numberofelements)+md.numberofelements2d;
+					 md.upperelements(end-md.numberofelements2d+1:end)=NaN;
+				 end
+				 if ~isfield(structmd,'lowerelements');
+					 md.lowerelements=transpose(1:md.numberofelements)-md.numberofelements2d;
+					 md.lowerelements(1:md.numberofelements2d)=NaN;
+				 end
+			 end
+		 end% }}}
+	 end
+	 methods
+		 function md=model(varargin) % {{{1
+
+			 switch nargin
+				 case 0
+
+					 %set default parameters
+					 md=setdefaultparameters(md);
+
+				 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
+		 end
+		 %}}}
+		 function md=activaterepository(md) % {{{1
+			 %ACTIVATEREPOSITORY - save the model fields separately
+			 %
+			 %   This function is required when a matlab session 
+			 %   is running low on memory, or when a model is very large. 
+			 %   By activating a repository (for example, 'repository.model'), 
+			 %   the model fields will each separately be saved and accessed when 
+			 %   needed. This will limit the memory use at the maximum; 
+			 %
+			 %   Usage:
+			 %      md=activaterepository(md)
+
+
+			 %for each field of the model, we create an equivalent variable, with the same name, which we 
+			 %save in the repository. 
+
+			 modelfields=fields(md);
+
+			 %then save each of the big fields in their respective files
+			 for i=1:length(modelfields),
+				 field=char(modelfields(i));
+				 eval([field '=md.' field ';']);
+
+				 %anything NaN, will be saved 
+
+				 if isnumeric(eval(field)),
+					 if isnan(eval(field)),
+						 if exist(['.' md.repository '.' field],'file');
+							 error(['Repository ' md.repository ' already exists! Delete all files starting with .' md.repository ' in your directory, before activating a new repository']);
+						 end
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 end
+				 end
+			 end
+		 end % }}}
+		 function disp(md) % {{{1
+			 disp(sprintf('\n%s = \n',inputname(1)));
+
+			 %name
+			 disp(sprintf('   Name:'));
+			 if isempty(md.name),
+				 disp(sprintf('      N/A'));
+			 else
+				 disp(sprintf('      %s',md.name));
+			 end
+
+			 %notes: md.notes is a cell array of strings
+			 disp(sprintf('   Notes:'));
+			 if isempty(md.notes),
+				 disp(sprintf('      N/A'));
+			 else
+				 if iscell(md.notes),
+					 strings=md.notes;
+					 for i=1:length(strings),
+						 string=char(strings(i));
+						 if i==1,
+							 disp(sprintf('      %s',string));
+						 else
+							 disp(sprintf('      %s',string));
+						 end
+					 end
+				 elseif ischar(md.notes),
+					 disp(sprintf('      %s',md.notes));
+				 else
+					 disp('display error message: ''notes'' for model is neither a cell array of string, nor a string. Cannot display.');
+				 end
+			 end
+
+			 %repository
+			 if ~isempty(md.repository),
+				 disp(sprintf('   Repository: %s',md.repository));
+			 end
+
+			 %diverse
+			 disp(sprintf('   Mesh: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Mesh: type ''',inputname(1),'.mesh'' to display'));
+			 disp(sprintf('   Parameters: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Boundary conditions: type ''',inputname(1),'.bc'' to display'));
+			 disp(sprintf('%s%s%s','      Observations: type ''',inputname(1),'.obs'' to display'));
+			 disp(sprintf('%s%s%s','      Materials: type ''',inputname(1),'.mat'' to display'));
+			 disp(sprintf('%s%s%s','      Parameters: type ''',inputname(1),'.par'' to display'));
+			 disp(sprintf('%s%s%s','      ExpPar: type ''',inputname(1),'.exppar'' to display'));
+			 disp(sprintf('   Solution parameters: %s',md.repository));
+			 disp(sprintf('%s%s%s','      Qmu: type ''',inputname(1),'.qmu'' to display'));
+			 disp(sprintf('%s%s%s','      FLAIM: type ''',inputname(1),'.flaim'' to display'));
+			 disp(sprintf('%s%s%s','      Diagnostic solution parameters: type ''',inputname(1),'.diagnostic'' to display'));
+			 disp(sprintf('%s%s%s','      Thermal solution parameters: type ''',inputname(1),'.thermal'' to display'));
+			 disp(sprintf('%s%s%s','      Prognostic solution parameters: type ''',inputname(1),'.prognostic'' to display'));
+			 disp(sprintf('%s%s%s','      Transient solution parameters: type ''',inputname(1),'.transient'' to display'));
+			 disp(sprintf('%s%s%s','      Control solution parameters: type ''',inputname(1),'.control'' to display'));
+			 disp(sprintf('   Parallel:'));
+			 disp(sprintf('%s%s%s','      Parallel options: type ''',inputname(1),'.parallel'' to display'));
+			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
+			 disp(sprintf('   Results:'));
+			 disp(sprintf('%s%s%s','      Results: type ''',inputname(1),'.res'' to display'));
+			 if(md.numrifts)disp(sprintf('      rifts: %i',md.numrifts));end
+		 end
+		 %}}}
+		 function mdstruct=struct(md) % {{{1
+			 mdstruct=struct();
+			 modelprops=properties(md);
+
+			 for i=1:length(modelprops),
+				 mdstruct.(modelprops{i})=md.(modelprops{i});
+			 end
+		 end
+		 %}}}
+		 function md=setdefaultparameters(md) % {{{1
+		 %SETDEFAULTPARAMETERS - plug default parameters onto model
+		 %
+		 %   Although the model parameterization should be done in
+		 %   the parameter file, some parameters are initialized here
+		 %   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;
+
+			 %firn density (kg/m^3)
+			 md.rho_firn=830;
+
+			 %ice heat capacity cp (J/kg/K)
+			 md.heatcapacity=2093;
+
+			 %ice latent heat of fusion L (J/kg)
+			 md.latentheat=3.34*10^5;
+
+			 %basal melting rate correction: 
+			 md.gl_melting_rate=0; 
+			 md.basal_melting_rate_correction_apply=0; %do not apply
+
+			 %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
+
+			 %mechanical residue convergence criterion norm(K(uold)uold - F)/norm(F)
+			 md.eps_res=10^-4; %from test233 in nightly runs
+
+			 %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;
+
+			 %maximum of non-linear iterations.
+			 md.max_nonlinear_iterations=100;
+
+			 %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=25;
+
+			 %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;
+
+			 %parameter used to print temporary results (convergence criterion,
+			 %current step,...)
+			 md.verbose=verbose('solution',true,'qmu',true,'control',true);
+
+			 %Stokes
+			 md.stokesreconditioning=10^13;
+			 md.shelf_dampening=0;
+
+			 %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;
+			 md.min_mechanical_constraints=0;
+
+			 %Transient parameters
+
+			 %time between 2 time steps
+			 md.dt=1/2;
+
+			 %number of time steps = md.ndt / md.dt
+			 md.ndt=10*md.dt;
+
+			 %time adaptation? 
+			 md.time_adapt=0;
+			 md.cfl_coefficient=1/2;
+
+			 %Hydrostatic adjustment
+			 md.hydrostatic_adjustment=AbsoluteEnum;
+			 %the artificial diffusivity is used in prognostic to avoid
+			 %numerical wiggles of the solution.
+			 md.artificial_diffusivity=1;
+
+			 %Discontinuous Galerkin is used in prognostic to avoid
+			 %numerical wiggles of the solution and conserve mass
+			 md.prognostic_DG=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;
+
+			 %Solution activated for transient runs. By default we do a
+			 %full analysis: Diagnostic, Prognostic and Thermal
+			 md.isprognostic=1;
+			 md.isdiagnostic=1;
+			 md.isthermal=1;
+
+			 %Control
+
+			 %parameter to be inferred by control methods (only
+			 %drag and B are supported yet)
+			 md.control_type=DragCoefficientEnum();
+
+			 %number of steps in the control methods
+			 md.nsteps=20;
+
+			 %maximum number of iteration in the optimization algorithm for
+			 %each step
+			 md.maxiter=20*ones(md.nsteps,1);
+
+			 %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);
+
+			 %several responses can be used:
+			 md.cm_responses=SurfaceAbsVelMisfitEnum*ones(md.nsteps,1);
+
+			 %cm_jump is used to speed up control method. When
+			 %misfit(1)/misfit(0) < md.cm_jump, we go directly to
+			 %the next step
+			 md.cm_jump=.7*ones(md.nsteps,1); %30 per cent decrement.
+
+			 %stop control solution at the gradient computation and return it? 
+			 md.cm_gradient=0;
+
+			 %eps_cm is a criteria to stop the control methods.
+			 %if J[n]-J[n-1]/J[n] < criteria, the control run stops
+			 %NaN if not applied
+			 md.eps_cm=NaN; %not activated
+
+			 %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);
+
+			 %grounding line migration: 
+			 md.gl_migration=NoneEnum;
+
+			 %How often to save results, default is 1 so save every step
+			 md.output_frequency=1;
+
+			 %Parallelisation parameters
+
+			 %cluster set as none for serial
+			 md.cluster=none;
+
+			 %this option can be activated to load automatically the results
+			 %onto the model after a parallel run by waiting for the lock file
+			 %N minutes that is generated once the solution has converged
+			 %0 to desactivate
+			 md.waitonlock=Inf;
+
+			 %prerun allocation cleanup
+			 md.alloc_cleanup=1;
+
+			 %set petsc options for different analysis
+			 md.petscoptions=petscoptions;
+			 md.petscoptions=addoptions(md.petscoptions,DiagnosticVertAnalysisEnum,mumpsoptions);
+
+			 %solution speed-up
+			 md.kff=1;
+
+			 %hydrology:  from Johnson's 2002 thesis, section 3.5.4
+			 md.hydro_n=.02;
+			 md.hydro_CR=2;
+			 md.hydro_p=2;
+			 md.hydro_q=1;
+			 md.hydro_gamma=1/(md.hydro_n*md.hydro_CR^md.hydro_p*(md.rho_water*md.g)^md.hydro_q);
+			 md.hydro_kn=0;
+
+			 %Rheology law: what is the temperature dependence of B with T
+			 %available: None, Paterson and Arrhenius
+			 md.rheology_law=PatersonEnum;
+
+			 %i/o:
+			 md.io_gather=1;
+		 end
+		 %}}}
+		 function md=subsasgn(md,index,val) % {{{1
+
+			 if ~isempty(md.repository),
+				 field=index(1).subs;
+				 if exist(['.' md.repository '.' field],'file'),
+					 if length(index)==1,
+						 %save val in corresponding repository
+						 eval([field '=val;']);
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 else
+						 %load field from file
+						 load(['.' md.repository '.' field],'-mat');
+						 fieldval=eval(field);
+						 eval([field '=builtin(''subsasgn'',fieldval,index(2),val);']);
+						 eval(['save .' md.repository '.' field ' ' field]);
+					 end
+				 else
+					 md=builtin('subsasgn',md,index,val);
+				 end
+			 else
+				 md=builtin('subsasgn',md,index,val);
+			 end
+	 end %}}}
+		 function result=subsref(md,index) % {{{1
+
+			 if length(index)==1,
+				 index1=index(1);
+				 if(strcmp(index1.subs,'mesh')), displaymesh(md);return; end
+				 if(strcmp(index1.subs,'bc')), displaybc(md);return; end
+				 if(strcmp(index1.subs,'mat')), displaymaterials(md);return; end
+				 if(strcmp(index1.subs,'par')), displayparameters(md);return; end
+				 if(strcmp(index1.subs,'exppar')), displayexppar(md);return; end
+				 if(strcmp(index1.subs,'res')), displayresults(md);return; end
+				 if(strcmp(index1.subs,'obs')), displayobservations(md);return; end
+				 if(strcmp(index1.subs,'qmu')), displayqmu(md);return; end
+				 if(strcmp(index1.subs,'flaim')), displayflaim(md);return; end
+				 if(strcmp(index1.subs,'diagnostic')), displaydiagnostic(md);return; end
+				 if(strcmp(index1.subs,'prognostic')), displayprognostic(md);return; end
+				 if(strcmp(index1.subs,'thermal')), displaythermal(md);return; end
+				 if(strcmp(index1.subs,'transient')), displaytransient(md);return; end
+				 if(strcmp(index1.subs,'control')), displaycontrol(md);return; end
+				 if(strcmp(index1.subs,'mesh2grid')), displaymesh2grid(md);return; end
+				 if(strcmp(index1.subs,'parallel')), displayparallel(md);return; end
+
+				 if ~isempty(md.repository),
+					 %load variable needed, if necessary
+					 if exist(['.' md.repository '.' index1.subs],'file'),
+						 structure=load(['.' md.repository '.' index1.subs],'-mat');
+						 eval(['result=structure.' index1.subs ';']);
+					 else
+						 %the file does not exist, just return built in value
+						 result=builtin('subsref',md,index);
+					 end
+				 else
+					 %return built in value
+					 result=builtin('subsref',md,index);
+				 end
+			 else
+				 if ~isempty(md.repository),
+					 if exist(['.' md.repository '.' index(1).subs],'file'),
+						 %load variable needed
+						 variable=load(['.' md.repository '.' index(1).subs],'-mat');
+
+						 %use temporary structure.
+						 result=builtin('subsref',variable,index);
+					 else
+						 result=builtin('subsref',md,index);
+					 end
+				 else
+					 result=builtin('subsref',md,index);
+				 end
+			 end
+		 end %}}}
+	 end
+ end
Index: /issm/trunk/startup.m
===================================================================
--- /issm/trunk/startup.m	(revision 8457)
+++ /issm/trunk/startup.m	(revision 8458)
@@ -61,4 +61,13 @@
 	addpath(genpath_ice([ISSM_TIER '/src/m/usr/' USERNAME]));
 end
+
+%version dependent stuff
+rmpath(genpath_ice([ISSM_TIER '/src/m/classes/version']));
+if verLessThan('matlab','7.6.0'),
+	addpath(genpath_ice([ISSM_TIER '/src/m/classes/version/7.6']));
+else
+	addpath(genpath_ice([ISSM_TIER '/src/m/classes/version/7.7']));
+end
+
 clear USERNAME ISSM_TIER;