Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26442)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26443)
@@ -405,10 +405,10 @@
 
 	IssmDouble  calvingrate[NUMVERTICES];
-	IssmDouble  vx,vy,vel;
+	IssmDouble  vx,vy;
 	IssmDouble  water_height, bed,Ho,thickness,surface;
 	IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
-	IssmDouble  B, strainparallel, straineffective,n;
+	IssmDouble  strainparallel, straineffective,B,n;
 	IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
-	int crevasse_opening_stress;
+	int         crevasse_opening_stress;
 
 	/*retrieve the type of crevasse_opening_stress*/
@@ -442,6 +442,5 @@
 		bed_input->GetInputValue(&bed,&gauss);
 		surface_input->GetInputValue(&surface,&gauss);
-		strainrateparallel_input->GetInputValue(&strainparallel,&gauss);
-		strainrateeffective_input->GetInputValue(&straineffective,&gauss);
+
 		vx_input->GetInputValue(&vx,&gauss);
 		vy_input->GetInputValue(&vy,&gauss);
@@ -450,41 +449,45 @@
 		s_xy_input->GetInputValue(&s_xy,&gauss);
 		s_yy_input->GetInputValue(&s_yy,&gauss);
-		B_input->GetInputValue(&B,&gauss);
-		n_input->GetInputValue(&n,&gauss);
-
-		vel=sqrt(vx*vx+vy*vy)+1.e-14;
-
-		s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
-		s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
-		if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
-
-		Ho = thickness - (rho_seawater/rho_ice) * (-bed);
-		if(Ho<0.)  Ho=0.;
-
-		if(crevasse_opening_stress==0){		/*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
-			surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / n)-1)) / (rho_ice * constant_g);
-			basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/n)-1)) / (rho_ice*constant_g) - Ho);
-		}
-		else if(crevasse_opening_stress==1){	 /*Benn2017,Todd2018: maximum principal stress */
-			surface_crevasse[iv] = s1 / (rho_ice*constant_g);
-			basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
-		}
-
-		/* some constraints */
-		if (surface_crevasse[iv]<0.) {
+
+      /*Get longitudinal or maximum Eigen stress*/
+		if(crevasse_opening_stress==0){
+         /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
+			strainrateparallel_input->GetInputValue(&strainparallel,&gauss);
+			strainrateeffective_input->GetInputValue(&straineffective,&gauss);
+			B_input->GetInputValue(&B,&gauss);
+			n_input->GetInputValue(&n,&gauss);
+         s1 =  B * strainparallel * pow(straineffective, (1./n)-1);
+		}
+		else if(crevasse_opening_stress==1){
+         /*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;
+         }
+		}
+      else{
+         _error_("not supported");
+      }
+
+      /*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;
+		if(surface_crevasse[iv]<0.){
 			surface_crevasse[iv]=0.;
 			water_height = 0.;
 		}
-		if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
-		if (bed>0.) basal_crevasse[iv] = 0.;
-
-		//if (surface_crevasse[iv]<water_height){
-		//	water_height = surface_crevasse[iv];
-		//}
-
-		/* add water in surface crevasse */
-		surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */
-		crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */
-
+
+      /*Basal crevasse: sigma'_xx - rho_i g (H-d) - rho_w g (b+d) = 0*/
+      if(bed>0.){
+         basal_crevasse[iv] = 0.;
+      }
+      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);
+         if(basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
+      }
+
+      /*Total crevasse depth (surface + basal)*/
+		crevasse_depth[iv]   = surface_crevasse[iv] + basal_crevasse[iv];
 	}