Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17393)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17394)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 101
+//#define FSANALYTICAL 4
 
 /*Model processing*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17393)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17394)
@@ -129,4 +129,36 @@
 void Element::DeleteMaterials(void){/*{{{*/
 	delete this->material;
+}/*}}}*/
+IssmDouble Element::Divergence(void){/*{{{*/
+	/*Compute element divergence*/
+
+	/*Intermediaries*/
+	IssmDouble Jdet;
+	IssmDouble divergence=0.;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble *xyz_list = NULL;
+
+	/*Get inputs and parameters*/
+	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
+	this->GetVerticesCoordinates(&xyz_list);
+
+	Gauss* gauss=this->NewGaussBase(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+		/*Get strain rate assuming that epsilon has been allocated*/
+		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+
+		divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet;
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return divergence;
 }/*}}}*/
 void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17393)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17394)
@@ -56,4 +56,5 @@
 		void       CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
 		void       DeleteMaterials(void);
+		IssmDouble Divergence(void);
 		void       FindParam(bool* pvalue,int paramenum);
 		void       FindParam(int* pvalue,int paramenum);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 17393)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 17394)
@@ -432,4 +432,5 @@
 	switch (response_descriptor_enum){
 
+		case DivergenceEnum:               this->Divergencex(responses); break;
 		case IceVolumeEnum:                this->IceVolumex(responses); break;
 		case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break;
@@ -504,4 +505,5 @@
 
 				/*Scalar output*/
+				case DivergenceEnum:               this->Divergencex(&double_result);               break;
 				case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
 				case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
@@ -1042,4 +1044,20 @@
 
 }/*}}}*/
+void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/
+
+	IssmDouble local_divergence=0;
+	IssmDouble total_divergence;
+
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		local_divergence+=element->Divergence();
+	}
+	ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+	/*Assign output pointers: */
+	*pdiv=total_divergence;
+
+}/*}}}*/
 void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 17393)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 17394)
@@ -72,4 +72,5 @@
 		void MinVzx(IssmDouble* presponse);
 		void TotalSmbx(IssmDouble* pSmb);
+		void Divergencex(IssmDouble* pdiv);
 		void IceVolumex(IssmDouble* pV);
 		void IceVolumeAboveFloatationx(IssmDouble* pV);
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17393)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17394)
@@ -533,4 +533,5 @@
 	StressTensoryzEnum,
 	StressTensorzzEnum,
+	DivergenceEnum,
 	GiaCrossSectionShapeEnum,
 	GiadWdtEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17393)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17394)
@@ -522,4 +522,5 @@
 		case StressTensoryzEnum : return "StressTensoryz";
 		case StressTensorzzEnum : return "StressTensorzz";
+		case DivergenceEnum : return "Divergence";
 		case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape";
 		case GiadWdtEnum : return "GiadWdt";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17393)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17394)
@@ -534,4 +534,5 @@
 	      else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
 	      else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+	      else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
 	      else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
 	      else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
 	      else if (strcmp(name,"Regular")==0) return RegularEnum;
-	      else if (strcmp(name,"Scaled")==0) return ScaledEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Separate")==0) return SeparateEnum;
+	      if (strcmp(name,"Scaled")==0) return ScaledEnum;
+	      else if (strcmp(name,"Separate")==0) return SeparateEnum;
 	      else if (strcmp(name,"Sset")==0) return SsetEnum;
 	      else if (strcmp(name,"Verbose")==0) return VerboseEnum;