Changeset 9597

Timestamp:

09/02/11 17:10:47 (14 years ago)

Author:

Mathieu Morlighem

Message:

Put yts g and referencetemperature in md.constants

Location:

issm/trunk

Files:

• src/c/EnumDefinitions/EnumDefinitions.h (modified) (3 diffs)
• src/c/modules/EnumToStringx/EnumToStringx.cpp (modified) (3 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (modified) (1 diff)
• src/c/modules/StringToEnumx/StringToEnumx.cpp (modified) (3 diffs)
• src/c/objects/Elements/Penta.cpp (modified) (4 diffs)
• src/c/objects/Elements/Tria.cpp (modified) (4 diffs)
• src/c/objects/Materials/Matpar.cpp (modified) (2 diffs)
• src/c/solutions/convergence.cpp (modified) (1 diff)
• src/c/solutions/transient_core.cpp (modified) (1 diff)
• src/m/classes/model/model.m (modified) (8 diffs)
• src/m/model/basalstress.m (modified) (1 diff)
• src/m/model/drivingstress.m (modified) (1 diff)
• src/m/model/effectivepressure.m (modified) (1 diff)
• src/m/model/extrude.m (modified) (1 diff)
• src/m/model/ismodelselfconsistent.m (modified) (3 diffs)
• src/m/model/marshall.m (modified) (3 diffs)
• src/m/model/mechanicalproperties.m (modified) (1 diff)
• src/m/model/misfit.m (modified) (2 diffs)
• src/m/model/plot/plot_basaldrag.m (modified) (1 diff)
• src/m/model/plot/plot_gridded.m (modified) (2 diffs)
• src/m/model/plot/plot_transient_field.m (modified) (1 diff)
• src/m/model/plot/plot_transient_movie.m (modified) (1 diff)
• src/m/model/sia.m (modified) (1 diff)
• src/m/solutions/convergence.m (modified) (1 diff)
• src/m/solutions/transient_core.m (modified) (1 diff)
• src/m/utils/Mesh/BuildAdaptedMesh.m (modified) (3 diffs)
• test/Miscellaneous/connectivity/Square.par (modified) (1 diff)
• test/NightlyRun/test1301.m (modified) (1 diff)
• test/NightlyRun/test1302.m (modified) (1 diff)
• test/NightlyRun/test1501.m (modified) (2 diffs)
• test/NightlyRun/test1502.m (modified) (2 diffs)
• test/NightlyRun/test233.m (modified) (1 diff)
• test/NightlyRun/test234.m (modified) (1 diff)
• test/NightlyRun/test621.m (modified) (1 diff)
• test/NightlyRun/test622.m (modified) (1 diff)
• test/NightlyRun/test623.m (modified) (1 diff)
• test/NightlyRun/test624.m (modified) (1 diff)
• test/NightlyRun/test625.m (modified) (1 diff)
• test/NightlyRun/test626.m (modified) (1 diff)
• test/Par/ISMIPC.par (modified) (1 diff)
• test/Par/ISMIPD.par (modified) (1 diff)
• test/Par/ISMIPF.par (modified) (1 diff)
• test/Par/RoundSheetEISMINT.par (modified) (1 diff)
• test/Par/SquareThermal.par (modified) (1 diff)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -363 +363 @@
         SsetEnum,
         GroundinglineMigrationEnum,
-        YtsEnum,
+        ConstantsYtsEnum,
         /*}}}*/
         /*Interpolation {{{1*/
@@ -412 +412 @@
         VerticesTypeEnum,
         HydroKnEnum,
-        GEnum,
+        ConstantsGEnum,
         SpcvzEnum,
         NumberOfNodes2DEnum,
@@ -446 +446 @@
         ThermalconductivityEnum,
         MeltingpointEnum,
-        ReferencetemperatureEnum,
+        ConstantsReferencetemperatureEnum,
         DhdtEnum,
         GeothermalfluxEnum,

issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -314 +314 @@
                 case SsetEnum : return "Sset";
                 case GroundinglineMigrationEnum : return "GroundinglineMigration";
-                case YtsEnum : return "Yts";
+                case ConstantsYtsEnum : return "ConstantsYts";
                 case TriangleInterpEnum : return "TriangleInterp";
                 case BilinearInterpEnum : return "BilinearInterp";
@@ -355 +355 @@
                 case VerticesTypeEnum : return "VerticesType";
                 case HydroKnEnum : return "HydroKn";
-                case GEnum : return "G";
+                case ConstantsGEnum : return "ConstantsG";
                 case SpcvzEnum : return "Spcvz";
                 case NumberOfNodes2DEnum : return "NumberOfNodes2D";
@@ -389 +389 @@
                 case ThermalconductivityEnum : return "Thermalconductivity";
                 case MeltingpointEnum : return "Meltingpoint";
-                case ReferencetemperatureEnum : return "Referencetemperature";
+                case ConstantsReferencetemperatureEnum : return "ConstantsReferencetemperature";
                 case DhdtEnum : return "Dhdt";
                 case GeothermalfluxEnum : return "Geothermalflux";

issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -35 +35 @@
         parameters->AddObject(iomodel->CopyConstantObject(MaxNonlinearIterationsEnum));
         parameters->AddObject(iomodel->CopyConstantObject(MaxSteadystateIterationsEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(YtsEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ConstantsYtsEnum));
         parameters->AddObject(iomodel->CopyConstantObject(DtEnum));
         parameters->AddObject(iomodel->CopyConstantObject(NdtEnum));
@@ -46 +46 @@
         parameters->AddObject(iomodel->CopyConstantObject(BetaEnum));
         parameters->AddObject(iomodel->CopyConstantObject(MeltingpointEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(ReferencetemperatureEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
         parameters->AddObject(iomodel->CopyConstantObject(LatentheatEnum));
         parameters->AddObject(iomodel->CopyConstantObject(HeatcapacityEnum));

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

@@ -41 +41 @@
         /*Fetch parameters: */
         iomodel->Constant(&dim,DimEnum);
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
-        iomodel->Constant(&g,GEnum);
+        iomodel->Constant(&g,ConstantsGEnum);
         iomodel->Constant(&rho_ice,RhoIceEnum);
         iomodel->Constant(&stokesreconditioning,StokesreconditioningEnum);

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp

@@ -28 +28 @@
 
         /*Fetch parameters: */
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
         iomodel->Constant(&ishutter,IshutterEnum);

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp

@@ -25 +25 @@
         /*Fetch parameters: */
         iomodel->Constant(&dim,DimEnum);
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
 

issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp

@@ -29 +29 @@
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
         iomodel->Constant(&heatcapacity,HeatcapacityEnum);
-        iomodel->Constant(&referencetemperature,ReferencetemperatureEnum);
+        iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
 
         /*Recover pointer: */

issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -312 +312 @@
         else if (strcmp(name,"Sset")==0) return SsetEnum;
         else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
-        else if (strcmp(name,"Yts")==0) return YtsEnum;
+        else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
         else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
         else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
@@ -353 +353 @@
         else if (strcmp(name,"VerticesType")==0) return VerticesTypeEnum;
         else if (strcmp(name,"HydroKn")==0) return HydroKnEnum;
-        else if (strcmp(name,"G")==0) return GEnum;
+        else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
         else if (strcmp(name,"Spcvz")==0) return SpcvzEnum;
         else if (strcmp(name,"NumberOfNodes2D")==0) return NumberOfNodes2DEnum;
@@ -387 +387 @@
         else if (strcmp(name,"Thermalconductivity")==0) return ThermalconductivityEnum;
         else if (strcmp(name,"Meltingpoint")==0) return MeltingpointEnum;
-        else if (strcmp(name,"Referencetemperature")==0) return ReferencetemperatureEnum;
+        else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
         else if (strcmp(name,"Dhdt")==0) return DhdtEnum;
         else if (strcmp(name,"Geothermalflux")==0) return GeothermalfluxEnum;

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -4846 +4846 @@
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
         iomodel->Constant(&numberofelements,NumberOfElementsEnum);
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
 
         /*Branch on type of vector: nodal or elementary: */
@@ -5153 +5153 @@
 
         /*Fetch parameters: */
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&control_analysis,ControlAnalysisEnum);
         iomodel->Constant(&num_control_type,NumControlTypeEnum);
@@ -7503 +7503 @@
 
         /*Fetch parameters: */
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&prognostic_DG,PrognosticDGEnum);
         iomodel->Constant(&dakota_analysis,DakotaAnalysisEnum);
@@ -7509 +7509 @@
         iomodel->Constant(&beta,BetaEnum);
         iomodel->Constant(&heatcapacity,HeatcapacityEnum);
-        iomodel->Constant(&referencetemperature,ReferencetemperatureEnum);
+        iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
         iomodel->Constant(&meltingpoint,MeltingpointEnum);
         iomodel->Constant(&latentheat,LatentheatEnum);

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -3273 +3273 @@
         iomodel->Constant(&control_analysis,ControlAnalysisEnum);
         iomodel->Constant(&num_control_type,NumControlTypeEnum);
-        iomodel->Constant(&yts,YtsEnum); 
+        iomodel->Constant(&yts,ConstantsYtsEnum); 
         iomodel->Constant(&num_cm_responses,NumCmResponsesEnum);
 
@@ -3900 +3900 @@
         iomodel->Constant(&numberofvertices,NumberOfVerticesEnum);
         iomodel->Constant(&numberofelements,NumberOfElementsEnum);
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
 
         /*Branch on type of vector: nodal or elementary: */
@@ -4581 +4581 @@
 
         /*recover parameters: */
-        parameters->FindParam(&yts,YtsEnum);
+        parameters->FindParam(&yts,ConstantsYtsEnum);
         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
 
@@ -5218 +5218 @@
 
         /*Fetch parameters: */
-        iomodel->Constant(&yts,YtsEnum);
+        iomodel->Constant(&yts,ConstantsYtsEnum);
         iomodel->Constant(&prognostic_DG,PrognosticDGEnum);
         iomodel->Constant(&dakota_analysis,DakotaAnalysisEnum);

issm/trunk/src/c/objects/Materials/Matpar.cpp

@@ -33 +33 @@
         iomodel->Constant(&this->beta,BetaEnum);
         iomodel->Constant(&this->meltingpoint,MeltingpointEnum);
-        iomodel->Constant(&this->referencetemperature,ReferencetemperatureEnum);
+        iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
         iomodel->Constant(&this->mixed_layer_capacity,MixedLayerCapacityEnum);
         iomodel->Constant(&this->thermal_exchange_velocity,ThermalExchangeVelocityEnum);
-        iomodel->Constant(&this->g,GEnum);
+        iomodel->Constant(&this->g,ConstantsGEnum);
         
         iomodel->Constant(&this->kn,HydroKnEnum);
@@ -255 +255 @@
                         break;
 
-                case  ReferencetemperatureEnum:
+                case  ConstantsReferencetemperatureEnum:
                         this->referencetemperature=constant;
                         break;

issm/trunk/src/c/solutions/convergence.cpp

@@ -41 +41 @@
         parameters->FindParam(&eps_rel,EpsRelEnum);
         parameters->FindParam(&eps_abs,EpsAbsEnum);
-        parameters->FindParam(&yts,YtsEnum);
+        parameters->FindParam(&yts,ConstantsYtsEnum);
 
         /*Display solver caracteristics*/

issm/trunk/src/c/solutions/transient_core.cpp

@@ -31 +31 @@
         femmodel->parameters->FindParam(&finaltime,NdtEnum);
         femmodel->parameters->FindParam(&dt,DtEnum);
-        femmodel->parameters->FindParam(&yts,YtsEnum);
+        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
         femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);

issm/trunk/src/m/classes/model/model.m

@@ -8 +8 @@
                  % {{{1
                  %Careful here: no other class should be used as default value this is a bug of matlab
-                 cluster = modelfield('default',0,'marshall',false);
-                 verbose = modelfield('default',0,'marshall',true,'preprocess','marshallverbose','format','Integer');
-                 results = modelfield('default',0,'marshall',false);
-                 solver  = modelfield('default',0,'marshall',false);
-                 debug   = modelfield('default',0,'marshall',false);
+                 cluster   = modelfield('default',0,'marshall',false);
+                 verbose   = modelfield('default',0,'marshall',true,'preprocess','marshallverbose','format','Integer');
+                 results   = modelfield('default',0,'marshall',false);
+                 solver    = modelfield('default',0,'marshall',false);
+                 debug     = modelfield('default',0,'marshall',false);
+                 constants = modelfield('default',0,'marshall',true);
 
                  %FIXME: all other fields should belong to other classes
@@ -108 +109 @@
                  thermalconductivity        = modelfield('default',0,'marshall',true,'format','Double');
                  meltingpoint               = modelfield('default',0,'marshall',true,'format','Double');
-                 referencetemperature       = modelfield('default',0,'marshall',true,'format','Double'); %for enthalpy
                  beta                       = modelfield('default',0,'marshall',true,'format','Double');
                  mixed_layer_capacity       = modelfield('default',0,'marshall',true,'format','Double');
@@ -117 +117 @@
 
                  %Physical parameters
-                 g                = modelfield('default',0,'marshall',true,'format','Double');
-                 yts              = modelfield('default',0,'marshall',true,'format','Double');
+
                  drag_coefficient = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
                  drag_p           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',2);
@@ -336 +335 @@
                                                  md = in;
                                          elseif (isa(in,'numeric') & in==0)
-                                                         %just requesting model, without default parameters, do nothing.
+                                                 %requesting template
+                                                 %FIXME: for now recursive call
+                                                 md2=model;
+                                                 modelprops=properties('model');
+                                                 for i=1:length(modelprops);
+                                                         fieldname=modelprops{i};
+                                                         if md.(fieldname).marshall & isobject(md2.(fieldname)) & isempty(md.(fieldname).preprocess),
+                                                                 eval(['md.' fieldname '=' fieldname '(0);']);
+                                                         end
+                                                 end
                                          else
                                                  error('model constructor error message: unknown constructor for ''model'' class object');
@@ -438 +446 @@
                          if isfield(structmd,'gridonboundary'), md.nodeonboundary=structmd.gridonboundary; end
                          if isfield(structmd,'petscoptions'), md.solver=structmd.petscoptions; end
+                         if isfield(structmd,'g'), md.constants.g=structmd.g; end
+                         if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end
                          if isfield(structmd,'dakota_analysis'), md.dakota_analysis=structmd.dakota_analysis; end
 
@@ -543 +553 @@
                          md.results=struct();
                          md.debug=debug;
+                         md.constants=constants;
 
 
@@ -568 +579 @@
                          md.meltingpoint=273.15;
 
-                         %the reference temperature for enthalpy model (cf Aschwanden)
-                         md.referencetemperature=223.15;
-
                          %rate of change of melting point with pressure (K/Pa)
                          md.beta=9.8*10^-8;
@@ -581 +589 @@
 
                          %Physical parameters
-
-                         %acceleration due to gravity (m/s^2)
-                         md.g=9.81;
-
-                         %converstion from year to seconds
-                         md.yts=365*24*3600;
 
                          %Solver parameters

issm/trunk/src/m/model/basalstress.m

@@ -13 +13 @@
 
 %compute horizontal velocity
-ub=sqrt(md.vx.^2+md.vy.^2)/md.yts;
-ubx=md.vx/md.yts;
-uby=md.vy/md.yts;
+ub=sqrt(md.vx.^2+md.vy.^2)/md.constants.yts;
+ubx=md.vx/md.constants.yts;
+uby=md.vy/md.constants.yts;
 
 %compute basal drag
-bx=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*ubx.^s;
-by=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*uby.^s;
-b=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*ub.^s;
+bx=(md.constants.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*ubx.^s;
+by=(md.constants.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*uby.^s;
+b=(md.constants.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag_coefficient).^2.*ub.^s;

issm/trunk/src/m/model/drivingstress.m

@@ -14 +14 @@
 thickness_bar=(md.thickness(md.elements(:,1))+md.thickness(md.elements(:,2))+md.thickness(md.elements(:,3)))/3;
 
-px=md.rho_ice*md.g*thickness_bar.*sx;
-py=md.rho_ice*md.g*thickness_bar.*sy;
+px=md.rho_ice*md.constants.g*thickness_bar.*sx;
+py=md.rho_ice*md.constants.g*thickness_bar.*sy;
 pmag=sqrt(px.^2+py.^2);

issm/trunk/src/m/model/effectivepressure.m

@@ -10 +10 @@
 
 
-Neff=md.rho_ice*md.g*md.thickness+md.rho_ice*md.g*md.bed;
+Neff=md.rho_ice*md.constants.g*md.thickness+md.rho_ice*md.constants.g*md.bed;
 
 

issm/trunk/src/m/model/extrude.m

@@ -217 +217 @@
 %Put lithostatic pressure is there is an existing pressure
 if ~isnan(md.pressure),
-        md.pressure=md.g*md.rho_ice*(md.surface-md.z);
+        md.pressure=md.constants.g*md.rho_ice*(md.surface-md.z);
 end
 

issm/trunk/src/m/model/ismodelselfconsistent.m

@@ -83 +83 @@
 %NO NAN {{{1
 fields={'numberofelements','numberofnodes','x','y','z','drag_coefficient','drag_p','drag_q',...
-        'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','nsteps','maxiter',...
+        'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','constants.g','lowmem','nsteps','maxiter',...
         'eps_res','max_nonlinear_iterations','rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface','elementconnectivity'};
 checknan(md,fields);
@@ -89 +89 @@
 %FIELDS >= 0 {{{1
 fields={'numberofelements','numberofnodes','elements','drag_coefficient','drag_p','drag_q',...
-        'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter',...
+        'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','constants.g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter',...
         'lowmem','rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface'};
 checkgreater(md,fields,0);
@@ -95 +95 @@
 %FIELDS > 0 {{{1
 fields={'numberofelements','numberofnodes','elements','drag_p',...
-        'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter'};
+        'rho_ice','rho_water','rheology_B','thickness','constants.g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter'};
 checkgreaterstrict(md,fields,0);
 %}}}

issm/trunk/src/m/model/marshall.m

@@ -18 +18 @@
 %automatically marshall model md, using template information available from an empty model class.
 template=model.template();
-fields=fieldnames(template);
-for i=1:length(fields),
-        fieldname=fields{i};
-        if template.(fieldname).marshall,
-                if ~isempty(template.(fieldname).preprocess),
-                        eval(['fieldvalue=' template.(fieldname).preprocess '(md.(fieldname));']);
-                else
-                        fieldvalue=md.(fieldname);
-                end
-                WriteData(fid,fieldname,fieldvalue,template.(fieldname));
-        end
-end
+MarshallObject(fid,template,md,'');
 
 %FIXME add some fields needed by parameters
@@ -40 +29 @@
 end
 
+function MarshallObject(fid,template,object,prefix) % {{{
+fields=fieldnames(template);
+for i=1:length(fields),
+        fieldname=fields{i};
+        if ~isa(template.(fieldname),'modelfield'),
+                %Recursive call
+                MarshallObject(fid,template.(fieldname),object.(fieldname),[prefix fieldname])
+        elseif template.(fieldname).marshall,
+                if ~isempty(template.(fieldname).preprocess),
+                        eval(['fieldvalue=' template.(fieldname).preprocess '(object.(fieldname));']);
+                else
+                        fieldvalue=object.(fieldname);
+                end
+
+                %Capitalize for enums
+                fieldnamecap=fieldname;
+                fieldnamecap(1)=upper(fieldname(1));
+                prefixcap=prefix;
+                if ~isempty(prefix), prefixcap(1)=upper(prefix(1)); end
+
+                %Write Data
+                WriteData(fid,[prefixcap fieldnamecap],fieldvalue,template.(fieldname));
+        end
+end
+        % }}}
 function WriteData(fid,fieldname,fieldvalue,fieldprop) % {{{
 %WRITEDATA - write model field in binary file
@@ -254 +268 @@
 
 %FIXME: delete this one
-if strcmp(string,'numberofnodes'),
-        string='numberofvertices';
+if strcmp(string,'Numberofnodes'),
+        string='Numberofvertices';
 end
 

issm/trunk/src/m/model/mechanicalproperties.m

@@ -37 +37 @@
 
 %compute shear
-vxlist=vx(index)/md.yts;
-vylist=vy(index)/md.yts;
+vxlist=vx(index)/md.constants.yts;
+vylist=vy(index)/md.constants.yts;
 ux=(vxlist.*alpha)*summation;
 uy=(vxlist.*beta)*summation;

issm/trunk/src/m/model/misfit.m

@@ -8 +8 @@
 %      J=misfit(md)
 %
-
 
 if md.dim==2,
@@ -32 +31 @@
 
 %compute delta v on elements
-deltav=1/2*(   (vx-vx_obs).^2+(vy-vy_obs).^2)/md.yts^2;
+deltav=1/2*(   (vx-vx_obs).^2+(vy-vy_obs).^2)/md.constants.yts^2;
 deltav_elem=deltav(elements)*[1;1;1]/3;
 

issm/trunk/src/m/model/plot/plot_basaldrag.m

@@ -22 +22 @@
 %compute horizontal velocity
 if strcmpi(type,'basal_drag')
-        ub=sqrt(md.vx.^2+md.vy.^2)/md.yts;
+        ub=sqrt(md.vx.^2+md.vy.^2)/md.constants.yts;
 elseif strcmpi(type,'basal_dragx')
-        ub=md.vx/md.yts;
+        ub=md.vx/md.constants.yts;
 elseif strcmpi(type,'basal_dragy')
-        ub=md.vy/md.yts;
+        ub=md.vy/md.constants.yts;
 end
 
 %compute basal drag
-drag=(max(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed),0)).^r.*(md.drag_coefficient).^2.*ub.^s/1000;
+drag=(max(md.constants.g*(md.rho_ice*md.thickness+md.rho_water*md.bed),0)).^r.*(md.drag_coefficient).^2.*ub.^s/1000;
 
 %Figure out if this is a Section plot

issm/trunk/src/m/model/plot/plot_gridded.m

@@ -37 +37 @@
         data_max=max(data_grid(:));
 end
-if whiteup,
-        white  =data_max + (data_max-data_min)/55;
-        options=changefieldvalue(options,'caxis',[data_min white]);
-else
-        white  =data_min - (data_max-data_min)/55;
-        options=changefieldvalue(options,'caxis',[white data_max]);
-end
-data_grid(isnan(data_grid))=white;
+%if whiteup,
+%       white  =data_max + (data_max-data_min)/55;
+%       options=changefieldvalue(options,'caxis',[data_min white]);
+%else
+%       white  =data_min - (data_max-data_min)/55;
+%       options=changefieldvalue(options,'caxis',[white data_max]);
+%end
+%data_grid(isnan(data_grid))=white;
 
 %Select plot area 
@@ -52 +52 @@
 h=imagesc(xlim,ylim,data_grid);set(gca,'YDir','normal');
 map=getfieldvalue(options,'colormap',jet);
-if whiteup,
-        map(end,:)=[1 1 1];
-else
-        map(1,:)=[1 1 1];
-end
-options=changefieldvalue(options,'colormap',map);
+%if whiteup,
+%       map(end,:)=[1 1 1];
+%else
+%       map(1,:)=[1 1 1];
+%end
+%options=changefieldvalue(options,'colormap',map);
 
 %last step: mesh gridded?

issm/trunk/src/m/model/plot/plot_transient_field.m

@@ -34 +34 @@
         %process data and change title if needed
         [data datatype]=processdata(md,transient(steps(i)).(field),options);
-        options=changefieldvalue(options,'title',[field ' at time ' num2str(transient(steps(i)).time/md.yts) ' a']);
+        options=changefieldvalue(options,'title',[field ' at time ' num2str(transient(steps(i)).time/md.constants.yts) ' a']);
 
         %create plot of step i

issm/trunk/src/m/model/plot/plot_transient_movie.m

@@ -24 +24 @@
                 [data datatype]=processdata(md,results(i).(field),options);
 
-                titlestring=[field ' at time ' num2str(results(i).time/md.yts) ' year'];
+                titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
                 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                 apply_options_movie(md,options,titlestring);

issm/trunk/src/m/model/sia.m

@@ -22 +22 @@
 Ael=Bel.^(-3);
 
-velx=-2*(md.rho_ice*md.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
-vely=-2*(md.rho_ice*md.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
+velx=-2*(md.rho_ice*md.constants.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
+vely=-2*(md.rho_ice*md.constants.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
 vel=sqrt(velx.^2+vely.^2);

issm/trunk/src/m/solutions/convergence.m

@@ -5 +5 @@
 
 %Get convergence options
-yts=parameters.Yts;
+yts=parameters.ConstantsYts;
 eps_res=parameters.EpsRes;
 eps_rel=parameters.EpsRel;

issm/trunk/src/m/solutions/transient_core.m

@@ -9 +9 @@
         ndt=femmodel.parameters.Ndt;
         dt=femmodel.parameters.Dt;
-        yts=femmodel.parameters.Yts;
+        yts=femmodel.parameters.ConstantsYts;
         control_analysis=femmodel.parameters.ControlAnalysis;
         solution_type=femmodel.parameters.SolutionType;

issm/trunk/src/m/utils/Mesh/BuildAdaptedMesh.m

@@ -38 +38 @@
         %interpolate velocities onto mesh
         disp('   interpolating velocities');
-        md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
-        md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
+        md.vx_obs=DataInterp(x_m,y_m,vx*md.constants.yts,md.x,md.y);
+        md.vy_obs=DataInterp(x_m,y_m,vy*md.constants.yts,md.x,md.y);
         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 
@@ -49 +49 @@
 disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
 
-md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
-md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
+md.vx_obs=DataInterp(x_m,y_m,vx*md.constants.yts,md.x,md.y);
+md.vy_obs=DataInterp(x_m,y_m,vy*md.constants.yts,md.x,md.y);
 md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 return;
@@ -59 +59 @@
 md=parameterize(md,parameterfile);
 
-%md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
-%md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
+%md.vx_obs=DataInterp(x_m,y_m,vx*md.constants.yts,md.x,md.y);
+%md.vy_obs=DataInterp(x_m,y_m,vy*md.constants.yts,md.x,md.y);
 %md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 %return;

issm/trunk/test/Miscellaneous/connectivity/Square.par

@@ -2 +2 @@
 
 %dynamics
-md.dt=1*md.yts; %1 year
+md.dt=1*md.constants.yts; %1 year
 md.ndt=md.dt*10; 
 md.artificial_diffusivity=1;

issm/trunk/test/NightlyRun/test1301.m

@@ -22 +22 @@
 %melting heat = geothermal flux
 %Mb*L*rho=G   => Mb=G/L*rho
-melting=md.geothermalflux/(md.rho_ice*md.latentheat)*md.yts;
+melting=md.geothermalflux/(md.rho_ice*md.latentheat)*md.constants.yts;
 
 %modeled  results

issm/trunk/test/NightlyRun/test1302.m

@@ -20 +20 @@
 %analytical results
 %d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
-alpha=0.1/md.yts*md.rho_ice*md.heatcapacity/md.thermalconductivity;   %alpha=w rho_ice c /k  and w=0.1m/an
+alpha=0.1/md.constants.yts*md.rho_ice*md.heatcapacity/md.thermalconductivity;   %alpha=w rho_ice c /k  and w=0.1m/an
 A=10/(exp(alpha*(-1000))-1);    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
 B=-A;

issm/trunk/test/NightlyRun/test1501.m

@@ -131 +131 @@
                 [data datatype]=processdata(md,results(i).(field),options);
 
-                titlestring=[field ' at time ' num2str(results(i).time/md.yts) ' year'];
+                titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
                 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                 options=changefieldvalue(options,'title',titlestring);
@@ -146 +146 @@
                 [data datatype]=processdata(md,results(i).(field),options);
 
-                titlestring=[field ' at time ' num2str(results(i).time/md.yts) ' year'];
+                titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
                 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                 options=changefieldvalue(options,'title',titlestring);

issm/trunk/test/NightlyRun/test1502.m

@@ -137 +137 @@
                 [data datatype]=processdata(md,results(i).(field),options);
 
-                titlestring=[field ' at time ' num2str(results(i).time/md.yts) ' year'];
+                titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
                 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                 options=changefieldvalue(options,'title',titlestring);
@@ -152 +152 @@
                 [data datatype]=processdata(md,results(i).(field),options);
 
-                titlestring=[field ' at time ' num2str(results(i).time/md.yts) ' year'];
+                titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
                 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                 options=changefieldvalue(options,'title',titlestring);

issm/trunk/test/NightlyRun/test233.m

@@ -61 +61 @@
 
 %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
-%ey_c=(md.rho_ice*md.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
-%vy_c=ey_c.*md.y*md.yts;
+%ey_c=(md.rho_ice*md.constants.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
+%vy_c=ey_c.*md.y*md.constants.yts;
 
 %Fields and tolerances to track changes

issm/trunk/test/NightlyRun/test234.m

@@ -61 +61 @@
 
 %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
-%ey_c=(md.rho_ice*md.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
-%vy_c=ey_c.*md.y*md.yts;
+%ey_c=(md.rho_ice*md.constants.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
+%vy_c=ey_c.*md.y*md.constants.yts;
 
 %Fields and tolerances to track changes

issm/trunk/test/NightlyRun/test621.m

@@ -14 +14 @@
 md.cm_responses=201*ones(md.nsteps,1);
 md.weights=ones(md.numberofnodes,1);
-md.optscal=10/md.yts*ones(md.nsteps,1);
+md.optscal=10/md.constants.yts*ones(md.nsteps,1);
 md.maxiter=4*ones(md.nsteps,1);
 md.cm_jump=0.99*ones(md.nsteps,1);

issm/trunk/test/NightlyRun/test622.m

@@ -14 +14 @@
 md.cm_responses=201*ones(md.nsteps,1);
 md.weights=ones(md.numberofnodes,1);
-md.optscal=10/md.yts*ones(md.nsteps,1);
+md.optscal=10/md.constants.yts*ones(md.nsteps,1);
 md.maxiter=4*ones(md.nsteps,1);
 md.cm_jump=0.99*ones(md.nsteps,1);

issm/trunk/test/NightlyRun/test623.m

@@ -15 +15 @@
 md.cm_responses=201*ones(md.nsteps,1);
 md.weights=ones(md.numberofnodes,1);
-md.optscal=10/md.yts*ones(md.nsteps,1);
+md.optscal=10/md.constants.yts*ones(md.nsteps,1);
 md.maxiter=4*ones(md.nsteps,1);
 md.cm_jump=0.99*ones(md.nsteps,1);

issm/trunk/test/NightlyRun/test624.m

@@ -15 +15 @@
 md.cm_responses=201*ones(md.nsteps,1);
 md.weights=ones(md.numberofnodes,1);
-md.optscal=10/md.yts*ones(md.nsteps,1);
+md.optscal=10/md.constants.yts*ones(md.nsteps,1);
 md.maxiter=4*ones(md.nsteps,1);
 md.cm_jump=0.99*ones(md.nsteps,1);

issm/trunk/test/NightlyRun/test625.m

@@ -17 +17 @@
 md.control_type={'Vx' 'Vy'};
 md.prognostic_DG=0;
-md.optscal=[10/md.yts*ones(md.nsteps,1) 10/md.yts*ones(md.nsteps,1)];
+md.optscal=[10/md.constants.yts*ones(md.nsteps,1) 10/md.constants.yts*ones(md.nsteps,1)];
 md.cm_min=[-2000*ones(md.numberofnodes,1) -2000*ones(md.numberofnodes,1)];
 md.cm_max=[+2000*ones(md.numberofnodes,1) +2000*ones(md.numberofnodes,1)];

issm/trunk/test/NightlyRun/test626.m

@@ -17 +17 @@
 md.control_type={'Vx' 'Vy'};
 md.prognostic_DG=0;
-md.optscal=[10/md.yts*ones(md.nsteps,1) 10/md.yts*ones(md.nsteps,1)];
+md.optscal=[10/md.constants.yts*ones(md.nsteps,1) 10/md.constants.yts*ones(md.nsteps,1)];
 md.cm_min=[-2000*ones(md.numberofnodes,1) -2000*ones(md.numberofnodes,1)];
 md.cm_max=[+2000*ones(md.numberofnodes,1) +2000*ones(md.numberofnodes,1)];

issm/trunk/test/Par/ISMIPC.par

@@ -7 +7 @@
 
 disp('      creating drag');
-%md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x/2)).*sin(md.y*2*pi/max(md.x/2)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
-md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x))));
+%md.drag_coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.x*2*pi/max(md.x/2)).*sin(md.y*2*pi/max(md.x/2)))./(md.constants.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
+md.drag_coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x))));
 %Take care of iceshelves: no basal drag
 pos=find(md.elementoniceshelf);

issm/trunk/test/Par/ISMIPD.par

@@ -7 +7 @@
 
 disp('      creating drag');
-md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x))));
+md.drag_coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.x*2*pi/max(md.x))));
 %Take care of iceshelves: no basal drag
 pos=find(md.elementoniceshelf);

issm/trunk/test/Par/ISMIPF.par

@@ -9 +9 @@
 
 disp('      creating drag');
-md.drag_coefficient=sqrt(md.yts/(2.140373*10^-7*1000))*ones(md.numberofnodes,1);
+md.drag_coefficient=sqrt(md.constants.yts/(2.140373*10^-7*1000))*ones(md.numberofnodes,1);
 md.drag_p=ones(md.numberofelements,1);
 md.drag_q=zeros(md.numberofelements,1);

issm/trunk/test/Par/RoundSheetEISMINT.par

@@ -56 +56 @@
 md.thermalconductivity=2.1;
 md.latentheat=3.35*10^5;
-md.beta=8.66*10^-4/(md.rho_ice*md.g); %conversion from K/m to K/Pa
-md.yts=31556926;
+md.beta=8.66*10^-4/(md.rho_ice*md.constants.g); %conversion from K/m to K/Pa
+md.constants.yts=31556926;

issm/trunk/test/Par/SquareThermal.par

@@ -30 +30 @@
 
 disp('      creating surface_mass_balance');
-md.surface_mass_balance=ones(md.numberofnodes,1)/md.yts; %1m/a
-md.basal_melting_rate=0*ones(md.numberofnodes,1)/md.yts; %1m/a
+md.surface_mass_balance=ones(md.numberofnodes,1)/md.constants.yts; %1m/a
+md.basal_melting_rate=0*ones(md.numberofnodes,1)/md.constants.yts; %1m/a
 
 %Deal with boundary conditions: