Changeset 22482


Ignore:
Timestamp:
02/27/18 19:33:05 (7 years ago)
Author:
schlegel
Message:

CHG: add option for forcing albedo only above a specific surface density

Location:
issm/trunk-jpl
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r22475 r22482  
    179179                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
    180180                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
     181                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
    181182                        break;
    182183                case SMBpddEnum:
  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r22475 r22482  
    26062606        IssmDouble init_scaling=0.0;
    26072607        IssmDouble thermo_scaling=1.0;
     2608        IssmDouble adThresh=1023.0;
    26082609
    26092610        /*}}}*/
     
    26682669        parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
    26692670        parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum);
     2671        parameters->FindParam(&adThresh,SmbAdThreshEnum);
    26702672
    26712673        /*}}}*/
     
    26942696        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
    26952697
    2696         if (aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);
    2697 
    26982698        zTop_input->GetInputValue(&zTop,gauss);
    26992699        dzTop_input->GetInputValue(&dzTop,gauss);
     
    27072707        Vz_input->GetInputValue(&Vz,gauss);
    27082708        teValue_input->GetInputValue(&teValue,gauss);
     2709        aValue_input->GetInputValue(&aValue,gauss);
    27092710        isinitialized_input->GetInputValue(&isinitialized);
    27102711        /*}}}*/
     
    28332834                pAir_input->GetInputValue(&pAir,gauss,t);  // screen level air pressure [Pa]
    28342835                teValue_input->GetInputValue(&teValue,gauss);  // screen level air pressure [Pa]
    2835                 if(aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);  // screen level air pressure [Pa]
     2836                aValue_input->GetInputValue(&aValue,gauss);  // screen level air pressure [Pa]
    28362837                //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
    28372838                /*}}}*/
     
    28412842
    28422843                /*Snow, firn and ice albedo:*/
    2843                 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());
     2844                if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    28442845
    28452846
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r22475 r22482  
    337337
    338338}  /*}}}*/
    339 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) { /*{{{*/
     339void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, 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:
     
    345345        //   3 : density and cloud amount (Greuell & Konzelmann, 1994)
    346346        //   4 : exponential time decay & wetness (Bougamont & Bamber, 2005)
    347         //   5 : ingest MODIS mode, direct input from aValue parameter applied to surface ice only
    348347
    349348        //// Inputs
    350349        // aIdx      = albedo method to use
    351 
    352         // Method 0 & 5
    353         //  aValue   = direct input value for albedo
     350       
     351        // Method 0
     352        //  aValue   = direct input value for albedo, override all changes to albedo
     353
     354        // adThresh
     355        //  Apply below method to all areas with densities below this value,
     356        //  or else apply direct input value, allowing albedo to be altered. 
     357        //  Default value is rho water (1023 kg m-3).
    354358
    355359        // Methods 1 & 2
     
    397401        const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
    398402
    399         if (aIdx==0){
    400                   for(int i=0;i<m;i++)a[i] = aValue;
    401         }
    402         else if(aIdx==1){
    403         //function of effective grain radius
    404        
    405         //convert effective radius to specific surface area [cm2 g-1]
    406         IssmDouble S = 3.0 / (0.091 * re[0]);
    407        
    408         //determine broadband albedo
    409         a[0]= 1.48 - pow(S,-0.07);
    410                   for(int i=1;i<m;i++)a[i]=a[0];
    411         }
    412         else if(aIdx==2){
    413                
    414         // Spectral fractions  (Lefebre et al., 2003)
    415         // [0.3-0.8um 0.8-1.5um 1.5-2.8um]
    416        
    417         IssmDouble sF[3] = {0.606, 0.301, 0.093};
    418        
    419         // convert effective radius to grain size in meters
    420         IssmDouble gsz = (re[0] * 2.0) / 1000.0;
    421        
    422         // spectral range:
    423         // 0.3 - 0.8um
    424         IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5));
    425         // 0.8 - 1.5um
    426         IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5));
    427         // 1.5 - 2.8um
    428         IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
    429        
    430         // broadband surface albedo
    431         a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
    432                   for(int i=1;i<m;i++)a[i]=a[0];
    433         }
    434         else if(aIdx==3){
    435                
    436         // a as a function of density
    437        
    438         // calculate albedo
    439         a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
    440                   for(int i=1;i<m;i++)a[i]=a[0];
    441         }
    442         else if(aIdx==4){
    443                
    444         // exponential time decay & wetness
    445        
    446         // change in albedo with time:
    447         //   (d_a) = (a - a_old)/(t0)
    448         // where: t0 = timescale for albedo decay
    449        
    450         dt = dt / dts;    // convert from [s] to [d]
    451        
    452         // initialize variables
    453         IssmDouble* t0=xNew<IssmDouble>(m);
    454         IssmDouble* T=xNew<IssmDouble>(m);
    455         IssmDouble* t0warm=xNew<IssmDouble>(m);
    456         IssmDouble* d_a=xNew<IssmDouble>(m);
    457        
    458         // specify constants
    459         // a_wet = 0.15;        // water albedo (0.15)
    460         // a_new = aSnow        // new snow albedo (0.64 - 0.89)
    461         // a_old = aIce;        // old snow/ice albedo (0.27-0.53)
    462         // t0_wet = t0wet;      // time scale for wet snow (15-21.9) [d]
    463         // t0_dry = t0dry;      // warm snow timescale [15] [d]
    464         // K = 7                // time scale temperature coef. (7) [d]
    465         // W0 = 300;            // 200 - 600 [mm]
    466         const IssmDouble z_snow = 15.0;            // 16 - 32 [mm]
    467        
    468         // determine timescale for albedo decay
    469         for(int i=0;i<m;i++)if(W[i]>0.0+Wtol)t0[i]=t0wet; // wet snow timescale
    470         for(int i=0;i<m;i++)T[i]=TK[i] - CtoK; // change T from K to degC
    471         for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
    472         for(int i=0;i<m;i++)if(fabs(W[i])<Wtol && T[i]>=-10.0-Ttol)t0[i]= t0warm[i];
    473         for(int i=0;i<m;i++)if(T[i]<-10.0-Ttol) t0[i] =  10.0 * K + t0dry; // 'cold' snow timescale
    474        
    475         // calculate new albedo
    476         for(int i=0;i<m;i++)d_a[i] = (a[i] - aIce) / t0[i] * dt;           // change in albedo
    477         for(int i=0;i<m;i++)a[i] -= d_a[i];                            // new albedo
    478        
    479         // modification of albedo due to thin layer of snow or solid
    480         // condensation (deposition) at the surface surface
    481        
    482         // check if condensation occurs & if it is deposited in solid phase
    483         if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]
    484        
    485         a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);
    486        
    487         //----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------
    488         // modification of albedo due to thin layer of water on the surface
    489         // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);
    490 
    491         /*Free ressources:*/
    492         xDelete<IssmDouble>(t0);
    493         xDelete<IssmDouble>(T);
    494         xDelete<IssmDouble>(t0warm);
    495         xDelete<IssmDouble>(d_a);
    496 
    497         }
    498         else if(aIdx==5){
    499                 for(int i=0;i<m;i++)if(dIce - d[i]<Dtol) a[i] = aValue;
    500         }
    501         else _error_("albedo method switch should range from 0 to 5!");
    502        
     403        if(aIdx==0 | (adThresh - d[0]<Dtol)){
     404                a[0] = aValue;
     405        }
     406        else{
     407                if(aIdx==1){
     408                        //function of effective grain radius
     409
     410                        //convert effective radius to specific surface area [cm2 g-1]
     411                        IssmDouble S = 3.0 / (0.091 * re[0]);
     412
     413                        //determine broadband albedo
     414                        a[0]= 1.48 - pow(S,-0.07);
     415                }
     416                else if(aIdx==2){
     417
     418                        // Spectral fractions  (Lefebre et al., 2003)
     419                        // [0.3-0.8um 0.8-1.5um 1.5-2.8um]
     420
     421                        IssmDouble sF[3] = {0.606, 0.301, 0.093};
     422
     423                        // convert effective radius to grain size in meters
     424                        IssmDouble gsz = (re[0] * 2.0) / 1000.0;
     425
     426                        // spectral range:
     427                        // 0.3 - 0.8um
     428                        IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5));
     429                        // 0.8 - 1.5um
     430                        IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5));
     431                        // 1.5 - 2.8um
     432                        IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
     433
     434                        // broadband surface albedo
     435                        a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
     436                }
     437                else if(aIdx==3){
     438
     439                        // a as a function of density
     440
     441                        // calculate albedo
     442                        a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
     443                }
     444                else if(aIdx==4){
     445
     446                        // exponential time decay & wetness
     447
     448                        // change in albedo with time:
     449                        //   (d_a) = (a - a_old)/(t0)
     450                        // where: t0 = timescale for albedo decay
     451
     452                        dt = dt / dts;    // convert from [s] to [d]
     453
     454                        // initialize variables
     455                        IssmDouble* t0=xNew<IssmDouble>(m);
     456                        IssmDouble* T=xNew<IssmDouble>(m);
     457                        IssmDouble* t0warm=xNew<IssmDouble>(m);
     458                        IssmDouble* d_a=xNew<IssmDouble>(m);
     459
     460                        // specify constants
     461                        // a_wet = 0.15;        // water albedo (0.15)
     462                        // a_new = aSnow        // new snow albedo (0.64 - 0.89)
     463                        // a_old = aIce;        // old snow/ice albedo (0.27-0.53)
     464                        // t0_wet = t0wet;      // time scale for wet snow (15-21.9) [d]
     465                        // t0_dry = t0dry;      // warm snow timescale [15] [d]
     466                        // K = 7                // time scale temperature coef. (7) [d]
     467                        // W0 = 300;            // 200 - 600 [mm]
     468                        const IssmDouble z_snow = 15.0;            // 16 - 32 [mm]
     469
     470                        // determine timescale for albedo decay
     471                        for(int i=0;i<m;i++)if(W[i]>0.0+Wtol)t0[i]=t0wet; // wet snow timescale
     472                        for(int i=0;i<m;i++)T[i]=TK[i] - CtoK; // change T from K to degC
     473                        for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
     474                        for(int i=0;i<m;i++)if(fabs(W[i])<Wtol && T[i]>=-10.0-Ttol)t0[i]= t0warm[i];
     475                        for(int i=0;i<m;i++)if(T[i]<-10.0-Ttol) t0[i] =  10.0 * K + t0dry; // 'cold' snow timescale
     476
     477                        // calculate new albedo
     478                        for(int i=0;i<m;i++)d_a[i] = (a[i] - aIce) / t0[i] * dt;           // change in albedo
     479                        for(int i=0;i<m;i++)a[i] -= d_a[i];                            // new albedo
     480
     481                        // modification of albedo due to thin layer of snow or solid
     482                        // condensation (deposition) at the surface surface
     483
     484                        // check if condensation occurs & if it is deposited in solid phase
     485                        if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]
     486
     487                        a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);
     488
     489                        //----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------
     490                        // modification of albedo due to thin layer of water on the surface
     491                        // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);
     492
     493                        /*Free ressources:*/
     494                        xDelete<IssmDouble>(t0);
     495                        xDelete<IssmDouble>(T);
     496                        xDelete<IssmDouble>(t0warm);
     497                        xDelete<IssmDouble>(d_a);
     498
     499                }
     500                else _error_("albedo method switch should range from 0 to 4!");
     501        }
     502
    503503        // Check for erroneous values
    504504        if (a[0] > 1) _printf_("albedo > 1.0\n");
     
    11761176
    11771177                                // adjust a, re, gdn & gsp
    1178                                 if(aIdx>0 | aIdx!=5 | (aIdx==5 & !(dIce - d[0]<Dtol)))a[0] = (aSnow * P + a[0] * mInit[0])/mass;
     1178                                if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
    11791179                                re[0] = (reNew * P + re[0] * mInit[0])/mass;
    11801180                                gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass;
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r22475 r22482  
    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 aValue, 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 adThresh, 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);
    2929void 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 thermo_scaling, IssmDouble dIce, int sid);
  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r22475 r22482  
    436436        SmbInitDensityScalingEnum,
    437437        SmbThermoDeltaTScalingEnum,
     438        SmbAdThreshEnum,
    438439        SmbTaEnum,
    439440        SmbVEnum,
  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r22475 r22482  
    438438                case SmbInitDensityScalingEnum : return "SmbInitDensityScaling";
    439439                case SmbThermoDeltaTScalingEnum : return "SmbThermoDeltaTScaling";
     440                case SmbAdThreshEnum : return "SmbAdThresh";
    440441                case SmbTaEnum : return "SmbTa";
    441442                case SmbVEnum : return "SmbV";
  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r22475 r22482  
    447447              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    448448              else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum;
     449              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
    449450              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    450451              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
     
    505506              else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
    506507              else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
    507               else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
     511              if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
     512              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    512513              else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
    513514              else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum;
     
    628629              else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    629630              else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    630               else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
     634              if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
     635              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
    635636              else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
    636637              else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
     
    751752              else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
    752753              else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
    753               else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
     757              if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
     758              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    758759              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    759760              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
     
    874875              else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum;
    875876              else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum;
    876               else if (strcmp(name,"LoveG0")==0) return LoveG0Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"LoveR0")==0) return LoveR0Enum;
     880              if (strcmp(name,"LoveG0")==0) return LoveG0Enum;
     881              else if (strcmp(name,"LoveR0")==0) return LoveR0Enum;
    881882              else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
    882883              else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
     
    997998              else if (strcmp(name,"Penta")==0) return PentaEnum;
    998999              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    999               else if (strcmp(name,"Vertex")==0) return VertexEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
     1003              if (strcmp(name,"Vertex")==0) return VertexEnum;
     1004              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    10041005              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    10051006              else if (strcmp(name,"Option")==0) return OptionEnum;
     
    11201121              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    11211122              else if (strcmp(name,"MINI")==0) return MINIEnum;
    1122               else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
     1126              if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
     1127              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    11271128              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    11281129              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
  • issm/trunk-jpl/src/m/classes/SMBgemb.m

    r22475 r22482  
    6060                                          % 3: density and cloud amount [Greuell & Konzelmann, 1994]
    6161                                          % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
    62                                           % 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only
    6362
    6463                swIdx  = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
     
    8988                t0dry = NaN; % warm snow timescale (30)
    9089                K     = NaN; % time scale temperature coef. (7)
     90                adThresh = NaN; %Apply aIdx method to all areas with densities below this value,
     91                                %or else apply direct input value from aValue, allowing albedo to be altered.
     92                                                         %Default value is rho water (1023 kg m-3).
    9193
    9294                %densities:
     
    125127                        self.pAir=project3d(md,'vector',self.pAir,'type','node');
    126128
    127                         if (aIdx == 0 | aIdx == 5) & ~isnan(self.aValue)
     129                        if (aIdx == 0) & ~isnan(self.aValue)
    128130                                self.aValue=project3d(md,'vector',self.aValue,'type','node');
    129131                        end
     
    169171                self.t0dry = 30;
    170172                self.K = 7;
     173                self.adThresh = 1023;
    171174
    172175                self.teValue = ones(mesh.numberofelements,1);
     
    212215                        md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    213216
    214                         md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4,5]);
     217                        md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]);
    215218                        md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
    216219                        md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);
     
    223226                        md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
    224227                        md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
     228                        md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0);
    225229
    226230                        switch self.aIdx,
     
    278282                        fielddisplay(self,'ThermoDeltaTScaling',{'scaling factor to multiply the thermal diffusion timestep (delta t)'});
    279283                        fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)');
     284                        fielddisplay(self,'adThresh',{'Apply aIdx method to all areas with densities below this value,','or else apply direct input value from aValue, allowing albedo to be altered.'});
    280285                        fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
    281286                                '0: direct input from aValue parameter',...
     
    283288                                '2: effective grain radius [Brun et al., 2009]',...
    284289                                '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
    285                                 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]',...
    286                                 '5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only'});
     290                                '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'})
    287291
    288292                        fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)');
     
    373377                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double');
    374378                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
     379                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double');
    375380
    376381                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r22475 r22482  
    9595                t0dry = float('NaN')    # warm snow timescale (30)
    9696                K     = float('NaN')    # time scale temperature coef. (7)
     97                adThresh = float('NaN') # Apply aIdx method to all areas with densities below this value,
     98                                        # or else apply direct input value from aValue, allowing albedo to be altered.
     99                                                                                # Default value is rho water (1023 kg m-3).
    97100
    98101                #densities:
     
    145148                string = "%s\n%s"%(string,fielddisplay(self,'ThermoDeltaTScaling','scaling factor to multiply the thermal diffusion timestep (delta t)'))
    146149                string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'))
     150                string = "%s\n%s"%(string,fielddisplay(self,'adThresh','Apply aIdx method to all areas with densities below this value,','or else apply direct input value from aValue, allowing albedo to be altered.'))
    147151                string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)',
    148152                                 '0: direct input from aValue parameter',
     
    150154                                                '2: effective grain radius [Brun et al., 2009]',
    151155                                                '3: density and cloud amount [Greuell & Konzelmann, 1994]',
    152                                                 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]',
    153                                                 '5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only']))
     156                                                '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
    154157
    155158                string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'))
     
    203206                self.pAir = project3d(md,'vector',self.pAir,'type','node')
    204207
    205                 if (aIdx == 0 or aIdx == 5) and np.isnan(self.aValue):
     208                if (aIdx == 0) and np.isnan(self.aValue):
    206209                        self.aValue=project3d(md,'vector',self.aValue,'type','node');
    207210                if np.isnan(self.teValue):
     
    245248                self.t0dry = 30
    246249                self.K = 7
     250                self.adThresh = 1023
    247251
    248252                self.teValue = np.ones((mesh.numberofelements,));
     
    288292                md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    289293
    290                 md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4,5])
     294                md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4])
    291295                md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1])
    292296                md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5])
     
    299303                md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1)
    300304                md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1)
     305                md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0)
    301306
    302307                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)):
     
    368373                WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double')
    369374                WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double')
     375                WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double');
    370376
    371377                WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
Note: See TracChangeset for help on using the changeset viewer.