Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19061)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19062)
@@ -694,4 +694,5 @@
 		case DivergenceEnum:               this->Divergencex(responses); break;
 		case MaxDivergenceEnum:            this->MaxDivergencex(responses); break;
+		case IceMassEnum:                  this->IceMassx(responses); break;
 		case IceVolumeEnum:                this->IceVolumex(responses); break;
 		case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break;
@@ -777,4 +778,5 @@
 					case DivergenceEnum:               this->Divergencex(&double_result);               break;
 					case MaxDivergenceEnum:            this->MaxDivergencex(&double_result);            break;
+					case IceMassEnum:                this->IceMassx(&double_result);                break;
 					case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
 					case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
@@ -1394,4 +1396,20 @@
 
 }/*}}}*/
+void FemModel::IceMassx(IssmDouble* pM){/*{{{*/
+
+	IssmDouble local_ice_mass = 0;
+	IssmDouble total_ice_mass;
+
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+		local_ice_mass+=element->IceMass();
+	}
+	ISSM_MPI_Reduce(&local_ice_mass,&total_ice_mass,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&total_ice_mass,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+	/*Assign output pointers: */
+	*pM=total_ice_mass;
+
+}/*}}}*/
 void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 19061)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 19062)
@@ -79,4 +79,5 @@
 		void Divergencex(IssmDouble* pdiv);
 		void MaxDivergencex(IssmDouble* pdiv);
+		void IceMassx(IssmDouble* pV);
 		void IceVolumex(IssmDouble* pV);
 		void IceVolumeAboveFloatationx(IssmDouble* pV);