Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26325)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26326)
@@ -4054,5 +4054,5 @@
 }
 /*}}}*/
-void	   Tria::MovingFrontalVelocity(void){/*{{{*/
+void	      Tria::MovingFrontalVelocity(void){/*{{{*/
 
 	int  dim, domaintype, calvinglaw, i;
@@ -4319,4 +4319,40 @@
 		movingfrontvy[iv] = w[1];		
 	}
+
+	/**************************************  MICI  START ************************************/
+	/*MICI from Crawford et al. 2021:
+	 * Tice = 20C; Bf  ->  I = 3.7e-16; α = 6.9
+	 * Tice = 20C; Bn  ->  I = 5.1e-14; α = 6.0
+	 * Tice = 20C; Bh  ->  I = 3.2e-17; α = 7.2
+	 * Tice = 10C; Bn  ->  I = 6.9e-17; α = 7.3
+	 * 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;
+	//	}
+	//}
+	/**************************************  END MICI  *************************************/
 
 	/*Add input*/