source: issm/oecreview/Archive/21724-22754/ISSM-22449-22450.diff@ 22755

../trunk-jpl/src/m/classes/SMBgemb.m

@@ -35 +35 @@
                 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 +54 @@
 
                 %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 +119 @@
                         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 +162 @@
                 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 +182 @@
                 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 +203 @@
                         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 +217 @@
                         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 +270 @@
                         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 +293 @@
             
                         %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 +362 @@
                         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);

../trunk-jpl/src/m/classes/SMBgemb.py

@@ -41 +41 @@
                 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 +60 @@
 
                 #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 +139 @@
                 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 +163 @@
                 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 +192 @@
                 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 +226 @@
                 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 +234 @@
                 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 +275 @@
                 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 +291 @@
                 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 +355 @@
                 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)

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

@@ -68 +68 @@
                         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");

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

@@ -24 +24 @@
 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);

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

@@ -336 +336 @@
         *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 +348 @@
         //// Inputs
         // aIdx      = albedo method to use
 
+        // Method 0
+        //  aValue   = direct input value for albedo
+
         // Methods 1 & 2
         //   re      = surface effective grain radius [mm]
 
@@ -391 +395 @@
         //some constants:
         const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
 
-        if(aIdx==1){ 
+        if (aIdx==0){
+                  for(int i=0;i<m;i++)a[i] = aValue;
+        }
+        else if(aIdx==1){ 
         //function of effective grain radius
         
         //convert effective radius to specific surface area [cm2 g-1]
@@ -399 +406 @@
         
         //determine broadband albedo
         a[0]= 1.48 - pow(S,-0.07);
+                  for(int i=1;i<m;i++)a[i]=a[0];
         }
         else if(aIdx==2){
                 
@@ -420 +428 @@
         
         // broadband surface albedo
         a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
-
+                  for(int i=1;i<m;i++)a[i]=a[0];
         }
         else if(aIdx==3){
                 
@@ -428 +436 @@
         
         // calculate albedo
         a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
+                  for(int i=1;i<m;i++)a[i]=a[0];
         }
         else if(aIdx==4){
                 
@@ -485 +494 @@
         xDelete<IssmDouble>(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 +505 @@
         *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 +816 @@
                 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 +1077 @@
         *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 +1170 @@
                                 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 +1795 @@
 } /*}}}*/ 
 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
 

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

@@ -2555 +2555 @@
         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 +2660 @@
         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 +2678 @@
         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 +2725 @@
                         W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
                         a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
                         T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
-                        /*            /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
-                                                          Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
+                        /*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) 
+                         *    because don't know Tmean yet when set default values.
+                         *    Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
 
                         //fixed lower temperature bounday condition - T is fixed
                         T_bottom=T[m-1];
@@ -2796 +2804 @@
                 P_input->GetInputValue(&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 +2813 @@
                 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 +2823 @@
                 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 +2830 @@
                 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 +2841 @@
 
                 /*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 +2849 @@
                 /*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 +2861 @@
                                                 << "gdn[" << cellsum(gdn,m)  << "] "
                                                 << "gsp[" << cellsum(gsp,m)  << "] "
                                                 << "swf[" << netSW << "] "
+                                                << "a[" << a << "] "
+                                                << "te[" << teValue << "] "
                                                 << "\n");
                 } /*}}}*/
 

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

@@ -458 +458 @@
         SmbWEnum,
         SmbAEnum,
         SmbTEnum,
+        SmbAValueEnum,
+        SmbTeValueEnum,
         SmbIsgraingrowthEnum,
         SmbIsalbedoEnum,
         SmbIsshortwaveEnum,

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

@@ -460 +460 @@
                 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";

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

@@ -469 +469 @@
               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 +505 @@
               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 +628 @@
               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 +751 @@
               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 +874 @@
               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 +997 @@
               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 +1120 @@
               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;