Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24690) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24691) @@ -3627,4 +3627,5 @@ IssmDouble sumR=0.0; IssmDouble sumM=0.0; + IssmDouble sumMsurf=0.0; IssmDouble sumEC=0.0; IssmDouble sumP=0.0; @@ -3659,4 +3660,5 @@ IssmDouble T_bottom = 0.0; IssmDouble M = 0.0; + IssmDouble Msurf = 0.0; IssmDouble R = 0.0; IssmDouble mAdd = 0.0; @@ -3815,5 +3817,5 @@ // initialize cumulative variables - sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; + sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; sumMsurf = 0; //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. @@ -3837,4 +3839,5 @@ Input2 *MassAdd_input = this->GetInput2(SmbMAddEnum); _assert_(MassAdd_input); Input2 *InitMass_input = this->GetInput2(SmbMInitnum); _assert_(InitMass_input); + Input2 *sumMsurf_input = this->GetInput2(SmbMSurfEnum); _assert_(sumMsurf_input); ULW_input->GetInputAverage(&meanULW); @@ -3851,4 +3854,6 @@ sumM_input->GetInputAverage(&sumM); sumM=sumM*dt*rho_ice; + sumMsurf_input->GetInputAverage(&sumMsurf); + sumMsurf=sumMsurf*dt*rho_ice; sumR_input->GetInputAverage(&sumR); sumR=sumR*dt*rho_ice; @@ -3886,5 +3891,5 @@ /*Snow, firn and ice albedo:*/ - 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()); + 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()); /*Distribution of absorbed short wave radation with depth:*/ @@ -3906,5 +3911,5 @@ /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ - 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()); + 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()); /*Allow non-melt densification and determine compaction [m]*/ @@ -3946,4 +3951,5 @@ sumMassAdd = mAdd + sumMassAdd; sumM = M + sumM; + sumMsurf = Msurf + sumMsurf; sumR = R + sumR; sumW = cellsum(W,m); @@ -3999,4 +4005,5 @@ this->SetElementInput(SmbMInitnum,initMass); this->SetElementInput(SmbMAddEnum,sumMassAdd/dt); + this->SetElementInput(SmbMSurfEnum,sumMsurf/dt/rho_ice); this->SetElementInput(SmbWAddEnum,sumW/dt); this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24690) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24691) @@ -400,5 +400,5 @@ } /*}}}*/ -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) { /*{{{*/ +void 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) { /*{{{*/ //// Calculates Snow, firn and ice albedo as a function of: @@ -463,4 +463,8 @@ //some constants: const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] + const IssmDouble dPHC = 830; //Pore closeoff density + const IssmDouble ai_max = 0.58; //maximum ice albedo, from Lefebre,2003 + const IssmDouble ai_min = aIce; //minimum ice albedo + const IssmDouble as_min = 0.65; //minimum snow albedo, from Alexander 2014 if(aIdx==0 || (adThresh - 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); + + } } else if(aIdx==3){ @@ -817,5 +857,6 @@ // calculated. The estimated enegy balance & melt are significanly // less when Ts is taken as the mean of the x top grid cells. - Ts = (T[0] + T[1])/2.0; + if(m>1) Ts = (T[0] + T[1])/2.0; + else Ts = T[0]; Ts = min(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) @@ -1329,5 +1370,5 @@ *pm=m; } /*}}}*/ -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){ /*{{{*/ +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 dIce, int sid){ /*{{{*/ //// MELT ROUTINE @@ -1388,4 +1429,5 @@ /*outputs:*/ + IssmDouble Msurf = 0.0; IssmDouble mAdd = 0.0; IssmDouble surplusE = 0.0; @@ -1438,4 +1480,5 @@ // specify irreducible water content saturation [fraction] const IssmDouble Swi = 0.07; // assumed constant after Colbeck, 1974 + const IssmDouble dPHC = 830; //Pore closeoff density //// REFREEZE PORE WATER @@ -1516,4 +1559,5 @@ M[i] = min(Mmax, m[i]); // melt } + Msurf = M[0]; sumM = cellsum(M,n); // total melt [kg] @@ -1537,4 +1581,5 @@ } + IssmDouble depthice=0.0; //// meltwater percolation for(int i=0;i= dPHC-Dtol){ + for(int l=i;(l=dPHC-Dtol);l++) depthice=depthice+dz[l]; + } + // break loop if there is no meltwater and if depth is > mw_depth if (fabs(inM) < Wtol && i > X){ break; } - // if reaches impermeable ice layer all liquid water runs off (R) - else if (d[i] >= dIce-Dtol){ // dPHC = pore hole close off [kg m-3] + else if (d[i] >= dIce-Dtol || (d[i] >= dPHC-Dtol && depthice>0.1)){ // dPHC = pore hole close off [kg m-3] // _printf_("ICE LAYER"); // no water freezes in this cell @@ -1900,4 +1949,5 @@ /*Assign output pointers:*/ + *pMs=Msurf; *pM=sumM; *pR=Rsum; @@ -2040,6 +2090,6 @@ c1arth = 0.03 * H; //ERA-5 - //M0 = max(2.3999 - (0.2610 * log(C)),0.25); - //M1 = max(2.7469 - (0.3228 * log(C)),0.25); + //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); @@ -2060,6 +2110,6 @@ c1arth = 0.03 * H; // ERA5 - //M0 = max(1.8920 - (0.1569 * log(C)),0.25); - //M1 = max(2.5662 - (0.2274 * log(C)),0.25); + //M0 = max(1.8554 - (0.1316 * log(C)),0.25); + //M1 = max(2.8901 - (0.3014 * log(C)),0.25); // RACMO M0 = max(1.6201 - (0.1450 * log(C)),0.25); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 24690) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 24691) @@ -29,9 +29,9 @@ IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); 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); -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); +void 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); void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); void 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); 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* 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); +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 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 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); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24690) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24691) @@ -739,4 +739,5 @@ syn keyword cConstant SmbMeltEnum syn keyword cConstant SmbMonthlytemperaturesEnum +syn keyword cConstant SmbMSurfEnum syn keyword cConstant SmbNetLWEnum syn keyword cConstant SmbNetSWEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24690) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24691) @@ -737,4 +737,5 @@ SmbMInitnum, SmbMonthlytemperaturesEnum, + SmbMSurfEnum, SmbNetLWEnum, SmbNetSWEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24690) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24691) @@ -741,4 +741,5 @@ case SmbMeltEnum : return "SmbMelt"; case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures"; + case SmbMSurfEnum : return "SmbMSurf"; case SmbNetLWEnum : return "SmbNetLW"; case SmbNetSWEnum : return "SmbNetSW"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24690) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24691) @@ -759,4 +759,5 @@ else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum; else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum; + else if (strcmp(name,"SmbMSurf")==0) return SmbMSurfEnum; else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum; else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; - else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; + if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; + else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; - else if (strcmp(name,"IntInput2")==0) return IntInput2Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + if (strcmp(name,"IntInput2")==0) return IntInput2Enum; + else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; - else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Incremental")==0) return IncrementalEnum; + if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; + else if (strcmp(name,"Incremental")==0) return IncrementalEnum; else if (strcmp(name,"Indexed")==0) return IndexedEnum; else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; - else if (strcmp(name,"Regular")==0) return RegularEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; + if (strcmp(name,"Regular")==0) return RegularEnum; + else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;