Index: issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 26463)
+++ issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 26490)
@@ -181,4 +181,5 @@
 	IssmDouble Jdet;
 	IssmDouble* xyz_list = NULL;
+	IssmDouble  gap,bed,thickness,head,g,rho_ice,rho_water,A,B,n,head_old;/***/
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -201,4 +202,16 @@
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 
+        /***Get all inputs and parameters*/
+        element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum);
+        element->FindParam(&rho_ice,MaterialsRhoIceEnum);
+        element->FindParam(&g,ConstantsGEnum);
+        Input* B_input = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        Input* gap_input = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
+        Input* thickness_input = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
+        Input* head_input = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
+        Input* base_input = element->GetInput(BaseEnum);                      _assert_(base_input);
+	Input* head_old_input = element->GetInput(HydrologyHeadOldEnum);	_assert_(head_old_input);
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(1);
@@ -209,8 +222,31 @@
 		element->NodalFunctions(basis,gauss);
 
+		/***/
+                base_input->GetInputValue(&bed,gauss);
+                thickness_input->GetInputValue(&thickness,gauss);
+                gap_input->GetInputValue(&gap,gauss);
+                head_input->GetInputValue(&head,gauss);
+  		head_old_input->GetInputValue(&head_old,gauss);
+
+                /***Get ice A parameter*/
+                B_input->GetInputValue(&B,gauss);
+                n_input->GetInputValue(&n,gauss);
+                A=pow(B,-n);
+
+                /***Get water and ice pressures*/
+                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
+                IssmDouble pressure_water = rho_water*g*(head-bed);
+                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
+
+                /***Get water pressure from previous time step to use in lagged creep term*/
+                IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
+                if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
+
+
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
 				Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
-				  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
+				  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]
+                			+gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j];
 			}
 		}
@@ -232,8 +268,9 @@
 	IssmDouble  Jdet,meltrate,G,dh[2],B,A,n;
 	IssmDouble  gap,bed,thickness,head,ieb,head_old;
-	IssmDouble  lr,br,vx,vy,beta;
+	IssmDouble  lr,br,vx,vy,beta,lc;
 	IssmDouble  alpha2,frictionheat;
    IssmDouble  PMPheat,dpressure_water[2],dbed[2];	
 	IssmDouble* xyz_list = NULL;
+        IssmDouble dgap[2]; /***/
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -261,4 +298,6 @@
 	Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
+//	Input* dgapxx_input	    = element->GetInput(HydrologyGapXXEnum); /***/
+//	Input* dgapyy_input	    = element->GetInput(HydrologyGapYYEnum); /***/
 
 	/*Get conductivity from inputs*/
@@ -289,5 +328,5 @@
 		lr_input->GetInputValue(&lr,gauss);
 		br_input->GetInputValue(&br,gauss);
-      headold_input->GetInputValue(&head_old,gauss);
+                headold_input->GetInputValue(&head_old,gauss);
 
 		/*Get ice A parameter*/
@@ -302,4 +341,12 @@
 		 beta = 0.;
 
+		/* TEST Compute lc term */
+//		if(gap<br)
+//		 lc = gap*(1-(br-gap)/br);
+//		else
+		 lc = gap; 
+//                lc=gap*tanh(pow(gap,3)/pow(br,3));
+//		lc=pow(gap,2)/lr;
+
 		/*Compute frictional heat flux*/
 		friction->GetAlpha2(&alpha2,gauss);
@@ -324,12 +371,35 @@
 		_assert_(meltrate>0.);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
-		 (
-		  meltrate*(1/rho_water-1/rho_ice)
-		  +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
-		  -beta*sqrt(vx*vx+vy*vy)
-		  +ieb
-		  +storage*head_old/dt
-		  )*basis[i];     	
+
+//		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
+//		 (
+//		  meltrate*(1/rho_water-1/rho_ice)
+//		  +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
+//		  -beta*sqrt(vx*vx+vy*vy)
+//		  +ieb
+//		  +storage*head_old/dt
+//		  )*basis[i];     
+
+		 /* With weighted creep term*/
+                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
+                  (
+                   meltrate*(1/rho_water-1/rho_ice)
+                   +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
+                   +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
+                   -beta*sqrt(vx*vx+vy*vy)
+                   +ieb
+                   +storage*head_old/dt
+                   )*basis[i];
+
+	
+//                 /*Test with linear creep rate*/
+ //                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
+  //                (
+   //                meltrate*(1/rho_water-1/rho_ice)
+    //               +(pressure_ice - pressure_water)*gap/pow(10,13)
+     //              -beta*sqrt(vx*vx+vy*vy)
+      //             +ieb
+       //            +storage*head_old/dt
+        //           )*basis[i];
 	}
 	/*Clean up and return*/
@@ -382,13 +452,13 @@
 		values[i]=solution[doflist[i]];
 
-		/*make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
-		if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
-			values[i] = rho_ice*thickness[i]/rho_water+bed[i];
-		}
+		/*overburden cap: make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
+//		if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
+//			values[i] = rho_ice*thickness[i]/rho_water+bed[i];
+//		}
 
 		/*Make sure that negative pressure is not allowed*/
-  //    if(values[i]<bed[i]){
-	//		values[i] = bed[i];
-	//	}
+//      if(values[i]<bed[i]){
+//			values[i] = bed[i];
+//		}
 
 		/*Under-relaxation*/
@@ -472,5 +542,5 @@
 	IssmDouble  Jdet,meltrate,G,dh[2],B,A,n,dt;
 	IssmDouble  gap,bed,thickness,head;
-	IssmDouble  lr,br,vx,vy,beta;
+	IssmDouble  lr,br,vx,vy,beta,lc;
 	IssmDouble  alpha2,frictionheat;
 	IssmDouble* xyz_list = NULL;
@@ -532,4 +602,12 @@
 		 beta = 0.;
 
+		/* TEST Compute lc term*/
+//		if(gap<br)
+//		 lc = gap*(1-(br-gap)/br);
+//		else
+		 lc = gap;
+//               lc=gap*tanh(pow(gap,3)/pow(br,3)); 
+//		lc = pow(gap,2)/lr;
+
 		/*Compute frictional heat flux*/
 		friction->GetAlpha2(&alpha2,gauss);
@@ -555,4 +633,20 @@
 					+beta*sqrt(vx*vx+vy*vy)
 					));
+
+		/* TEST with gap height "diffusivity" for melting walls - not yet completed*/
+//                newgap += gauss->weight*Jdet*(gap+dt*(
+//                                         meltrate/rho_ice
+//                                         -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*lc
+//                                         +beta*sqrt(vx*vx+vy*vy)
+//					 +(meltrate*gap/rho_ice*(dgapxx+dgapyy))
+//                                         ));
+
+
+                 /*TEST with linear creep rate*/
+//                 newgap += gauss->weight*Jdet*(gap+dt*(
+ //                                        meltrate/rho_ice
+  //                                       -(pressure_ice-pressure_water)*gap/pow(10,13)
+   //                                      +beta*sqrt(vx*vx+vy*vy)
+    //                                     ));
 		totalweights +=gauss->weight*Jdet;
 
