Changeset 22448

Timestamp:

02/21/18 18:19:30 (7 years ago)

Author:

schlegel

Message:

CHG: add pdd factor input field

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SmbAnalysis.cpp (modified) (7 diffs)
• c/classes/Elements/Element.cpp (modified) (4 diffs)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (3 diffs)
• c/shared/Elements/PddSurfaceMassBalance.cpp (modified) (4 diffs)
• c/shared/Elements/elements.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (8 diffs)
• m/classes/SMBd18opdd.m (modified) (8 diffs)
• m/classes/SMBd18opdd.py (modified) (8 diffs)
• m/classes/SMBpdd.m (modified) (2 diffs)
• m/classes/SMBpdd.py (modified) (2 diffs)

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

@@ -21 +21 @@
         
         int    smb_model;
-        bool   isdelta18o,ismungsm,isd18opd;
+        bool   isdelta18o,ismungsm,isd18opd,issetpddfac;
         
         /*Update elements: */
@@ -90 +90 @@
                         iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
                         iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
+                        iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
                         iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
                         iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
@@ -97 +98 @@
                                 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
                         }
-
+                        if(issetpddfac){
+                                iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
+                                iomodel->FetchDataToInput(elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
+                        }
                         break;
                 case SMBgradientsEnum:
@@ -137 +141 @@
         int     numoutputs;
         char**  requestedoutputs = NULL;
-        bool    isdelta18o,ismungsm,isd18opd,interp;
+        bool    isdelta18o,ismungsm,isd18opd,issetpddfac,interp;
         int     smb_model;
         IssmDouble *temp = NULL;
@@ -176 +180 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
@@ -208 +213 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
@@ -213 +219 @@
                         iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
                         iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
+                        iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
                         if(isd18opd){
                                 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);

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

@@ -2128 +2128 @@
 }
 /*}}}*/
-void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
+void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
 
         int  numvertices = this->GetNumberOfVertices();
@@ -2144 +2144 @@
         IssmDouble* s0p=xNew<IssmDouble>(numvertices);
         IssmDouble* s0t=xNew<IssmDouble>(numvertices);
+        IssmDouble pddsnowfac = -1.;
+        IssmDouble pddicefac = -1.;
         IssmDouble rho_water,rho_ice,desfac,rlaps,rlapslgm;
         IssmDouble PfacTime,TdiffTime,sealevTime;
@@ -2190 +2192 @@
         }
 
+        /*Recover pdd factors at time t. 
+         *     This parameter is set, if the user wants to define the 
+         *     pdd factors regionally, if issetpddfac==1 in the d18opdd method */
+        if (issetpddfac==1){
+                input=this->GetInput(SmbPddfacSnowEnum); _assert_(input);
+                input2=this->GetInput(SmbPddfacIceEnum); _assert_(input2);
+        }
+
         /*Recover info at the vertices: */
         GetInputListOnVertices(&h[0],ThicknessEnum);
@@ -2198 +2208 @@
         /*measure the surface mass balance*/
         for (int iv = 0; iv<numvertices; iv++){
+                gauss->GaussVertex(iv);
+                pddsnowfac=-1.;
+                pddicefac=-1.;
+                input->GetInputValue(&pddsnowfac,gauss);
+                input2->GetInputValue(&pddicefac,gauss);
                 agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
                                         pdds, pds, &melt[iv], &accu[iv], signorm, yts, h[iv], s[iv],
                                         desfac, s0t[iv], s0p[iv],rlaps,rlapslgm,TdiffTime,sealevTime,
-                                        rho_water,rho_ice);
+                                        pddsnowfac,pddicefac,rho_water,rho_ice);
                 /*Get yearlytemperatures */
                 for(int month=0;month<12;month++) {

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

@@ -140 +140 @@
                 ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
                 ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
-                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
+                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
                 IssmDouble         PureIceEnthalpy(IssmDouble pressure);
                 void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);

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

@@ -188 +188 @@
 
         bool ismungsm;
+        bool issetpddfac;
 
         IssmDouble *pdds    = NULL;
@@ -198 +199 @@
         // Get ismungsm parameter
         femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
+
+        // Get issetpddfac parameter
+        femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
 
         /* initialize PDD (creation of a lookup table)*/
@@ -253 +257 @@
         for(i=0;i<femmodel->elements->Size();i++){
                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                element->PositiveDegreeDay(pdds,pds,signorm,ismungsm);
+                element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
         }
         /*free ressouces: */

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

@@ -4 +4 @@
  */
 
+#include "../io/io.h"
 #include "./elements.h"
 #include "../Numerics/numerics.h"
@@ -11 +12 @@
                                  IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,
                                  IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm,
-                                 IssmDouble TdiffTime,IssmDouble sealevTime,
+                                 IssmDouble TdiffTime,IssmDouble sealevTime, IssmDouble pddsnowfac,IssmDouble pddicefac,
                                  IssmDouble rho_water,IssmDouble rho_ice){
 
@@ -67 +68 @@
   IssmDouble fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
   IssmDouble pddtj, hmx2;
+  IssmDouble pddsnowfac0=4.3, pddicefac0=8.3;
+  IssmDouble snowfac, icefac;
 
   sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
@@ -128 +131 @@
   prect = qmp;     // total precipitation during 1 year taking into account des. ef.
   Tsum=Tsum/3;
+
+  snowfac=pddsnowfac0;
+  icefac=pddicefac0;
+  if (pddsnowfac>=0) {
+          if (pddsnowfac<1.65) {
+                  _printf0_("WARNING: Pdd snow factor input, " << pddsnowfac << ", results in a negative value. It will be ignored. \n");
+          }
+          else{
+                snowfac=pddsnowfac;
+          }
+  }
+  if (pddicefac>=0) {
+          if (pddicefac>17.22) {
+                  _printf0_("WARNING: Pdd ice factor input, " << pddsnowfac << ", results in a negative value. It will be ignored. \n");
+          }
+          else{
+            icefac=pddicefac;
+          }
+  }
   
   /***** 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
+    snwmf=(2.65+snowfac-pddsnowfac0)*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;
+    snwmf = (0.15*(Tsum+1) + (2.65+snowfac-pddsnowfac0))*0.001;
+    smf = (((icefac-pddicefac0)/(Tsum+1))*pow((10.-Tsum),3) + pddicefac0)*0.001;
   }
   else{
-    snwmf=4.3*0.001;
-    smf=8.3*0.001;
+    snwmf=snowfac*0.001;
+    smf=icefac*0.001;
   }
   snwmf=0.95*snwmf;

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

@@ -22 +22 @@
                                  IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
                                  IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
-                                 IssmDouble TdiffTime,IssmDouble sealevTime,
+                                 IssmDouble TdiffTime,IssmDouble sealevTime,IssmDouble pddsnowfac,IssmDouble pddicefac,
                                  IssmDouble rho_water, IssmDouble rho_ice);
 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,

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

@@ -474 +474 @@
         SmbDelta18oSurfaceEnum,
         SmbIsdelta18oEnum,
+        SmbIssetpddfacEnum,
         SmbIsmungsmEnum,
         SmbIsd18opdEnum,
@@ -481 +482 @@
         SmbTemperaturesLgmEnum,
         SmbPrecipitationEnum,
+        SmbPddfacSnowEnum,
+        SmbPddfacIceEnum,
         SmbDesfacEnum,
         SmbS0pEnum,

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

@@ -475 +475 @@
                 case SmbDelta18oSurfaceEnum : return "SmbDelta18oSurface";
                 case SmbIsdelta18oEnum : return "SmbIsdelta18o";
+                case SmbIssetpddfacEnum : return "SmbIssetpddfac";
                 case SmbIsmungsmEnum : return "SmbIsmungsm";
                 case SmbIsd18opdEnum : return "SmbIsd18opd";
@@ -482 +483 @@
                 case SmbTemperaturesLgmEnum : return "SmbTemperaturesLgm";
                 case SmbPrecipitationEnum : return "SmbPrecipitation";
+                case SmbPddfacSnowEnum : return "SmbPddfacSnow";
+                case SmbPddfacIceEnum : return "SmbPddfacIce";
                 case SmbDesfacEnum : return "SmbDesfac";
                 case SmbS0pEnum : return "SmbS0p";

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

@@ -484 +484 @@
               else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
               else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum;
+              else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum;
               else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
               else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum;
@@ -491 +492 @@
               else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
               else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
+              else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum;
+              else if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum;
               else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
               else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
@@ -503 +506 @@
               else if (strcmp(name,"SmbF")==0) return SmbFEnum;
               else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
-              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
-              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
-              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
+              if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
+              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
+              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
+              else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
               else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum;
               else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
@@ -626 +629 @@
               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
-              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
-              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
-              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"J")==0) return JEnum;
+              if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
+              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
+              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
+              else if (strcmp(name,"J")==0) return JEnum;
               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
               else if (strcmp(name,"Step")==0) return StepEnum;
@@ -749 +752 @@
               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
-              else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
-              else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum;
-              else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum;
+              if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
+              else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum;
+              else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum;
+              else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum;
               else if (strcmp(name,"AugmentedLagrangianTheta")==0) return AugmentedLagrangianThetaEnum;
               else if (strcmp(name,"None")==0) return NoneEnum;
@@ -872 +875 @@
               else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
               else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
-              else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
-              else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
-              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
+              if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
+              else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
+              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
+              else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
               else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum;
               else if (strcmp(name,"EsaStrainrateyy")==0) return EsaStrainrateyyEnum;
@@ -995 +998 @@
               else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
               else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
-              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
-              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
-              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
+              if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
+              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
+              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
+              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
               else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
               else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
@@ -1118 +1121 @@
               else if (strcmp(name,"Loads")==0) return LoadsEnum;
               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
-              else if (strcmp(name,"Nodes")==0) return NodesEnum;
-              else if (strcmp(name,"Contours")==0) return ContoursEnum;
-              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Vertices")==0) return VerticesEnum;
+              if (strcmp(name,"Nodes")==0) return NodesEnum;
+              else if (strcmp(name,"Contours")==0) return ContoursEnum;
+              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
+              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
               else if (strcmp(name,"Results")==0) return ResultsEnum;
               else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;

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

@@ -17 +17 @@
                 ismungsm                  = 0;
                 isd18opd                  = 0;
+                issetpddfac               = 0;
                 delta18o                  = NaN;
                 delta18o_surface          = NaN;
                 temperatures_presentday   = NaN;
                 precipitations_presentday = NaN;
+                pddfac_snow               = NaN;
+                pddfac_ice                = NaN;
                 requested_outputs      = {};
         end
@@ -35 +38 @@
                         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
+                        if(self.issetpddfac), self.pddfac_snow=project3d(md,'vector',self.pddfac_snow,'type','node');end
+                        if(self.issetpddfac), self.pddfac_ice=project3d(md,'vector',self.pddfac_ice,'type','node');end
                         self.s0p=project3d(md,'vector',self.s0p,'type','node');
                         self.s0t=project3d(md,'vector',self.s0t,'type','node');
 
                 end % }}}
-                        function list = defaultoutputs(self,md) % {{{
+         function list = defaultoutputs(self,md) % {{{
 
                         list = {''};
@@ -65 +70 @@
                   self.dpermil    = 2.4;
                   self.f          = 0.169;
+                  self.issetpddfac = 0;
                   
                 end % }}}
@@ -86 +92 @@
                                    md = checkfield(md,'fieldname','smb.f','>=',0,'numel',1);
                                 end
+                                if(self.issetpddfac==1)
+                                        md = checkfield(md,'fieldname','smb.pddfac_snow','>=',0,'NaN',1,'Inf',1);
+                                        md = checkfield(md,'fieldname','smb.pddfac_ice','>=',0,'NaN',1,'Inf',1);
+                                end
                         end
                         md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
@@ -94 +104 @@
                         disp(sprintf('\n   PDD and deltaO18 parameters:'));
                         fielddisplay(self,'isd18opd','is delta18o parametrisation from present day temperature and precipitation activated (0 or 1, default is 0)');
+                        fielddisplay(self,'issetpddfac','is user passing in defined pdd factors at each vertex (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]');
@@ -104 +115 @@
                                 fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated');                            
                            fielddisplay(self,'f','precip/temperature scaling factor, required if d18opd is activated');
+                        end
+                        if(self.issetpddfac==1)
+                                fielddisplay(self,'pddfac_snow','Pdd factor for snow, at each vertex [mm ice equiv/day/degree C]');
+                                fielddisplay(self,'pddfac_ice','Pdd factor for ice, at each vertex [mm ice equiv/day/degree C]');
                         end
                         fielddisplay(self,'requested_outputs','additional outputs requested');
@@ -118 +133 @@
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismungsm','format','Boolean');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','isd18opd','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','issetpddfac','format','Boolean');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0p','format','DoubleMat','mattype',1);
@@ -133 +149 @@
                            WriteData(fid,prefix,'object',self,'class','smb','fieldname','f','format','Double');
                         end
+                        if self.issetpddfac==1
+                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_snow','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_ice','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+                        end
                         
                         %process requested outputs

issm/trunk-jpl/src/m/classes/SMBd18opdd.py

@@ -25 +25 @@
                 self.ismungsm                  = 0
                 self.isd18opd                  = 0
+                self.setpddfac                 = 0
                 self.delta18o                  = float('NaN')
                 self.delta18o_surface          = float('NaN')
                 self.temperatures_presentday   = float('NaN')
                 self.precipitations_presentday = float('NaN')
+                self.pddfac_snow               = float('NaN')
+                self.pddfac_ice                = float('NaN')
 
                 #set defaults
@@ -38 +41 @@
 
                 string="%s\n%s"%(string,fielddisplay(self,'isd18opd','is delta18o parametrisation from present day temperature and precipitation activated (0 or 1, default is 0)'))
+                string="%s\n%s"%(string,fielddisplay(self,'issetpddfac','is user passing in defined pdd factors at each vertex (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,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
@@ -47 +51 @@
                         string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o [per mil], required if pdd is activated and delta18o activated'))
                         string="%s\n%s"%(string,fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated'))
+
+                if self.issetpddfac==1:
+                        string="%s\n%s"%(string,fielddisplay(self,'pddfac_snow','Pdd factor for snow, at each vertex [mm ice equiv/day/degree C]'))
+                        string="%s\n%s"%(string,fielddisplay(self,'pddfac_ice','Pdd factor for ice, at each vertex [mm ice equiv/day/degree C]'))
                 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
 
@@ -55 +63 @@
                 if self.isd18opd: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node')
                 if self.isd18opd: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node')
+                if self.issetpddfac: self.pddfac_snow=project3d(md,'vector',self.pddfac_snow,'type','node')
+                if self.issetpddfac: self.pddfac_ice=project3d(md,'vector',self.pddfac_ice,'type','node')
                 self.s0p=project3d(md,'vector',self.s0p,'type','node')
                 self.s0t=project3d(md,'vector',self.s0t,'type','node')
@@ -85 +95 @@
                 self.dpermil    = 2.4
                 self.f          = 0.169
+                self.issetpddfac = 0
                 return self
         #}}}
@@ -107 +118 @@
                                 md = checkfield(md,'fieldname','smb.f','>=',0,'numel',[1])
 
+                        if self.issetpddfac:
+                                md = checkfield(md,'fieldname','smb.pddfac_snow','>=',0,'NaN',1,'Inf',1)
+                                md = checkfield(md,'fieldname','smb.pddfac_ice','>=',0,'NaN',1,'Inf',1)
+
                 md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1)
 
@@ -119 +134 @@
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismungsm','format','Boolean')
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','isd18opd','format','Boolean')
+                WriteData(fid,prefix,'object',self,'class','smb','fieldname','issetpddfac','format','Boolean');
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double')
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0p','format','DoubleMat','mattype',1);
@@ -133 +149 @@
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','dpermil','format','Double')
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','f','format','Double')
+
+                if self.issetpddfac:
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_snow','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_ice','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+
                 #process requested outputs
                 outputs = self.requested_outputs

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

@@ -18 +18 @@
                 isdelta18o                = 0;
                 ismungsm                  = 0;
+                issetpddfac               = 0;
                 delta18o                  = NaN;
                 delta18o_surface          = NaN;
@@ -148 +149 @@
 
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdelta18o','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','issetpddfac','format','Boolean');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismungsm','format','Boolean');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');

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

@@ -26 +26 @@
                 self.isdelta18o                = 0
                 self.ismungsm                  = 0
+                self.issetpddfac               = 0
                 self.delta18o                  = float('NaN')
                 self.delta18o_surface          = float('NaN')
@@ -156 +157 @@
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdelta18o','format','Boolean')
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismungsm','format','Boolean')
+                WriteData(fid,prefix,'object',self,'class','smb','fieldname','issetpddfac','format','Boolean');
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double')
                 WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0p','format','DoubleMat','mattype',1);