Index: issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 27937)
+++ issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 28050)
@@ -151,4 +151,6 @@
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.hydrology.rheology_B_base",HydrologyRheologyBBaseEnum);
+
 	if(iomodel->domaintype==Domain2DhorizontalEnum){
 		iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
@@ -249,12 +251,12 @@
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
 	IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
-	Input* hr_input      = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
+	Input* hr_input  = element->GetInput(HydrologyBumpHeightEnum);       _assert_(hr_input);
 	Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
-	Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
-	Input* h_input   = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input);
-	Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
-	Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
-	Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
-	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+	Input* phi_input = element->GetInput(HydraulicPotentialEnum);        _assert_(phi_input);
+	Input* h_input   = element->GetInput(HydrologySheetThicknessEnum);   _assert_(h_input);
+	Input* H_input   = element->GetInput(ThicknessEnum);                 _assert_(H_input);
+	Input* b_input   = element->GetInput(BedEnum);                       _assert_(b_input);
+	Input* B_input   = element->GetInput(HydrologyRheologyBBaseEnum);    _assert_(B_input);
+	Input* n_input   = element->GetInput(MaterialsRheologyNEnum);        _assert_(n_input);	
 
 	/* Start  looping on the number of gaussian points: */
@@ -275,7 +277,4 @@
 		b_input->GetInputValue(&b,gauss);
 		H_input->GetInputValue(&H,gauss);
-
-		/*Hard code B*/
-		B = Cuffey(273.15-2);
 		
 		/*Get norm of gradient of hydraulic potential and make sure it is >0*/
@@ -364,15 +363,15 @@
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
 	IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
-	Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
-	Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
-	Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
-	Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
-	Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
+	Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);                _assert_(hr_input);
+	Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);            _assert_(h_input);
+	Input* H_input      = element->GetInput(ThicknessEnum);                          _assert_(H_input);
+	Input* b_input      = element->GetInput(BedEnum);                                _assert_(b_input);
+	Input* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);        _assert_(G_input);
 	Input* melt_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);_assert_(melt_input);
 	Input* RO_input     = NULL;
-	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);
-	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
+	Input* B_input      = element->GetInput(HydrologyRheologyBBaseEnum);             _assert_(B_input);
+	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);                 _assert_(n_input);
+	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);              _assert_(phiold_input);
+	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);                 _assert_(phi_input);
 
 	/*Build friction element, needed later: */
@@ -397,7 +396,4 @@
 		melt_input->GetInputValue(&melt,gauss);
 
-		/*Hard code B*/
-		B = Cuffey(273.15-2);
-
 		/*Get basal velocity*/
 		friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss);
@@ -457,9 +453,9 @@
 	element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
 
-	/*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/
+	/*Compute Hydrology Vx and Vy for time stepping purposes, and Sheet Discharge as an optional output (These inputs do not affect GlaDS)*/
 
 	/*Intermediaries*/
    IssmDouble  dphi[3],h,k,phi;
-   	IssmDouble  h_r;
+	IssmDouble  h_r;
 	IssmDouble  oceanLS,iceLS;
 	IssmDouble* xyz_list = NULL;
@@ -468,7 +464,8 @@
 	int numvertices = element->GetNumberOfVertices();
 
-	/*Initialize new sheet thickness*/
+	/*Initialize water sheet velocity and discharge*/
 	IssmDouble* vx = xNew<IssmDouble>(numvertices);
 	IssmDouble* vy = xNew<IssmDouble>(numvertices);
+	IssmDouble* d = xNew<IssmDouble>(numvertices);
 
 	/*Set to 0 if inactive element*/
@@ -476,8 +473,11 @@
 		for(int iv=0;iv<numvertices;iv++) vx[iv] = 0.;
 		for(int iv=0;iv<numvertices;iv++) vy[iv] = 0.;
+		for(int iv=0;iv<numvertices;iv++) d[iv] = 0.;
 		element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
 		element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+		element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
 		xDelete<IssmDouble>(vx);
 		xDelete<IssmDouble>(vy);
+		xDelete<IssmDouble>(d);
 		return;
 	}
@@ -499,5 +499,5 @@
 	Input *iceLS_input   = element->GetInput(MaskIceLevelsetEnum);            _assert_(iceLS_input);
 
-	/* Start  looping on the number of gaussian points: */
+	/* Start looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss();
 	for(int iv=0;iv<numvertices;iv++){
@@ -508,13 +508,14 @@
       phi_input->GetInputValue(&phi,gauss);
       h_input->GetInputValue(&h,gauss);
-      	hr_input->GetInputValue(&h_r,gauss); 
+      hr_input->GetInputValue(&h_r,gauss); 
       k_input->GetInputValue(&k,gauss);
 		oceanLS_input->GetInputValue(&oceanLS,gauss);
 		iceLS_input->GetInputValue(&iceLS,gauss);
 
-		/*Set sheet thickness to zero if floating or no ice*/
+		/*Set to zero if floating or no ice*/
 		if(oceanLS<0. || iceLS>0.){
 			vx[iv] = 0.;
          vy[iv] = 0.;
+			d[iv] = 0.;
 		}
 		else{
@@ -526,17 +527,18 @@
          /*If omega is zero, use standard model, otherwise transition model*/
          IssmDouble nu = mu_water/rho_water;
-	IssmDouble coeff;
-	if(istransition==1 && omega>=AEPS){
-		IssmDouble hratio = fabs(h/h_r);
-		IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
-		coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi/max(AEPS,h);  // divide by h to get speed instead of discharge
-	}
-	else {
-		coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h);  // divide by h to get speed instead of discharge
-	}
-
-
-			vx[iv] = -coeff*dphi[0];
-			vy[iv] = -coeff*dphi[1];
+			IssmDouble coeff;
+			if(istransition==1 && omega>=AEPS){
+				IssmDouble hratio = fabs(h/h_r);
+				IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
+				coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi;  // coeff gives discharge; divide by h to get speed instead of discharge
+			}
+			else {
+			coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);  // coeff gives discharge; divide by h to get speed instead of discharge
+			}
+
+			vx[iv] = -coeff/max(AEPS,h)*dphi[0];
+			vy[iv] = -coeff/max(AEPS,h)*dphi[1];
+				
+			d[iv] = coeff*normgradphi;
 		}
 	}
@@ -544,4 +546,5 @@
 	element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
 	element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+	element->AddInput(HydrologySheetDischargeEnum,d,P1DGEnum);
 
 	/*Clean up and return*/
@@ -549,4 +552,5 @@
 	xDelete<IssmDouble>(vx);
 	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(d);
 	delete gauss;
 }/*}}}*/
@@ -646,16 +650,15 @@
 	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
-	Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
-	Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input);
-	Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input);
-	Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
-	Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
-	Input* hold_input  = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
-	Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
-	Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+	Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);         _assert_(hr_input);
+	Input* vx_input = element->GetInput(VxBaseEnum);                      _assert_(vx_input);
+	Input* vy_input = element->GetInput(VyBaseEnum);                      _assert_(vy_input);
+	Input* H_input = element->GetInput(ThicknessEnum);                    _assert_(H_input);
+	Input* b_input = element->GetInput(BedEnum);                          _assert_(b_input);
+	Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
+	Input* B_input = element->GetInput(HydrologyRheologyBBaseEnum);       _assert_(B_input);
+	Input* n_input = element->GetInput(MaterialsRheologyNEnum);           _assert_(n_input);
 	Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
-	Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
-	Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
-
+	Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);      _assert_(oceanLS_input);
+	Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum);          _assert_(iceLS_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -677,7 +680,4 @@
 		iceLS_input->GetInputValue(&iceLS,gauss);
 
-		/*Hard code B*/
-		B = Cuffey(273.15-2);
-
 		/*Set sheet thickness to zero if floating or no ice*/
 		if(oceanLS<0. || iceLS>0.){
