Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 27224)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp	(revision 27225)
@@ -242,5 +242,4 @@
 		element->NodalFunctions(basis,gauss);
 
-		/***/
                 base_input->GetInputValue(&bed,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
@@ -288,6 +287,4 @@
    IssmDouble  PMPheat,dissipation,dpressure_water[2],dbed[2];	
 	IssmDouble* xyz_list = NULL;
-//        IssmDouble dgapxx; /***/
-//	IssmDouble dgapyy; /***/
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -315,6 +312,4 @@
 	Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
-//	Input* dgapxx_input	    = element->GetInput(HydrologyGapHeightXXEnum); /***/
-//	Input* dgapyy_input	    = element->GetInput(HydrologyGapHeightYYEnum); /***/
 
 	/*Get conductivity from inputs*/
@@ -346,6 +341,4 @@
 		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*/
@@ -365,8 +358,4 @@
 		frictionheat=alpha2*(vx*vx+vy*vy);
 
-		/* Take out frictional heat for large gap height */
-		if(gap>br)
-		 frictionheat=0.;
-
 		/*Get water and ice pressures*/
 		IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
@@ -379,10 +368,8 @@
 
 		/*Compute change in sensible heat due to changes in pressure melting point*/
-   	dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
+   		dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
 		dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
-		PMPheat=CT*CW*conductivity*rho_water*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
-		PMPheat=0; /*** TEST no PMPheat***/
-
-   	meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
+
+   	meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
 
                   for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
@@ -396,16 +383,4 @@
                     )*basis[i];
 
-                /* Test with experimental diffusivity 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
-        //             +0.e-5*(dgapxx+dgapyy)
-         //            )*basis[i];
-
 	
 	}
@@ -458,14 +433,4 @@
 		values[i]=solution[doflist[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];
-//		}
-
 		/*Under-relaxation*/
 	   values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
@@ -549,9 +514,7 @@
 	IssmDouble  alpha2,frictionheat;
 	IssmDouble* xyz_list = NULL;
-   IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
+  	IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
 	IssmDouble q = 0.;
-   IssmDouble channelization = 0.;
-//	IssmDouble dgapxx; /***/
-//	IssmDouble dgapyy; /***/
+   	IssmDouble channelization = 0.;
 
 	/*Retrieve all inputs and parameters*/
@@ -571,6 +534,4 @@
 	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*/
@@ -598,6 +559,4 @@
 		lr_input->GetInputValue(&lr,gauss);
 		br_input->GetInputValue(&br,gauss);
-//		dgapxx_input->GetInputValue(&dgapxx,gauss); /***/
-//		dgapyy_input->GetInputValue(&dgapyy,gauss); /***/
 
 		/*Get ice A parameter*/
@@ -617,8 +576,4 @@
 		frictionheat=alpha2*(vx*vx+vy*vy);
 
-		/* Take out frictional heat for large gap height */
-		if(gap>br)
-		 frictionheat=0.;
-
 		/*Get water and ice pressures*/
 		IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
@@ -629,9 +584,7 @@
 	   dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
 		dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
-		PMPheat=CT*CW*conductivity*rho_water*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
-PMPheat=0; /*** TEST no PMPheat***/
 		dissipation=rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]);
 
-		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
+		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
 
 		element->AddInput(DummyEnum,&meltrate,P0Enum);
@@ -647,12 +600,4 @@
 					));
 
-		/* TEST with experimental gap height "diffusivity" for melting walls */
-//                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)
-//					 +0.e-5*(dgapxx+dgapyy))
- //                                        );
-
 
 		totalweights +=gauss->weight*Jdet;
@@ -667,10 +612,8 @@
 	/*Divide by connectivity*/
 	newgap = newgap/totalweights;
-	IssmDouble mingap = 1e-6;
+	IssmDouble mingap = 1e-3;
 	if(newgap<mingap) newgap=mingap;
 
-	/*Limit gap height to grow to surface*/
-//	if(newgap>thickness)
-//	 newgap = thickness;
+	/*Limit gap height*/
 	if(newgap>1)
 	 newgap = 1;