Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 25941)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 25942)
@@ -5,4 +5,6 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+//#define INWOOVZ 1 //For Inwoo Park 0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID
 
 /*Model processing*/
@@ -112,5 +114,7 @@
 	}
 	iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
+	#ifdef INWOOVZ
 	iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
+	#endif
 
 
@@ -154,6 +158,4 @@
 
 	/* specific parameters*/
-	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvzsolution",StressbalanceIsVzSolutionEnum));
-	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvelbc",StressbalanceIsVelBCEnum));
 
 }/*}}}*/
@@ -518,7 +520,4 @@
 	IssmDouble   rho_ice,g;
 
-	int          isvzsolution;
-	IssmDouble   isvelbc, rheology_n, eta;
-
 	/*Get the approximation and do nothing if the element in FS or None*/
 	element->GetInputValue(&approximation,ApproximationEnum);
@@ -526,8 +525,4 @@
 		return;
 	}
-
-	/*Get parameters*/
-	element->FindParam(&isvzsolution,StressbalanceIsVzSolutionEnum);
-	element->FindParam(&isvelbc,StressbalanceIsVelBCEnum);
 
 	/*Get dof list and vertices coordinates: */
@@ -544,7 +539,4 @@
 	IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
 	IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
-	IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
-	IssmDouble*  base	     = xNew<IssmDouble>(numnodes);
-	IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector vz: */
@@ -588,17 +580,23 @@
 	}
 
+
+	#ifdef INWOOVZ
+	IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
+	IssmDouble*  base	     = xNew<IssmDouble>(numnodes);
+	IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
+
+	IssmDouble rheology_n=element->material->GetN();
+	IssmDouble isvelbc   = 10./(365.25*24*3600); /*10 m/yr*/
 	/*Set analytical vertical velocity field at slow flow region.*/
-	if (isvzsolution == 1){
-		rheology_n=element->material->GetN();
+	if(INWOOVZ == 1){
 		element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
 		element->GetInputListOnNodes(&base[0],BaseEnum,0.);
 		element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
 		for(i=0;i<numnodes;i++){
-			eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+			IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
 			vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
 		}	
 	}
-	else if (isvzsolution == 2) {
-		rheology_n=element->material->GetN();
+	else if(INWOOVZ== 2) {
 		element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
 		element->GetInputListOnNodes(&base[0],BaseEnum,0.);
@@ -606,9 +604,21 @@
 		for(i=0;i<numnodes;i++){
 			if ( vx[i]*vx[i]+vy[i]*vy[i] < isvelbc*isvelbc ){
-				eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+				IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
 				vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
 			}
 		}	
 	}
+	else if(INWOOVZ==0){
+		/*nothing to add*/
+	}
+	else{
+		_error_("not supported yet");
+	}
+
+	/*Cleanup*/
+	xDelete<IssmDouble>(smb);
+	xDelete<IssmDouble>(base);
+	xDelete<IssmDouble>(thickness);
+	#endif
 
 
@@ -637,7 +647,4 @@
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(smb);
-	xDelete<IssmDouble>(base);
-	xDelete<IssmDouble>(thickness);
 	xDelete<IssmDouble>(surface);
 	xDelete<IssmDouble>(pressure);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25941)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25942)
@@ -3015,10 +3015,19 @@
 	IssmDouble* values    = xNew<IssmDouble>(numnodes);
 
+		if(this->Id()==104 || this->Id()==266 || this->Id()==264 || this->Id()==247 || this->Id()==263 || this->Id()==249 || this->Id()==250){
+			if(enum_type==MaskIceLevelsetEnum) printf("-----------------\n");
+		}
+
 	/*Use the dof list to index into the solution vector: */
 	for(int i=0;i<numnodes;i++){
 		values[i]=solution[doflist[i]];
+		if(this->Id()==104 || this->Id()==266 || this->Id()==264 || this->Id()==247 || this->Id()==263 || this->Id()==249 || this->Id()==250){
+			if(enum_type==MaskIceLevelsetEnum) printf("%g\n",values[i]);
+			if(enum_type==MaskIceLevelsetEnum && values[i]>0.5 && values[i]<1.5 )this->nodes[i]->DeepEcho();
+		}
 		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
 		if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector, SID = " << this->Sid());
 	}
+
 
 	/*Add input to the element: */
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 25941)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 25942)
@@ -1399,4 +1399,5 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1405,6 +1406,6 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1419,6 +1420,6 @@
 syn keyword cType ElementInput
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1427,9 +1428,10 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1446,4 +1448,5 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1456,4 +1459,5 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1468,6 +1472,6 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1483,11 +1487,11 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType SealevelMasks
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1508,8 +1512,4 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis
@@ -1530,6 +1530,6 @@
 syn keyword cType FreeSurfaceBaseAnalysis
 syn keyword cType FreeSurfaceTopAnalysis
+syn keyword cType GiaAnalysis
 syn keyword cType GLheightadvectionAnalysis
-syn keyword cType GiaAnalysis
 syn keyword cType HydrologyDCEfficientAnalysis
 syn keyword cType HydrologyDCInefficientAnalysis
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25941)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25942)
@@ -441,6 +441,4 @@
 	StressbalanceRiftPenaltyThresholdEnum,
 	StressbalanceShelfDampeningEnum,
-	StressbalanceIsVzSolutionEnum,
-	StressbalanceIsVelBCEnum,
 	ThermalIsdrainicecolumnEnum,
 	ThermalIsdynamicbasalspcEnum,
Index: /issm/trunk-jpl/src/m/classes/stressbalance.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 25941)
+++ /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 25942)
@@ -23,6 +23,4 @@
 		loadingforce           = NaN;
 		requested_outputs      = {};
-		isvzsolution           = NaN;
-		isvelbc                = NaN;	
 	end
 	methods
@@ -77,8 +75,4 @@
 			 self.rift_penalty_lock=10;
 
-			 %Set vertical velocity option
-			 self.isvzsolution = 0;
-			 self.isvelbc      = 10.;
-
 			 %output default:
 			 self.requested_outputs={'default'};
@@ -102,6 +96,4 @@
 			md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices 6]);
 			md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices 3]);
-			md = checkfield(md,'fieldname','stressbalance.isvzsolution','numel',[1],'values',[0 1 2]);
-			md = checkfield(md,'fieldname','stressbalance.isvelbc','numel',[1],'>=',0.0);
 			md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
 			if ~any(isnan(md.stressbalance.vertex_pairing)),
@@ -166,8 +158,4 @@
 			fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
 			fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized');
-
-			disp(sprintf('\n      %s','Estimation of vertical velocity field options:'));
-			fielddisplay(self,'isvzsolution','0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID');
-			fielddisplay(self,'isvelbc','If turned on isvzsolution as 2, region (vel < isvelbc) calculates for ID.');
 
 			disp(sprintf('\n      %s','Other:'));
@@ -199,6 +187,4 @@
 			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
 			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvzsolution','format','Integer','numel',1,'values',[0 1 2]);
-			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvelbc','format','Double','numel',1,'scale',1./yts);
 
 			if size(self.loadingforce,2)==3,