Index: ../trunk-jpl/test/NightlyRun/test236.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test236.py	(revision 18967)
+++ ../trunk-jpl/test/NightlyRun/test236.py	(revision 18968)
@@ -15,6 +15,7 @@
 # Use of ispdd and isdelta18o methods
 md.surfaceforcings = SMBpdd();
 md.surfaceforcings.isdelta18o=1
+md.surfaceforcings.ismungsm=0
 
 # Add temperature, precipitation and delta18o needed to measure the surface mass balance
 # creating delta18o
@@ -49,14 +50,20 @@
 for imonth in xrange(0,12):
     md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
     md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
+    md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
+    md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
 
+# Interpolation factors
+md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5;
+md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5;
+# Year of each data point
+md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
+md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
+
 # time steps and resolution
 md.timestepping.time_step=20.
 md.timestepping.final_time=60.
 
-md.surfaceforcings.Pfac=[1,1]
-md.surfaceforcings.Tdiff=[1,1]
-md.surfaceforcings.sealev=[1,1]
 
 # 
 md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures']
Index: ../trunk-jpl/test/NightlyRun/test237.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test237.py	(revision 18967)
+++ ../trunk-jpl/test/NightlyRun/test237.py	(revision 18968)
@@ -11,20 +11,21 @@
 
 md=triangle(model(),'../Exp/Square.exp',600000)    #180000
 md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelf.py')
 
 # Use of ispdd and isdelta18o methods
 md.surfaceforcings = SMBpdd();
-md.surfaceforcings.isdelta18o=1
+md.surfaceforcings.isdelta18o=0
+md.surfaceforcings.ismungsm=1
 
-md=parameterize(md,'../Par/SquareShelf.py')
 
-# Add temperature, precipitation and delta18o needed to measure the surface mass balance
-# creating delta18o
-delta18o=numpy.loadtxt('../Data/delta18o.data')
-md.surfaceforcings.delta18o=delta18o
-# creating delta18oSurface
-md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1)))
-md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:]
+# # Add temperature, precipitation and delta18o needed to measure the surface mass balance
+# # creating delta18o
+# delta18o=numpy.loadtxt('../Data/delta18o.data')
+# md.surfaceforcings.delta18o=delta18o
+# # creating delta18oSurface
+# md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1)))
+# md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:]
 
 # creating Present day and lgm temperatures
 # Same temperature over the all region:
@@ -51,16 +52,23 @@
 for imonth in xrange(0,12):
     md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
     md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
+    md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
+    md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
 
+# Interpolation factors
+md.surfaceforcings.Pfac[1,1:md.timestepping.final_time]=0.5;
+md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5;
+md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5;
+# Year of each data point
+md.surfaceforcings.Pfac[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
+md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
+md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
+
 # time steps and resolution
 md.timestepping.time_step=20.
 md.settings.output_frequency=1
 md.timestepping.final_time=60.
 
-md.surfaceforcings.Pfac=[1,1]
-md.surfaceforcings.Tdiff=[1,1]
-md.surfaceforcings.sealev=[1,1]
-
 #
 md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures']
 md.extrude(3,1.)
Index: ../trunk-jpl/test/NightlyRun/test236.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/test236.m	(revision 18967)
+++ ../trunk-jpl/test/NightlyRun/test236.m	(revision 18968)
@@ -2,12 +2,18 @@
 md=setmask(md,'all','');
 md=parameterize(md,'../Par/SquareShelf.par');
 
+%md.verbose=verbose('all');
+
 % Use of ispdd and isdelta18o methods
 md.surfaceforcings = SMBpdd();
 md.surfaceforcings.isdelta18o=1;
+md.surfaceforcings.ismungsm=0;
 
+%md.surfaceforcings.precipitation(1:md.mesh.numberofvertices,1:12)=0;
+%md.surfaceforcings.monthlytemperatures(1:md.mesh.numberofvertices,1:12)=273;
+
 % Add temperature, precipitation and delta18o needed to measure the surface mass balance
-% creating delta18o
+%  creating delta18o
 load '../Data/delta18o.data'
 md.surfaceforcings.delta18o=delta18o;
 % creating delta18oSurface
@@ -25,7 +31,8 @@
     md.surfaceforcings.temperatures_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
 end
 
-% creating initialization and spc temperatures initialization and spc
+% creating initialization and spc temperatures initialization and
+% spc
 md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1);
 md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step);
 itemp=0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step;
@@ -36,21 +43,28 @@
 % creating precipitation
 for imonth=0:11
     md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
+    md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
+    % Time for the last line:
     md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
+    md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
 end
 
+% Interpolation factors
+md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5;
+md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5;
+% Year of each data point
+md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
+md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
+
 % time steps and resolution
 md.timestepping.time_step=20;
+md.settings.output_frequency=1;
 md.timestepping.final_time=60;
 
-md.surfaceforcings.Pfac=[1;1];
-md.surfaceforcings.Tdiff=[1;1];
-md.surfaceforcings.sealev=[1;1];
-
 % 
 md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'};
 md=setflowequation(md,'SSA','all');
-md.cluster=generic('name',oshostname(),'np',3);
+md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
 md=solve(md,TransientSolutionEnum());
 
 %Fields and tolerances to track changes
Index: ../trunk-jpl/test/NightlyRun/test237.m
===================================================================
--- ../trunk-jpl/test/NightlyRun/test237.m	(revision 18967)
+++ ../trunk-jpl/test/NightlyRun/test237.m	(revision 18968)
@@ -1,19 +1,27 @@
 md=triangle(model(),'../Exp/Square.exp',600000.);%180000
 md=setmask(md,'all','');
+md=parameterize(md,'../Par/SquareShelf.par');
 
-% Use of ispdd and isdelta18o methods
+%md.verbose=verbose('all');
+
+% Use of ispdd methods
 md.surfaceforcings = SMBpdd();
-md.surfaceforcings.isdelta18o=1;
+md.surfaceforcings.isdelta18o=0;
+md.surfaceforcings.ismungsm=1;
 
-md=parameterize(md,'../Par/SquareShelf.par');
+if md.surfaceforcings.isdelta18o==0 & md.surfaceforcings.ismungsm==0
+    md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
+    md.surfaceforcings.monthlytemperatures=273*ones(md.mesh.numberofvertices,1);
+end
+    
 
 % Add temperature, precipitation and delta18o needed to measure the surface mass balance
-% creating delta18o
-load '../Data/delta18o.data'
-md.surfaceforcings.delta18o=delta18o;
-% creating delta18oSurface
-md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0;
-md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:);
+% % creating delta18o
+% load '../Data/delta18o.data'
+% md.surfaceforcings.delta18o=delta18o;
+% % creating delta18oSurface
+% md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0;
+% md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:);
 
 % creating Present day and lgm temperatures
 % Same temperature over the all region:
@@ -37,12 +45,19 @@
 % creating precipitation
 for imonth=0:11
     md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
+    md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
+    % Time for the last line:
     md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
+    md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
 end
 
-md.surfaceforcings.Pfac=[1;1];
-md.surfaceforcings.Tdiff=[1;1];
-md.surfaceforcings.sealev=[1;1];
+md.surfaceforcings.Pfac(1,1:md.timestepping.final_time)=0.5;
+md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5;
+md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5;
+% Year of each data point
+md.surfaceforcings.Pfac(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
+md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
+md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
 
 % time steps and resolution
 md.timestepping.time_step=20;
@@ -51,10 +66,11 @@
 
 %
 md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'};
-md=extrude(md,3,1.);
+md=extrude(md,3,1);
+
 md=setflowequation(md,'SSA','all');
-md.cluster=generic('name',oshostname(),'np',1);
-md=solve(md,TransientSolutionEnum());
+md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
+md=solve(md,TransientSolutionEnum);
 
 %Fields and tolerances to track changes
 field_names     ={'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', ...
Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18968)
@@ -119,7 +119,7 @@
 
 	int    stabilization,finiteelement,smb_model;
 	bool   dakota_analysis;
-	bool   isdelta18o;
+	bool   isdelta18o,ismungsm;
 	bool   isgroundingline;
 	bool   islevelset;
 
@@ -178,17 +178,19 @@
 			break;
 		case SMBpddEnum:
 			iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
+			iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
 			iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
-			if(isdelta18o){
+			if(isdelta18o || ismungsm){
 				iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum);
 			}
 			else{
-				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
+			        iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum);
 			}
+
 			break;
 		case SMBgradientsEnum:
 			iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum);
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18967)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18968)
@@ -343,6 +343,7 @@
 	SurfaceforcingsDelta18oEnum,
 	SurfaceforcingsDelta18oSurfaceEnum,
 	SurfaceforcingsIsdelta18oEnum,
+	SurfaceforcingsIsmungsmEnum,
 	SurfaceforcingsPrecipitationsPresentdayEnum,
 	SurfaceforcingsPrecipitationsLgmEnum,
 	SurfaceforcingsTemperaturesPresentdayEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18968)
@@ -349,6 +349,7 @@
 		case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o";
 		case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface";
 		case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o";
+		case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm";
 		case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday";
 		case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm";
 		case SurfaceforcingsTemperaturesPresentdayEnum : return "SurfaceforcingsTemperaturesPresentday";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18968)
@@ -355,6 +355,7 @@
 	      else if (strcmp(name,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum;
 	      else if (strcmp(name,"SurfaceforcingsDelta18oSurface")==0) return SurfaceforcingsDelta18oSurfaceEnum;
 	      else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum;
+	      else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum;
 	      else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum;
 	      else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum;
 	      else if (strcmp(name,"SurfaceforcingsTemperaturesPresentday")==0) return SurfaceforcingsTemperaturesPresentdayEnum;
@@ -381,11 +382,11 @@
 	      else if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum;
 	      else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
 	      else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum;
-	      else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
+	      if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum;
+	      else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
 	      else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
 	      else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
 	      else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
@@ -504,11 +505,11 @@
 	      else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
 	      else if (strcmp(name,"Masscon")==0) return MassconEnum;
 	      else if (strcmp(name,"MassconName")==0) return MassconNameEnum;
-	      else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
+	      if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
+	      else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
 	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
 	      else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
 	      else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
 	      else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
 	      else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
-	      else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
+	      if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
+	      else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
 	      else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
 	      else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
 	      else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
@@ -750,11 +751,11 @@
 	      else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
 	      else if (strcmp(name,"Gset")==0) return GsetEnum;
 	      else if (strcmp(name,"Index")==0) return IndexEnum;
-	      else if (strcmp(name,"Indexed")==0) return IndexedEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"Intersect")==0) return IntersectEnum;
+	      if (strcmp(name,"Indexed")==0) return IndexedEnum;
+	      else if (strcmp(name,"Intersect")==0) return IntersectEnum;
 	      else if (strcmp(name,"Nodal")==0) return NodalEnum;
 	      else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
 	      else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
Index: ../trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18967)
+++ ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18968)
@@ -12,18 +12,21 @@
 IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n);
 IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat);
 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
-				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
-				IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
-				IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
-				IssmDouble h, IssmDouble s, IssmDouble desfac,
-				IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
-				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime);
+IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
+				 IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
+				 IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
+				 IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
+				 IssmDouble TdiffTime,IssmDouble sealevTime,
+				 IssmDouble rho_water, IssmDouble rho_ice);
 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
-				     IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, 
-				     IssmDouble* PrecipitationsPresentday,
-				     IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
-					  IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
+					     IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime,
+					     IssmDouble* PrecipitationsPresentday,
+					     IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
+					     IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
+void ComputeMungsmTemperaturePrecipitation(IssmDouble TdiffTime, IssmDouble PfacTime,
+					   IssmDouble* PrecipitationsLgm,IssmDouble* PrecipitationsPresentday,
+					   IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
+					   IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
 IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
 
Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18968)
@@ -1,22 +1,18 @@
 /* file:  PddSurfaceMassBlance.cpp
    Calculating the surface mass balance using the positive degree day method.
+   Updating the precipitation and temperature to the new elevation
  */
 
 #include "./elements.h"
 #include "../Numerics/numerics.h"
 
-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, 
-				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
-				IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
-				IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, 
-				IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, 
-				IssmDouble s0t,
-				IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, 
-				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ 
+IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
+				 IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm,
+				 IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,
+				 IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm,
+				 IssmDouble TdiffTime,IssmDouble sealevTime,
+				 IssmDouble rho_water,IssmDouble rho_ice){
 
-  // CURENTLY monthlytemperatures and monthlyprec ARE NOT USED HERE.
-  // THEY ARE REPLACE BY tdiffh No usage of deltO18 anymore
-
   // output:
   IssmDouble B;    // surface mass balance, melt+accumulation
   int    iqj,imonth;
@@ -27,15 +23,15 @@
   IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
   IssmDouble water; //water=rain + snowmelt 
   IssmDouble runoff; //meltwater only, does not include rain 
+  IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
   //IssmDouble  sealev=0.;         // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
   //IssmDouble  Pfac=0.5,Tdiff=0.5;
-  IssmDouble  rtlaps;
+  IssmDouble rtlaps;
   // IssmDouble lapser=6.5         // lapse rate
   // IssmDouble desfac = 0.3;      // desert elevation factor
   // IssmDouble s0p=0.;            // should be set to elevation from precip source
   // IssmDouble s0t=0.;         // should be set to elevation from temperature source
-  IssmDouble tdiffh;  
   IssmDouble st;             // elevation between altitude of the temp record and current altitude
   IssmDouble sp;             // elevation between altitude of the prec record and current altitude
   IssmDouble deselcut=1.0;
@@ -48,7 +44,6 @@
   IssmDouble DT = 0.02;
   IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
 
-  IssmDouble premt;     // monthly precipitation (m of ice equivalent)
   IssmDouble q, qmpt;   // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
   IssmDouble qm = 0.;   // snow part of the precipitation 
   IssmDouble qmt = 0.;  // precipitation without desertification effect adjustment
@@ -72,6 +67,7 @@
   IssmDouble fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
   IssmDouble pddtj, hmx2;
 
+  sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
 
   /*PDD constant*/
   siglim = 2.5*signorm; 
@@ -89,8 +85,8 @@
     rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
     
     /*********compute Surface temperature *******/
-    tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
-    tstar = tdiffh + TemperaturesPresentday[imonth] - rtlaps *max(st,sealevTime*0.001);
+    monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
+    tstar = monthlytemperatures[imonth];
     Tsurf = tstar*deltm+Tsurf;        
 
       /*********compute PD ****************/
@@ -105,10 +101,8 @@
       if (sp>0.0){q = exp(-desfac*sp);}
       else {q = 1.0;}
 
-      premt =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime));   // [m/month]
-
-      qmt= qmt + premt;  
-      qmpt= q*premt;           
+      qmt= qmt + monthlyprec[imonth]*sconv;  //*sconv to convert in m of ice equivalent per month  
+      qmpt= q*monthlyprec[imonth]*sconv;           
       qmp= qmp + qmpt;
       qm= qm + qmpt*pd;
 
@@ -123,6 +117,10 @@
 	pdd = pdd + pddsig*deltm;
 	frzndd = frzndd - (tstar-pddsig)*deltm;}
       else{frzndd = frzndd - tstar*deltm; }
+
+      /*Assign output pointer*/
+      *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
+      *(monthlyprec+imonth) = monthlyprec[imonth];      
   } // end of seasonal loop 
   //******************************************************************
 
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18968)
@@ -21,7 +21,7 @@
 	int         numoutputs,materialtype,smb_model,basalforcing_model;
 	char**      requestedoutputs = NULL;
 	IssmDouble  time;
-	bool        isdelta18o;
+	bool        isdelta18o,ismungsm;
 
 	/*parameters for mass flux:*/
 	int          mass_flux_num_profiles     = 0;
@@ -107,27 +107,31 @@
 			break;
 		case SMBpddEnum:
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0tEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapsEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapslgmEnum));
 			iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
+			iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
 
-			iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2);
-			for(i=0;i<M;i++) temp[M+i]=temp[M+i];
-			parameters->AddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M));
-			iomodel->DeleteData(temp,SurfaceforcingsPfacEnum);
+			if(ismungsm){
+			  iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2);
+			  for(i=0;i<M;i++) temp[M+i]=temp[M+i];
+			  parameters->AddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M));
+			  iomodel->DeleteData(temp,SurfaceforcingsPfacEnum);
 			
-			iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2);
-			for(i=0;i<M;i++) temp[M+i]=temp[M+i];
-			parameters->AddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M));
-			iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum);
+			  iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2);
+			  for(i=0;i<M;i++) temp[M+i]=temp[M+i];
+			  parameters->AddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M));
+			  iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum);
 
-			iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2);
-			for(i=0;i<M;i++) temp[M+i]=temp[M+i];
-			parameters->AddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M));
-			iomodel->DeleteData(temp,SurfaceforcingsSealevEnum);
+			  iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2);
+			  for(i=0;i<M;i++) temp[M+i]=temp[M+i];
+			  parameters->AddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M));
+			  iomodel->DeleteData(temp,SurfaceforcingsSealevEnum);
+			}
 
 			if(isdelta18o){
 				iomodel->Constant(&yts,ConstantsYtsEnum);
Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 18967)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 18968)
@@ -11,6 +11,7 @@
 void SurfaceMassBalancex(FemModel* femmodel);
 void SmbGradientsx(FemModel* femmodel);
 void Delta18oParameterizationx(FemModel* femmodel);
+void MungsmtpParameterizationx(FemModel* femmodel);
 void PositiveDegreeDayx(FemModel* femmodel);
 void SmbHenningx(FemModel* femmodel);
 void SmbComponentsx(FemModel* femmodel);
Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 18968)
@@ -10,7 +10,7 @@
 
 	/*Intermediaties*/
 	int  smb_model;
-	bool isdelta18o;
+	bool isdelta18o,ismungsm;
 
 	/*First, get SMB model from parameters*/
 	femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum);
@@ -22,10 +22,15 @@
 			break;
 		case SMBpddEnum:
 			femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
+			femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum);
 			if(isdelta18o){
-				if(VerboseSolution()) _printf0_("   call Delta18oParametrization module\n");
+				if(VerboseSolution()) _printf0_("   call Delta18oParameterization module\n");
 				Delta18oParameterizationx(femmodel);
 			} 
+			if(ismungsm){
+				if(VerboseSolution()) _printf0_("   call MungsmtpParameterization module\n");
+				MungsmtpParameterizationx(femmodel);
+			} 
 			if(VerboseSolution()) _printf0_("   call positive degree day module\n");
 			PositiveDegreeDayx(femmodel);
 			break;
@@ -117,6 +122,14 @@
 	}
 
 }/*}}}*/
+void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
+
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->MungsmtpParameterization();
+	}
+
+}/*}}}*/
 void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
 
 	// void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
@@ -143,6 +156,8 @@
 	IssmDouble PDup, PDCUT = 2.0;    // PDcut: rain/snow cutoff temperature (C)
 	IssmDouble tstar; // monthly mean surface temp
 
+	bool ismungsm;
+
 	IssmDouble *pdds    = NULL;
 	IssmDouble *pds     = NULL;
 	Element    *element = NULL;
@@ -150,6 +165,9 @@
 	pdds=xNew<IssmDouble>(NPDMAX+1); 
 	pds=xNew<IssmDouble>(NPDCMAX+1); 
 
+	// Get ismungsm parameter
+	femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum);
+
 	/* initialize PDD (creation of a lookup table)*/
 	tstep    = 0.1;
 	tsint    = tstep*0.5;
@@ -204,9 +222,8 @@
 
 	for(i=0;i<femmodel->elements->Size();i++){
 		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->PositiveDegreeDay(pdds,pds,signorm);
+		element->PositiveDegreeDay(pdds,pds,signorm,ismungsm);
 	}
-
 	/*free ressouces: */
 	xDelete<IssmDouble>(pdds);
 	xDelete<IssmDouble>(pds);
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am	(revision 18967)
+++ ../trunk-jpl/src/c/Makefile.am	(revision 18968)
@@ -219,6 +219,7 @@
 					./shared/Elements/PrintArrays.cpp\
 					./shared/Elements/PddSurfaceMassBalance.cpp\
 					./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\
+					./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\
 					./shared/Elements/DrainageFunctionWaterfraction.cpp\
 					./shared/String/sharedstring.h\
 					./shared/String/DescriptorIndex.cpp\
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18968)
@@ -179,6 +179,7 @@
 		virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
 		virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
 		virtual void       Delta18oParameterization(void)=0;
+		virtual void       MungsmtpParameterization(void)=0;
 		virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
 		virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
 		virtual int        FiniteElement(void)=0;
@@ -250,7 +251,7 @@
 		virtual void       NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
 		virtual int        NumberofNodesPressure(void)=0;
 		virtual int        NumberofNodesVelocity(void)=0;
-		virtual void       PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0;
+		virtual void       PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm)=0;
 		virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
 		virtual int        PressureInterpolation()=0;
 		virtual void       ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18968)
@@ -608,7 +608,6 @@
 			input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
 			input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
 			input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
-			PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec
 		}
 	}
 
@@ -650,6 +649,69 @@
 	/*clean-up*/
 	delete gauss;
 }
+void       Tria::MungsmtpParameterization(void){/*{{{*/
+	/*Are we on the base? If not, return*/
+	if(!IsOnBase()) return;
+
+	int        i;
+	IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+	IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+	IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
+	IssmDouble tmp[NUMVERTICES];
+	IssmDouble TdiffTime,PfacTime;
+	IssmDouble time,yts;
+	this->parameters->FindParam(&time,TimeEnum);
+	this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+	/*Recover present day temperature and precipitation*/
+	Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
+	Input*     input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
+	Input*     input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
+	Input*     input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);        _assert_(input4);
+	GaussTria* gauss=new GaussTria();
+	for(int month=0;month<12;month++) {
+		for(int iv=0;iv<NUMVERTICES;iv++) {
+			gauss->GaussVertex(iv);
+			input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
+			input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
+			input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+			input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
+		}
+	}
+
+	/*Recover interpolation parameters at time t*/
+	this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
+	this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
+
+	/*Compute the temperature and precipitation*/
+	for(int iv=0;iv<NUMVERTICES;iv++){
+	  ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
+					&PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0], 
+					&TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
+					&monthlytemperatures[iv][0], &monthlyprec[iv][0]);
+	}
+
+	/*Update inputs*/ 
+	TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
+	TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
+	for (int imonth=0;imonth<12;imonth++) {
+		for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
+		TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
+		NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
+
+		for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
+		TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
+		NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
+	}
+	NewTemperatureInput->Configure(this->parameters);
+	NewPrecipitationInput->Configure(this->parameters);
+
+	this->inputs->AddInput(NewTemperatureInput);
+	this->inputs->AddInput(NewPrecipitationInput);
+
+	/*clean-up*/
+	delete gauss;
+}
 /*}}}*/
 int        Tria::EdgeOnBaseIndex(void){/*{{{*/
 
@@ -2441,83 +2503,94 @@
 	return TriaRef::NumberofNodes(this->VelocityInterpolation());
 }
 /*}}}*/
-void       Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/
-
+void       Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
+ 
+   int        i;
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
-   IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
-   IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
+   IssmDouble tmp[NUMVERTICES];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
-   IssmDouble sconv; //rhow_rain/rhoi / 12 months
    
    /*Get material parameters :*/
    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-   //rho_ice=matpar->GetRhoIce();
-   //rho_water=matpar->GetRhoFreshwater();
 
-   sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
+  /*Get some pdd parameters*/
+  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
+  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
+  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
+  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
+  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
 
    /*Recover monthly temperatures and precipitation and Present day and LGm ones*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
-   Input*     input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input3);
-   Input*     input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input4);
-   Input*     input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum);  _assert_(input5);
-   Input*     input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);   _assert_(input6);
    GaussTria* gauss=new GaussTria();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
    this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+
    for(int month=0;month<12;month++) {
      for(int iv=0;iv<NUMVERTICES;iv++) {
        gauss->GaussVertex(iv);
        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
-       monthlyprec[iv][month]=monthlyprec[iv][month]*sconv*yts; // convertion in m/y
-       input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
-       TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
-       input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
-       TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
-       input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
-       PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; 
-       input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
-       PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];  
      }
    }
 
-  /*Recover Pfac, Tdiff and sealev at time t:*/
-  this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
-  this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
-  this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
+  /*Recover Pfac, Tdiff and sealev at time t:
+    This parameters are used to interpolate the temperature 
+    and precipitaton between PD and LGM when ismungsm==1 */ 
+  if (ismungsm==1){  
+    this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
+    this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
+  }
+  else {
+    TdiffTime=0;
+    sealevTime=0;  
+  }
 
   /*Recover info at the vertices: */
   GetInputListOnVertices(&h[0],ThicknessEnum);
   GetInputListOnVertices(&s[0],SurfaceEnum);
-
-  /*Get other pdd parameters*/
-  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
-  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
-  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
-  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
-  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
      
    /*measure the surface mass balance*/
-   for (int iv = 0; iv<NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
-				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
-				  &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], 
+  for (int iv = 0; iv<NUMVERTICES; iv++){
+     agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
 				  pdds, pds, signorm, yts, h[iv], s[iv],
-				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
+				  desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
+				  rho_water,rho_ice);
    }
 
    /*Update inputs*/    
+   // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
+   // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
+   // for (int imonth=0;imonth<12;imonth++) {
+   //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
+   //   TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
+   //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
+   // 
+   //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
+   //   TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
+   //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
+   // }
+   // NewTemperatureInput->Configure(this->parameters);
+   // NewPrecipitationInput->Configure(this->parameters);
+
+
    this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
+   // this->inputs->AddInput(NewTemperatureInput);
+   // this->inputs->AddInput(NewPrecipitationInput);
    // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
 
+   //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
+   // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
+   // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
+
 	/*clean-up*/
 	delete gauss;
 }
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18968)
@@ -62,6 +62,7 @@
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
 		void        Delta18oParameterization(void);
+		void        MungsmtpParameterization(void);
 		int         EdgeOnBaseIndex();
 		void        EdgeOnBaseIndices(int* pindex1,int* pindex);
 		int         EdgeOnSurfaceIndex();
@@ -108,7 +109,7 @@
 		int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
 		int         NumberofNodesPressure(void);
 		int         NumberofNodesVelocity(void);
-		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
+		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		int         PressureInterpolation();
 		void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18968)
@@ -607,7 +607,6 @@
 			input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
 			input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
 			input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
-			PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec
 		}
 	}
 
@@ -652,6 +651,72 @@
 	/*clean-up*/
 	delete gauss;
 }
+void       Penta::MungsmtpParameterization(void){/*{{{*/
+	/*Are we on the base? If not, return*/
+	if(!IsOnBase()) return;
+
+	int        i;
+	IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+	IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+	IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
+	IssmDouble tmp[NUMVERTICES];
+	IssmDouble TdiffTime,PfacTime;
+	IssmDouble time,yts;
+	this->parameters->FindParam(&time,TimeEnum);
+	this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+	/*Recover present day temperature and precipitation*/
+	Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
+	Input*     input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
+	Input*     input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
+	Input*     input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);        _assert_(input4);
+	GaussPenta* gauss=new GaussPenta();
+	for(int month=0;month<12;month++) {
+		for(int iv=0;iv<NUMVERTICES;iv++) {
+			gauss->GaussVertex(iv);
+			input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
+			input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
+			input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+			input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
+		}
+	}
+
+	/*Recover interpolation parameters at time t*/
+	this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
+	this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
+
+	/*Compute the temperature and precipitation*/
+	for(int iv=0;iv<NUMVERTICES;iv++){
+	  ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
+					&PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0], 
+					&TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
+					&monthlytemperatures[iv][0], &monthlyprec[iv][0]);
+	}
+
+	/*Update inputs*/ 
+	TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
+	TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
+	for (int imonth=0;imonth<12;imonth++) {
+		for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
+		PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
+		NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
+	
+		for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
+		PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
+		NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
+	}
+	NewTemperatureInput->Configure(this->parameters);
+	NewPrecipitationInput->Configure(this->parameters);
+
+	this->inputs->AddInput(NewTemperatureInput);
+	this->inputs->AddInput(NewPrecipitationInput);
+	
+	this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
+	this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
+
+	/*clean-up*/
+	delete gauss;
+}
 /*}}}*/
 void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
 
@@ -2081,32 +2146,30 @@
 
 }
 /*}}}*/
-void       Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/
+void       Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
 
+   int        i;
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
-   IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
-   IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
+   IssmDouble tmp[NUMVERTICES];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
-   IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
    /*Get material parameters :*/
    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-   //rho_ice=matpar->GetRhoIce();
-   //rho_water=matpar->GetRhoFreshwater();
 
-   sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
+  /*Get some pdd parameters*/
+  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
+  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
+  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
+  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
+  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
 
    /*Recover monthly temperatures and precipitation*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
-   Input*     input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);  _assert_(input3);
-   Input*     input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);   _assert_(input4);
-   Input*     input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum);  _assert_(input5);
-   Input*     input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);   _assert_(input6);
    GaussPenta* gauss=new GaussPenta();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
@@ -2118,50 +2181,61 @@
        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
-       monthlyprec[iv][month]=monthlyprec[iv][month]*sconv; // convertion to m/yr
-       input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
-       TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
-       input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
-       TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
-       input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
-       PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; 
-       input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
-       PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];  
      }
    } 
 
-  /*Recover Pfac, Tdiff and sealev at time t:*/
-  this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
-  this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
-  this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
+  /*Recover Pfac, Tdiff and sealev at time t:
+    This parameters are used to interpolate the temperature 
+    and precipitaton between PD and LGM when ismungsm==1 */ 
+  if (ismungsm==1){  
+    this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
+    this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
+  }
+  else {
+    TdiffTime=0;
+    sealevTime=0;  
+  }
 
   /*Recover info at the vertices: */
   GetInputListOnVertices(&h[0],ThicknessEnum);
   GetInputListOnVertices(&s[0],SurfaceEnum); 
 
-  /*Get other pdd parameters*/
-  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
-  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
-  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
-  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
-  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
 
    /*measure the surface mass balance*/
    for (int iv = 0; iv < NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],  
-				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
-				  &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], 
+     agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],  
 				  pdds,pds, signorm, yts, h[iv], s[iv],
-				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
+				  desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
+				  rho_water,rho_ice);
    }
 
    /*Update inputs*/    
+   // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
+   // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
+   // for (int imonth=0;imonth<12;imonth++) {
+   //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
+   //   PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
+   //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
+   // 
+   //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
+   //   PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
+   //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
+   // }
+   // NewTemperatureInput->Configure(this->parameters);
+   // NewPrecipitationInput->Configure(this->parameters);
+
+
+
    this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
-   //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
-   this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
+   // this->inputs->AddInput(NewTemperatureInput);
+   // this->inputs->AddInput(NewPrecipitationInput);
+   // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
+    this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
+   // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
+   // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
 
-	/*clean-up*/
-	delete gauss;
+   /*clean-up*/
+   delete gauss;
 }
 /*}}}*/
 void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18968)
@@ -58,6 +58,7 @@
 		void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
 		ElementMatrix* CreateBasalMassMatrix(void);
 		void           Delta18oParameterization(void);
+		void           MungsmtpParameterization(void);
 		void           ElementResponse(IssmDouble* presponse,int response_enum);
 		void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
 		int            FiniteElement(void);
@@ -136,7 +137,7 @@
 		int            NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
 		int            NumberofNodesPressure(void);
 		int            NumberofNodesVelocity(void);
-		void           PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
+		void           PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
 		void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		int            PressureInterpolation();
 		void           ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18968)
@@ -53,6 +53,7 @@
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
 		void        Delta18oParameterization(void){_error_("not implemented yet");};
+		void        MungsmtpParameterization(void){_error_("not implemented yet");};
 		void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
 		void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
 		void        FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){_error_("not implemented yet");};
@@ -128,7 +129,7 @@
 		void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
 		int         NumberofNodesPressure(void){_error_("not implemented yet");};
 		int         NumberofNodesVelocity(void){_error_("not implemented yet");};
-		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
+		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
 		int         PressureInterpolation(void){_error_("not implemented yet");};
 		void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18967)
+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18968)
@@ -53,6 +53,7 @@
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
 		void        Delta18oParameterization(void){_error_("not implemented yet");};
+		void        MungsmtpParameterization(void){_error_("not implemented yet");};
 		IssmDouble  DragCoefficientAbsGradient(void){_error_("not implemented yet");};
 		void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
 		void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
@@ -134,7 +135,7 @@
 		void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
 		int         NumberofNodesPressure(void);
 		int         NumberofNodesVelocity(void);
-		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
+		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
 		int         PressureInterpolation(void);
 		void        ResetFSBasalBoundaryCondition(void);
Index: ../trunk-jpl/src/m/classes/SMBpdd.py
===================================================================
--- ../trunk-jpl/src/m/classes/SMBpdd.py	(revision 18967)
+++ ../trunk-jpl/src/m/classes/SMBpdd.py	(revision 18968)
@@ -25,11 +25,13 @@
 		self.Tdiff                     = float('NaN')
 		self.sealev                    = float('NaN')
 		self.isdelta18o                = 0
+		self.ismungsm                  = 0
 		self.delta18o                  = float('NaN')
 		self.delta18o_surface          = float('NaN')
 		self.temperatures_presentday   = float('NaN')
 		self.temperatures_lgm          = float('NaN')
 		self.precipitations_presentday = float('NaN')
+		self.precipitations_lgm        = float('NaN')
 
 		#set defaults
 		self.setdefaultparameters()
@@ -37,46 +39,63 @@
 	def __repr__(self): # {{{
 		string="   surface forcings parameters:"
 
-		string="%s\n%s"%(string,fielddisplay(self,'precipitation','surface precipitation [m/yr water eq]'))
+		string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'))
+		string="%s\n%s"%(string,fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)'))
 		string="%s\n%s"%(string,fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]'))
-		string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'))
 		string="%s\n%s"%(string,fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
 		string="%s\n%s"%(string,fielddisplay(self,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
 		string="%s\n%s"%(string,fielddisplay(self,'rlaps','present day lapse rate [degree/km]'))
 		string="%s\n%s"%(string,fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]'))
-		string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D (year)'))
-		string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D (year)'))
-		string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year)'))
-		string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated'))
-		string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated'))
-		string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated'))
-		string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'))
-		string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated'))
-		string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'))
-
+		if not (self.isdelta18o and self.ismungsm):
+			string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']))
+			string="%s\n%s"%(string,fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']))
+			if self.isdelta18o:
+				string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'))
+			if self.ismungsm:
+				string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'))
+				string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'))
 		return string
 		#}}}
 	def extrude(self,md): # {{{
 
-		self.precipitation=project3d(md,'vector',self.precipitation,'type','node');
-		self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');
+		if not (self.isdelta18o and self.ismungsm):
+			self.precipitation=project3d(md,'vector',self.precipitation,'type','node')
+			self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node')
 		if self.isdelta18o: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node')
 		if self.isdelta18o: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node')
 		if self.isdelta18o: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node')
+		if self.isdelta18o: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node')
+		if self.ismungsm: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node')
+		if self.ismungsm: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node')
+		if self.ismungsm: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node')
+		if self.ismungsm: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node')
 		return self
 	#}}}
 	def initialize(self,md): # {{{
 
-		if numpy.all(numpy.isnan(self.precipitation)):
-			self.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-			print "      no SMBpdd.precipitation specified: values set as zero"
-
-		return self
+		# if numpy.all(numpy.isnan(self.precipitation)):
+		# 	self.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
+		# 	print "      no SMBpdd.precipitation specified: values set as zero"
+		# 
+		# return self
 	#}}}
 	def setdefaultparameters(self): # {{{
 		  
 		#pdd method not used in default mode
 		self.isdelta18o = 0
+		self.ismungsm   = 0
 		self.desfac     = 0.5
 		self.s0p        = 0.
 		self.s0t        = 0.
@@ -88,23 +107,32 @@
 	def checkconsistency(self,md,solution,analyses):    # {{{
 
 		if MasstransportAnalysisEnum() in analyses:
-			md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]);
-			md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]);
-			md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]);
-			md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]);
-			md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]);
-			md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]);
-			md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]);
-			md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]);
-			if not self.isdelta18o:
-				md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1)
-				md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1)
-			else:
+			md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1])
+			md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1])
+			md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1])
+			md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1])
+			md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1])
+
+			if not (self.isdelta18o and self.ismungsm):
+				md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','NaN',1)
+				md = checkfield(md,'fieldname','surfaceforcings.precipitation','NaN',1)
+			elif self.isdelta18o:
 				md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1)
 				md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1)
 				md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
 				md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1)
 				md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
+				md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)                                       
+				md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan])
+				md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan])
+			elif self.ismungsm:
+				md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1)
+				md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)
+				md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1)
+				md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)                                       
+				md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan])
+				md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan])
+				md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan])
 
 		return md
 	# }}}
@@ -112,24 +140,34 @@
 
 		yts=365.0*24.0*3600.0
 
-		WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer');
+		WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer')
 
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
 		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean')
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double');
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double');
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double');
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double');
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double');
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
-		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
-		if self.isdelta18o:
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean')
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double')
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double')
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double')
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double')
+		WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double')
+
+		if not (self.isdelta18o and self.ismungsm):
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
+		elif self.isdelta18o:
 			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1)
 			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1)
 			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1)
 			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1)
 			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1)
-		else:
-			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1)			
+		elif self.ismungsm:
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1)
+			WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1)
 	# }}}
Index: ../trunk-jpl/src/m/classes/SMBpdd.m
===================================================================
--- ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18967)
+++ ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18968)
@@ -16,6 +16,7 @@
                 Tdiff                     = NaN;
                 sealev                    = NaN;
 		isdelta18o                = 0;
+		ismungsm                  = 0;
 		delta18o                  = NaN;
 		delta18o_surface          = NaN;
 		temperatures_presentday   = NaN;
@@ -33,27 +34,30 @@
 			end
 		end % }}}
 		function self = extrude(self,md) % {{{
-
-			self.precipitation=project3d(md,'vector',self.precipitation,'type','node');
-			self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');
+			if(self.isdelta18o==0 & self.ismungsm==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end
+			if(self.isdelta18o==0 & self.ismungsm==0),self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');end
 			if(self.isdelta18o),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end
 			if(self.isdelta18o),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
 			if(self.isdelta18o),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
 			if(self.isdelta18o),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end
+			if(self.ismungsm),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end
+			if(self.ismungsm),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
+			if(self.ismungsm),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
+			if(self.ismungsm),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end
 
-
 		end % }}}
 		function self = initialize(self,md) % {{{
+                    
+			% if isnan(self.precipitation),
+			% 	self.precipitation=zeros(md.mesh.numberofvertices,1);
+			% 	disp('      no SMBpdd.precipitation specified: values set as zero');
+			% end
 
-			if isnan(self.precipitation),
-				self.precipitation=zeros(md.mesh.numberofvertices,1);
-				disp('      no SMBpdd.precipitation specified: values set as zero');
-			end
-
 		end % }}}
 		function obj = setdefaultparameters(obj) % {{{
 
 		  obj.isdelta18o = 0;
+		  obj.ismungsm   = 0;
 		  obj.desfac     = 0.5;
 		  obj.s0p        = 0;
 		  obj.s0t        = 0;
@@ -69,44 +73,61 @@
 				md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1);
 				md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1);
-            md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]);
-            md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,NaN]);
-            md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,NaN]);
-				if(obj.isdelta18o==0)
-					md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1);
-				else
-					md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
-					md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+				if(obj.isdelta18o==0 & obj.ismungsm==0)
+                                    md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1);
+                                elseif(obj.isdelta18o==1) 
+                                    md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);                                       
+                                    md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
+                                elseif(obj.ismungsm==1) 
+                                    md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);                                       
+                                    md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]);
+                                    md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
+                                    md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
 				end
-			end
+                        end
 		end % }}}
 		function disp(obj) % {{{
 			disp(sprintf('   surface forcings parameters:'));
 
 			disp(sprintf('\n   PDD and deltaO18 parameters:'));
 			fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)');
+			fielddisplay(obj,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)');
 			fielddisplay(obj,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]');
 			fielddisplay(obj,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]');
 			fielddisplay(obj,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]');
 			fielddisplay(obj,'rlaps','present day lapse rate [degree/km]');
 			fielddisplay(obj,'rlapslgm','LGM lapse rate [degree/km]');
-			fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year)');
-			fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year)');
-			fielddisplay(obj,'sealev','sea level [m], 1D(year)');
-			fielddisplay(obj,'monthlytemperatures',['CURRENTLY NOT USED monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
-			fielddisplay(obj,'precipitation',['CURRENTLY NOT USED monthly surface precipitation [m/yr water eq]']);
-			fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated');
-			fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated');
-			fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
-			fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
-			fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated');
-			fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated');
-
+                        if(obj.isdelta18o==0 & obj.ismungsm==0)
+                            fielddisplay(obj,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
+                            fielddisplay(obj,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']);
+                        elseif(obj.isdelta18o==1)
+                            fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated');
+                            fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated');
+                            fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
+                            fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
+                            fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated');
+                        elseif(obj.ismungsm==1)
+                            fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
+                            fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated');
+                            fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
+                            fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated');
+                        end
 		end % }}}
 		function marshall(obj,md,fid) % {{{
 
@@ -114,25 +135,35 @@
 
 			WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer');
 
-			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','ismungsm','format','Boolean');
 			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','desfac','format','Double');
 			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0p','format','Double');
 			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0t','format','Double');
 			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlaps','format','Double');
 			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlapslgm','format','Double');
-			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean');
-			if obj.isdelta18o
+
+                        if(obj.isdelta18o==0 & obj.ismungsm==0)
+                            %WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
+                            WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+                            WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+			elseif obj.isdelta18o
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1);
-			else
-				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+                                WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
+			elseif obj.ismungsm
+				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
+				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
+				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
+				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
 			end
 		end % }}}
 	end
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18967)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18968)
@@ -341,6 +341,7 @@
 def SurfaceforcingsDelta18oEnum(): return StringToEnum("SurfaceforcingsDelta18o")[0]
 def SurfaceforcingsDelta18oSurfaceEnum(): return StringToEnum("SurfaceforcingsDelta18oSurface")[0]
 def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0]
+def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0]
 def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0]
 def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0]
 def SurfaceforcingsTemperaturesPresentdayEnum(): return StringToEnum("SurfaceforcingsTemperaturesPresentday")[0]
