Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27904)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27905)
@@ -29,5 +29,5 @@
 #define NUMVERTICES   3
 #define NUMVERTICES1D 2
-//#define MICI          1 //1 = DeConto & Pollard, 2 = DOMINOS
+#define MICI          2 //1 = DeConto & Pollard, 2 = Anna Crawford DOMINOS
 
 /*Constructors/destructor/copy*/
@@ -2772,13 +2772,12 @@
 IssmDouble Tria::GetIcefrontArea(){/*{{{*/
 
-	IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
+	IssmDouble  bed[NUMVERTICES];
 	IssmDouble	Haverage,frontarea;
 	IssmDouble  x1,y1,x2,y2,distance;
 	IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
 	int* indices=NULL;
-	IssmDouble* H=NULL;;
-	int nrfrontbed,numiceverts;
-
-	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
+	/*Return if no ice front present*/
+	if(!this->IsIcefront()) return 0.;
 
 	/*Retrieve all inputs and parameters*/
@@ -2788,91 +2787,87 @@
 	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
 
-	nrfrontbed=0;
-	for(int i=0;i<NUMVERTICES;i++){
-		/*Find if bed<0*/
-		if(bed[i]<0.) nrfrontbed++;
-	}
-
-	if(nrfrontbed==3){
-		/*2. Find coordinates of where levelset crosses 0*/
-		int         numiceverts;
-		IssmDouble  s[2],x[2],y[2];
-		this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
-		_assert_(numiceverts);
-
-		/*3 Write coordinates*/
-		IssmDouble  xyz_list[NUMVERTICES][3];
-		::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-		int counter = 0;
-		if((numiceverts>0) && (numiceverts<NUMVERTICES)){
-			for(int i=0;i<numiceverts;i++){
-				for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
-					x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
-					y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
-					counter++;
+	/*Only continue if all 3 vertices are below sea level*/
+	for(int i=0;i<NUMVERTICES;i++) if(bed[i]>=0.) return 0.;
+
+	/*2. Find coordinates of where levelset crosses 0*/
+	int         numiceverts;
+	IssmDouble  s[2],x[2],y[2];
+	this->GetLevelsetIntersection(&indices, &numiceverts, &s[0],MaskIceLevelsetEnum,0.);
+	_assert_(numiceverts);
+	if(numiceverts>2){
+		Input* ls_input = this->GetInput(MaskIceLevelsetEnum);
+		ls_input->Echo();
+	}
+
+	/*3 Write coordinates*/
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+	int counter = 0;
+	if((numiceverts>0) && (numiceverts<NUMVERTICES)){
+		for(int i=0;i<numiceverts;i++){
+			for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+				x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+				y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+				counter++;
+			}
+		}
+	}
+	else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+
+		for(int i=0;i<NUMVERTICES;i++){
+			if(lsf[indices[i]]==0.){
+				x[counter]=xyz_list[indices[i]][0];
+				y[counter]=xyz_list[indices[i]][1];
+				counter++;
+			}
+			if(counter==2) break;
+		}
+		if(counter==1){
+			/*We actually have only 1 vertex on levelset, write a single point as a segment*/
+			x[counter]=x[0];
+			y[counter]=y[0];
+			counter++;
+		}
+	}
+	else{
+		_error_("not sure what's going on here...");
+	}
+	x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
+	distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
+	if(distance<1e-3) return 0.;
+
+	IssmDouble H[4];
+	for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
+	xDelete<int>(indices);
+
+	switch(numiceverts){
+		case 1: // average over triangle
+			H[0]=Haux[0];
+			H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+			H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+			Haverage=(H[1]+H[2])/2;
+			break;
+		case 2: // average over quadrangle
+			H[0]=Haux[0];
+			H[1]=Haux[1];
+			H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+			H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+			Haverage=(H[2]+H[3])/2;
+			break;
+		case 3:
+			if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex
+			else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both
+				Haverage = 0;
+				for(int i=0;i<NUMVERTICES;i++){
+					if(lsf[indices[i]]==0.) Haverage -= Haux[indices[i]]/2;
+					if(Haverage<Haux[indices[i]]/2-1e-3) break; //done with the two vertices
 				}
 			}
-		}
-		else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
-
-			for(int i=0;i<NUMVERTICES;i++){
-				if(lsf[indices[i]]==0.){
-					x[counter]=xyz_list[indices[i]][0];
-					y[counter]=xyz_list[indices[i]][1];
-					counter++;
-				}
-				if(counter==2) break;
-			}
-			if(counter==1){
-				/*We actually have only 1 vertex on levelset, write a single point as a segment*/
-				x[counter]=x[0];
-				y[counter]=y[0];
-				counter++;
-			}
-		}
-		else{
-			_error_("not sure what's going on here...");
-		}
-		x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
-		distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
-
-		int numthk=numiceverts+2;
-		H=xNew<IssmDouble>(numthk);
-		for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
-
-		switch(numiceverts){
-			case 1: // average over triangle
-				H[0]=Haux[0];
-				H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
-				H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
-				Haverage=(H[1]+H[2])/2;
-				break;
-			case 2: // average over quadrangle
-				H[0]=Haux[0];
-				H[1]=Haux[1];
-				H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
-				H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
-				Haverage=(H[2]+H[3])/2;
-				break;
-			case 3:
-				if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex
-				else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both
-					Haverage = 0;
-					for(int i=0;i<NUMVERTICES;i++){
-						if(lsf[indices[i]]==0.) Haverage -= Haux[indices[i]]/2;
-						if(Haverage<Haux[indices[i]]/2-1e-3) break; //done with the two vertices
-					}
-				}
-				break;
-			default:
-				_error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
-				break;
-		}
-		frontarea=distance*Haverage;
-	}
-	else return 0;
-
-	xDelete<int>(indices);
-	xDelete<IssmDouble>(H);
+			break;
+		default:
+			_error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
+			break;
+	}
+	frontarea=distance*Haverage;
 
 	_assert_(frontarea>0);
@@ -4733,6 +4728,11 @@
 
 		/*Do we assume that the calving front does not move if MICI is not engaged?*/
-		movingfrontvx[iv] = 0.;
-		movingfrontvy[iv] = 0.;
+		bool regrowth = false;
+		bool apply_as_retreat = true;
+		if(!regrowth){
+			movingfrontvx[iv] = 0.;
+			movingfrontvy[iv] = 0.;
+		}
+
 		//movingfrontvx[iv] = -2000./(365*24*3600.)*dlsf[0]/norm_dlsf;
 		//movingfrontvy[iv] = -2000./(365*24*3600.)*dlsf[1]/norm_dlsf;
@@ -4746,4 +4746,8 @@
 		}
 		else if (MICI==2 && Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
+
+			/*if 1: RETREAT rate
+			 *if 0: calving rate*/
+			if(0) v[0]=0.; v[1]=0.;
 
 			/*5C Bn (worst case scenario)*/
@@ -4755,4 +4759,10 @@
 			movingfrontvx[iv] = v[0] -C*dlsf[0]/norm_dlsf;
 			movingfrontvy[iv] = v[1] -C*dlsf[1]/norm_dlsf;
+
+			/*disable regrowth if calving rate is too low*/
+			if(!regrowth && C<vel){
+				movingfrontvx[iv] = 0.;
+				movingfrontvy[iv] = 0.;
+			}
 		}
 	}
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 27904)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 27905)
@@ -30,8 +30,8 @@
 	/*retrieve mat_type option: */
 	#if PETSC_VERSION_LT(3,7,0)
+	PetscOptionsGetString(PETSC_NULL,"-mat_type",&option[0],100,&flag);
+	#elif PETSC_VERSION_LT(3,19,0)
 	PetscOptionsGetString(NULL,PETSC_NULL,"-mat_type",&option[0],100,&flag);
-	#elif PETSC_VERSION_LT(3,19,0)
-	PetscOptionsGetString(PETSC_NULL,"-mat_type",&option[0],100,&flag);
-	#else
+	#else /*newest version*/
 	PetscOptionsGetString(NULL,PETSC_NULLPTR,"-mat_type",&option[0],100,&flag);
 	#endif