Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 24145)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 24146)
@@ -60,4 +60,14 @@
 	iomodel->DeleteData(1,"md.flowequation.element_equation");
 
+   /*First, reset all F to 0 */
+	for(int i=0;i<elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+		int numvertices = element->GetNumberOfVertices();
+		IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
+		element->AddInput(DamageFEnum,values,P1Enum);
+		xDelete<IssmDouble>(values);
+	}
+
+
 	/*What input do I need to run my damage evolution model?*/
 	iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
@@ -97,5 +107,7 @@
 		parameters->AddObject(iomodel->CopyConstantObject("md.damage.c2",DamageC2Enum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.damage.c3",DamageC3Enum));
+		parameters->AddObject(iomodel->CopyConstantObject("md.damage.c4",DamageC4Enum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_threshold",DamageStressThresholdEnum));
+		parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_ubound",DamageStressUBoundEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.damage.kappa",DamageKappaEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.damage.healing",DamageHealingEnum));
@@ -130,4 +142,89 @@
 	xDelete<IssmDouble>(f);
 }/*}}}*/
+void           DamageEvolutionAnalysis::CreateDamageFInputArctan(Element* element){/*{{{*/
+	IssmDouble c1, c2, stress_threshold, stress_ubound;
+	IssmDouble damage;
+	IssmDouble yts;
+	IssmDouble principalDevStress1, principalDevStress2;
+	IssmDouble tensileStress, compressiveStress;
+	
+	int equivstress, domaintype, dim;
+
+	/*Fetch number of vertices and allocate output*/
+	int numnodes = element->GetNumberOfNodes();
+	IssmDouble* f   = xNew<IssmDouble>(numnodes);
+
+	/*retrieve parameters:*/
+	element->FindParam(&c1,DamageC1Enum);
+	element->FindParam(&c2,DamageC2Enum);
+	element->FindParam(&yts,ConstantsYtsEnum);
+	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
+	element->FindParam(&stress_ubound,DamageStressUBoundEnum);
+	element->FindParam(&domaintype,DomainTypeEnum);
+
+	/*Get problem dimension*/
+	switch(domaintype){
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 3; break;
+		default: _error_("not implemented");
+	}
+	/*Compute stress tensor and Stress Max Principal: */
+	element->ComputeDeviatoricStressTensor();
+
+	Input* principalDevStress1_input = element->GetInput(DeviatoricStress1Enum);     _assert_(principalDevStress1_input);
+	Input* principalDevStress2_input = element->GetInput(DeviatoricStress2Enum);     _assert_(principalDevStress2_input);
+
+	Input* damage_input = NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		damage_input = element->GetInput(DamageDbarEnum); 	_assert_(damage_input);
+	}
+	else{
+		damage_input = element->GetInput(DamageDEnum);   _assert_(damage_input);
+	}
+
+	/*Calculate damage evolution source term */
+	Gauss* gauss=element->NewGauss();
+
+	/* To keep arctan output (bounded by -pi/2 and pi/2) within the specified boundaries */
+	c1 /= (PI/2);
+	c2 /= (PI/2);
+	/* To have per second output with per annum parameters */
+	c1 /= yts;
+	c2 /= yts;
+
+	for (int i=0;i<numnodes;i++){
+		f[i] = 0;
+		
+		gauss->GaussNode(element->GetElementType(),i);
+
+		damage_input->GetInputValue(&damage,gauss);
+		principalDevStress1_input->GetInputValue(&principalDevStress1,gauss);
+		principalDevStress2_input->GetInputValue(&principalDevStress2,gauss);
+
+		tensileStress     = sqrt(1.5*(pow(max(principalDevStress1, 0.), 2) + pow(max(principalDevStress2, 0.), 2)));
+		compressiveStress = sqrt(1.5*(pow(min(principalDevStress1, 0.), 2) + pow(min(principalDevStress2, 0.), 2)));
+
+		/* Calculate principal effective stresses */
+		if(dim==2){
+			f[i] = 0;
+			if(tensileStress > stress_threshold)
+				f[i] += c1*atan((tensileStress/stress_threshold - 1)/(1-damage));
+
+			if(compressiveStress < stress_ubound)
+				f[i] += c2*atan((compressiveStress/stress_ubound - 1)/(1-damage));
+		}
+		else{
+			_error_("Only 2D is implemented.");
+		}
+	}
+
+	/*Add input*/
+	element->AddInput(DamageFEnum,f,element->GetElementType());
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(f);
+	delete gauss;
+}/*}}}*/
+
 void           DamageEvolutionAnalysis::CreateDamageFInputExp(Element* element){/*{{{*/
 
@@ -191,4 +288,8 @@
 		}
 		else f[i]=0;
+
+		/*Edits from MM*/
+		if(f[i]>10.) f[i]=10.;
+		if(f[i]<-10.) f[i]=-10.;
 	}
 
@@ -358,4 +459,5 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	//printf("dt %f\n", dt);
 	element->FindParam(&stabilization,DamageStabilizationEnum);
 	Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
@@ -525,4 +627,7 @@
 			this->CreateDamageFInputExp(element);
 			break;
+		case 3:
+			this->CreateDamageFInputArctan(element);
+			break;
 		default:
 			_error_("not implemented yet");
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 24145)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 24146)
@@ -23,4 +23,5 @@
 		void           Core(FemModel* femmodel);
 		void           CreateDamageFInput(Element* element);
+		void           CreateDamageFInputArctan(Element* element);
 		void           CreateDamageFInputExp(Element* element);
 		void           CreateDamageFInputPralong(Element* element);
Index: /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp	(revision 24145)
+++ /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp	(revision 24146)
@@ -24,4 +24,8 @@
 			/* Damage calculated using source term in DamageEvolutionAnalysis */
 			break;
+		case 3:
+			if(VerboseModule()) _printf0_("   computing damage using source term in advection scheme\n");
+			/* Damage calculated using source term in DamageEvolutionAnalysis */
+			break;
 		default:
 			_error_("Damage law "<<EnumToStringx(damagelaw)<<" not implemented yet");
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 24145)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 24146)
@@ -109,4 +109,5 @@
 	DamageC2Enum,
 	DamageC3Enum,
+	DamageC4Enum,
 	DamageEnum,
 	DamageEquivStressEnum,
@@ -119,4 +120,5 @@
 	DamageStabilizationEnum,
 	DamageStressThresholdEnum,
+	DamageStressUBoundEnum,
 	DebugProfilingEnum,
 	DomainDimensionEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 24145)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 24146)
@@ -117,4 +117,5 @@
 		case DamageC2Enum : return "DamageC2";
 		case DamageC3Enum : return "DamageC3";
+		case DamageC4Enum : return "DamageC4";
 		case DamageEnum : return "Damage";
 		case DamageEquivStressEnum : return "DamageEquivStress";
@@ -127,4 +128,5 @@
 		case DamageStabilizationEnum : return "DamageStabilization";
 		case DamageStressThresholdEnum : return "DamageStressThreshold";
+		case DamageStressUBoundEnum : return "DamageStressUBound";
 		case DebugProfilingEnum : return "DebugProfiling";
 		case DomainDimensionEnum : return "DomainDimension";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 24145)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 24146)
@@ -117,4 +117,5 @@
 	      else if (strcmp(name,"DamageC2")==0) return DamageC2Enum;
 	      else if (strcmp(name,"DamageC3")==0) return DamageC3Enum;
+	      else if (strcmp(name,"DamageC4")==0) return DamageC4Enum;
 	      else if (strcmp(name,"Damage")==0) return DamageEnum;
 	      else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
@@ -127,4 +128,5 @@
 	      else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
 	      else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
+	      else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum;
 	      else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
 	      else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
@@ -135,10 +137,10 @@
 	      else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
 	      else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
-	      else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
-	      else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
          else stage=2;
    }
    if(stage==2){
-	      if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum;
+	      if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
+	      else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
+	      else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum;
 	      else if (strcmp(name,"FlowequationFeFS")==0) return FlowequationFeFSEnum;
 	      else if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum;
@@ -258,10 +260,10 @@
 	      else if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum;
 	      else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum;
-	      else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum;
-	      else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
+	      if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum;
+	      else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
+	      else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
 	      else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
 	      else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
@@ -381,10 +383,10 @@
 	      else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
 	      else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
-	      else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
-	      else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum;
+	      if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
+	      else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
+	      else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum;
 	      else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum;
 	      else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
@@ -504,10 +506,10 @@
 	      else if (strcmp(name,"DamageD")==0) return DamageDEnum;
 	      else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
-	      else if (strcmp(name,"DamageF")==0) return DamageFEnum;
-	      else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
+	      if (strcmp(name,"DamageF")==0) return DamageFEnum;
+	      else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
+	      else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
 	      else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
 	      else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
@@ -627,10 +629,10 @@
 	      else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
 	      else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum;
-	      else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
-	      else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
+	      if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
+	      else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
+	      else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
 	      else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
 	      else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
@@ -750,10 +752,10 @@
 	      else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
 	      else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
-	      else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
-	      else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
+	      if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
+	      else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
+	      else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
 	      else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
 	      else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
@@ -873,10 +875,10 @@
 	      else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
 	      else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
-	      else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
-	      else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
+	      if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
+	      else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
+	      else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
 	      else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
 	      else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
@@ -996,10 +998,10 @@
 	      else if (strcmp(name,"Dense")==0) return DenseEnum;
 	      else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
-	      else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
-	      else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"Divergence")==0) return DivergenceEnum;
+	      if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
+	      else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
+	      else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
 	      else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
 	      else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
@@ -1119,10 +1121,10 @@
 	      else if (strcmp(name,"Masscon")==0) return MassconEnum;
 	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
-	      else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
-	      else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+	      if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
+	      else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
+	      else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
 	      else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
 	      else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
@@ -1242,10 +1244,10 @@
 	      else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
 	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
-	      else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
-	      else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
+	      if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
+	      else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
+	      else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
 	      else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
 	      else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
Index: /issm/trunk-jpl/src/m/classes/damage.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/damage.m	(revision 24145)
+++ /issm/trunk-jpl/src/m/classes/damage.m	(revision 24146)
@@ -20,4 +20,5 @@
 		%general parameters for evolution law: 
 		stress_threshold    = 0;
+		stress_ubound       = 0;
 		kappa               = 0;
 		c1                  = 0;
@@ -88,5 +89,5 @@
 			md = checkfield(md,'fieldname','damage.isdamage','values',[1,0]);
 			if self.isdamage,
-				md = checkfield(md,'fieldname','damage.law','numel',[1],'values',[0,1,2]);
+				md = checkfield(md,'fieldname','damage.law','numel',[1],'values',[0,1,2,3]);
 				md = checkfield(md,'fieldname','damage.D','>=',0,'<=',self.max_damage,'size',[md.mesh.numberofvertices 1]);
 				md = checkfield(md,'fieldname','damage.spcdamage','Inf',1,'timeseries',1);
@@ -96,4 +97,5 @@
 				md = checkfield(md,'fieldname','damage.elementinterp','values',{'P1','P2'});
 				md = checkfield(md,'fieldname','damage.stress_threshold','>=',0);
+				md = checkfield(md,'fieldname','damage.stress_ubound','>=',0);
 				md = checkfield(md,'fieldname','damage.kappa','>',1);
 				md = checkfield(md,'fieldname','damage.healing','>=',0);
@@ -131,5 +133,6 @@
 				fielddisplay(self,'maxiter','maximum number of non linear iterations');
 				fielddisplay(self,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}');
-				fielddisplay(self,'stress_threshold','stress threshold for damage initiation [Pa]');
+				fielddisplay(self,'stress_threshold','stress threshold for damage initiation (Pa)');
+				fielddisplay(self,'stress_ubound','stress upper bound for damage healing (Pa), arctan law');
 				fielddisplay(self,'kappa','ductility parameter for stress softening and damage');
 				fielddisplay(self,'c1','damage parameter 1');
@@ -156,4 +159,5 @@
 				WriteData(fid,prefix,'name','md.damage.elementinterp','data',self.elementinterp,'format','String');
 				WriteData(fid,prefix,'object',self,'fieldname','stress_threshold','format','Double');
+				WriteData(fid,prefix,'object',self,'fieldname','stress_ubound','format','Double');
 				WriteData(fid,prefix,'object',self,'fieldname','kappa','format','Double');
 				WriteData(fid,prefix,'object',self,'fieldname','c1','format','Double');
Index: /issm/trunk-jpl/src/m/classes/damage.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/damage.py	(revision 24145)
+++ /issm/trunk-jpl/src/m/classes/damage.py	(revision 24146)
@@ -27,4 +27,5 @@
 		#general parameters for evolution law:
 		self.stress_threshold  = float('NaN')
+                self.stress_ubound     = float('NaN')
 		self.kappa             = float('NaN')
 		self.c1                = float('NaN')
@@ -53,5 +54,6 @@
 			s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations")
 			s+="%s\n" %	fielddisplay(self,"elementinterp","interpolation scheme for finite elements [''P1'',''P2'']")
-			s+="%s\n" % fielddisplay(self,"stress_threshold","stress threshold for damage initiation [Pa]")
+			s+="%s\n" % fielddisplay(self,"stress_threshold","stress threshold for damage initiation (Pa)")
+			s+="%s\n" % fielddisplay(self,"stress_ubound","stress upper bound for damage healing (Pa)")
 			s+="%s\n" % fielddisplay(self,"kappa","ductility parameter for stress softening and damage [>1]")
 			s+="%s\n" % fielddisplay(self,"c1","damage parameter 1 ")
@@ -123,4 +125,5 @@
 			md = checkfield(md,'fieldname','damage.elementinterp','values',['P1','P2'])
 			md = checkfield(md,'fieldname','damage.stress_threshold','>=',0)
+			md = checkfield(md,'fieldname','damage.stress_ubound','>=',0)
 			md = checkfield(md,'fieldname','damage.kappa','>',1)
 			md = checkfield(md,'fieldname','damage.healing','>=',0)
@@ -150,4 +153,5 @@
 			WriteData(fid,prefix,'name','md.damage.elementinterp','data',self.elementinterp,'format','String')
 			WriteData(fid,prefix,'object',self,'fieldname','stress_threshold','format','Double')
+			WriteData(fid,prefix,'object',self,'fieldname','stress_ubound','format','Double')
 			WriteData(fid,prefix,'object',self,'fieldname','kappa','format','Double')
 			WriteData(fid,prefix,'object',self,'fieldname','c1','format','Double')