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

Last change on this file was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 34.0 KB

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

     
    3535                C     = NaN; %mean annual snow accumulation [kg m-2 yr-1]
    3636                Tz    = NaN; %height above ground at which temperature (T) was sampled [m]
    3737                Vz    = NaN; %height above ground at which wind (V) eas sampled [m]
     38
     39                %optional inputs:
     40                aValue = NaN; %Albedo forcing at every element.  Used only if aIdx == 0.
     41                teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
    3842       
    3943                % Initialization of snow properties
    4044                Dzini = NaN; %cell depth (m)
    … …  
    5054
    5155                %settings:
    5256                aIdx   = NaN; %method for calculating albedo and subsurface absorption (default is 1)
    53                               % 1: effective grain radius [Gardner & Sharp, 2009]
     57                           % 0: direct input from aValue parameter
     58                           % 1: effective grain radius [Gardner & Sharp, 2009]
    5459                                          % 2: effective grain radius [Brun et al., 2009]
    5560                                          % 3: density and cloud amount [Greuell & Konzelmann, 1994]
    5661                                          % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
    … …  
    114119                        self.eAir=project3d(md,'vector',self.eAir,'type','node');
    115120                        self.pAir=project3d(md,'vector',self.pAir,'type','node');
    116121
     122                        if aIdx == 0 & ~isnan(self.aValue)
     123                                self.aValue=project3d(md,'vector',self.aValue,'type','node');
     124                        end
     125                        if ~isnan(self.teValue)
     126                                self.teValue=project3d(md,'vector',self.teValue,'type','node');
     127                        end
     128
     129
    117130                end % }}}
    118131                function list = defaultoutputs(self,md) % {{{
    119132                        list = {'SmbMassBalance'};
    … …  
    149162                self.t0wet = 15;
    150163                self.t0dry = 30;
    151164                self.K = 7;
     165
     166                self.teValue = ones(mesh.numberofelements,1);
     167                self.aValue = self.aSnow*ones(mesh.numberofelements,1);
    152168       
    153169                self.Dzini=0.05*ones(mesh.numberofelements,2);
    154170                self.Dini=910.0*ones(mesh.numberofelements,2);
    … …  
    166182                end % }}}
    167183                function md = checkconsistency(self,md,solution,analyses) % {{{
    168184
    169 
    170185                        md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
    171186                        md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
    172187                        md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
    … …  
    188203                        md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000);
    189204                        md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000);
    190205
    191                         md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]);
     206                        md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
     207
     208                        md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]);
    192209                        md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
    193210                        md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);
    194211
    … …  
    200217                        md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
    201218
    202219                        switch self.aIdx,
     220                                case 0
     221                                        md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    203222                                case {1 2}
    204223                                        md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
    205224                                        md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
    … …  
    251270                        fielddisplay(self,'InitDensityScaling',{'initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'});
    252271                        fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)');
    253272                        fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
     273                                               '0: direct input from aValue parameter',...
    254274                                                                        '1: effective grain radius [Gardner & Sharp, 2009]',...
    255275                                                                        '2: effective grain radius [Brun et al., 2009]',...
    256276                                                                        '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
    257277                                                                        '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'});
     278
     279                        fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)');
    258280                               
    259281                        %snow properties init
    260                         fielddisplay(self,'Dzini','Initial cel depth when restart [m]');
     282                        fielddisplay(self,'Dzini','Initial cell depth when restart [m]');
    261283                        fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]');
    262284                        fielddisplay(self,'Reini','Initial grain size when restart [mm]');
    263285                        fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]');
    … …  
    271293           
    272294                        %additional albedo parameters:
    273295                        switch self.aIdx
     296                        case 0
     297                                fielddisplay(self,'aValue','Albedo forcing at every element.  Used only if aIdx == 0.');
    274298                        case {1 2}
    275299                                fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)');
    276300                                fielddisplay(self,'aIce','albedo of ice (0.27-0.58)');
    … …  
    338362                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double');
    339363                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double');
    340364                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
     365
     366                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     367                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    341368           
    342369                        %snow properties init
    343370                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3);

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

     
    4141                C     = float('NaN')    #mean annual snow accumulation [kg m-2 yr-1]
    4242                Tz    = float('NaN')    #height above ground at which temperature (T) was sampled [m]
    4343                Vz    = float('NaN')    #height above ground at which wind (V) eas sampled [m]
     44
     45                #optional inputs:
     46                aValue  = float('NaN') #Albedo forcing at every element.  Used only if aIdx == 0.
     47                teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
    4448       
    4549                # Initialization of snow properties
    4650                Dzini = float('NaN')    #cell depth (m)
    … …  
    5660
    5761                #settings:
    5862                aIdx   = float('NaN')   #method for calculating albedo and subsurface absorption (default is 1)
    59                                         # 1: effective grain radius [Gardner & Sharp, 2009]
     63                           # 0: direct input from aValue parameter
     64                                          # 1: effective grain radius [Gardner & Sharp, 2009]
    6065                                          # 2: effective grain radius [Brun et al., 2009]
    6166                                          # 3: density and cloud amount [Greuell & Konzelmann, 1994]
    6267                                          # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
    … …  
    134139                string = "%s\n%s"%(string,fielddisplay(self,'InitDensityScaling',['initial scaling factor multiplying the density of ice','which describes the density of the snowpack.']))
    135140                string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'))
    136141                string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)',
     142                                 '0: direct input from aValue parameter',
    137143                                                '1: effective grain radius [Gardner & Sharp, 2009]',
    138144                                                '2: effective grain radius [Brun et al., 2009]',
    139145                                                '3: density and cloud amount [Greuell & Konzelmann, 1994]',
    140146                                                '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
     147
     148                string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'))
    141149                               
    142150                #snow properties init
    143                 string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cel depth when restart [m]'))
     151                string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cell depth when restart [m]'))
    144152                string = "%s\n%s"%(string,fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]'))
    145153                string = "%s\n%s"%(string,fielddisplay(self,'Reini','Initial grain size when restart [mm]'))
    146154                string = "%s\n%s"%(string,fielddisplay(self,'Gdnini','Initial grain dricity when restart [-]'))
    … …  
    155163                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)):
    156164                        string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'))
    157165                        string = "%s\n%s"%(string,fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'))
     166                elif self.aIdx == 0:
     167                        string = "%s\n%s"%(string,fielddisplay(self,'aValue','Albedo forcing at every element.  Used only if aIdx == 0.'))
    158168                elif self.aIdx == 3:
    159169                        string = "%s\n%s"%(string,fielddisplay(self,'cldFrac','average cloud amount'))
    160170                elif self.aIdx == 4:
    … …  
    182192                self.P = project3d(md,'vector',self.P,'type','node')
    183193                self.eAir = project3d(md,'vector',self.eAir,'type','node')
    184194                self.pAir = project3d(md,'vector',self.pAir,'type','node')
     195
     196                if aIdx == 0 and np.isnan(self.aValue):
     197                        self.aValue=project3d(md,'vector',self.aValue,'type','node');
     198                if np.isnan(self.teValue):
     199                        self.teValue=project3d(md,'vector',self.teValue,'type','node');
     200
    185201                return self
    186202        #}}}
    187203        def defaultoutputs(self,md): # {{{
    … …  
    210226                self.zMin = 130*np.ones((mesh.numberofelements,))
    211227                self.zY = 1.10*np.ones((mesh.numberofelements,))
    212228                self.outputFreq = 30
    213                
     229
    214230                #additional albedo parameters
    215231                self.aSnow = 0.85
    216232                self.aIce = 0.48
    … …  
    218234                self.t0wet = 15
    219235                self.t0dry = 30
    220236                self.K = 7
     237
     238                self.teValue = np.ones((mesh.numberofelements,));
     239                self.aValue = self.aSnow*np.ones(mesh.numberofelements,);
    221240       
    222241                self.Dzini = 0.05*np.ones((mesh.numberofelements,2))
    223242                self.Dini = 910.0*np.ones((mesh.numberofelements,2))
    … …  
    256275                md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000)
    257276                md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000)
    258277
    259                 md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4])
     278                md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
     279
     280                md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4])
    260281                md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1])
    261282                md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5])
    262283
    … …  
    270291                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)):
    271292                        md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'> = ',.64,'< = ',.89)
    272293                        md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'> = ',.27,'< = ',.58)
     294                elif self.aIdx == 0:
     295                        md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    273296                elif self.aIdx == 3:
    274297                        md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'> = ',0,'< = ',1)
    275298                elif self.aIdx == 4:
    … …  
    332355                WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double')
    333356                WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double')
    334357                WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double')
     358
     359                WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     360                WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    335361           
    336362                #snow properties init
    337363                WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3)

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

     
    6868                        iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum);
    6969                        iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum);
    7070                        iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum);
     71                        iomodel->FetchDataToInput(elements,"md.smb.aValue",SmbAValueEnum);
     72                        iomodel->FetchDataToInput(elements,"md.smb.teValue",SmbTeValueEnum);
    7173                        break;
    7274                case SMBpddEnum:
    7375                        iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");

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

     
    2424void       GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY);
    2525IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
    2626void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
    27 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);
     27void 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);
    2828void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
    29 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);
    30 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);
     29void 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);
     30void 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);
    3131void 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);
    3232void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
    3333void 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

     
    336336        *pgsp=gsp;
    337337
    338338}  /*}}}*/
    339 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) { /*{{{*/
     339void 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) { /*{{{*/
    340340
    341341        //// Calculates Snow, firn and ice albedo as a function of:
     342        //   0 : direct input from aValue parameter
    342343        //   1 : effective grain radius (Gardner & Sharp, 2009)
    343344        //   2 : effective grain radius (Brun et al., 2009)
    344345        //   3 : density and cloud amount (Greuell & Konzelmann, 1994)
    … …  
    347348        //// Inputs
    348349        // aIdx      = albedo method to use
    349350
     351        // Method 0
     352        //  aValue   = direct input value for albedo
     353
    350354        // Methods 1 & 2
    351355        //   re      = surface effective grain radius [mm]
    352356
    … …  
    391395        //some constants:
    392396        const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
    393397
    394         if(aIdx==1){
     398        if (aIdx==0){
     399                  for(int i=0;i<m;i++)a[i] = aValue;
     400        }
     401        else if(aIdx==1){
    395402        //function of effective grain radius
    396403       
    397404        //convert effective radius to specific surface area [cm2 g-1]
    … …  
    399406       
    400407        //determine broadband albedo
    401408        a[0]= 1.48 - pow(S,-0.07);
     409                  for(int i=1;i<m;i++)a[i]=a[0];
    402410        }
    403411        else if(aIdx==2){
    404412               
    … …  
    420428       
    421429        // broadband surface albedo
    422430        a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
    423 
     431                  for(int i=1;i<m;i++)a[i]=a[0];
    424432        }
    425433        else if(aIdx==3){
    426434               
    … …  
    428436       
    429437        // calculate albedo
    430438        a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
     439                  for(int i=1;i<m;i++)a[i]=a[0];
    431440        }
    432441        else if(aIdx==4){
    433442               
    … …  
    485494        xDelete<IssmDouble>(d_a);
    486495
    487496        }
    488         else _error_("albedo method switch should range from 1 to 4!");
     497        else _error_("albedo method switch should range from 0 to 4!");
    489498       
    490499        // Check for erroneous values
    491500        if (a[0] > 1) _printf_("albedo > 1.0\n");
    … …  
    496505        *pa=a;
    497506
    498507}  /*}}}*/
    499 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) { /*{{{*/
     508void 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) { /*{{{*/
    500509
    501510        /* ENGLACIAL THERMODYNAMICS*/
    502511         
    … …  
    807816                dT_turb = turb  / TCs;
    808817
    809818                // upward longwave contribution
    810                 ulw = - SB * pow(Ts,4.0) * dt;
     819                ulw = - SB * pow(Ts,4.0)* teValue * dt;
    811820                dT_ulw = ulw / TCs;
    812821               
    813822                // new grid point temperature
    … …  
    10681077        *pswf=swf;
    10691078
    10701079} /*}}}*/
    1071 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){ /*{{{*/
     1080void 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){ /*{{{*/
    10721081
    10731082        // Adds precipitation and deposition to the model grid
    10741083
    … …  
    11611170                                T[0] = (T_air * P + T[0] * mInit[0])/mass;
    11621171
    11631172                                // adjust a, re, gdn & gsp
    1164                                 a[0] = (aSnow * P + a[0] * mInit[0])/mass;
     1173                                if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
    11651174                                re[0] = (reNew * P + re[0] * mInit[0])/mass;
    11661175                                gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass;
    11671176                                gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass;
    … …  
    17861795} /*}}}*/
    17871796void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/
    17881797
    1789         //// 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???]
     1798        //// 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???]
    17901799
    17911800        //// FUNCTION INFO
    17921801

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

     
    25552555        IssmDouble P=0.0;
    25562556        IssmDouble eAir=0.0;
    25572557        IssmDouble pAir=0.0;
     2558        IssmDouble teValue=1.0;
     2559        IssmDouble aValue=0.0;
    25582560        int        aIdx=0;
    25592561        int        denIdx=0;
    25602562        int        swIdx=0;
    … …  
    26582660        Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input);
    26592661        Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input);
    26602662        Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input);
     2663        Input* teValue_input=this->GetInput(SmbTeValueEnum); _assert_(teValue_input);
     2664        Input* aValue_input=this->GetInput(SmbAValueEnum); _assert_(aValue_input);
    26612665        Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input);
    26622666        /*Retrieve input values:*/
    26632667        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
    26642668
     2669        if (aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);
     2670
    26652671        zTop_input->GetInputValue(&zTop,gauss);
    26662672        dzTop_input->GetInputValue(&dzTop,gauss);
    26672673        dzMin_input->GetInputValue(&dzMin,gauss);
    … …  
    26722678        C_input->GetInputValue(&C,gauss);
    26732679        Tz_input->GetInputValue(&Tz,gauss);
    26742680        Vz_input->GetInputValue(&Vz,gauss);
     2681        teValue_input->GetInputValue(&teValue,gauss);
    26752682        isinitialized_input->GetInputValue(&isinitialized);
    26762683        /*}}}*/
    26772684
    … …  
    27182725                        W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
    27192726                        a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
    27202727                        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]
    2721                         /*            /!\ 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.
    2722                                                           Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
     2728                        /*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties)
     2729                         *    because don't know Tmean yet when set default values.
     2730                         *    Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
    27232731
    27242732                        //fixed lower temperature bounday condition - T is fixed
    27252733                        T_bottom=T[m-1];
    … …  
    27962804                P_input->GetInputValue(&P,gauss,t);        //precipitation [kg m-2]
    27972805                eAir_input->GetInputValue(&eAir,gauss,t);  //screen level vapor pressure [Pa]
    27982806                pAir_input->GetInputValue(&pAir,gauss,t);  // screen level air pressure [Pa]
     2807                teValue_input->GetInputValue(&teValue,gauss);  // screen level air pressure [Pa]
     2808                if(aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);  // screen level air pressure [Pa]
    27992809                //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
    28002810                /*}}}*/
    28012811
    … …  
    28032813                if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
    28042814
    28052815                /*Snow, firn and ice albedo:*/
    2806                 if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
     2816                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());
    28072817
    28082818
    28092819                /*Distribution of absorbed short wave radation with depth:*/
    … …  
    28132823                netSW = cellsum(swf,m);
    28142824
    28152825                /*Thermal profile computation:*/
    2816                 if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid());
     2826                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());
    28172827
    28182828                /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
    28192829                 * need to fix this in case all or more of cell evaporates */
    … …  
    28202830                dz[0] = dz[0] + EC / d[0];
    28212831
    28222832                /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
    2823                 if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
     2833                if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
    28242834
    28252835                /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    28262836                 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    … …  
    28312841
    28322842                /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
    28332843                 * sub-time step in thermo equations*/
    2834                 ulw = 5.67E-8 * pow(T[0],4.0);
     2844                ulw = 5.67E-8 * pow(T[0],4.0) * teValue;
    28352845
    28362846                /*Calculate net longwave [W m-2]*/
    28372847                netLW = dlw - ulw;
    … …  
    28392849                /*Calculate turbulent heat fluxes [W m-2]*/
    28402850                if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
    28412851
    2842                 /*Verbose some resuls in debug mode: {{{*/
     2852                /*Verbose some results in debug mode: {{{*/
    28432853                if(VerboseSmb() && 0){
    28442854                        _printf_("smb log: count[" << count << "] m[" << m << "] "
    28452855                                                << setprecision(16)   << "T[" << cellsum(T,m)  << "] "
    … …  
    28512861                                                << "gdn[" << cellsum(gdn,m)  << "] "
    28522862                                                << "gsp[" << cellsum(gsp,m)  << "] "
    28532863                                                << "swf[" << netSW << "] "
     2864                                                << "a[" << a << "] "
     2865                                                << "te[" << teValue << "] "
    28542866                                                << "\n");
    28552867                } /*}}}*/
    28562868

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

     
    458458        SmbWEnum,
    459459        SmbAEnum,
    460460        SmbTEnum,
     461        SmbAValueEnum,
     462        SmbTeValueEnum,
    461463        SmbIsgraingrowthEnum,
    462464        SmbIsalbedoEnum,
    463465        SmbIsshortwaveEnum,

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

     
    460460                case SmbWEnum : return "SmbW";
    461461                case SmbAEnum : return "SmbA";
    462462                case SmbTEnum : return "SmbT";
     463                case SmbAValueEnum : return "SmbAValue";
     464                case SmbTeValueEnum : return "SmbTeValue";
    463465                case SmbIsgraingrowthEnum : return "SmbIsgraingrowth";
    464466                case SmbIsalbedoEnum : return "SmbIsalbedo";
    465467                case SmbIsshortwaveEnum : return "SmbIsshortwave";

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

     
    469469              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
    470470              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    471471              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
     472              else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
     473              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    472474              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
    473475              else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
    474476              else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum;
    … …  
    503505              else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
    504506              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    505507              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    506               else if (strcmp(name,"SmbF")==0) return SmbFEnum;
    507               else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
     511              if (strcmp(name,"SmbF")==0) return SmbFEnum;
     512              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
     513              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
    512514              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
    513515              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    514516              else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
    … …  
    626628              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
    627629              else if (strcmp(name,"GiaW")==0) return GiaWEnum;
    628630              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    629               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    630               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     634              if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     635              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     636              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    635637              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    636638              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
    637639              else if (strcmp(name,"J")==0) return JEnum;
    … …  
    749751              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    750752              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    751753              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
    752               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    753               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
     757              if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
     758              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     759              else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
    758760              else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum;
    759761              else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum;
    760762              else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum;
    … …  
    872874              else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    873875              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    874876              else if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum;
    875               else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
    876               else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
     880              if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
     881              else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
     882              else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
    881883              else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
    882884              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
    883885              else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
    … …  
    995997              else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
    996998              else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
    997999              else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
    998               else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    999               else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
     1003              if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
     1004              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
     1005              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    10041006              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    10051007              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
    10061008              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
    … …  
    11181120              else if (strcmp(name,"Water")==0) return WaterEnum;
    11191121              else if (strcmp(name,"DataSet")==0) return DataSetEnum;
    11201122              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    1121               else if (strcmp(name,"Loads")==0) return LoadsEnum;
    1122               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Nodes")==0) return NodesEnum;
     1126              if (strcmp(name,"Loads")==0) return LoadsEnum;
     1127              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
     1128              else if (strcmp(name,"Nodes")==0) return NodesEnum;
    11271129              else if (strcmp(name,"Contours")==0) return ContoursEnum;
    11281130              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
    11291131              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
Note: See TracBrowser for help on using the repository browser.