Index: ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp	(revision 18604)
@@ -47,7 +47,7 @@
 
 	/* Intermediaries */
 	int         domaintype;
-	IssmDouble  Jdet,thickness,l=8.;
+	IssmDouble  Jdet,thickness,l;
 	IssmDouble *xyz_list = NULL;
 
 	/*Check dimension*/
@@ -67,6 +67,7 @@
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
+	element->FindParam(&l,SmoothThicknessMultiplierEnum); _assert_(l>0.);
 	element->GetVerticesCoordinates(&xyz_list);
 	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
 
Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18604)
@@ -194,27 +194,46 @@
 }/*}}}*/
 void AdjointBalancethickness2Analysis::GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
+	/*Intermediaries*/
+	IssmDouble lambda,Jdet; 
+	IssmDouble *xyz_list= NULL;
+
 	/*Fetch number of vertices for this finite element*/
 	int numvertices = element->GetNumberOfVertices();
 
 	/*Initialize some vectors*/
+	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
 	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
-	IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
 	int*        vertexpidlist = xNew<int>(numvertices);
 
 	/*Retrieve all inputs we will be needing: */
+	element->GetVerticesCoordinates(&xyz_list);
 	element->GradientIndexing(&vertexpidlist[0],control_index);
-	element->GetInputListOnVertices(lambda,AdjointEnum);
-	for(int i=0;i<numvertices;i++){
-		ge[i]= - lambda[i];
-		_assert_(!xIsNan<IssmDouble>(ge[i]));
+	Input* adjoint_input = element->GetInput(AdjointEnum);            _assert_(adjoint_input);
+
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsP1(basis,gauss);
+
+		adjoint_input->GetInputValue(&lambda,gauss);
+
+		/*Build gradient vector (actually -dJ/da): */
+		for(int i=0;i<numvertices;i++){
+			ge[i]+= -Jdet*gauss->weight*basis[i]*lambda;
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
 	}
-	gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(ge);
-	xDelete<IssmDouble>(lambda);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
 	xDelete<int>(vertexpidlist);
+	delete gauss;
 }/*}}}*/
 void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	element->InputUpdateFromSolutionOneDof(solution,AdjointEnum);
Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18604)
@@ -81,10 +81,10 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
 	Input* D0_input    = element->GetInput(BalancethicknessD0Enum);
-	if(!D0_input){
+	//if(!D0_input){
 		this->CreateD0(element);
 		D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
-	}
+	//}
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -219,7 +219,7 @@
 void Balancethickness2Analysis::CreateD0(Element* element){/*{{{*/
 
 	/*Intermediaries */
-	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b;
+	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b,normds;
 	const int        n = 3;
 	const IssmDouble Hstar = 500.;
 	const IssmDouble Lstar = 500.e+3;
@@ -231,9 +231,11 @@
 
 	/*retrieve what we need: */
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
-	Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
-	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+	Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum);
+	Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum);
+	Input* surface_input       = element->GetInput(SurfaceEnum);            _assert_(surface_input);
+	Input* B_input             = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+	IssmDouble rhog            = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
 	/*Calculate damage evolution source term: */
 	Gauss* gauss=element->NewGauss();
@@ -241,8 +243,25 @@
 		gauss->GaussNode(element->GetElementType(),i);
 		
 		B_input->GetInputValue(&B,gauss);
-		surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+		if(surfaceslopex_input && surfaceslopey_input){
+			surfaceslopex_input->GetInputValue(&ds[0],gauss);
+			surfaceslopey_input->GetInputValue(&ds[1],gauss);
+		}
+		else{
+			surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+		}
 
+		/*check slopes*/
+		normds = sqrt(ds[0]*ds[0]+ds[1]*ds[1]);
+		if (normds==0.){
+			_error_("surface slope is zero");
+		}
+		if(normds<1.e-5){
+			ds[0] = ds[0]/normds*1.e+5;
+			ds[1] = ds[1]/normds*1.e+5;
+			normds = 1.e-5;
+		}
+
 		mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
 		Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
 
@@ -251,6 +270,10 @@
 
 	/*Add input*/
 	element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType());
+	//if(surfaceslopex_input && surfaceslopey_input){
+	//	element->DeleteInput(SurfaceSlopeXEnum);
+	//	element->DeleteInput(SurfaceSlopeYEnum);
+	//}
 	
 	/*Clean up and return*/
 	xDelete<IssmDouble>(D0);
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18603)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18604)
@@ -449,6 +449,7 @@
 	InputToL2ProjectEnum,
 	InputToDepthaverageEnum,
 	InputToSmoothEnum,
+	SmoothThicknessMultiplierEnum,
 	IntParamEnum,
 	IntVecParamEnum,
 	TransientParamEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18604)
@@ -446,6 +446,7 @@
 		case InputToL2ProjectEnum : return "InputToL2Project";
 		case InputToDepthaverageEnum : return "InputToDepthaverage";
 		case InputToSmoothEnum : return "InputToSmooth";
+		case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier";
 		case IntParamEnum : return "IntParam";
 		case IntVecParamEnum : return "IntVecParam";
 		case TransientParamEnum : return "TransientParam";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18604)
@@ -455,6 +455,7 @@
 	      else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum;
 	      else if (strcmp(name,"InputToDepthaverage")==0) return InputToDepthaverageEnum;
 	      else if (strcmp(name,"InputToSmooth")==0) return InputToSmoothEnum;
+	      else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
 	      else if (strcmp(name,"IntParam")==0) return IntParamEnum;
 	      else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
 	      else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
@@ -504,11 +505,11 @@
 	      else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
 	      else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
 	      else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
-	      else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"Boundary")==0) return BoundaryEnum;
+	      if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
+	      else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
 	      else if (strcmp(name,"Converged")==0) return ConvergedEnum;
 	      else if (strcmp(name,"Fill")==0) return FillEnum;
 	      else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"Step")==0) return StepEnum;
 	      else if (strcmp(name,"Time")==0) return TimeEnum;
 	      else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
-	      else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
+	      if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
+	      else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
 	      else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
 	      else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
 	      else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum;
@@ -750,11 +751,11 @@
 	      else if (strcmp(name,"BasalforcingsOceanVx")==0) return BasalforcingsOceanVxEnum;
 	      else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
 	      else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
-	      else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
+	      if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
+	      else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
 	      else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
 	      else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;
 	      else if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18604)
@@ -28,7 +28,7 @@
 	iomodel->FetchDataToInput(elements,InversionVxObsEnum,0.);
 	iomodel->FetchDataToInput(elements,InversionVyObsEnum,0.); 
 	iomodel->FetchDataToInput(elements,InversionThicknessObsEnum,0.);
-	//iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.);
+	iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.);
 
 	iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum);
 
Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18604)
@@ -12,7 +12,7 @@
 
 	int         solution_type,n;
 	int         num_responses;
-	IssmDouble  j0,j;
+	IssmDouble  j0,j,yts;
 	IssmDouble  Ialpha,exponent,alpha;
 	IssmDouble* jlist = NULL;
 	IssmDouble *G = NULL;
@@ -27,6 +27,7 @@
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
 	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
 
 	/*Get initial guess*/
 	Vector<IssmDouble> *Xpetsc = NULL;
@@ -63,7 +64,7 @@
 	/*Range of tests*/
 	IssmDouble exp0 = 0.;
 	IssmDouble incr = -0.2;
-	IssmDouble exp1 = -10.;
+	IssmDouble exp1 = -18.;
 	int        num  = reCast<int,IssmDouble>((exp1-exp0)/incr);
 
 	/*Allocate output*/
@@ -87,7 +88,7 @@
 		femmodel->CostFunctionx(&j,NULL,NULL);
 
 		IssmDouble Den = 0.;
-		for(int i=0;i<n;i++) Den += alpha * G[i] * 1.;
+		for(int i=0;i<n;i++) Den += alpha* G[i] * 1.;
 		Ialpha = fabs((j - j0)/Den - 1.);
 
 		_printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
Index: ../trunk-jpl/src/c/cores/balancevelocity_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/balancevelocity_core.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/cores/balancevelocity_core.cpp	(revision 18604)
@@ -17,7 +17,8 @@
 	/*recover parameters: */
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-	if(VerboseSolution()) _printf0_("computing smoothed slopes:\n");
+	if(VerboseSolution()) _printf0_("computing smooth driving stress:\n");
+	femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum);
 	femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
 	femmodel->parameters->SetParam(DrivingStressXEnum,InputToSmoothEnum);
 	solutionsequence_linear(femmodel);
Index: ../trunk-jpl/src/c/cores/balancethickness2_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18603)
+++ ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18604)
@@ -14,20 +14,28 @@
 	/*parameters: */
 	bool save_results;
 
-	/*activate formulation: */
-	femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
-
 	/*recover parameters: */
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
+	if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n");
+	femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum);
+	femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
+	femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToSmoothEnum);
+	solutionsequence_linear(femmodel);
+	femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum);
+	solutionsequence_linear(femmodel);
+
 	if(VerboseSolution()) _printf0_("call computational core:\n");
+	femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
 	solutionsequence_linear(femmodel);
 	//solutionsequence_nonlinear(femmodel,false);
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		const int numoutputs = 4;
-		int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
+		const int numoutputs = 6;
+		int outputs[numoutputs] = {SurfaceEnum,SurfaceSlopeXEnum,SurfaceSlopeYEnum,VxEnum,VyEnum,VelEnum};
+		//const int numoutputs = 4;
+		//int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
 		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
 	}
 
Index: ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m	(revision 18603)
+++ ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m	(revision 18604)
@@ -1,11 +0,0 @@
-function macro=SmoothedSurfaceSlopeXAnalysisEnum()
-%SMOOTHEDSURFACESLOPEXANALYSISENUM - Enum of SmoothedSurfaceSlopeXAnalysis
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
-%            Please read src/c/shared/Enum/README for more information
-%
-%   Usage:
-%      macro=SmoothedSurfaceSlopeXAnalysisEnum()
-
-macro=StringToEnum('SmoothedSurfaceSlopeXAnalysis');
Index: ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m	(revision 18603)
+++ ../trunk-jpl/src/m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m	(revision 18604)
@@ -1,11 +0,0 @@
-function macro=SmoothedSurfaceSlopeYAnalysisEnum()
-%SMOOTHEDSURFACESLOPEYANALYSISENUM - Enum of SmoothedSurfaceSlopeYAnalysis
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
-%            Please read src/c/shared/Enum/README for more information
-%
-%   Usage:
-%      macro=SmoothedSurfaceSlopeYAnalysisEnum()
-
-macro=StringToEnum('SmoothedSurfaceSlopeYAnalysis');
Index: ../trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m	(revision 18604)
@@ -0,0 +1,11 @@
+function macro=SmoothThicknessMultiplierEnum()
+%SMOOTHTHICKNESSMULTIPLIERENUM - Enum of SmoothThicknessMultiplier
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SmoothThicknessMultiplierEnum()
+
+macro=StringToEnum('SmoothThicknessMultiplier');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18603)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18604)
@@ -438,6 +438,7 @@
 def InputToL2ProjectEnum(): return StringToEnum("InputToL2Project")[0]
 def InputToDepthaverageEnum(): return StringToEnum("InputToDepthaverage")[0]
 def InputToSmoothEnum(): return StringToEnum("InputToSmooth")[0]
+def SmoothThicknessMultiplierEnum(): return StringToEnum("SmoothThicknessMultiplier")[0]
 def IntParamEnum(): return StringToEnum("IntParam")[0]
 def IntVecParamEnum(): return StringToEnum("IntVecParam")[0]
 def TransientParamEnum(): return StringToEnum("TransientParam")[0]
