Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 27164)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 27165)
@@ -344,5 +344,21 @@
 
 	/*Get effective pressure*/
-	IssmDouble Neff = EffectivePressure(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,FrictionCoulombWaterPressureEnum);
+		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 vmag = VelMag(gauss);
 
Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27164)
+++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27165)
@@ -207,5 +207,5 @@
 							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]               = p_water_deterministic[i] + noisefield[dimensionid]; 
 						}
 						element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
@@ -228,5 +228,5 @@
 							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]               = p_water_deterministic[i] + noisefield[dimensionid]; 
 						}
 						element->AddInput(FrictionSchoofWaterPressureEnum,p_water,P1DGEnum);
@@ -234,5 +234,27 @@
 						delete friction;
 					}
-					break;				default:
+					break;
+				case FrictionCoulombWaterPressureEnum:
+					/*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]; 
+						}
+						element->AddInput(FrictionCoulombWaterPressureEnum,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 27164)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27165)
@@ -747,4 +747,5 @@
 syn keyword cConstant FrictionCoefficientEnum
 syn keyword cConstant FrictionCoefficientcoulombEnum
+syn keyword cConstant FrictionCoulombWaterPressureEnum
 syn keyword cConstant FrictionEffectivePressureEnum
 syn keyword cConstant FrictionMEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27164)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27165)
@@ -743,4 +743,5 @@
 	FrictionCoefficientEnum,
 	FrictionCoefficientcoulombEnum,
+	FrictionCoulombWaterPressureEnum,
 	FrictionEffectivePressureEnum,
 	FrictionMEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27164)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27165)
@@ -749,4 +749,5 @@
 		case FrictionCoefficientEnum : return "FrictionCoefficient";
 		case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
+		case FrictionCoulombWaterPressureEnum : return "FrictionCoulombWaterPressure";
 		case FrictionEffectivePressureEnum : return "FrictionEffectivePressure";
 		case FrictionMEnum : return "FrictionM";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27164)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27165)
@@ -740,4 +740,5 @@
 syn keyword juliaConstC FrictionCoefficientEnum
 syn keyword juliaConstC FrictionCoefficientcoulombEnum
+syn keyword juliaConstC FrictionCoulombWaterPressureEnum
 syn keyword juliaConstC FrictionEffectivePressureEnum
 syn keyword juliaConstC FrictionMEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27164)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27165)
@@ -767,4 +767,5 @@
 	      else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
 	      else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
+	      else if (strcmp(name,"FrictionCoulombWaterPressure")==0) return FrictionCoulombWaterPressureEnum;
 	      else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
 	      else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
 	      else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
-	      else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
+	      if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
+	      else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
 	      else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
 	      else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum;
 	      else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
-	      else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
+	      if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
+	      else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
 	      else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
 	      else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
@@ -1120,9 +1121,9 @@
 	      else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
 	      else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
-	      else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
+	      if (strcmp(name,"Thickness")==0) return ThicknessEnum;
+	      else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
 	      else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
 	      else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
@@ -1243,9 +1244,9 @@
 	      else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
 	      else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
-	      else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
+	      if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
+	      else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
 	      else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
 	      else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
@@ -1366,9 +1367,9 @@
 	      else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
 	      else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
-	      else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
          else stage=12;
    }
    if(stage==12){
-	      if (strcmp(name,"Free")==0) return FreeEnum;
+	      if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
+	      else if (strcmp(name,"Free")==0) return FreeEnum;
 	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
 	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
@@ -1489,9 +1490,9 @@
 	      else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
 	      else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
-	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
          else stage=13;
    }
    if(stage==13){
-	      if (strcmp(name,"MeshY")==0) return MeshYEnum;
+	      if (strcmp(name,"MeshX")==0) return MeshXEnum;
+	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
 	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
 	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
@@ -1612,9 +1613,9 @@
 	      else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
 	      else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
-	      else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
          else stage=14;
    }
    if(stage==14){
-	      if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
+	      if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
+	      else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
 	      else if (strcmp(name,"Tetra")==0) return TetraEnum;
 	      else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27164)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27165)
@@ -99,5 +99,5 @@
                   error('md.friction does not agree with stochasticforcing field %s', char(field));
                end
-					if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof'))
+					if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof') || strcmp(class(md.friction),'frictioncoulomb'))
                   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));
@@ -271,4 +271,5 @@
 		'FloatingMeltRate',...
 		'FrictionWaterPressure',...
+		'FrictionCoulombWaterPressure',...
 		'FrictionSchoofWaterPressure',...
 		'FrontalForcingsRignotAutoregression',...
@@ -282,4 +283,5 @@
 		'basalforcings',...
 		'friction',...
+		'frictioncoulomb',...
 		'frictionschoof',...
 		'frontalforcingsrignotautoregression',...
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27164)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27165)
@@ -107,5 +107,5 @@
                 if (type(md.friction).__name__ != mdname):
                     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 (type(md.friction).__name__=='friction' or type(md.friction).__name__=='frictionschoof' or type(md.friction).__name__=='frictioncoulomb'):
                     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))
@@ -234,4 +234,5 @@
                      'FloatingMeltRate': 'basalforcings',
                      'FrictionWaterPressure': 'friction',
+                     'FrictionCoulombWaterPressure': 'frictioncoulomb',
                      'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',