Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27013)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27014)
@@ -646,10 +646,56 @@
 	int  calvinglaw;
 	femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
+	IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum);
+	IssmDouble dt      = femmodel->parameters->FindParam(TimesteppingTimeStepEnum);
+
+
+   /*Get current distance to terminus*/
+   InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
+   femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
+
 
 	if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum){
 
 		/*Intermediaries*/
-		IssmDouble thickness,bed,sealevel;
+		IssmDouble thickness,bed,sealevel,distance;
 		IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
+
+		/*Loop over all elements of this partition*/
+		for(Object* & object : femmodel->elements->objects){
+			Element* element  = xDynamicCast<Element*>(object);
+
+			int      numnodes = element->GetNumberOfNodes();
+			Gauss*   gauss    = element->NewGauss();
+			Input*   H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
+			Input*   b_input  = element->GetInput(BedEnum); _assert_(b_input);
+			Input*   sl_input = element->GetInput(SealevelEnum); _assert_(sl_input);
+			Input*   dis_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
+
+			/*Potentially constrain nodes of this element*/
+			for(int in=0;in<numnodes;in++){
+				gauss->GaussNode(element->GetElementType(),in);
+				Node* node=element->GetNode(in);
+				if(!node->IsActive()) continue;
+
+				H_input->GetInputValue(&thickness,gauss);
+				b_input->GetInputValue(&bed,gauss);
+				sl_input->GetInputValue(&sealevel,gauss);
+				dis_input->GetInputValue(&distance,gauss);
+				if(thickness<min_thickness && bed<sealevel && distance<mig_max*dt){
+					node->ApplyConstraint(0,+1.);
+				}
+				else {
+					/* no ice, set no spc */
+					node->DofInFSet(0);
+				}
+			}
+			delete gauss;
+		}
+	}
+	else if(calvinglaw==CalvingParameterizationEnum){
+
+      /*Intermediaries*/
+      IssmDouble thickness,bed,sealevel,distance;
+      IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
 
 		/*Loop over all elements of this partition*/
@@ -662,4 +708,5 @@
 			Input*   b_input = element->GetInput(BedEnum); _assert_(b_input);
 			Input*   sl_input = element->GetInput(SealevelEnum); _assert_(sl_input);
+			Input*   dis_input           = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
 
 			/*Potentially constrain nodes of this element*/
@@ -672,5 +719,6 @@
 				b_input->GetInputValue(&bed,gauss);
 				sl_input->GetInputValue(&sealevel,gauss);
-				if(thickness<min_thickness && bed<sealevel){
+				dis_input->GetInputValue(&distance,gauss);
+				if(thickness<min_thickness && bed<sealevel && distance<mig_max*dt){
 					node->ApplyConstraint(0,+1.);
 				}
@@ -683,9 +731,8 @@
 		}
 	}
-	else if(calvinglaw==CalvingParameterizationEnum){
-
-      /*Intermediaries*/
-      IssmDouble thickness,bed,sealevel;
-      IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
+   else if(calvinglaw==CalvingHabEnum){
+
+		/*Intermediaries*/
+		IssmDouble  thickness,water_depth,distance,hab_fraction;
 
 		/*Loop over all elements of this partition*/
@@ -693,9 +740,13 @@
 			Element* element  = xDynamicCast<Element*>(object);
 
-			int      numnodes = element->GetNumberOfNodes();
-			Gauss*   gauss    = element->NewGauss();
-			Input*   H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
-			Input*   b_input = element->GetInput(BedEnum); _assert_(b_input);
-			Input*   sl_input = element->GetInput(SealevelEnum); _assert_(sl_input);
+			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
+
+			int      numnodes           = element->GetNumberOfNodes();
+			Gauss*   gauss              = element->NewGauss();
+			Input*   H_input            = element->GetInput(ThicknessEnum); _assert_(H_input);
+			Input*   bed_input          = element->GetInput(BedEnum); _assert_(bed_input);
+			Input*   hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input);
+			Input*   dis_input           = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
 
 			/*Potentially constrain nodes of this element*/
@@ -706,52 +757,9 @@
 
 				H_input->GetInputValue(&thickness,gauss);
-				b_input->GetInputValue(&bed,gauss);
-				sl_input->GetInputValue(&sealevel,gauss);
-				if(thickness<min_thickness && bed<sealevel){
-					node->ApplyConstraint(0,+1.);
-				}
-				else {
-					/* no ice, set no spc */
-					node->DofInFSet(0);
-				}
-			}
-			delete gauss;
-		}
-	}
-   else if(calvinglaw==CalvingHabEnum){
-
-		/*Intermediaries*/
-		IssmDouble  thickness,water_depth,levelset,hab_fraction;
-
-		/*Get the fraction of the flotation thickness at the terminus*/
-		InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
-		femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
-
-		/*Loop over all elements of this partition*/
-		for(Object* & object : femmodel->elements->objects){
-			Element* element  = xDynamicCast<Element*>(object);
-
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
-			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
-
-			int      numnodes           = element->GetNumberOfNodes();
-			Gauss*   gauss              = element->NewGauss();
-			Input*   H_input            = element->GetInput(ThicknessEnum); _assert_(H_input);
-			Input*   bed_input          = element->GetInput(BedEnum); _assert_(bed_input);
-			Input*   hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input);
-			Input*   ls_input           = element->GetInput(DistanceToCalvingfrontEnum); _assert_(ls_input);
-
-			/*Potentially constrain nodes of this element*/
-			for(int in=0;in<numnodes;in++){
-				gauss->GaussNode(element->GetElementType(),in);
-				Node* node=element->GetNode(in);
-				if(!node->IsActive()) continue;
-
-				H_input->GetInputValue(&thickness,gauss);
 				bed_input->GetInputValue(&water_depth,gauss);
-				ls_input->GetInputValue(&levelset,gauss);
+				dis_input->GetInputValue(&distance,gauss);
 				hab_fraction_input->GetInputValue(&hab_fraction,gauss);
 
-				if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300. && levelset<0.){
+				if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && distance<mig_max*dt){
 					node->ApplyConstraint(0,+1.);
 				}
@@ -780,6 +788,5 @@
 
 		IssmDouble crevasse_threshold = femmodel->parameters->FindParam(CalvingCrevasseThresholdEnum);
-		IssmDouble mig_max            = femmodel->parameters->FindParam(MigrationMaxEnum);
-		IssmDouble dt                 = femmodel->parameters->FindParam(TimesteppingTimeStepEnum);
+
 
 		for(Object* & object : femmodel->elements->objects){
@@ -856,4 +863,5 @@
 						surface_input->GetInputValue(&surface,gauss);
 
+                  /*FIXME: not sure about levelset<0. && levelset>-mig_max*dt! SHould maybe be distance<mig_max*dt*/
                   if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0. && levelset<0. && levelset>-mig_max*dt && constraint_nodes[node->Lid()]==0.){
 							local_nflipped++;