Changeset 19172

Timestamp:

03/04/15 11:09:39 (10 years ago)

Author:

lemorzad

Message:

new method to interpolate present day temp and prec from delta O18

Location:

issm/trunk-jpl

Files:

• m4/analyses.m4 (modified) (32 diffs)
• src/c/Makefile.am (modified) (1 diff)
• src/c/analyses/MasstransportAnalysis.cpp (modified) (3 diffs)
• src/c/analyses/ThermalAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/classes/Elements/Penta.cpp (modified) (1 diff)
• src/c/classes/Elements/Penta.h (modified) (1 diff)
• src/c/classes/Elements/Seg.h (modified) (1 diff)
• src/c/classes/Elements/Tetra.h (modified) (1 diff)
• src/c/classes/Elements/Tria.cpp (modified) (1 diff)
• src/c/classes/Elements/Tria.h (modified) (1 diff)
• src/c/classes/Materials/Matpar.cpp (modified) (5 diffs)
• src/c/classes/Materials/Matpar.h (modified) (1 diff)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (5 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (4 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (modified) (1 diff)
• src/c/shared/Elements/PddSurfaceMassBalance.cpp (modified) (1 diff)
• src/c/shared/Elements/elements.h (modified) (1 diff)
• src/c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (7 diffs)
• src/m/classes/SMBpdd.m (modified) (14 diffs)
• src/m/classes/clusters/acenet.m (modified) (1 diff)
• src/m/enum/EnumDefinitions.py (modified) (2 diffs)

issm/trunk-jpl/m4/analyses.m4

@@ -19 +19 @@
 if test "x$ADJOINTBALANCETHICKNESS" = "xyes"; then
         HAVE_ADJOINTBALANCETHICKNESS=yes
-        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethicknesscapability])
+        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethickness capability])
 fi
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes])
@@ -33 +33 @@
 if test "x$ADJOINTBALANCETHICKNESS2" = "xyes"; then
         HAVE_ADJOINTBALANCETHICKNESS2=yes
-        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2capability])
+        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2 capability])
 fi
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes])
@@ -47 +47 @@
 if test "x$ADJOINTHORIZ" = "xyes"; then
         HAVE_ADJOINTHORIZ=yes
-        AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHorizcapability])
+        AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHoriz capability])
 fi
 AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes])
@@ -61 +61 @@
 if test "x$BALANCETHICKNESS" = "xyes"; then
         HAVE_BALANCETHICKNESS=yes
-        AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethicknesscapability])
+        AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethickness capability])
 fi
 AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes])
@@ -75 +75 @@
 if test "x$BALANCETHICKNESS2" = "xyes"; then
         HAVE_BALANCETHICKNESS2=yes
-        AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2capability])
+        AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2 capability])
 fi
 AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes])
@@ -89 +89 @@
 if test "x$BALANCETHICKNESSSOFT" = "xyes"; then
         HAVE_BALANCETHICKNESSSOFT=yes
-        AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoftcapability])
+        AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoft capability])
 fi
 AM_CONDITIONAL([BALANCETHICKNESSSOFT], [test x$HAVE_BALANCETHICKNESSSOFT = xyes])
@@ -103 +103 @@
 if test "x$BALANCEVELOCITY" = "xyes"; then
         HAVE_BALANCEVELOCITY=yes
-        AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocitycapability])
+        AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocity capability])
 fi
 AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes])
@@ -117 +117 @@
 if test "x$L2PROJECTIONEPL" = "xyes"; then
         HAVE_L2PROJECTIONEPL=yes
-        AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPLcapability])
+        AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPL capability])
 fi
 AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes])
@@ -131 +131 @@
 if test "x$L2PROJECTIONBASE" = "xyes"; then
         HAVE_L2PROJECTIONBASE=yes
-        AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBasecapability])
+        AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBase capability])
 fi
 AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes])
@@ -145 +145 @@
 if test "x$DAMAGEEVOLUTION" = "xyes"; then
         HAVE_DAMAGEEVOLUTION=yes
-        AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolutioncapability])
+        AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolution capability])
 fi
 AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes])
@@ -159 +159 @@
 if test "x$STRESSBALANCE" = "xyes"; then
         HAVE_STRESSBALANCE=yes
-        AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalancecapability])
+        AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalance capability])
 fi
 AM_CONDITIONAL([STRESSBALANCE], [test x$HAVE_STRESSBALANCE = xyes])
@@ -173 +173 @@
 if test "x$STRESSBALANCESIA" = "xyes"; then
         HAVE_STRESSBALANCESIA=yes
-        AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIAcapability])
+        AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIA capability])
 fi
 AM_CONDITIONAL([STRESSBALANCESIA], [test x$HAVE_STRESSBALANCESIA = xyes])
@@ -187 +187 @@
 if test "x$STRESSBALANCEVERTICAL" = "xyes"; then
         HAVE_STRESSBALANCEVERTICAL=yes
-        AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVerticalcapability])
+        AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVertical capability])
 fi
 AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes])
@@ -201 +201 @@
 if test "x$ENTHALPY" = "xyes"; then
         HAVE_ENTHALPY=yes
-        AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpycapability])
+        AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpy capability])
 fi
 AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes])
@@ -215 +215 @@
 if test "x$HYDROLOGYSHREVE" = "xyes"; then
         HAVE_HYDROLOGYSHREVE=yes
-        AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShrevecapability])
+        AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShreve capability])
 fi
 AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes])
@@ -229 +229 @@
 if test "x$HYDROLOGYDCINEFFICIENT" = "xyes"; then
         HAVE_HYDROLOGYDCINEFFICIENT=yes
-        AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficientcapability])
+        AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficient capability])
 fi
 AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes])
@@ -243 +243 @@
 if test "x$HYDROLOGYDCEFFICIENT" = "xyes"; then
         HAVE_HYDROLOGYDCEFFICIENT=yes
-        AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficientcapability])
+        AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficient capability])
 fi
 AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes])
@@ -257 +257 @@
 if test "x$MELTING" = "xyes"; then
         HAVE_MELTING=yes
-        AC_DEFINE([_HAVE_MELTING_],[1],[with Meltingcapability])
+        AC_DEFINE([_HAVE_MELTING_],[1],[with Melting capability])
 fi
 AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes])
@@ -271 +271 @@
 if test "x$MASSTRANSPORT" = "xyes"; then
         HAVE_MASSTRANSPORT=yes
-        AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransportcapability])
+        AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransport capability])
 fi
 AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
@@ -285 +285 @@
 if test "x$FREESURFACEBASE" = "xyes"; then
         HAVE_FREESURFACEBASE=yes
-        AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBasecapability])
+        AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBase capability])
 fi
 AM_CONDITIONAL([FREESURFACEBASE], [test x$HAVE_FREESURFACEBASE = xyes])
@@ -299 +299 @@
 if test "x$FREESURFACETOP" = "xyes"; then
         HAVE_FREESURFACETOP=yes
-        AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTopcapability])
+        AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTop capability])
 fi
 AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes])
@@ -313 +313 @@
 if test "x$EXTRUDEFROMBASE" = "xyes"; then
         HAVE_EXTRUDEFROMBASE=yes
-        AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBasecapability])
+        AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBase capability])
 fi
 AM_CONDITIONAL([EXTRUDEFROMBASE], [test x$HAVE_EXTRUDEFROMBASE = xyes])
@@ -327 +327 @@
 if test "x$EXTRUDEFROMTOP" = "xyes"; then
         HAVE_EXTRUDEFROMTOP=yes
-        AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTopcapability])
+        AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTop capability])
 fi
 AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes])
@@ -341 +341 @@
 if test "x$DEPTHAVERAGE" = "xyes"; then
         HAVE_DEPTHAVERAGE=yes
-        AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAveragecapability])
+        AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAverage capability])
 fi
 AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes])
@@ -355 +355 @@
 if test "x$SMOOTH" = "xyes"; then
         HAVE_SMOOTH=yes
-        AC_DEFINE([_HAVE_SMOOTH_],[1],[with Smoothcapability])
+        AC_DEFINE([_HAVE_SMOOTH_],[1],[with Smooth capability])
 fi
 AM_CONDITIONAL([SMOOTH], [test x$HAVE_SMOOTH = xyes])
@@ -369 +369 @@
 if test "x$THERMAL" = "xyes"; then
         HAVE_THERMAL=yes
-        AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermalcapability])
+        AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermal capability])
 fi
 AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes])
@@ -383 +383 @@
 if test "x$UZAWAPRESSURE" = "xyes"; then
         HAVE_UZAWAPRESSURE=yes
-        AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressurecapability])
+        AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressure capability])
 fi
 AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes])
@@ -397 +397 @@
 if test "x$GIA" = "xyes"; then
         HAVE_GIA=yes
-        AC_DEFINE([_HAVE_GIA_],[1],[with Giacapability])
+        AC_DEFINE([_HAVE_GIA_],[1],[with Gia capability])
 fi
 AM_CONDITIONAL([GIA], [test x$HAVE_GIA = xyes])
 AC_MSG_RESULT($HAVE_GIA)
-dnl }}}
-dnl with-Seaice{{{
-AC_ARG_WITH([Seaice],
-        AS_HELP_STRING([--with-Seaice = YES], [compile with Seaice capabilities (default is yes)]),
-        [SEAICE=$withval],[SEAICE=yes])
-AC_MSG_CHECKING(for Seaice capability compilation)
-
-HAVE_SEAICE=no 
-if test "x$SEAICE" = "xyes"; then
-        HAVE_SEAICE=yes
-        AC_DEFINE([_HAVE_SEAICE_],[1],[with Seaicecapability])
-fi
-AM_CONDITIONAL([SEAICE], [test x$HAVE_SEAICE = xyes])
-AC_MSG_RESULT($HAVE_SEAICE)
 dnl }}}
 dnl with-Meshdeformation{{{
@@ -425 +411 @@
 if test "x$MESHDEFORMATION" = "xyes"; then
         HAVE_MESHDEFORMATION=yes
-        AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformationcapability])
+        AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformation capability])
 fi
 AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes])
@@ -439 +425 @@
 if test "x$LEVELSET" = "xyes"; then
         HAVE_LEVELSET=yes
-        AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelsetcapability])
+        AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelset capability])
 fi
 AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes])
@@ -453 +439 @@
 if test "x$EXTRAPOLATION" = "xyes"; then
         HAVE_EXTRAPOLATION=yes
-        AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolationcapability])
+        AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolation capability])
 fi
 AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes])
@@ -467 +453 @@
 if test "x$LSFREINITIALIZATION" = "xyes"; then
         HAVE_LSFREINITIALIZATION=yes
-        AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitializationcapability])
+        AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitialization capability])
 fi
 AM_CONDITIONAL([LSFREINITIALIZATION], [test x$HAVE_LSFREINITIALIZATION = xyes])

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

@@ -114 +114 @@
                                         ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\
                                         ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\
+                                        ./shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp\
                                         ./shared/Elements/DrainageFunctionWaterfraction.cpp\
                                         ./shared/String/DescriptorIndex.cpp\

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

@@ -120 +120 @@
         int    stabilization,finiteelement,smb_model;
         bool   dakota_analysis;
-        bool   isdelta18o,ismungsm;
+        bool   isdelta18o,ismungsm,isd18opd;
         bool   isgroundingline;
         bool   islevelset;
@@ -180 +180 @@
                         iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
                         iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
+                        iomodel->Constant(&isd18opd,SurfaceforcingsIsd18opdEnum);
                         iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
                         if(isdelta18o || ismungsm){
@@ -187 +188 @@
                                 iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum);
                         }
+                        else if (isd18opd){
+                                iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum);
+                                iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum);
+                        }
                         else{
                                 iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
                                 iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum);
                         }
-
                         break;
                 case SMBgradientsEnum:

issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

@@ -35 +35 @@
         }
         else{
-                IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement);
+                IoModelToDynamicConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement);
         }
 

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

@@ -183 +183 @@
                 virtual void       Delta18oParameterization(void)=0;
                 virtual void       MungsmtpParameterization(void)=0;
+                virtual void       Delta18opdParameterization(void)=0;
                 virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
                 virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;

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

@@ -713 +713 @@
         this->inputs->AddInput(NewPrecipitationInput);
         
+        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
+        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+void       Penta::Delta18opdParameterization(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];
+        IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
+        IssmDouble tmp[NUMVERTICES];
+        IssmDouble Delta18oTime;
+        IssmDouble dpermil;
+        IssmDouble time,yts;
+        this->parameters->FindParam(&time,TimeEnum);
+        this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+        /*Get some pdd parameters*/
+        dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum);
+
+        /*Recover present day temperature and precipitation*/
+        Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
+        Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
+        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(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+                }
+        }
+
+        /*Recover interpolation parameters at time t*/
+        this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
+
+        /*Compute the temperature and precipitation*/
+        for(int iv=0;iv<NUMVERTICES;iv++){
+          ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
+                                        &PrecipitationsPresentday[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);

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

@@ -60 +60 @@
                 void           Delta18oParameterization(void);
                 void           MungsmtpParameterization(void);
+                void           Delta18opdParameterization(void);
                 void           ElementResponse(IssmDouble* presponse,int response_enum);
                 void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);

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

@@ -55 +55 @@
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
                 void        MungsmtpParameterization(void){_error_("not implemented yet");};
+                void        Delta18opdParameterization(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");};

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

@@ -55 +55 @@
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
                 void        MungsmtpParameterization(void){_error_("not implemented yet");};
+                void        Delta18opdParameterization(void){_error_("not implemented yet");};
                 IssmDouble  DragCoefficientAbsGradient(void){_error_("not implemented yet");};
                 void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};

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

@@ -750 +750 @@
                                         &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;
+}
+/*}}}*/
+void       Tria::Delta18opdParameterization(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];
+        IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
+        IssmDouble tmp[NUMVERTICES];
+        IssmDouble Delta18oTime;
+        IssmDouble dpermil;
+        IssmDouble time,yts;
+        this->parameters->FindParam(&time,TimeEnum);
+        this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+        /*Get some pdd parameters*/
+        dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum);
+                        
+        /*Recover present day temperature and precipitation*/
+        Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
+        Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
+        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(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+                }
+        }
+
+        /*Recover interpolation parameters at time t*/
+        this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
+
+        /*Compute the temperature and precipitation*/
+        for(int iv=0;iv<NUMVERTICES;iv++){
+          ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
+                                        &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0], 
                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
         }

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

@@ -65 +65 @@
                 void        Delta18oParameterization(void);
                 void        MungsmtpParameterization(void);
+                void        Delta18opdParameterization(void);
                 int         EdgeOnBaseIndex();
                 void        EdgeOnBaseIndices(int* pindex1,int* pindex);

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -55 +55 @@
                                         iomodel->Constant(&this->rlaps,SurfaceforcingsRlapsEnum);
                                         iomodel->Constant(&this->rlapslgm,SurfaceforcingsRlapslgmEnum);
+                                        iomodel->Constant(&this->dpermil,SurfaceforcingsDpermilEnum);
                                         break;
                                 case SMBgradientsEnum:
@@ -134 +135 @@
         _printf_("   rlaps: " << rlaps << "\n");
         _printf_("   rlapslgm: " << rlapslgm << "\n");
+        _printf_("   dpermil: " << dpermil << "\n");
         return;
 }
@@ -179 +181 @@
         matpar->rlaps=this->rlaps;
         matpar->rlapslgm=this->rlapslgm;
+        matpar->dpermil=this->dpermil;
 
         matpar->sediment_compressibility=this->sediment_compressibility;
@@ -273 +276 @@
                 case SurfaceforcingsRlapslgmEnum:
                         this->rlapslgm=constant;
+                        break;
+                case SurfaceforcingsDpermilEnum:
+                        this->dpermil=constant;
                         break;
                 default: 
@@ -400 +406 @@
                 case SurfaceforcingsRlapsEnum:               return this->rlaps;
                 case SurfaceforcingsRlapslgmEnum:            return this->rlapslgm;
+                case SurfaceforcingsDpermilEnum:             return this->dpermil;
                 case MaterialsLithosphereShearModulusEnum:   return this->lithosphere_shear_modulus;
                 case MaterialsLithosphereDensityEnum:        return this->lithosphere_density;

issm/trunk-jpl/src/c/classes/Materials/Matpar.h

@@ -36 +36 @@
                 IssmDouble  rlaps;
                 IssmDouble  rlapslgm;
+                IssmDouble  dpermil;
 
                 /*hydrology Dual Porous Continuum: */    

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

@@ -22 +22 @@
         char**      requestedoutputs = NULL;
         IssmDouble  time;
-        bool        isdelta18o,ismungsm;
+        bool        isdelta18o,ismungsm,isd18opd;
 
         /*parameters for mass flux:*/
@@ -107 +107 @@
                         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
                         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsd18opdEnum));
                         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum));
                         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum));
@@ -114 +115 @@
                         iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
                         iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
+                        iomodel->Constant(&isd18opd,SurfaceforcingsIsd18opdEnum);
 
                         if(ismungsm){
@@ -131 +133 @@
                           iomodel->DeleteData(temp,SurfaceforcingsSealevEnum);
                         }
-
                         if(isdelta18o){
                                 iomodel->Constant(&yts,ConstantsYtsEnum);
@@ -144 +145 @@
                                 parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oSurfaceEnum,&temp[0],&temp[M],M));
                                 iomodel->DeleteData(temp,SurfaceforcingsDelta18oSurfaceEnum);
+                        }
+                        if(isd18opd){
+                                iomodel->Constant(&yts,ConstantsYtsEnum);
+
+                                iomodel->FetchData(&temp,&N,&M,SurfaceforcingsDelta18oEnum); _assert_(N==2);
+                                for(i=0;i<M;i++) temp[M+i]=yts*temp[M+i];
+                                parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oEnum,&temp[0],&temp[M],M));
+                                iomodel->DeleteData(temp,SurfaceforcingsDelta18oEnum);
+                                
+                                parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDpermilEnum));
                         }
                         break;

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

@@ -11 +11 @@
         /*Intermediaties*/
         int  smb_model;
-        bool isdelta18o,ismungsm;
+        bool isdelta18o,ismungsm,isd18opd;
 
         /*First, get SMB model from parameters*/
@@ -24 +24 @@
                         femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
                         femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum);
+                        femmodel->parameters->FindParam(&isd18opd,SurfaceforcingsIsd18opdEnum);
                         if(isdelta18o){
                                 if(VerboseSolution()) _printf0_("   call Delta18oParameterization module\n");
@@ -31 +32 @@
                                 if(VerboseSolution()) _printf0_("   call MungsmtpParameterization module\n");
                                 MungsmtpParameterizationx(femmodel);
+                        } 
+                        if(isd18opd){
+                                if(VerboseSolution()) _printf0_("   call Delta18opdParameterization module\n");
+                                Delta18opdParameterizationx(femmodel);
                         } 
                         if(VerboseSolution()) _printf0_("   call positive degree day module\n");
@@ -128 +133 @@
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 element->MungsmtpParameterization();
+        }
+
+}/*}}}*/
+void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
+
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element->Delta18opdParameterization();
         }
 

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

@@ -13 +13 @@
 void Delta18oParameterizationx(FemModel* femmodel);
 void MungsmtpParameterizationx(FemModel* femmodel);
+void Delta18opdParameterizationx(FemModel* femmodel);
 void PositiveDegreeDayx(FemModel* femmodel);
 void SmbHenningx(FemModel* femmodel);

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

@@ -90 +90 @@
     Tsurf = tstar*deltm+Tsurf;        
 
-      /*********compute PD ****************/
-      if (tstar < PDup){
-        pd = 1.;
-        if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
-      else { 
-        pd = 0.;}
-
-      /******exp des/elev precip reduction*******/
-      sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
-      if (sp>0.0){q = exp(-desfac*sp);}
-      else {q = 1.0;}
-
-      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;
-
-      /*********compute PDD************/
-      // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
-      // gaussian=T_m, so ndd=-(Tsurf-pdd)
-      if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;} 
-
-      if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
-      else if (tstar> -siglim){
-        pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
-        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];      
+    /*********compute PD ****************/
+    if (tstar < PDup){
+      pd = 1.;
+      if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
+    else { 
+      pd = 0.;}
+    
+    /******exp des/elev precip reduction*******/
+    sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
+    if (sp>0.0){q = exp(-desfac*sp);}
+    else {q = 1.0;}
+    
+    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;
+
+    /*********compute PDD************/
+    // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
+    // gaussian=T_m, so ndd=-(Tsurf-pdd)
+    if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;} 
+    
+    if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
+    else if (tstar> -siglim){
+      pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
+      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 
   //******************************************************************
 
-    saccu = qm;
-    prect = qmp;     // total precipitation during 1 year taking into account des. ef.
-    Tsum=Tsum/3;
-
-    /***** determine PDD factors *****/
-    if(Tsum< -1.) {
-      snwmf=2.65*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
-      smf=17.22*0.001;    //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
-    } 
-    else if(Tsum< 10){
-      snwmf = (0.15*Tsum + 2.8)*0.001;
-      smf = (0.0067*pow((10.-Tsum),3) + 8.3)*0.001;
-    }
-    else{
-      snwmf=4.3*0.001;
-      smf=8.3*0.001;
-    }
-    snwmf=0.95*snwmf;
-    smf=0.95*smf;
-
-    /*****  compute PDD ablation and refreezing *****/
-    pddt = pdd *365;
-    snwm = snwmf*pddt;       // snow that could have been melted in a year
-    hmx2 = min(h,dfrz);   // refreeze active layer max depth: dfrz
-
-    if(snwm < saccu) {
-      water=prect-saccu + snwm; //water=rain + snowmelt
-      //     l 2.2= capillary factor
-      //     Should refreezing be controlled by frzndd or by mean annual Tsurf?
-      //     dCovrLm concept is of warming of active layer (thickness =d@=1-
-      //     >2m)
-      //     problem with water seepage into ice: should be sealed after
-      //     refreezing
-      //     so everything needs to be predicated on 1 year scale, except for
-      //     thermal
-      //     conductivity through ice
-      //     also, need to account that melt season has low accum, so what's
-      //     going to
-      //     hold the meltwater around for refreezing? And melt-time will have
-      //     low seasonal frzndd
-
-      //      Superimposed ice :  Pfeffer et al. 1991, Tarasov 2002
-
-      supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice
-      supcap=min(2.2*(saccu-snwm),water);
-      runoff=snwm - supice;  //meltwater only, does not include rain
-    }
-    else {  //all snow melted
-      supice= min(hmx2*CovrLm*frzndd, prect );
-      runoff= saccu + smf*(pddt-saccu/snwmf) - supice;
-      supcap=0;
-    }
-    //     pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
-    //     except pdd melt heat source is atmosphere, while refreeze is
-    //     ground/ice stored interim
-    //     assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
-    //     assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
-    //     <0
-    //     assume ndd>pdd, little melt => little supice 
-    //     bottom line: compare for Tsurf<0 : supice and no supice case,
-    //     expect Tsurf difference
-    //     except some of cooling flux comes from atmosphere//
-    //     1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
-    //     does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
-    //     < 0.1 
-
-    //     make more sense to just use residual pdd-ndd except that pdd
-    //     residual not clear yet
-    //     frzndd should not be used up by refreezing in snow, so stick in
-    //     supcap.
-    diffndd=0;
-    if (frzndd>0) {
-      diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd);
-      frzndd=frzndd-diffndd;
-    }
-    if(runoff<0){
-      saccu= saccu -runoff;
-      smelt = 0;
-      precrunoff=prect-saccu;
-      //here assume pdd residual is 0, => 
-      Tsurf= max(Tsurf,-frzndd);
-    }
-    else {
-      smelt = runoff;
-      precrunoff=prect-max(0.,supice)-saccu;}
-    //here really need pdd balance, try 0.5 fudge factor?
-    //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
-    //yet from site plots, can be ice free with Tsurf=-5.5C
-    if(Tsurf<0) {
-      Tsurf= min(Tsurf+fsupT*diffndd , 0.);}
-
-    B = -smelt+saccu;
-    B = B/yts;
-    pddtj=pddt;
+  saccu = qm;
+  prect = qmp;     // total precipitation during 1 year taking into account des. ef.
+  Tsum=Tsum/3;
+  
+  /***** determine PDD factors *****/
+  if(Tsum< -1.) {
+    snwmf=2.65*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
+    smf=17.22*0.001;    //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
+  } 
+  else if(Tsum< 10){
+    snwmf = (0.15*Tsum + 2.8)*0.001;
+    smf = (0.0067*pow((10.-Tsum),3) + 8.3)*0.001;
+  }
+  else{
+    snwmf=4.3*0.001;
+    smf=8.3*0.001;
+  }
+  snwmf=0.95*snwmf;
+  smf=0.95*smf;
+  
+  /*****  compute PDD ablation and refreezing *****/
+  pddt = pdd *365;
+  snwm = snwmf*pddt;       // snow that could have been melted in a year
+  hmx2 = min(h,dfrz);   // refreeze active layer max depth: dfrz
+  
+  if(snwm < saccu) {
+    water=prect-saccu + snwm; //water=rain + snowmelt
+    //     l 2.2= capillary factor
+    //     Should refreezing be controlled by frzndd or by mean annual Tsurf?
+    //     dCovrLm concept is of warming of active layer (thickness =d@=1-
+    //     >2m)
+    //     problem with water seepage into ice: should be sealed after
+    //     refreezing
+    //     so everything needs to be predicated on 1 year scale, except for
+    //     thermal
+    //     conductivity through ice
+    //     also, need to account that melt season has low accum, so what's
+    //     going to
+    //     hold the meltwater around for refreezing? And melt-time will have
+    //     low seasonal frzndd
+
+    //      Superimposed ice :  Pfeffer et al. 1991, Tarasov 2002
+
+    supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice
+    supcap=min(2.2*(saccu-snwm),water);
+    runoff=snwm - supice;  //meltwater only, does not include rain
+  }
+  else {  //all snow melted
+    supice= min(hmx2*CovrLm*frzndd, prect );
+    runoff= saccu + smf*(pddt-saccu/snwmf) - supice;
+    supcap=0;
+  }
+  //     pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
+  //     except pdd melt heat source is atmosphere, while refreeze is
+  //     ground/ice stored interim
+  //     assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
+  //     assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
+  //     <0
+  //     assume ndd>pdd, little melt => little supice 
+  //     bottom line: compare for Tsurf<0 : supice and no supice case,
+  //     expect Tsurf difference
+  //     except some of cooling flux comes from atmosphere//
+  //     1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
+  //     does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
+  //     < 0.1 
+
+  //     make more sense to just use residual pdd-ndd except that pdd
+  //     residual not clear yet
+  //     frzndd should not be used up by refreezing in snow, so stick in
+  //     supcap.
+  diffndd=0;
+  if (frzndd>0) {
+    diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd);
+    frzndd=frzndd-diffndd;
+  }
+  if(runoff<0){
+    saccu= saccu -runoff;
+    smelt = 0;
+    precrunoff=prect-saccu;
+    //here assume pdd residual is 0, => 
+    Tsurf= max(Tsurf,-frzndd);
+  }
+  else {
+    smelt = runoff;
+    precrunoff=prect-max(0.,supice)-saccu;}
+  //here really need pdd balance, try 0.5 fudge factor?
+  //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
+  //yet from site plots, can be ice free with Tsurf=-5.5C
+  if(Tsurf<0) {
+    Tsurf= min(Tsurf+fsupT*diffndd , 0.);}
+
+  B = -smelt+saccu;
+  B = B/yts;
+  pddtj=pddt;
 
   return B;

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

@@ -28 +28 @@
                                            IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
                                            IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
+void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,
+                                               IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
+                                               IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);  
 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
 IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);

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

@@ -348 +348 @@
         SurfaceforcingsIsdelta18oEnum,
         SurfaceforcingsIsmungsmEnum,
+        SurfaceforcingsIsd18opdEnum,
         SurfaceforcingsPrecipitationsPresentdayEnum,
         SurfaceforcingsPrecipitationsLgmEnum,
@@ -361 +362 @@
         SurfaceforcingsTdiffEnum,
         SurfaceforcingsSealevEnum,
+        SurfaceforcingsDpermilEnum,
         SMBgradientsEnum,
         SurfaceforcingsMonthlytemperaturesEnum,

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

@@ -354 +354 @@
                 case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o";
                 case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm";
+                case SurfaceforcingsIsd18opdEnum : return "SurfaceforcingsIsd18opd";
                 case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday";
                 case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm";
@@ -367 +368 @@
                 case SurfaceforcingsTdiffEnum : return "SurfaceforcingsTdiff";
                 case SurfaceforcingsSealevEnum : return "SurfaceforcingsSealev";
+                case SurfaceforcingsDpermilEnum : return "SurfaceforcingsDpermil";
                 case SMBgradientsEnum : return "SMBgradients";
                 case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures";

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

@@ -360 +360 @@
               else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum;
               else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum;
+              else if (strcmp(name,"SurfaceforcingsIsd18opd")==0) return SurfaceforcingsIsd18opdEnum;
               else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum;
               else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum;
@@ -373 +374 @@
               else if (strcmp(name,"SurfaceforcingsTdiff")==0) return SurfaceforcingsTdiffEnum;
               else if (strcmp(name,"SurfaceforcingsSealev")==0) return SurfaceforcingsSealevEnum;
+              else if (strcmp(name,"SurfaceforcingsDpermil")==0) return SurfaceforcingsDpermilEnum;
               else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
               else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
@@ -381 +383 @@
               else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
               else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
-              else if (strcmp(name,"SurfaceforcingsAccumulation")==0) return SurfaceforcingsAccumulationEnum;
-              else if (strcmp(name,"SurfaceforcingsEvaporation")==0) return SurfaceforcingsEvaporationEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum;
+              if (strcmp(name,"SurfaceforcingsAccumulation")==0) return SurfaceforcingsAccumulationEnum;
+              else if (strcmp(name,"SurfaceforcingsEvaporation")==0) return SurfaceforcingsEvaporationEnum;
+              else if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum;
               else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
               else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum;
@@ -504 +506 @@
               else if (strcmp(name,"Penpair")==0) return PenpairEnum;
               else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
-              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
-              else if (strcmp(name,"Masscon")==0) return MassconEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"MassconName")==0) return MassconNameEnum;
+              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 if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
@@ -627 +629 @@
               else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
               else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
-              else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
-              else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
+              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 if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
@@ -750 +752 @@
               else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
               else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
-              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
-              else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
+              if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
+              else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
+              else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
               else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
               else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
@@ -873 +875 @@
               else if (strcmp(name,"XY")==0) return XYEnum;
               else if (strcmp(name,"XYZ")==0) return XYZEnum;
-              else if (strcmp(name,"Dense")==0) return DenseEnum;
-              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
+              if (strcmp(name,"Dense")==0) return DenseEnum;
+              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
+              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
               else if (strcmp(name,"Seq")==0) return SeqEnum;
               else if (strcmp(name,"Mpi")==0) return MpiEnum;

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

@@ -12 +12 @@
                 s0t                       = 0;
                 rlaps                     = 0;
-                rlapslgm                  = 0;                
+                rlapslgm                  = 0; 
+                dpermil                   = 0; 
                 Pfac                      = NaN;
                 Tdiff                     = NaN;
@@ -18 +19 @@
                 isdelta18o                = 0;
                 ismungsm                  = 0;
+                isd18opd                  = 0;
                 delta18o                  = NaN;
                 delta18o_surface          = NaN;
@@ -35 +37 @@
                 end % }}}
                 function self = extrude(self,md) % {{{
-                        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==0 & self.ismungsm==0 & self.isd18opd==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end
+                        if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==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
@@ -45 +47 @@
                         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
+                        if(self.isd18opd),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
+                        if(self.isd18opd),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
 
                 end % }}}
@@ -59 +63 @@
                   self.isdelta18o = 0;
                   self.ismungsm   = 0;
+                  self.isd18opd   = 0;
                   self.desfac     = 0.5;
                   self.s0p        = 0;
@@ -64 +69 @@
                   self.rlaps      = 6.5;
                   self.rlapslgm   = 6.5;
+                  self.dpermil    = 2.4;
                   
                 end % }}}
@@ -74 +80 @@
                                 md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1);
                                 md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1);
-                                if(self.isdelta18o==0 & self.ismungsm==0)
+                                if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0)
                                         md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','timeseries',1,'NaN',1);
                                         md = checkfield(md,'fieldname','surfaceforcings.precipitation','timeseries',1,'NaN',1);
@@ -94 +100 @@
                                         md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
                                         md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
+                                elseif(self.isd18opd==1) 
+                                        md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','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.delta18o','NaN',1);
+                                        md = checkfield(md,'fieldname','surfaceforcings.dpermil','>=',0,'numel',1);
                                 end
                         end
@@ -103 +114 @@
                         fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)');
                         fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)');
+                        fielddisplay(self,'isd18opd','is delta18o parametrisation from present day temperature and precipitation activated (0 or 1, default is 0)');
                         fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]');
                         fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]');
@@ -108 +120 @@
                         fielddisplay(self,'rlaps','present day lapse rate [degree/km]');
                         fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]');
-                        if(self.isdelta18o==0 & self.ismungsm==0)
+                        if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0)
                             fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
                             fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']);
@@ -116 +128 @@
                             fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
                             fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
-                            fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
-                            fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated');
+                            fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm is activated');
                             fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
                             fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated');
                         elseif(self.ismungsm==1)
-                            fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
+                            fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated');
                             fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
-                            fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
-                            fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
+                            fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated');
+                            fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm is activated');
                             fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated');
                             fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
                             fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated');
+                        elseif(self.isd18opd==1) 
+                            fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated');
+                            fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated');
+                            fielddisplay(self,'delta18o','delta18o, required if pdd is activated and d18opd activated');  
+                            fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated');                            
                         end
                 end % }}}
@@ -138 +155 @@
                         WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean');
                         WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean');
+                        WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isd18opd','format','Boolean');
                         WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double');
                         WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double');
@@ -144 +162 @@
                         WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double');
 
-                        if(self.isdelta18o==0 & self.ismungsm==0)
+                        if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0)
                                 %WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
                                 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
@@ -165 +183 @@
                                 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
                                 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
+                        elseif self.isd18opd
+                                WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1);
+                                WriteData(fid,'object',self,'class','surfaceforcings','fieldname','dpermil','format','Double');
                         end
                 end % }}}

issm/trunk-jpl/src/m/classes/clusters/acenet.m

@@ -68 +68 @@
                          fprintf(fid,'#$ -cwd\n');
           fprintf(fid,'#$ -N issm\n');
-          %fprintf(fid,'#$ -l h_rt=25:00:0\n');
-          %fprintf(fid,'#$ -l h_rt=47:59:00\n');
-          %fprintf(fid,'#$ -l h_rt=72:00:0\n');
-          fprintf(fid,'#$ -l h_rt=96:00:0\n');
-          fprintf(fid,'#$ -l h_vmem=4G\n');
-          fprintf(fid,'#$ -pe ompi* %i\n',cluster.np);
-          fprintf(fid,'#$ -j y\n');
-          fprintf(fid,'#$ -l h=cl27*|cl28*|cl29*|cl30*|cl31*|cl320|cl267|cl268|cl269|cl338 \n');
+          % fprintf(fid,'#$ -l h_rt=25:00:0\n');
+          fprintf(fid,'#$ -l h_rt=47:59:00\n');
+          % fprintf(fid,'#$ -l h_rt=72:00:0\n');
+          % fprintf(fid,'#$ -l h_rt=96:00:0\n');
+          % fprintf(fid,'#$ -l h_rt=336:00:0\n');
+
+          fprintf(fid,'#$ -l h_vmem=3G\n');
+          % if cluster.np>10
+          %     fprintf(fid,'#$ -l h_vmem=3G\n');
+          % else
+          %     fprintf(fid,'#$ -l h_vmem=3G\n');
+          % end
+
+          % ---- Which queue to use ----
+          %fprintf(fid,'#$ -q !tarasov.q\n'); %
+          %fprintf(fid,'#$ -q medium.q@*,short.q@*\n');
+          fprintf(fid,'#$ -q short.q@*\n');
+
+          % ---- Which node are selected ----
+          % fprintf(fid,'#$ -l h=cl27*|cl28*|cl29*|cl30*|cl31*|cl320|cl267|cl268|cl269|cl338 \n');
+          %fprintf(fid,'#$ -l h=cl0* \n');
           %fprintf(fid,'#$ -l h=cl338 \n');
-          %fprintf(fid,'#$ -pe openmp 20 \n');
-          %fprintf(fid,'#$ -q !tarasov.q\n'); %
-          fprintf(fid,'#$ -pe openmp 8\n');
+          % Acenet nodes with 16cpus and more than 60G mem
+          % fprintf(fid,'#$ -l h=cl001|cl002|cl003|cl004|cl005|cl006|cl007|cl008|cl009|cl010|cl011|cl012|cl021|cl022|cl023|cl024 \n');
+  
+          % ---- cpus on different nodes ----
+          fprintf(fid,'#$ -pe ompi %i\n',cluster.np); % To avoid green acenet that does not have InfiniBand
+          % fprintf(fid,'#$ -pe ompi* %i\n',cluster.np);
+          % -------- All cpus in the same node --------          
+          % fprintf(fid,'#$ -pe openmp %i\n',cluster.np);
+          
+          % ---- misc ----
+          fprintf(fid,'#$ -j y\n');        
+                  
           fprintf(fid,'module purge\n');
           %fprintf(fid,'module load gcc openmpi/gcc\n');

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

@@ -346 +346 @@
 def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0]
 def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0]
+def SurfaceforcingsIsd18opdEnum(): return StringToEnum("SurfaceforcingsIsd18opd")[0]
 def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0]
 def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0]
@@ -359 +360 @@
 def SurfaceforcingsTdiffEnum(): return StringToEnum("SurfaceforcingsTdiff")[0]
 def SurfaceforcingsSealevEnum(): return StringToEnum("SurfaceforcingsSealev")[0]
+def SurfaceforcingsDpermilEnum(): return StringToEnum("SurfaceforcingsDpermil")[0]
 def SMBgradientsEnum(): return StringToEnum("SMBgradients")[0]
 def SurfaceforcingsMonthlytemperaturesEnum(): return StringToEnum("SurfaceforcingsMonthlytemperatures")[0]