Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23950)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23951)
@@ -112,4 +112,18 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
 
+	/*Deal with friction parameters*/
+	int frictionlaw;
+	iomodel->FindConstant(&frictionlaw,"md.friction.law");
+	if(frictionlaw==4 || frictionlaw==6){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+	}
+	if(frictionlaw==3 || frictionlaw==1 || frictionlaw==7){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+	}
+	if(frictionlaw==9){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+		parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
+	}
+
   /*Requested outputs*/
   iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
@@ -153,7 +167,7 @@
 	/*Get all inputs and parameters*/
 	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
-	IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
 	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
+	IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
 	Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
 	Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
@@ -179,8 +193,7 @@
 		IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);
 
+		/*Diffusive term*/
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
-
-				/*Diffusive term*/
 				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
 							coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
@@ -189,7 +202,14 @@
 		}
 
-		/*Transient term if dt*/
+		/*Transient term if dt>0*/
 		if(dt>0.){
+			/*Diffusive term*/
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*basis[i]*basis[j];
+				}
+			}
 		}
+
 	}
 
@@ -209,5 +229,5 @@
 	IssmDouble  Jdet,w,v,vx,vy,ub,h,N,h_r;
 	IssmDouble  G,m,frictionheat,alpha2;
-	IssmDouble  A,B,n;
+	IssmDouble  A,B,n,phi_old;
 	IssmDouble* xyz_list = NULL;
 
@@ -221,15 +241,20 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	IssmDouble  L       = element->FindParam(MaterialsLatentheatEnum);
-	IssmDouble  rho_ice = element->FindParam(MaterialsRhoIceEnum);
-	IssmDouble l_r      = element->FindParam(HydrologyCavitySpacingEnum);
-	Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
-	Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
-	Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
-	Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
-	Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
-	Input*  G_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
-	Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
-	Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
+	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+	IssmDouble l_r       = element->FindParam(HydrologyCavitySpacingEnum);
+	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+	IssmDouble g         = element->FindParam(ConstantsGEnum);
+	IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
+	Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
+	Input* vx_input     = element->GetInput(VxEnum);_assert_(vx_input);
+	Input* vy_input     = element->GetInput(VyEnum);_assert_(vy_input);
+	Input* N_input      = element->GetInput(EffectivePressureEnum); _assert_(N_input);
+	Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+	Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
+	Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
 
 	/*Build friction element, needed later: */
@@ -273,4 +298,10 @@
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i];
+
+		/*Transient term if dt>0*/
+		if(dt>0.){
+			phiold_input->GetInputValue(&phi_old,gauss);
+			for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
+		}
 	}