Changeset 22482
- Timestamp:
- 02/27/18 19:33:05 (7 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp
r22475 r22482 179 179 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum)); 180 180 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum)); 181 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum)); 181 182 break; 182 183 case SMBpddEnum: -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r22475 r22482 2606 2606 IssmDouble init_scaling=0.0; 2607 2607 IssmDouble thermo_scaling=1.0; 2608 IssmDouble adThresh=1023.0; 2608 2609 2609 2610 /*}}}*/ … … 2668 2669 parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); 2669 2670 parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum); 2671 parameters->FindParam(&adThresh,SmbAdThreshEnum); 2670 2672 2671 2673 /*}}}*/ … … 2694 2696 Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); 2695 2697 2696 if (aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);2697 2698 2698 zTop_input->GetInputValue(&zTop,gauss); 2699 2699 dzTop_input->GetInputValue(&dzTop,gauss); … … 2707 2707 Vz_input->GetInputValue(&Vz,gauss); 2708 2708 teValue_input->GetInputValue(&teValue,gauss); 2709 aValue_input->GetInputValue(&aValue,gauss); 2709 2710 isinitialized_input->GetInputValue(&isinitialized); 2710 2711 /*}}}*/ … … 2833 2834 pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa] 2834 2835 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] 2836 2837 //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); 2837 2838 /*}}}*/ … … 2841 2842 2842 2843 /*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()); 2844 2845 2845 2846 -
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
r22475 r22482 337 337 338 338 } /*}}}*/ 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) { /*{{{*/339 void 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) { /*{{{*/ 340 340 341 341 //// Calculates Snow, firn and ice albedo as a function of: … … 345 345 // 3 : density and cloud amount (Greuell & Konzelmann, 1994) 346 346 // 4 : exponential time decay & wetness (Bougamont & Bamber, 2005) 347 // 5 : ingest MODIS mode, direct input from aValue parameter applied to surface ice only348 347 349 348 //// Inputs 350 349 // 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). 354 358 355 359 // Methods 1 & 2 … … 397 401 const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] 398 402 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 503 503 // Check for erroneous values 504 504 if (a[0] > 1) _printf_("albedo > 1.0\n"); … … 1176 1176 1177 1177 // 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; 1179 1179 re[0] = (reNew * P + re[0] * mInit[0])/mass; 1180 1180 gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass; -
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
r22475 r22482 25 25 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); 26 26 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); 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);27 void 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); 28 28 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); 29 29 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 thermo_scaling, IssmDouble dIce, int sid); -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r22475 r22482 436 436 SmbInitDensityScalingEnum, 437 437 SmbThermoDeltaTScalingEnum, 438 SmbAdThreshEnum, 438 439 SmbTaEnum, 439 440 SmbVEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r22475 r22482 438 438 case SmbInitDensityScalingEnum : return "SmbInitDensityScaling"; 439 439 case SmbThermoDeltaTScalingEnum : return "SmbThermoDeltaTScaling"; 440 case SmbAdThreshEnum : return "SmbAdThresh"; 440 441 case SmbTaEnum : return "SmbTa"; 441 442 case SmbVEnum : return "SmbV"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r22475 r22482 447 447 else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; 448 448 else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum; 449 else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; 449 450 else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; 450 451 else if (strcmp(name,"SmbV")==0) return SmbVEnum; … … 505 506 else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum; 506 507 else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; 507 else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;508 508 else stage=5; 509 509 } 510 510 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; 512 513 else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; 513 514 else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum; … … 628 629 else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; 629 630 else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; 630 else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;631 631 else stage=6; 632 632 } 633 633 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; 635 636 else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; 636 637 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; … … 751 752 else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; 752 753 else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; 753 else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;754 754 else stage=7; 755 755 } 756 756 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; 758 759 else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; 759 760 else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; … … 874 875 else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; 875 876 else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; 876 else if (strcmp(name,"LoveG0")==0) return LoveG0Enum;877 877 else stage=8; 878 878 } 879 879 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; 881 882 else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; 882 883 else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; … … 997 998 else if (strcmp(name,"Penta")==0) return PentaEnum; 998 999 else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; 999 else if (strcmp(name,"Vertex")==0) return VertexEnum;1000 1000 else stage=9; 1001 1001 } 1002 1002 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; 1004 1005 else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; 1005 1006 else if (strcmp(name,"Option")==0) return OptionEnum; … … 1120 1121 else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; 1121 1122 else if (strcmp(name,"MINI")==0) return MINIEnum; 1122 else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;1123 1123 else stage=10; 1124 1124 } 1125 1125 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; 1127 1128 else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; 1128 1129 else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; -
issm/trunk-jpl/src/m/classes/SMBgemb.m
r22475 r22482 60 60 % 3: density and cloud amount [Greuell & Konzelmann, 1994] 61 61 % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] 62 % 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only63 62 64 63 swIdx = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) … … 89 88 t0dry = NaN; % warm snow timescale (30) 90 89 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). 91 93 92 94 %densities: … … 125 127 self.pAir=project3d(md,'vector',self.pAir,'type','node'); 126 128 127 if (aIdx == 0 | aIdx == 5) & ~isnan(self.aValue)129 if (aIdx == 0) & ~isnan(self.aValue) 128 130 self.aValue=project3d(md,'vector',self.aValue,'type','node'); 129 131 end … … 169 171 self.t0dry = 30; 170 172 self.K = 7; 173 self.adThresh = 1023; 171 174 172 175 self.teValue = ones(mesh.numberofelements,1); … … 212 215 md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); 213 216 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]); 215 218 md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]); 216 219 md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]); … … 223 226 md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1); 224 227 md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'>=',0,'<=',1); 228 md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0); 225 229 226 230 switch self.aIdx, … … 278 282 fielddisplay(self,'ThermoDeltaTScaling',{'scaling factor to multiply the thermal diffusion timestep (delta t)'}); 279 283 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.'}); 280 285 fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',... 281 286 '0: direct input from aValue parameter',... … … 283 288 '2: effective grain radius [Brun et al., 2009]',... 284 289 '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]'}) 287 291 288 292 fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'); … … 373 377 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double'); 374 378 WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double'); 379 WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double'); 375 380 376 381 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 95 95 t0dry = float('NaN') # warm snow timescale (30) 96 96 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). 97 100 98 101 #densities: … … 145 148 string = "%s\n%s"%(string,fielddisplay(self,'ThermoDeltaTScaling','scaling factor to multiply the thermal diffusion timestep (delta t)')) 146 149 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.')) 147 151 string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)', 148 152 '0: direct input from aValue parameter', … … 150 154 '2: effective grain radius [Brun et al., 2009]', 151 155 '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]'])) 154 157 155 158 string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)')) … … 203 206 self.pAir = project3d(md,'vector',self.pAir,'type','node') 204 207 205 if (aIdx == 0 or aIdx == 5) and np.isnan(self.aValue):208 if (aIdx == 0) and np.isnan(self.aValue): 206 209 self.aValue=project3d(md,'vector',self.aValue,'type','node'); 207 210 if np.isnan(self.teValue): … … 245 248 self.t0dry = 30 246 249 self.K = 7 250 self.adThresh = 1023 247 251 248 252 self.teValue = np.ones((mesh.numberofelements,)); … … 288 292 md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); 289 293 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]) 291 295 md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]) 292 296 md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]) … … 299 303 md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1) 300 304 md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1) 305 md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0) 301 306 302 307 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)): … … 368 373 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double') 369 374 WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double') 375 WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double'); 370 376 371 377 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.