Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18550)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18551)
@@ -33,6 +33,7 @@
 	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
 	iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
+	iomodel->FetchDataToInput(elements,BalancethicknessCmuEnum);
 }/*}}}*/
 void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
@@ -143,7 +144,7 @@
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
 					(ms-mb-dhdt)*basis[i]
-					-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
+					//-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
 					);
 	}
 
@@ -154,13 +155,63 @@
 	return pe;
 }/*}}}*/
 void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-	   _error_("not implemented yet");
+		element->GetSolutionFromInputsOneDof(solution,SurfaceEnum);
 }/*}}}*/
 void Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
 void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+
+	/*Intermediaries*/
+	IssmDouble  ds[2],s,b,D;
+	IssmDouble* xyz_list = NULL;
+
+	//element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+	element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
+
+	/*Fetch number of vertices and allocate velocity vectors*/
+	int numvertices = element->GetNumberOfVertices();
+	IssmDouble* vel_list = xNew<IssmDouble>(numvertices);
+	IssmDouble* vx_list  = xNew<IssmDouble>(numvertices);
+	IssmDouble* vy_list  = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);
+	Input* H_input   = element->GetInput(ThicknessEnum);                            _assert_(H_input);
+	Input* s_input   = element->GetInput(SurfaceEnum);                              _assert_(s_input);
+	Input* b_input   = element->GetInput(BaseEnum);                                 _assert_(b_input);
+
+	/*Calculate velocities*/
+	Gauss* gauss=element->NewGauss();
+	for(int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		if(D_input){
+			D_input->GetInputValue(&D,gauss);
+		}
+		else{
+			D = 0.;
+		}
+		b_input->GetInputValue(&b,gauss);
+		s_input->GetInputValue(&s,gauss);
+		s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+
+		vx_list[iv] = -1./(s-b)*D*ds[0];
+		vy_list[iv] = -1./(s-b)*D*ds[1];
+		vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2));
+	}
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddInput(VxEnum,vx_list,P1Enum);
+	element->AddInput(VyEnum,vy_list,P1Enum);
+	element->AddInput(VelEnum,vel_list,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(vy_list);
+	xDelete<IssmDouble>(vx_list);
+	xDelete<IssmDouble>(vel_list);
+	xDelete<IssmDouble>(xyz_list);
 }/*}}}*/
 void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
@@ -171,7 +222,7 @@
 void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/
 
 	/*Intermediaries */
-	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k;
+	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b;
 	const int        n = 3;
 	const IssmDouble Hstar = 500.;
 	const IssmDouble Lstar = 500.e+3;
@@ -184,10 +235,11 @@
 
 	/*retrieve what we need: */
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* thickness_input  = element->GetInput(ThicknessEnum);          _assert_(thickness_input);
+	Input* bed_input        = element->GetInput(BaseEnum);               _assert_(bed_input);
 	Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
 	Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
 	Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
+	Input* Cmu_input        = element->GetInput(BalancethicknessCmuEnum);_assert_(Cmu_input);
 	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
 	/*Calculate damage evolution source term: */
@@ -197,14 +249,18 @@
 		
 		B_input->GetInputValue(&B,gauss);
 		k_input->GetInputValue(&k,gauss);
-		thickness_input->GetInputValue(&h,gauss);
+		//thickness_input->GetInputValue(&h,gauss);
 		surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+		surface_input->GetInputValue(&s,gauss);
+		bed_input->GetInputValue(&b,gauss);
+		h = s-b;
 
 		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);
-		Cmu   = mu0/Hstar/(k*k);
+		//Cmu   = pow(mu0,n)/Hstar/(k*k);
+		Cmu_input->GetInputValue(&Cmu,gauss);
 
-		D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu);
+		D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h + Cmu);
 	}
 
 	/*Add input*/
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18550)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18551)
@@ -304,6 +304,7 @@
 	BalancethicknessApparentMassbalanceEnum,
 	Balancethickness2MisfitEnum,
 	BalancethicknessDiffusionCoefficientEnum,
+	BalancethicknessCmuEnum,
 	/*}}}*/
 	/*Surfaceforcings{{{*/
 	SurfaceforcingsEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18550)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18551)
@@ -312,6 +312,7 @@
 		case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
 		case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
 		case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient";
+		case BalancethicknessCmuEnum : return "BalancethicknessCmu";
 		case SurfaceforcingsEnum : return "Surfaceforcings";
 		case SMBEnum : return "SMB";
 		case SurfaceforcingsMassBalanceEnum : return "SurfaceforcingsMassBalance";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18550)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18551)
@@ -318,6 +318,7 @@
 	      else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
 	      else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
 	      else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
+	      else if (strcmp(name,"BalancethicknessCmu")==0) return BalancethicknessCmuEnum;
 	      else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
 	      else if (strcmp(name,"SMB")==0) return SMBEnum;
 	      else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
@@ -381,11 +382,11 @@
 	      else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
 	      else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
 	      else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
-	      else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
+	      if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
 	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
 	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
 	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
@@ -504,11 +505,11 @@
 	      else if (strcmp(name,"Fill")==0) return FillEnum;
 	      else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
 	      else if (strcmp(name,"Friction")==0) return FrictionEnum;
-	      else if (strcmp(name,"Internal")==0) return InternalEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
+	      if (strcmp(name,"Internal")==0) return InternalEnum;
+	      else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
 	      else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
 	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
 	      else if (strcmp(name,"Pressure")==0) return PressureEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
 	      else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
 	      else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
-	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
+	      if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
 	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
@@ -750,11 +751,11 @@
 	      else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
 	      else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
 	      else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
-	      else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
+	      if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
+	      else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
 	      else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
 	      else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
 	      else if (strcmp(name,"MaterialsCompressionCoef")==0) return MaterialsCompressionCoefEnum;
Index: ../trunk-jpl/src/c/cores/balancethickness2_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18550)
+++ ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18551)
@@ -21,12 +21,15 @@
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
 	if(VerboseSolution()) _printf0_("call computational core:\n");
-	solutionsequence_linear(femmodel);
+	//solutionsequence_linear(femmodel);
+	solutionsequence_nonlinear(femmodel,false);
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		int outputs[1] = {ThicknessEnum};
-		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
+		const int numoutputs = 5;
+		//int outputs[numoutputs] = {ThicknessEnum};
+		int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum,BalancethicknessDiffusionCoefficientEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
 	}
 
 }
Index: ../trunk-jpl/src/m/plot/plot_googlemaps.m
===================================================================
--- ../trunk-jpl/src/m/plot/plot_googlemaps.m	(revision 18550)
+++ ../trunk-jpl/src/m/plot/plot_googlemaps.m	(revision 18551)
@@ -35,6 +35,8 @@
 			[linspace(ylim(1),ylim(1),100) linspace(ylim(1),ylim(2),100) linspace(ylim(2),ylim(2),100) linspace(ylim(2),ylim(1),100)],...
 			-1,0,71);
 	else
+		latlist = md.mesh.lat; %That might work?
+		lonlist = md.mesh.long;
 		error('field hemisphere should either be ''n'' or ''s''');
 	end
 
Index: ../trunk-jpl/src/m/plot/googlemaps.m
===================================================================
--- ../trunk-jpl/src/m/plot/googlemaps.m	(revision 18550)
+++ ../trunk-jpl/src/m/plot/googlemaps.m	(revision 18551)
@@ -53,6 +53,8 @@
 	EPSGlocal = 'EPSG:3413'; % UPS Greenland  http://www.spatialreference.org/ref/epsg/3413/
 elseif strcmpi(md.mesh.hemisphere,'s'),
 	EPSGlocal = 'EPSG:3031'; % UPS Antarctica http://www.spatialreference.org/ref/epsg/3031/
+elseif strcmpi(md.mesh.hemisphere,'UTM18N'),
+	EPSGlocal = 'EPSG:32218'; % UTM 18 N Patagonia
 else
 	error('field hemisphere should either be ''n'' or ''s''');
 end
Index: ../trunk-jpl/src/m/classes/balancethickness.m
===================================================================
--- ../trunk-jpl/src/m/classes/balancethickness.m	(revision 18550)
+++ ../trunk-jpl/src/m/classes/balancethickness.m	(revision 18551)
@@ -8,22 +8,10 @@
 		spcthickness      = NaN;
 		thickening_rate   = NaN;
 		stabilization     = 0;
+
+		Cmu               = NaN;
 	end
 	methods
-        function createxml(obj,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- balancethickness -->');
-            
-			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcthickness" type="',   	class(obj.spcthickness),'" default="',        convert2str(obj.spcthickness),'">',	'     <section name="balancethickness" />','     <help> thickness constraints (NaN means no constraint) [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="thickening_rate" type="', 	class(obj.thickening_rate),'" default="',            convert2str(obj.thickening_rate),'">',     '     <section name="balancethickness" />','     <help> ice thickening rate used in the mass conservation (dh/dt) [m/yr] </help>','</parameter>');
-            
-            % balancethickness drop-down (1,2, or 3)
-            fprintf(fid,'%s\n%s\n%s\n',        '<parameter key ="stabilization" type="alternative" optional="false">','     <section name="balancethickness" />','     <help> 0: None, 1: SU, 2: SSAs artificial diffusivity, 3:DG </help>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="2" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="3" type="string" default="false"> </option>','</parameter>');
-
-        end % }}}
 		function obj = balancethickness(varargin) % {{{
 			switch nargin
 				case 0
@@ -45,6 +33,8 @@
 			md = checkfield(md,'fieldname','balancethickness.spcthickness');
 			md = checkfield(md,'fieldname','balancethickness.thickening_rate','size',[md.mesh.numberofvertices 1],'NaN',1);
 			md = checkfield(md,'fieldname','balancethickness.stabilization','size',[1 1],'values',[0 1 2 3]);
+
+			md = checkfield(md,'fieldname','balancethickness.Cmu','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0);
 		end % }}}
 		function disp(obj) % {{{
 			disp(sprintf('   balance thickness solution parameters:'));
@@ -61,6 +51,8 @@
 			WriteData(fid,'object',obj,'fieldname','spcthickness','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','thickening_rate','format','DoubleMat','mattype',1,'scale',1./yts);
 			WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
+
+			WriteData(fid,'object',obj,'fieldname','Cmu','format','DoubleMat','mattype',1);
 		end % }}}
 	end
 end
Index: ../trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BalancethicknessCmuEnum.m	(revision 18551)
@@ -0,0 +1,11 @@
+function macro=BalancethicknessCmuEnum()
+%BALANCETHICKNESSCMUENUM - Enum of BalancethicknessCmu
+%
+%   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=BalancethicknessCmuEnum()
+
+macro=StringToEnum('BalancethicknessCmu');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18550)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18551)
@@ -304,6 +304,7 @@
 def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
 def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
 def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0]
+def BalancethicknessCmuEnum(): return StringToEnum("BalancethicknessCmu")[0]
 def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
 def SMBEnum(): return StringToEnum("SMB")[0]
 def SurfaceforcingsMassBalanceEnum(): return StringToEnum("SurfaceforcingsMassBalance")[0]
