Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18717)
@@ -330,6 +330,12 @@
 	SurfaceforcingsPrecipitationEnum,
 	SurfaceforcingsDesfacEnum,
 	SurfaceforcingsS0pEnum,
+	SurfaceforcingsS0tEnum,
+	SurfaceforcingsRlapsEnum,
+	SurfaceforcingsRlapslgmEnum,
+	SurfaceforcingsPfacEnum,
+	SurfaceforcingsTdiffEnum,
+	SurfaceforcingsSealevEnum,
 	SMBgradientsEnum,
 	SurfaceforcingsMonthlytemperaturesEnum,
 	SurfaceforcingsHrefEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18717)
@@ -336,6 +336,12 @@
 		case SurfaceforcingsPrecipitationEnum : return "SurfaceforcingsPrecipitation";
 		case SurfaceforcingsDesfacEnum : return "SurfaceforcingsDesfac";
 		case SurfaceforcingsS0pEnum : return "SurfaceforcingsS0p";
+		case SurfaceforcingsS0tEnum : return "SurfaceforcingsS0t";
+		case SurfaceforcingsRlapsEnum : return "SurfaceforcingsRlaps";
+		case SurfaceforcingsRlapslgmEnum : return "SurfaceforcingsRlapslgm";
+		case SurfaceforcingsPfacEnum : return "SurfaceforcingsPfac";
+		case SurfaceforcingsTdiffEnum : return "SurfaceforcingsTdiff";
+		case SurfaceforcingsSealevEnum : return "SurfaceforcingsSealev";
 		case SMBgradientsEnum : return "SMBgradients";
 		case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures";
 		case SurfaceforcingsHrefEnum : return "SurfaceforcingsHref";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18717)
@@ -342,6 +342,12 @@
 	      else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
 	      else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
 	      else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
+	      else if (strcmp(name,"SurfaceforcingsS0t")==0) return SurfaceforcingsS0tEnum;
+	      else if (strcmp(name,"SurfaceforcingsRlaps")==0) return SurfaceforcingsRlapsEnum;
+	      else if (strcmp(name,"SurfaceforcingsRlapslgm")==0) return SurfaceforcingsRlapslgmEnum;
+	      else if (strcmp(name,"SurfaceforcingsPfac")==0) return SurfaceforcingsPfacEnum;
+	      else if (strcmp(name,"SurfaceforcingsTdiff")==0) return SurfaceforcingsTdiffEnum;
+	      else if (strcmp(name,"SurfaceforcingsSealev")==0) return SurfaceforcingsSealevEnum;
 	      else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
 	      else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
 	      else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum;
@@ -379,9 +385,6 @@
 	      else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
 	      else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
 	      else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
-	      else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
-	      else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
-	      else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
          else stage=4;
    }
    if(stage==4){
@@ -502,9 +505,6 @@
 	      else if (strcmp(name,"Ice")==0) return IceEnum;
 	      else if (strcmp(name,"Melange")==0) return MelangeEnum;
 	      else if (strcmp(name,"Water")==0) return WaterEnum;
-	      else if (strcmp(name,"Closed")==0) return ClosedEnum;
-	      else if (strcmp(name,"Free")==0) return FreeEnum;
-	      else if (strcmp(name,"Open")==0) return OpenEnum;
          else stage=5;
    }
    if(stage==5){
@@ -625,9 +625,6 @@
 	      else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
 	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
-	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
-	      else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
-	      else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
          else stage=6;
    }
    if(stage==6){
@@ -748,9 +745,6 @@
 	      else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
 	      else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
 	      else if (strcmp(name,"SeaiceCoriolisFactor")==0) return SeaiceCoriolisFactorEnum;
-	      else if (strcmp(name,"BasalforcingsRhoOcean")==0) return BasalforcingsRhoOceanEnum;
-	      else if (strcmp(name,"BasalforcingsOceanCoef")==0) return BasalforcingsOceanCoefEnum;
-	      else if (strcmp(name,"BasalforcingsOceanLinDragCoef")==0) return BasalforcingsOceanLinDragCoefEnum;
          else stage=7;
    }
    if(stage==7){
Index: ../trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18717)
@@ -12,9 +12,12 @@
 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* pdds, IssmDouble* pds, 
-				IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s,
-				IssmDouble rho_ice, IssmDouble rho_water, IssmDouble desfac, IssmDouble s0p);
+IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
+				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
+				IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
+				IssmDouble h, IssmDouble s, IssmDouble rho_ice,IssmDouble rho_water, IssmDouble desfac,
+				IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
+				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime);
 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
 				     IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, 
 				     IssmDouble* PrecipitationsPresentday,
Index: ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp	(revision 18717)
@@ -19,7 +19,7 @@
   glacialindex = 0.;//(Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))
   //  /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent));
 
-  for (int imonth = 0; imonth<12; imonth++){  
+  for (int imonth = 0; imonth<12; imonth++){
     monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1.-glacialindex)*TemperaturesPresentday[imonth];
     //monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-273.15-0)/10)); //equation from rob's paper
     monthlyprectmp[imonth] = PrecipitationsPresentday[imonth];
Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18717)
@@ -5,7 +5,13 @@
 #include "./elements.h"
 #include "../Numerics/numerics.h"
 
-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water, IssmDouble desfac, IssmDouble s0p){ 
+IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, 
+				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
+				IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, 
+				IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, 
+				IssmDouble rho_water, IssmDouble desfac, IssmDouble s0t,
+				IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, 
+				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ 
 
   // output:
   IssmDouble B;    // surface mass balance, melt+accumulation
@@ -19,10 +25,14 @@
   IssmDouble runoff; //meltwater only, does not include rain 
   IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
-  IssmDouble lapser=6.5, sealev=0.;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
-  // IssmDouble desfac = 0.5;                 // 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  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 lapser=6.5         // lapse rate
+  // IssmDouble desfac = 0.5;      // 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
 
@@ -67,13 +77,18 @@
   PDup = siglimc+PDCUT;
 
   // seasonal loop
-    for (iqj = 0; iqj < 12; iqj++){
-      imonth =  ismon[iqj];
+  for (iqj = 0; iqj < 12; iqj++){
+    imonth =  ismon[iqj];
+    
+    /*********compute lapse rate ****************/
+    st=(s-s0t)/1000.;
+    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);
+    Tsurf = tstar*deltm+Tsurf;        
 
-      st=(s-s0t)/1000.;
-      tstar = monthlytemperatures[imonth] - lapser *max(st,sealev);
-      Tsurf = tstar*deltm+Tsurf;        
-
       /*********compute PD ****************/
       if (tstar < PDup){
 	pd = 1.;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18717)
@@ -36,6 +36,10 @@
 	IssmDouble  *matrix                     = NULL;
 	int          count;
 
+	IssmDouble *temp = NULL;
+	IssmDouble  yts;
+	int         N,M;
+
 	/*Make sure current dataset is empty*/
 	_assert_(parameters->Size()==0); 
 
@@ -103,12 +107,29 @@
 			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
 			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->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,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){
-				IssmDouble *temp = NULL;
-				IssmDouble  yts;
-				int         N,M;
 				iomodel->Constant(&yts,ConstantsYtsEnum);
+
 				iomodel->FetchData(&temp,&N,&M,SurfaceforcingsDelta18oEnum); _assert_(N==2);
 				for(i=0;i<M;i++) temp[M+i]=yts*temp[M+i];
 				parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oEnum,&temp[0],&temp[M],M));
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18717)
@@ -51,6 +51,9 @@
 				case SMBpddEnum:
 					iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
 					iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
+					iomodel->Constant(&this->s0t,SurfaceforcingsS0tEnum);
+					iomodel->Constant(&this->rlaps,SurfaceforcingsRlapsEnum);
+					iomodel->Constant(&this->rlapslgm,SurfaceforcingsRlapslgmEnum);
 					break;
 				case SMBgradientsEnum:
 					/*Nothing to add*/
@@ -141,6 +144,9 @@
 	_printf_("   g: " << g << "\n");
 	_printf_("   desfac: " << desfac << "\n");
 	_printf_("   s0p: " << s0p << "\n");
+	_printf_("   s0t: " << s0t << "\n");
+	_printf_("   rlaps: " << rlaps << "\n");
+	_printf_("   rlapslgm: " << rlapslgm << "\n");
 	return;
 }
 /*}}}*/
@@ -183,6 +189,9 @@
 	matpar->g=this->g;
 	matpar->desfac=this->desfac;
 	matpar->s0p=this->s0p;
+	matpar->s0t=this->s0t;
+	matpar->rlaps=this->rlaps;
+	matpar->rlapslgm=this->rlapslgm;
 
 	matpar->sediment_compressibility=this->sediment_compressibility;
 	matpar->sediment_porosity=this->sediment_porosity;
@@ -269,6 +278,15 @@
 		case SurfaceforcingsS0pEnum:
 			this->s0p=constant;
 			break;
+		case SurfaceforcingsS0tEnum:
+			this->s0t=constant;
+			break;
+		case SurfaceforcingsRlapsEnum:
+			this->rlaps=constant;
+			break;
+		case SurfaceforcingsRlapslgmEnum:
+			this->rlapslgm=constant;
+			break;
 		default: 
 			break;
 	}
@@ -366,6 +384,18 @@
 	return s0p;
 }
 /*}}}*/
+IssmDouble Matpar::GetS0t(){/*{{{*/
+	return s0t;
+}
+/*}}}*/
+IssmDouble Matpar::GetRlaps(){/*{{{*/
+	return rlaps;
+}
+/*}}}*/
+IssmDouble Matpar::GetRlapslgm(){/*{{{*/
+	return rlapslgm;
+}
+/*}}}*/
 IssmDouble Matpar::TMeltingPoint(IssmDouble pressure){/*{{{*/
 	return meltingpoint-beta*pressure;
 }
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18716)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18717)
@@ -32,6 +32,9 @@
 		IssmDouble  omega;
 		IssmDouble  desfac;
 		IssmDouble  s0p;
+		IssmDouble  s0t;
+		IssmDouble  rlaps;
+		IssmDouble  rlapslgm;
 
 		/*hydrology Dual Porous Continuum: */	 
 		IssmDouble  sediment_compressibility;
@@ -124,6 +127,9 @@
 		void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
 		IssmDouble GetDesFac();
 		IssmDouble GetS0p(); 
+		IssmDouble GetS0t(); 
+		IssmDouble GetRlaps(); 
+		IssmDouble GetRlapslgm(); 
 		IssmDouble GetMaterialParameter(int in_enum); 
 		/*}}}*/
 
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18717)
@@ -1720,12 +1720,16 @@
 
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+   IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
-   IssmDouble rho_water,rho_ice,desfac,s0p;
+   IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
+   IssmDouble PfacTime,TdiffTime,sealevTime;
 
-   /*Recover monthly temperatures and precipitation*/
+   /*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);
    GaussTria* gauss=new GaussTria();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
@@ -1737,9 +1741,18 @@
        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]*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
      }
    }
 
+  /*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 info at the vertices: */
   GetInputListOnVertices(&h[0],ThicknessEnum);
   GetInputListOnVertices(&s[0],SurfaceEnum);
@@ -1748,13 +1761,19 @@
   rho_ice=matpar->GetRhoIce();
   rho_water=matpar->GetRhoFreshwater();
 
-  /*Get desertification effect parameters*/
+  /*Get other pdd parameters*/
   desfac=matpar->GetDesFac();
   s0p=matpar->GetS0p();
-
+  s0t=matpar->GetS0t();
+  rlaps=matpar->GetRlaps();
+  rlapslgm=matpar->GetRlapslgm();
+     
    /*measure the surface mass balance*/
    for (int iv = 0; iv<NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water, desfac, s0p);
+     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
+				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
+				  pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water, 
+				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
    }
 
    /*Update inputs*/    
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18716)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18717)
@@ -1915,12 +1915,16 @@
 
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
-   IssmDouble h[NUMVERTICES],s[NUMVERTICES]; // ,b
-   IssmDouble rho_water,rho_ice,desfac,s0p;
+   IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+   IssmDouble h[NUMVERTICES],s[NUMVERTICES];
+   IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
+   IssmDouble PfacTime,TdiffTime,sealevTime;
 
    /*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);
    GaussPenta* gauss=new GaussPenta();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
@@ -1932,9 +1936,18 @@
        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]*yts; // 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
      }
    } 
 
+  /*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 info at the vertices: */
   GetInputListOnVertices(&h[0],ThicknessEnum);
   GetInputListOnVertices(&s[0],SurfaceEnum);
@@ -1943,14 +1956,19 @@
   rho_ice=matpar->GetRhoIce();
   rho_water=matpar->GetRhoFreshwater();
 
-  /*Get desertification effect parameters*/
+  /*Get other pdd parameters*/
   desfac=matpar->GetDesFac();
   s0p=matpar->GetS0p();
+  s0t=matpar->GetS0t();
+  rlaps=matpar->GetRlaps();
+  rlapslgm=matpar->GetRlapslgm();
 
    /*measure the surface mass balance*/
    for (int iv = 0; iv < NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, 
-				  signorm, yts, h[iv], s[iv], rho_ice, rho_water, desfac, s0p);
+     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],  
+				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
+				  pdds,pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water,
+				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
    }
 
    /*Update inputs*/    
Index: ../trunk-jpl/src/m/classes/SMBpdd.m
===================================================================
--- ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18716)
+++ ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18717)
@@ -9,6 +9,12 @@
 		monthlytemperatures       = NaN;
 		desfac                    = 0.5;
 		s0p                       = 0;
+		s0t                       = 0;
+                rlaps                     = 6.5;
+                rlapslgm                  = 6.5;                
+                Pfac                      = NaN;
+                Tdiff                     = NaN;
+                sealev                    = NaN;
 		isdelta18o                = 0;
 		delta18o                  = NaN;
 		delta18o_surface          = NaN;
@@ -48,12 +54,22 @@
 		  obj.isdelta18o = 0;
 		  obj.desfac     = 0.5;
 		  obj.s0p        = 0;
+		  obj.s0t        = 0;
+		  obj.rlaps      = 6.5;
+		  obj.rlapslgm   = 6.5;
+                  
 		end % }}}
 		function md = checkconsistency(obj,md,solution,analyses) % {{{
 
 			if ismember(MasstransportAnalysisEnum(),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);
+                                md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
+                                md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
 				if(obj.isdelta18o==0)
 					md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
 					md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1);
@@ -73,7 +89,13 @@
 			fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o 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,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated');
+			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','monthly surface temperatures [K], required if pdd is activated and delta18o not activated');
 			fielddisplay(obj,'precipitation','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');
@@ -91,6 +113,12 @@
 			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','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
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
Index: ../trunk-jpl/src/m/classes/clusters/acenet.m
===================================================================
--- ../trunk-jpl/src/m/classes/clusters/acenet.m	(revision 18716)
+++ ../trunk-jpl/src/m/classes/clusters/acenet.m	(revision 18717)
@@ -8,8 +8,8 @@
 classdef acenet
     properties (SetAccess=public) 
 		 % {{{
-		 name='glacdyn.ace-net.ca'
-		 %name='placentia.ace-net.ca'
+		 %name='glacdyn.ace-net.ca'
+		 name='placentia.ace-net.ca'
 		 %name='brasdor.ace-net.ca'
 		 login='klemorza';
 		 np=10;
@@ -18,7 +18,9 @@
 		 time=10;
 		 % codepath='/usr/local/issm-r11321/bin'; % this one is for issm on acenet global
 		 codepath='/home/klemorza/issm/trunk-jpl/bin'; % this one is for issm on my acenet directory
-		 executionpath='/home/klemorza/issm/trunk-jpl/execution';
+		 %executionpath='/home/klemorza/issm/trunk-jpl/execution';
+		 %executionpath='/home/klemorza/scratch/issmres.dir';
+		 executionpath='/net/glacdyn-data/glacdyn/1/klemorza/issm.dir';
 		 %}}}
 	 end
 	 methods
@@ -49,7 +51,7 @@
 		 function md = checkconsistency(cluster,md,solution,analyses) % {{{
 
 			 available_queues={'debug','shortq','longq'};
-			 queue_requirements_time=[60*1 60*3 60*17];
+			 queue_requirements_time=[48*1 48*2 48*5];
 			 queue_requirements_np=[32 128 256];
 
 			 QueueRequirements(available_queues,queue_requirements_time,queue_requirements_np,cluster.queue,cluster.np,cluster.time)
@@ -65,17 +67,32 @@
 			 fprintf(fid,'#!/bin/bash\n');
 			 fprintf(fid,'#$ -cwd\n');
           fprintf(fid,'#$ -N issm\n');
-          fprintf(fid,'#$ -l h_rt=96:0:0\n');
+          %fprintf(fid,'#$ -l h_rt=25:00:0\n');
+          %fprintf(fid,'#$ -l h_rt=47:59:00\n');
+          %fprintf(fid,'#$ -l h_rt=72:00:0\n');
+          fprintf(fid,'#$ -l h_rt=96:00:0\n');
+	  fprintf(fid,'#$ -l h_vmem=4G\n');
           fprintf(fid,'#$ -pe ompi* %i\n',cluster.np);
           fprintf(fid,'#$ -j y\n');
+          fprintf(fid,'#$ -l h=cl27*|cl28*|cl29*|cl30*|cl31*|cl320|cl267|cl268|cl269|cl338 \n');
+          %fprintf(fid,'#$ -l h=cl338 \n');
+          %fprintf(fid,'#$ -pe openmp 20 \n');
+          %fprintf(fid,'#$ -q !tarasov.q\n'); %
+          fprintf(fid,'#$ -pe openmp 8\n');
           fprintf(fid,'module purge\n');
-          fprintf(fid,'module load gcc openmpi/gcc\n');
-          fprintf(fid,'module load issm\n');
-			 fprintf(fid,'export ISSM_DIR="%s/../"\n',cluster.codepath); %FIXME
-			 fprintf(fid,'source $ISSM_DIR/etc/environment.sh\n');       %FIXME
+          %fprintf(fid,'module load gcc openmpi/gcc\n');
+          %fprintf(fid,'module unload openmpi\n');
+          fprintf(fid,'module load intel/12.1.7.367\n');
+          fprintf(fid,'module load openmpi/intel/1.2.9\n');
+
+          fprintf(fid,'module load gsl\n');
+          %fprintf(fid,'module load issm\n');
+          fprintf(fid,'export ISSM_DIR="%s/../"\n',cluster.codepath); %FIXME
+          fprintf(fid,'source $ISSM_DIR/etc/environment.sh\n');       %FIXME
           fprintf(fid,'\n');
-          fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog ',...
-                   cluster.np,cluster.codepath,EnumToString(solution),[cluster.executionpath '/' dirname],modelname,modelname,modelname);
+          fprintf(fid,'mpiexec %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog\n',...
+                   cluster.codepath,EnumToString(solution),[cluster.executionpath '/' dirname],modelname,modelname,modelname);
+          %fprintf(fid,'echo $HOSTNAME >>%s.outlog',modelname);
 			 fclose(fid);
 
 		 end
@@ -94,7 +111,7 @@
 
 			 disp('launching solution sequence on remote cluster');
 			 launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' dirname ' && mkdir ' dirname ...
-				 ' && cd ' dirname ' && mv ../' dirname '.tar.gz ./ && tar -zxf ' dirname '.tar.gz  && qsub ' modelname '.queue '];
+				 ' && cd ' dirname ' && mv ../' dirname '.tar.gz ./ && tar -zxf ' dirname '.tar.gz && qsub ' modelname '.queue '];
 			 issmssh(cluster.name,cluster.login,cluster.port,launchcommand);
 		 end %}}}
 		 function Download(cluster,dirname,filelist)% {{{
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsTdiffEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsTdiffEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsTdiffEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsTdiffEnum()
+%SURFACEFORCINGSTDIFFENUM - Enum of SurfaceforcingsTdiff
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsTdiffEnum()
+
+macro=StringToEnum('SurfaceforcingsTdiff');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsRlapslgmEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsRlapslgmEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsRlapslgmEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsRlapslgmEnum()
+%SURFACEFORCINGSRLAPSLGMENUM - Enum of SurfaceforcingsRlapslgm
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsRlapslgmEnum()
+
+macro=StringToEnum('SurfaceforcingsRlapslgm');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsRlapsEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsRlapsEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsRlapsEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsRlapsEnum()
+%SURFACEFORCINGSRLAPSENUM - Enum of SurfaceforcingsRlaps
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsRlapsEnum()
+
+macro=StringToEnum('SurfaceforcingsRlaps');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsS0tEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsS0tEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsS0tEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsS0tEnum()
+%SURFACEFORCINGSS0TENUM - Enum of SurfaceforcingsS0t
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsS0tEnum()
+
+macro=StringToEnum('SurfaceforcingsS0t');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsPfacEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsPfacEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsPfacEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsPfacEnum()
+%SURFACEFORCINGSPFACENUM - Enum of SurfaceforcingsPfac
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsPfacEnum()
+
+macro=StringToEnum('SurfaceforcingsPfac');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsSealevEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsSealevEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsSealevEnum.m	(revision 18717)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsSealevEnum()
+%SURFACEFORCINGSSEALEVENUM - Enum of SurfaceforcingsSealev
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsSealevEnum()
+
+macro=StringToEnum('SurfaceforcingsSealev');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18716)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18717)
@@ -328,6 +328,12 @@
 def SurfaceforcingsPrecipitationEnum(): return StringToEnum("SurfaceforcingsPrecipitation")[0]
 def SurfaceforcingsDesfacEnum(): return StringToEnum("SurfaceforcingsDesfac")[0]
 def SurfaceforcingsS0pEnum(): return StringToEnum("SurfaceforcingsS0p")[0]
+def SurfaceforcingsS0tEnum(): return StringToEnum("SurfaceforcingsS0t")[0]
+def SurfaceforcingsRlapsEnum(): return StringToEnum("SurfaceforcingsRlaps")[0]
+def SurfaceforcingsRlapslgmEnum(): return StringToEnum("SurfaceforcingsRlapslgm")[0]
+def SurfaceforcingsPfacEnum(): return StringToEnum("SurfaceforcingsPfac")[0]
+def SurfaceforcingsTdiffEnum(): return StringToEnum("SurfaceforcingsTdiff")[0]
+def SurfaceforcingsSealevEnum(): return StringToEnum("SurfaceforcingsSealev")[0]
 def SMBgradientsEnum(): return StringToEnum("SMBgradients")[0]
 def SurfaceforcingsMonthlytemperaturesEnum(): return StringToEnum("SurfaceforcingsMonthlytemperatures")[0]
 def SurfaceforcingsHrefEnum(): return StringToEnum("SurfaceforcingsHref")[0]
