Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22449) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22450) @@ -69,4 +69,6 @@ iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum); iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum); + iomodel->FetchDataToInput(elements,"md.smb.aValue",SmbAValueEnum); + iomodel->FetchDataToInput(elements,"md.smb.teValue",SmbTeValueEnum); break; case SMBpddEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22449) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22450) @@ -2556,4 +2556,6 @@ IssmDouble eAir=0.0; IssmDouble pAir=0.0; + IssmDouble teValue=1.0; + IssmDouble aValue=0.0; int aIdx=0; int denIdx=0; @@ -2659,7 +2661,11 @@ Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input); Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input); + Input* teValue_input=this->GetInput(SmbTeValueEnum); _assert_(teValue_input); + Input* aValue_input=this->GetInput(SmbAValueEnum); _assert_(aValue_input); Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input); /*Retrieve input values:*/ Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); + + if (aIdx == 0) aValue_input->GetInputValue(&aValue,gauss); zTop_input->GetInputValue(&zTop,gauss); @@ -2673,4 +2679,5 @@ Tz_input->GetInputValue(&Tz,gauss); Vz_input->GetInputValue(&Vz,gauss); + teValue_input->GetInputValue(&teValue,gauss); isinitialized_input->GetInputValue(&isinitialized); /*}}}*/ @@ -2719,6 +2726,7 @@ a = xNewZeroInit(m); for(int i=0;i(m); for(int i=0;iGetInputValue(&eAir,gauss,t); //screen level vapor pressure [Pa] pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa] + teValue_input->GetInputValue(&teValue,gauss); // screen level air pressure [Pa] + if(aIdx == 0) aValue_input->GetInputValue(&aValue,gauss); // screen level air pressure [Pa] //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); /*}}}*/ @@ -2804,5 +2814,5 @@ /*Snow, firn and ice albedo:*/ - if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid()); + if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce, aSnow,aValue,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid()); @@ -2814,5 +2824,5 @@ /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid()); + if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid()); /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. @@ -2821,5 +2831,5 @@ /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ - if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,rho_ice,this->Sid()); + if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, Ta, P, dzMin, aSnow,rho_ice,this->Sid()); /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K @@ -2832,5 +2842,5 @@ /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every * sub-time step in thermo equations*/ - ulw = 5.67E-8 * pow(T[0],4.0); + ulw = 5.67E-8 * pow(T[0],4.0) * teValue; /*Calculate net longwave [W m-2]*/ @@ -2840,5 +2850,5 @@ if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid()); - /*Verbose some resuls in debug mode: {{{*/ + /*Verbose some results in debug mode: {{{*/ if(VerboseSmb() && 0){ _printf_("smb log: count[" << count << "] m[" << m << "] " @@ -2852,4 +2862,6 @@ << "gsp[" << cellsum(gsp,m) << "] " << "swf[" << netSW << "] " + << "a[" << a << "] " + << "te[" << teValue << "] " << "\n"); } /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22449) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22450) @@ -337,7 +337,8 @@ } /*}}}*/ -void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/ +void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/ //// Calculates Snow, firn and ice albedo as a function of: + // 0 : direct input from aValue parameter // 1 : effective grain radius (Gardner & Sharp, 2009) // 2 : effective grain radius (Brun et al., 2009) @@ -347,4 +348,7 @@ //// Inputs // aIdx = albedo method to use + + // Method 0 + // aValue = direct input value for albedo // Methods 1 & 2 @@ -392,5 +396,8 @@ const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] - if(aIdx==1){ + if (aIdx==0){ + for(int i=0;i0)a[0] = (aSnow * P + a[0] * mInit[0])/mass; re[0] = (reNew * P + re[0] * mInit[0])/mass; gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass; @@ -1787,5 +1796,5 @@ void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/ - //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???] + //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPACTION IS COMPENSATED FOR BY TRACES OF SNOW???] //// FUNCTION INFO Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22449) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22450) @@ -25,8 +25,8 @@ IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid); -void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid); +void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid); void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); -void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid); -void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); +void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid); +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx,IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid); void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22449) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22450) @@ -459,4 +459,6 @@ SmbAEnum, SmbTEnum, + SmbAValueEnum, + SmbTeValueEnum, SmbIsgraingrowthEnum, SmbIsalbedoEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22449) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22450) @@ -461,4 +461,6 @@ case SmbAEnum : return "SmbA"; case SmbTEnum : return "SmbT"; + case SmbAValueEnum : return "SmbAValue"; + case SmbTeValueEnum : return "SmbTeValue"; case SmbIsgraingrowthEnum : return "SmbIsgraingrowth"; case SmbIsalbedoEnum : return "SmbIsalbedo"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22449) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22450) @@ -470,4 +470,6 @@ else if (strcmp(name,"SmbA")==0) return SmbAEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; + else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; + else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; - else if (strcmp(name,"SmbF")==0) return SmbFEnum; - else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum; + 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; @@ -627,10 +629,10 @@ else if (strcmp(name,"GiaW")==0) return GiaWEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; - else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; - else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; + 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; @@ -750,10 +752,10 @@ else if (strcmp(name,"VyObs")==0) return VyObsEnum; else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; - else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; - else if (strcmp(name,"Incremental")==0) return IncrementalEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum; + 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; @@ -873,10 +875,10 @@ else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; else if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum; - else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; - else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum; + 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; @@ -996,10 +998,10 @@ else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum; else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum; - else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; - else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; + 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; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; - else if (strcmp(name,"Loads")==0) return LoadsEnum; - else if (strcmp(name,"Materials")==0) return MaterialsEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Nodes")==0) return NodesEnum; + 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; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 22449) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 22450) @@ -36,4 +36,8 @@ Tz = NaN; %height above ground at which temperature (T) was sampled [m] Vz = NaN; %height above ground at which wind (V) eas sampled [m] + + %optional inputs: + aValue = NaN; %Albedo forcing at every element. Used only if aIdx == 0. + teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) % Initialization of snow properties @@ -51,5 +55,6 @@ %settings: aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1) - % 1: effective grain radius [Gardner & Sharp, 2009] + % 0: direct input from aValue parameter + % 1: effective grain radius [Gardner & Sharp, 2009] % 2: effective grain radius [Brun et al., 2009] % 3: density and cloud amount [Greuell & Konzelmann, 1994] @@ -114,4 +119,12 @@ self.eAir=project3d(md,'vector',self.eAir,'type','node'); self.pAir=project3d(md,'vector',self.pAir,'type','node'); + + if aIdx == 0 & ~isnan(self.aValue) + self.aValue=project3d(md,'vector',self.aValue,'type','node'); + end + if ~isnan(self.teValue) + self.teValue=project3d(md,'vector',self.teValue,'type','node'); + end + end % }}} @@ -150,4 +163,7 @@ self.t0dry = 30; self.K = 7; + + self.teValue = ones(mesh.numberofelements,1); + self.aValue = self.aSnow*ones(mesh.numberofelements,1); self.Dzini=0.05*ones(mesh.numberofelements,2); @@ -167,5 +183,4 @@ function md = checkconsistency(self,md,solution,analyses) % {{{ - md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]); md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]); @@ -189,5 +204,7 @@ md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); - md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]); + md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); + + md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]); md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]); md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]); @@ -201,4 +218,6 @@ switch self.aIdx, + case 0 + md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); case {1 2} md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89); @@ -252,11 +271,14 @@ fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'); fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',... + '0: direct input from aValue parameter',... '1: effective grain radius [Gardner & Sharp, 2009]',... '2: effective grain radius [Brun et al., 2009]',... '3: density and cloud amount [Greuell & Konzelmann, 1994]',... '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}); + + fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'); %snow properties init - fielddisplay(self,'Dzini','Initial cel depth when restart [m]'); + fielddisplay(self,'Dzini','Initial cell depth when restart [m]'); fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]'); fielddisplay(self,'Reini','Initial grain size when restart [mm]'); @@ -272,4 +294,6 @@ %additional albedo parameters: switch self.aIdx + case 0 + fielddisplay(self,'aValue','Albedo forcing at every element. Used only if aIdx == 0.'); case {1 2} fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'); @@ -339,4 +363,7 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double'); + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); %snow properties init Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 22449) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 22450) @@ -42,4 +42,8 @@ Tz = float('NaN') #height above ground at which temperature (T) was sampled [m] Vz = float('NaN') #height above ground at which wind (V) eas sampled [m] + + #optional inputs: + aValue = float('NaN') #Albedo forcing at every element. Used only if aIdx == 0. + teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) # Initialization of snow properties @@ -57,5 +61,6 @@ #settings: aIdx = float('NaN') #method for calculating albedo and subsurface absorption (default is 1) - # 1: effective grain radius [Gardner & Sharp, 2009] + # 0: direct input from aValue parameter + # 1: effective grain radius [Gardner & Sharp, 2009] # 2: effective grain radius [Brun et al., 2009] # 3: density and cloud amount [Greuell & Konzelmann, 1994] @@ -135,11 +140,14 @@ string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)')) string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)', + '0: direct input from aValue parameter', '1: effective grain radius [Gardner & Sharp, 2009]', '2: effective grain radius [Brun et al., 2009]', '3: density and cloud amount [Greuell & Konzelmann, 1994]', '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) + + string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)')) #snow properties init - string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cel depth when restart [m]')) + string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cell depth when restart [m]')) string = "%s\n%s"%(string,fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]')) string = "%s\n%s"%(string,fielddisplay(self,'Reini','Initial grain size when restart [mm]')) @@ -156,4 +164,6 @@ string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)')) string = "%s\n%s"%(string,fielddisplay(self,'aIce','albedo of ice (0.27-0.58)')) + elif self.aIdx == 0: + string = "%s\n%s"%(string,fielddisplay(self,'aValue','Albedo forcing at every element. Used only if aIdx == 0.')) elif self.aIdx == 3: string = "%s\n%s"%(string,fielddisplay(self,'cldFrac','average cloud amount')) @@ -183,4 +193,10 @@ self.eAir = project3d(md,'vector',self.eAir,'type','node') self.pAir = project3d(md,'vector',self.pAir,'type','node') + + if aIdx == 0 and np.isnan(self.aValue): + self.aValue=project3d(md,'vector',self.aValue,'type','node'); + if np.isnan(self.teValue): + self.teValue=project3d(md,'vector',self.teValue,'type','node'); + return self #}}} @@ -211,5 +227,5 @@ self.zY = 1.10*np.ones((mesh.numberofelements,)) self.outputFreq = 30 - + #additional albedo parameters self.aSnow = 0.85 @@ -219,4 +235,7 @@ self.t0dry = 30 self.K = 7 + + self.teValue = np.ones((mesh.numberofelements,)); + self.aValue = self.aSnow*np.ones(mesh.numberofelements,); self.Dzini = 0.05*np.ones((mesh.numberofelements,2)) @@ -257,5 +276,7 @@ md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) - md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]) + md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); + + md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]) md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]) md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]) @@ -271,4 +292,6 @@ md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'> = ',.64,'< = ',.89) md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'> = ',.27,'< = ',.58) + elif self.aIdx == 0: + md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); elif self.aIdx == 3: md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'> = ',0,'< = ',1) @@ -333,4 +356,7 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double') WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double') + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); #snow properties init