source: issm/oecreview/Archive/21724-22754/ISSM-22449-22450.diff

Last change on this file was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 34.0 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/m/classes/SMBgemb.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/SMBgemb.m (revision 22449)
4+++ ../trunk-jpl/src/m/classes/SMBgemb.m (revision 22450)
5@@ -35,6 +35,10 @@
6 C = NaN; %mean annual snow accumulation [kg m-2 yr-1]
7 Tz = NaN; %height above ground at which temperature (T) was sampled [m]
8 Vz = NaN; %height above ground at which wind (V) eas sampled [m]
9+
10+ %optional inputs:
11+ aValue = NaN; %Albedo forcing at every element. Used only if aIdx == 0.
12+ teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
13
14 % Initialization of snow properties
15 Dzini = NaN; %cell depth (m)
16@@ -50,7 +54,8 @@
17
18 %settings:
19 aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1)
20- % 1: effective grain radius [Gardner & Sharp, 2009]
21+ % 0: direct input from aValue parameter
22+ % 1: effective grain radius [Gardner & Sharp, 2009]
23 % 2: effective grain radius [Brun et al., 2009]
24 % 3: density and cloud amount [Greuell & Konzelmann, 1994]
25 % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
26@@ -114,6 +119,14 @@
27 self.eAir=project3d(md,'vector',self.eAir,'type','node');
28 self.pAir=project3d(md,'vector',self.pAir,'type','node');
29
30+ if aIdx == 0 & ~isnan(self.aValue)
31+ self.aValue=project3d(md,'vector',self.aValue,'type','node');
32+ end
33+ if ~isnan(self.teValue)
34+ self.teValue=project3d(md,'vector',self.teValue,'type','node');
35+ end
36+
37+
38 end % }}}
39 function list = defaultoutputs(self,md) % {{{
40 list = {'SmbMassBalance'};
41@@ -149,6 +162,9 @@
42 self.t0wet = 15;
43 self.t0dry = 30;
44 self.K = 7;
45+
46+ self.teValue = ones(mesh.numberofelements,1);
47+ self.aValue = self.aSnow*ones(mesh.numberofelements,1);
48
49 self.Dzini=0.05*ones(mesh.numberofelements,2);
50 self.Dini=910.0*ones(mesh.numberofelements,2);
51@@ -166,7 +182,6 @@
52 end % }}}
53 function md = checkconsistency(self,md,solution,analyses) % {{{
54
55-
56 md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
57 md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
58 md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
59@@ -188,7 +203,9 @@
60 md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000);
61 md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000);
62
63- md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]);
64+ md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
65+
66+ md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]);
67 md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
68 md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);
69
70@@ -200,6 +217,8 @@
71 md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
72
73 switch self.aIdx,
74+ case 0
75+ md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
76 case {1 2}
77 md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
78 md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
79@@ -251,13 +270,16 @@
80 fielddisplay(self,'InitDensityScaling',{'initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'});
81 fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)');
82 fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
83+ '0: direct input from aValue parameter',...
84 '1: effective grain radius [Gardner & Sharp, 2009]',...
85 '2: effective grain radius [Brun et al., 2009]',...
86 '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
87 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'});
88+
89+ fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)');
90
91 %snow properties init
92- fielddisplay(self,'Dzini','Initial cel depth when restart [m]');
93+ fielddisplay(self,'Dzini','Initial cell depth when restart [m]');
94 fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]');
95 fielddisplay(self,'Reini','Initial grain size when restart [mm]');
96 fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]');
97@@ -271,6 +293,8 @@
98
99 %additional albedo parameters:
100 switch self.aIdx
101+ case 0
102+ fielddisplay(self,'aValue','Albedo forcing at every element. Used only if aIdx == 0.');
103 case {1 2}
104 fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)');
105 fielddisplay(self,'aIce','albedo of ice (0.27-0.58)');
106@@ -338,6 +362,9 @@
107 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double');
108 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double');
109 WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
110+
111+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
112+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
113
114 %snow properties init
115 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3);
116Index: ../trunk-jpl/src/m/classes/SMBgemb.py
117===================================================================
118--- ../trunk-jpl/src/m/classes/SMBgemb.py (revision 22449)
119+++ ../trunk-jpl/src/m/classes/SMBgemb.py (revision 22450)
120@@ -41,6 +41,10 @@
121 C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1]
122 Tz = float('NaN') #height above ground at which temperature (T) was sampled [m]
123 Vz = float('NaN') #height above ground at which wind (V) eas sampled [m]
124+
125+ #optional inputs:
126+ aValue = float('NaN') #Albedo forcing at every element. Used only if aIdx == 0.
127+ teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
128
129 # Initialization of snow properties
130 Dzini = float('NaN') #cell depth (m)
131@@ -56,7 +60,8 @@
132
133 #settings:
134 aIdx = float('NaN') #method for calculating albedo and subsurface absorption (default is 1)
135- # 1: effective grain radius [Gardner & Sharp, 2009]
136+ # 0: direct input from aValue parameter
137+ # 1: effective grain radius [Gardner & Sharp, 2009]
138 # 2: effective grain radius [Brun et al., 2009]
139 # 3: density and cloud amount [Greuell & Konzelmann, 1994]
140 # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
141@@ -134,13 +139,16 @@
142 string = "%s\n%s"%(string,fielddisplay(self,'InitDensityScaling',['initial scaling factor multiplying the density of ice','which describes the density of the snowpack.']))
143 string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'))
144 string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)',
145+ '0: direct input from aValue parameter',
146 '1: effective grain radius [Gardner & Sharp, 2009]',
147 '2: effective grain radius [Brun et al., 2009]',
148 '3: density and cloud amount [Greuell & Konzelmann, 1994]',
149 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
150+
151+ string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'))
152
153 #snow properties init
154- string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cel depth when restart [m]'))
155+ string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cell depth when restart [m]'))
156 string = "%s\n%s"%(string,fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]'))
157 string = "%s\n%s"%(string,fielddisplay(self,'Reini','Initial grain size when restart [mm]'))
158 string = "%s\n%s"%(string,fielddisplay(self,'Gdnini','Initial grain dricity when restart [-]'))
159@@ -155,6 +163,8 @@
160 if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)):
161 string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'))
162 string = "%s\n%s"%(string,fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'))
163+ elif self.aIdx == 0:
164+ string = "%s\n%s"%(string,fielddisplay(self,'aValue','Albedo forcing at every element. Used only if aIdx == 0.'))
165 elif self.aIdx == 3:
166 string = "%s\n%s"%(string,fielddisplay(self,'cldFrac','average cloud amount'))
167 elif self.aIdx == 4:
168@@ -182,6 +192,12 @@
169 self.P = project3d(md,'vector',self.P,'type','node')
170 self.eAir = project3d(md,'vector',self.eAir,'type','node')
171 self.pAir = project3d(md,'vector',self.pAir,'type','node')
172+
173+ if aIdx == 0 and np.isnan(self.aValue):
174+ self.aValue=project3d(md,'vector',self.aValue,'type','node');
175+ if np.isnan(self.teValue):
176+ self.teValue=project3d(md,'vector',self.teValue,'type','node');
177+
178 return self
179 #}}}
180 def defaultoutputs(self,md): # {{{
181@@ -210,7 +226,7 @@
182 self.zMin = 130*np.ones((mesh.numberofelements,))
183 self.zY = 1.10*np.ones((mesh.numberofelements,))
184 self.outputFreq = 30
185-
186+
187 #additional albedo parameters
188 self.aSnow = 0.85
189 self.aIce = 0.48
190@@ -218,6 +234,9 @@
191 self.t0wet = 15
192 self.t0dry = 30
193 self.K = 7
194+
195+ self.teValue = np.ones((mesh.numberofelements,));
196+ self.aValue = self.aSnow*np.ones(mesh.numberofelements,);
197
198 self.Dzini = 0.05*np.ones((mesh.numberofelements,2))
199 self.Dini = 910.0*np.ones((mesh.numberofelements,2))
200@@ -256,7 +275,9 @@
201 md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000)
202 md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000)
203
204- md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4])
205+ md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
206+
207+ md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4])
208 md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1])
209 md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5])
210
211@@ -270,6 +291,8 @@
212 if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)):
213 md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'> = ',.64,'< = ',.89)
214 md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'> = ',.27,'< = ',.58)
215+ elif self.aIdx == 0:
216+ md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
217 elif self.aIdx == 3:
218 md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'> = ',0,'< = ',1)
219 elif self.aIdx == 4:
220@@ -332,6 +355,9 @@
221 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double')
222 WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double')
223 WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double')
224+
225+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
226+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
227
228 #snow properties init
229 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3)
230Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp
231===================================================================
232--- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22449)
233+++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22450)
234@@ -68,6 +68,8 @@
235 iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum);
236 iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum);
237 iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum);
238+ iomodel->FetchDataToInput(elements,"md.smb.aValue",SmbAValueEnum);
239+ iomodel->FetchDataToInput(elements,"md.smb.teValue",SmbTeValueEnum);
240 break;
241 case SMBpddEnum:
242 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
243Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
244===================================================================
245--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22449)
246+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 22450)
247@@ -24,10 +24,10 @@
248 void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY);
249 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
250 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);
251-void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
252+void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
253 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
254-void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);
255-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid);
256+void thermo(IssmDouble* pEC, IssmDouble** T, 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 dIce, int sid);
257+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx,IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid);
258 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);
259 void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
260 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);
261Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
262===================================================================
263--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22449)
264+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22450)
265@@ -336,9 +336,10 @@
266 *pgsp=gsp;
267
268 } /*}}}*/
269-void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
270+void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
271
272 //// Calculates Snow, firn and ice albedo as a function of:
273+ // 0 : direct input from aValue parameter
274 // 1 : effective grain radius (Gardner & Sharp, 2009)
275 // 2 : effective grain radius (Brun et al., 2009)
276 // 3 : density and cloud amount (Greuell & Konzelmann, 1994)
277@@ -347,6 +348,9 @@
278 //// Inputs
279 // aIdx = albedo method to use
280
281+ // Method 0
282+ // aValue = direct input value for albedo
283+
284 // Methods 1 & 2
285 // re = surface effective grain radius [mm]
286
287@@ -391,7 +395,10 @@
288 //some constants:
289 const IssmDouble dSnow = 300; // density of fresh snow [kg m-3]
290
291- if(aIdx==1){
292+ if (aIdx==0){
293+ for(int i=0;i<m;i++)a[i] = aValue;
294+ }
295+ else if(aIdx==1){
296 //function of effective grain radius
297
298 //convert effective radius to specific surface area [cm2 g-1]
299@@ -399,6 +406,7 @@
300
301 //determine broadband albedo
302 a[0]= 1.48 - pow(S,-0.07);
303+ for(int i=1;i<m;i++)a[i]=a[0];
304 }
305 else if(aIdx==2){
306
307@@ -420,7 +428,7 @@
308
309 // broadband surface albedo
310 a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
311-
312+ for(int i=1;i<m;i++)a[i]=a[0];
313 }
314 else if(aIdx==3){
315
316@@ -428,6 +436,7 @@
317
318 // calculate albedo
319 a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
320+ for(int i=1;i<m;i++)a[i]=a[0];
321 }
322 else if(aIdx==4){
323
324@@ -485,7 +494,7 @@
325 xDelete<IssmDouble>(d_a);
326
327 }
328- else _error_("albedo method switch should range from 1 to 4!");
329+ else _error_("albedo method switch should range from 0 to 4!");
330
331 // Check for erroneous values
332 if (a[0] > 1) _printf_("albedo > 1.0\n");
333@@ -496,7 +505,7 @@
334 *pa=a;
335
336 } /*}}}*/
337-void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/
338+void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/
339
340 /* ENGLACIAL THERMODYNAMICS*/
341
342@@ -807,7 +816,7 @@
343 dT_turb = turb / TCs;
344
345 // upward longwave contribution
346- ulw = - SB * pow(Ts,4.0) * dt;
347+ ulw = - SB * pow(Ts,4.0)* teValue * dt;
348 dT_ulw = ulw / TCs;
349
350 // new grid point temperature
351@@ -1068,7 +1077,7 @@
352 *pswf=swf;
353
354 } /*}}}*/
355-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/
356+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/
357
358 // Adds precipitation and deposition to the model grid
359
360@@ -1161,7 +1170,7 @@
361 T[0] = (T_air * P + T[0] * mInit[0])/mass;
362
363 // adjust a, re, gdn & gsp
364- a[0] = (aSnow * P + a[0] * mInit[0])/mass;
365+ if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
366 re[0] = (reNew * P + re[0] * mInit[0])/mass;
367 gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass;
368 gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass;
369@@ -1786,7 +1795,7 @@
370 } /*}}}*/
371 void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/
372
373- //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???]
374+ //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPACTION IS COMPENSATED FOR BY TRACES OF SNOW???]
375
376 //// FUNCTION INFO
377
378Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
379===================================================================
380--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22449)
381+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22450)
382@@ -2555,6 +2555,8 @@
383 IssmDouble P=0.0;
384 IssmDouble eAir=0.0;
385 IssmDouble pAir=0.0;
386+ IssmDouble teValue=1.0;
387+ IssmDouble aValue=0.0;
388 int aIdx=0;
389 int denIdx=0;
390 int swIdx=0;
391@@ -2658,10 +2660,14 @@
392 Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input);
393 Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input);
394 Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input);
395+ Input* teValue_input=this->GetInput(SmbTeValueEnum); _assert_(teValue_input);
396+ Input* aValue_input=this->GetInput(SmbAValueEnum); _assert_(aValue_input);
397 Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input);
398 /*Retrieve input values:*/
399 Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
400
401+ if (aIdx == 0) aValue_input->GetInputValue(&aValue,gauss);
402+
403 zTop_input->GetInputValue(&zTop,gauss);
404 dzTop_input->GetInputValue(&dzTop,gauss);
405 dzMin_input->GetInputValue(&dzMin,gauss);
406@@ -2672,6 +2678,7 @@
407 C_input->GetInputValue(&C,gauss);
408 Tz_input->GetInputValue(&Tz,gauss);
409 Vz_input->GetInputValue(&Vz,gauss);
410+ teValue_input->GetInputValue(&teValue,gauss);
411 isinitialized_input->GetInputValue(&isinitialized);
412 /*}}}*/
413
414@@ -2718,8 +2725,9 @@
415 W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0]; //set water content to zero [kg m-2]
416 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0]; //set albedo equal to fresh snow [fraction]
417 T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean; //set initial grid cell temperature to the annual mean temperature [K]
418- /* /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
419- Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
420+ /*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties)
421+ * because don't know Tmean yet when set default values.
422+ * Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
423
424 //fixed lower temperature bounday condition - T is fixed
425 T_bottom=T[m-1];
426@@ -2796,6 +2804,8 @@
427 P_input->GetInputValue(&P,gauss,t); //precipitation [kg m-2]
428 eAir_input->GetInputValue(&eAir,gauss,t); //screen level vapor pressure [Pa]
429 pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa]
430+ teValue_input->GetInputValue(&teValue,gauss); // screen level air pressure [Pa]
431+ if(aIdx == 0) aValue_input->GetInputValue(&aValue,gauss); // screen level air pressure [Pa]
432 //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
433 /*}}}*/
434
435@@ -2803,7 +2813,7 @@
436 if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
437
438 /*Snow, firn and ice albedo:*/
439- if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
440+ if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce, aSnow,aValue,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
441
442
443 /*Distribution of absorbed short wave radation with depth:*/
444@@ -2813,7 +2823,7 @@
445 netSW = cellsum(swf,m);
446
447 /*Thermal profile computation:*/
448- if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid());
449+ if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid());
450
451 /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell.
452 * need to fix this in case all or more of cell evaporates */
453@@ -2820,7 +2830,7 @@
454 dz[0] = dz[0] + EC / d[0];
455
456 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
457- if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
458+ if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
459
460 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
461 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
462@@ -2831,7 +2841,7 @@
463
464 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
465 * sub-time step in thermo equations*/
466- ulw = 5.67E-8 * pow(T[0],4.0);
467+ ulw = 5.67E-8 * pow(T[0],4.0) * teValue;
468
469 /*Calculate net longwave [W m-2]*/
470 netLW = dlw - ulw;
471@@ -2839,7 +2849,7 @@
472 /*Calculate turbulent heat fluxes [W m-2]*/
473 if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
474
475- /*Verbose some resuls in debug mode: {{{*/
476+ /*Verbose some results in debug mode: {{{*/
477 if(VerboseSmb() && 0){
478 _printf_("smb log: count[" << count << "] m[" << m << "] "
479 << setprecision(16) << "T[" << cellsum(T,m) << "] "
480@@ -2851,6 +2861,8 @@
481 << "gdn[" << cellsum(gdn,m) << "] "
482 << "gsp[" << cellsum(gsp,m) << "] "
483 << "swf[" << netSW << "] "
484+ << "a[" << a << "] "
485+ << "te[" << teValue << "] "
486 << "\n");
487 } /*}}}*/
488
489Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
490===================================================================
491--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22449)
492+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22450)
493@@ -458,6 +458,8 @@
494 SmbWEnum,
495 SmbAEnum,
496 SmbTEnum,
497+ SmbAValueEnum,
498+ SmbTeValueEnum,
499 SmbIsgraingrowthEnum,
500 SmbIsalbedoEnum,
501 SmbIsshortwaveEnum,
502Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
503===================================================================
504--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22449)
505+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22450)
506@@ -460,6 +460,8 @@
507 case SmbWEnum : return "SmbW";
508 case SmbAEnum : return "SmbA";
509 case SmbTEnum : return "SmbT";
510+ case SmbAValueEnum : return "SmbAValue";
511+ case SmbTeValueEnum : return "SmbTeValue";
512 case SmbIsgraingrowthEnum : return "SmbIsgraingrowth";
513 case SmbIsalbedoEnum : return "SmbIsalbedo";
514 case SmbIsshortwaveEnum : return "SmbIsshortwave";
515Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
516===================================================================
517--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22449)
518+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22450)
519@@ -469,6 +469,8 @@
520 else if (strcmp(name,"SmbW")==0) return SmbWEnum;
521 else if (strcmp(name,"SmbA")==0) return SmbAEnum;
522 else if (strcmp(name,"SmbT")==0) return SmbTEnum;
523+ else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
524+ else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
525 else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
526 else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
527 else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum;
528@@ -503,12 +505,12 @@
529 else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
530 else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
531 else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
532- else if (strcmp(name,"SmbF")==0) return SmbFEnum;
533- else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
534 else stage=5;
535 }
536 if(stage==5){
537- if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
538+ if (strcmp(name,"SmbF")==0) return SmbFEnum;
539+ else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
540+ else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
541 else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
542 else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
543 else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
544@@ -626,12 +628,12 @@
545 else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
546 else if (strcmp(name,"GiaW")==0) return GiaWEnum;
547 else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
548- else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
549- else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
550 else stage=6;
551 }
552 if(stage==6){
553- if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
554+ if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
555+ else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
556+ else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
557 else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
558 else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
559 else if (strcmp(name,"J")==0) return JEnum;
560@@ -749,12 +751,12 @@
561 else if (strcmp(name,"VxObs")==0) return VxObsEnum;
562 else if (strcmp(name,"VyObs")==0) return VyObsEnum;
563 else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
564- else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
565- else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
566 else stage=7;
567 }
568 if(stage==7){
569- if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
570+ if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
571+ else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
572+ else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
573 else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum;
574 else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum;
575 else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum;
576@@ -872,12 +874,12 @@
577 else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
578 else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
579 else if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum;
580- else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
581- else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
582 else stage=8;
583 }
584 if(stage==8){
585- if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
586+ if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
587+ else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
588+ else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum;
589 else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
590 else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
591 else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum;
592@@ -995,12 +997,12 @@
593 else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
594 else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
595 else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
596- else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
597- else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
598 else stage=9;
599 }
600 if(stage==9){
601- if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
602+ if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
603+ else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
604+ else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
605 else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
606 else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
607 else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
608@@ -1118,12 +1120,12 @@
609 else if (strcmp(name,"Water")==0) return WaterEnum;
610 else if (strcmp(name,"DataSet")==0) return DataSetEnum;
611 else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
612- else if (strcmp(name,"Loads")==0) return LoadsEnum;
613- else if (strcmp(name,"Materials")==0) return MaterialsEnum;
614 else stage=10;
615 }
616 if(stage==10){
617- if (strcmp(name,"Nodes")==0) return NodesEnum;
618+ if (strcmp(name,"Loads")==0) return LoadsEnum;
619+ else if (strcmp(name,"Materials")==0) return MaterialsEnum;
620+ else if (strcmp(name,"Nodes")==0) return NodesEnum;
621 else if (strcmp(name,"Contours")==0) return ContoursEnum;
622 else if (strcmp(name,"Parameters")==0) return ParametersEnum;
623 else if (strcmp(name,"Vertices")==0) return VerticesEnum;
Note: See TracBrowser for help on using the repository browser.