Index: issm/trunk-jpl/src/m/classes/basalforcings.py
===================================================================
--- issm/trunk-jpl/src/m/classes/basalforcings.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/basalforcings.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class basalforcings:
@@ -28,2 +31,23 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		if PrognosticAnalysisEnum in analyses and not (solution==TransientSolutionEnum and not md.transient.isprognostic):
+			md = checkfield(md,'basalforcings.melting_rate','NaN',1,'forcing',1)
+
+		if BalancethicknessAnalysisEnum in analyses:
+			md = checkfield(md,'basalforcings.melting_rate','NaN',1,'size',[md.mesh.numberofvertices])
+
+		if ThermalAnalysisEnum in analyses and not (solution==TransientSolutionEnum and not md.transient.isthermal):
+			md = checkfield(md,'basalforcings.melting_rate','NaN',1,'forcing',1)
+			md = checkfield(md,'basalforcings.geothermalflux','NaN',1,'forcing',1,'>=',0)
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','melting_rate','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','melting_rate_correction','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','geothermalflux','format','DoubleMat','mattype',1)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/constants.py
===================================================================
--- issm/trunk-jpl/src/m/classes/constants.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/constants.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class constants:
@@ -40,2 +43,17 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		md = checkfield(md,'constants.g','>',0,'size',[1])
+		md = checkfield(md,'constants.yts','>',0,'size',[1])
+		md = checkfield(md,'constants.referencetemperature','size',[1])
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','g','format','Double')
+		WriteData(fid,'object',self,'fieldname','yts','format','Double')
+		WriteData(fid,'object',self,'fieldname','referencetemperature','format','Double')
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- issm/trunk-jpl/src/m/classes/flowequation.m	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/flowequation.m	(revision 12949)
@@ -44,5 +44,5 @@
 					md = checkfield(md,'flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[0:7]);
 				end
-				if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0),
+				if ~(md.flowequation.ismacayealpattyn || md.flowequation.ishutter || md.flowequation.isstokes),
 					md = checkmessage(md,['no element types set for this model. At least one of ismacayealpattyn, ishutter or isstokes need to be =1']);
 				end
Index: issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- issm/trunk-jpl/src/m/classes/flowequation.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/flowequation.py	(revision 12949)
@@ -1,4 +1,8 @@
 #module imports
+import numpy
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class flowequation:
@@ -40,2 +44,58 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		if DiagnosticHorizAnalysisEnum in analyses:
+			md = checkfield(md,'flowequation.ismacayealpattyn','numel',1,'values',[0,1])
+			md = checkfield(md,'flowequation.ishutter','numel',1,'values',[0,1])
+			md = checkfield(md,'flowequation.isstokes','numel',1,'values',[0,1])
+			md = checkfield(md,'flowequation.bordermacayeal','size',[md.mesh.numberofvertices],'values',[0,1])
+			md = checkfield(md,'flowequation.borderpattyn','size',[md.mesh.numberofvertices],'values',[0,1])
+			md = checkfield(md,'flowequation.borderstokes','size',[md.mesh.numberofvertices],'values',[0,1])
+			if md.mesh.dimension==2:
+				md = checkfield(md,'flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[1,2])
+				md = checkfield(md,'flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2])
+			else:
+				md = checkfield(md,'flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',range(0,7+1))
+				md = checkfield(md,'flowequation.element_equation','size',[md.mesh.numberofelements],'values',range(0,7+1))
+			if not (md.flowequation.ismacayealpattyn or md.flowequation.ishutter or md.flowequation.isstokes):
+				md.checkmessage("no element types set for this model. At least one of ismacayealpattyn, ishutter or isstokes need to be =1")
+
+		if DiagnosticHutterAnalysisEnum in analyses:
+			if any(md.flowequation.element_equation==1):
+				if any(numpy.logical_and(md.flowequation.element_equation,md.mask.elementonfloatingice)):
+					print "\n !!! Warning: Hutter's model is not consistent on ice shelves !!!\n"
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','ismacayealpattyn','format','Boolean')
+		WriteData(fid,'object',self,'fieldname','ishutter','format','Boolean')
+		WriteData(fid,'object',self,'fieldname','isstokes','format','Boolean')
+		WriteData(fid,'object',self,'fieldname','bordermacayeal','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','borderpattyn','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','borderstokes','format','DoubleMat','mattype',1)
+		#convert approximations to enums
+		data=self.vertex_equation
+		data[[i for i,item in enumerate(data) if item==0]]=NoneApproximationEnum
+		data[[i for i,item in enumerate(data) if item==1]]=HutterApproximationEnum
+		data[[i for i,item in enumerate(data) if item==2]]=MacAyealApproximationEnum
+		data[[i for i,item in enumerate(data) if item==3]]=PattynApproximationEnum
+		data[[i for i,item in enumerate(data) if item==4]]=StokesApproximationEnum
+		data[[i for i,item in enumerate(data) if item==5]]=MacAyealPattynApproximationEnum
+		data[[i for i,item in enumerate(data) if item==6]]=MacAyealStokesApproximationEnum
+		data[[i for i,item in enumerate(data) if item==7]]=PattynStokesApproximationEnum
+		WriteData(fid,'data',data,'enum',FlowequationVertexEquationEnum,'format','DoubleMat','mattype',1)
+		data=self.element_equation
+		data[[i for i,item in enumerate(data) if item==0]]=NoneApproximationEnum
+		data[[i for i,item in enumerate(data) if item==1]]=HutterApproximationEnum
+		data[[i for i,item in enumerate(data) if item==2]]=MacAyealApproximationEnum
+		data[[i for i,item in enumerate(data) if item==3]]=PattynApproximationEnum
+		data[[i for i,item in enumerate(data) if item==4]]=StokesApproximationEnum
+		data[[i for i,item in enumerate(data) if item==5]]=MacAyealPattynApproximationEnum
+		data[[i for i,item in enumerate(data) if item==6]]=MacAyealStokesApproximationEnum
+		data[[i for i,item in enumerate(data) if item==7]]=PattynStokesApproximationEnum
+		WriteData(fid,'data',data,'enum',FlowequationElementEquationEnum,'format','DoubleMat','mattype',2)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/friction.py
===================================================================
--- issm/trunk-jpl/src/m/classes/friction.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/friction.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class friction:
@@ -28,2 +31,21 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		#Early return
+		if not DiagnosticHorizAnalysisEnum in analyses and not ThermalAnalysisEnum in analyses:
+			return md
+
+		md = checkfield(md,'friction.coefficient','NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'friction.q','NaN',1,'size',[md.mesh.numberofelements])
+		md = checkfield(md,'friction.p','NaN',1,'size',[md.mesh.numberofelements])
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','p','format','DoubleMat','mattype',2)
+		WriteData(fid,'object',self,'fieldname','q','format','DoubleMat','mattype',2)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/geometry.py
===================================================================
--- issm/trunk-jpl/src/m/classes/geometry.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/geometry.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class geometry:
@@ -34,2 +37,24 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		md = checkfield(md,'geometry.surface'  ,'NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'geometry.bed'      ,'NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'geometry.thickness','NaN',1,'size',[md.mesh.numberofvertices],'>',0)
+		if any((self.thickness-self.surface+self.bed)>10**-9):
+			md.checkmessage("equality thickness=surface-bed violated")
+		end 
+		if solution==TransientSolutionEnum and md.transient.isgroundingline:
+			md = checkfield(md,'geometry.bathymetry','NaN',1,'size',[md.mesh.numberofvertices])
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'data',self.surface,'format','DoubleMat','mattype',1,'enum',SurfaceEnum)
+		WriteData(fid,'data',self.thickness,'format','DoubleMat','mattype',1,'enum',ThicknessEnum)
+		WriteData(fid,'data',self.bed,'format','DoubleMat','mattype',1,'enum',BedEnum)
+		WriteData(fid,'data',self.bathymetry,'format','DoubleMat','mattype',1,'enum',BathymetryEnum)
+		WriteData(fid,'object',self,'fieldname','hydrostatic_ratio','format','DoubleMat','mattype',1)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/initialization.m
===================================================================
--- issm/trunk-jpl/src/m/classes/initialization.m	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/initialization.m	(revision 12949)
@@ -29,5 +29,5 @@
 		function md = checkconsistency(obj,md,solution,analyses) % {{{
 			if ismember(DiagnosticHorizAnalysisEnum,analyses)
-				if ~isnan(md.initialization.vx) & ~isnan(md.initialization.vy),
+				if ~(isnan(md.initialization.vx) | isnan(md.initialization.vy)),
 					md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
 					md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
Index: issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- issm/trunk-jpl/src/m/classes/initialization.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/initialization.py	(revision 12949)
@@ -1,4 +1,8 @@
 #module imports
+import numpy
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class initialization:
@@ -41,2 +45,40 @@
 	#}}}
 
+	def checkconsistency(obj,md,solution,analyses):    # {{{
+		if DiagnosticHorizAnalysisEnum in analyses:
+			if not any(numpy.isnan(md.initialization.vx) or numpy.isnan(md.initialization.vy)):
+				md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
+				md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
+		if PrognosticAnalysisEnum in analyses:
+			md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
+			md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
+		if HydrologyAnalysisEnum in analyses:
+			md = checkfield(md,'initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices])
+		if BalancethicknessAnalysisEnum in analyses:
+			md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
+			md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
+			#Triangle with zero velocity
+			if any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements.astype(int)-1]),axis=1)==0,\
+			                         numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements.astype(int)-1]),axis=1)==0)):
+				md.checkmessage("at least one triangle has all its vertices with a zero velocity")
+		if ThermalAnalysisEnum in analyses:
+			md = checkfield(md,'initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
+			md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
+			md = checkfield(md,'initialization.vz','NaN',1,'size',[md.mesh.numberofvertices])
+			md = checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices])
+		if (EnthalpyAnalysisEnum in analyses and md.thermal.isenthalpy) or solution==EnthalpySolutionEnum:
+			md = checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
+
+		return md
+	# }}}
+
+	def marshall(obj,fid):    # {{{
+		WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum)
+		WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum)
+		WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum)
+		WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum)
+		WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum)
+		WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum)
+		WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/mask.py
===================================================================
--- issm/trunk-jpl/src/m/classes/mask.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/mask.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class mask:
@@ -35,2 +38,23 @@
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		md = checkfield(md,'mask.elementonfloatingice','size',[md.mesh.numberofelements],'values',[0,1])
+		md = checkfield(md,'mask.elementongroundedice','size',[md.mesh.numberofelements],'values',[0,1])
+		md = checkfield(md,'mask.elementonwater'      ,'size',[md.mesh.numberofelements],'values',[0,1])
+		md = checkfield(md,'mask.vertexonfloatingice' ,'size',[md.mesh.numberofvertices],'values',[0,1])
+		md = checkfield(md,'mask.vertexongroundedice' ,'size',[md.mesh.numberofvertices],'values',[0,1])
+		md = checkfield(md,'mask.vertexonwater'       ,'size',[md.mesh.numberofvertices],'values',[0,1])
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','elementonfloatingice','format','BooleanMat','mattype',2)
+		WriteData(fid,'object',self,'fieldname','elementongroundedice','format','BooleanMat','mattype',2)
+		WriteData(fid,'object',self,'fieldname','elementonwater','format','BooleanMat','mattype',2)
+		WriteData(fid,'object',self,'fieldname','vertexonfloatingice','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','vertexongroundedice','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','vertexonwater','format','DoubleMat','mattype',1)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/materials.py
===================================================================
--- issm/trunk-jpl/src/m/classes/materials.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/materials.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class materials:
@@ -6,17 +9,17 @@
 	def __init__(self):
 		# {{{ Properties
-		self.rho_ice                    = 0;
-		self.rho_water                  = 0;
-		self.mu_water                   = 0;
-		self.heatcapacity               = 0;
-		self.latentheat                 = 0;
-		self.thermalconductivity        = 0;
-		self.meltingpoint               = 0;
-		self.beta                       = 0;
-		self.mixed_layer_capacity       = 0;
-		self.thermal_exchange_velocity  = 0;
+		self.rho_ice                    = 0
+		self.rho_water                  = 0
+		self.mu_water                   = 0
+		self.heatcapacity               = 0
+		self.latentheat                 = 0
+		self.thermalconductivity        = 0
+		self.meltingpoint               = 0
+		self.beta                       = 0
+		self.mixed_layer_capacity       = 0
+		self.thermal_exchange_velocity  = 0
 		self.rheology_B   = float('NaN')
 		self.rheology_n   = float('NaN')
-		self.rheology_law = "";
+		self.rheology_law = ""
 
 		self.setdefaultparameters()
@@ -79,2 +82,31 @@
 		return obj
 		#}}}
+
+	def checkconsistency(self,md,solution,analyses):    # {{{
+		md = checkfield(md,'materials.rho_ice','>',0)
+		md = checkfield(md,'materials.rho_water','>',0)
+		md = checkfield(md,'materials.rho_freshwater','>',0)
+		md = checkfield(md,'materials.mu_water','>',0)
+		md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements])
+		md = checkfield(md,'materials.rheology_law','values',['None','Paterson','Arrhenius'])
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','rho_ice','format','Double')
+		WriteData(fid,'object',self,'fieldname','rho_water','format','Double')
+		WriteData(fid,'object',self,'fieldname','rho_freshwater','format','Double')
+		WriteData(fid,'object',self,'fieldname','mu_water','format','Double')
+		WriteData(fid,'object',self,'fieldname','heatcapacity','format','Double')
+		WriteData(fid,'object',self,'fieldname','latentheat','format','Double')
+		WriteData(fid,'object',self,'fieldname','thermalconductivity','format','Double')
+		WriteData(fid,'object',self,'fieldname','meltingpoint','format','Double')
+		WriteData(fid,'object',self,'fieldname','beta','format','Double')
+		WriteData(fid,'object',self,'fieldname','mixed_layer_capacity','format','Double')
+		WriteData(fid,'object',self,'fieldname','thermal_exchange_velocity','format','Double')
+		WriteData(fid,'object',self,'fieldname','rheology_B','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','rheology_n','format','DoubleMat','mattype',2)
+		WriteData(fid,'data',StringToEnum(self.rheology_law),'enum',MaterialsRheologyLawEnum,'format','Integer')
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/mesh.py
===================================================================
--- issm/trunk-jpl/src/m/classes/mesh.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/mesh.py	(revision 12949)
@@ -2,6 +2,6 @@
 import numpy
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
 from checkfield import *
-from EnumDefinitions import *
 from MatlabFuncs import *
 
@@ -131,5 +131,5 @@
 			md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,6])
 		if any(numpy.logical_not(ismember(range(1,md.mesh.numberofvertices+1),md.mesh.elements))):
-			md = checkmessage(md,"orphan nodes have been found. Check the mesh outline")
+			md.checkmessage("orphan nodes have been found. Check the mesh outline")
 		md = checkfield(md,'mesh.dimension','values',[2,3])
 		md = checkfield(md,'mesh.numberoflayers','>=',0)
Index: issm/trunk-jpl/src/m/classes/private.py
===================================================================
--- issm/trunk-jpl/src/m/classes/private.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/private.py	(revision 12949)
@@ -6,8 +6,8 @@
 	def __init__(self):
 		# {{{ Properties
-		self.isconsistent = True;
+		self.isconsistent = True
 		self.runtimename  = ''
 		self.bamg         = {}
-		self.solution     = '';
+		self.solution     = ''
 
 		#set defaults
Index: issm/trunk-jpl/src/m/classes/surfaceforcings.m
===================================================================
--- issm/trunk-jpl/src/m/classes/surfaceforcings.m	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/surfaceforcings.m	(revision 12949)
@@ -11,5 +11,5 @@
 		issmbgradients = 0;
 		isdelta18o = 0;
-	   hc = NaN;
+		hc = NaN;
 		smb_pos_max = NaN;
 		smb_pos_min = NaN;
@@ -46,24 +46,24 @@
 			if ismember(PrognosticAnalysisEnum,analyses),
 				md = checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0 1]);
-				checkfield(md,'surfaceforcings.issmbgradients','numel',1,'values',[0 1]);
+				md = checkfield(md,'surfaceforcings.issmbgradients','numel',1,'values',[0 1]);
 				if(obj.ispdd)
-				        if(obj.isdelta18o==0)
-				        	md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
-				        	md = checkfield(md,'surfaceforcings.precipitation','forcing',1,'NaN',1);
-			                else
-				                md = checkfield(md,'surfaceforcings.delta18o','NaN',1);
-				        	md = checkfield(md,'surfaceforcings.delta18o_surface','NaN',1);
-				        	md = checkfield(md,'surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-				        	md = checkfield(md,'surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-				        	md = checkfield(md,'surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-				        end
+					if(obj.isdelta18o==0)
+						md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
+						md = checkfield(md,'surfaceforcings.precipitation','forcing',1,'NaN',1);
+					else
+						md = checkfield(md,'surfaceforcings.delta18o','NaN',1);
+						md = checkfield(md,'surfaceforcings.delta18o_surface','NaN',1);
+						md = checkfield(md,'surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+						md = checkfield(md,'surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+						md = checkfield(md,'surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+					end
 				elseif(obj.issmbgradients)
-					checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.a_pos','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.b_pos','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.a_neg','forcing',1,'NaN',1);
-					checkfield(md,'surfaceforcings.b_neg','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.a_pos','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.b_pos','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.a_neg','forcing',1,'NaN',1);
+					md = checkfield(md,'surfaceforcings.b_neg','forcing',1,'NaN',1);
 				else
 					md = checkfield(md,'surfaceforcings.mass_balance','forcing',1,'NaN',1);
@@ -103,14 +103,14 @@
 			WriteData(fid,'object',obj,'fieldname','isdelta18o','format','Boolean');
 			if obj.ispdd,
-			  if obj.isdelta18o
-			        WriteData(fid,'object',obj,'fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
-			        WriteData(fid,'object',obj,'fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
-			        WriteData(fid,'object',obj,'fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
-			        WriteData(fid,'object',obj,'fieldname','delta18o_surface','format','DoubleMat','mattype',1);
-			        WriteData(fid,'object',obj,'fieldname','delta18o','format','DoubleMat','mattype',1);
-			  else
-				WriteData(fid,'object',obj,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
-				WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1);
-			  end
+				if obj.isdelta18o
+					WriteData(fid,'object',obj,'fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
+					WriteData(fid,'object',obj,'fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
+					WriteData(fid,'object',obj,'fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
+					WriteData(fid,'object',obj,'fieldname','delta18o_surface','format','DoubleMat','mattype',1);
+					WriteData(fid,'object',obj,'fieldname','delta18o','format','DoubleMat','mattype',1);
+				else
+					WriteData(fid,'object',obj,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
+					WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1);
+				end
 			end
 			WriteData(fid,'object',obj,'fieldname','issmbgradients','format','Boolean');
Index: issm/trunk-jpl/src/m/classes/surfaceforcings.py
===================================================================
--- issm/trunk-jpl/src/m/classes/surfaceforcings.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/surfaceforcings.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class surfaceforcings:
@@ -8,4 +11,20 @@
 		self.precipitation = float('NaN')
 		self.mass_balance  = float('NaN')
+		self.ispdd = 0
+		self.issmbgradients = 0
+		self.isdelta18o = 0
+		self.hc = float('NaN')
+		self.smb_pos_max = float('NaN')
+		self.smb_pos_min = float('NaN')
+		self.a_pos = float('NaN')
+		self.b_pos = float('NaN')
+		self.a_neg = float('NaN')
+		self.b_neg = float('NaN')
+		self.monthlytemperatures = float('NaN')
+		self.delta18o = float('NaN')
+		self.delta18o_surface = float('NaN')
+		self.temperatures_presentday = float('NaN')
+		self.temperatures_lgm = float('NaN')
+		self.precipitations_presentday = float('NaN')
 
 		#set defaults
@@ -19,4 +38,21 @@
 		string="%s\n\n%s"%(string,fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]'))
 		string="%s\n%s"%(string,fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]'))
+		string="%s\n%s"%(string,fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)'))
+		string="%s\n%s"%(string,fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'))
+		string="%s\n%s"%(string,fielddisplay(obj,'monthlytemperatures','monthly surface temperatures [Kelvin], required if pdd is activated and delta18o not activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]'))
+		string="%s\n%s"%(string,fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [Kelvin], required if pdd is activated and delta18o activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [Kelvin], required if pdd is activated and delta18o activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)'))
+		string="%s\n%s"%(string,fielddisplay(obj,'hc',' elevation of intersection between accumulation and ablation regime required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'smb_pos_max',' maximum value of positive smb required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'smb_pos_min',' minimum value of positive smb required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'a_pos',' intercept of hs - smb regression line for accumulation regime required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'b_pos',' slope of hs - smb regression line for accumulation regime required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'a_neg',' intercept of hs - smb regression line for ablation regime required if smb gradients is activated'))
+		string="%s\n%s"%(string,fielddisplay(obj,'b_neg',' slope of hs - smb regression line for ablation regime required if smb gradients is activated'))
 
 		return string
@@ -25,5 +61,72 @@
 	def setdefaultparameters(obj):
 		# {{{setdefaultparameters
+		  
+		#pdd method not used in default mode
+		obj.ispdd=0
+		obj.issmbgradients=0
+		obj.isdelta18o=0
+
 		return obj
 	#}}}
 
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		if PrognosticAnalysisEnum in analyses:
+			md = checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0,1])
+			md = checkfield(md,'surfaceforcings.issmbgradients','numel',1,'values',[0,1])
+			if   self.ispdd:
+				if not self.isdelta18o:
+					md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1)
+					md = checkfield(md,'surfaceforcings.precipitation','forcing',1,'NaN',1)
+				else:
+					md = checkfield(md,'surfaceforcings.delta18o','NaN',1)
+					md = checkfield(md,'surfaceforcings.delta18o_surface','NaN',1)
+					md = checkfield(md,'surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
+					md = checkfield(md,'surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1)
+					md = checkfield(md,'surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
+			elif self.issmbgradients:
+				md = checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.a_pos','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.b_pos','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.a_neg','forcing',1,'NaN',1)
+				md = checkfield(md,'surfaceforcings.b_neg','forcing',1,'NaN',1)
+			else:
+				md = checkfield(md,'surfaceforcings.mass_balance','forcing',1,'NaN',1)
+
+		if BalancethicknessAnalysisEnum in analyses:
+			md = checkfield(md,'surfaceforcings.mass_balance','size',[md.mesh.numberofvertices],'NaN',1)
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','precipitation','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','mass_balance','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','ispdd','format','Boolean')
+		WriteData(fid,'object',self,'fieldname','isdelta18o','format','Boolean')
+
+		if self.ispdd:
+			if self.isdelta18o:
+				WriteData(fid,'object',self,'fieldname','temperatures_presentday','format','DoubleMat','mattype',1)
+				WriteData(fid,'object',self,'fieldname','temperatures_lgm','format','DoubleMat','mattype',1)
+				WriteData(fid,'object',self,'fieldname','precipitations_presentday','format','DoubleMat','mattype',1)
+				WriteData(fid,'object',self,'fieldname','delta18o_surface','format','DoubleMat','mattype',1)
+				WriteData(fid,'object',self,'fieldname','delta18o','format','DoubleMat','mattype',1)
+			else:
+				WriteData(fid,'object',self,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1)
+				WriteData(fid,'object',self,'fieldname','precipitation','format','DoubleMat','mattype',1)
+
+		WriteData(fid,'object',self,'fieldname','issmbgradients','format','Boolean')
+
+		if self.issmbgradients:
+			WriteData(fid,'object',self,'fieldname','hc','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','smb_pos_max','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','smb_pos_min','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','a_pos','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','b_pos','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','a_neg','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'fieldname','b_neg','format','DoubleMat','mattype',1)
+	# }}}
+
Index: issm/trunk-jpl/src/m/classes/timestepping.py
===================================================================
--- issm/trunk-jpl/src/m/classes/timestepping.py	(revision 12945)
+++ issm/trunk-jpl/src/m/classes/timestepping.py	(revision 12949)
@@ -1,4 +1,7 @@
 #module imports
 from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from WriteData import *
 
 class timestepping:
@@ -6,6 +9,7 @@
 	def __init__(self):
 		# {{{ Properties
+		self.start_time      = 0
+		self.final_time      = 0
 		self.time_step       = 0
-		self.final_time      = 0
 		self.time_adapt      = 0
 		self.cfl_coefficient = 0
@@ -18,6 +22,7 @@
 		# {{{ Display
 		string="   timestepping parameters:"
-		string="%s\n\n%s"%(string,fielddisplay(obj,"time_step","length of time steps [yrs]"))
+		string="%s\n\n%s"%(string,fielddisplay(obj,"start_time","simulation starting time [yrs]"))
 		string="%s\n%s"%(string,fielddisplay(obj,"final_time","final time to stop the simulation [yrs]"))
+		string="%s\n%s"%(string,fielddisplay(obj,"time_step","length of time steps [yrs]"))
 		string="%s\n%s"%(string,fielddisplay(obj,"time_adapt","use cfl condition to define time step ? (0 or 1) "))
 		string="%s\n%s"%(string,fielddisplay(obj,"cfl_coefficient","coefficient applied to cfl condition"))
@@ -40,2 +45,24 @@
 		return obj
 	#}}}
+
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		md = checkfield(md,'timestepping.start_time','numel',1,'NaN',1)
+		md = checkfield(md,'timestepping.final_time','numel',1,'NaN',1)
+		md = checkfield(md,'timestepping.time_step','numel',1,'>=',0,'NaN',1)
+		md = checkfield(md,'timestepping.time_adapt','numel',1,'values',[0,1])
+		md = checkfield(md,'timestepping.cfl_coefficient','numel',1,'>',0,'<=',1)
+		if self.final_time-self.start_time<0:
+			md.checkmessage("timestepping.final_time should be larger than timestepping.start_time")
+
+		return md
+	# }}}
+
+	def marshall(self,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','start_time','format','Double')
+		WriteData(fid,'object',self,'fieldname','final_time','format','Double')
+		WriteData(fid,'object',self,'fieldname','time_step','format','Double')
+		WriteData(fid,'object',self,'fieldname','time_adapt','format','Boolean')
+		WriteData(fid,'object',self,'fieldname','cfl_coefficient','format','Double')
+	# }}}
+
Index: issm/trunk-jpl/src/m/model/setflowequation.m
===================================================================
--- issm/trunk-jpl/src/m/model/setflowequation.m	(revision 12945)
+++ issm/trunk-jpl/src/m/model/setflowequation.m	(revision 12949)
@@ -35,9 +35,9 @@
 %Flag the elements that have not been flagged as filltype
 if strcmpi(filltype,'hutter'),
-	hutterflag(find(~macayealflag & ~pattynflag))=1;
+	hutterflag(find(~(macayealflag | pattynflag)))=1;
 elseif strcmpi(filltype,'macayeal'),
-	macayealflag(find(~hutterflag & ~pattynflag & ~stokesflag))=1;
+	macayealflag(find(~(hutterflag | pattynflag | ~stokesflag)))=1;
 elseif strcmpi(filltype,'pattyn'),
-	pattynflag(find(~hutterflag & ~macayealflag & ~stokesflag))=1;
+	pattynflag(find(~(hutterflag | ~macayealflag | ~stokesflag)))=1;
 end
 
Index: issm/trunk-jpl/src/m/model/solve.py
===================================================================
--- issm/trunk-jpl/src/m/model/solve.py	(revision 12945)
+++ issm/trunk-jpl/src/m/model/solve.py	(revision 12949)
@@ -16,14 +16,14 @@
  
 	   solution types available comprise:
-		  - DiagnosticSolutionEnum
-	 	  - PrognosticSolutionEnum
-	 	  - ThermalSolutionEnum
-	 	  - SteadystateSolutionEnum
-	 	  - TransientSolutionEnum...
-	 	  - BalancethicknessSolutionEnum
-	 	  - BedSlopeSolutionEnum
-	 	  - SurfaceSlopeSolutionEnum
-	 	  - HydrologySolutionEnum
-	 	  - FlaimSolutionEnum
+	      - DiagnosticSolutionEnum
+	      - PrognosticSolutionEnum
+	      - ThermalSolutionEnum
+	      - SteadystateSolutionEnum
+	      - TransientSolutionEnum...
+	      - BalancethicknessSolutionEnum
+	      - BedSlopeSolutionEnum
+	      - SurfaceSlopeSolutionEnum
+	      - HydrologySolutionEnum
+	      - FlaimSolutionEnum
  
 	   extra options:
