Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 13720)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 13721)
@@ -521,7 +521,5 @@
 					     ./modules/IceVolumex/IceVolumex.cpp\
 					     ./modules/ElementResponsex/ElementResponsex.h\
-					     ./modules/ElementResponsex/ElementResponsex.cpp\
-					     ./modules/MassFluxx/MassFluxx.cpp\
-					     ./modules/MassFluxx/MassFluxx.h
+					     ./modules/ElementResponsex/ElementResponsex.cpp
 #}}}
 #Slope sources  {{{
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13720)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13721)
@@ -228,4 +228,47 @@
 }
 /*}}}*/
+/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
+void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
+
+	/*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
+	 * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
+	 * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
+	 * Slope configuration.*/
+
+	int found=-1;
+	for(int i=0;i<nummodels;i++){
+		if (analysis_type_list[i]==configuration_type){
+			found=i;
+			break;
+		}
+	}
+	if(found!=-1) analysis_counter=found;
+	else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
+
+	/*Now, plug analysis_counter and analysis_type inside the parameters: */
+	this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
+	this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
+	this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
+
+	/*configure elements, loads and nodes, for this new analysis: */
+	this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
+	this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
+	this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
+
+	#ifdef _HAVE_PETSC_
+	/*take care of petsc options, that depend on this analysis type (present only after model processor)*/
+	if(this->parameters->Exist(PetscOptionsStringsEnum)){
+		PetscOptionsFromAnalysis(this->parameters,analysis_type);
+		if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
+	}
+	#endif
+
+}
+/*}}}*/
+/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
+void FemModel::SetCurrentConfiguration(int configuration_type){
+	this->SetCurrentConfiguration(configuration_type,configuration_type);
+}
+/*}}}*/
 /*FUNCTION FemModel::Solve {{{*/
 void FemModel::Solve(void){
@@ -280,51 +323,6 @@
 /*}}}*/
 
-/*Numerics: */
-/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
-void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
-
-	/*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
-	 * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
-	 * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
-	 * Slope configuration.*/
-
-	int found=-1;
-	for(int i=0;i<nummodels;i++){
-		if (analysis_type_list[i]==configuration_type){
-			found=i;
-			break;
-		}
-	}
-	if(found!=-1) analysis_counter=found;
-	else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
-
-	/*Now, plug analysis_counter and analysis_type inside the parameters: */
-	this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
-	this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
-	this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
-
-	/*configure elements, loads and nodes, for this new analysis: */
-	this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
-	this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
-	this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
-
-	#ifdef _HAVE_PETSC_
-	/*take care of petsc options, that depend on this analysis type (present only after model processor)*/
-	if(this->parameters->Exist(PetscOptionsStringsEnum)){
-		PetscOptionsFromAnalysis(this->parameters,analysis_type);
-		if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
-	}
-	#endif
-
-}
-/*}}}*/
-/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
-void FemModel::SetCurrentConfiguration(int configuration_type){
-	this->SetCurrentConfiguration(configuration_type,configuration_type);
-}
-/*}}}*/
-
 /*Modules:*/
-int FemModel::UpdateVertexPositionsx(void){ /*{{{*/
+int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
 	int     i;
@@ -378,36 +376,4 @@
 	NodesDofx(nodes,parameters,analysis_type);
 
-}
-/*}}}*/
-void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
-
-	int      i;
-
-	/*output: */
-	IssmDouble   dt;
-
-	/*intermediary: */
-	Element *element     = NULL;
-	IssmDouble   min_dt      = 0;
-	IssmDouble   node_min_dt = 0;
-
-	/*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
-	element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
-
-	for (i=1;i<elements->Size();i++){
-		element=(Element*)elements->GetObjectByOffset(i);
-		dt=element->TimeAdapt();
-		if(dt<min_dt)min_dt=dt;
-	}
-
-	/*Figure out minimum across the cluster: */
-	#ifdef _HAVE_MPI_
-	MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
-	MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
-	min_dt=node_min_dt;
-	#endif
-
-	/*Assign output pointers:*/
-	*pdt=min_dt;
 }
 /*}}}*/
@@ -440,5 +406,5 @@
 		case MaxVzEnum:                  MaxVzx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
 		case MaxAbsVzEnum:               MaxAbsVzx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
-		case MassFluxEnum:               MassFluxx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
+		case MassFluxEnum:               this->MassFluxx(          responses,process_units); break;
 		case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
 		case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
@@ -543,4 +509,98 @@
 }
 /*}}}*/
+void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
+
+	int      i;
+
+	/*output: */
+	IssmDouble   dt;
+
+	/*intermediary: */
+	Element *element     = NULL;
+	IssmDouble   min_dt      = 0;
+	IssmDouble   node_min_dt = 0;
+
+	/*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
+	element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
+
+	for (i=1;i<elements->Size();i++){
+		element=(Element*)elements->GetObjectByOffset(i);
+		dt=element->TimeAdapt();
+		if(dt<min_dt)min_dt=dt;
+	}
+
+	/*Figure out minimum across the cluster: */
+#ifdef _HAVE_MPI_
+	MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
+	MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
+	min_dt=node_min_dt;
+#endif
+
+	/*Assign output pointers:*/
+	*pdt=min_dt;
+}
+/*}}}*/
+#ifdef _HAVE_RESPONSES_
+void FemModel::MassFluxx(IssmDouble* pmass_flux,bool process_units){/*{{{*/
+
+	int          i,j;
+	Element     *element       = NULL;
+	int          element_id;
+	bool         ispresent     = false;
+	IssmDouble   mass_flux     = 0;
+	IssmDouble   all_mass_flux = 0;
+	int          counter;
+	IssmDouble **array         = NULL;
+	int          M;
+	int         *mdims_array   = NULL;
+	int         *ndims_array   = NULL;
+	IssmDouble  *segments      = NULL;
+	int          num_segments;
+
+	/*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */
+	this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum);
+	if(!ispresent)_error_("no mass flux segments available!");
+	this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum);
+
+	/*Retrieve index of segments being used for MassFlux computation: */
+	parameters->FindParam(&counter,IndexEnum);
+
+	/*retrieve segments from array: */
+	segments     = array[counter-1]; //matlab to "C" indexing
+	num_segments = mdims_array[counter-1];
+
+	/*Go through segments, and then elements, and figure out which elements belong to a segment. 
+	 * When we find one, use the element to compute the mass flux on the segment: */
+	for(i=0;i<num_segments;i++){
+		element_id=reCast<int,IssmDouble>(*(segments+5*i+4));
+		for(j=0;j<elements->Size();j++){
+			element=(Element*)this->elements->GetObjectByOffset(j);
+			if (element->Id()==element_id){
+				/*We found the element which owns this segment, use it to compute the mass flux: */
+				mass_flux+=element->MassFlux(segments+5*i+0,process_units);
+				break;
+			}
+		}
+	}
+
+#ifdef _HAVE_MPI_
+	MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
+	mass_flux=all_mass_flux;
+#endif
+
+	/*Free ressources:*/
+	for(i=0;j<M;i++){
+		IssmDouble* matrix=array[j];
+		xDelete<IssmDouble>(matrix);
+	}
+	xDelete<int>(mdims_array);
+	xDelete<int>(ndims_array);
+	xDelete<IssmDouble*>(array);
+
+	/*Assign output pointers: */
+	*pmass_flux=mass_flux;
+
+}/*}}}*/
+#endif
 #ifdef _HAVE_CONTROL_
 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/
@@ -601,5 +661,4 @@
 /*}}}*/
 #endif
-
 #ifdef  _HAVE_DAKOTA_
 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13720)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13721)
@@ -50,21 +50,21 @@
 		FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
 		~FemModel();
+
+		/*Methods:*/
+		void Echo();
 		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
-
-		/*Methods: {{{*/
-		void Echo();
 		void Solve(void);
 		void OutputResults(void);
 		void SetStaticComm();
-		/*}}}*/
-		/*Fem: {{{*/
-		void  SetCurrentConfiguration(int configuration_type);
-		void  SetCurrentConfiguration(int configuration_type,int analysis_type);
-		/*}}}*/
-		/*Modules: {{{*/
-		 #ifdef  _HAVE_DAKOTA_
+		void SetCurrentConfiguration(int configuration_type);
+		void SetCurrentConfiguration(int configuration_type,int analysis_type);
+
+		/*Modules*/ 
+		#ifdef _HAVE_RESPONSES_
+		void MassFluxx(IssmDouble* pmass_flux,bool process_units);
+		#endif
+		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
-
 		void RequestedOutputsx(int* requested_outputs, int numoutputs);
 		void RequestedDependentsx(void);
@@ -78,8 +78,4 @@
 		void UpdateConstraintsx(void);
 		int  UpdateVertexPositionsx(void);
-
-
-		/*}}}*/
-
 };
 
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 13720)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 13721)
@@ -63,5 +63,4 @@
 #include "./Krigingx/Krigingx.h"
 #include "./Shp2Kmlx/Shp2Kmlx.h"
-#include "./MassFluxx/MassFluxx.h"
 #include "./MaxAbsVxx/MaxAbsVxx.h"
 #include "./MaxAbsVyx/MaxAbsVyx.h"