Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22448) @@ -21,5 +21,5 @@ int smb_model; - bool isdelta18o,ismungsm,isd18opd; + bool isdelta18o,ismungsm,isd18opd,issetpddfac; /*Update elements: */ @@ -90,4 +90,5 @@ 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,5 +98,8 @@ 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,5 +141,5 @@ int numoutputs; char** requestedoutputs = NULL; - bool isdelta18o,ismungsm,isd18opd,interp; + bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp; int smb_model; IssmDouble *temp = NULL; @@ -176,4 +180,5 @@ 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,4 +213,5 @@ 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,4 +219,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22448) @@ -2128,5 +2128,5 @@ } /*}}}*/ -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,4 +2144,6 @@ IssmDouble* s0p=xNew(numvertices); IssmDouble* s0t=xNew(numvertices); + IssmDouble pddsnowfac = -1.; + IssmDouble pddicefac = -1.; IssmDouble rho_water,rho_ice,desfac,rlaps,rlapslgm; IssmDouble PfacTime,TdiffTime,sealevTime; @@ -2190,4 +2192,12 @@ } + /*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,8 +2208,13 @@ /*measure the surface mass balance*/ for (int iv = 0; ivGaussVertex(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++) { Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22447) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22448) @@ -140,5 +140,5 @@ 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); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 22448) @@ -188,4 +188,5 @@ bool ismungsm; + bool issetpddfac; IssmDouble *pdds = NULL; @@ -198,4 +199,7 @@ // Get ismungsm parameter femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum); + + // Get issetpddfac parameter + femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum); /* initialize PDD (creation of a lookup table)*/ @@ -253,5 +257,5 @@ for(i=0;ielements->Size();i++){ element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->PositiveDegreeDay(pdds,pds,signorm,ismungsm); + element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac); } /*free ressouces: */ Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 22448) @@ -4,4 +4,5 @@ */ +#include "../io/io.h" #include "./elements.h" #include "../Numerics/numerics.h" @@ -11,5 +12,5 @@ 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,4 +68,6 @@ 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,17 +131,36 @@ 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; Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 22447) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 22448) @@ -22,5 +22,5 @@ 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, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22447) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22448) @@ -474,4 +474,5 @@ SmbDelta18oSurfaceEnum, SmbIsdelta18oEnum, + SmbIssetpddfacEnum, SmbIsmungsmEnum, SmbIsd18opdEnum, @@ -481,4 +482,6 @@ SmbTemperaturesLgmEnum, SmbPrecipitationEnum, + SmbPddfacSnowEnum, + SmbPddfacIceEnum, SmbDesfacEnum, SmbS0pEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22448) @@ -475,4 +475,5 @@ case SmbDelta18oSurfaceEnum : return "SmbDelta18oSurface"; case SmbIsdelta18oEnum : return "SmbIsdelta18o"; + case SmbIssetpddfacEnum : return "SmbIssetpddfac"; case SmbIsmungsmEnum : return "SmbIsmungsm"; case SmbIsd18opdEnum : return "SmbIsd18opd"; @@ -482,4 +483,6 @@ case SmbTemperaturesLgmEnum : return "SmbTemperaturesLgm"; case SmbPrecipitationEnum : return "SmbPrecipitation"; + case SmbPddfacSnowEnum : return "SmbPddfacSnow"; + case SmbPddfacIceEnum : return "SmbPddfacIce"; case SmbDesfacEnum : return "SmbDesfac"; case SmbS0pEnum : return "SmbS0p"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22447) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22448) @@ -484,4 +484,5 @@ 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,4 +492,6 @@ 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,11 +506,11 @@ 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,11 +629,11 @@ 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,11 +752,11 @@ 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,11 +875,11 @@ 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,11 +998,11 @@ 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,11 +1121,11 @@ 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; Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22447) +++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22448) @@ -17,8 +17,11 @@ 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,9 +38,11 @@ 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,4 +70,5 @@ self.dpermil = 2.4; self.f = 0.169; + self.issetpddfac = 0; end % }}} @@ -86,4 +92,8 @@ 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,4 +104,5 @@ 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,4 +115,8 @@ 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,4 +133,5 @@ 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,4 +149,8 @@ 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 Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 22447) +++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 22448) @@ -25,8 +25,11 @@ 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,4 +41,5 @@ 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,4 +51,8 @@ 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,4 +63,6 @@ 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,4 +95,5 @@ self.dpermil = 2.4 self.f = 0.169 + self.issetpddfac = 0 return self #}}} @@ -107,4 +118,8 @@ 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,4 +134,5 @@ 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,4 +149,9 @@ 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 Index: /issm/trunk-jpl/src/m/classes/SMBpdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 22447) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 22448) @@ -18,4 +18,5 @@ isdelta18o = 0; ismungsm = 0; + issetpddfac = 0; delta18o = NaN; delta18o_surface = NaN; @@ -148,4 +149,5 @@ 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'); Index: /issm/trunk-jpl/src/m/classes/SMBpdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 22447) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 22448) @@ -26,4 +26,5 @@ self.isdelta18o = 0 self.ismungsm = 0 + self.issetpddfac = 0 self.delta18o = float('NaN') self.delta18o_surface = float('NaN') @@ -156,4 +157,5 @@ 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);