Index: ../trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- ../trunk-jpl/src/m/classes/SMBgemb.m (revision 22449) +++ ../trunk-jpl/src/m/classes/SMBgemb.m (revision 22450) @@ -35,6 +35,10 @@ C = NaN; %mean annual snow accumulation [kg m-2 yr-1] 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 Dzini = NaN; %cell depth (m) @@ -50,7 +54,8 @@ %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] % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] @@ -114,6 +119,14 @@ 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 % }}} function list = defaultoutputs(self,md) % {{{ list = {'SmbMassBalance'}; @@ -149,6 +162,9 @@ self.t0wet = 15; 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); self.Dini=910.0*ones(mesh.numberofelements,2); @@ -166,7 +182,6 @@ end % }}} 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]); md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]); @@ -188,7 +203,9 @@ md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 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]); @@ -200,6 +217,8 @@ md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1); 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); md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58); @@ -251,13 +270,16 @@ fielddisplay(self,'InitDensityScaling',{'initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'}); 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]'); fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]'); @@ -271,6 +293,8 @@ %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)'); fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'); @@ -338,6 +362,9 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double'); 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 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3); Index: ../trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- ../trunk-jpl/src/m/classes/SMBgemb.py (revision 22449) +++ ../trunk-jpl/src/m/classes/SMBgemb.py (revision 22450) @@ -41,6 +41,10 @@ C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1] 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 Dzini = float('NaN') #cell depth (m) @@ -56,7 +60,8 @@ #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] # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] @@ -134,13 +139,16 @@ string = "%s\n%s"%(string,fielddisplay(self,'InitDensityScaling',['initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'])) 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]')) string = "%s\n%s"%(string,fielddisplay(self,'Gdnini','Initial grain dricity when restart [-]')) @@ -155,6 +163,8 @@ if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)): 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')) elif self.aIdx == 4: @@ -182,6 +192,12 @@ self.P = project3d(md,'vector',self.P,'type','node') 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 #}}} def defaultoutputs(self,md): # {{{ @@ -210,7 +226,7 @@ self.zMin = 130*np.ones((mesh.numberofelements,)) self.zY = 1.10*np.ones((mesh.numberofelements,)) self.outputFreq = 30 - + #additional albedo parameters self.aSnow = 0.85 self.aIce = 0.48 @@ -218,6 +234,9 @@ self.t0wet = 15 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)) self.Dini = 910.0*np.ones((mesh.numberofelements,2)) @@ -256,7 +275,9 @@ md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) 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]) @@ -270,6 +291,8 @@ if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)): 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) elif self.aIdx == 4: @@ -332,6 +355,9 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double') 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 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3) Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22449) +++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22450) @@ -68,6 +68,8 @@ iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum); 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: iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o"); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22449) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22450) @@ -24,10 +24,10 @@ void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY); 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); void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22449) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22450) @@ -336,9 +336,10 @@ *pgsp=gsp; } /*}}}*/ -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) // 3 : density and cloud amount (Greuell & Konzelmann, 1994) @@ -347,6 +348,9 @@ //// Inputs // aIdx = albedo method to use + // Method 0 + // aValue = direct input value for albedo + // Methods 1 & 2 // re = surface effective grain radius [mm] @@ -391,7 +395,10 @@ //some constants: const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] - if(aIdx==1){ + if (aIdx==0){ + for(int i=0;i(d_a); } - else _error_("albedo method switch should range from 1 to 4!"); + else _error_("albedo method switch should range from 0 to 4!"); // Check for erroneous values if (a[0] > 1) _printf_("albedo > 1.0\n"); @@ -496,7 +505,7 @@ *pa=a; } /*}}}*/ -void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/ +void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/ /* ENGLACIAL THERMODYNAMICS*/ @@ -807,7 +816,7 @@ dT_turb = turb / TCs; // upward longwave contribution - ulw = - SB * pow(Ts,4.0) * dt; + ulw = - SB * pow(Ts,4.0)* teValue * dt; dT_ulw = ulw / TCs; // new grid point temperature @@ -1068,7 +1077,7 @@ *pswf=swf; } /*}}}*/ -void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, 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 T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/ // Adds precipitation and deposition to the model grid @@ -1161,7 +1170,7 @@ T[0] = (T_air * P + T[0] * mInit[0])/mass; // adjust a, re, gdn & gsp - a[0] = (aSnow * P + a[0] * mInit[0])/mass; + if(aIdx>0)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; gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass; @@ -1786,7 +1795,7 @@ } /*}}}*/ 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: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22449) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22450) @@ -2555,6 +2555,8 @@ IssmDouble P=0.0; IssmDouble eAir=0.0; IssmDouble pAir=0.0; + IssmDouble teValue=1.0; + IssmDouble aValue=0.0; int aIdx=0; int denIdx=0; int swIdx=0; @@ -2658,10 +2660,14 @@ Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input); 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); dzTop_input->GetInputValue(&dzTop,gauss); dzMin_input->GetInputValue(&dzMin,gauss); @@ -2672,6 +2678,7 @@ C_input->GetInputValue(&C,gauss); Tz_input->GetInputValue(&Tz,gauss); Vz_input->GetInputValue(&Vz,gauss); + teValue_input->GetInputValue(&teValue,gauss); isinitialized_input->GetInputValue(&isinitialized); /*}}}*/ @@ -2718,8 +2725,9 @@ W = xNewZeroInit(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;iGetInputValue(&P,gauss,t); //precipitation [kg m-2] eAir_input->GetInputValue(&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"); /*}}}*/ @@ -2803,7 +2813,7 @@ if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid()); /*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()); /*Distribution of absorbed short wave radation with depth:*/ @@ -2813,7 +2823,7 @@ netSW = cellsum(swf,m); /*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. * need to fix this in case all or more of cell evaporates */ @@ -2820,7 +2830,7 @@ dz[0] = dz[0] + EC / d[0]; /*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 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ @@ -2831,7 +2841,7 @@ /*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]*/ netLW = dlw - ulw; @@ -2839,7 +2849,7 @@ /*Calculate turbulent heat fluxes [W m-2]*/ 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 << "] " << setprecision(16) << "T[" << cellsum(T,m) << "] " @@ -2851,6 +2861,8 @@ << "gdn[" << cellsum(gdn,m) << "] " << "gsp[" << cellsum(gsp,m) << "] " << "swf[" << netSW << "] " + << "a[" << a << "] " + << "te[" << teValue << "] " << "\n"); } /*}}}*/ Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22449) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22450) @@ -458,6 +458,8 @@ SmbWEnum, SmbAEnum, SmbTEnum, + SmbAValueEnum, + SmbTeValueEnum, SmbIsgraingrowthEnum, SmbIsalbedoEnum, SmbIsshortwaveEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22449) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22450) @@ -460,6 +460,8 @@ case SmbWEnum : return "SmbW"; case SmbAEnum : return "SmbA"; case SmbTEnum : return "SmbT"; + case SmbAValueEnum : return "SmbAValue"; + case SmbTeValueEnum : return "SmbTeValue"; case SmbIsgraingrowthEnum : return "SmbIsgraingrowth"; case SmbIsalbedoEnum : return "SmbIsalbedo"; case SmbIsshortwaveEnum : return "SmbIsshortwave"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22449) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22450) @@ -469,6 +469,8 @@ else if (strcmp(name,"SmbW")==0) return SmbWEnum; 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; else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; @@ -503,12 +505,12 @@ else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum; 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; else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; @@ -626,12 +628,12 @@ else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; 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; else if (strcmp(name,"J")==0) return JEnum; @@ -749,12 +751,12 @@ else if (strcmp(name,"VxObs")==0) return VxObsEnum; 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; else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum; @@ -872,12 +874,12 @@ else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; 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; else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum; @@ -995,12 +997,12 @@ else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum; 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; else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; @@ -1118,12 +1120,12 @@ else if (strcmp(name,"Water")==0) return WaterEnum; 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; else if (strcmp(name,"Vertices")==0) return VerticesEnum;