Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24188) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24189) @@ -4047,5 +4047,7 @@ /*Check bottom grid cell T is unchanged:*/ - if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); + if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){ + if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); + } /*Free ressources: */ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24188) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24189) @@ -248,61 +248,74 @@ if (gdn[i]>0.0+Gdntol){ - //_printf_("Dendritic dry snow metamorphism\n"); - - //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1 - if(dT[i]<5.0-Ttol){ - //determine coefficients - IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt; - IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt; - //new dentricity and sphericity for dT < 5 degC m-1 - gdn[i] += A; - gsp[i] += B; + if(W[i]<=0.0+Wtol){ + //_printf_("Dendritic dry snow metamorphism\n"); + //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1 + if(dT[i]<=5.0+Ttol){ + //determine coefficients + IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt; + IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt; + //new dentricity and sphericity for dT < 5 degC m-1 + gdn[i] += A; + gsp[i] += B; + } + else{ + // new dendricity and sphericity for dT >= 5 degC m-1 + + //determine coeff0icients + IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); + gdn[i] +=C; + gsp[i] +=C; + } } else{ - // new dendricity and sphericity for dT >= 5 degC m-1 - - //determine coefficients - IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); - gdn[i] +=C; - gsp[i] +=C; - } - - // wet snow metamorphism - if(W[i]>0.0+Wtol){ - - //_printf_("D}ritic wet snow metamorphism\n"); + // wet snow metamorphism + //_printf_("Dendritic wet snow metamorphism\n"); //determine coefficient IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt; - // new dendricity and sphericity for wet snow + // new dendricity for wet snow gdn[i] -= D; + // new sphericity for wet snow gsp[i] += D; } - - // dendricity and sphericity can not be > 1 or < 0 - if (gdn[i]<0.0+Wtol)gdn[i]=0.0; - if (gsp[i]<0.0+Wtol)gsp[i]=0.0; - if (gdn[i]>1.0-Wtol)gdn[i]=1.0; - if (gsp[i]>1.0-Wtol)gsp[i]=1.0; - - // determine new grain size (mm) - gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1; + // dendricity and sphericity can not be > 1 or < 0 + if (gdn[i]<0.0+Wtol)gdn[i]=0.0; + if (gsp[i]<0.0+Wtol)gsp[i]=0.0; + if (gdn[i]>1.0-Wtol)gdn[i]=1.0; + if (gsp[i]>1.0-Wtol)gsp[i]=1.0; + + // determine new grain size (mm) + gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1; } else{ - /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficinets + //When the state of "faceted crystals" (gsp==0) is fully reached, + // snow evolves towards depth hoar if the gradient is + // higher than 15 degC m-1 (Brun et al., 1992) + //When wet-snow grains (class 6) are submitted to a + // temperature gradient higher than degC m-1, their sphericity + // decreases according to Equations (4). When sphericity + // reaches 0, their size increases according to the functions + // determined by Marbouty. (Brun et al., 1992) + if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol && (dT[i]>15.0+Ttol || (T[i]>5.0+Ttol && W[i]>0.0+Wtol)) ){ + //determine coeff0icients + IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); + gsp[i] +=C; + if (gsp[i]<0.0+Wtol)gsp[i]=0.0; + } + /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents *from Marbouty, 1980: Figure 9*/ - - //_printf_("Nond}ritic snow metamorphism\n"); - Q = Marbouty(T[i],d[i],dT[i]); - - // calculate grain growth - gsz[i] += (Q*dt); - + else if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && dT[i]>15.0+Ttol) || (gsp[i]<=0.0+Gdntol && T[i]>5.0+Ttol && W[i]>0.0+Wtol)){ + //_printf_("Nondendritic snow metamorphism\n"); + Q = Marbouty(T[i],d[i],dT[i]); + + // calculate grain growth + gsz[i] += (Q*dt); + } //Wet snow metamorphism (Brun, 1989) - if(W[i]>0.0+Wtol){ - //_printf_("Nond}ritic wet snow metamorphism\n"); + else{ + //_printf_("Nondendritic wet snow metamorphism\n"); //wet rate of change coefficient IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts); // [mm^3 s^-1] @@ -310,5 +323,4 @@ // calculate change in grain volume and convert to grain size gsz[i] = 2.0 * pow(3.0/(Pi * 4.0)*((4.0/ 3.0)*Pi*pow(gsz[i]/2.0,3.0) + E),1.0/3.0); - } @@ -317,8 +329,8 @@ // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992) - if (fabs(gsp[i]-1.0)>Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0; - } - - //convert grain size back to effective grain radius: + if (fabs(gsp[i]-1.0)>=Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0; + } + + //convert grain size back to effective grain radius: re[i]=gsz[i]/2.0; } @@ -1092,5 +1104,5 @@ adjusts the properties of the top grid cell to account for accumulation T_air & T = Air and top grid cell temperatures [K] - Tmean = average surface temperature [K] + Vmean = average wind velocity [m s-1] V = wind velocity [m s-1] C = average accumulation rate [kg m-2 yr-1] @@ -1105,5 +1117,5 @@ // specify constants IssmDouble dSnow = 150; // density of snow [kg m-3] - const IssmDouble reNew = 0.1; // new snow grain size [mm] + const IssmDouble reNew = 0.05; // new snow grain size [mm] const IssmDouble gdnNew = 1.0; // new snow dendricity const IssmDouble gspNew = 0.5; // new snow sphericity @@ -1156,5 +1168,5 @@ //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W // 7.36x10-2 1.06x10-3 6.69x10-2 4.77x10-3 - dSnow=(7.362e-2 + 1.06e-3*Tmean + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000; + dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000; break; @@ -1175,4 +1187,12 @@ IssmDouble z_snow = P/dSnow; // depth of snow + IssmDouble dfall = gdnNew; + IssmDouble sfall = gspNew; + IssmDouble refall = reNew; + + //From Vionnet et al., 2012 (Crocus) + dfall = min(max(1.29 - 0.17*V,0.20),1.0); + sfall = min(max(0.08*V + 0.38,0.5),0.9); + refall = (0.1 + (1.0-dfall)*0.25 + (0.5-sfall)*0.1)/2.0; // if snow depth is greater than specified min dz, new cell created @@ -1184,7 +1204,7 @@ newcell(&W,0.0,top,m); //new cell W newcell(&a,aSnow,top,m); //new cell a - newcell(&re,reNew,top,m); //new cell grain size - newcell(&gdn,gdnNew,top,m); //new cell grain dendricity - newcell(&gsp,gspNew,top,m); //new cell grain sphericity + newcell(&re,refall,top,m); //new cell grain size + newcell(&gdn,dfall,top,m); //new cell grain dendricity + newcell(&gsp,sfall,top,m); //new cell grain sphericity m=m+1; } @@ -1202,7 +1222,7 @@ // adjust a, re, gdn & gsp if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass; - re[0] = (reNew * P + re[0] * mInit[0])/mass; - gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass; - gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass; + re[0] = (refall * P + re[0] * mInit[0])/mass; + gdn[0] = (dfall * P + gdn[0] * mInit[0])/mass; + gsp[0] = (sfall * P + gsp[0] * mInit[0])/mass; } }