Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18758)
@@ -19,10 +19,9 @@
 }
 /*}}}*/
 void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-	int  finiteelement;
 
 	/*Finite element type*/
-	finiteelement = P1Enum;
+	int finiteelement = P1Enum;
 
 	/*Update elements: */
 	int counter=0;
@@ -37,7 +36,20 @@
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 	iomodel->FetchDataToInput(elements,VxEnum);
 	iomodel->FetchDataToInput(elements,VyEnum);
-	iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
+
+	/*Get calving parameters*/
+	int calvinglaw;
+	iomodel->Constant(&calvinglaw,CalvingLawEnum);
+	switch(calvinglaw){
+		case DefaultCalvingEnum:
+			iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
+			break;
+		case CalvingLevermannEnum:
+			iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
+			break;
+		default:
+			_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+	}
 }
 /*}}}*/
 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -99,19 +111,21 @@
 	Element* basalelement = element->SpawnBasalElement();
 
 	/*Intermediaries */
-	int  dim, domaintype;
-	bool iscalvingrate;
+	int  stabilization=2;
+	int  dim, domaintype, calvinglaw;
+	bool iscalving;
 	int i, row, col;
 	IssmDouble kappa;
 	IssmDouble Jdet, dt, D_scalar;
 	IssmDouble h,hx,hy,hz;
 	IssmDouble vel;
-//	IssmDouble norm_dlsf;
+	IssmDouble norm_dlsf, calvingrate;
 	IssmDouble* xyz_list = NULL;
 
 	/*Get problem dimension and whether there is calving or not*/
-	basalelement->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
+	basalelement->FindParam(&iscalving,TransientIscalvingEnum);
 	basalelement->FindParam(&domaintype,DomainTypeEnum);
+	basalelement->FindParam(&calvinglaw,CalvingLawEnum);
 	switch(domaintype){
 		case Domain2DverticalEnum:   dim = 1; break;
 		case Domain2DhorizontalEnum: dim = 2; break;
@@ -130,45 +144,63 @@
 	IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
 	IssmDouble*    v        = xNew<IssmDouble>(dim);
 	IssmDouble*    w        = xNew<IssmDouble>(dim);
-	IssmDouble*    c        = xNew<IssmDouble>(dim);
-	//IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
+	IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
+	IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
 
 	/*Retrieve all inputs and parameters*/
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	Input* calvingratex_input=NULL;
-	Input* calvingratey_input=NULL;
+	Input* vx_input           = NULL;
+	Input* vy_input           = NULL;
+	Input* calvingratex_input = NULL;
+	Input* calvingratey_input = NULL;
+	Input* lsf_slopex_input   = NULL;
+	Input* lsf_slopey_input   = NULL;
+	Input* calvingrate_input  = NULL;
+
+	/*Load velocities*/
 	if(domaintype==Domain2DhorizontalEnum){
 		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
 		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-		if(iscalvingrate){
-			calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-			calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-		}
 	}
 	else{
 		if(dim==1){
 			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-			if(iscalvingrate){
-				calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-			}
 		}
 		if(dim==2){
 			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
 			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-			if(iscalvingrate){
-				calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-				calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-			}
 		}
 	}
 
-//	Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-//	Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-	//Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
-	
+	/*Load calving inputs*/
+	if(iscalving){
+		switch(calvinglaw){
+			case DefaultCalvingEnum:
+				lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+				if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
+				calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
+				break;
+			case CalvingLevermannEnum:
+				if(domaintype==Domain2DhorizontalEnum){
+					calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+					calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+				}
+				else{
+					if(dim==1){
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+					}
+					if(dim==2){
+						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+					}
+				}
+				break;
+			default:
+				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+		}
+	}
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -192,34 +224,45 @@
 		GetBprime(Bprime,basalelement,xyz_list,gauss); 
 		vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
 		vy_input->GetInputValue(&v[1],gauss); 
-		if(iscalvingrate){
-			calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
-			calvingratey_input->GetInputValue(&c[1],gauss); 
-			for(i=0;i<dim;i++) w[i]=v[i]-c[i];
+
+		/*Get calving speed*/
+		if(iscalving){
+			switch(calvinglaw){
+				case DefaultCalvingEnum:
+					lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+					if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+					calvingrate_input->GetInputValue(&calvingrate,gauss);
+
+					norm_dlsf=0.;
+					for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
+					norm_dlsf=sqrt(norm_dlsf);
+
+					if(norm_dlsf>1.e-10)
+					 for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
+					else
+					 for(i=0;i<dim;i++) c[i]=0.;
+					break;
+				case CalvingLevermannEnum:
+					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
+					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
+					break;
+				default:
+					_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+			}
 		}
-		else{
-			for(i=0;i<dim;i++) w[i]=v[i];
-		}
-		//lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-		//lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-		//calvingrate_input->GetInputValue(&calvingrate,gauss);
 
-		//norm_dlsf=0.;
-		//for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-		//norm_dlsf=sqrt(norm_dlsf);
+		/*Levelset speed is ice velocity - calving rate*/
+		for(i=0;i<dim;i++) w[i]=v[i]-c[i];
 
-		//if(norm_dlsf>1.e-10)
-		//	for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
-		//else
-		//for(i=0;i<dim;i++) c[i]=0.;
-		
-
-		for(row=0;row<dim;row++)
-			for(col=0;col<dim;col++)
+		/*Compute D*/
+		for(row=0;row<dim;row++){
+			for(col=0;col<dim;col++){
 				if(row==col)
-					D[row*dim+col]=D_scalar*w[row];
+				 D[row*dim+col]=D_scalar*w[row];
 				else
-				   D[row*dim+col]=0.;
+				 D[row*dim+col]=0.;
+			}
+		}
 
 		TripleMultiply(B,dim,numnodes,1,
 					D,dim,dim,0,
@@ -227,9 +270,8 @@
 					&Ke->values[0],1);
 
 		/* Stabilization */
-		int stabilization=2;
 		vel=0.;
-		for(i=0;i<dim;i++) vel+=w[i]*w[i];
+		for(i=0;i<dim;i++) vel+=v[i]*v[i];
 		vel=sqrt(vel)+1.e-14;
 		switch(stabilization){
 			case 0:
@@ -255,10 +297,10 @@
 			case 2:
 				/* Streamline Upwinding */
 				basalelement->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
 				for(row=0;row<dim;row++) 
 					for(col=0;col<dim;col++) 
-						D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
+						D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
 
 				TripleMultiply(Bprime,dim,numnodes,1,
 							D,dim,dim,0,
@@ -279,7 +321,7 @@
 	xDelete<IssmDouble>(v);
 	xDelete<IssmDouble>(w);
 	xDelete<IssmDouble>(c);
-	//xDelete<IssmDouble>(dlsf);
+	xDelete<IssmDouble>(dlsf);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18757)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18758)
@@ -208,6 +208,12 @@
 	DamageEnum,
 	NewDamageEnum,
 	StressIntensityFactorEnum,
+	CalvingLawEnum,
+	CalvingCalvingrateEnum,
+	CalvingLevermannEnum,
+	DefaultCalvingEnum,
+	CalvingRequestedOutputsEnum,
+	CalvinglevermannCoeffEnum,
 	CalvingratexEnum,
 	CalvingrateyEnum,
 	CalvingratexAverageEnum,
@@ -249,13 +255,10 @@
 	Domain3DEnum,
 	MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
 	MasstransportHydrostaticAdjustmentEnum,
-	MasstransportIscalvingrateEnum,
-	MasstransportLevermannCalvingCoeffEnum,
 	MasstransportIsfreesurfaceEnum,
 	MasstransportMinThicknessEnum,
 	MasstransportPenaltyFactorEnum,
 	MasstransportSpcthicknessEnum,
-	MasstransportCalvingrateEnum,
 	MasstransportStabilizationEnum,
 	MasstransportVertexPairingEnum,
 	MasstransportNumRequestedOutputsEnum,
@@ -313,6 +316,7 @@
 	TransientIsgiaEnum,
 	TransientIsdamageevolutionEnum,
 	TransientIshydrologyEnum,
+	TransientIscalvingEnum,
 	TransientNumRequestedOutputsEnum,
 	TransientRequestedOutputsEnum,
 	PotentialEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18758)
@@ -216,6 +216,12 @@
 		case DamageEnum : return "Damage";
 		case NewDamageEnum : return "NewDamage";
 		case StressIntensityFactorEnum : return "StressIntensityFactor";
+		case CalvingLawEnum : return "CalvingLaw";
+		case CalvingCalvingrateEnum : return "CalvingCalvingrate";
+		case CalvingLevermannEnum : return "CalvingLevermann";
+		case DefaultCalvingEnum : return "DefaultCalving";
+		case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
+		case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
 		case CalvingratexEnum : return "Calvingratex";
 		case CalvingrateyEnum : return "Calvingratey";
 		case CalvingratexAverageEnum : return "CalvingratexAverage";
@@ -257,13 +263,10 @@
 		case Domain3DEnum : return "Domain3D";
 		case MiscellaneousNameEnum : return "MiscellaneousName";
 		case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
-		case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate";
-		case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff";
 		case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
 		case MasstransportMinThicknessEnum : return "MasstransportMinThickness";
 		case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
 		case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness";
-		case MasstransportCalvingrateEnum : return "MasstransportCalvingrate";
 		case MasstransportStabilizationEnum : return "MasstransportStabilization";
 		case MasstransportVertexPairingEnum : return "MasstransportVertexPairing";
 		case MasstransportNumRequestedOutputsEnum : return "MasstransportNumRequestedOutputs";
@@ -321,6 +324,7 @@
 		case TransientIsgiaEnum : return "TransientIsgia";
 		case TransientIsdamageevolutionEnum : return "TransientIsdamageevolution";
 		case TransientIshydrologyEnum : return "TransientIshydrology";
+		case TransientIscalvingEnum : return "TransientIscalving";
 		case TransientNumRequestedOutputsEnum : return "TransientNumRequestedOutputs";
 		case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
 		case PotentialEnum : return "Potential";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18758)
@@ -219,6 +219,12 @@
 	      else if (strcmp(name,"Damage")==0) return DamageEnum;
 	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
 	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
+	      else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
+	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
+	      else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
+	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
+	      else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
+	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
 	      else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
 	      else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
 	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
@@ -253,23 +259,20 @@
 	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
 	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
 	      else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
-	      else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
 	      else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
 	      else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
 	      else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
 	      else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
 	      else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
-	      else if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum;
-	      else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum;
+	      else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
 	      else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
 	      else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
 	      else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
 	      else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
-	      else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
 	      else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
 	      else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
 	      else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
@@ -327,6 +330,7 @@
 	      else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
 	      else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
 	      else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
+	      else if (strcmp(name,"TransientIscalving")==0) return TransientIscalvingEnum;
 	      else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
 	      else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
 	      else if (strcmp(name,"Potential")==0) return PotentialEnum;
@@ -378,14 +382,14 @@
 	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
-	      else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
 	      else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
+	      else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
 	      else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
 	      else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
@@ -501,14 +505,14 @@
 	      else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
 	      else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
 	      else if (strcmp(name,"StringParam")==0) return StringParamEnum;
-	      else if (strcmp(name,"Seg")==0) return SegEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"Seg")==0) return SegEnum;
 	      else if (strcmp(name,"SegInput")==0) return SegInputEnum;
 	      else if (strcmp(name,"Tria")==0) return TriaEnum;
 	      else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"Tetra")==0) return TetraEnum;
+	      else if (strcmp(name,"Tetra")==0) return TetraEnum;
 	      else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
 	      else if (strcmp(name,"Penta")==0) return PentaEnum;
 	      else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
@@ -624,14 +628,14 @@
 	      else if (strcmp(name,"P2")==0) return P2Enum;
 	      else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
 	      else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
-	      else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
 	      else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
 	      else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
 	      else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"P1P1")==0) return P1P1Enum;
+	      else if (strcmp(name,"P1P1")==0) return P1P1Enum;
 	      else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
 	      else if (strcmp(name,"MINI")==0) return MINIEnum;
 	      else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
@@ -747,14 +751,14 @@
 	      else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
 	      else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
 	      else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
-	      else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
 	      else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
+	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
 	      else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
 	      else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
 	      else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18758)
@@ -66,6 +66,7 @@
 	parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(CalvingLawEnum));
 	if(solution_type==SeaiceSolutionEnum){
 	}
 	else{
@@ -82,8 +83,7 @@
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
-		parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum));
-		parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIscalvingEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
 
Index: ../trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18758)
@@ -20,11 +20,11 @@
 
 	/*parameters: */
 	IssmDouble starttime,finaltime,dt,yts;
-	bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate;
+	bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
 	bool       save_results,dakota_analysis;
 	bool       time_adapt=false;
 	int        output_frequency;
-	int        domaintype,groundingline_migration;
+	int        domaintype,groundingline_migration,calvinglaw;
 	int        numoutputs         = 0;
 	Analysis  *analysis = NULL;
 	char**     requested_outputs = NULL;
@@ -52,7 +52,8 @@
 	femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
 	femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
 	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
-	femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
+	femmodel->parameters->FindParam(&iscalving,TransientIscalvingEnum);
+	femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
 	if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
 	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
 	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
@@ -101,7 +102,7 @@
 		}
 
 		if(islevelset){
-			if(iscalvingrate){
+			if(iscalving && calvinglaw==CalvingLevermannEnum){
 				if(VerboseSolution()) _printf0_("   computing calving rate\n");
 				femmodel->StrainRateparallelx();
 				femmodel->StrainRateperpendicularx();
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18758)
@@ -438,7 +438,7 @@
 	Input* vy_input=inputs->GetInput(VyEnum);																		_assert_(vy_input);
 	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);								_assert_(strainparallel_input);
 	Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);					_assert_(strainperpendicular_input);
-	this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
+	Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
 
 	/* Start looping on the number of vertices: */
 	gauss=new GaussTria();
@@ -451,6 +451,7 @@
 		vel=vx*vx+vy*vy;
 		strainparallel_input->GetInputValue(&strainparallel,gauss);
 		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
+		levermanncoeff_input->GetInputValue(&propcoeff,gauss);
 
 		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
 		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
@@ -464,7 +465,7 @@
 	/*Add input*/
 	this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
 	this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
 
 	/*Clean up and return*/
 	delete gauss;
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18758)
@@ -366,8 +366,8 @@
 	Input* vx_input=inputs->GetInput(VxEnum);																		_assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);																		_assert_(vy_input);
 	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);								_assert_(strainparallel_input);
-	Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);             _assert_(strainperpendicular_input);
-	this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
+	Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
+	Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
 
 	/* Start looping on the number of vertices: */
 	gauss=new GaussPenta();
@@ -380,6 +380,7 @@
 		vel=vx*vx+vy*vy;
 		strainparallel_input->GetInputValue(&strainparallel,gauss);
 		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
+		levermanncoeff_input->GetInputValue(&propcoeff,gauss);
 
 		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
 		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;	
@@ -393,7 +394,7 @@
 	/*Add input*/
 	this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
 	this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
 
 	/*Clean up and return*/
 	delete gauss;
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18757)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18758)
@@ -1091,7 +1091,7 @@
 				name==LevelsetfunctionSlopeXEnum ||
 				name==LevelsetfunctionSlopeYEnum ||
 				name==LevelsetfunctionPicardEnum ||
-				name==MasstransportCalvingrateEnum ||
+				//name==CalvingCalvingrateEnum ||
 				name==GradientEnum ||
 				name==OldGradientEnum  ||
 				name==ConvergedEnum || 
@@ -1211,7 +1211,7 @@
 				break;
 			case CalvingratexEnum:
 			case CalvingrateyEnum:
-			case MasstransportCalvingrateEnum:
+			case CalvingCalvingrateEnum:
 				this->StrainRateparallel();
 				this->StrainRateperpendicular();
 				this->CalvingRateLevermann();
