Changeset 24691


Ignore:
Timestamp:
04/03/20 17:56:25 (5 years ago)
Author:
schlegel
Message:

CHG: add pore closeoff and ice albedo to GEMB and Archives

Location:
issm/trunk-jpl
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r24683 r24691  
    36273627   IssmDouble sumR=0.0;
    36283628   IssmDouble sumM=0.0;
     3629        IssmDouble sumMsurf=0.0;
    36293630   IssmDouble sumEC=0.0;
    36303631   IssmDouble sumP=0.0;
    … …  
    36593660        IssmDouble  T_bottom = 0.0;
    36603661        IssmDouble  M = 0.0;
     3662        IssmDouble  Msurf = 0.0;
    36613663        IssmDouble  R = 0.0;
    36623664        IssmDouble  mAdd = 0.0;
    … …  
    38153817
    38163818        // initialize cumulative variables
    3817         sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
     3819        sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; sumMsurf = 0;
    38183820
    38193821        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
    … …  
    38373839                Input2 *MassAdd_input       = this->GetInput2(SmbMAddEnum);  _assert_(MassAdd_input);
    38383840                Input2 *InitMass_input      = this->GetInput2(SmbMInitnum);  _assert_(InitMass_input);
     3841                Input2 *sumMsurf_input         = this->GetInput2(SmbMSurfEnum);  _assert_(sumMsurf_input);
    38393842
    38403843                ULW_input->GetInputAverage(&meanULW);
    … …  
    38513854                sumM_input->GetInputAverage(&sumM);
    38523855                sumM=sumM*dt*rho_ice;
     3856                sumMsurf_input->GetInputAverage(&sumMsurf);
     3857      sumMsurf=sumMsurf*dt*rho_ice;
    38533858                sumR_input->GetInputAverage(&sumR);
    38543859                sumR=sumR*dt*rho_ice;
    … …  
    38863891
    38873892        /*Snow, firn and ice albedo:*/
    3888         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());
     3893        if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    38893894
    38903895        /*Distribution of absorbed short wave radation with depth:*/
    … …  
    39063911        /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    39073912         * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    3908         if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
     3913        if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
    39093914
    39103915        /*Allow non-melt densification and determine compaction [m]*/
    … …  
    39463951        sumMassAdd = mAdd + sumMassAdd;
    39473952        sumM = M + sumM;
     3953        sumMsurf = Msurf + sumMsurf;
    39483954        sumR = R + sumR;
    39493955        sumW = cellsum(W,m);
    … …  
    39994005        this->SetElementInput(SmbMInitnum,initMass);
    40004006        this->SetElementInput(SmbMAddEnum,sumMassAdd/dt);
     4007        this->SetElementInput(SmbMSurfEnum,sumMsurf/dt/rho_ice);
    40014008        this->SetElementInput(SmbWAddEnum,sumW/dt);
    40024009        this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters

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

    r24689 r24691  
    400400
    401401}  /*}}}*/
    402 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) { /*{{{*/
     402void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
    403403
    404404        //// Calculates Snow, firn and ice albedo as a function of:
    … …  
    463463        //some constants:
    464464        const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
     465        const IssmDouble dPHC = 830;  //Pore closeoff density
     466        const IssmDouble ai_max = 0.58;  //maximum ice albedo, from Lefebre,2003
     467        const IssmDouble ai_min = aIce;  //minimum ice albedo
     468        const IssmDouble as_min = 0.65;  //minimum snow albedo, from Alexander 2014
    465469
    466470        if(aIdx==0 || (adThresh - d[0])<Dtol){
    … …  
    481485                        // Spectral fractions  (Lefebre et al., 2003)
    482486                        // [0.3-0.8um 0.8-1.5um 1.5-2.8um]
    483 
    484                         IssmDouble sF[3] = {0.606, 0.301, 0.093};
    485 
    486                         // convert effective radius to grain size in meters
    487                         IssmDouble gsz = (re[0] * 2.0) / 1000.0;
    488 
    489                         // spectral range:
    490                         // 0.3 - 0.8um
    491                         IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5));
    492                         // 0.8 - 1.5um
    493                         IssmDouble a1 = max(0., 0.95 - 15.4 *pow(gsz,0.5));
    494                         // 1.5 - 2.8um
    495                         IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
    496 
    497                         // broadband surface albedo
    498                         a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
     487                        if (d[0]<dPHC-Dtol){
     488
     489                                IssmDouble sF[3] = {0.606, 0.301, 0.093};
     490
     491                                // convert effective radius to grain size in meters
     492                                IssmDouble gsz = (re[0] * 2.0) / 1000.0;
     493
     494                                // spectral range:
     495                                // 0.3 - 0.8um
     496                                IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5));
     497                                // 0.8 - 1.5um
     498                                IssmDouble a1 = max(0., 0.95 - 15.4 *pow(gsz,0.5));
     499                                // 1.5 - 2.8um
     500                                IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
     501
     502                                // broadband surface albedo
     503                                a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
     504
     505                                // In a layer < 10cm, account for mix of ice and snow,
     506                                // after P. Alexander et al., 2014
     507                                IssmDouble depthsnow=0.0;
     508                                IssmDouble aice=0.0;
     509                                int lice=0;
     510                                for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
     511                                        depthsnow=depthsnow+dz[l];
     512                                        lice=l+1;
     513                                }
     514                                if (depthsnow<=0.1+Dtol & lice<m & d[lice]>=dPHC-Dtol){
     515                                        aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce);
     516                                        a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1);
     517                                }
     518                        }
     519                        else if (d[0]<dIce-Dtol){ //For continuity of albedo in firn i.e. P. Alexander et al., 2014
     520
     521                                //ai=ai_max + (as_min - ai_max)*(dI-dIce)/(dC-dIce)
     522                                //dC is pore close off (830 kg m^-3)
     523                                //dI is density of the upper firn layer
     524
     525                                a[0] = ai_max + (as_min - ai_max)*(d[0]-dIce)/(dPHC-dIce);
     526
     527                        }
     528                        else{ //surface layer is density of ice
     529
     530                                //When density is > dIce (typically 910 kg m^-3, 920 is used by Alexander in MAR),
     531                                //ai=ai_min + (ai_max - ai_min)*e^(-1*(Msw(t)/K))
     532                                //K is a scale factor (set to 200 kg m^-2)
     533                                //Msw(t) is the time-dependent accumulated amount of excessive surface meltwater
     534                                //  before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly)
     535                                IssmDouble M = Msurf+W[0];
     536                                a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min);
     537
     538                        }
    499539                }
    500540                else if(aIdx==3){
    … …  
    817857                // calculated.  The estimated enegy balance & melt are significanly
    818858                // less when Ts is taken as the mean of the x top grid cells.
    819                 Ts = (T[0] + T[1])/2.0;
     859                if(m>1) Ts = (T[0] + T[1])/2.0;
     860                else Ts = T[0];
    820861                Ts = min(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
    821862
    … …  
    13291370        *pm=m;
    13301371} /*}}}*/
    1331 void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid){ /*{{{*/
     1372void melt(IssmDouble* pM, IssmDouble* pMs, 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){ /*{{{*/
    13321373
    13331374        //// MELT ROUTINE
    … …  
    13881429
    13891430        /*outputs:*/
     1431        IssmDouble  Msurf = 0.0;
    13901432        IssmDouble  mAdd = 0.0;
    13911433        IssmDouble  surplusE = 0.0;
    … …  
    14381480        // specify irreducible water content saturation [fraction]
    14391481        const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
     1482        const IssmDouble dPHC = 830;                     //Pore closeoff density
    14401483
    14411484        //// REFREEZE PORE WATER
    … …  
    15161559                        M[i] = min(Mmax, m[i]);  // melt
    15171560                }
     1561                Msurf = M[0];
    15181562                sumM = cellsum(M,n);                                                       // total melt [kg]
    15191563
    … …  
    15371581                }
    15381582
     1583                IssmDouble depthice=0.0;
    15391584                //// meltwater percolation
    15401585                for(int i=0;i<n;i++){
    … …  
    15421587                        IssmDouble inM = M[i]+ flxDn[i];
    15431588
     1589                        depthice=0.0;
     1590                        if (d[i] >= dPHC-Dtol){
     1591                                for(int l=i;(l<n & d[l]>=dPHC-Dtol);l++) depthice=depthice+dz[l];
     1592                        }
     1593
    15441594                        // break loop if there is no meltwater and if depth is > mw_depth
    15451595                        if (fabs(inM) < Wtol && i > X){
    15461596                                break;
    15471597                        }
    1548 
    15491598                        // if reaches impermeable ice layer all liquid water runs off (R)
    1550                         else if (d[i] >= dIce-Dtol){  // dPHC = pore hole close off [kg m-3]
     1599                        else if (d[i] >= dIce-Dtol || (d[i] >= dPHC-Dtol && depthice>0.1)){  // dPHC = pore hole close off [kg m-3]
    15511600                                // _printf_("ICE LAYER");
    15521601                                // no water freezes in this cell
    … …  
    19001949
    19011950        /*Assign output pointers:*/
     1951        *pMs=Msurf;
    19021952        *pM=sumM;
    19031953        *pR=Rsum;
    … …  
    20402090                                c1arth = 0.03 * H;
    20412091                                //ERA-5
    2042                                 //M0 = max(2.3999 - (0.2610 * log(C)),0.25);
    2043                                 //M1 = max(2.7469 - (0.3228 * log(C)),0.25);
     2092                                //M0 = max(2.3128 - (0.2480 * log(C)),0.25);
     2093                                //M1 = max(2.7950 - (0.3318 * log(C)),0.25);
    20442094                                //RACMO
    20452095                                M0 = max(1.6599 - (0.1724 * log(C)),0.25);
    … …  
    20602110                                c1arth = 0.03 * H;
    20612111                                // ERA5
    2062                                 //M0 = max(1.8920 - (0.1569 * log(C)),0.25);
    2063                                 //M1 = max(2.5662 - (0.2274 * log(C)),0.25);
     2112                                //M0 = max(1.8554 - (0.1316 * log(C)),0.25);
     2113                                //M1 = max(2.8901 - (0.3014 * log(C)),0.25);
    20642114                                // RACMO
    20652115                                M0 = max(1.6201 - (0.1450 * log(C)),0.25);

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

    r23846 r24691  
    2929IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
    3030void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
    31 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);
     31void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
    3232void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
    3333void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, 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);
    3434void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
    35 void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid);
     35void melt(IssmDouble* pM, IssmDouble* pMs, 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);
    3636void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
    3737void 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);

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24683 r24691  
    739739syn keyword cConstant SmbMeltEnum
    740740syn keyword cConstant SmbMonthlytemperaturesEnum
     741syn keyword cConstant SmbMSurfEnum
    741742syn keyword cConstant SmbNetLWEnum
    742743syn keyword cConstant SmbNetSWEnum

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

    r24683 r24691  
    737737        SmbMInitnum,
    738738        SmbMonthlytemperaturesEnum,
     739        SmbMSurfEnum,
    739740        SmbNetLWEnum,
    740741        SmbNetSWEnum,

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

    r24683 r24691  
    741741                case SmbMeltEnum : return "SmbMelt";
    742742                case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures";
     743                case SmbMSurfEnum : return "SmbMSurf";
    743744                case SmbNetLWEnum : return "SmbNetLW";
    744745                case SmbNetSWEnum : return "SmbNetSW";

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

    r24683 r24691  
    759759              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    760760              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
     761              else if (strcmp(name,"SmbMSurf")==0) return SmbMSurfEnum;
    761762              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
    762763              else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    875876              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    876               else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
     880              if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
     881              else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    881882              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    882883              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    … …  
    997998              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    998999              else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
    999               else if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     1003              if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
     1004              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    10041005              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    10051006              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    … …  
    11201121              else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
    11211122              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
    1122               else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     1126              if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
     1127              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    11271128              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    11281129              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    … …  
    12431244              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
    12441245              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
    1245               else if (strcmp(name,"Regular")==0) return RegularEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
     1249              if (strcmp(name,"Regular")==0) return RegularEnum;
     1250              else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    12501251              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    12511252              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
Note: See TracChangeset for help on using the changeset viewer.