Index: issm/trunk/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 25836)
+++ issm/trunk/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 26744)
@@ -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*/
@@ -105,5 +107,5 @@
 	iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
 	iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
-	iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
 	iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
 	if(iomodel->domaintype!=Domain2DhorizontalEnum){
@@ -112,13 +114,20 @@
 	}
 	iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
-	//iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
-
+	#ifdef INWOOVZ
+	iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
+	#endif
 
 	/*Add basal forcings to compute melt rate*/
+	bool isstochastic;
 	int basalforcing_model;
 	iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
+	iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
 	switch(basalforcing_model){
 		case FloatingMeltRateEnum:
 			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
+			if(isstochastic){
+            iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum);
+         }
 			break;
 		case LinearFloatingMeltRateEnum:
@@ -153,5 +162,5 @@
 void StressbalanceVerticalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
-	/*No specific parameters*/
+	/* specific parameters*/
 
 }/*}}}*/
@@ -163,4 +172,7 @@
 		femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
 		solutionsequence_linear(femmodel);
+}/*}}}*/
+void           StressbalanceVerticalAnalysis::PreCore(FemModel* femmodel){/*{{{*/
+	_error_("not implemented");
 }/*}}}*/
 ElementVector* StressbalanceVerticalAnalysis::CreateDVector(Element* element){/*{{{*/
@@ -576,4 +588,45 @@
 	}
 
+	#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(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++){
+			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== 2) {
+		element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+		element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+		element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+		for(i=0;i<numnodes;i++){
+			if ( vx[i]*vx[i]+vy[i]*vy[i] < isvelbc*isvelbc ){
+				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
+
 	/*Now Compute vel*/
 	for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
