Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26463)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26464)
@@ -29,4 +29,5 @@
 #define NUMVERTICES   3
 #define NUMVERTICES1D 2
+//#define ISMICI        1
 
 /*Constructors/destructor/copy*/
@@ -406,9 +407,24 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 	IssmDouble  vx,vy;
-	IssmDouble  water_height, bed,Ho,thickness,surface;
+	IssmDouble  water_height, bed,Hab,thickness,surface;
 	IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
 	IssmDouble  strainparallel, straineffective,B,n;
 	IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
 	int         crevasse_opening_stress;
+
+
+	/*reset if no ice in element*/
+	if(!this->IsIceInElement()){
+		for(int i=0;i<NUMVERTICES;i++){
+			surface_crevasse[i] = 0.;
+			basal_crevasse[i] = 0.;
+			crevasse_depth[i] = 0.;
+		}
+		this->AddInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1DGEnum);
+		this->AddInput(BasalCrevasseEnum,&basal_crevasse[0],P1DGEnum);
+		this->AddInput(CrevasseDepthEnum,&crevasse_depth[0],P1DGEnum);
+		return;
+	}
+
 
 	/*retrieve the type of crevasse_opening_stress*/
@@ -462,7 +478,5 @@
          /*Benn2017,Todd2018: maximum principal stress */
 			Matrix2x2Eigen(&s1,&s2,NULL,NULL,s_xx,s_xy,s_yy);
-			if(fabs(s2)>fabs(s1)){
-            stmp=s2; s2=s1; s1=stmp;
-         }
+			s1 = max(0.,max(s1,s2));
 		}
       else{
@@ -471,5 +485,5 @@
 
       /*Surface crevasse: sigma'_xx - rho_i g d + rho_fw g d_w = 0*/
-      surface_crevasse[iv] = s1 / (rho_ice*constant_g) + (rho_freshwater/rho_ice)*water_height;
+      surface_crevasse[iv] = 2*s1 / (rho_ice*constant_g) + (rho_freshwater/rho_ice)*water_height;
 		if(surface_crevasse[iv]<0.){
 			surface_crevasse[iv]=0.;
@@ -482,7 +496,7 @@
       }
       else{
-         Ho = thickness - (rho_seawater/rho_ice) * (-bed);
-         if(Ho<0.)  Ho=0.;
-         basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
+         Hab = thickness - (rho_seawater/rho_ice) * (-bed);
+         if(Hab<0.)  Hab=0.;
+         basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (2*s1/ (rho_ice*constant_g)-Hab);
          if(basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
       }
@@ -4323,4 +4337,5 @@
 	}
 
+	#ifdef MICI
 	/**************************************  MICI  START ************************************/
 	/*MICI from Crawford et al. 2021:
@@ -4331,30 +4346,45 @@
 	 * Tice = 5C;  Bn  ->  I = 1.9e-16; α = 7.3
 	 */
-	//Input* surface_input = this->GetInput(SurfaceEnum); _assert_(surface_input);
-	//Input* bed_input     = this->GetInput(BedEnum);     _assert_(bed_input);
-	//Input* ls_input      = this->GetInput(MaskIceLevelsetEnum);   _assert_(ls_input);
-	//IssmDouble Hc,bed,ls;
-	//for(int iv=0;iv<NUMVERTICES;iv++){
-	//	gauss.GaussVertex(iv);
-
-	//	surface_input->GetInputValue(&Hc,&gauss);
-	//	bed_input->GetInputValue(&bed,&gauss);
-	//	ls_input->GetInputValue(&ls,&gauss);
-
-	//	if(Hc>135. && bed<0. && fabs(ls)<100.e3){
-
-	//		lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
-	//		norm_dlsf=0.;
-	//		for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-	//		norm_dlsf=sqrt(norm_dlsf);
-
-	//		/*5C Bn (worst case scenario)*/
-	//		IssmDouble I     = 1.9e-16;
-	//		IssmDouble alpha = 7.3;
-	//		IssmDouble C = I*pow(Hc,alpha)/(24*3600.); /*convert from m/day to m/s*/
-	//		movingfrontvx[iv] = -C*dlsf[0]/norm_dlsf;
-	//		movingfrontvy[iv] = -C*dlsf[1]/norm_dlsf;
-	//	}
-	//}
+	Input* surface_input = this->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* bed_input     = this->GetInput(BedEnum);     _assert_(bed_input);
+	Input* ls_input      = this->GetInput(MaskIceLevelsetEnum);   _assert_(ls_input);
+	IssmDouble Hc,bed,ls;
+	for(int iv=0;iv<NUMVERTICES;iv++){
+		gauss.GaussVertex(iv);
+
+		surface_input->GetInputValue(&Hc,&gauss);
+		bed_input->GetInputValue(&bed,&gauss);
+		ls_input->GetInputValue(&ls,&gauss);
+
+		/*Do we assume that the calving front does not move?*/
+		//movingfrontvx[iv] = 0.;
+		//movingfrontvy[iv] = 0.;
+
+		//if(Hc>80. && bed<0. && fabs(ls)<100.e3){ //Pollard & De Conto
+		if(Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
+
+			lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
+			norm_dlsf=0.;
+			for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
+			norm_dlsf=sqrt(norm_dlsf);
+
+			/*use vel direction instead of LSF*/
+			vx_input->GetInputValue(&v[0],&gauss);
+			vy_input->GetInputValue(&v[1],&gauss);
+			vel=sqrt(v[0]*v[0] + v[1]*v[1]);
+			norm_dlsf = max(vel,1.e-10);
+			dlsf[0] = v[0];
+			dlsf[1] = v[1];
+
+			/*5C Bn (worst case scenario)*/
+			IssmDouble I     = 1.9e-16;
+			IssmDouble alpha = 7.3;
+			IssmDouble C = min(2000.,I*pow(Hc,alpha))/(24*3600.); /*convert from m/day to m/s*/
+			//IssmDouble C = (min(max(Hc,80.),100.) - 80.)/20. * 10./(24*3600.); /*Original MICI! convert from m/day to m/s*/
+			movingfrontvx[iv] = -C*dlsf[0]/norm_dlsf;
+			movingfrontvy[iv] = -C*dlsf[1]/norm_dlsf;
+		}
+	}
+	#endif
 	/**************************************  END MICI  *************************************/