Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 26490)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 26491)
@@ -272,5 +272,6 @@
    IssmDouble  PMPheat,dpressure_water[2],dbed[2];	
 	IssmDouble* xyz_list = NULL;
-        IssmDouble dgap[2]; /***/
+        IssmDouble dgapxx; /***/
+	IssmDouble dgapyy; /***/
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -298,6 +299,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); /***/
+	Input* dgapxx_input	    = element->GetInput(HydrologyGapHeightXXEnum); /***/
+	Input* dgapyy_input	    = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
 	/*Get conductivity from inputs*/
@@ -329,4 +330,6 @@
 		br_input->GetInputValue(&br,gauss);
                 headold_input->GetInputValue(&head_old,gauss);
+		dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
+		dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
 		/*Get ice A parameter*/
@@ -382,13 +385,25 @@
 
 		 /* 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];
+//                 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];
+
+                 /* 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
+		    +meltrate*gap/rho_ice*(dgapxx+dgapyy)
+                    )*basis[i];
 
 	
@@ -548,4 +563,6 @@
 	IssmDouble q = 0.;
    IssmDouble channelization = 0.;
+	IssmDouble dgapxx; /***/
+	IssmDouble dgapyy; /***/
 
 	/*Retrieve all inputs and parameters*/
@@ -565,4 +582,6 @@
 	Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
 	Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
+        Input* dgapxx_input         = element->GetInput(HydrologyGapHeightXXEnum); /***/
+        Input* dgapyy_input         = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
 	/*Get conductivity from inputs*/
@@ -590,4 +609,6 @@
 		lr_input->GetInputValue(&lr,gauss);
 		br_input->GetInputValue(&br,gauss);
+		dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
+		dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
 		/*Get ice A parameter*/
@@ -628,17 +649,17 @@
 		_assert_(meltrate>0.);
 
-		newgap += gauss->weight*Jdet*(gap+dt*(
-					meltrate/rho_ice
-					-A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
-					+beta*sqrt(vx*vx+vy*vy)
-					));
+//		newgap += gauss->weight*Jdet*(gap+dt*(
+//					meltrate/rho_ice
+//					-A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
+//					+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))
-//                                         ));
+                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))
+                                         ));