Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16775)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16776)
@@ -26,6 +26,7 @@
 		void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16775)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16776)
@@ -952,7 +952,7 @@
 			InputUpdateFromSolutionHOFS(solution,element);
 			return;
 		case SSAFSApproximationEnum:
-			_error_("not implemented yet");
+			InputUpdateFromSolutionSSAFS(solution,element);
 			return;
 		default:
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
@@ -1423,6 +1423,112 @@
 	xDelete<int>(doflist);
 	if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
 }/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i;
+	IssmDouble  rho_ice,g,FSreconditioning;
+	int*        doflistSSA  = NULL;
+	int*        doflistFSv = NULL;
+	int*        doflistFSp = NULL;
+
+	/*we have to add results of this element for FS and results from the element
+	 * at base for SSA, so we need to have the pointer toward the basal element*/
+	Element* basalelement=element->GetBasalElement();
+	if(basalelement->ObjectEnum()!=PentaEnum){
+		_error_("Coupling not supported for "<<EnumToStringx(basalelement->ObjectEnum()));
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes  = 6;
+	int numdof2d  = numnodes;
+	int numdofSSA = 6*2;
+	int numdofFSv = 6*3;
+	int numdofFSp = 6;
+
+	/*Fetch dof list and allocate solution vectors*/
+	element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+	element->GetDofListPressure(&doflistFSp,GsetEnum);
+	basalelement->GetDofList(&doflistSSA, SSAApproximationEnum, GsetEnum);
+	IssmDouble* SSAvalues  = xNew<IssmDouble>(numdofSSA);
+	IssmDouble* FSvalues  = xNew<IssmDouble>(numdofFSv+numdofFSp);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vzSSA      = xNew<IssmDouble>(numnodes);
+	IssmDouble* vzFS      = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(2*numnodes);
+	for(i=0;i<numnodes;i++) cs_list[i]          = XYZEnum;
+	for(i=0;i<numnodes;i++) cs_list[numnodes+i] = PressureEnum;
+
+	/*Use the dof list to index into the solution vector: */
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	for(i=0;i<numdof2d;i++){
+		SSAvalues[i]          = solution[doflistSSA[i]];
+		SSAvalues[i+numdof2d] = solution[doflistSSA[i]];
+	}
+	for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
+	for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
+
+	/*Transform solution in Cartesian Space*/
+	element->TransformSolutionCoord(FSvalues,2*numnodes,cs_list);
+	element->TransformSolutionCoord(SSAvalues,numnodes,XYEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+
+	for(i=0;i<numnodes;i++){
+		vx[i]       = FSvalues[i*3+0]+SSAvalues[i*2+0];
+		vy[i]       = FSvalues[i*3+1]+SSAvalues[i*2+1];
+		vzFS[i]     = FSvalues[i*3+2];
+		pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get Vz and compute vel*/
+	element->GetInputListOnVertices(vzSSA,VzSSAEnum);
+	for(i=0;i<numnodes;i++){
+		vz[i] = vzSSA[i]+vzFS[i];
+		vel[i]= sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+	}
+
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+	element->InputChangeName(VzEnum,VzPicardEnum);
+	element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+	/*Add vx and vy as inputs to element: */
+	element->AddInput(VxEnum,vx,P1Enum);
+	element->AddInput(VyEnum,vy,P1Enum);
+	element->AddInput(VzEnum,vz,P1Enum);
+	element->AddInput(VzFSEnum,vzFS,P1Enum);
+	element->AddInput(VelEnum,vel,P1Enum);
+	element->AddInput(PressureEnum,pressure,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vzSSA);
+	xDelete<IssmDouble>(vzFS);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(FSvalues);
+	xDelete<IssmDouble>(SSAvalues);
+	xDelete<int>(doflistFSp);
+	xDelete<int>(doflistFSv);
+	xDelete<int>(doflistSSA);
+	xDelete<int>(cs_list);
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,meshtype;