source: issm/oecreview/Archive/24307-24683/ISSM-24507-24508.diff@ 24684

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

CHG: added new review
File size: 5.4 KB
RevLine 
[24684]1Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24507)
4+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24508)
5@@ -493,7 +493,7 @@
6 for(int i=0;i<m;i++)a[i] -= d_a[i]; // new albedo
7
8 // modification of albedo due to thin layer of snow or solid
9- // condensation (deposition) at the surface surface
10+ // condensation (deposition) at the surface
11
12 // check if condensation occurs & if it is deposited in solid phase
13 if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0; // add cond to precip [mm]
14@@ -1339,6 +1339,8 @@
15
16 /*outputs:*/
17 IssmDouble mAdd = 0.0;
18+ IssmDouble surplusE = 0.0;
19+ IssmDouble surplusT = 0.0;
20 IssmDouble dz_add = 0.0;
21 IssmDouble Rsum = 0.0;
22 IssmDouble* T=*pT;
23@@ -1432,15 +1434,25 @@
24
25 int i = 0;
26 while (cellsum(surpE,n) > 0.0+Ttol && i<n){
27- // use surplus energy to increase the temperature of lower cell
28- T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
29
30- exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
31- T[i+1] = min(CtoK, T[i+1]);
32+ if (i<n-1){
33+ // use surplus energy to increase the temperature of lower cell
34+ T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
35
36- surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
37- surpE[i+1] = surpT[i+1] * CI * m[i+1];
38+ exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
39+ T[i+1] = min(CtoK, T[i+1]);
40
41+ surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
42+ surpE[i+1] = surpT[i+1] * CI * m[i+1];
43+ }
44+ else{
45+ surplusT=max(0.0, exsT[i] - LF/CI);
46+ surplusE=surpE[i];
47+ if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0){
48+ _printf0_(" WARNING: surplus energy at the base of GEMB column\n");
49+ }
50+ }
51+
52 // adjust current cell properties (again 159.1342 is the max T)
53 exsT[i] = LF/CI;
54 surpE[i] = 0.0;
55@@ -1648,7 +1660,7 @@
56 W[i+1] = W[i+1] + W[i]; // combine liquid water
57 m[i+1] = m_new; // combine top masses
58
59- // set cell to 99999 for deletion
60+ // set cell to -99999 for deletion
61 m[i] = Delflag;
62 }
63 }
64@@ -1655,7 +1667,7 @@
65
66 //If last cell has to be merged
67 if(lastCellFlag){
68- //find closest cell to merge with
69+ //find closest cell to merge with
70 for(int i=n-2;i>=0;i--){
71 if(m[i]!=Delflag){
72 X2=i;
73@@ -1678,7 +1690,7 @@
74 W[X1] = W[X1] + W[X2]; // combine liquid water
75 m[X1] = m_new; // combine top masses
76
77- // set cell to 99999 for deletion
78+ // set cell to -99999 for deletion
79 m[X2] = Delflag;
80 }
81
82@@ -1703,7 +1715,7 @@
83
84 // check if any of the top 10 cell depths are too large
85 X=0;
86- for(int i=9;i>=0;i--){
87+ for(int i=min(9,n-1);i>=0;i--){
88 if(dz[i]> 2.0*dzMin+Dtol){
89 X=i;
90 break;
91@@ -1714,21 +1726,21 @@
92 while(j<=X){
93 if (dz[j] > dzMin*2.0+Dtol){
94
95- // _printf_("dz > dzMin * 2");
96- // split in two
97- cellsplit(&dz, n, j,.5);
98- cellsplit(&W, n, j,.5);
99- cellsplit(&m, n, j,.5);
100- cellsplit(&T, n, j,1.0);
101- cellsplit(&d, n, j,1.0);
102- cellsplit(&a, n, j,1.0);
103- cellsplit(&EI, n, j,.5);
104- cellsplit(&EW, n, j,.5);
105- cellsplit(&re, n, j,1.0);
106- cellsplit(&gdn, n, j,1.0);
107- cellsplit(&gsp, n, j,1.0);
108- n++;
109- X=X+1;
110+ // _printf_("dz > dzMin * 2");
111+ // split in two
112+ cellsplit(&dz, n, j,.5);
113+ cellsplit(&W, n, j,.5);
114+ cellsplit(&m, n, j,.5);
115+ cellsplit(&T, n, j,1.0);
116+ cellsplit(&d, n, j,1.0);
117+ cellsplit(&a, n, j,1.0);
118+ cellsplit(&EI, n, j,.5);
119+ cellsplit(&EW, n, j,.5);
120+ cellsplit(&re, n, j,1.0);
121+ cellsplit(&gdn, n, j,1.0);
122+ cellsplit(&gsp, n, j,1.0);
123+ n++;
124+ X=X+1;
125 }
126 else j++;
127 }
128@@ -1816,7 +1828,7 @@
129 /*only in forward mode! avoid round in AD mode as it is not differentiable: */
130 #ifndef _HAVE_AD_
131 dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
132- dE = round(sumE0 - sumE1 - sumER + addE);
133+ dE = round(sumE0 - sumE1 - sumER + addE - surplusE);
134 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
135 << "dm: " << dm << " dE: " << dE << "\n");
136 #endif
137@@ -1977,8 +1989,14 @@
138 H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
139 c0arth = 0.07 * H;
140 c1arth = 0.03 * H;
141+ //ERA-5
142+ //M0 = max(2.3999 - (0.2610 * log(C)),0.25);
143+ //M1 = max(2.7469 - (0.3228 * log(C)),0.25);
144+ //RACMO
145 M0 = max(1.6599 - (0.1724 * log(C)),0.25);
146 M1 = max(2.0102 - (0.2458 * log(C)),0.25);
147+ //From Ligtenberg
148+ //H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
149 //M0 = max(1.435 - (0.151 * log(C)),0.25);
150 //M1 = max(2.366 - (0.293 * log(C)),0.25);
151 c0 = M0*c0arth;
152@@ -1991,10 +2009,15 @@
153 H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
154 c0arth = 0.07 * H;
155 c1arth = 0.03 * H;
156+ // ERA5
157+ //M0 = max(1.8920 - (0.1569 * log(C)),0.25);
158+ //M1 = max(2.5662 - (0.2274 * log(C)),0.25);
159+ // RACMO
160+ M0 = max(1.6201 - (0.1450 * log(C)),0.25);
161+ M1 = max(2.5577 - (0.2899 * log(C)),0.25);
162+ // From Kuipers Munneke
163 //M0 = max(1.042 - (0.0916 * log(C)),0.25);
164 //M1 = max(1.734 - (0.2039 * log(C)),0.25);
165- M0 = max(1.6201 - (0.1450 * log(C)),0.25);
166- M1 = max(2.5577 - (0.2899 * log(C)),0.25);
167 c0 = M0*c0arth;
168 c1 = M1*c1arth;
169 break;
Note: See TracBrowser for help on using the repository browser.