Index: /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp	(revision 27491)
@@ -10,4 +10,5 @@
 
 #define FINITEELEMENT P1Enum
+#define EPS 1e-14
 
 /*Model processing*/
@@ -59,5 +60,5 @@
 
 	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
-	::CreateNodes(nodes,iomodel,DebrisAnalysisEnum,FINITEELEMENT);
+	::CreateNodes(nodes,iomodel,DebrisAnalysisEnum,FINITEELEMENT,isamr);
 	iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 
@@ -127,85 +128,9 @@
 void           DebrisAnalysis::Core(FemModel* femmodel){/*{{{*/
 
-
-	Element* element= NULL;
-	for(Object* & object : femmodel->elements->objects){
-		element=xDynamicCast<Element*>(object);
-
-
-		int numvertices = element->GetNumberOfNodes();
-
-		IssmDouble* vx = xNew<IssmDouble>(numvertices);
-		IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
-		IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
-		IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
-		IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
-
-		element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
-		element->GetInputListOnVertices(&vx[0],VxEnum);
-		element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
-		element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
-		element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
-
-		IssmDouble slope,rad2deg=180./M_PI; //=57.2958
-		IssmDouble vslipx,rhod=1900.;
-		IssmDouble gravity=element->FindParam(ConstantsGEnum);
-		IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
-		IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
-
-		int step;
-		IssmDouble dt, maxv;
-		IssmDouble yts=31536000.;
-		femmodel->parameters->FindParam(&step,StepEnum);
-		femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-
-		bool isminthicknessinelement=false;
-		for(int i=0;i<numvertices;i++){
-			if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
-		}
-		if(isminthicknessinelement){
-			//do nothing
-			for(int i=0;i<numvertices;i++){
-				if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;                         
-			}
-		}else{
-			for(int i=0;i<numvertices;i++){
-				//if(onsurface[i]>.5){
-				slope=fabs(slopex[i]);
-				if((atan(slope)*rad2deg)>25.){
-					//if(debristhickness[i]>0.01){
-					vslipx=slope_threshold/yts;
-					//maxv=10.0/2./dt;
-					//vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts;
-					vx[i]=vx[i]+vslipx;
-					//debristhickness[i]=debristhickness[i];
-					//if(vx[i]>maxv) vx[i]=maxv;
-					//}
-				} 
-				//}
-			}
-		}
-		//if(step%100==0)   
-		element->AddInput(VxDebrisEnum,vx,P1Enum);
-		//element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
-
-		/* Free resources */
-		xDelete<IssmDouble>(debristhickness);
-		xDelete<IssmDouble>(icethickness);
-		xDelete<IssmDouble>(vx);
-		xDelete<IssmDouble>(slopex);
-		xDelete<IssmDouble>(onsurface);
-	}
-
-	//if(step%7==0) 
 	//PreProcessing(femmodel);
 	//femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
 	//extrudefromtop_core(femmodel);
-
 	femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);        
 	solutionsequence_linear(femmodel);
-
-	int step;
-	femmodel->parameters->FindParam(&step,StepEnum);
-	//if(step%7==0) PreProcessing(femmodel);
 	PostProcessing(femmodel);
 
@@ -230,5 +155,6 @@
 	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble xi,tau;
+	IssmDouble yts=31536000.;
+	IssmDouble tau;
 	IssmDouble dvx[2],dvy[2];
 	Element*    topelement = NULL;
@@ -271,5 +197,5 @@
 	Input* vy_input=NULL;
 	if(dim>1){vy_input = topelement->GetInput(VyDebrisEnum); _assert_(vy_input);}
-	h = topelement->CharacteristicLength();
+	h=topelement->CharacteristicLength();
 
 	/* Start  looping on the number of gaussian points: */
@@ -302,6 +228,5 @@
 				}
 			}
-		}
-		else{
+		}else{
 			dvxdx=dvx[0];
 			for(int i=0;i<numnodes;i++){
@@ -313,101 +238,51 @@
 		}
 
-		/*******************************************************************/
-		/* Diffusion */
-		bool isdisplacement=false;
-		int step;
-		topelement->FindParam(&step,StepEnum);
-		IssmDouble slope_threshold;
-		topelement->FindParam(&slope_threshold,DebrisRemovalSlopeThresholdEnum);
-		IssmDouble kappa,f,smb,debristhickness,slopex;
-		IssmDouble Diff,fraction,M=1,C;
-		IssmDouble rad2deg=180./M_PI;
-		Diff=3.2/3e7;
-		Input* slopex_input=topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
-		Input* smb_input=topelement->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
-		Input* debristhickness_input=topelement->GetInput(DebrisThicknessEnum); _assert_(debristhickness_input);
-
-		if(isdisplacement){
-
-			slopex_input->GetInputValue(&slopex, gauss);
-			smb_input->GetInputValue(&smb, gauss);
-			debristhickness_input->GetInputValue(&debristhickness, gauss);
-			if((atan(fabs(slopex))*rad2deg)>30.){
-				f=1.;
-			}else{
-				f=0.;
-			}
-			//f=1;
-			//kappa=-5.6e16*smb*debristhickness*f;
-			//kappa=debristhickness/h*4e9*f;
-			//kappa=14.2809e8*f; // 25°
-			kappa=slope_threshold*1e8*f;
-			if(dim==2){
-				for(int i=0;i<numnodes;i++){
-					for(int j=0;j<numnodes;j++){
-						Ke->values[i*numnodes+j] +=  D_scalar*kappa*(
-								dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
-								);
-					}
-				}
-			}else{
-				for(int i=0;i<numnodes;i++){
-					for(int j=0;j<numnodes;j++){
-						Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i]);
-					}
-				}
-			}                
-		}
-
-		/*******************************************************************/                
-
 		IssmDouble rho;
 		if(FINITEELEMENT==P1Enum){
-			rho=2;
+			rho=2.;
 		}else if(FINITEELEMENT==P2Enum){
 			rho=4.;
 		}
-		if(stabilization==2){
+
+		for(int i=0;i<(dim*dim);i++) D[i]=0.;
+		if(stabilization==1){
+			/*SSA*/
+			if(dim==1){
+				vx_input->GetInputValue(&vx,gauss);
+				D[0]=h/rho*fabs(vx);
+			}else{
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel=sqrt(vx*vx+vy*vy);
+				D[0*dim+0]=h/rho*fabs(vx);
+				D[1*dim+1]=h/rho*fabs(vy);
+			}
+		}else if(stabilization==2){  
 			/*Streamline upwinding*/
 			if(dim==1){
 				vx_input->GetInputValue(&vx,gauss);
-				vel=fabs(vx)+1.e-10;
+				vel=fabs(vx)+EPS;
 				D[0] = h/(rho*vel)*vx*vx;
-			}
-			else{
-				vx_input->GetInputAverage(&vx);
-				vy_input->GetInputAverage(&vy);
-				vel=sqrt(vx*vx+vy*vy)+1.e-10;
+			}else{
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel=sqrt(vx*vx+vy*vy)+EPS;
 				D[0*dim+0]=h/(rho*vel)*vx*vx;
 				D[1*dim+0]=h/(rho*vel)*vy*vx;
 				D[0*dim+1]=h/(rho*vel)*vx*vy;
 				D[1*dim+1]=h/(rho*vel)*vy*vy;
-			}
-		}
-		else if(stabilization==1){  
-			/*SSA*/
-			if(dim==1){
-				vx_input->GetInputAverage(&vx);
-				D[0]=h/rho*fabs(vx);
-			}
-			else{
-				vx_input->GetInputAverage(&vx);
-				vy_input->GetInputAverage(&vy);
-				D[0*dim+0]=h/rho*fabs(vx);
-				D[1*dim+1]=h/rho*fabs(vy);
-			}
-		}
-		else if(stabilization==3){  
+			}		
+		}else if(stabilization==3){  
 			/*SUPG*/
 			if(dim==1){
-				vx_input->GetInputAverage(&vx);
-				tau=h/(rho*fabs(vx)+1e-10);
-			}
-			else{
-				vx_input->GetInputAverage(&vx);
-				vy_input->GetInputAverage(&vy);
-				tau=1*h/(rho*sqrt(vx*vx+vy*vy)+1e-10);
-			}
-		}
+				vx_input->GetInputValue(&vx,gauss);
+				tau=h/(rho*max(fabs(vx),EPS));
+			}else{
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				tau=h/(rho*sqrt(vx*vx+vy*vy)+EPS);
+			}
+		}
+
 		if(stabilization==1 || stabilization==2){
 			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
@@ -417,15 +292,13 @@
 						Ke->values[i*numnodes+j] += (
 								dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
-								dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
-								);
+								dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]));
 					}   
 				}
-			}
-			else{
+			}else{
 				for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
 			}
 		}else if(stabilization==3){ 
-			/*Mass matrix - part 2*/
 			if(dim==1){
+				/*Mass matrix - part 2*/
 				for(int i=0;i<numnodes;i++){
 					for(int j=0;j<numnodes;j++){
@@ -454,4 +327,5 @@
 				}
 
+
 				/*Advection matrix - part 2, B*/
 				for(int i=0;i<numnodes;i++){
@@ -464,5 +338,44 @@
 				for(int i=0;i<numnodes;i++){
 					for(int j=0;j<numnodes;j++){
-						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx);
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);;
+					}
+				}
+			}else if(dim==2){
+				/*Mass matrix - part 2*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+					}
+				}
+				/*Mass matrix - part 3*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(basis[i]*dvxdx+basis[i]*dvydy);
+					}
+				}
+
+				/*Advection matrix - part 2, A*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+					}
+				}
+				/*Advection matrix - part 3, A*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(basis[i]*dvxdx+basis[i]*dvydy);
+					}
+				}
+
+				/*Advection matrix - part 2, B*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+					}
+				}
+				/*Advection matrix - part 3, B*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);
 					}
 				}
@@ -487,6 +400,7 @@
 	int	stabilization,dim,domaintype;
 	IssmDouble  Jdet,dt;
-	IssmDouble  smb,thickness;
-	IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau,pf;
+	IssmDouble  smb,thickness,psi;
+	IssmDouble  vx,vy,vel,dvxdx,dvydy,h,tau,pf;
+	IssmDouble yts=31536000.;
 	IssmDouble  dvx[2],dvy[2];
 	IssmDouble* xyz_list = NULL;
@@ -533,4 +447,5 @@
 		vy_input=topelement->GetInput(VyDebrisEnum); _assert_(vy_input);
 	}
+	h=topelement->CharacteristicLength();
 
 	IssmDouble rho;
@@ -550,35 +465,46 @@
 		smb_input->GetInputValue(&smb,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
-
 		if(smb>0.){
-			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness-0.*dt*smb*pf)*basis[i];
-		} else {
-			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness-dt*smb*pf)*basis[i]; // take the negative of melt, because it is a debris production term here
-		}
+			psi=thickness-0.*dt*smb*pf;
+		}else{
+			psi=thickness-dt*smb*pf;
+		}
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*psi*basis[i]; 
 
 		if(stabilization==3){
 			/*SUPG*/
 			topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+			dvxdx=dvx[0];
 			if(dim==1){
-
 				vx_input->GetInputValue(&vx,gauss);
-				vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-				dvxdx=dvx[0];
-				tau=h/(rho*fabs(vx)+1e-10);
-				IssmDouble psi;
-				if(smb>0.){
-					psi=thickness;
-				} else {
-					psi=thickness-dt*smb*pf;
-				}
+				tau=h/(rho*max(fabs(vx),EPS));
+
 				/*Force vector - part 2*/
 				for(int i=0;i<numnodes;i++){
-					pe->values[i]+=Jdet*gauss->weight*psi*tau*(vx*dbasis[0*numnodes+i]);
+					pe->values[i]+=Jdet*gauss->weight*psi*tau*vx*dbasis[0*numnodes+i];
 				}
 				/*Force vector - part 3*/
 				for(int i=0;i<numnodes;i++){
-					pe->values[i]+=Jdet*gauss->weight*psi*tau*(basis[i]*dvxdx);
-				}
-
+					pe->values[i]+=Jdet*gauss->weight*psi*tau*basis[i]*dvxdx;
+				}
+
+			}else if(dim==2){
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+				vel=sqrt(vx*vx+vy*vy);
+				dvydy=dvy[1];
+				tau=h/(rho*vel+EPS);
+
+				/*Force vector - part 2*/
+				for(int i=0;i<numnodes;i++){
+					pe->values[i]+=Jdet*gauss->weight*psi*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+				}
+				/*Force vector - part 3*/
+				for(int i=0;i<numnodes;i++){
+					pe->values[i]+=Jdet*gauss->weight*psi*tau*(basis[i]*dvxdx+basis[i]*dvydy);
+				}
 			}
 		}
@@ -600,22 +526,12 @@
 }/*}}}*/
 void           DebrisAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	//element->InputUpdateFromSolutionOneDof(solution,DebrisThicknessEnum);
-	int*         ddoflist=NULL;
+
+	int *ddoflist=NULL;
 
 	int numnodes = element->GetNumberOfNodes();
-	IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
-	IssmDouble* thicknessold  = xNew<IssmDouble>(numnodes);
 	IssmDouble* newthickness  = xNew<IssmDouble>(numnodes);
-	IssmDouble* icethickness  = xNew<IssmDouble>(numnodes);
-	IssmDouble* bedslopex     = xNew<IssmDouble>(numnodes);
-	IssmDouble* surfaceslopex = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
 	IssmDouble minthickness = element->FindParam(DebrisMinThicknessEnum);
-	IssmDouble iceminthickness = element->FindParam(MasstransportMinThicknessEnum);
-	element->GetInputListOnVertices(&thickness[0],DebrisThicknessEnum);   
-	element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
-	element->GetInputListOnVertices(&bedslopex[0],BedSlopeXEnum);
-	element->GetInputListOnVertices(&surfaceslopex[0],SurfaceSlopeXEnum);
 	element->GetDofListLocal(&ddoflist,NoneApproximationEnum,GsetEnum);
 
@@ -625,26 +541,50 @@
 		if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
 
-		/* check for thickness<minthickness */
-		if(thickness[i]<minthickness) newthickness[i]=minthickness;
-
-		/* Carlos model sets all values below Hmin to zero */
-		if(icethickness[i]<=(iceminthickness+0.0001) & fabs((fabs(surfaceslopex[i])-fabs(bedslopex[i])))<1e-3 ) newthickness[i]=0;
-		//if(icethickness[i]<=(iceminthickness+0.01)) newthickness[i]=0;
-	}
-
-	/* update inputs */
+		// check for thickness<minthickness
+		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
+	}
+
+	// update inputs
 	element->AddInput(DebrisThicknessEnum,newthickness,P1Enum);
 
-	/* Free resources */
+	// Free resources
 	xDelete<IssmDouble>(newthickness);
-	xDelete<IssmDouble>(thickness);
-	xDelete<IssmDouble>(icethickness);
-	xDelete<IssmDouble>(bedslopex);
-	xDelete<IssmDouble>(surfaceslopex);
 	xDelete<int>(ddoflist);
+	//*/
 }/*}}}*/
 void           DebrisAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	//        return;
-	SetActiveNodesLSMx(femmodel);
+	//SetActiveNodesLSMx(femmodel);
+
+	// Update active elements based on ice levelset and ocean levelset*/
+	GetMaskOfIceVerticesLSMx(femmodel,false,true); //FIXME?
+	SetActiveNodesLSMx(femmodel,false,true); //FIXME?
+
+	/*Constrain all nodes that are grounded and unconstrain the ones that float*/
+	for(Object* & object : femmodel->elements->objects){
+		Element    *element  = xDynamicCast<Element*>(object);
+		int         numnodes  = element->GetNumberOfNodes();
+		IssmDouble *mask      = xNew<IssmDouble>(numnodes);
+		IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
+
+		element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
+		element->GetInputListOnNodes(&ls_active[0],DebrisMaskNodeActivationEnum);
+
+		for(int in=0;in<numnodes;in++){
+			Node* node=element->GetNode(in);
+			if(mask[in]>0. && ls_active[in]==1.){
+				// Do nothing
+				node->Activate(); //Not sure if we need this!
+			}
+			else{
+				IssmDouble phi=0;
+				node->Deactivate();// Not sure if we need this
+				node->ApplyConstraint(0,phi);
+			}
+		}
+		xDelete<IssmDouble>(mask);
+		xDelete<IssmDouble>(ls_active);
+	}
+	//*/
+	return;	
 }/*}}}*/
 void           DebrisAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
@@ -676,7 +616,10 @@
 			IssmDouble* slopey	   = xNew<IssmDouble>(numnodes); 
 			IssmDouble* onsurface	   = xNew<IssmDouble>(numnodes); 
+			IssmDouble* ls_active      = xNew<IssmDouble>(numnodes); 
 			element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
 			element->GetInputListOnNodes(icethickness,ThicknessEnum);
 			element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
+			element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
+
 			dim=1;
 			element->GetInputListOnNodes(slopex,SurfaceSlopeXEnum);
@@ -688,4 +631,5 @@
 			bool isminthicknessinelement=false;
 			bool remove_debris=false;
+			bool isactive=false;
 
 			IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);                        
@@ -693,72 +637,71 @@
 			switch(removalmodel){
 				case 1:{
-					IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
-
-					for(k=0; k<numnodes;k++){
-						if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
-					}
-					isminthicknessinelement=true;
-					if(isminthicknessinelement){
-						for(k=0; k<numnodes;k++){
-							if(onsurface[k]>0.5){
-								slope=fabs(slopex[k]);
-								if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
-								if((atan(slope)*rad2deg)>slope_threshold) debristhickness[k]=remove_debris=true;
-							}
-						}
-						if(remove_debris){
-							for(k=0; k<numnodes;k++){
-								if(icethickness[k]<=(iceminthickness+0.01)) debristhickness[k]=0.;
-							}
-						}
-					}
-					//int finite_element = element->GetElementType(); 
-					//element->AddInput(DebrisThicknessEnum,debristhickness,FINITEELEMENT);
-					element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
-
-					xDelete<IssmDouble>(debristhickness);
-					xDelete<IssmDouble>(icethickness);
-					xDelete<IssmDouble>(slopex);
-					xDelete<IssmDouble>(slopey);
-					break;
-						 }
+					       IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
+					       int kk=0;
+					       for(k=0; k<numnodes;k++){
+						       if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
+						       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
+					       }
+					       isminthicknessinelement=true;
+					       if(kk<numnodes && isminthicknessinelement){
+						       for(k=0; k<numnodes;k++){
+							       slope=fabs(slopex[k]);
+							       if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
+							       //slope_mean=slope_mean+slope;
+							       if((atan(slope)*rad2deg)>slope_threshold) remove_debris=true;
+							       //if((atan(slope)*rad2deg)>slope_threshold) debristhickness[k]=0.;
+						       }
+						       //if((atan(slope_mean)*rad2deg)>slope_threshold) remove_debris=true;
+						       if(remove_debris){
+							       for(k=0; k<numnodes;k++){
+								       debristhickness[k]=0.;
+							       }
+						       }
+					       }
+					       element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
+
+					       xDelete<IssmDouble>(debristhickness);
+					       xDelete<IssmDouble>(icethickness);
+					       xDelete<IssmDouble>(slopex);
+					       xDelete<IssmDouble>(slopey);
+					       break;
+				       }
 				case 2:{
-					IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
-					IssmDouble gravity = element->FindParam(ConstantsGEnum);
-					IssmDouble stress,rhod=1900.;
-
-					for(k=0; k<numnodes;k++){
-						if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
-					}
-					isminthicknessinelement=true;
-					if(isminthicknessinelement){
-						//stress=0;
-						int kk=0;
-						for(k=0; k<numnodes;k++){
-							if(onsurface[k]>0.5){
-								slope=fabs(slopex[k]);
-								if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
-								stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
-								if(stress>stress_threshold) debristhickness[k]=0.;
-								//kk++;
-							}
-						}
-						/*if((stress/double(kk))>stress_threshold) remove_debris=true;
-						  if(remove_debris){
-						  for(k=0; k<numnodes;k++){
-						  debristhickness[k]=0.;
-						  }
-						  }*/
-					}
-					//int finite_element = element->GetElementType(); 
-					//element->AddInput(DebrisThicknessEnum,debristhickness,FINITEELEMENT);
-					element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
-
-					xDelete<IssmDouble>(debristhickness);
-					xDelete<IssmDouble>(icethickness);
-					xDelete<IssmDouble>(slopex);
-					xDelete<IssmDouble>(slopey);
-					break;
-						 }
+					       IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
+					       IssmDouble gravity = element->FindParam(ConstantsGEnum);
+					       IssmDouble stress,rhod=1900.;
+					       int kk=0;
+					       for(k=0; k<numnodes;k++){
+						       if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
+						       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
+					       }
+					       isminthicknessinelement=true;
+					       if(kk<numnodes && isminthicknessinelement){
+						       //stress=0;
+						       IssmDouble stress_sum=0.;
+						       for(k=0; k<numnodes;k++){
+							       slope=fabs(slopex[k]);
+							       if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
+							       stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
+							       //stress_sum=stress_sum+stress;
+							       if(stress>stress_threshold) remove_debris=true;
+							       //if(stress>stress_threshold) debristhickness[k]=0.;
+						       }
+						       //if((stress_sum/double(kk))>stress_threshold) remove_debris=true;
+						       if(remove_debris){
+							       for(k=0; k<numnodes;k++){
+								       debristhickness[k]=0.;
+							       }
+						       }
+					       }
+					       element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
+
+					       xDelete<IssmDouble>(debristhickness);
+					       xDelete<IssmDouble>(icethickness);
+					       xDelete<IssmDouble>(slopex);
+					       xDelete<IssmDouble>(slopey);
+					       xDelete<IssmDouble>(ls_active);
+					       break;
+				       }
 				default: _error_("removalmodel "<<EnumToStringx(removalmodel)<<" not implemented yet");
 			}
@@ -767,138 +710,74 @@
 
 }/*}}}*/
-void           DebrisAnalysis::PreProcessing(FemModel* femmodel){/*{{{*/
+void DebrisAnalysis::PreProcessing(FemModel* femmodel){/*{{{*/
 
 	if(VerboseSolution()) _printf0_("   Debris preprocessing\n");
 
-	/*Intermediaries*/
-	bool isdebrisdisplacement=true;
-	int displacementmodel=1;
-	int k,numnodes;
-	int domaintype,dim;
-	femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
-	//femmodel->parameters->FindParam(&displacementmodel,DebrisDisplacementmodelEnum);
 	Element* element= NULL;
-	IssmDouble *xyz_list = NULL;
-	//IssmDouble top_normal[2];
-
-	if(isdebrisdisplacement){
-
-		//if(displacementmodel==0){
-		// no displacement, do nothing
-		//}else{
-		// deterministic or random displacement
-
-		for(Object* & object : femmodel->elements->objects){
-			element=xDynamicCast<Element*>(object);
-
-			numnodes=element->GetNumberOfNodes();
-			IssmDouble* debristhickness= xNew<IssmDouble>(numnodes);
-			IssmDouble* icethickness= xNew<IssmDouble>(numnodes);
-			IssmDouble* slopex	   = xNew<IssmDouble>(numnodes);
-			IssmDouble* slopey	   = xNew<IssmDouble>(numnodes); 
-			IssmDouble* vx 		   = xNew<IssmDouble>(numnodes);
-			IssmDouble* vy 		   = xNew<IssmDouble>(numnodes);
-			IssmDouble* surface        = xNew<IssmDouble>(numnodes);
-			IssmDouble* onsurface  	   = xNew<IssmDouble>(numnodes);
-			element->GetInputListOnNodes(vx,VxDebrisEnum);
-			element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
-			element->GetInputListOnNodes(icethickness,ThicknessEnum);
-			element->GetInputListOnNodes(surface,SurfaceEnum);
-			element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
-			element->GetVerticesCoordinates(&xyz_list);
-
-			dim=1;
-			element->GetInputListOnNodes(slopex,SurfaceSlopeXEnum);
-			if(domaintype!=Domain2DverticalEnum){
-				element->GetInputListOnNodes(slopey,SurfaceSlopeYEnum);
-				element->GetInputListOnNodes(vy,VyDebrisEnum);
-				dim=2;
-			}
-			IssmDouble slope,rad2deg=180./M_PI; //=57.2958
-			IssmDouble h=10.,f;
-			IssmDouble debrissum;
-			IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
-			bool displacedebris;
-
-			switch(displacementmodel){
-				case 1:{
-					IssmDouble slope_threshold = element->FindParam(DebrisRemovalSlopeThresholdEnum);
-					IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
-
-					bool isminthicknessinelement=false;
-					for(k=0; k<numnodes;k++){
-						if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
-					}
-					if(isminthicknessinelement){
-						//do nothing
-					}else{
-						debrissum=0.;
-						int test;
-						test=0;
-						for(k=0; k<numnodes;k++){
-							if(onsurface[k]>.5){
-								slope=pow(slopex[k]*slopex[k],0.5);
-								if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
-								if((atan(slope)*rad2deg)>30.){
-									f=0.5;
-									test=test+0;
-									debrissum=debrissum+debristhickness[k]*f;
-									//displacedebris=true;
-									//if(debristhickness[k]>1.e-6) vx[k]=vx[k]+10./31536000.;//*vx[k]/pow(pow(vx[k],2),0.5);
-									//debristhickness[k]=debristhickness[k]*(1.-f);
-								}
-							}
-						}
-						if(test>1){
-							test=test;
-						}else{
-							test=1;
-						}
-						//if(displacedebris){
-						int index=-1;
-						IssmDouble min=1e14;
-						for(k=0; k<numnodes;k++){
-							if(onsurface[k]>.5){
-								if(surface[k]<min){
-									index=k;
-									min=surface[k];
-								} 
-							}
-						}
-						for(k=0; k<numnodes;k++){
-							if(onsurface[k]>.5){
-								if(k==index){
-									debristhickness[k]=debristhickness[k]+debrissum;
-								}else{
-									debristhickness[k]=debristhickness[k]-debrissum;
-									if(debristhickness[k]<=0) debristhickness[k]=0;
-								}
-								//if(debristhickness[k]>10.) debristhickness[k]=10.;
-							}
-						}
-						//}
-					}
-
-					int finite_element = element->GetElementType(); 
-					//element->AddInput(DebrisThicknessEnum,debristhickness,finite_element);
-					//element->AddInput(VxDebrisEnum,vx,P1Enum);
-					element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
-
-					xDelete<IssmDouble>(debristhickness);
-					xDelete<IssmDouble>(icethickness);
-					xDelete<IssmDouble>(vx);
-					xDelete<IssmDouble>(vy); 
-					xDelete<IssmDouble>(slopex);
-					xDelete<IssmDouble>(slopey);
-					xDelete<IssmDouble>(surface);
-					break;
-						 }
-				case 2:
-					// Do nothing
-
-				default: _error_("Debris displacement model "<<EnumToStringx(displacementmodel)<<" not implemented yet");
-			}
-		}
-	}
-	//}
-}/*}}}*/        
+	for(Object* & object : femmodel->elements->objects){
+		element=xDynamicCast<Element*>(object);
+
+		int numvertices = element->GetNumberOfVertices();
+
+		IssmDouble* vx = xNew<IssmDouble>(numvertices);
+		IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
+		IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
+		IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
+		IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
+
+		element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
+		element->GetInputListOnVertices(&vx[0],VxDebrisEnum);
+		element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
+		element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
+		element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
+
+		IssmDouble slope,rad2deg=180./M_PI; //=57.2958
+		IssmDouble vslipx,rhod=1900.;
+		IssmDouble gravity=element->FindParam(ConstantsGEnum);
+		IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
+		IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
+
+		int step;
+		IssmDouble dt, maxv;
+		IssmDouble yts=31536000.;
+		femmodel->parameters->FindParam(&step,StepEnum);
+		femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+		bool isminthicknessinelement=false;
+		for(int i=0;i<numvertices;i++){
+			if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
+		}
+		if(isminthicknessinelement){
+			//do nothing
+			for(int i=0;i<numvertices;i++){
+				if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;
+			}
+		}else{
+			for(int i=0;i<numvertices;i++){
+				//if(onsurface[i]>.5){
+				slope=fabs(slopex[i]);
+				//if((atan(slope)*rad2deg)>25.){
+				//if(debristhickness[i]>0.01){
+				vslipx=1.0/yts;
+				//maxv=10.0/2./dt;
+				//vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts;
+				vx[i]=vx[i]+vslipx;
+				//debristhickness[i]=debristhickness[i];
+				//if(vx[i]>maxv) vx[i]=maxv;
+				//}
+				//} 
+				//}
+			}
+		}
+		//if(step%100==0)   
+		element->AddInput(VxDebrisEnum,vx,P1Enum);
+
+		/* Free resources */
+		xDelete<IssmDouble>(debristhickness);
+		xDelete<IssmDouble>(icethickness);
+		xDelete<IssmDouble>(vx);
+		xDelete<IssmDouble>(slopex);
+		xDelete<IssmDouble>(onsurface);
+	}
+
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27491)
@@ -2360,4 +2360,17 @@
 	Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
 	return (input->GetInputMax()<=0.);
+}
+/*}}}*/
+bool       Element::IsAllMinThicknessInElement(){/*{{{*/
+
+	IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum);
+
+	Input* input=this->GetInput(ThicknessEnum); _assert_(input);
+	if(input->GetInputMax()<=(minthickness+0.00000001)){
+		return true;
+	}
+	else{
+                return false;
+        }
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27490)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27491)
@@ -155,4 +155,5 @@
 		bool               IsOceanInElement();
 		bool               IsOceanOnlyInElement();
+		bool		   IsAllMinThicknessInElement();
 		bool               IsLandInElement();
 		void               Ismip6FloatingiceMeltingRate();
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 27491)
@@ -921,4 +921,7 @@
 			}
 			if(isdebris){
+				analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
+				analyses_temp[numanalyses++]=SmbAnalysisEnum;
+				analyses_temp[numanalyses++]=ExtrudeFromTopAnalysisEnum;
 				analyses_temp[numanalyses++]=DebrisAnalysisEnum;
 			}
Index: /issm/trunk-jpl/src/c/cores/debris_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/debris_core.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/cores/debris_core.cpp	(revision 27491)
@@ -34,28 +34,30 @@
 	if(VerboseSolution()) _printf0_("   computing debris transport\n");
 
-	// We need surface and bed slopes for removal model
+	// We need surface slopes for removal model
 	surfaceslope_core(femmodel);
-	bedslope_core(femmodel);
 
 	/*Transport Debris*/
 	if(VerboseSolution()) _printf0_("   call computational core\n");
-	//InputDuplicatex(femmodel,VxEnum,VxDebrisEnum);
 	femmodel->inputs->DuplicateInput(VxEnum,VxDebrisEnum);
-	//InputDuplicatex(femmodel,VyEnum,VyDebrisEnum);
+	if(domaintype!=Domain2DverticalEnum){
+		femmodel->inputs->DuplicateInput(VyEnum,VyDebrisEnum);	
+	}
+	femmodel->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
+	femmodel->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
+	depthaverage_core(femmodel);
+	if(domaintype!=Domain2DverticalEnum){
+		femmodel->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
+		femmodel->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
+		depthaverage_core(femmodel);
+	}
+
 	debris_analysis = new DebrisAnalysis();
-	//debris_analysis->PreCore(femmodel);
-	//femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
-	//extrudefromtop_core(femmodel);	
-
 	debris_analysis->Core(femmodel);
-	delete debris_analysis;
+	delete debris_analysis;	
 
 	femmodel->parameters->SetParam(DebrisThicknessEnum,InputToExtrudeEnum);
 	extrudefromtop_core(femmodel);	
-	//femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
-	//extrudefromtop_core(femmodel);
 
 	if(save_results) femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
-
 	if(solution_type==DebrisSolutionEnum)femmodel->RequestedDependentsx();
 
Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 27491)
@@ -10,5 +10,5 @@
 #include "../modules.h"
 
-void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
+void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
 	/* activate/deactivate nodes for levelset method according to IceMaskNodeActivation */
 
@@ -16,4 +16,5 @@
 	int nodeactivationmask = IceMaskNodeActivationEnum;
 	if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
+	if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
 
 
@@ -59,9 +60,10 @@
 }/*}}}*/
 
-void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology){/*{{{*/
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
 
 	/*Determine which node activation to construct*/
 	int nodeactivationmask = IceMaskNodeActivationEnum;
 	if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
+	if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
 
 	/*Initialize vector with number of vertices*/
@@ -83,6 +85,15 @@
 			}
 		}
-	}
-	else{
+	}else if(isdebris){
+		for(Object* & object : femmodel->elements->objects){
+                        Element* element=xDynamicCast<Element*>(object);
+                        if(element->IsIceInElement() && !element->IsAllMinThicknessInElement()){
+                                int nbv = element->GetNumberOfVertices();
+                                for(int iv=0;iv<nbv;iv++){
+                                        vec_mask_ice->SetValue(element->vertices[iv]->Pid(),1.,INS_VAL);
+                                }
+                        }
+                }
+	}else{
 		for(Object* & object : femmodel->elements->objects){
 			Element* element=xDynamicCast<Element*>(object);
@@ -101,5 +112,5 @@
 	delete vec_mask_ice;
 }/*}}}*/
-void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
+void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
 
 	/*Set configuration to levelset*/
@@ -107,6 +118,7 @@
 		/*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/
 		femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
-	}
-	else{
+	}else if(isdebris){
+		femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);
+	}else{
 		femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
 	}
@@ -115,4 +127,5 @@
 	int nodeactivationmask = IceMaskNodeActivationEnum;
 	if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
+	if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
 
 	/*Create vector on gset*/
@@ -137,6 +150,16 @@
 				xDelete<int>(glist_local);
 			}
-		}
-		else{
+		}else if(isdebris){
+			if(element->IsIceInElement() && !element->IsAllMinThicknessInElement()){
+                                int numnodes = element->GetNumberOfNodes();
+                                int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                IssmDouble* ones = xNew<IssmDouble>(gsize_local);
+                                for(int n=0;n<gsize_local;n++) ones[n] = 1.;
+                                vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
+                                xDelete<IssmDouble>(ones);
+                                xDelete<int>(glist_local);
+			}
+		}else{
 			if(element->IsIceInElement()){
 				int numnodes = element->GetNumberOfNodes();
Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h	(revision 27490)
+++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h	(revision 27491)
@@ -8,6 +8,6 @@
 #include "../../classes/classes.h"
 
-void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false);
-void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false);
-void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false);
+void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
+void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
 #endif  /* _SETACTIVENODESLSMX_H*/
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 27491)
@@ -522,5 +522,5 @@
 	IssmDouble AccG=0.1;                    // m w.e. /100m
 	IssmDouble AccMax=1.;                    // m w.e.
-	IssmDouble ReferenceElevation=2200.;     // m M&L
+	IssmDouble ReferenceElevation; 
 	IssmDouble AblationDays=120.;            //
 
@@ -537,4 +537,5 @@
         IssmDouble Alphad=0.07;            //              debris albedo                 0.07
         IssmDouble Alphai=0.4;             //              ice ablbedo
+        IssmDouble Alphaeff;
         IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
         IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
@@ -552,4 +553,7 @@
 	IssmDouble smb,yts,z,debris;
 	IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
+	bool isdebris;
+	int domaintype;
+	femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
 
 	/*Get material parameters and constants */
@@ -558,5 +562,5 @@
 	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); 
 	PhiD=0.;
-	//if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
+	if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
 
 	/* Loop over all the elements of this partition */
@@ -573,46 +577,56 @@
 		/* Get inputs */
 		element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
+		element->FindParam(&domaintype,DomainTypeEnum);		
 
 		/*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
 		for(int v=0;v<numvertices;v++){
 
-                        /*Get vertex elevation */
-                        z=surfacelist[v];
-
-                        /* Get debris cover */
-                        debris=debriscover[v];
-
-                        /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
-                        IssmDouble n=debris/dk;
-                        IssmDouble nmax=1000;
-                        IssmDouble Alphaeff;
-                        if(n>nmax){
-                                Alphaeff=Alphad;
-                        } else {
-                                Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
-                        }*/
-
-                        // M&L
-                       IssmDouble Alphaeff=Alphad;
-
-                        /* compute smb */
-                        for (int ismb=0;ismb<2;ismb++){
-                                if(ismb==0){
-                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
-                                        debris=0.;
-                                }else{
-                                	// only in the meltregime debris develops
-                                        if(-MassBalanceCmDayDebris<0) debris=debriscover[v]; 
-                                }
-                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
-                                    (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
-                                    (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
-                                    ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
-                                    ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
-                                    K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
-                                    ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
-                                    Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
+			/*Get vertex elevation */
+			z=surfacelist[v];
+
+			/*Get top element*/
+			//if(domaintype==Domain3DEnum){
+
+			//}else{
+			//	Alphaeff=Alphad;
+			//	ReferenceElevation=2200.;     // m M&L                        	
+			//}
+
+			/* compute smb */
+			for (int ismb=0;ismb<2;ismb++){
+				if(ismb==0){
+					// calc a reference smb to identify accum and melt region; debris only develops in ablation area
+					debris=0.;
+					PhiD=0.;
+				}else{
+					// only in the meltregime debris develops
+					if(-MassBalanceCmDayDebris<1e-14) debris=debriscover[v]; 
+				}
+				if(debris<=0.) debris=0.;
+				IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
+				IssmDouble n=debris/dk;
+				IssmDouble nmax=1000;
+				IssmDouble Alphaeff;
+				if(n>nmax){
+					Alphaeff=Alphad;
+				} else {
+					Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
+				}
+				ReferenceElevation=3200.;     // m HEF
+
+
+				Alphaeff=Alphad;
+				ReferenceElevation=2200.;     // m M&L  
+
+				MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
+							(Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
+							(Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+									WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
+										ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
+									((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+											WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
+									K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
+														ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
+											Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
                                     -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
                         }
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27490)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27491)
@@ -211,4 +211,5 @@
 syn keyword cConstant FrictionLinearizeEnum
 syn keyword cConstant FrictionPseudoplasticityExponentEnum
+syn keyword cConstant FrictionU0Enum
 syn keyword cConstant FrictionThresholdSpeedEnum
 syn keyword cConstant FrictionVoidRatioEnum
@@ -868,4 +869,5 @@
 syn keyword cConstant HydrologyWaterVyEnum
 syn keyword cConstant HydrologyMaskNodeActivationEnum
+syn keyword cConstant DebrisMaskNodeActivationEnum
 syn keyword cConstant IceEnum
 syn keyword cConstant IceMaskNodeActivationEnum
@@ -1724,4 +1726,5 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1730,6 +1733,6 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1744,6 +1747,6 @@
 syn keyword cType ElementInput
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1752,9 +1755,10 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1773,4 +1777,5 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1783,4 +1788,5 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1795,6 +1801,6 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1810,11 +1816,11 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType SealevelGeometry
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1835,8 +1841,4 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27490)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27491)
@@ -865,4 +865,5 @@
 	HydrologyWaterVyEnum,
 	HydrologyMaskNodeActivationEnum,
+	DebrisMaskNodeActivationEnum,
 	IceEnum,
 	IceMaskNodeActivationEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27491)
@@ -213,4 +213,5 @@
 		case FrictionLinearizeEnum : return "FrictionLinearize";
 		case FrictionPseudoplasticityExponentEnum : return "FrictionPseudoplasticityExponent";
+		case FrictionU0Enum : return "FrictionU0";
 		case FrictionThresholdSpeedEnum : return "FrictionThresholdSpeed";
 		case FrictionVoidRatioEnum : return "FrictionVoidRatio";
@@ -870,4 +871,5 @@
 		case HydrologyWaterVyEnum : return "HydrologyWaterVy";
 		case HydrologyMaskNodeActivationEnum : return "HydrologyMaskNodeActivation";
+		case DebrisMaskNodeActivationEnum : return "DebrisMaskNodeActivation";
 		case IceEnum : return "Ice";
 		case IceMaskNodeActivationEnum : return "IceMaskNodeActivation";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27490)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27491)
@@ -204,4 +204,5 @@
 syn keyword juliaConstC FrictionLinearizeEnum
 syn keyword juliaConstC FrictionPseudoplasticityExponentEnum
+syn keyword juliaConstC FrictionU0Enum
 syn keyword juliaConstC FrictionThresholdSpeedEnum
 syn keyword juliaConstC FrictionVoidRatioEnum
@@ -861,4 +862,5 @@
 syn keyword juliaConstC HydrologyWaterVyEnum
 syn keyword juliaConstC HydrologyMaskNodeActivationEnum
+syn keyword juliaConstC DebrisMaskNodeActivationEnum
 syn keyword juliaConstC IceEnum
 syn keyword juliaConstC IceMaskNodeActivationEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27490)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27491)
@@ -216,4 +216,5 @@
 	      else if (strcmp(name,"FrictionLinearize")==0) return FrictionLinearizeEnum;
 	      else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum;
+	      else if (strcmp(name,"FrictionU0")==0) return FrictionU0Enum;
 	      else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
 	      else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"HydrologyarmaNumBreaks")==0) return HydrologyarmaNumBreaksEnum;
 	      else if (strcmp(name,"HydrologyarmaNumParams")==0) return HydrologyarmaNumParamsEnum;
-	      else if (strcmp(name,"Hydrologyarmapolyparams")==0) return HydrologyarmapolyparamsEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum;
+	      if (strcmp(name,"Hydrologyarmapolyparams")==0) return HydrologyarmapolyparamsEnum;
+	      else if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum;
 	      else if (strcmp(name,"HydrologyAveraging")==0) return HydrologyAveragingEnum;
 	      else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
 	      else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
-	      else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
+	      if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
+	      else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
 	      else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
 	      else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
 	      else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
-	      else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
+	      if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
+	      else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
 	      else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
 	      else if (strcmp(name,"StochasticForcingIsEffectivePressure")==0) return StochasticForcingIsEffectivePressureEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
 	      else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
-	      else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
+	      if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
+	      else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
 	      else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
 	      else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum;
@@ -751,9 +752,9 @@
 	      else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
 	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
-	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"CalvingAblationrate")==0) return CalvingAblationrateEnum;
+	      if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
+	      else if (strcmp(name,"CalvingAblationrate")==0) return CalvingAblationrateEnum;
 	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
 	      else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
 	      else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
-	      else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
+	      if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
+	      else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
 	      else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
 	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
@@ -891,4 +892,5 @@
 	      else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
 	      else if (strcmp(name,"HydrologyMaskNodeActivation")==0) return HydrologyMaskNodeActivationEnum;
+	      else if (strcmp(name,"DebrisMaskNodeActivation")==0) return DebrisMaskNodeActivationEnum;
 	      else if (strcmp(name,"Ice")==0) return IceEnum;
 	      else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
@@ -996,10 +998,10 @@
 	      else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
 	      else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
-	      else if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
-	      else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
+	      if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
+	      else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
+	      else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
 	      else if (strcmp(name,"SealevelchangeGrot")==0) return SealevelchangeGrotEnum;
 	      else if (strcmp(name,"SealevelchangeGSatGravirot")==0) return SealevelchangeGSatGravirotEnum;
@@ -1119,10 +1121,10 @@
 	      else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
 	      else if (strcmp(name,"SmbT")==0) return SmbTEnum;
-	      else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
-	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
+	      if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
+	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
+	      else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
 	      else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
 	      else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
@@ -1242,10 +1244,10 @@
 	      else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
 	      else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
-	      else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
-	      else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
+	      if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
+	      else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
+	      else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
 	      else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
 	      else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
@@ -1365,10 +1367,10 @@
 	      else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
 	      else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
-	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
-	      else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
          else stage=12;
    }
    if(stage==12){
-	      if (strcmp(name,"IntInput")==0) return IntInputEnum;
+	      if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
+	      else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
+	      else if (strcmp(name,"IntInput")==0) return IntInputEnum;
 	      else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
 	      else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
@@ -1488,10 +1490,10 @@
 	      else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
 	      else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
-	      else if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
-	      else if (strcmp(name,"HydrologySubTime")==0) return HydrologySubTimeEnum;
          else stage=13;
    }
    if(stage==13){
-	      if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
+	      if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
+	      else if (strcmp(name,"HydrologySubTime")==0) return HydrologySubTimeEnum;
+	      else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
 	      else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum;
 	      else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum;
@@ -1611,10 +1613,10 @@
 	      else if (strcmp(name,"P1DG")==0) return P1DGEnum;
 	      else if (strcmp(name,"P1P1")==0) return P1P1Enum;
-	      else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
-	      else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
          else stage=14;
    }
    if(stage==14){
-	      if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
+	      if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
+	      else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
+	      else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
 	      else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
 	      else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
@@ -1734,10 +1736,10 @@
 	      else if (strcmp(name,"XYZ")==0) return XYZEnum;
 	      else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
-	      else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
-	      else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
          else stage=15;
    }
    if(stage==15){
-	      if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
+	      if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
+	      else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
+	      else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
 	      else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
 	      else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;