Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 27160)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 27161)
@@ -749,5 +749,5 @@
 	 *               C |u_b|^(m-1)
 	 * alpha2= __________________________
-	 *          (1+(C/(Cmax N))^1/m |u_b| )^m
+	 *          (1+(C/(Cmax Neff))^1/m |u_b| )^m
 	 *
 	 * */
@@ -764,14 +764,30 @@
 	C = coeff*coeff;
 
-	/*Get effective pressure and velocity magnitude*/
-	IssmDouble N  = EffectivePressure(gauss);
+	/*Get effective pressure*/
+	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,FrictionSchoofWaterPressureEnum);
+		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);
+	}
+
+	/*Get velocity magnitude*/
 	IssmDouble ub = VelMag(gauss);
 
 	/*Compute alpha^2*/
-	if((ub<1e-10) ||(N==0.0)){
+	if((ub<1e-10) ||(Neff==0.0)){
 		alpha2 = 0.;
 	}
 	else{
-		alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
+		alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*Neff),1./m)*ub,m);
 	}
 
Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27160)
+++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27161)
@@ -214,5 +214,25 @@
 					}
 					break;
-				default:
+				case FrictionSchoofWaterPressureEnum:
+					/*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)
+						}
+						element->AddInput(FrictionSchoofWaterPressureEnum,p_water,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 27160)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27161)
@@ -754,5 +754,7 @@
 syn keyword cConstant FrictionSedimentCompressibilityCoefficientEnum
 syn keyword cConstant FrictionTillFrictionAngleEnum
+syn keyword cConstant FrictionSchoofWaterPressureEnum
 syn keyword cConstant FrictionWaterLayerEnum
+syn keyword cConstant FrictionWaterPressureEnum
 syn keyword cConstant FrictionfEnum
 syn keyword cConstant FrontalForcingsBasinIdEnum
@@ -1123,6 +1125,4 @@
 syn keyword cConstant WaterfractionEnum
 syn keyword cConstant WaterheightEnum
-syn keyword cConstant FrictionWaterPressureEnum
-syn keyword cConstant FrictionWaterPressureNoiseEnum
 syn keyword cConstant WeightsLevelsetObservationEnum
 syn keyword cConstant WeightsSurfaceObservationEnum
@@ -1639,4 +1639,5 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1645,6 +1646,6 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1659,6 +1660,6 @@
 syn keyword cType ElementInput
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1667,9 +1668,10 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1688,4 +1690,5 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1698,4 +1701,5 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1710,6 +1714,6 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1725,11 +1729,11 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType SealevelGeometry
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1750,8 +1754,4 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27160)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27161)
@@ -750,5 +750,7 @@
 	FrictionSedimentCompressibilityCoefficientEnum,
 	FrictionTillFrictionAngleEnum,
+	FrictionSchoofWaterPressureEnum,
 	FrictionWaterLayerEnum,
+	FrictionWaterPressureEnum,
 	FrictionfEnum,
 	FrontalForcingsBasinIdEnum,
@@ -1120,6 +1122,4 @@
 	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 27160)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27161)
@@ -756,5 +756,7 @@
 		case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient";
 		case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle";
+		case FrictionSchoofWaterPressureEnum : return "FrictionSchoofWaterPressure";
 		case FrictionWaterLayerEnum : return "FrictionWaterLayer";
+		case FrictionWaterPressureEnum : return "FrictionWaterPressure";
 		case FrictionfEnum : return "Frictionf";
 		case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
@@ -1125,6 +1127,4 @@
 		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 27160)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27161)
@@ -747,5 +747,7 @@
 syn keyword juliaConstC FrictionSedimentCompressibilityCoefficientEnum
 syn keyword juliaConstC FrictionTillFrictionAngleEnum
+syn keyword juliaConstC FrictionSchoofWaterPressureEnum
 syn keyword juliaConstC FrictionWaterLayerEnum
+syn keyword juliaConstC FrictionWaterPressureEnum
 syn keyword juliaConstC FrictionfEnum
 syn keyword juliaConstC FrontalForcingsBasinIdEnum
@@ -1116,6 +1118,4 @@
 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 27160)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27161)
@@ -774,5 +774,7 @@
 	      else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
 	      else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
+	      else if (strcmp(name,"FrictionSchoofWaterPressure")==0) return FrictionSchoofWaterPressureEnum;
 	      else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
+	      else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
 	      else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
 	      else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
@@ -873,10 +875,10 @@
 	      else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
 	      else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
-	      else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
-	      else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
+	      if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
+	      else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
+	      else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
 	      else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
 	      else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
@@ -996,10 +998,10 @@
 	      else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
 	      else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
-	      else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
-	      else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
+	      if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
+	      else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
+	      else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
 	      else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
 	      else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
@@ -1119,10 +1121,10 @@
 	      else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
 	      else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
-	      else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
-	      else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
+	      if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
+	      else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
+	      else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
 	      else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
 	      else if (strcmp(name,"Vel")==0) return VelEnum;
@@ -1152,6 +1154,4 @@
 	      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;
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27160)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27161)
@@ -95,14 +95,18 @@
 				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));
+					mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
+					if~(isequal(class(md.friction),char(mdname)))
+                  error('md.friction does not agree with stochasticforcing field %s', 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
+					if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof'))
+                  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
+					end
+					if(strcmp(class(md.friction),'friction'))
+                  if(any(md.friction.q==0))
+                     error('stochasticforcing field %s requires non-zero q exponent',char(field));
+                  end
+					end
             end
 				%Checking for specific dimensions
@@ -267,4 +271,5 @@
 		'FloatingMeltRate',...
 		'FrictionWaterPressure',...
+		'FrictionSchoofWaterPressure',...
 		'FrontalForcingsRignotAutoregression',...
 		'SMBautoregression',...
@@ -277,4 +282,5 @@
 		'basalforcings',...
 		'friction',...
+		'frictionschoof',...
 		'frontalforcingsrignotautoregression',...
 		'SMBautoregression',...
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27160)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27161)
@@ -106,9 +106,11 @@
                 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 not in[0, 1, 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))
+                    raise TypeError('md.friction does not agree with stochasticforcing field {}'.format(field))
+                if (type(md.friction).__name__=='friction' or type(md.friction).__name__=='frictionschoof'):
+                    if md.friction.coupling not in[0, 1, 2]:
+                        raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
+                if (type(md.friction).__name__=='friction'):
+                    if (np.any(md.friction.q == 0)):
+                        raise TypeError('stochasticforcing field {} requires non-zero q exponent'.format(field))
 
             # Checking for specific dimensions
@@ -232,4 +234,5 @@
                      'FloatingMeltRate': 'basalforcings',
                      'FrictionWaterPressure': 'friction',
+                     'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',
                      'SMBautoregression': 'SMBautoregression',