Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 26675)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 26676)
@@ -545,6 +545,21 @@
 
 	/*Get effective pressure and basal velocity*/
-	IssmDouble Neff = EffectivePressure(gauss);
 	IssmDouble vmag = VelMag(gauss);
+
+	bool ispwStochastic;
+	IssmDouble Neff;
+	element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
+	if(ispwStochastic){
+		/*Retrieve stochastic water pressure and compute ice pressure*/
+		IssmDouble p_ice,p_water,Neff_limit;
+		element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum);
+		element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
+		p_ice = IcePressure(gauss);
+		Neff  = max(Neff_limit*p_ice, p_ice - p_water);
+	}	
+	else{
+		/*Compute effective pressure directly*/
+		Neff = EffectivePressure(gauss);
+	}
 
 	/*Check to prevent dividing by zero if vmag==0*/
@@ -845,14 +860,14 @@
 	element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
 
+	/*Compute ice pressure*/
+	p_ice = IcePressure(gauss);
+
 	/*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
 	switch(coupled_flag){
 		case 0:{
-			element->GetInputValue(&thickness, gauss,ThicknessEnum);
 			element->GetInputValue(&base, gauss,BaseEnum);
 			element->GetInputValue(&sealevel, gauss,SealevelEnum);
 			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
 			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
-			p_ice   = gravity*rho_ice*thickness;
 			p_water = rho_water*gravity*(sealevel-base);
 			Neff = p_ice - p_water;
@@ -860,8 +875,4 @@
 			break;
 		case 1:{
-			element->GetInputValue(&thickness, gauss,ThicknessEnum);
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
-			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
-			p_ice   = gravity*rho_ice*thickness;
 			p_water = 0.;
 			Neff = p_ice - p_water;
@@ -869,11 +880,8 @@
 			break;
 		case 2:{
-			element->GetInputValue(&thickness, gauss,ThicknessEnum);
 			element->GetInputValue(&base, gauss,BaseEnum);
 			element->GetInputValue(&sealevel, gauss,SealevelEnum);
 			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
 			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
-			p_ice   = gravity*rho_ice*thickness;
 			p_water = max(0.,rho_water*gravity*(sealevel-base));
 			Neff = p_ice - p_water;
@@ -882,16 +890,8 @@
 		case 3:{
 			element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
-			element->GetInputValue(&thickness, gauss,ThicknessEnum);
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
-			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
-			p_ice   = gravity*rho_ice*thickness;
 		}
 			break;
 		case 4:{
 			element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
-			element->GetInputValue(&thickness, gauss,ThicknessEnum);
-			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
-			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
-			p_ice   = gravity*rho_ice*thickness;
 		}
 			break;
@@ -905,4 +905,66 @@
 	/*Return effective pressure*/
 	return Neff;
+
+}/*}}}*/
+IssmDouble Friction::IcePressure(Gauss* gauss){/*{{{*/
+	/*Get ice pressure*/
+
+	IssmDouble  thickness,p_ice;
+	/*Recover Inputs and Parameters*/
+	element->GetInputValue(&thickness, gauss,ThicknessEnum);
+	IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
+	IssmDouble gravity = element->FindParam(ConstantsGEnum);
+
+	/*Compute*/
+	p_ice = gravity*rho_ice*thickness;
+
+	/*Return ice pressure*/
+	return p_ice;
+
+}/*}}}*/
+IssmDouble Friction::SubglacialWaterPressure(Gauss* gauss){/*{{{*/
+	/*Get water pressure as a function of  flag */
+
+	int         coupled_flag;
+	IssmDouble  base,sealevel,p_water;
+
+	/*Recover parameters: */
+	element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
+
+	switch(coupled_flag){
+		case 0:{
+			element->GetInputValue(&base, gauss,BaseEnum);
+			element->GetInputValue(&sealevel, gauss,SealevelEnum);
+			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
+			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
+			p_water = rho_water*gravity*(sealevel-base);
+		}
+			break;
+		case 1:{
+			p_water = 0.;
+		}
+			break;
+		case 2:{
+			element->GetInputValue(&base, gauss,BaseEnum);
+			element->GetInputValue(&sealevel, gauss,SealevelEnum);
+			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
+			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
+			p_water = max(0.,rho_water*gravity*(sealevel-base));
+		}
+			break;
+		case 3:{
+			_error_("water pressure not computed for coupling==3 in friction law");
+		}
+			break;
+		case 4:{
+			_error_("water pressure not computed for coupling==4 in friction law");
+		}
+			break;
+		default:
+			_error_("not supported");
+	}
+
+	/*Return water pressure*/
+	return p_water;
 
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.h	(revision 26675)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h	(revision 26676)
@@ -52,4 +52,6 @@
 
 		IssmDouble EffectivePressure(Gauss* gauss);
+		IssmDouble IcePressure(Gauss* gauss);
+		IssmDouble SubglacialWaterPressure(Gauss* gauss);
 		IssmDouble VelMag(Gauss* gauss);
 		void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss);
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26675)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26676)
@@ -501,4 +501,6 @@
 	bool isstochasticforcing;
    parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
+	/*Stochastic Effective Pressure false by default*/
+	parameters->AddObject(new BoolParam(StochasticForcingIsWaterPressureEnum,false));
    if(isstochasticforcing){
       int num_fields,stochastic_dim;
Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26675)
+++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26676)
@@ -4,4 +4,5 @@
 
 #include "./StochasticForcingx.h"
+#include "../../classes/Loads/Friction.h"
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
@@ -141,4 +142,50 @@
 					}
 					break;
+				//VV working(29Nov2021)
+				case FrictionWaterPressureEnum:
+					/*Specify that WaterPressure is stochastic*/ 
+					femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
+					for(Object* &object:femmodel->elements->objects){
+                  Element* element = xDynamicCast<Element*>(object);
+                  int numvertices  = element->GetNumberOfVertices();
+                  IssmDouble p_water_deterministic[numvertices];
+                  IssmDouble p_water[numvertices];
+						element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
+						Gauss* gauss=element->NewGauss();
+						Friction* friction = new Friction(element);
+						for(int i=0;i<numvertices;i++){
+							gauss->GaussVertex(i);
+							p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
+							p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
+							p_water[i]               = max(0.0,p_water[i]);
+						}
+						element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
+						delete gauss;
+						delete friction;
+					}
+					break;
+
+				/*
+				case FrictionEffectivePressureEnum:
+					femmodel->parameters->SetParam(true,StochasticForcingIsEffectivePressureEnum);
+					for(Object* &object:femmodel->elements->objects){
+                  Element* element = xDynamicCast<Element*>(object);
+                  int numvertices  = element->GetNumberOfVertices();
+                  IssmDouble Neff[numvertices];
+						element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
+						Gauss* gauss=element->NewGauss();
+						Friction* friction = new Friction(element);
+						for(int i=0;i<numvertices;i++){
+							gauss->GaussVertex(i);
+							Neff[i] = friction->EffectivePressure(gauss);
+							Neff[i] = Neff[i]+noisefield[dimensionid];
+						}
+						element->AddInput(FrictionEffectivePressureEnum,Neff,P1DGEnum);
+						delete gauss;
+						delete friction;
+					}
+					break;
+				*/
+
 				default:
 					_error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26675)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26676)
@@ -401,5 +401,7 @@
 syn keyword cConstant StochasticForcingDimensionsEnum
 syn keyword cConstant StochasticForcingFieldsEnum
+syn keyword cConstant StochasticForcingIsEffectivePressureEnum
 syn keyword cConstant StochasticForcingIsStochasticForcingEnum
+syn keyword cConstant StochasticForcingIsWaterPressureEnum
 syn keyword cConstant StochasticForcingNumFieldsEnum
 syn keyword cConstant StochasticForcingRandomflagEnum
@@ -603,4 +605,5 @@
 syn keyword cConstant BaselineBasalforcingsFloatingiceMeltingRateEnum
 syn keyword cConstant BaselineCalvingCalvingrateEnum
+syn keyword cConstant BaselineFrictionEffectivePressureEnum
 syn keyword cConstant BedEnum
 syn keyword cConstant BedGRDEnum
@@ -1073,4 +1076,6 @@
 syn keyword cConstant WaterfractionEnum
 syn keyword cConstant WaterheightEnum
+syn keyword cConstant FrictionWaterPressureEnum
+syn keyword cConstant FrictionWaterPressureNoiseEnum
 syn keyword cConstant WeightsLevelsetObservationEnum
 syn keyword cConstant WeightsSurfaceObservationEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26675)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26676)
@@ -395,5 +395,7 @@
 	StochasticForcingDimensionsEnum,
 	StochasticForcingFieldsEnum,
+	StochasticForcingIsEffectivePressureEnum,
 	StochasticForcingIsStochasticForcingEnum,
+	StochasticForcingIsWaterPressureEnum,
 	StochasticForcingNumFieldsEnum,
 	StochasticForcingRandomflagEnum,
@@ -599,4 +601,5 @@
 	BaselineBasalforcingsFloatingiceMeltingRateEnum,
 	BaselineCalvingCalvingrateEnum,
+	BaselineFrictionEffectivePressureEnum,
 	BedEnum,
 	BedGRDEnum,
@@ -1070,4 +1073,6 @@
 	WaterfractionEnum,
 	WaterheightEnum,
+	FrictionWaterPressureEnum,
+	FrictionWaterPressureNoiseEnum,
 	WeightsLevelsetObservationEnum,
 	WeightsSurfaceObservationEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26675)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26676)
@@ -403,5 +403,7 @@
 		case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
 		case StochasticForcingFieldsEnum : return "StochasticForcingFields";
+		case StochasticForcingIsEffectivePressureEnum : return "StochasticForcingIsEffectivePressure";
 		case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
+		case StochasticForcingIsWaterPressureEnum : return "StochasticForcingIsWaterPressure";
 		case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
 		case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
@@ -605,4 +607,5 @@
 		case BaselineBasalforcingsFloatingiceMeltingRateEnum : return "BaselineBasalforcingsFloatingiceMeltingRate";
 		case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
+		case BaselineFrictionEffectivePressureEnum : return "BaselineFrictionEffectivePressure";
 		case BedEnum : return "Bed";
 		case BedGRDEnum : return "BedGRD";
@@ -1075,4 +1078,6 @@
 		case WaterfractionEnum : return "Waterfraction";
 		case WaterheightEnum : return "Waterheight";
+		case FrictionWaterPressureEnum : return "FrictionWaterPressure";
+		case FrictionWaterPressureNoiseEnum : return "FrictionWaterPressureNoise";
 		case WeightsLevelsetObservationEnum : return "WeightsLevelsetObservation";
 		case WeightsSurfaceObservationEnum : return "WeightsSurfaceObservation";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26675)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26676)
@@ -394,5 +394,7 @@
 syn keyword juliaConstC StochasticForcingDimensionsEnum
 syn keyword juliaConstC StochasticForcingFieldsEnum
+syn keyword juliaConstC StochasticForcingIsEffectivePressureEnum
 syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
+syn keyword juliaConstC StochasticForcingIsWaterPressureEnum
 syn keyword juliaConstC StochasticForcingNumFieldsEnum
 syn keyword juliaConstC StochasticForcingRandomflagEnum
@@ -596,4 +598,5 @@
 syn keyword juliaConstC BaselineBasalforcingsFloatingiceMeltingRateEnum
 syn keyword juliaConstC BaselineCalvingCalvingrateEnum
+syn keyword juliaConstC BaselineFrictionEffectivePressureEnum
 syn keyword juliaConstC BedEnum
 syn keyword juliaConstC BedGRDEnum
@@ -1066,4 +1069,6 @@
 syn keyword juliaConstC WaterfractionEnum
 syn keyword juliaConstC WaterheightEnum
+syn keyword juliaConstC FrictionWaterPressureEnum
+syn keyword juliaConstC FrictionWaterPressureNoiseEnum
 syn keyword juliaConstC WeightsLevelsetObservationEnum
 syn keyword juliaConstC WeightsSurfaceObservationEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26675)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26676)
@@ -412,5 +412,7 @@
 	      else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
 	      else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
+	      else if (strcmp(name,"StochasticForcingIsEffectivePressure")==0) return StochasticForcingIsEffectivePressureEnum;
 	      else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
+	      else if (strcmp(name,"StochasticForcingIsWaterPressure")==0) return StochasticForcingIsWaterPressureEnum;
 	      else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
 	      else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
@@ -504,10 +506,10 @@
 	      else if (strcmp(name,"Step")==0) return StepEnum;
 	      else if (strcmp(name,"Steps")==0) return StepsEnum;
-	      else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
-	      else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
+	      if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
+	      else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
+	      else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
 	      else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
 	      else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
@@ -617,4 +619,5 @@
 	      else if (strcmp(name,"BaselineBasalforcingsFloatingiceMeltingRate")==0) return BaselineBasalforcingsFloatingiceMeltingRateEnum;
 	      else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
+	      else if (strcmp(name,"BaselineFrictionEffectivePressure")==0) return BaselineFrictionEffectivePressureEnum;
 	      else if (strcmp(name,"Bed")==0) return BedEnum;
 	      else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
@@ -626,11 +629,11 @@
 	      else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
 	      else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum;
-	      else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
-	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
-	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
+	      if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
+	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
+	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
+	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
 	      else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
 	      else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
@@ -749,11 +752,11 @@
 	      else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
 	      else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
-	      else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
-	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
-	      else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
+	      if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
+	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
+	      else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
+	      else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
 	      else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
 	      else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
@@ -872,11 +875,11 @@
 	      else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
 	      else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
-	      else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
-	      else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
-	      else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
+	      if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
+	      else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
+	      else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
+	      else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
 	      else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
 	      else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
@@ -995,11 +998,11 @@
 	      else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
 	      else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
-	      else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
-	      else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
-	      else if (strcmp(name,"SmbT")==0) return SmbTEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
+	      if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
+	      else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
+	      else if (strcmp(name,"SmbT")==0) return SmbTEnum;
+	      else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
 	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
 	      else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
@@ -1099,4 +1102,6 @@
 	      else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
 	      else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
+	      else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
+	      else if (strcmp(name,"FrictionWaterPressureNoise")==0) return FrictionWaterPressureNoiseEnum;
 	      else if (strcmp(name,"WeightsLevelsetObservation")==0) return WeightsLevelsetObservationEnum;
 	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
@@ -1116,13 +1121,13 @@
 	      else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
 	      else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
-	      else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
 	      else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
 	      else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
 	      else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
 	      else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
+	      else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
 	      else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
 	      else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
@@ -1239,13 +1244,13 @@
 	      else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
 	      else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
-	      else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
 	      else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
 	      else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
 	      else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
 	      else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
+	      else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
 	      else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
 	      else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
@@ -1362,13 +1367,13 @@
 	      else if (strcmp(name,"Indexed")==0) return IndexedEnum;
 	      else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
-	      else if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
 	      else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
 	      else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
 	      else if (strcmp(name,"IntParam")==0) return IntParamEnum;
 	      else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"Inputs")==0) return InputsEnum;
+	      else if (strcmp(name,"Inputs")==0) return InputsEnum;
 	      else if (strcmp(name,"Internal")==0) return InternalEnum;
 	      else if (strcmp(name,"Intersect")==0) return IntersectEnum;
@@ -1485,13 +1490,13 @@
 	      else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
 	      else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
-	      else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
 	      else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
 	      else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum;
 	      else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
 	      else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
+	      else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
 	      else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
 	      else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26675)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26676)
@@ -72,4 +72,16 @@
 					end
 				end
+				if(contains(field,'WaterPressure'))
+               mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
+               if~(isequal(class(md.friction),char(mdname)))
+                  error('stochasticforcing field %s is only implemented for default friction', char(field));
+               end
+               if(md.friction.coupling~=0 && md.friction.coupling~=1 && md.friction.coupling~=2)
+                  error('stochasticforcing field %s is only implemented for cases md.friction.coupling 0 or 1 or 2', char(field));
+               end
+               if(any(md.friction.q==0))
+                  error('stochasticforcing field %s requires non-zero q exponent',char(field));
+               end
+            end
 				%Checking for specific dimensions
 				if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
@@ -185,4 +197,5 @@
 		'DefaultCalving',...
 		'FloatingMeltRate',...
+		'FrictionWaterPressure',...
 		'FrontalForcingsRignotAutoregression',...
 		'SMBautoregression'
@@ -191,4 +204,5 @@
 		'calving',...
 		'basalforcings',...
+		'friction',...
 		'frontalforcingsrignotautoregression',...
 		'SMBautoregression'
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 26675)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 26676)
@@ -70,17 +70,26 @@
                 mdname = structstoch[field]
                 if (type(md.smb).__name__ != mdname):
-                    raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(mdname))
+                    raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(field))
             if 'FrontalForcings' in field:
                 mdname = structstoch[field]
                 if (type(md.frontalforcings).__name__ != mdname):
-                    raise TypeError('md.frontalforcings does not agree with stochasticforcing field {}'.format(mdname))
+                    raise TypeError('md.frontalforcings does not agree with stochasticforcing field {}'.format(field))
             if 'Calving' in field:
                 mdname = structstoch[field]
                 if (type(md.calving).__name__ != mdname):
-                    raise TypeError('md.calving does not agree with stochasticforcing field {}'.format(mdname))
+                    raise TypeError('md.calving does not agree with stochasticforcing field {}'.format(field))
             if 'BasalforcingsFloatingice' in field:
                 mdname = structstoch[field]
                 if (type(md.basalforcings).__name__ != mdname):
-                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname))
+                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(field))
+            if 'WaterPressure' in field:
+                mdname = structstoch[field]
+                if (type(md.friction).__name__ != mdname):
+                    raise TypeError('stochasticforcing field {} is only implemented for default friction'.format(field))
+                if (md.friction.coupling!=0 and md.friction.coupling!=1 and md.friction.coupling!=2):
+                    raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
+                if (np.any(md.friction.q==0):
+                        raise TypeError('stochasticforcing field {} requires non-zero q exponent'.format(field))
+
             # Checking for specific dimensions
             if not (field == 'SMBautoregression' or field == 'FrontalForcingsRignotAutoregression'):
@@ -173,4 +182,5 @@
         structure = {'DefaultCalving': 'calving',
                      'FloatingMeltRate': 'basalforcings',
+                     'FrictionWaterPressure': 'friction',
                      'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',
                      'SMBautoregression': 'SMBautoregression'}
Index: /issm/trunk-jpl/test/NightlyRun/test621.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test621.m	(revision 26676)
+++ /issm/trunk-jpl/test/NightlyRun/test621.m	(revision 26676)
@@ -0,0 +1,77 @@
+%Test Name: 79NorthStochFrictionWaterPressure
+md=triangle(model(),'../Exp/79North.exp',6000.);
+md=setmask(md,'../Exp/79NorthShelf.exp','');
+md=parameterize(md,'../Par/79North.par');
+md=setflowequation(md,'SSA','all');
+
+%Default friction
+md.friction             = friction();
+md.friction.coefficient = 30*ones(md.mesh.numberofvertices,1);
+md.friction.p           = 1*ones(md.mesh.numberofelements,1);
+md.friction.q           = 1*ones(md.mesh.numberofelements,1);
+
+% Basin separation default
+idb_df = zeros(md.mesh.numberofelements,1);
+iid1   = find(md.mesh.y<=-1.08e6);
+for ii=1:md.mesh.numberofelements
+    for vertex=1:3
+        if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1
+            idb_df(ii) = 1;
+        end
+    end
+    if idb_df(ii)==0 %no vertex was found in basin 1
+        idb_df(ii) = 2;
+    end
+end
+%Covariance matrix
+covPw      = 0.5e10*eye(2);
+covPw(1,1) = 1.5*covPw(1,1);
+
+%Stochastic forcing
+md.stochasticforcing.isstochasticforcing = 1;
+md.stochasticforcing.fields              = [{'FrictionWaterPressure'}];
+md.stochasticforcing.defaultdimension    = 2;
+md.stochasticforcing.default_id          = idb_df;
+md.stochasticforcing.covariance          = covPw; %global covariance
+md.stochasticforcing.randomflag          = 0; %determines true/false randomness
+
+md.transient.issmb              = 0;
+md.transient.ismasstransport    = 1;
+md.transient.isstressbalance    = 1;
+md.transient.isthermal          = 0;
+md.transient.isgroundingline    = 0;
+
+md.transient.requested_outputs = {'default','FrictionWaterPressure'};
+md.timestepping.start_time = 0;
+md.timestepping.time_step  = 1;
+md.timestepping.final_time = 5;
+md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Transient');
+
+
+%Fields and tolerances to track changes
+field_names      = {'Vx1','Vy1','Vel1','Thickness1','FrictionWaterPressure1',...
+                    'Vx2','Vy2','Vel2','Thickness2','FrictionWaterPressure2',...
+                    'Vx5','Vy5','Vel5','Thickness5','FrictionWaterPressure5'};
+field_tolerances={2e-10,2e-10,2e-10,2e-10,2e-10,...
+                  4e-10,4e-10,4e-10,4e-10,4e-10,...
+                  8e-10,8e-10,8e-10,8e-10,8e-10};
+              
+field_values={...
+	(md.results.TransientSolution(1).Vx),...
+    (md.results.TransientSolution(1).Vy),...
+    (md.results.TransientSolution(1).Vel),...
+    (md.results.TransientSolution(1).Thickness),...
+    (md.results.TransientSolution(1).FrictionWaterPressure),...
+    (md.results.TransientSolution(2).Vx),...
+    (md.results.TransientSolution(2).Vy),...
+    (md.results.TransientSolution(2).Vel),...
+    (md.results.TransientSolution(2).Thickness),...
+    (md.results.TransientSolution(2).FrictionWaterPressure),...
+    (md.results.TransientSolution(5).Vx),...
+    (md.results.TransientSolution(5).Vy),...
+    (md.results.TransientSolution(5).Vel),...
+    (md.results.TransientSolution(5).Thickness),...
+    (md.results.TransientSolution(5).FrictionWaterPressure),...
+    };
+
Index: /issm/trunk-jpl/test/NightlyRun/test621.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test621.py	(revision 26676)
+++ /issm/trunk-jpl/test/NightlyRun/test621.py	(revision 26676)
@@ -0,0 +1,91 @@
+#Test Name: 79NorthStochFrictionWaterPressure
+import numpy as np
+
+from socket import gethostname
+from model import *
+from parameterize import *
+from setflowequation import *
+from setmask import *
+from solve import *
+from triangle import *
+
+
+md = triangle(model(), '../Exp/79North.exp', 6000)
+md = setmask(md, '../Exp/79NorthShelf.exp')
+md = parameterize(md, '../Par/79North.py')
+md = setflowequation(md, 'SSA', 'all')
+
+#Default friction
+md.friction         = friction()
+md.friction.coefficient = 30*np.ones(md.mesh.numberofvertices)
+md.friction.p           = 1*np.ones((md.mesh.numberofelements))
+md.friction.q           = 1*np.ones((md.mesh.numberofelements))
+
+# Basin separation default
+idb_df = np.zeros((md.mesh.numberofelements,))
+iid1 = np.where(md.mesh.y<=-1.08e6)[0]
+for ii in range(md.mesh.numberofelements):
+    for vertex in range(3):
+        if md.mesh.elements[ii][vertex] - 1 in iid1:  # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing
+            idb_df[ii] = 1
+    if idb_df[ii] == 0:  # no vertex was found in basin 1
+        for vertex in range(3):
+            idb_df[ii] = 2
+#Covariance matrix
+covPw = np.array([[0.75e10,0.0],[0.0,0.5e10]])
+
+# Stochastic forcing
+md.stochasticforcing.isstochasticforcing = 1
+md.stochasticforcing.fields              = ['FrontalForcingsRignotAutoregression','DefaultCalving','FloatingMeltRate']
+md.stochasticforcing.defaultdimension    = 2
+md.stochasticforcing.default_id          = idb_df-1 #NOTE: offset because of 1-based vertex indexing
+md.stochasticforcing.covariance          = covPw # global covariance
+md.stochasticforcing.randomflag          = 0 # determines true/false randomness
+
+md.transient.issmb              = 0;
+md.transient.ismasstransport    = 1;
+md.transient.isstressbalance    = 1;
+md.transient.isthermal          = 0;
+md.transient.isgroundingline    = 0;
+
+md.transient.requested_outputs = ['default', 'FrictionWaterPressure']
+md.timestepping.start_time = 0
+md.timestepping.time_step  = 1
+md.timestepping.final_time = 5
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md, 'Transient')
+
+# Fields and tolerances to track changes
+field_names = [
+    'Vx1','Vy1','Vel1','Thickness1','FrictionWaterPressure1',
+    'Vx2','Vy2','Vel2','Thickness2','FrictionWaterPressure2',
+    'Vx10','Vy10','Vel10','Thickness10','FrictionWaterPressure10'
+    ]
+
+field_tolerances = [
+    2e-10,2e-10,2e-10,2e-10,2e-10,
+    4e-10,4e-10,4e-10,4e-10,4e-10,
+    8e-10,8e-10,8e-10,8e-10,8e-10
+    ]
+
+field_values = [
+    md.results.TransientSolution[0].Vx,
+    md.results.TransientSolution[0].Vy,
+    md.results.TransientSolution[0].Vel,
+    md.results.TransientSolution[0].Thickness,
+    md.results.TransientSolution[0].FrictionWaterPressure,
+    md.results.TransientSolution[1].Vx,
+    md.results.TransientSolution[1].Vy,
+    md.results.TransientSolution[1].Vel,
+    md.results.TransientSolution[1].Thickness,
+    md.results.TransientSolution[1].FrictionWaterPressure,
+    md.results.TransientSolution[4].Vx,
+    md.results.TransientSolution[4].Vy,
+    md.results.TransientSolution[4].Vel,
+    md.results.TransientSolution[4].Thickness,
+    md.results.TransientSolution[4].FrictionWaterPressure
+    ]
+    
+
+
+