Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22572)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22573)
@@ -873,5 +873,4 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.FSreconditioning",StressbalanceFSreconditioningEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.shelf_dampening",StressbalanceShelfDampeningEnum));
-	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.viscosity_overshoot",StressbalanceViscosityOvershootEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
 
@@ -1443,6 +1442,5 @@
 	/*Intermediaries*/
 	int         dim,domaintype;
-	IssmDouble  viscosity,newviscosity,oldviscosity;
-	IssmDouble  viscosity_overshoot,thickness,Jdet;
+	IssmDouble  viscosity,thickness,Jdet;
 	IssmDouble *xyz_list = NULL;
 
@@ -1468,12 +1466,8 @@
 	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
-	Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
 	Input* vy_input    = NULL;
-	Input* vyold_input = NULL;
 	if(dim==2){
 		vy_input    = element->GetInput(VyEnum);       _assert_(vy_input);
-		vyold_input = element->GetInput(VyPicardEnum); _assert_(vyold_input);
-	}
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -1487,20 +1481,18 @@
 		thickness_input->GetInputValue(&thickness, gauss);
 		element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
-		element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 
 		if(dim==2){
 			for(int i=0;i<numnodes;i++){
 				for(int j=0;j<numnodes;j++){
-					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*thickness*(
+					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*thickness*(
 								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
 								);
-					Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*thickness*(
+					Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*thickness*(
 								2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
 								);
-					Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*thickness*(
+					Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*thickness*(
 								2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
 								);
-					Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*thickness*(
+					Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*thickness*(
 								dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
 								);
@@ -1511,5 +1503,5 @@
 			for(int i=0;i<numnodes;i++){
 				for(int j=0;j<numnodes;j++){
-					Ke->values[i*numnodes+j] += gauss->weight*Jdet*newviscosity*thickness*(
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*viscosity*thickness*(
 								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
 								);
@@ -1908,9 +1900,4 @@
 	}
 
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	element->InputChangeName(VxEnum,VxPicardEnum);
-	if(dim==2)element->InputChangeName(VyEnum,VyPicardEnum);
-
 	/*Add vx and vy as inputs to the tria element: */
 	element->AddBasalInput(VxEnum,vx,element->GetElementType());
@@ -2236,9 +2223,4 @@
 	for(i=0;i<numnodes;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);
-
 	/*Add vx and vy as inputs to the tria element: */
 	element->AddBasalInput(VxEnum,vx,element->GetElementType());
@@ -2428,6 +2410,5 @@
 	/*Intermediaries*/
 	int         dim,bsize;
-	IssmDouble  viscosity,newviscosity,oldviscosity;
-	IssmDouble  viscosity_overshoot,thickness,Jdet;
+	IssmDouble  viscosity,thickness,Jdet;
 	IssmDouble *xyz_list = NULL;
 
@@ -2446,12 +2427,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	Input* vx_input    = element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input);
 	Input* vy_input    = NULL;
-	Input* vyold_input = NULL;
 	if(dim==3){
 		vy_input=element->GetInput(VyEnum);          _assert_(vy_input);
-		vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
-	}
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -2463,21 +2440,18 @@
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
-		element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-
-		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 
 		if(dim==3){
 			for(int i=0;i<numnodes;i++){
 				for(int j=0;j<numnodes;j++){
-					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
 								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
 								);
-					Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+					Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
 								2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
 								);
-					Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+					Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
 								2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
 								);
-					Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+					Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
 								dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
 								);
@@ -2488,5 +2462,5 @@
 			for(int i=0;i<numnodes;i++){
 				for(int j=0;j<numnodes;j++){
-					Ke->values[i*numnodes+j] += gauss->weight*Jdet*newviscosity*(
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*viscosity*(
 								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
 								);
@@ -2905,9 +2879,4 @@
 		for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[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);
-	if(dim==3)element->InputChangeName(VyEnum,VyPicardEnum);
 
 	/*Add vx and vy as inputs to the element: */
@@ -5085,15 +5054,8 @@
 	else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
 
-	/*Now, we have to move the previous inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	element->InputChangeName(VxEnum,VxPicardEnum);
-	element->InputChangeName(VyEnum,VyPicardEnum);
-	if(pnumdof>0) element->InputChangeName(PressureEnum,PressurePicardEnum);
-	if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
-
 	/*Add vx and vy as inputs to the tria element: */
-	element->AddInput(VxEnum,      vx,      element->VelocityInterpolation());
-	element->AddInput(VyEnum,      vy,      element->VelocityInterpolation());
-	element->AddInput(VelEnum,     vel,     element->VelocityInterpolation());
+	element->AddInput(VxEnum, vx, element->VelocityInterpolation());
+	element->AddInput(VyEnum, vy, element->VelocityInterpolation());
+	element->AddInput(VelEnum,vel,element->VelocityInterpolation());
 	if(pnumdof>0) element->AddInput(PressureEnum,pressure,element->PressureInterpolation());
 	if(dim==3) element->AddInput(VzEnum,vz, element->VelocityInterpolation());
@@ -5879,5 +5841,5 @@
 	/*Intermediaries */
 	int         i,j;
-	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
+	IssmDouble  Jdet,viscosity;
 	IssmDouble  *xyz_list      = NULL;
 	IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
@@ -5909,9 +5871,6 @@
 	/* Get node coordinates and dof list: */
 	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
 	Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
-	Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
-	Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -5927,8 +5886,6 @@
 		this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 
 		element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
-		element->material->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
-
-		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2*newviscosity*gauss->weight*Jdet;
+
+		D_scalar=2*viscosity*gauss->weight*Jdet;
 		for (i=0;i<3;i++) D[i][i]=D_scalar;
 
@@ -6051,5 +6008,5 @@
 	int         i,j,approximation;
 	int         dim=3;
-	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
+	IssmDouble  Jdet,viscosity;
 	IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
 	IssmDouble  epsilons[6];                    //6 for FS
@@ -6071,9 +6028,6 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
 	Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
-	Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
-	Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
 	Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
 
@@ -6092,13 +6046,11 @@
 		if(approximation==SSAHOApproximationEnum){
 			element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
-			element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-			newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 		}
 		else if (approximation==SSAFSApproximationEnum){
-			element->material->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+			element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
 		}
 		else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
 
-		D_scalar=2*newviscosity*gauss->weight*Jdet;
+		D_scalar=2*viscosity*gauss->weight*Jdet;
 		for (i=0;i<3;i++) D[i][i]=D_scalar;
 
@@ -7018,11 +6970,4 @@
 	}
 
-	/*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);
@@ -7128,11 +7073,4 @@
 	}
 
-	/*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);
@@ -7228,10 +7166,4 @@
 	for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
 
-	/*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(PressureEnum,PressurePicardEnum);
-
 	/*Add vx and vy as inputs to element: */
 	element->AddInput(VxEnum,vx,P1Enum);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22572)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22573)
@@ -545,5 +545,4 @@
 	int found=-1;
 	for(int i=0;i<nummodels;i++){
-	
 		if (analysis_type_list[i]==configuration_type){
 			found=i;
@@ -566,9 +565,8 @@
 	if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){
 		ToolkitsOptionsFromAnalysis(this->parameters,analysis_type);
-		if(VerboseSolver()) _printf0_("      toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n");
-	}
-
-}
-/*}}}*/
+		if(VerboseSolver()) _printf0_("      toolkits Options set for analysis: " << EnumToStringx(analysis_type) << "\n");
+	}
+
+}/*}}}*/
 void FemModel::SetCurrentConfiguration(int configuration_type){/*{{{*/
 	this->SetCurrentConfiguration(configuration_type,configuration_type);