source: issm/oecreview/Archive/23390-24306/ISSM-24188-24189.diff

Last change on this file was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 8.6 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24188)
4+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24189)
5@@ -247,69 +247,81 @@
6 for(int i=0;i<m;i++){
7 if (gdn[i]>0.0+Gdntol){
8
9- //_printf_("Dendritic dry snow metamorphism\n");
10+ if(W[i]<=0.0+Wtol){
11+ //_printf_("Dendritic dry snow metamorphism\n");
12+ //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
13+ if(dT[i]<=5.0+Ttol){
14+ //determine coefficients
15+ IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
16+ IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt;
17+ //new dentricity and sphericity for dT < 5 degC m-1
18+ gdn[i] += A;
19+ gsp[i] += B;
20+ }
21+ else{
22+ // new dendricity and sphericity for dT >= 5 degC m-1
23
24- //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
25- if(dT[i]<5.0-Ttol){
26- //determine coefficients
27- IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
28- IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt;
29- //new dentricity and sphericity for dT < 5 degC m-1
30- gdn[i] += A;
31- gsp[i] += B;
32+ //determine coeff0icients
33+ IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
34+ gdn[i] +=C;
35+ gsp[i] +=C;
36+ }
37 }
38 else{
39- // new dendricity and sphericity for dT >= 5 degC m-1
40+ // wet snow metamorphism
41+ //_printf_("Dendritic wet snow metamorphism\n");
42
43- //determine coefficients
44- IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
45- gdn[i] +=C;
46- gsp[i] +=C;
47- }
48-
49- // wet snow metamorphism
50- if(W[i]>0.0+Wtol){
51-
52- //_printf_("D}ritic wet snow metamorphism\n");
53-
54 //determine coefficient
55 IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt;
56
57- // new dendricity and sphericity for wet snow
58+ // new dendricity for wet snow
59 gdn[i] -= D;
60+ // new sphericity for wet snow
61 gsp[i] += D;
62 }
63+ // dendricity and sphericity can not be > 1 or < 0
64+ if (gdn[i]<0.0+Wtol)gdn[i]=0.0;
65+ if (gsp[i]<0.0+Wtol)gsp[i]=0.0;
66+ if (gdn[i]>1.0-Wtol)gdn[i]=1.0;
67+ if (gsp[i]>1.0-Wtol)gsp[i]=1.0;
68
69- // dendricity and sphericity can not be > 1 or < 0
70- if (gdn[i]<0.0+Wtol)gdn[i]=0.0;
71- if (gsp[i]<0.0+Wtol)gsp[i]=0.0;
72- if (gdn[i]>1.0-Wtol)gdn[i]=1.0;
73- if (gsp[i]>1.0-Wtol)gsp[i]=1.0;
74+ // determine new grain size (mm)
75+ gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
76
77- // determine new grain size (mm)
78- gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
79-
80 }
81 else{
82
83- /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficinets
84+ //When the state of "faceted crystals" (gsp==0) is fully reached,
85+ // snow evolves towards depth hoar if the gradient is
86+ // higher than 15 degC m-1 (Brun et al., 1992)
87+ //When wet-snow grains (class 6) are submitted to a
88+ // temperature gradient higher than degC m-1, their sphericity
89+ // decreases according to Equations (4). When sphericity
90+ // reaches 0, their size increases according to the functions
91+ // determined by Marbouty. (Brun et al., 1992)
92+ 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)) ){
93+ //determine coeff0icients
94+ IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
95+ gsp[i] +=C;
96+ if (gsp[i]<0.0+Wtol)gsp[i]=0.0;
97+ }
98+ /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents
99 *from Marbouty, 1980: Figure 9*/
100+ 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)){
101+ //_printf_("Nondendritic snow metamorphism\n");
102+ Q = Marbouty(T[i],d[i],dT[i]);
103
104- //_printf_("Nond}ritic snow metamorphism\n");
105- Q = Marbouty(T[i],d[i],dT[i]);
106-
107- // calculate grain growth
108- gsz[i] += (Q*dt);
109-
110+ // calculate grain growth
111+ gsz[i] += (Q*dt);
112+ }
113 //Wet snow metamorphism (Brun, 1989)
114- if(W[i]>0.0+Wtol){
115- //_printf_("Nond}ritic wet snow metamorphism\n");
116+ else{
117+ //_printf_("Nondendritic wet snow metamorphism\n");
118 //wet rate of change coefficient
119 IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts); // [mm^3 s^-1]
120
121 // calculate change in grain volume and convert to grain size
122 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);
123-
124 }
125
126 // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
127@@ -316,10 +328,10 @@
128 if (fabs(gsp[i]-1.0)<Wtol && gsz[i]>2.0-Wtol) gsz[i]=2.0;
129
130 // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
131- if (fabs(gsp[i]-1.0)>Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;
132+ if (fabs(gsp[i]-1.0)>=Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;
133 }
134
135- //convert grain size back to effective grain radius:
136+ //convert grain size back to effective grain radius:
137 re[i]=gsz[i]/2.0;
138 }
139
140@@ -1091,7 +1103,7 @@
141 /* Description:
142 adjusts the properties of the top grid cell to account for accumulation
143 T_air & T = Air and top grid cell temperatures [K]
144- Tmean = average surface temperature [K]
145+ Vmean = average wind velocity [m s-1]
146 V = wind velocity [m s-1]
147 C = average accumulation rate [kg m-2 yr-1]
148 dz = topgrid cell length [m]
149@@ -1104,7 +1116,7 @@
150 // MAIN FUNCTION
151 // specify constants
152 IssmDouble dSnow = 150; // density of snow [kg m-3]
153- const IssmDouble reNew = 0.1; // new snow grain size [mm]
154+ const IssmDouble reNew = 0.05; // new snow grain size [mm]
155 const IssmDouble gdnNew = 1.0; // new snow dendricity
156 const IssmDouble gspNew = 0.5; // new snow sphericity
157
158@@ -1155,7 +1167,7 @@
159 case 3: //Surface snow accumulation density from Kaspers et al., 2004, Antarctica
160 //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W
161 // 7.36x10-2 1.06x10-3 6.69x10-2 4.77x10-3
162- dSnow=(7.362e-2 + 1.06e-3*Tmean + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
163+ dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
164 break;
165
166 case 4: // Kuipers Munneke and others (2015), Greenland
167@@ -1174,7 +1186,15 @@
168 if (T_air <= CtoK+Ttol){ // if snow
169
170 IssmDouble z_snow = P/dSnow; // depth of snow
171+ IssmDouble dfall = gdnNew;
172+ IssmDouble sfall = gspNew;
173+ IssmDouble refall = reNew;
174
175+ //From Vionnet et al., 2012 (Crocus)
176+ dfall = min(max(1.29 - 0.17*V,0.20),1.0);
177+ sfall = min(max(0.08*V + 0.38,0.5),0.9);
178+ refall = (0.1 + (1.0-dfall)*0.25 + (0.5-sfall)*0.1)/2.0;
179+
180 // if snow depth is greater than specified min dz, new cell created
181 if (z_snow > dzMin+Dtol){
182
183@@ -1183,9 +1203,9 @@
184 newcell(&d,dSnow,top,m); //new cell d
185 newcell(&W,0.0,top,m); //new cell W
186 newcell(&a,aSnow,top,m); //new cell a
187- newcell(&re,reNew,top,m); //new cell grain size
188- newcell(&gdn,gdnNew,top,m); //new cell grain dendricity
189- newcell(&gsp,gspNew,top,m); //new cell grain sphericity
190+ newcell(&re,refall,top,m); //new cell grain size
191+ newcell(&gdn,dfall,top,m); //new cell grain dendricity
192+ newcell(&gsp,sfall,top,m); //new cell grain sphericity
193 m=m+1;
194 }
195 else { // if snow depth is less than specified minimum dz snow
196@@ -1201,9 +1221,9 @@
197
198 // adjust a, re, gdn & gsp
199 if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
200- re[0] = (reNew * P + re[0] * mInit[0])/mass;
201- gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass;
202- gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass;
203+ re[0] = (refall * P + re[0] * mInit[0])/mass;
204+ gdn[0] = (dfall * P + gdn[0] * mInit[0])/mass;
205+ gsp[0] = (sfall * P + gsp[0] * mInit[0])/mass;
206 }
207 }
208 else{ // if rain
209Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
210===================================================================
211--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24188)
212+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24189)
213@@ -4046,7 +4046,9 @@
214 #endif
215
216 /*Check bottom grid cell T is unchanged:*/
217- if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
218+ if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
219+ if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
220+ }
221
222 /*Free ressources: */
223 xDelete<IssmDouble>(swf);
224Index: ../trunk-jpl/test/Archives/Archive243.arch
225===================================================================
226Cannot display: file marked as a binary type.
227svn:mime-type = application/octet-stream
228Index: ../trunk-jpl/test/Archives/Archive244.arch
229===================================================================
230Cannot display: file marked as a binary type.
231svn:mime-type = application/octet-stream
Note: See TracBrowser for help on using the repository browser.