Index: /issm/trunk/src/m/classes/model.m
===================================================================
--- /issm/trunk/src/m/classes/model.m	(revision 7117)
+++ /issm/trunk/src/m/classes/model.m	(revision 7117)
@@ -0,0 +1,771 @@
+%MODEL class definition
+%
+%   Usage:
+%      md = model(varargin)
+
+classdef model
+    properties (SetAccess=public) 
+		 % {{{1
+
+		 %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=NaN;
+		 dim=0;
+		 numberofelements=0;
+		 numberofgrids=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;
+		 numberofgrids2d=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
+		 gridonhutter=NaN;
+		 gridonmacayeal=NaN;
+		 gridonpattyn=NaN;
+		 gridonstokes=NaN;
+		 borderstokes=NaN;
+
+		 %Stokes mesh
+		 stokesreconditioning=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
+		 uppergrids=NaN;
+		 upperelements=NaN;
+		 lowerelements=NaN;
+		 lowergrids=NaN;
+
+		 %Extrusion
+		 numlayers=0;
+		 extrusionexponent=0;
+		 elementonbed=NaN;
+		 elementonsurface=NaN;
+		 gridonbed=NaN;
+		 gridonsurface=NaN;
+		 minh=0;
+		 firn_layer=NaN;
+
+		 %Extraction
+		 extractedgrids=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;
+
+		 %Geometrical parameters
+		 elementoniceshelf=NaN;
+		 elementonicesheet=NaN;
+		 elementonwater=NaN;
+		 gridoniceshelf=NaN;
+		 gridonicesheet=NaN;
+		 gridonwater=NaN;
+		 surface=NaN;
+		 thickness=NaN;
+		 thickness_coeff=NaN;
+		 bed=NaN;
+		 bathymetry=NaN;
+
+		 %Boundary conditions
+		 gridonboundary=NaN;
+		 pressureload=NaN;
+		 spcvelocity=NaN;
+		 spctemperature=NaN;
+		 spcthickness=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;
+		 accumulation_rate=NaN;
+		 dhdt=NaN;
+		 geothermalflux=NaN;
+		 observed_temperature=NaN;
+		 thickness_obs=NaN;
+
+		 %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;
+
+		 %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=0;
+		 cm_max=0;
+		 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
+		 melting_rate=NaN;
+		 pressure=NaN;
+
+		 %Parallelisation parameters
+		 cluster=none;
+		 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;
+		 petscoptions_analyses=[];
+		 petscoptions_strings={};
+
+		 %Analysis and sub_analysis
+		 analysis_type='';
+		 sub_analysis_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;
+
+		 %grounding line migration: 
+		 gl_migration=NaN;
+
+		 %solution parameters
+		 kff=0;
+
+		 %partitioner:
+		 adjacency=[];
+		 vwgt=[];
+		 %}}}
+	 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','      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 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;
+
+			 %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.3;
+
+			 %parameter used to print temporary results (convergence criterion,
+			 %current step,...)
+			 md.verbose=verbose('solution',true,'qmu',true,'control',true);
+
+			 %Stokes mesh
+			 md.stokesreconditioning=10^13;
+
+			 %Penalties
+
+			 %the penalty offset is the power used to compute the value
+			 %of the penalty as follows
+			 %kappa=max(stiffness matrix)*10^penalty_offset
+			 md.penalty_offset=3;
+
+			 %penalty_melting is the penalty used to keep T<Tpmp
+			 %kappa=penalty_melting
+			 md.penalty_melting=10^7;
+
+			 %in some solutions, it might be needed to stop a run when only
+			 %a few constraints remain unstable. For thermal computation, this
+			 %parameter is often used.
+			 md.penalty_lock=10;
+
+			 %minimum thickness to avoid stiffness singularity, used in
+			 %the parameter file
+			 md.minh=1;
+
+			 %in some cases, it might be needed to stop a run when only
+			 %a few constraints remain unstable. For thermal computation, this
+			 %parameter is often used.
+			 md.min_thermal_constraints=0;
+			 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;
+
+			 %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.
+
+			 %cm_min is the minimum acceptable value of the
+			 %inversed parameter (B>0 for example)
+			 md.cm_min=0;
+
+			 %cm_max is the maximum acceptable value of the
+			 %inversed parameter (drag<200 for example)
+			 md.cm_max=200;
+
+			 %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=0;
+
+			 %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;
+
+			 %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,'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