Index: issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 21626)
@@ -61,5 +61,4 @@
 			break;
 		case CalvingDevEnum:
-			iomodel->FetchDataToInput(elements,"md.calving.coeff",CalvingdevCoeffEnum);
 			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			break;
@@ -75,4 +74,5 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.levelset.stabilization",LevelsetStabilizationEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.levelset.calving_max",CalvingMaxEnum));
 	int  calvinglaw;
 	iomodel->FindConstant(&calvinglaw,"md.calving.law");
@@ -80,5 +80,8 @@
 		case DefaultCalvingEnum:
 		case CalvingLevermannEnum:
+			break;
 		case CalvingDevEnum:
+			parameters->AddObject(iomodel->CopyConstantObject("md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum));
+			parameters->AddObject(iomodel->CopyConstantObject("md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum));
 			break;
 		case CalvingMinthicknessEnum:
@@ -138,5 +141,6 @@
 	IssmDouble h,hx,hy,hz;
 	IssmDouble vel;
-	IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate;
+	IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
+	IssmDouble calvingmax;
 	IssmDouble* xyz_list = NULL;
 
@@ -152,4 +156,6 @@
 	}
 
+	/*Calving threshold*/
+	
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes    = basalelement->GetNumberOfNodes();
@@ -170,4 +176,5 @@
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	basalelement->FindParam(&calvingmax,CalvingMaxEnum);
 	Input* vx_input           = NULL;
 	Input* vy_input           = NULL;
@@ -178,4 +185,5 @@
 	Input* calvingrate_input  = NULL;
 	Input* meltingrate_input  = NULL;
+	Input* gr_input           = NULL;
 
 	/*Load velocities*/
@@ -187,8 +195,10 @@
 			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
 			vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
+			gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
 			break;
 		case Domain3DEnum:
 			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
 			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
+			gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
 			break;
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
@@ -253,4 +263,5 @@
 		vx_input->GetInputValue(&v[0],gauss);
 		vy_input->GetInputValue(&v[1],gauss); 
+		gr_input->GetInputValue(&groundedice,gauss);
 
 		/*Get calving speed*/
@@ -262,4 +273,9 @@
 				calvingrate_input->GetInputValue(&calvingrate,gauss);
 				meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+				/*Limit calving rate to c <= v + 3 km/yr */
+				vel=sqrt(v[0]*v[0] + v[1]*v[1]);
+				if(calvingrate>calvingmax+vel) calvingrate = vel+calvingmax;
+				if(groundedice<0) meltingrate = 0;
 
 				norm_dlsf=0.;
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 21626)
@@ -186,5 +186,6 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 	IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-	IssmDouble  sigma_vm,sigma_max,epse_2,groundedice;
+	IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
+	IssmDouble  epse_2,groundedice;
 
 	/* Get node coordinates and dof list: */
@@ -202,4 +203,6 @@
 	IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
 	IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
+	this->parameters->FindParam(&sigma_max_floating,CalvingStressThresholdFloatingiceEnum);
+	this->parameters->FindParam(&sigma_max_grounded,CalvingStressThresholdGroundediceEnum);
 
 	/* Start looping on the number of vertices: */
@@ -229,7 +232,10 @@
 		epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
 		sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
-		sigma_max = 1000.e+3;
-
-		if(groundedice<0) sigma_max=200.e+3;
+
+		/*Tensile stress threshold*/
+		if(groundedice<0)
+		 sigma_max = sigma_max_floating;
+		else
+		 sigma_max = sigma_max_grounded;
 
 		/*Assign values*/
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21626)
@@ -216,5 +216,6 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 	IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-	IssmDouble  sigma_vm,sigma_max,epse_2,groundedice;
+	IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
+	IssmDouble  epse_2,groundedice;
 
 	/* Get node coordinates and dof list: */
@@ -227,4 +228,6 @@
 	IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
 	IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
+	this->parameters->FindParam(&sigma_max_floating,CalvingStressThresholdFloatingiceEnum);
+	this->parameters->FindParam(&sigma_max_grounded,CalvingStressThresholdGroundediceEnum);
 
 	/* Start looping on the number of vertices: */
@@ -254,17 +257,15 @@
 		epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
 		sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
-		//sigma_max = 125.e+3;
-		sigma_max = 350.e+3;
-		sigma_max = 450.e+3;
-		sigma_max = 800.e+3; //too much
-		//sigma_max = 700.e+3;
-		//sigma_max = 670.e+3;
-		//sigma_max = 550.e+3;
-		sigma_max = 750.e+3; //too high
-		sigma_max = 850.e+3; //too low
-		sigma_max = 800.e+3; //IUGG previous test
-		sigma_max = 1000.e+3; //850 seems small
-
-		if(groundedice<0) sigma_max=200.e+3;
+
+		/*OLD (keep for a little bit)*/
+		//sigma_max = 800.e+3; //IUGG previous test
+		//sigma_max = 1000.e+3; //GRL
+		//if(groundedice<0) sigma_max=150.e+3;
+
+		/*Tensile stress threshold*/
+		if(groundedice<0)
+		 sigma_max = sigma_max_floating;
+		else
+		 sigma_max = sigma_max_grounded;
 
 		/*Assign values*/
@@ -294,5 +295,4 @@
 	IssmDouble  calvingratey[NUMVERTICES];
 	IssmDouble  calvingrate[NUMVERTICES];
-
 
 	/* Get node coordinates and dof list: */
Index: issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp	(revision 21626)
@@ -33,5 +33,5 @@
 	}
 	if(domaintype==Domain2DverticalEnum){
-	      femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToExtrudeEnum);
+		femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToExtrudeEnum);
 		extrudefrombase_core(femmodel);
 	}
Index: issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21625)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21626)
@@ -259,4 +259,7 @@
 	CalvingratexAverageEnum,
 	CalvingrateyAverageEnum,
+	CalvingStressThresholdGroundediceEnum,
+	CalvingStressThresholdFloatingiceEnum,
+	CalvingMaxEnum,
 	StrainRateparallelEnum,
 	StrainRateperpendicularEnum,
Index: issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21626)
@@ -265,4 +265,7 @@
 		case CalvingratexAverageEnum : return "CalvingratexAverage";
 		case CalvingrateyAverageEnum : return "CalvingrateyAverage";
+		case CalvingStressThresholdGroundediceEnum : return "CalvingStressThresholdGroundedice";
+		case CalvingStressThresholdFloatingiceEnum : return "CalvingStressThresholdFloatingice";
+		case CalvingMaxEnum : return "CalvingMax";
 		case StrainRateparallelEnum : return "StrainRateparallel";
 		case StrainRateperpendicularEnum : return "StrainRateperpendicular";
Index: issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21625)
+++ issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21626)
@@ -271,4 +271,7 @@
 	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
 	      else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
+	      else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
+	      else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
+	      else if (strcmp(name,"CalvingMax")==0) return CalvingMaxEnum;
 	      else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
 	      else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
@@ -380,11 +383,11 @@
 	      else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
 	      else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
-	      else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
-	      else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
-	      else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
+	      if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
+	      else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
+	      else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
+	      else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
 	      else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
 	      else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
@@ -503,11 +506,11 @@
 	      else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
 	      else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
-	      else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
-	      else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
-	      else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
+	      if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
+	      else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
+	      else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
+	      else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
 	      else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
 	      else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
@@ -626,11 +629,11 @@
 	      else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
 	      else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
-	      else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
-	      else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
-	      else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
+	      if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
+	      else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
+	      else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
+	      else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
 	      else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
 	      else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
@@ -749,11 +752,11 @@
 	      else if (strcmp(name,"Scaled")==0) return ScaledEnum;
 	      else if (strcmp(name,"Separate")==0) return SeparateEnum;
-	      else if (strcmp(name,"Sset")==0) return SsetEnum;
-	      else if (strcmp(name,"Dense")==0) return DenseEnum;
-	      else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
+	      if (strcmp(name,"Sset")==0) return SsetEnum;
+	      else if (strcmp(name,"Dense")==0) return DenseEnum;
+	      else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
+	      else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
 	      else if (strcmp(name,"Seq")==0) return SeqEnum;
 	      else if (strcmp(name,"Mpi")==0) return MpiEnum;
@@ -872,11 +875,11 @@
 	      else if (strcmp(name,"Penta")==0) return PentaEnum;
 	      else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
-	      else if (strcmp(name,"Vertex")==0) return VertexEnum;
-	      else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
-	      else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"Option")==0) return OptionEnum;
+	      if (strcmp(name,"Vertex")==0) return VertexEnum;
+	      else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
+	      else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
+	      else if (strcmp(name,"Option")==0) return OptionEnum;
 	      else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
 	      else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
@@ -995,11 +998,11 @@
 	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
 	      else if (strcmp(name,"Closed")==0) return ClosedEnum;
-	      else if (strcmp(name,"Free")==0) return FreeEnum;
-	      else if (strcmp(name,"Open")==0) return OpenEnum;
-	      else if (strcmp(name,"Air")==0) return AirEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"Ice")==0) return IceEnum;
+	      if (strcmp(name,"Free")==0) return FreeEnum;
+	      else if (strcmp(name,"Open")==0) return OpenEnum;
+	      else if (strcmp(name,"Air")==0) return AirEnum;
+	      else if (strcmp(name,"Ice")==0) return IceEnum;
 	      else if (strcmp(name,"Melange")==0) return MelangeEnum;
 	      else if (strcmp(name,"Water")==0) return WaterEnum;
Index: issm/trunk-jpl/src/m/classes/calvingdev.m
===================================================================
--- issm/trunk-jpl/src/m/classes/calvingdev.m	(revision 21625)
+++ issm/trunk-jpl/src/m/classes/calvingdev.m	(revision 21626)
@@ -6,5 +6,6 @@
 classdef calvingdev
 	properties (SetAccess=public) 
-		coeff         = NaN;
+		stress_threshold_groundedice = 0.;
+		stress_threshold_floatingice = 0.;
 		meltingrate   = NaN;
 	end
@@ -29,5 +30,4 @@
 		end % }}}
 		function self = extrude(self,md) % {{{
-			self.coeff=project3d(md,'vector',self.coeff,'type','node');
 			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
@@ -35,5 +35,9 @@
 
 			%Proportionality coefficient in Pi model
-			self.coeff=2e13;
+			%self.coeff=2e13;
+
+			%Default sigma max
+			self.stress_threshold_groundedice = 1e6;
+			self.stress_threshold_floatingice = 150e3;
 		end % }}}
 		function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -41,10 +45,12 @@
 			if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end
 
-			md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'numel',1,'nan',1,'Inf',1);
+			md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'numel',1,'nan',1,'Inf',1);
 			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Calving Pi parameters:'));
-			fielddisplay(self,'coeff','proportionality coefficient in Pi model');
+			fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]');
+			fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]');
 			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 
@@ -53,6 +59,7 @@
 			yts=md.constants.yts;
 			WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer');
-			WriteData(fid,prefix,'object',self,'fieldname','coeff','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
+			WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
Index: issm/trunk-jpl/src/m/classes/calvinglevermann.m
===================================================================
--- issm/trunk-jpl/src/m/classes/calvinglevermann.m	(revision 21625)
+++ issm/trunk-jpl/src/m/classes/calvinglevermann.m	(revision 21626)
@@ -42,5 +42,5 @@
 
 			md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
-			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>=',0);
+			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 		end % }}}
 		function disp(self) % {{{
@@ -54,5 +54,5 @@
 			WriteData(fid,prefix,'name','md.calving.law','data',3,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','coeff','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
+			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
Index: issm/trunk-jpl/src/m/classes/calvingminthickness.m
===================================================================
--- issm/trunk-jpl/src/m/classes/calvingminthickness.m	(revision 21625)
+++ issm/trunk-jpl/src/m/classes/calvingminthickness.m	(revision 21626)
@@ -53,5 +53,5 @@
 			WriteData(fid,prefix,'name','md.calving.law','data',4,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double');
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
+			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
Index: issm/trunk-jpl/src/m/classes/levelset.m
===================================================================
--- issm/trunk-jpl/src/m/classes/levelset.m	(revision 21625)
+++ issm/trunk-jpl/src/m/classes/levelset.m	(revision 21626)
@@ -9,4 +9,5 @@
 		spclevelset			= NaN;
 		reinit_frequency	= 5;
+		calving_max       = 0.;
 	end
 	methods
@@ -36,6 +37,7 @@
 
 			%stabilization = 2 by default
-			self.stabilization = 2;
+			self.stabilization    = 2;
 			self.reinit_frequency = 5;
+			self.calving_max      = 3000.;
 
 		end % }}}
@@ -46,4 +48,5 @@
 			md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1);
 			md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2]);
+			md = checkfield(md,'fieldname','levelset.calving_max','numel',1,'NaN',1,'Inf',1,'>',0);
 		end % }}}
 		function disp(self) % {{{
@@ -52,9 +55,14 @@
 			fielddisplay(self,'spclevelset','Levelset constraints (NaN means no constraint)');
 			fielddisplay(self,'reinit_frequency','Amount of time steps after which the levelset function in re-initialized');
+			fielddisplay(self,'calving_max','maximum allowed calving rate (m/a)');
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
+
+			yts=md.constants.yts;
+
 			WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','spclevelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','reinit_frequency','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','calving_max','format','Double','scale',1./yts);
 		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
@@ -63,4 +71,5 @@
 			writejs1Darray(fid,[modelname '.levelset.spclevelset'],self.spclevelset);
 			writejs1Darray(fid,[modelname '.levelset.reinit_frequency'],self.reinit_frequency);
+			writejsdouble(fid,[modelname '.levelset.calving_max'],self.calving_max);
 
 		end % }}}
Index: issm/trunk-jpl/src/m/classes/levelset.py
===================================================================
--- issm/trunk-jpl/src/m/classes/levelset.py	(revision 21625)
+++ issm/trunk-jpl/src/m/classes/levelset.py	(revision 21626)
@@ -14,7 +14,8 @@
 	def __init__(self): # {{{
 
-		self.stabilization = 0
-		self.spclevelset   = float('NaN')
+		self.stabilization    = 0
+		self.spclevelset      = float('NaN')
 		self.reinit_frequency = 0
+		self.calving_max      = 0.
 
 		#set defaults
@@ -27,4 +28,5 @@
 		string="%s\n%s"%(string,fielddisplay(self,'spclevelset','levelset constraints (NaN means no constraint)'))
 		string="%s\n%s"%(string,fielddisplay(self,'reinit_frequency','Amount of time steps after which the levelset function in re-initialized'))
+		string="%s\n%s"%(string,fielddisplay(self,'calving_max','maximum allowed calving rate (m/a)'))
 
 		return string
@@ -39,4 +41,5 @@
 		self.stabilization = 2
 		self.reinit_frequency = 5
+		self.calving_max      = 3000
 
 		return self
@@ -50,4 +53,5 @@
 		md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1)
 		md = checkfield(md,'fieldname','levelset.stabilization','values',[0,1,2]);
+		md = checkfield(md,'fieldname','levelset.calving_max','numel',1,'NaN',1,'Inf',1,'>',0);
 
 		return md
@@ -55,6 +59,9 @@
 	def marshall(self,prefix,md,fid):    # {{{
 
+		yts=md.constants.yts;
+
 		WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer');
 		WriteData(fid,prefix,'object',self,'fieldname','spclevelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 		WriteData(fid,prefix,'object',self,'fieldname','reinit_frequency','format','Integer');
+		WriteData(fid,prefix,'object',self,'fieldname','calving_max','format','Double','scale',1./yts);
 	# }}}
