Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 27297)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 27298)
@@ -284,4 +284,5 @@
 	./cores/smb_core.cpp \
 	./cores/bmb_core.cpp \
+	./cores/debris_core.cpp \
 	./solutionsequences/solutionsequence_thermal_nonlinear.cpp \
 	./modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp \
@@ -443,4 +444,7 @@
 issm_sources += ./modules/Damagex/Damagex.cpp
 endif
+if DEBRIS
+issm_sources += ./analyses/DebrisAnalysis.cpp
+endif
 if STRESSBALANCE
 issm_sources += ./analyses/StressbalanceAnalysis.cpp
Index: /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp	(revision 27298)
+++ /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp	(revision 27298)
@@ -0,0 +1,905 @@
+#include "./DebrisAnalysis.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../classes/Inputs/TransientInput.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../cores/cores.h"
+#include <math.h>
+
+#define FINITEELEMENT P1Enum
+
+/*Model processing*/
+void DebrisAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+
+	/*add constraints*/
+	IoModelToConstraintsx(constraints,iomodel,"md.debris.spcthickness",DebrisAnalysisEnum,FINITEELEMENT);
+
+}/*}}}*/
+void DebrisAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+
+	/*Intermediaries*/
+	int penpair_ids[2];
+	int count=0;
+	int numvertex_pairing;
+
+	/*Create Penpair for vertex_pairing: */
+	IssmDouble *vertex_pairing=NULL;
+	IssmDouble *nodeonsurface=NULL;
+	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.debris.vertex_pairing");
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,"md.mesh.vertexonsurface");
+	for(int i=0;i<numvertex_pairing;i++){
+
+		if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
+
+			/*In debugging mode, check that the second node is in the same cpu*/
+			_assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
+
+			/*Skip if one of the two is not on the bed*/
+			if(iomodel->domaintype!=Domain2DhorizontalEnum){
+				if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
+			}
+
+			/*Get node ids*/
+			penpair_ids[0]=reCast<int>(vertex_pairing[2*i+0]);
+			penpair_ids[1]=reCast<int>(vertex_pairing[2*i+1]);
+
+			/*Create Load*/
+			loads->AddObject(new Penpair(count+1,&penpair_ids[0]));
+			count++;
+		}
+	}
+
+	/*free ressources: */
+	iomodel->DeleteData(vertex_pairing,"md.debris.vertex_pairing");
+	iomodel->DeleteData(nodeonsurface,"md.mesh.vertexonsurface");
+}/*}}}*/
+void DebrisAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
+
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
+	::CreateNodes(nodes,iomodel,DebrisAnalysisEnum,FINITEELEMENT);
+	iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
+
+}/*}}}*/
+int  DebrisAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 1;
+}/*}}}*/
+void DebrisAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	/*Now, is the model 3d? otherwise, do nothing: */
+	if (iomodel->domaintype==Domain2DhorizontalEnum)return;
+
+	int smb_model;
+
+	/*Fetch data needed: */
+	iomodel->FindConstant(&smb_model,"md.smb.model");
+
+	/*Update elements: */
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(inputs,i,iomodel,analysis_counter,analysis_type,FINITEELEMENT);
+			counter++;
+		}
+	}
+
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxDebrisEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyDebrisEnum);
+	if(iomodel->domaintype!=Domain2DhorizontalEnum){
+		iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
+		iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
+	}
+
+	switch(smb_model){
+		case SMBforcingEnum:
+			iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
+			break;
+		default:
+			/*Nothing for now*/
+			;
+	}
+
+}/*}}}*/
+void DebrisAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+
+	int     numoutputs;
+	char**  requestedoutputs = NULL;
+
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.stabilization",DebrisStabilizationEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.removalmodel",DebrisRemovalmodelEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.displacementmodel",DebrisDisplacementmodelEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.min_thickness",DebrisMinThicknessEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.packingfraction",DebrisPackingFractionEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_slope_threshold",DebrisRemovalSlopeThresholdEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_stress_threshold",DebrisRemovalStressThresholdEnum));
+
+	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.debris.requested_outputs");
+	parameters->AddObject(new IntParam(DebrisNumRequestedOutputsEnum,numoutputs));
+	if(numoutputs)parameters->AddObject(new StringArrayParam(DebrisRequestedOutputsEnum,requestedoutputs,numoutputs));
+	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.debris.requested_outputs");
+
+}/*}}}*/
+
+/*Finite Element Analysis*/
+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]){
+				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);
+
+}/*}}}*/
+void           DebrisAnalysis::PreCore(FemModel* femmodel){/*{{{*/
+	_error_("not implemented");
+}/*}}}*/
+ElementVector* DebrisAnalysis::CreateDVector(Element* element){/*{{{*/
+	/*Default, return NULL*/
+	return NULL;
+}/*}}}*/
+ElementMatrix* DebrisAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
+	_error_("Not implemented");
+}/*}}}*/
+ElementMatrix* DebrisAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries */
+	int        stabilization,domaintype,dim;
+	IssmDouble Jdet,D_scalar,dt,h;
+	IssmDouble vel,vx,vy,dvxdx,dvydy;
+	IssmDouble xi,tau;
+	IssmDouble dvx[2],dvy[2];
+	Element*    topelement = NULL;
+	IssmDouble* xyz_list = NULL;
+
+	/*Get top element*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			topelement = element;
+			dim = 2;
+			break;
+		case Domain2DverticalEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 1;
+			break;
+		case Domain3DEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = topelement->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = topelement->NewElementMatrix(NoneApproximationEnum);
+	IssmDouble*   	basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble* 	dbasis = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*	D      = xNew<IssmDouble>(dim*dim);
+
+	/*Retrieve all inputs and parameters*/
+	topelement->GetVerticesCoordinates(&xyz_list);
+	topelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	topelement->FindParam(&stabilization,DebrisStabilizationEnum);
+	Input* vx_input=topelement->GetInput(VxDebrisEnum); _assert_(vx_input);
+	Input* vy_input=NULL;
+	if(dim>1){vy_input = topelement->GetInput(VyDebrisEnum); _assert_(vy_input);}
+	h = topelement->CharacteristicLength();
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=topelement->NewGauss(2);
+	while(gauss->next()){
+
+		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		topelement->NodalFunctions(basis,gauss);
+		topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+		if(dim==2) vy_input->GetInputValue(&vy,gauss);
+
+		D_scalar=gauss->weight*Jdet;
+		for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
+
+		/*Advection terms*/
+		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		D_scalar=dt*gauss->weight*Jdet;
+		if(dim==2){
+			vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+			dvxdx=dvx[0];
+			dvydy=dvy[1];
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					/*\phi_i \phi_j \nabla\cdot v*/
+					Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+					/*\phi_i v\cdot\nabla\phi_j*/
+					Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+				}
+			}
+		}
+		else{
+			dvxdx=dvx[0];
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D_scalar*dvxdx*basis[i]*basis[j];
+					Ke->values[i*numnodes+j] += D_scalar*vx*dbasis[0*numnodes+j]*basis[i];
+				}
+			}
+		}
+
+		/*******************************************************************/
+		/* 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;
+		}else if(FINITEELEMENT==P2Enum){
+			rho=4.;
+		}
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			if(dim==1){
+				vx_input->GetInputValue(&vx,gauss);
+				vel=fabs(vx)+1.e-10;
+				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;
+				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){  
+			/*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);
+			}
+		}
+		if(stabilization==1 || stabilization==2){
+			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
+			if(dim==2){
+				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*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]) 
+								);
+					}   
+				}
+			}
+			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){
+				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];
+					}
+				}
+				/*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;
+					}
+				}
+
+				/*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])*(vx*dbasis[0*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])*(basis[i]*dvxdx);
+					}
+				}
+
+				/*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)*(vx*dbasis[0*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);
+					}
+				}
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	delete gauss;
+	if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
+	return Ke;
+}/*}}}*/
+ElementVector* DebrisAnalysis::CreatePVector(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries */
+	int	stabilization,dim,domaintype;
+	IssmDouble  Jdet,dt;
+	IssmDouble  smb,thickness;
+	IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau,pf;
+	IssmDouble  dvx[2],dvy[2];
+	IssmDouble* xyz_list = NULL;
+	Element*    topelement = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			topelement = element;
+			dim = 2; 
+			break;
+		case Domain2DverticalEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 1; 
+			break;
+		case Domain3DEnum:           
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 2; break;
+		default: 
+			_error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = topelement->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = topelement->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	topelement->GetVerticesCoordinates(&xyz_list);
+	topelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	topelement->FindParam(&pf,DebrisPackingFractionEnum);
+	topelement->FindParam(&stabilization,DebrisStabilizationEnum);
+	Input* smb_input      	= topelement->GetInput(SmbMassBalanceEnum);  _assert_(smb_input);
+	Input* thickness_input  = topelement->GetInput(DebrisThicknessEnum); _assert_(thickness_input);
+	Input* vx_input  = topelement->GetInput(VxDebrisEnum);  _assert_(vx_input);
+	Input* vy_input=NULL;
+	if(dim==2){
+		vy_input=topelement->GetInput(VyDebrisEnum); _assert_(vy_input);
+	}
+
+	IssmDouble rho;
+	if(FINITEELEMENT==P1Enum){
+		rho=2.;
+	}else if(FINITEELEMENT==P2Enum){
+		rho=4.;
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=topelement->NewGauss(2);
+	while(gauss->next()){
+
+		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		topelement->NodalFunctions(basis,gauss);
+
+		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
+		}
+
+		if(stabilization==3){
+			/*SUPG*/
+			topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			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;
+				}
+				/*Force vector - part 2*/
+				for(int i=0;i<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);
+				}
+
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	delete gauss;
+	if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
+	return pe;
+}/*}}}*/
+void           DebrisAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void           DebrisAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
+	_error_("Not implemented yet");
+}/*}}}*/
+void           DebrisAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+	//element->InputUpdateFromSolutionOneDof(solution,DebrisThicknessEnum);
+	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);
+
+	for(int i=0;i<numnodes;i++){
+		newthickness[i] = solution[ddoflist[i]];
+		if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
+		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 */
+	element->AddInput(DebrisThicknessEnum,newthickness,P1Enum);
+
+	/* 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);
+}/*}}}*/
+void           DebrisAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
+
+	if(VerboseSolution()) _printf0_("   Debris postprocessing\n");
+
+	/*Intermediaries*/
+	int removalmodel;
+	int k,numnodes;
+	int domaintype,dim;
+	femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
+	femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum);
+	Element* element= NULL;
+	Element*    topelement = NULL;
+
+	if(removalmodel==0){
+		// no removal, do nothing
+	}else{
+		// slope or driving stress removal
+
+		for(Object* & object : femmodel->elements->objects){
+			element=xDynamicCast<Element*>(object);
+
+			numnodes=element->GetNumberOfNodes();
+			element->FindParam(&domaintype,DomainTypeEnum);
+			IssmDouble* icethickness   = xNew<IssmDouble>(numnodes);
+			IssmDouble* debristhickness= xNew<IssmDouble>(numnodes);
+			IssmDouble* slopex	   = xNew<IssmDouble>(numnodes);
+			IssmDouble* slopey	   = xNew<IssmDouble>(numnodes); 
+			IssmDouble* onsurface	   = xNew<IssmDouble>(numnodes); 
+			element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
+			element->GetInputListOnNodes(icethickness,ThicknessEnum);
+			element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
+			dim=1;
+			element->GetInputListOnNodes(slopex,SurfaceSlopeXEnum);
+			if(domaintype!=Domain2DverticalEnum){
+				element->GetInputListOnNodes(slopey,SurfaceSlopeYEnum);
+				dim=2;
+			}
+			IssmDouble slope,rad2deg=180./M_PI; //=57.2958
+			bool isminthicknessinelement=false;
+			bool remove_debris=false;
+
+			IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);                        
+
+			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]){
+								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;
+
+				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]){
+								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;
+
+				default: _error_("removalmodel "<<EnumToStringx(removalmodel)<<" not implemented yet");
+			}
+		}
+	}
+
+}/*}}}*/
+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]){
+								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]){
+								if(surface[k]<min){
+									index=k;
+									min=surface[k];
+								} 
+							}
+						}
+						for(k=0; k<numnodes;k++){
+							if(onsurface[k]){
+								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");
+			}
+		}
+	}
+	//}
+}/*}}}*/        
Index: /issm/trunk-jpl/src/c/analyses/analyses.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 27297)
+++ /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 27298)
@@ -17,4 +17,5 @@
 #include "./BalancevelocityAnalysis.h"
 #include "./DamageEvolutionAnalysis.h"
+#include "./DebrisAnalysis.h"
 #include "./DepthAverageAnalysis.h"
 #include "./EnthalpyAnalysis.h"
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 27298)
@@ -851,5 +851,5 @@
 		case TransientSolutionEnum:{
 			/*We have multiple analyses here, process one by one*/
-			bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isoceantransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslc,isesa,issampling;
+			bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isoceantransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslc,isesa,isdebris,issampling;
 			iomodel->FindConstant(&isthermal,"md.transient.isthermal");
 			iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
@@ -863,4 +863,5 @@
 			iomodel->FindConstant(&isslc,"md.transient.isslc");
 			iomodel->FindConstant(&isesa,"md.transient.isesa");
+			iomodel->FindConstant(&isdebris,"md.transient.isdebris");
 			iomodel->FindConstant(&issampling,"md.transient.issampling");
       int* analyses_iter     = NULL;
@@ -910,4 +911,7 @@
 				analyses_temp[numanalyses++]=EsaAnalysisEnum;
 			}
+			if(isdebris){
+				analyses_temp[numanalyses++]=DebrisAnalysisEnum;
+			}
 			if(issampling){
 				analyses_temp[numanalyses++]=SamplingAnalysisEnum;
@@ -993,4 +997,5 @@
 		if(profiler->Used(GROUNDINGLINECORE)) _printf0_("   "<<setw(40)<<left<<"Groundingline migration core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(GROUNDINGLINECORE) << " sec\n");
 		if(profiler->Used(ESACORE)) _printf0_("   "<<setw(40)<<left<<"ESA core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(ESACORE) << " sec\n");
+		if(profiler->Used(DEBRISCORE)) _printf0_("   "<<setw(40)<<left<<"DEBRIS core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(SLRCORE) << " sec\n");
 		if(profiler->Used(SLRCORE)) _printf0_("   "<<setw(40)<<left<<"SLR core elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(SLRCORE) << " sec\n");
 		if(profiler->Used(MPISERIAL)) _printf0_("   "<<setw(40)<<left<<"MPISERIAL elapsed time:"<<setw(7)<<setprecision(6)<<profiler->TotalTime(MPISERIAL) << " sec\n");
Index: /issm/trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 27298)
@@ -159,5 +159,6 @@
 	}
 	if(
-				analysis_enum==FreeSurfaceTopAnalysisEnum
+			analysis_enum==FreeSurfaceTopAnalysisEnum ||
+			analysis_enum==DebrisAnalysisEnum
 	  ){
 		if(iomodel->domaintype!=Domain2DhorizontalEnum){
Index: /issm/trunk-jpl/src/c/classes/Profiler.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Profiler.h	(revision 27297)
+++ /issm/trunk-jpl/src/c/classes/Profiler.h	(revision 27298)
@@ -27,14 +27,15 @@
 #define ESACORE				14 /*Profiling ESA */
 #define SLRCORE				15 /*Profiling SLR */
-#define SAMPLINGCORE	16 /*Profiling SAMPLING */
-#define MPISERIAL				17 /*Profiling MPISerial */
-#define SEDLOOP				18 /*Profiling MPISerial */
-#define SEDMatrix				19 /*Profiling MPISerial */
-#define SEDUpdate				20 /*Profiling MPISerial */
-#define EPLLOOP				21 /*Profiling MPISerial */
-#define EPLMasking			22 /*Profiling MPISerial */
-#define EPLMatrices			23 /*Profiling MPISerial */
-#define EPLUpdate				24 /*Profiling MPISerial */
-#define MAXPROFSIZE			25 /*Used to initialize static arrays*/
+#define DEBRISCORE                              16 /*Profiling DEBRIS */
+#define SAMPLINGCORE	17 /*Profiling SAMPLING */
+#define MPISERIAL				18 /*Profiling MPISerial */
+#define SEDLOOP				19 /*Profiling MPISerial */
+#define SEDMatrix				20 /*Profiling MPISerial */
+#define SEDUpdate				21 /*Profiling MPISerial */
+#define EPLLOOP				22 /*Profiling MPISerial */
+#define EPLMasking			23 /*Profiling MPISerial */
+#define EPLMatrices			24 /*Profiling MPISerial */
+#define EPLUpdate				25 /*Profiling MPISerial */
+#define MAXPROFSIZE			26 /*Used to initialize static arrays*/
 
 class Profiler: public Object{
Index: /issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- /issm/trunk-jpl/src/c/cores/cores.h	(revision 27297)
+++ /issm/trunk-jpl/src/c/cores/cores.h	(revision 27298)
@@ -57,4 +57,5 @@
 void damage_core(FemModel* femmodel);
 void sampling_core(FemModel* femmodel);
+void debris_core(FemModel* femmodel);
 
 /*sealevel change cores:*/
Index: /issm/trunk-jpl/src/c/cores/debris_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/debris_core.cpp	(revision 27298)
+++ /issm/trunk-jpl/src/c/cores/debris_core.cpp	(revision 27298)
@@ -0,0 +1,67 @@
+/*!\file: debris_core.cpp
+ * \brief: core of the debris solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../classes/Inputs/TransientInput.h"
+
+void debris_core(FemModel* femmodel){ /*{{{*/
+
+	/*Start profiler*/
+	femmodel->profiler->Start(DEBRISCORE);
+
+	/*parameters: */
+	int    numoutputs,domaintype;
+	bool   save_results;
+	int    solution_type,stabilization;
+	char** requested_outputs = NULL;
+	DebrisAnalysis * debris_analysis = NULL;
+
+	/*activate configuration*/
+	femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&numoutputs,DebrisNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DebrisRequestedOutputsEnum);
+
+	if(VerboseSolution()) _printf0_("   computing debris transport\n");
+
+	// We need surface and bed 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);
+	debris_analysis = new DebrisAnalysis();
+	//debris_analysis->PreCore(femmodel);
+	//femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
+	//extrudefromtop_core(femmodel);	
+
+	debris_analysis->Core(femmodel);
+	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();
+
+	/*Free ressources:*/
+	if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
+
+	/*profiler*/
+	femmodel->profiler->Stop(DEBRISCORE);
+} /*}}}*/
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 27298)
@@ -78,5 +78,5 @@
 		if(VerboseSolution()){
 			_printf0_("iteration " << step << "/" << ceil((finaltime-time)/dt)+step << \
-						"  time [yr]: " <<std::fixed<<setprecision(2)<< time/yts << " (time step: " << dt/yts << ")\n");
+					"  time [yr]: " <<std::fixed<<setprecision(2)<< time/yts << " (time step: " << dt/yts << ")\n");
 		}
 		bool save_results=false;
@@ -101,5 +101,5 @@
 		if(amr_frequency){
 
-			#if !defined(_HAVE_AD_)
+#if !defined(_HAVE_AD_)
 			if(save_results) femmodel->WriteMeshInResults();
 			if(step%amr_frequency==0 && time<finaltime){
@@ -108,7 +108,7 @@
 			}
 
-			#else
+#else
 			_error_("AMR not suppored with AD");
-			#endif
+#endif
 		}
 
@@ -132,4 +132,5 @@
 	bool isstressbalance,ismasstransport,isage,isoceantransport,issmb,isthermal,isgroundingline,isesa,issampling;
 	bool isslc,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,isstochasticforcing,save_results;
+	bool isdebris;
 	int  step,sb_coupling_frequency;
 	int  domaintype,numoutputs;
@@ -153,11 +154,12 @@
 	femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
 	femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
+	femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
 	femmodel->parameters->FindParam(&issampling,TransientIssamplingEnum);
 	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
 	femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
 
-	#if defined(_HAVE_OCEAN_ )
+#if defined(_HAVE_OCEAN_ )
 	if(isoceancoupling) OceanExchangeDatax(femmodel,false);
-	#endif
+#endif
 
 	if(isstochasticforcing){
@@ -208,4 +210,8 @@
 	}
 
+	if(isdebris){
+		debris_core(femmodel);
+	}
+
 	/* from here on, prepare geometry for next time step*/
 
@@ -229,9 +235,9 @@
 	/*Sea level change: */
 	if(isslc){
-		#ifdef _HAVE_SEALEVELCHANGE_
+#ifdef _HAVE_SEALEVELCHANGE_
 		sealevelchange_core(femmodel);
-		#else
+#else
 		_error_("Compiled with SeaLevelChange capability");
-		#endif
+#endif
 	}
 
@@ -263,18 +269,18 @@
 	femmodel->parameters->FindParam(&isslc,TransientIsslcEnum);
 
-	#if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
+#if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
 	if(amr_frequency){
 		femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
 		if(amr_restart) femmodel->ReMesh();
 	}
-	#endif
-
-	#if defined(_HAVE_OCEAN_ )
+#endif
+
+#if defined(_HAVE_OCEAN_ )
 	if(isoceancoupling) OceanExchangeDatax(femmodel,true);
-	#endif
-
-	#if defined(_HAVE_SEALEVELCHANGE_)
+#endif
+
+#if defined(_HAVE_SEALEVELCHANGE_)
 	if(isslc) sealevelchange_initialgeometry(femmodel);
-	#endif
+#endif
 }/*}}}*/
 
@@ -335,5 +341,5 @@
 		if(VerboseSolution()){
 			_printf0_("iteration " << step << "/" << ceil((finaltime-time)/dt)+step << \
-						"  time [yr]: " <<std::fixed<<setprecision(2)<< time/yts << " (time step: " << dt/yts << ")\n");
+					"  time [yr]: " <<std::fixed<<setprecision(2)<< time/yts << " (time step: " << dt/yts << ")\n");
 		}
 
@@ -383,11 +389,11 @@
 
 	/*Start tracing*/
-	#if _CODIPACK_MAJOR_==2
+#if _CODIPACK_MAJOR_==2
 	auto& tape_codi = IssmDouble::getTape();
-	#elif _CODIPACK_MAJOR_==1
+#elif _CODIPACK_MAJOR_==1
 	auto& tape_codi = IssmDouble::getGlobalTape();
-	#else
-	#error "_CODIPACK_MAJOR_ not supported"
-	#endif
+#else
+#error "_CODIPACK_MAJOR_ not supported"
+#endif
 	tape_codi.setActive();
 
@@ -410,15 +416,15 @@
 
 		tape_codi.registerOutput(J);
-		#if _CODIPACK_MAJOR_==2
+#if _CODIPACK_MAJOR_==2
 		codi_global.output_indices.push_back(J.getIdentifier());
-		#elif _CODIPACK_MAJOR_==1
+#elif _CODIPACK_MAJOR_==1
 		codi_global.output_indices.push_back(J.getGradientData());
-		#else
-		#error "_CODIPACK_MAJOR_ not supported"
-		#endif
+#else
+#error "_CODIPACK_MAJOR_ not supported"
+#endif
 
 		/*Keep track of output for printing*/
 		Jlist[count] = output_value.getValue();
-      count++;
+		count++;
 	}
 	Jlist[count] = J.getValue();
@@ -473,5 +479,5 @@
 			if(VerboseSolution()){
 				_printf0_("step "<<thisstep<<" ("<<ii+1<<"/"<<checkpoint_frequency<<") time [yr]: "\
-							<<std::fixed<<std::setprecision(2)<<thistime/yts<< " (time step: " << thisdt/yts << ")\n");
+						<<std::fixed<<std::setprecision(2)<<thistime/yts<< " (time step: " << thisdt/yts << ")\n");
 			}
 
@@ -524,5 +530,5 @@
 	xDelete<double>(Yb);
 	xDelete<int>(Yin);
-   return J.getValue();
+	return J.getValue();
 }/*}}}*/
 #endif
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 27298)
@@ -618,5 +618,5 @@
 		xDelete<int>(vertex_pairing);
 	}
-	if(!isamr && (analysis==MasstransportAnalysisEnum || analysis==FreeSurfaceBaseAnalysisEnum || analysis==FreeSurfaceTopAnalysisEnum)){
+	if(!isamr && (analysis==MasstransportAnalysisEnum || analysis==FreeSurfaceBaseAnalysisEnum || analysis==FreeSurfaceTopAnalysisEnum || analysis==DebrisAnalysisEnum)){
 		int *vertex_pairing = NULL;
 		int  numvertex_pairing;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 27298)
@@ -90,4 +90,5 @@
 		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));
+		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isdebris",TransientIsdebrisEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.transient.issampling",TransientIssamplingEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.isstochasticforcing",StochasticForcingIsStochasticForcingEnum));
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27297)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27298)
@@ -156,4 +156,14 @@
 	DamageStressThresholdEnum,
 	DamageStressUBoundEnum,
+	DebrisMinThicknessEnum,
+        DebrisNumRequestedOutputsEnum,
+        DebrisPenaltyFactorEnum,
+        DebrisRequestedOutputsEnum,
+        DebrisStabilizationEnum,
+        DebrisRemovalmodelEnum,
+        DebrisDisplacementmodelEnum,
+        DebrisRemovalSlopeThresholdEnum,
+        DebrisRemovalStressThresholdEnum,
+        DebrisPackingFractionEnum,
 	DebugProfilingEnum,
 	DebrisThicknessEnum,
@@ -593,4 +603,5 @@
 	TransientIsageEnum,
 	TransientIsdamageevolutionEnum,
+	TransientIsdebrisEnum,
 	TransientIsesaEnum,
 	TransientIsgiaEnum,
@@ -1121,4 +1132,5 @@
 	VxAverageEnum,
 	VxBaseEnum,
+	VxDebrisEnum,
 	VxEnum,
 	VxMeshEnum,
@@ -1128,4 +1140,5 @@
 	VyAverageEnum,
 	VyBaseEnum,
+	VyDebrisEnum,
 	VyEnum,
 	VyMeshEnum,
@@ -1318,4 +1331,6 @@
 	DataSetParamEnum,
 	DatasetInputEnum,
+	DebrisAnalysisEnum,
+	DebrisSolutionEnum,
 	DefaultAnalysisEnum,
 	DefaultCalvingEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27298)
@@ -164,4 +164,14 @@
 		case DamageStressThresholdEnum : return "DamageStressThreshold";
 		case DamageStressUBoundEnum : return "DamageStressUBound";
+		case DebrisMinThicknessEnum : return "DebrisMinThickness";
+		case DebrisNumRequestedOutputsEnum : return "DebrisNumRequestedOutputs";
+		case DebrisPenaltyFactorEnum : return "DebrisPenaltyFactor";
+		case DebrisRequestedOutputsEnum : return "DebrisRequestedOutputs";
+		case DebrisStabilizationEnum : return "DebrisStabilization";
+		case DebrisRemovalmodelEnum : return "DebrisRemovalmodel";
+		case DebrisDisplacementmodelEnum : return "DebrisDisplacementmodel";
+		case DebrisRemovalSlopeThresholdEnum : return "DebrisRemovalSlopeThreshold";
+		case DebrisRemovalStressThresholdEnum : return "DebrisRemovalStressThreshold";
+		case DebrisPackingFractionEnum : return "DebrisPackingFraction";
 		case DebugProfilingEnum : return "DebugProfiling";
 		case DebrisThicknessEnum : return "DebrisThickness";
@@ -601,4 +611,5 @@
 		case TransientIsageEnum : return "TransientIsage";
 		case TransientIsdamageevolutionEnum : return "TransientIsdamageevolution";
+		case TransientIsdebrisEnum : return "TransientIsdebris";
 		case TransientIsesaEnum : return "TransientIsesa";
 		case TransientIsgiaEnum : return "TransientIsgia";
@@ -1126,4 +1137,5 @@
 		case VxAverageEnum : return "VxAverage";
 		case VxBaseEnum : return "VxBase";
+		case VxDebrisEnum : return "VxDebris";
 		case VxEnum : return "Vx";
 		case VxMeshEnum : return "VxMesh";
@@ -1133,4 +1145,5 @@
 		case VyAverageEnum : return "VyAverage";
 		case VyBaseEnum : return "VyBase";
+		case VyDebrisEnum : return "VyDebris";
 		case VyEnum : return "Vy";
 		case VyMeshEnum : return "VyMesh";
@@ -1321,4 +1334,6 @@
 		case DataSetParamEnum : return "DataSetParam";
 		case DatasetInputEnum : return "DatasetInput";
+		case DebrisAnalysisEnum : return "DebrisAnalysis";
+		case DebrisSolutionEnum : return "DebrisSolution";
 		case DefaultAnalysisEnum : return "DefaultAnalysis";
 		case DefaultCalvingEnum : return "DefaultCalving";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27297)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27298)
@@ -167,4 +167,14 @@
 	      else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
 	      else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum;
+	      else if (strcmp(name,"DebrisMinThickness")==0) return DebrisMinThicknessEnum;
+	      else if (strcmp(name,"DebrisNumRequestedOutputs")==0) return DebrisNumRequestedOutputsEnum;
+	      else if (strcmp(name,"DebrisPenaltyFactor")==0) return DebrisPenaltyFactorEnum;
+	      else if (strcmp(name,"DebrisRequestedOutputs")==0) return DebrisRequestedOutputsEnum;
+	      else if (strcmp(name,"DebrisStabilization")==0) return DebrisStabilizationEnum;
+	      else if (strcmp(name,"DebrisRemovalmodel")==0) return DebrisRemovalmodelEnum;
+	      else if (strcmp(name,"DebrisDisplacementmodel")==0) return DebrisDisplacementmodelEnum;
+	      else if (strcmp(name,"DebrisRemovalSlopeThreshold")==0) return DebrisRemovalSlopeThresholdEnum;
+	      else if (strcmp(name,"DebrisRemovalStressThreshold")==0) return DebrisRemovalStressThresholdEnum;
+	      else if (strcmp(name,"DebrisPackingFraction")==0) return DebrisPackingFractionEnum;
 	      else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
 	      else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum;
@@ -250,5 +260,8 @@
 	      else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
 	      else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
-	      else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
 	      else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
 	      else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
@@ -260,8 +273,5 @@
 	      else if (strcmp(name,"HydrologydcSedimentThickness")==0) return HydrologydcSedimentThicknessEnum;
 	      else if (strcmp(name,"HydrologyStepAdapt")==0) return HydrologyStepAdaptEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
+	      else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
 	      else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum;
 	      else if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum;
@@ -373,5 +383,8 @@
 	      else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
 	      else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
-	      else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
 	      else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
 	      else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
@@ -383,8 +396,5 @@
 	      else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
 	      else if (strcmp(name,"QmuVariableDescriptors")==0) return QmuVariableDescriptorsEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum;
+	      else if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum;
 	      else if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum;
 	      else if (strcmp(name,"QmuVariablePartitionsNt")==0) return QmuVariablePartitionsNtEnum;
@@ -496,5 +506,8 @@
 	      else if (strcmp(name,"SmbARMAInitialTime")==0) return SmbARMAInitialTimeEnum;
 	      else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum;
-	      else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum;
 	      else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum;
 	      else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
@@ -506,8 +519,5 @@
 	      else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum;
 	      else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
+	      else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
 	      else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
 	      else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
@@ -613,4 +623,5 @@
 	      else if (strcmp(name,"TransientIsage")==0) return TransientIsageEnum;
 	      else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
+	      else if (strcmp(name,"TransientIsdebris")==0) return TransientIsdebrisEnum;
 	      else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
 	      else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
@@ -618,5 +629,8 @@
 	      else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
 	      else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
-	      else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
 	      else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
 	      else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum;
@@ -629,8 +643,5 @@
 	      else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
 	      else if (strcmp(name,"Velocity")==0) return VelocityEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"Xxe")==0) return XxeEnum;
+	      else if (strcmp(name,"Xxe")==0) return XxeEnum;
 	      else if (strcmp(name,"Yye")==0) return YyeEnum;
 	      else if (strcmp(name,"Zze")==0) return ZzeEnum;
@@ -741,5 +752,8 @@
 	      else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
 	      else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
-	      else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
 	      else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
 	      else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
@@ -752,8 +766,5 @@
 	      else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
 	      else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum;
+	      else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum;
 	      else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum;
 	      else if (strcmp(name,"Dummy")==0) return DummyEnum;
@@ -864,5 +875,8 @@
 	      else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
 	      else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
-	      else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum;
 	      else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum;
 	      else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
@@ -875,8 +889,5 @@
 	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
 	      else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"Misfit")==0) return MisfitEnum;
+	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
 	      else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
 	      else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
@@ -987,5 +998,8 @@
 	      else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum;
 	      else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
-	      else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
 	      else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
 	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
@@ -998,8 +1012,5 @@
 	      else if (strcmp(name,"SmbC")==0) return SmbCEnum;
 	      else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
+	      else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
 	      else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
 	      else if (strcmp(name,"SmbD")==0) return SmbDEnum;
@@ -1110,5 +1121,8 @@
 	      else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
 	      else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
-	      else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
 	      else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
 	      else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
@@ -1121,8 +1135,5 @@
 	      else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
 	      else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Surface")==0) return SurfaceEnum;
+	      else if (strcmp(name,"Surface")==0) return SurfaceEnum;
 	      else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
 	      else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
@@ -1153,4 +1164,5 @@
 	      else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
 	      else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
+	      else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum;
 	      else if (strcmp(name,"Vx")==0) return VxEnum;
 	      else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
@@ -1160,4 +1172,5 @@
 	      else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
 	      else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
+	      else if (strcmp(name,"VyDebris")==0) return VyDebrisEnum;
 	      else if (strcmp(name,"Vy")==0) return VyEnum;
 	      else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
@@ -1231,5 +1244,8 @@
 	      else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
 	      else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
-	      else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
 	      else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
 	      else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
@@ -1244,8 +1260,5 @@
 	      else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum;
 	      else if (strcmp(name,"Outputdefinition67")==0) return Outputdefinition67Enum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum;
+	      else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum;
 	      else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
 	      else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
@@ -1351,6 +1364,11 @@
 	      else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
 	      else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
+	      else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum;
+	      else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum;
 	      else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
-	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
 	      else if (strcmp(name,"Dense")==0) return DenseEnum;
 	      else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
@@ -1367,8 +1385,5 @@
 	      else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
 	      else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
+	      else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
 	      else if (strcmp(name,"Element")==0) return ElementEnum;
 	      else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
@@ -1475,5 +1490,8 @@
 	      else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
 	      else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
-	      else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
 	      else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
 	      else if (strcmp(name,"LoveLt")==0) return LoveLtEnum;
@@ -1490,8 +1508,5 @@
 	      else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
 	      else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"Masscon")==0) return MassconEnum;
+	      else if (strcmp(name,"Masscon")==0) return MassconEnum;
 	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
 	      else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
@@ -1598,5 +1613,8 @@
 	      else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
 	      else if (strcmp(name,"Scaled")==0) return ScaledEnum;
-	      else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
+         else stage=14;
+   }
+   if(stage==14){
+	      if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
 	      else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
 	      else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum;
@@ -1613,8 +1631,5 @@
 	      else if (strcmp(name,"Separate")==0) return SeparateEnum;
 	      else if (strcmp(name,"Seq")==0) return SeqEnum;
-         else stage=14;
-   }
-   if(stage==14){
-	      if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
+	      else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
 	      else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum;
 	      else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
Index: /issm/trunk-jpl/src/m/classes/debris.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/debris.m	(revision 27298)
+++ /issm/trunk-jpl/src/m/classes/debris.m	(revision 27298)
@@ -0,0 +1,172 @@
+%DEBRIS class definition
+%
+%   Usage:
+%      debris=debris();
+
+classdef debris
+	properties (SetAccess=public)
+		 spcthickness           = NaN;
+		 min_thickness          = 0;
+		 stabilization          = 0;
+		 packingfraction	= 0;
+		 removalmodel		= 0;
+		 displacementmodel	= 0;
+		 removal_slope_threshold= 0;
+		 removal_stress_threshold=0;
+		 vertex_pairing         = NaN;
+		 requested_outputs      = {};
+	end
+	methods (Static)
+		function self = loadobj(self) % {{{
+			% This function is directly called by matlab when a model object is
+			% loaded. If the input is a struct it is an old version of this class and
+			% old fields must be recovered (make sure they are in the deprecated
+			% model properties)
+
+
+
+			if verLessThan('matlab','7.9'),
+				disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
+				disp('         if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');
+
+				% This is a Matlab bug: all the fields of md have their default value
+				% Example of error message:
+				% Warning: Error loading an object of class 'model':
+				% Undefined function or method 'exist' for input arguments of type 'cell'
+				%
+				% This has been fixed in MATLAB 7.9 (R2009b) and later versions
+			end
+
+			if isstruct(self)
+				disp('Recovering debris from older version');
+				self = structtoobj(debris(),self);
+			end
+		end% }}}
+	end
+	methods
+		function self = debris(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				case 1
+					inputstruct=varargin{1};
+					list1 = properties('debris');
+					list2 = fieldnames(inputstruct);
+					for i=1:length(list1)
+						fieldname = list1{i};
+						if ismember(fieldname,list2),
+							self.(fieldname) = inputstruct.(fieldname);
+						end
+					end
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.spcthickness=project3d(md,'vector',self.spcthickness,'type','node');
+		end % }}}
+		function list = defaultoutputs(self,md) % {{{
+
+			list = {'DebrisThickness'};
+
+		end % }}}
+		function self = setdefaultparameters(self) % {{{
+
+			%Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
+			self.stabilization=2;
+
+			%Minimum debris thickness that can be used
+			self.min_thickness=0;
+
+			%Fraction of debris covered in the ice
+			self.packingfraction=0.01;
+
+			%Type of frontal debris removal
+			self.removalmodel=0;
+
+			%Type of debris displacement
+			self.displacementmodel=0;
+
+			%Slope threshold for removalmodel (1)
+                 	self.removal_slope_threshold=0;
+                 	
+			%Stress threshold for removalmodel (2)
+			self.removal_stress_threshold=0;
+
+			%default output
+			self.requested_outputs={'default'};
+		end % }}}
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+
+			%Early return
+			if ~ismember('MasstransportAnalysis',analyses) |  (strcmp(solution,'TransientSolution') & md.transient.isdebris==0), return; end
+
+			md = checkfield(md,'fieldname','debris.spcthickness','Inf',1,'timeseries',1);
+			md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3]);
+			md = checkfield(md,'fieldname','debris.min_thickness','>=',0);
+			md = checkfield(md,'fieldname','debris.packingfraction','>=',0);
+			md = checkfield(md,'fieldname','debris.removalmodel','values',[0 1 2]);
+			md = checkfield(md,'fieldname','debris.displacementmodel','values',[0 1 2]);
+			md = checkfield(md,'fieldname','debris.removal_slope_threshold','>=',0);
+			md = checkfield(md,'fieldname','debris.removal_stress_threshold','>=',0);
+			md = checkfield(md,'fieldname','debris.requested_outputs','stringrow',1);
+			if ~any(isnan(md.stressbalance.vertex_pairing)),
+				md = checkfield(md,'fieldname','stressbalance.vertex_pairing','>',0);
+			end
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   Debris solution parameters:'));
+			fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
+			fielddisplay(self,'min_thickness','minimum ice thickness allowed [m]');
+			fielddisplay(self,'packingfraction','fraction of debris covered in the ice');
+			fielddisplay(self,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: streamline upwind Petrov-Galerkin (SUPG)');
+			fielddisplay(self,'removalmodel','frontal removal of debris. 0: no removal, 1: Slope-triggered debris removal, 2: driving-stress triggered debris removal');
+			fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: ...');
+			fielddisplay(self,'removal_slope_threshold','critical slope (degrees) for removalmodel (1)');
+			fielddisplay(self,'removal_stress_threshold','critical stress (Pa) for removalmodel (2)');
+
+			disp(sprintf('\n      %s','Penalty options:'));
+			fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized');
+			fielddisplay(self,'requested_outputs','additional outputs requested');
+
+		end % }}}
+		function marshall(self,prefix,md,fid) % {{{
+
+			yts=md.constants.yts;
+
+			WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','removalmodel','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','displacementmodel','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','removal_slope_threshold','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','removal_stress_threshold','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','packingfraction','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','vertex_pairing','format','DoubleMat','mattype',3);
+
+			%process requested outputs
+			outputs = self.requested_outputs;
+			pos  = find(ismember(outputs,'default'));
+			if ~isempty(pos),
+				outputs(pos) = [];                         %remove 'default' from outputs
+				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
+			end
+			WriteData(fid,prefix,'data',outputs,'name','md.debris.requested_outputs','format','StringArray');
+		end % }}}
+		function savemodeljs(self,fid,modelname) % {{{
+
+			writejs1Darray(fid,[modelname '.debris.spcthickness'],self.spcthickness);
+			writejsdouble(fid,[modelname '.debris.min_thickness'],self.min_thickness);
+			writejsdouble(fid,[modelname '.debris.stabilization'],self.stabilization);
+			writejsdouble(fid,[modelname '.debris.removalmodel'],self.removalmodel);
+			writejsdouble(fid,[modelname '.debris.displacementmodel'],self.displacementmodel);
+			writejsdouble(fid,[modelname '.debris.removal_slope_threshold'],self.removal_slope_threshold);
+			writejsdouble(fid,[modelname '.debris.removal_stress_threshold'],self.removal_stress_threshold);
+			writejsdouble(fid,[modelname '.debris.packingfraction'],self.packingfraction);
+			writejs2Darray(fid,[modelname '.debris.vertex_pairing'],self.vertex_pairing);
+			writejsdouble(fid,[modelname '.debris.penalty_factor'],self.penalty_factor);
+			writejscellstring(fid,[modelname '.debris.requested_outputs'],self.requested_outputs);
+
+		end % }}}
+	end
+end