![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/23390-24306/ISSM-24188-24189.diff
| Last change on this file was 24307, checked in by , 5 years ago | |
|---|---|
| File size: 8.6 KB | |
-
../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
247 247 for(int i=0;i<m;i++){ 248 248 if (gdn[i]>0.0+Gdntol){ 249 249 250 //_printf_("Dendritic dry snow metamorphism\n"); 250 if(W[i]<=0.0+Wtol){ 251 //_printf_("Dendritic dry snow metamorphism\n"); 252 //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1 253 if(dT[i]<=5.0+Ttol){ 254 //determine coefficients 255 IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt; 256 IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt; 257 //new dentricity and sphericity for dT < 5 degC m-1 258 gdn[i] += A; 259 gsp[i] += B; 260 } 261 else{ 262 // new dendricity and sphericity for dT >= 5 degC m-1 251 263 252 //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1 253 if(dT[i]<5.0-Ttol){ 254 //determine coefficients 255 IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt; 256 IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt; 257 //new dentricity and sphericity for dT < 5 degC m-1 258 gdn[i] += A; 259 gsp[i] += B; 264 //determine coeff0icients 265 IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); 266 gdn[i] +=C; 267 gsp[i] +=C; 268 } 260 269 } 261 270 else{ 262 // new dendricity and sphericity for dT >= 5 degC m-1 271 // wet snow metamorphism 272 //_printf_("Dendritic wet snow metamorphism\n"); 263 273 264 //determine coefficients265 IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);266 gdn[i] +=C;267 gsp[i] +=C;268 }269 270 // wet snow metamorphism271 if(W[i]>0.0+Wtol){272 273 //_printf_("D}ritic wet snow metamorphism\n");274 275 274 //determine coefficient 276 275 IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt; 277 276 278 // new dendricity and sphericityfor wet snow277 // new dendricity for wet snow 279 278 gdn[i] -= D; 279 // new sphericity for wet snow 280 280 gsp[i] += D; 281 281 } 282 // dendricity and sphericity can not be > 1 or < 0 283 if (gdn[i]<0.0+Wtol)gdn[i]=0.0; 284 if (gsp[i]<0.0+Wtol)gsp[i]=0.0; 285 if (gdn[i]>1.0-Wtol)gdn[i]=1.0; 286 if (gsp[i]>1.0-Wtol)gsp[i]=1.0; 282 287 283 // dendricity and sphericity can not be > 1 or < 0 284 if (gdn[i]<0.0+Wtol)gdn[i]=0.0; 285 if (gsp[i]<0.0+Wtol)gsp[i]=0.0; 286 if (gdn[i]>1.0-Wtol)gdn[i]=1.0; 287 if (gsp[i]>1.0-Wtol)gsp[i]=1.0; 288 // determine new grain size (mm) 289 gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1; 288 290 289 // determine new grain size (mm)290 gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1;291 292 291 } 293 292 else{ 294 293 295 /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficinets 294 //When the state of "faceted crystals" (gsp==0) is fully reached, 295 // snow evolves towards depth hoar if the gradient is 296 // higher than 15 degC m-1 (Brun et al., 1992) 297 //When wet-snow grains (class 6) are submitted to a 298 // temperature gradient higher than degC m-1, their sphericity 299 // decreases according to Equations (4). When sphericity 300 // reaches 0, their size increases according to the functions 301 // determined by Marbouty. (Brun et al., 1992) 302 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)) ){ 303 //determine coeff0icients 304 IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); 305 gsp[i] +=C; 306 if (gsp[i]<0.0+Wtol)gsp[i]=0.0; 307 } 308 /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents 296 309 *from Marbouty, 1980: Figure 9*/ 310 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)){ 311 //_printf_("Nondendritic snow metamorphism\n"); 312 Q = Marbouty(T[i],d[i],dT[i]); 297 313 298 //_printf_("Nond}ritic snow metamorphism\n"); 299 Q = Marbouty(T[i],d[i],dT[i]); 300 301 // calculate grain growth 302 gsz[i] += (Q*dt); 303 314 // calculate grain growth 315 gsz[i] += (Q*dt); 316 } 304 317 //Wet snow metamorphism (Brun, 1989) 305 if(W[i]>0.0+Wtol){306 //_printf_("Nond }ritic wet snow metamorphism\n");318 else{ 319 //_printf_("Nondendritic wet snow metamorphism\n"); 307 320 //wet rate of change coefficient 308 321 IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts); // [mm^3 s^-1] 309 322 310 323 // calculate change in grain volume and convert to grain size 311 324 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); 312 313 325 } 314 326 315 327 // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992) … … 316 328 if (fabs(gsp[i]-1.0)<Wtol && gsz[i]>2.0-Wtol) gsz[i]=2.0; 317 329 318 330 // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992) 319 if (fabs(gsp[i]-1.0)> Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;331 if (fabs(gsp[i]-1.0)>=Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0; 320 332 } 321 333 322 //convert grain size back to effective grain radius: 334 //convert grain size back to effective grain radius: 323 335 re[i]=gsz[i]/2.0; 324 336 } 325 337 … … 1091 1103 /* Description: 1092 1104 adjusts the properties of the top grid cell to account for accumulation 1093 1105 T_air & T = Air and top grid cell temperatures [K] 1094 Tmean = average surface temperature [K]1106 Vmean = average wind velocity [m s-1] 1095 1107 V = wind velocity [m s-1] 1096 1108 C = average accumulation rate [kg m-2 yr-1] 1097 1109 dz = topgrid cell length [m] … … 1104 1116 // MAIN FUNCTION 1105 1117 // specify constants 1106 1118 IssmDouble dSnow = 150; // density of snow [kg m-3] 1107 const IssmDouble reNew = 0. 1; // new snow grain size [mm]1119 const IssmDouble reNew = 0.05; // new snow grain size [mm] 1108 1120 const IssmDouble gdnNew = 1.0; // new snow dendricity 1109 1121 const IssmDouble gspNew = 0.5; // new snow sphericity 1110 1122 … … 1155 1167 case 3: //Surface snow accumulation density from Kaspers et al., 2004, Antarctica 1156 1168 //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W 1157 1169 // 7.36x10-2 1.06x10-3 6.69x10-2 4.77x10-3 1158 dSnow=(7.36 2e-2 + 1.06e-3*Tmean+ 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;1170 dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000; 1159 1171 break; 1160 1172 1161 1173 case 4: // Kuipers Munneke and others (2015), Greenland … … 1174 1186 if (T_air <= CtoK+Ttol){ // if snow 1175 1187 1176 1188 IssmDouble z_snow = P/dSnow; // depth of snow 1189 IssmDouble dfall = gdnNew; 1190 IssmDouble sfall = gspNew; 1191 IssmDouble refall = reNew; 1177 1192 1193 //From Vionnet et al., 2012 (Crocus) 1194 dfall = min(max(1.29 - 0.17*V,0.20),1.0); 1195 sfall = min(max(0.08*V + 0.38,0.5),0.9); 1196 refall = (0.1 + (1.0-dfall)*0.25 + (0.5-sfall)*0.1)/2.0; 1197 1178 1198 // if snow depth is greater than specified min dz, new cell created 1179 1199 if (z_snow > dzMin+Dtol){ 1180 1200 … … 1183 1203 newcell(&d,dSnow,top,m); //new cell d 1184 1204 newcell(&W,0.0,top,m); //new cell W 1185 1205 newcell(&a,aSnow,top,m); //new cell a 1186 newcell(&re,re New,top,m); //new cell grain size1187 newcell(&gdn, gdnNew,top,m); //new cell grain dendricity1188 newcell(&gsp, gspNew,top,m); //new cell grain sphericity1206 newcell(&re,refall,top,m); //new cell grain size 1207 newcell(&gdn,dfall,top,m); //new cell grain dendricity 1208 newcell(&gsp,sfall,top,m); //new cell grain sphericity 1189 1209 m=m+1; 1190 1210 } 1191 1211 else { // if snow depth is less than specified minimum dz snow … … 1201 1221 1202 1222 // adjust a, re, gdn & gsp 1203 1223 if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass; 1204 re[0] = (re New* P + re[0] * mInit[0])/mass;1205 gdn[0] = ( gdnNew* P + gdn[0] * mInit[0])/mass;1206 gsp[0] = ( gspNew* P + gsp[0] * mInit[0])/mass;1224 re[0] = (refall * P + re[0] * mInit[0])/mass; 1225 gdn[0] = (dfall * P + gdn[0] * mInit[0])/mass; 1226 gsp[0] = (sfall * P + gsp[0] * mInit[0])/mass; 1207 1227 } 1208 1228 } 1209 1229 else{ // if rain -
../trunk-jpl/src/c/classes/Elements/Element.cpp
4046 4046 #endif 4047 4047 4048 4048 /*Check bottom grid cell T is unchanged:*/ 4049 if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); 4049 if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){ 4050 if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); 4051 } 4050 4052 4051 4053 /*Free ressources: */ 4052 4054 xDelete<IssmDouble>(swf);
Note:
See TracBrowser
for help on using the repository browser.
