Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 26789)
@@ -108,7 +108,8 @@
 			break;
 		case SpatialLinearFloatingMeltRateEnum:
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
 			break;
 		case BasalforcingsPicoEnum:
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 26789)
@@ -182,7 +182,9 @@
 			break;
 		case SpatialLinearFloatingMeltRateEnum:
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.perturbation_melting_rate",BasalforcingsPerturbationMeltingRateEnum,0.);
 			break;
 		case BasalforcingsPicoEnum:
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26789)
@@ -836,7 +836,8 @@
 				break;
 			case SpatialLinearFloatingMeltRateEnum:
-				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
-				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
-				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
+				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
+				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
+				iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
 				break;
 			case BasalforcingsPicoEnum:
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 26789)
@@ -139,7 +139,8 @@
 			break;
 		case SpatialLinearFloatingMeltRateEnum:
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
-			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
+			iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
 			break;
 		case BasalforcingsPicoEnum:
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26789)
@@ -2298,19 +2298,32 @@
 	const int NUM_VERTICES = this->GetNumberOfVertices();
 
+	IssmDouble alpha;
 	IssmDouble deepwatermelt[MAXVERTICES];
-	IssmDouble deepwaterel[MAXVERTICES];;
+	IssmDouble deepwaterel[MAXVERTICES];
+	IssmDouble upperwatermelt[MAXVERTICES];
 	IssmDouble upperwaterel[MAXVERTICES];
+	IssmDouble perturbation[MAXVERTICES];
 	IssmDouble base[MAXVERTICES];
 	IssmDouble values[MAXVERTICES];
 
 	this->GetInputListOnVertices(&base[0],BaseEnum);
-	this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsDeepwaterMeltingRateEnum);
-	this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsDeepwaterElevationEnum);
-	this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsUpperwaterElevationEnum);
+	this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsSpatialDeepwaterMeltingRateEnum);
+	this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsSpatialDeepwaterElevationEnum);
+	this->GetInputListOnVertices(&upperwatermelt[0],BasalforcingsSpatialUpperwaterMeltingRateEnum);
+	this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsSpatialUpperwaterElevationEnum);
+   this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum);
 
 	for(int i=0;i<NUM_VERTICES;i++){
-		if(base[i]>upperwaterel[i])      values[i]=0;
-		else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
-		else values[i]=deepwatermelt[i]*(base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
+		if(base[i]>=upperwaterel[i]){
+			values[i]=upperwatermelt[i];
+		}
+		else if (base[i]<deepwaterel[i]){
+			values[i]=deepwatermelt[i];
+		}
+		else{
+			alpha = (base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
+			values[i]=deepwatermelt[i]*alpha+(1.-alpha)*upperwatermelt[i];
+		}
+		values[i] += perturbation[i];
 	}
 
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26788)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26789)
@@ -582,4 +582,8 @@
 syn keyword cConstant BasalforcingsGroundediceMeltingRateEnum
 syn keyword cConstant BasalforcingsPerturbationMeltingRateEnum
+syn keyword cConstant BasalforcingsSpatialDeepwaterElevationEnum
+syn keyword cConstant BasalforcingsSpatialDeepwaterMeltingRateEnum
+syn keyword cConstant BasalforcingsSpatialUpperwaterElevationEnum
+syn keyword cConstant BasalforcingsSpatialUpperwaterMeltingRateEnum
 syn keyword cConstant BasalforcingsIsmip6BasinIdEnum
 syn keyword cConstant BasalforcingsIsmip6TfEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26788)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26789)
@@ -578,4 +578,8 @@
 	BasalforcingsGroundediceMeltingRateEnum,
 	BasalforcingsPerturbationMeltingRateEnum,
+	BasalforcingsSpatialDeepwaterElevationEnum,
+	BasalforcingsSpatialDeepwaterMeltingRateEnum,
+	BasalforcingsSpatialUpperwaterElevationEnum,
+	BasalforcingsSpatialUpperwaterMeltingRateEnum,
 	BasalforcingsIsmip6BasinIdEnum,
 	BasalforcingsIsmip6TfEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26788)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26789)
@@ -575,4 +575,8 @@
 syn keyword juliaConstC BasalforcingsGroundediceMeltingRateEnum
 syn keyword juliaConstC BasalforcingsPerturbationMeltingRateEnum
+syn keyword juliaConstC BasalforcingsSpatialDeepwaterElevationEnum
+syn keyword juliaConstC BasalforcingsSpatialDeepwaterMeltingRateEnum
+syn keyword juliaConstC BasalforcingsSpatialUpperwaterElevationEnum
+syn keyword juliaConstC BasalforcingsSpatialUpperwaterMeltingRateEnum
 syn keyword juliaConstC BasalforcingsIsmip6BasinIdEnum
 syn keyword juliaConstC BasalforcingsIsmip6TfEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26788)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26789)
@@ -595,4 +595,8 @@
 	      else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
+	      else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
+	      else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
+	      else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
+	      else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
 	      else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
@@ -625,12 +629,12 @@
 	      else if (strcmp(name,"BedNorth")==0) return BedNorthEnum;
 	      else if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum;
-	      else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
 	      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 stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
+	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
 	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
 	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
@@ -748,12 +752,12 @@
 	      else if (strcmp(name,"HydrologyGapHeightY")==0) return HydrologyGapHeightYEnum;
 	      else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
-	      else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
 	      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 stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
+	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
 	      else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
 	      else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
@@ -871,12 +875,12 @@
 	      else if (strcmp(name,"SealevelGNrotm2")==0) return SealevelGNrotm2Enum;
 	      else if (strcmp(name,"SealevelGNrotm3")==0) return SealevelGNrotm3Enum;
-	      else if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
 	      else if (strcmp(name,"SealevelGErotm2")==0) return SealevelGErotm2Enum;
 	      else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum;
 	      else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
+	      else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
 	      else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
 	      else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
@@ -994,12 +998,12 @@
 	      else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
 	      else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
-	      else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
 	      else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
 	      else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
 	      else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
+	      else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
 	      else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
 	      else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
@@ -1117,12 +1121,12 @@
 	      else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
 	      else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
-	      else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
 	      else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
 	      else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
 	      else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
+	      else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
 	      else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
 	      else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
@@ -1240,12 +1244,12 @@
 	      else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
 	      else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
-	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
 	      else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
 	      else if (strcmp(name,"IntInput")==0) return IntInputEnum;
 	      else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
+	      else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
 	      else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
 	      else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
@@ -1363,12 +1367,12 @@
 	      else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
 	      else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
-	      else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
 	      else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
 	      else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
 	      else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"Incremental")==0) return IncrementalEnum;
+	      else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
 	      else if (strcmp(name,"Indexed")==0) return IndexedEnum;
 	      else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
@@ -1486,12 +1490,12 @@
 	      else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
 	      else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
-	      else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
 	      else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
 	      else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
 	      else if (strcmp(name,"Regular")==0) return RegularEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
+	      else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
 	      else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
 	      else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
Index: /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m	(revision 26788)
+++ /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m	(revision 26789)
@@ -9,5 +9,7 @@
 		deepwater_melting_rate    = NaN;
 		deepwater_elevation       = NaN;
+		upperwater_melting_rate   = NaN;
 		upperwater_elevation      = NaN;
+		perturbation_melting_rate = NaN;
 		geothermalflux            = NaN;
 	end
@@ -25,5 +27,7 @@
 						self.deepwater_elevation=lb.deepwater_elevation*ones(nvertices,1);
 						self.deepwater_melting_rate=lb.deepwater_melting_rate*ones(nvertices,1);
+						self.upperwater_melting_rate=lb.upperwater_melting_rate*ones(nvertices,1);
 						self.upperwater_elevation=lb.upperwater_elevation*ones(nvertices,1);
+						self.perturbation_melting_rate=lb.perturbation_melting_rate*ones(nvertices,1);
 					else 
 						self=structtoobj(spatiallinearbasalforcings(),varargin{1});
@@ -37,6 +41,8 @@
 			self.deepwater_melting_rate=project3d(md,'vector',self.deepwater_melting_rate,'type','node','layer',1); 
 			self.deepwater_elevation=project3d(md,'vector',self.deepwater_elevation,'type','node','layer',1); 
+			self.upperwater_melting_rate=project3d(md,'vector',self.upperwater_melting_rate,'type','node','layer',1); 
 			self.upperwater_elevation=project3d(md,'vector',self.upperwater_elevation,'type','node','layer',1); 
 			self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
+			self.perturbation_melting_rate=project3d(md,'vector',self.perturbation_melting_rate,'type','node','layer',1); 
 		end % }}}
 		function self = initialize(self,md) % {{{
@@ -58,8 +64,13 @@
 		function md = checkconsistency(self,md,solution,analyses) % {{{
 
+			if numel(md.basalforcings.perturbation_melting_rate)>1
+            md = checkfield(md,'fieldname','basalforcings.perturbation_melting_rate','NaN',1,'Inf',1,'timeseries',1);
+         end
+
 			if ismember('MasstransportAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.ismasstransport==0),
 				md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
 				md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 				md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'timeseries',1);
+				md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 				md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'timeseries',1,'<',0);
 			end
@@ -69,4 +80,5 @@
 				md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1);
+				md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<=',0,'numel',1);
 			end
@@ -76,4 +88,5 @@
 				md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1);
+				md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<',0,'numel',1);
 				md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0);
@@ -86,5 +99,7 @@
 			fielddisplay(self,'deepwater_melting_rate','basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]');
 			fielddisplay(self,'deepwater_elevation','elevation of ocean deepwater [m]');
+			fielddisplay(self,'upperwater_melting_rate','basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]');
 			fielddisplay(self,'upperwater_elevation','elevation of ocean upperwater [m]');
+			fielddisplay(self,'perturbation_melting_rate','basal melting rate perturbation added to computed melting rate (positive if melting) [m/yr]');
 			fielddisplay(self,'geothermalflux','geothermal heat flux [W/m^2]');
 
@@ -94,16 +109,12 @@
 			yts=md.constants.yts;
 
-			floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
-			pos=find(md.geometry.base<=md.basalforcings.deepwater_elevation);
-			floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos);
-			pos=find(md.geometry.base>md.basalforcings.deepwater_elevation & md.geometry.base<md.basalforcings.upperwater_elevation);
-			floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos).*(md.geometry.base(pos)-md.basalforcings.upperwater_elevation(pos))./(md.basalforcings.deepwater_elevation(pos)-md.basalforcings.upperwater_elevation(pos));
 			WriteData(fid,prefix,'name','md.basalforcings.model','data',6,'format','Integer');
-			WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','deepwater_melting_rate','format','DoubleMat','name','md.basalforcings.deepwater_melting_rate','scale',1./yts,'mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','upperwater_melting_rate','format','DoubleMat','name','md.basalforcings.upperwater_melting_rate','scale',1./yts,'mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','perturbation_melting_rate','format','DoubleMat','name','md.basalforcings.perturbation_melting_rate','scale',1./yts,'mattype',1);
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.py	(revision 26788)
+++ /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.py	(revision 26789)
@@ -21,6 +21,8 @@
             self.deepwater_melting_rate     = np.nan
             self.deepwater_elevation        = np.nan
+            self.upperwater_melting_rate    = np.nan
             self.upperwater_elevation       = np.nan
             self.geothermalflux             = np.nan
+            self.perturbation_melting_rate  = np.nan
 
             self.setdefaultparameters()
@@ -30,8 +32,10 @@
                 nvertices = len(lb.groundedice_melting_rate)
                 self.groundedice_melting_rate   = lb.groundedice_melting_rate
-                self.deepwater_melting_rate     = lb.geothermalflux
+                self.geothermalflux             = lb.geothermalflux
                 self.deepwater_elevation        = lb.deepwater_elevation * np.ones((nvertices, ))
-                self.upperwater_elevation       = lb.deepwater_melting_rate * np.ones((nvertices, ))
-                self.geothermalflux             = lb.upperwater_elevation * np.ones((nvertices, ))
+                self.deepwater_melting_rate     = lb.deepwater_melting_rate * np.ones((nvertices, ))
+                self.upperwater_melting_rate    = lb.upperwater_melting_rate * np.ones((nvertices, ))
+                self.upperwater_elevation       = lb.upperwater_elevation * np.ones((nvertices, ))
+                self.perturbation_melting_rate  = lb.perturbation_melting_rate * np.ones((nvertices, ))
             else:
                 # TODO: This has not been tested
@@ -46,5 +50,7 @@
         s += '{}\n'.format(fielddisplay(self, 'deepwater_melting_rate', 'basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]'))
         s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'elevation of ocean deepwater [m]'))
+        s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'))
         s += '{}\n'.format(fielddisplay(self, 'upperwater_elevation', 'elevation of ocean upperwater [m]'))
+        s += '{}\n'.format(fielddisplay(self, 'perturbation_melting_rate', 'basal melting rate perturbation added to computed melting rate (positive if melting) [m/yr]'))
         s += '{}\n'.format(fielddisplay(self, 'geothermalflux', 'geothermal heat flux [W/m^2]'))
         return s
@@ -55,6 +61,8 @@
         self.deepwater_melting_rate = project3d(md, 'vector', self.deepwater_melting_rate, 'type', 'node', 'layer', 1) 
         self.deepwater_elevation = project3d(md, 'vector', self.deepwater_elevation, 'type', 'node', 'layer', 1)
+        self.upperwater_melting_rate = project3d(md, 'vector', self.upperwater_melting_rate, 'type', 'node', 'layer', 1) 
         self.upperwater_elevation = project3d(md, 'vector', self.upperwater_elevation, 'type', 'node', 'layer', 1) 
         self.geothermalflux = project3d(md, 'vector', self.geothermalflux, 'type', 'node', 'layer', 1) # Bedrock only gets geothermal flux
+        self.perturbation_melting_rate = project3d(md, 'vector', self.upperwater_melting_rate, 'type', 'node', 'layer', 1) 
         return self
     #}}}
@@ -72,8 +80,11 @@
 
     def checkconsistency(self, md, solution, analyses): #{{{
+        if not np.all(np.isnan(self.perturbation_melting_rate)):
+            md = checkfield(md, 'fieldname', 'basalforcings.perturbation_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
         if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
             md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
             md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
             md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
             md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '<', 0)
         if 'BalancethicknessAnalysis' in analyses:
@@ -82,4 +93,5 @@
             md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'numel', 1)
             md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'numel', 1)
+            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', '>=', 0, 'numel', 1)
             md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<=', 0, 'numel', 1)
         if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal:
@@ -89,4 +101,5 @@
             md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'numel', 1)
             md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<', 0, 'numel', 1)
+            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'numel', 1)
             md = checkfield(md, 'fieldname', 'basalforcings.geothermalflux', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
         return md
@@ -96,15 +109,11 @@
         yts = md.constants.yts
 
-        floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
-        pos = np.where(md.geometry.base <= md.basalforcings.deepwater_elevation)[0]
-        floatingice_melting_rate[pos] = md.basalforcings.deepwater_melting_rate[pos]
-        pos = np.where(np.logical_and.reduce(md.geometry.base > md.basalforcings.deepwater_elevation, md.geometry.base < md.basalforcings.upperwater_elevation))
-        floatingice_melting_rate[pos] = md.basalforcings.deepwater_melting_rate[pos] * (md.geometry.base[pos] - md.basalforcings.upperwater_elevation[pos]) / (md.basalforcings.deepwater_elevation[pos] - md.basalforcings.upperwater_elevation[pos])
         WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 6, 'format', 'Integer')
-        WriteData(fid, prefix, 'data', floatingice_melting_rate, 'format', 'DoubleMat', 'name', 'md.basalforcings.floatingice_melting_rate', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.groundedice_melting_rate', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1,'yts', md.constants.yts)
         WriteData(fid, prefix,'object', self, 'fieldname', 'geothermalflux', 'name', 'md.basalforcings.geothermalflux', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_melting_rate', 'scale', 1. / yts, 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1. / yts, 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_elevation', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'perturbation_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.perturbation_melting_rate', 'scale', 1. / yts, 'mattype', 1)
     #}}}