source: issm/oecreview/Archive/18296-19100/ISSM-18967-18968.diff@ 19102

../trunk-jpl/test/NightlyRun/test236.py

@@ -15 +15 @@
 # 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 +50 @@
 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']

../trunk-jpl/test/NightlyRun/test237.py

@@ -11 +11 @@
 
 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 +52 @@
 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.)

../trunk-jpl/test/NightlyRun/test236.m

@@ -2 +2 @@
 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 +31 @@
     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 +43 @@
 % 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

../trunk-jpl/test/NightlyRun/test237.m

@@ -1 +1 @@
 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 +45 @@
 % 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 +66 @@
 
 %
 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', ...

../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -119 +119 @@
 
         int    stabilization,finiteelement,smb_model;
         bool   dakota_analysis;
-        bool   isdelta18o;
+        bool   isdelta18o,ismungsm;
         bool   isgroundingline;
         bool   islevelset;
 
@@ -178 +178 @@
                         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);

../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -343 +343 @@
         SurfaceforcingsDelta18oEnum,
         SurfaceforcingsDelta18oSurfaceEnum,
         SurfaceforcingsIsdelta18oEnum,
+        SurfaceforcingsIsmungsmEnum,
         SurfaceforcingsPrecipitationsPresentdayEnum,
         SurfaceforcingsPrecipitationsLgmEnum,
         SurfaceforcingsTemperaturesPresentdayEnum,

../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -349 +349 @@
                 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";

../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -355 +355 @@
               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 +382 @@
               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 +505 @@
               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 +628 @@
               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 +751 @@
               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;

../trunk-jpl/src/c/shared/Elements/elements.h

@@ -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* 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);
 

../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp

@@ -1 +1 @@
 /* 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 +23 @@
   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 +44 @@
   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 +67 @@
   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 +85 @@
     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 +101 @@
       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 +117 @@
         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 
   //******************************************************************
 

../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -21 +21 @@
         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 +107 @@
                         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);

../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

@@ -11 +11 @@
 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);

../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

@@ -10 +10 @@
 
         /*Intermediaties*/
         int  smb_model;
-        bool isdelta18o;
+        bool isdelta18o,ismungsm;
 
         /*First, get SMB model from parameters*/
         femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum);
@@ -22 +22 @@
                         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 +122 @@
         }
 
 }/*}}}*/
+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 +156 @@
         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 +165 @@
         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 +222 @@
 
         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);

../trunk-jpl/src/c/Makefile.am

@@ -219 +219 @@
                                         ./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\

../trunk-jpl/src/c/classes/Elements/Element.h

@@ -179 +179 @@
                 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 +251 @@
                 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;

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -608 +608 @@
                         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 +649 @@
         /*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 +2503 @@
         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;
 }

../trunk-jpl/src/c/classes/Elements/Tria.h

@@ -62 +62 @@
                 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 +109 @@
                 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);

../trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -607 +607 @@
                         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 +651 @@
         /*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 +2146 @@
 
 }
 /*}}}*/
-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 +2181 @@
        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){/*{{{*/

../trunk-jpl/src/c/classes/Elements/Penta.h

@@ -58 +58 @@
                 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 +137 @@
                 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);

../trunk-jpl/src/c/classes/Elements/Seg.h

@@ -53 +53 @@
                 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 +129 @@
                 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");};

../trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -53 +53 @@
                 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 +135 @@
                 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);

../trunk-jpl/src/m/classes/SMBpdd.py

@@ -25 +25 @@
                 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 +39 @@
         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 +107 @@
         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 +140 @@
 
                 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)
         # }}}

../trunk-jpl/src/m/classes/SMBpdd.m

@@ -16 +16 @@
                 Tdiff                     = NaN;
                 sealev                    = NaN;
                 isdelta18o                = 0;
+                ismungsm                  = 0;
                 delta18o                  = NaN;
                 delta18o_surface          = NaN;
                 temperatures_presentday   = NaN;
@@ -33 +34 @@
                         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 +73 @@
                                 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 +135 @@
 
                         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

../trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -341 +341 @@
 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]