Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25395) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25396) @@ -3915,5 +3915,5 @@ /*Allow non-melt densification and determine compaction [m]*/ - if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); + if(isdensification)densification(&d,&dz, T, re, denIdx, aIdx, swIdx, adThresh, C, smb_dt, Tmean,rho_ice,m,this->Sid()); /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 25395) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 25396) @@ -487,55 +487,33 @@ // Spectral fractions (Lefebre et al., 2003) // [0.3-0.8um 0.8-1.5um 1.5-2.8um] - if (d[0]=dPHC-Dtol){ - aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce); - a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1); - } - } - else if (d[0] dIce (typically 910 kg m^-3, 920 is used by Alexander in MAR), - //ai=ai_min + (ai_max - ai_min)*e^(-1*(Msw(t)/K)) - //K is a scale factor (set to 200 kg m^-2) - //Msw(t) is the time-dependent accumulated amount of excessive surface meltwater - // before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly) - IssmDouble M = Msurf+W[0]; - a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min); - + + IssmDouble sF[3] = {0.606, 0.301, 0.093}; + + // convert effective radius to grain size in meters + IssmDouble gsz = (re[0] * 2.0) / 1000.0; + + // spectral range: + // 0.3 - 0.8um + IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5)); + // 0.8 - 1.5um + IssmDouble a1 = max(0., 0.95 - 15.4 *pow(gsz,0.5)); + // 1.5 - 2.8um + IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); + + // broadband surface albedo + a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2; + + // In a layer < 10cm, account for mix of ice and snow, + // after P. Alexander et al., 2014 + IssmDouble depthsnow=0.0; + IssmDouble aice=0.0; + int lice=0; + for(int l=0;(l=dPHC-Dtol){ + aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce); + a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1); } } @@ -604,4 +582,44 @@ } else _error_("albedo method switch should range from 0 to 4!"); + + //If we do not have fresh snow + if (aIdx<3 && aIdx>0){ + // In a layer < 10cm, account for mix of ice and snow, + // after P. Alexander et al., 2014 + IssmDouble depthsnow=0.0; + IssmDouble aice=0.0; + int lice=0; + for(int l=0;(l=dPHC-Dtol){ + aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce); + a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1); + } + + if (d[0]>=dPHC-Dtol){ + if (d[0] dIce (typically 910 kg m^-3, 920 is used by Alexander in MAR), + //ai=ai_min + (ai_max - ai_min)*e^(-1*(Msw(t)/K)) + //K is a scale factor (set to 200 kg m^-2) + //Msw(t) is the time-dependent accumulated amount of excessive surface meltwater + // before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly) + IssmDouble M = Msurf+W[0]; + a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min); + + } + } + } } @@ -872,9 +890,11 @@ // calculate the Bulk Richardson Number (Ri) Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0)); + //Ri = (2.0*9.81*(Ta - Ts)) / ((Tz - z0)*(Ta + Ts)* pow(V/(Vz - z0),2.0)); // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' - // do not allow Ri to exceed 0.19 - Ri = min(Ri, 0.19); + // do not allow Ri to exceed 0.16 + // Ri = min(Ri, 0.19); + Ri = min(Ri, 0.16); //Ohmura, 1982 // calculate momentum 'coefM' stability factor @@ -929,5 +949,7 @@ // upward longwave contribution - ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here + //ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here + IssmDouble deltatest=0; + ulw = - (SB * pow(Ts,4.0)* teValue - deltatest) * dt; // - deltatest here ulwrf = ulwrf - ulw/dt0; @@ -996,6 +1018,6 @@ // swIdx = 1 : absorbed SW is distributed with depth as a function of: - // 1 : snow density (taken from Bassford, 2004) - // 2 : grain size in 3 spectral bands (Brun et al., 1992) + // default : snow density (taken from Bassford, 2004) + // if aIdx=2 : grain size in 3 spectral bands (Brun et al., 1992) // Inputs @@ -1254,4 +1276,5 @@ case 1: // Density of Antarctica snow dSnow = 350; + //dSnow = 360; //FirnMICE Lundin et al., 2017 break; @@ -1930,4 +1953,5 @@ xDelete(Zcum); xDelete(dzMin2); + xDelete(dzMax2); /*Assign output pointers:*/ @@ -1949,5 +1973,5 @@ } /*}}}*/ -void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/ +void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/ //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPACTION IS COMPENSATED FOR BY TRACES OF SNOW???] @@ -2072,14 +2096,38 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - //ERA-5 + //ERA-5 old //M0 = max(2.3128 - (0.2480 * log(C)),0.25); //M1 = max(2.7950 - (0.3318 * log(C)),0.25); - //RACMO - M0 = max(1.6599 - (0.1724 * log(C)),0.25); - M1 = max(2.0102 - (0.2458 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0 + if (aIdx==1 && swIdx==0){ + if (abs(adThresh - 820) < Dtol){ + // ERA5 new aIdx=1, swIdx=0, MODIS 820 + M0 = max(1.8785 - (0.1811 * log(C)),0.25); + M1 = max(2.0005 - (0.2346 * log(C)),0.25); + } + else{ + // ERA5 new aIdx=1, swIdx=0 + //M0 = max(1.8785 - (0.1811 * log(C)),0.25); + //M1 = max(2.0005 - (0.2346 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0, bare ice + M0 = max(1.8785 - (0.1811 * log(C)),0.25); + M1 = max(2.0005 - (0.2346 * log(C)),0.25); + } + } + // ERA5 new aIdx=2, swIdx=1 + else if (aIdx<3 && swIdx>0){ + M0 = max(2.2191 - (0.2301 * log(C)),0.25); + M1 = max(2.2917 - (0.2710 * log(C)),0.25); + } + // ERA5 new aIdx=2, swIdx=0 + //else if (aIdx==2){ + //} //From Ligtenberg //H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81); //M0 = max(1.435 - (0.151 * log(C)),0.25); //M1 = max(2.366 - (0.293 * log(C)),0.25); + //RACMO + M0 = max(1.6599 - (0.1724 * log(C)),0.25); + M1 = max(2.0102 - (0.2458 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; @@ -2092,13 +2140,37 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - // ERA5 + // ERA5 old //M0 = max(1.8554 - (0.1316 * log(C)),0.25); //M1 = max(2.8901 - (0.3014 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0 + if (aIdx==1 && swIdx==0){ + if (abs(adThresh - 820) < Dtol){ + // ERA5 new aIdx=1, swIdx=0, MODIS 820 + M0 = max(1.4174 - (0.1037 * log(C)),0.25); + M1 = max(2.2010 - (0.2460 * log(C)),0.25); + } + else{ + // ERA5 new aIdx=1, swIdx=0 + //M0 = max(1.4574 - (0.1123 * log(C)),0.25); + //M1 = max(2.0238 - (0.2070 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0, bare ice + M0 = max(1.4318 - (0.1055 * log(C)),0.25); + M1 = max(2.0453 - (0.2137 * log(C)),0.25); + } + } + // ERA5 new aIdx=2, swIdx=1 + else if (aIdx<3 && swIdx>0){ + M0 = max(1.7834 - (0.1409 * log(C)),0.25); + M1 = max(1.9260 - (0.1527 * log(C)),0.25); + } + // ERA5 new aIdx=2, swIdx=0 + //else if (aIdx==2){ + //} + // From Kuipers Munneke + //M0 = max(1.042 - (0.0916 * log(C)),0.25); + //M1 = max(1.734 - (0.2039 * log(C)),0.25); // RACMO M0 = max(1.6201 - (0.1450 * log(C)),0.25); M1 = max(2.5577 - (0.2899 * log(C)),0.25); - // From Kuipers Munneke - //M0 = max(1.042 - (0.0916 * log(C)),0.25); - //M1 = max(1.734 - (0.2039 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; @@ -2190,9 +2262,11 @@ // calculate the Bulk Richardson Number (Ri) Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0)); + //Ri = (2.0*9.81*(Ta - Ts)) / ((Tz - z0)*(Ta + Ts)* pow(V/(Vz - z0),2.0)); // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' - // do not allow Ri to exceed 0.19 - Ri = min(Ri, 0.19); + // do not allow Ri to exceed 0.16 + // Ri = min(Ri, 0.19); + Ri = min(Ri, 0.16); //Ohmura, 1982 // calculate momentum 'coefM' stability factor Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 25395) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 25396) @@ -34,5 +34,5 @@ void 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); void 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 zY, IssmDouble dIce, int sid); -void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); +void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid); #endif /* _SurfaceMassBalancex_H*/ Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 25395) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 25396) @@ -60,9 +60,9 @@ % 0: direct input from aValue parameter % 1: effective grain radius [Gardner & Sharp, 2009] - % 2: effective grain radius [Brun et al., 2009] + % 2: effective grain radius [Brun et al., 2009], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992) % 3: density and cloud amount [Greuell & Konzelmann, 1994] % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] - swIdx = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) + swIdx = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, with snow density (taken from Bassford, 2004)) denIdx = NaN; %densification model to use (default is 2): @@ -227,5 +227,5 @@ md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400,'size',size(self.Ta)); md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'size',size(self.Ta)); - md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100,'size',size(self.Ta)); + md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',200,'size',size(self.Ta)); md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1,'size',size(self.Ta)); @@ -330,5 +330,5 @@ '0: direct input from aValue parameter',... '1: effective grain radius [Gardner & Sharp, 2009]',... - '2: effective grain radius [Brun et al., 2009]',... + '2: effective grain radius [Brun et al., 2009], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)',... '3: density and cloud amount [Greuell & Konzelmann, 1994]',... '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}) @@ -364,5 +364,5 @@ end - fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]'); + fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1, with snow density (taken from Bassford, 2004)]'); fielddisplay(self,'denIdx',{'densification model to use (default is 2):',... '1 = emperical model of Herron and Langway (1980)',... Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25395) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25396) @@ -64,7 +64,27 @@ #settings: self.aIdx = np.nan # method for calculating albedo and subsurface absorption (default is 1) - self.swIdx = np.nan # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) + # 0: direct input from aValue parameter + # 1: effective grain radius [Gardner & Sharp, 2009] + # 2: effective grain radius [Brun et al., 2009], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992) + # 3: density and cloud amount [Greuell & Konzelmann, 1994] + # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + + self.swIdx = np.nan # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, with snow density (taken from Bassford, 2004)) self.denIdx = np.nan # densification model to use (default is 2): + # 1 = emperical model of Herron and Langway (1980) + # 2 = semi-emperical model of Anthern et al. (2010) + # 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010) + # 4 = DO NOT USE: emperical model of Li and Zwally (2004) + # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) + # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) + # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) + self.dsnowIdx = np.nan # model for fresh snow accumulation density (default is 1): + # 0 = Original GEMB value, 150 kg/m^3 + # 1 = Antarctica value of fresh snow density, 350 kg/m^3 + # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) + # 3 = Antarctica model of Kaspers et al. (2004) + # 4 = Greenland model of Kuipers Munneke et al. (2015) + self.zTop = np.nan # depth over which grid length is constant at the top of the snopack (default 10) [m] self.dzTop = np.nan # initial top vertical grid spacing (default .05) [m] @@ -155,5 +175,5 @@ '0: direct input from aValue parameter', '1: effective grain radius [Gardner & Sharp, 2009]', - '2: effective grain radius [Brun et al., 2009]', + '2: effective grain radius [Brun et al., 2009], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)', '3: density and cloud amount [Greuell & Konzelmann, 1994]', '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) @@ -186,5 +206,5 @@ string = "%s\n%s" % (string, fielddisplay(self, 'K', 'time scale temperature coef. (7) [d]')) - string = "%s\n%s" % (string, fielddisplay(self, 'swIdx', 'apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]')) + string = "%s\n%s" % (string, fielddisplay(self, 'swIdx', 'apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1, with snow density (taken from Bassford, 2004)]')) string = "%s\n%s" % (string, fielddisplay(self, 'denIdx', ['densification model to use (default is 2):', '1 = emperical model of Herron and Langway (1980)', @@ -305,5 +325,5 @@ md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta)) md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta)) - md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 100, 'size', np.shape(self.Ta)) + md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 200, 'size', np.shape(self.Ta)) md = checkfield(md, 'fieldname', 'smb.eAir', 'timeseries', 1, 'NaN', 1, 'Inf', 1, 'size', np.shape(self.Ta))