source: issm/trunk/src/c/modules/SurfaceMassBalancex/Gembx.cpp@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 58.1 KB
Line 
1/*!\file GEMB module from Alex Gardner.
2 * \brief: calculates SMB
3 */
4
5#include "./SurfaceMassBalancex.h"
6#include "../../shared/shared.h"
7#include "../../toolkits/toolkits.h"
8
9const double Pi = 3.141592653589793;
10
11void Gembx(FemModel* femmodel){ /*{{{*/
12
13 for(int i=0;i<femmodel->elements->Size();i++){
14 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
15 element->SmbGemb();
16 }
17
18} /*}}}*/
19void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY){ /*{{{*/
20
21 /* This file sets up the initial grid spacing and total grid depth. The
22 grid structure is set as constant grid length 'dzTop' for the top
23 'zTop' meters of the model grid. Bellow 'zTop' the gid length increases
24 linearly with depth */
25
26 /*intermediary:*/
27 IssmDouble dgpTop;
28 int gpTop, gpBottom;
29 int i;
30 IssmDouble gp0,z0;
31 IssmDouble* dzT=NULL;
32 IssmDouble* dzB=NULL;
33
34 /*output: */
35 int m;
36 IssmDouble* dz=NULL;
37
38 //----------------------Calculate Grid Lengths------------------------------
39 //calculate number of top grid points
40 dgpTop = zTop/dzTop;
41
42 //check to see if the top grid cell structure length (dzTop) goes evenly
43 //into specified top structure depth (zTop). Also make sure top grid cell
44 //structure length (dzTop) is greater than 5 cm
45 #ifndef _HAVE_ADOLC_ //avoid the round operation check!
46 if (dgpTop != round(dgpTop)){
47 _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
48 }
49 #endif
50 if(dzTop < 0.05){
51 _printf_("initial top grid cell length (dzTop) is < 0.05 m");
52 }
53 gpTop=reCast<int,IssmDouble>(dgpTop);
54
55 //initialize top grid depth vector
56 dzT = xNew<IssmDouble>(gpTop);
57 for (i=0;i<gpTop;i++)dzT[i]=dzTop;
58
59 //bottom grid cell depth = x*zY^(cells from to structure)
60 //figure out the number of grid points in the bottom vector (not known a priori)
61 gp0 = dzTop;
62 z0 = zTop;
63 gpBottom = 0;
64 while (zMax > z0){
65 gp0= gp0 * zY;
66 z0 = z0 + gp0;
67 gpBottom++;
68 }
69 //initialize bottom vectors
70 dzB = xNewZeroInit<IssmDouble>(gpBottom);
71 gp0 = dzTop;
72 z0 = zTop;
73 for(i=0;i<gpBottom;i++){
74 gp0=gp0*zY;
75 dzB[i]=gp0;
76 }
77
78 //combine top and bottom dz vectors
79 dz = xNew<IssmDouble>(gpTop+gpBottom);
80 for(i=0;i<gpTop;i++){
81 dz[i]=dzT[i];
82 }
83 for(i=0;i<gpBottom;i++){
84 dz[gpTop+i]=dzB[i];
85 }
86
87 /*Free ressouces:*/
88 xDelete(dzT);
89 xDelete(dzB);
90
91 //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
92
93 /*assign ouput pointers: */
94 *pdz=dz;
95 *psize=gpTop+gpBottom;
96} /*}}}*/
97IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT){ /*{{{*/
98
99 // calculates grain growth according to Fig. 9 of Marbouty, 1980
100 // ------NO GRAIN GROWTH FOR d > 400 kg m-3 because (H is set to zero)------
101 // ---------------this is a major limitation of the model-------------------
102
103 // initialize
104 IssmDouble F = 0, H=0, G=0;
105 const IssmDouble E = 0.09; //[mm d-1] model time growth constant E
106 // convert T from K to ºC
107 T = T - 273.15;
108
109 // temperature coefficient F
110 if(T>-6.0) F = 0.7 + ((T/-6.0) * 0.3);
111 if(T<=-6.0 && T>-22.0) F = 1 - ((T+6.0)/-16.0 * 0.8);
112 if(T<=-22.0 && T>-40.0) F = 0.2 - ((T+22.0)/-18.0 * 0.2);
113
114 // density coefficient F
115 if(d<150.0) H=1.0;
116
117 if(d>=150.0 & d<400.0) H = 1 - ((d-150.0)/250.0);
118
119 // temperature gradient coefficient G
120 if(dT >= 0.16 && dT < 0.25) G = ((dT - 0.16)/0.09) * 0.1;
121 if(dT >= 0.25 && dT < 0.4) G = 0.1 + (((dT - 0.25)/0.15) * 0.57);
122 if(dT >= 0.4 && dT < 0.5) G = 0.67 + (((dT - 0.4)/0.1) * 0.23);
123 if(dT >= 0.5 && dT < 0.7) G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
124 if(dT >= .7 ) G = 1.0;
125
126 // grouped coefficient Q
127 return F*H*G*E;
128
129} /*}}}*/
130void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/
131
132 /*Created by: Alex S. Gardner, University of Alberta
133
134 *Description*: models the effective snow grain size
135
136 *Reference:*
137 DENDRITIC SNOW METAMORPHISM:
138 Brun, E., P. David, M. Sudul, and G. Brunot, 1992: A numerical model to
139 simulate snow-cover stratigraphy for operational avalanche forecasting.
140 Journal of Glaciology, 38, 13-22.
141
142 NONDENDRITIC SNOW METAMORPHISM:
143 Dry snow metamorphism:
144 Marbouty, D., 1980: An experimental study of temperature-gradient
145 metamorphism. Journal of Glaciology, 26, 303-312.
146
147 WET SNOW METAMORPHISM:
148 Brun, E., 1989: Investigation on wet-snow metamorphism in respect of
149 liquid-water content. Annals of Glaciology, 13, 22-26.
150
151 INPUTS
152 * T: grid cell temperature [k]
153 * dz: grid cell depth [m]
154 * d: grid cell density [kg m-3]
155 * W: water content [kg/m^2]
156 * re: effective grain size [mm]
157 * gdn: grain dentricity
158 * gsp: grain sphericity
159 * dt: time step of input data [s]
160
161 OUTPUTS
162 * re: effective grain size [mm]
163 * gdn: grain dentricity
164 * gsp: grain sphericity*/
165
166 /*intermediary: */
167 IssmDouble dt;
168 IssmDouble* gsz=NULL;
169 IssmDouble* dT=NULL;
170 IssmDouble* zGPC=NULL;
171 IssmDouble* lwc=NULL;
172 IssmDouble Q=0;
173
174 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" grain growth module\n");
175
176 /*only when aIdx = 1 or 2 do we run grainGrowth: */
177 if(aIdx!=1 & aIdx!=2){
178 /*come out as we came in: */
179 return;
180 }
181
182 /*Figure out grain size from effective grain radius: */
183 gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2.0;
184
185 /*Convert dt from seconds to day: */
186 dt=smb_dt/86400.0;
187
188 /*Determine liquid-water content in percentage: */
189 lwc=xNew<IssmDouble>(m); for(int i=0;i<m;i++)lwc[i]= W[i] / (d[i]*dz[i])*100.0;
190
191 //set maximum water content by mass to 9 percent (Brun, 1980)
192 for(int i=0;i<m;i++)if(lwc[i]>9.0) lwc[i]=9.0;
193
194
195 /* Calculate temperature gradiant across grid cells
196 * Returns the avereage gradient across the upper and lower grid cell */
197
198 //initialize
199 dT=xNewZeroInit<IssmDouble>(m);
200
201 //depth of grid point center from surface
202 zGPC=xNewZeroInit<IssmDouble>(m);
203 for(int i=0;i<m;i++){
204 for (int j=0;j<=i;j++) zGPC[i]+=dz[j];
205 zGPC[i]-=dz[i]/2;
206 }
207
208 // Take forward differences on left and right edges
209 dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
210 dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
211
212 //Take centered differences on interior points
213 for(int i=1;i<m-1;i++) dT[i] = (T[i+1]-T[i-1])/(zGPC[i+1]-zGPC[i-1]);
214
215 // take absolute value of temperature gradient
216 for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
217
218 /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
219 for(int i=0;i<m;i++){
220 if (gdn[i]>0){
221
222 //_printf_("Dendritic dry snow metamorphism\n");
223
224 //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
225 if(dT[i]<5.0){
226 //determine coefficients
227 IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
228 IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt;
229 //new dentricity and sphericity for dT < 5 degC m-1
230 gdn[i] += A;
231 gsp[i] += B;
232 }
233 else{
234 // new dendricity and sphericity for dT >= 5 degC m-1
235
236 //determine coefficients
237 IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
238 gdn[i] +=C;
239 gsp[i] +=C;
240 }
241
242 // wet snow metamorphism
243 if(W[i]>0.0){
244
245 //_printf_("D}ritic wet snow metamorphism\n");
246
247 //determine coefficient
248 IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt;
249
250 // new dendricity and sphericity for wet snow
251 gdn[i] -= D;
252 gsp[i] += D;
253 }
254
255 // dendricity and sphericity can not be > 1 or < 0
256 if (gdn[i]<0.0)gdn[i]=0.0;
257 if (gsp[i]<0.0)gsp[i]=0.0;
258 if (gdn[i]>1.0)gdn[i]=1.0;
259 if (gsp[i]>1.0)gsp[i]=1.0;
260
261 // determine new grain size (mm)
262 gsz[i] = 0.1 + (1-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
263
264 }
265 else{
266
267 /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficinets
268 *from Marbouty, 1980: Figure 9*/
269
270 //_printf_("Nond}ritic snow metamorphism\n");
271 Q = Marbouty(T[i],d[i],dT[i]);
272
273 // calculate grain growth
274 gsz[i] += Q* dt;
275
276 //Wet snow metamorphism (Brun, 1989)
277 if(W[i]>0.0){
278 //_printf_("Nond}ritic wet snow metamorphism\n");
279 //wet rate of change coefficient
280 IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400.0); // [mm^3 s^-1]
281
282 // calculate change in grain volume and convert to grain size
283 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);
284
285 }
286
287 // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
288 if (gsp[i]==1.0 && gsz[i]>2.0) gsz[i]=2.0;
289
290 // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
291 if (gsp[i]!=1.0 && gsz[i]>5.0) gsz[i]=5.0;
292 }
293
294 //convert grain size back to effective grain radius:
295 re[i]=gsz[i]/2.0;
296 }
297
298 /*Free ressources:*/
299 xDelete<IssmDouble>(gsz);
300 xDelete<IssmDouble>(dT);
301 xDelete<IssmDouble>(zGPC);
302 xDelete<IssmDouble>(lwc);
303
304} /*}}}*/
305void albedo(IssmDouble* a, 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, int m,int sid) { /*{{{*/
306
307 //// Calculates Snow, firn and ice albedo as a function of:
308 // 1 : effective grain radius (Gardner & Sharp, 2009)
309 // 2 : effective grain radius (Brun et al., 2009)
310 // 3 : density and cloud amount (Greuell & Konzelmann, 1994)
311 // 4 : exponential time decay & wetness (Bougamont & Bamber, 2005)
312
313 //// Inputs
314 // aIdx = albedo method to use
315
316 // Methods 1 & 2
317 // re = surface effective grain radius [mm]
318
319 // Methods 3
320 // d = snow surface density [kg m-3]
321 // n = cloud amount
322 // aIce = albedo of ice
323 // aSnow = albedo of fresh snow
324
325 // Methods 4
326 // aIce = albedo of ice
327 // aSnow = albedo of fresh snow
328 // a = grid cell albedo from prevous time step;
329 // T = grid cell temperature [k]
330 // W = pore water [kg]
331 // P = precipitation [mm w.e.] or [kg m-3]
332 // EC = surface evaporation (-) condensation (+) [kg m-2]
333 // t0wet = time scale for wet snow (15-21.9) [d]
334 // t0dry = warm snow timescale [15] [d]
335 // K = time scale temperature coef. (7) [d]
336 // dt = time step of input data [s]
337
338 //// Output
339 // a = grid cell albedo
340
341 //// Usage
342 // Method 1
343 // a = albedo(1, 0.1);
344
345 // Method 4
346 // a = albedo(4, [], [], [], 0.48, 0.85, [0.8 0.5 ... 0.48], ...
347 // [273 272.5 ... 265], [0 0.001 ... 0], 0, 0.01, 15, 15, 7, 3600)
348
349 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" albedo module\n");
350
351 //some constants:
352 const IssmDouble dSnow = 300; // density of fresh snow [kg m-3]
353 const IssmDouble dIce = 910; // density of ice [kg m-3]
354
355 if(aIdx==1){
356 //function of effective grain radius
357
358 //convert effective radius to specific surface area [cm2 g-1]
359 IssmDouble S = 3.0 / (.091 * re[0]);
360
361 //determine broadband albedo
362 a[0]= 1.48 - pow(S,-.07);
363 }
364 else if(aIdx==2){
365
366 // Spectral fractions (Lefebre et al., 2003)
367 // [0.3-0.8um 0.8-1.5um 1.5-2.8um]
368
369 IssmDouble sF[3] = {0.606, 0.301, 0.093};
370
371 // convert effective radius to grain size in meters
372 IssmDouble gsz = (re[0] * 2) / 1000.0;
373
374 // spectral range:
375 // 0.3 - 0.8um
376 IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5));
377 // 0.8 - 1.5um
378 IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5));
379 // 1.5 - 2.8um
380 IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
381
382 // broadband surface albedo
383 a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
384
385 }
386 else if(aIdx==3){
387
388 // a as a function of density
389
390 // calculate albedo
391 a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
392 }
393 else if(aIdx==4){
394
395 // exponential time decay & wetness
396
397 // change in albedo with time:
398 // (d_a) = (a - a_old)/(t0)
399 // where: t0 = timescale for albedo decay
400
401 dt = dt / 86400; // convert from [s] to [d]
402
403 // initialize variables
404 IssmDouble* t0=xNew<IssmDouble>(m);
405 IssmDouble* T=xNew<IssmDouble>(m);
406 IssmDouble* t0warm=xNew<IssmDouble>(m);
407 IssmDouble* d_a=xNew<IssmDouble>(m);
408
409 // specify constants
410 // a_wet = 0.15; // water albedo (0.15)
411 // a_new = aSnow // new snow albedo (0.64 - 0.89)
412 // a_old = aIce; // old snow/ice albedo (0.27-0.53)
413 // t0_wet = t0wet; // time scale for wet snow (15-21.9) [d]
414 // t0_dry = t0dry; // warm snow timescale [15] [d]
415 // K = 7 // time scale temperature coef. (7) [d]
416 // W0 = 300; // 200 - 600 [mm]
417 const IssmDouble z_snow = 15; // 16 - 32 [mm]
418
419 // determine timescale for albedo decay
420 for(int i=0;i<m;i++)if(W[i]>0)t0[i]=t0wet; // wet snow timescale
421 for(int i=0;i<m;i++)T[i]=TK[i] - 273.15; // change T from K to degC
422 for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
423 for(int i=0;i<m;i++)if(W[i]==0.0 && T[i]>=-10)t0[i]= t0warm[i];
424 for(int i=0;i<m;i++)if(T[i]<-10) t0[i] = 10 * K + t0dry; // 'cold' snow timescale
425
426 // calculate new albedo
427 for(int i=0;i<m;i++)d_a[i] = (a[i] - aIce) / t0[i] * dt; // change in albedo
428 for(int i=0;i<m;i++)a[i] -= d_a[i]; // new albedo
429
430 // modification of albedo due to thin layer of snow or solid
431 // condensation (deposition) at the surface surface
432
433 // check if condensation occurs & if it is deposited in solid phase
434 if ( EC > 0 && T[0] < 0) P = P + (EC/dSnow) * 1000; // add cond to precip [mm]
435
436 a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);
437
438 //----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------
439 // modification of albedo due to thin layer of water on the surface
440 // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);
441
442 /*Free ressources:*/
443 xDelete<IssmDouble>(t0);
444 xDelete<IssmDouble>(T);
445 xDelete<IssmDouble>(t0warm);
446 xDelete<IssmDouble>(d_a);
447
448 }
449 else _error_("albedo method switch should range from 1 to 4!");
450
451 // Check for erroneous values
452 if (a[0] > 1) _printf_("albedo > 1.0\n");
453 else if (a[0] < 0) _printf_("albedo is negative\n");
454 else if (xIsNan(a[0])) _error_("albedo == NAN\n");
455} /*}}}*/
456void thermo(IssmDouble* pEC, IssmDouble* T, 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,int sid) { /*{{{*/
457
458 /* ENGLACIAL THERMODYNAMICS*/
459
460 /* Description:
461 computes new temperature profile accounting for energy absorption and
462 thermal diffusion.*/
463
464 // INPUTS
465 // T: grid cell temperature [k]
466 // dz: grid cell depth [m]
467 // d: grid cell density [kg m-3]
468 // swf: shortwave radiation fluxes [W m-2]
469 // dlwrf: downward longwave radiation fluxes [W m-2]
470 // Ta: 2 m air temperature
471 // V: wind velocity [m s-1]
472 // eAir: screen level vapor pressure [Pa]
473 // Ws: surface water content [kg]
474 // dt0: time step of input data [s]
475 // elev: surface elevation [m a.s.l.]
476 // Vz: air temperature height above surface [m]
477 // Tz: wind height above surface [m]
478
479 // OUTPUTS
480 // T: grid cell temperature [k]
481 // EC: evaporation/condensation [kg]
482
483 /*intermediary: */
484 IssmDouble* K = NULL;
485 IssmDouble* KU = NULL;
486 IssmDouble* KD = NULL;
487 IssmDouble* KP = NULL;
488 IssmDouble* Au = NULL;
489 IssmDouble* Ad = NULL;
490 IssmDouble* Ap = NULL;
491 IssmDouble* Nu = NULL;
492 IssmDouble* Nd = NULL;
493 IssmDouble* Np = NULL;
494 IssmDouble* dzU = NULL;
495 IssmDouble* dzD = NULL;
496 IssmDouble* sw = NULL;
497 IssmDouble* dT_sw = NULL;
498 IssmDouble* lw = NULL;
499 IssmDouble* T0 = NULL;
500 IssmDouble* Tu = NULL;
501 IssmDouble* Td = NULL;
502
503 IssmDouble z0;
504 IssmDouble dt;
505 IssmDouble max_fdt=0;
506 IssmDouble Ts=0;
507 IssmDouble L;
508 IssmDouble eS;
509 IssmDouble Ri=0;
510 IssmDouble coefM;
511 IssmDouble coefH;
512 IssmDouble An;
513 IssmDouble C;
514 IssmDouble shf;
515 IssmDouble SB;
516 IssmDouble CI;
517 IssmDouble ds;
518 IssmDouble dAir;
519 IssmDouble TCs;
520 IssmDouble lhf;
521 IssmDouble EC_day;
522 IssmDouble dT_turb;
523 IssmDouble turb;
524 IssmDouble ulw;
525 IssmDouble dT_ulw;
526 IssmDouble dlw;
527 IssmDouble dT_dlw;
528
529 /*outputs:*/
530 IssmDouble EC;
531
532 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" thermal module\n");
533
534 // INITIALIZE
535 CI = 2102; // heat capacity of snow/ice (J kg-1 k-1)
536 // CA = 1005; // heat capacity of air (J kg-1 k-1)
537 // LF = 0.3345E6; // latent heat of fusion(J kg-1)
538 // LV = 2.495E6; // latent heat of vaporization(J kg-1)
539 // dIce = 910; // density of ice [kg m-3]
540 // dSnow = 300; // density of snow [kg m-3]
541 SB = 5.67E-8; // Stefan-Boltzmann constant [W m-2 K-4]
542
543 ds = d[0]; // density of top grid cell
544
545 // calculated air density [kg/m3]
546 dAir = 0.029 * pAir /(8.314 * Ta);
547
548 // thermal capacity of top grid cell [J/k]
549 TCs = d[0]*dz[0]*CI;
550
551 //initialize Evaporation - Condenstation
552 EC = 0;
553
554 // check if all SW applied to surface or distributed throught subsurface
555 // swIdx = length(swf) > 1
556
557 // SURFACE ROUGHNESS (Bougamont, 2006)
558 // wind/temperature surface roughness height [m]
559 if (ds < 910 && Ws == 0) z0 = 0.00012; // 0.12 mm for dry snow
560 else if (ds >= 910) z0 = 0.0032; // 3.2 mm for ice
561 else z0 = 0.0013; // 1.3 mm for wet snow
562
563 // if V = 0 goes to infinity therfore if V = 0 change
564 if(V<.01)V=.01;
565
566 // Bulk-transfer coefficient for turbulent fluxes
567 An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient
568 C = An * dAir * V; // shf & lhf common coefficient
569
570 // THERMAL CONDUCTIVITY (Sturm, 1997: J. Glaciology)
571 // calculate new K profile [W m-1 K-1]
572
573 // initialize conductivity
574 K= xNewZeroInit<IssmDouble>(m);
575
576 // for snow and firn (density < 910 kg m-3) (Sturn et al, 1997)
577 for(int i=0;i<m;i++) if(d[i]<910) K[i] = 0.138 - 1.01E-3 * d[i] + 3.233E-6 * (pow(d[i],2));
578
579 // for ice (density >= 910 kg m-3)
580 for(int i=0;i<m;i++) if(d[i]>=910) K[i] = 9.828 * exp(-5.7E-3*T[i]);
581
582 // THERMAL DIFFUSION COEFFICIENTS
583
584 // A discretization scheme which truncates the Taylor-Series expansion
585 // after the 3rd term is used. See Patankar 1980, Ch. 3&4
586
587 // discretized heat equation:
588
589 // Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap
590
591 // where neighbor coefficients Au, Ap, & Ad are
592
593 // Au = [dz_u/2KP + dz_p/2KE]^-1
594 // Ad = [dz_d/2KP + dz_d/2KD]^-1
595 // Ap = d*CI*dz/Dt
596
597 // and u & d represent grid points up and down from the center grid point
598 // p and // u & d represent grid points up and down from the center grid
599 // point p and ° identifies previous time step values. S is a source term.
600
601 // u, d, and p conductivities
602 KU = xNew<IssmDouble>(m);
603 KD = xNew<IssmDouble>(m);
604 KP = xNew<IssmDouble>(m);
605
606 KU[0] = UNDEF;
607 KD[m-1] = UNDEF;
608 for(int i=1;i<m;i++) KU[i]= K[i-1];
609 for(int i=0;i<m-1;i++) KD[i] = K[i+1];
610 for(int i=0;i<m;i++) KP[i] = K[i];
611
612 // determine u, d & p cell widths
613 dzU = xNew<IssmDouble>(m);
614 dzD = xNew<IssmDouble>(m);
615 dzU[0]=UNDEF;
616 dzD[m-1]=UNDEF;
617
618 for(int i=1;i<m;i++) dzU[i]= dz[i-1];
619 for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
620
621 // determine minimum acceptable delta t (diffusion number > 1/2) [s]
622 dt=1e12;
623 for(int i=0;i<m;i++)dt = fmin(dt,CI * pow(dz[i],2) * d[i] / (3 * K[i]));
624
625 // smallest possible even integer of 60 min where diffusion number > 1/2
626 // must go evenly into one hour or the data frequency if it is smaller
627
628 // all integer factors of the number of second in a day (86400 [s])
629 int f[45] = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 30, 36, 40, 45, 48, 50, 60,
630 72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600};
631
632 // return the min integer factor that is < dt
633 max_fdt=f[0];
634 for(int i=0;i<45;i++){
635 if (f[i]<dt)if(f[i]>=max_fdt)max_fdt=f[i];
636 }
637 dt=max_fdt;
638
639 // determine mean (harmonic mean) of K/dz for u, d, & p
640 Au = xNew<IssmDouble>(m);
641 Ad = xNew<IssmDouble>(m);
642 Ap = xNew<IssmDouble>(m);
643 for(int i=0;i<m;i++){
644 Au[i] = pow((dzU[i]/2/KP[i] + dz[i]/2/KU[i]),-1);
645 Ad[i] = pow((dzD[i]/2/KP[i] + dz[i]/2/KD[i]),-1);
646 Ap[i] = (d[i]*dz[i]*CI)/dt;
647 }
648
649 // create "neighbor" coefficient matrix
650 Nu = xNew<IssmDouble>(m);
651 Nd = xNew<IssmDouble>(m);
652 Np = xNew<IssmDouble>(m);
653 for(int i=0;i<m;i++){
654 Nu[i] = Au[i] / Ap[i];
655 Nd[i] = Ad[i] / Ap[i];
656 Np[i]= 1 - Nu[i] - Nd[i];
657 }
658
659 // specify boundary conditions: constant flux at bottom
660 Nu[m-1] = 0;
661 Np[m-1] = 1;
662
663 // zero flux at surface
664 Np[0] = 1 - Nd[0];
665
666 // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
667 Nu[0] = 0;
668 Nd[m-1] = 0;
669
670 /* RADIATIVE FLUXES*/
671
672 // energy supplied by shortwave radiation [J]
673 sw = xNew<IssmDouble>(m);
674 for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
675
676 // temperature change due to SW
677 dT_sw = xNew<IssmDouble>(m);
678 for(int i=0;i<m;i++) dT_sw[i]= sw[i] / (CI * d[i] * dz[i]);
679
680 // Upward longwave radiation flux is calculated from the snow surface
681 // temperature which is set equal to the average temperature of the
682 // top grid cells.
683
684 // energy supplied by downward longwave radiation to the top grid cell [J]
685 dlw = dlwrf * dt;
686
687 // temperature change due to dlw_surf
688 dT_dlw = dlw / TCs;
689
690 // PREALLOCATE ARRAYS BEFORE LOOP FOR IMPROVED PERFORMANCE
691 T0 = xNewZeroInit<IssmDouble>(m+2);
692 Tu=xNew<IssmDouble>(m);
693 Td=xNew<IssmDouble>(m);
694
695 /* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
696 for (IssmDouble i=1;i<=dt0;i+=dt){
697
698 // PART OF ENERGY CONSERVATION CHECK
699 // store initial temperature
700 //T_init = T;
701
702 // calculate temperature of snow surface (Ts)
703 // when incoming SW radition is allowed to penetrate the surface,
704 // the modeled energy balance becomes very sensitive to how Ts is
705 // calculated. The estimated enegy balance & melt are significanly
706 // less when Ts is taken as the mean of the x top grid cells.
707 Ts = (T[0] + T[1])/2;
708 Ts = fmin(273.15,Ts); // don't allow Ts to exceed 273.15 K (0°C)
709
710 //TURBULENT HEAT FLUX
711
712 // Monin–Obukhov Stability Correction
713 // Reference:
714 // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
715 // Journal of Climatology, 2, 65-84.
716
717 // calculate the Bulk Richardson Number (Ri)
718 Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
719
720 // calculate Monin-–Obukhov stability factors 'coefM' and 'coefH'
721
722 // do not allow Ri to exceed 0.19
723 Ri = fmin(Ri, 0.19);
724
725 // calculate momentum 'coefM' stability factor
726 if (Ri > 0){
727 // if stable
728 coefM = 1/(1-5.2*Ri);
729 }
730 else {
731 coefM =pow (1-18*Ri,-0.25);
732 }
733
734 // calculate heat/wind 'coef_H' stability factor
735 if (Ri < -0.03) coefH = 1.3 * coefM;
736 else coefH = coefM;
737
738 //// Sensible Heat
739 // calculate the sensible heat flux [W m-2](Patterson, 1998)
740 shf = C * 1005 * (Ta - Ts);
741
742 // adjust using Monin–Obukhov stability theory
743 shf = shf / (coefM * coefH);
744
745 //// Latent Heat
746 // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
747 if (Ts >= 273.15){
748 L = 2.495E6;
749
750 // for an ice surface Murphy and Koop, 2005 [Equation 7]
751 eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
752 }
753 else{
754 L = 2.8295E6; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C)
755 // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
756 eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15));
757 }
758
759 // Latent heat flux [W m-2]
760 lhf = C * L * (eAir - eS) * 0.622 / pAir;
761
762 // adjust using Monin–Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface,
763 // if '-' then there is mass and energy loss at the surface.
764 lhf = lhf / (coefM * coefH);
765
766 //mass loss (-)/acreation(+) due to evaporation/condensation [kg]
767 EC_day = lhf * 86400 / L;
768
769 // temperature change due turbulent fluxes
770 turb = (shf + lhf)* dt;
771 dT_turb = turb / TCs;
772
773 // upward longwave contribution
774 ulw = - SB * pow(Ts,4.0) * dt;
775 dT_ulw = ulw / TCs;
776
777 // new grid point temperature
778
779 //SW penetrates surface
780 for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
781 T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
782
783 // temperature diffusion
784 for(int j=0;j<m;j++)T0[1+j]=T[j];
785 for(int j=0;j<m;j++) Tu[j] = T0[j];
786 for(int j=0;j<m;j++) Td[j] = T0[2+j];
787 for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
788
789 // calculate cumulative evaporation (+)/condensation(-)
790 EC = EC + (EC_day/86400)*dt;
791
792 /* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
793 //energy flux across lower boundary (energy supplied by underling ice)
794 base_flux = Ad(-1)*(T_init()-T_init(-1)) * dt;
795
796 E_used = sum((T - T_init) * (d*dz*CI));
797 E_sup = ((sum(swf) * dt) + dlw + ulw + turb + base_flux);
798
799 E_diff = E_used - E_sup;
800
801 if abs(E_diff) > 1E-6 || isnan(E_diff)
802 disp(T(1))
803 _error_("energy not conserved in thermodynamics equations");
804 */
805 }
806
807 /*Free ressources:*/
808 xDelete<IssmDouble>(K);
809 xDelete<IssmDouble>(KU);
810 xDelete<IssmDouble>(KD);
811 xDelete<IssmDouble>(KP);
812 xDelete<IssmDouble>(Au);
813 xDelete<IssmDouble>(Ad);
814 xDelete<IssmDouble>(Ap);
815 xDelete<IssmDouble>(Nu);
816 xDelete<IssmDouble>(Nd);
817 xDelete<IssmDouble>(Np);
818 xDelete<IssmDouble>(dzU);
819 xDelete<IssmDouble>(dzD);
820 xDelete<IssmDouble>(sw);
821 xDelete<IssmDouble>(dT_sw);
822 xDelete<IssmDouble>(lw);
823 xDelete<IssmDouble>(T0);
824 xDelete<IssmDouble>(Tu);
825 xDelete<IssmDouble>(Td);
826
827
828 /*Assign output pointers:*/
829 *pEC=EC;
830
831} /*}}}*/
832void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
833
834 // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
835
836 // swIdx = 0 : all absorbed SW energy is assigned to the top grid cell
837
838 // swIdx = 1 : absorbed SW is distributed with depth as a function of:
839 // 1 : snow density (taken from Bassford, 2004)
840 // 2 : grain size in 3 spectral bands (Brun et al., 1992)
841
842 // Inputs
843 // swIdx = shortwave allowed to penetrate surface (0 = No, 1 = Yes)
844 // aIdx = method for calculating albedo (1-4)
845 // dsw = downward shortwave radiative flux [w m-2]
846 // as = surface albedo
847 // d = grid cell density [kg m-3]
848 // dz = grid cell depth [m]
849 // re = grid cell effective grain radius [mm]
850
851 // Outputs
852 // swf = absorbed shortwave radiation [W m-2]
853 //
854
855 /*outputs: */
856 IssmDouble* swf=NULL;
857
858 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" shortwave module\n");
859
860 /*Initialize and allocate: */
861 swf=xNewZeroInit<IssmDouble>(m);
862
863
864 // SHORTWAVE FUNCTION
865 if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
866
867 // calculate surface shortwave radiation fluxes [W m-2]
868 swf[0] = (1 - as) * dsw;
869 }
870 else{ // sw radation is absorbed at depth within the glacier
871
872 if (aIdx == 2){ // function of effective radius (3 spectral bands)
873
874 IssmDouble * gsz=NULL;
875 IssmDouble * B1_cum=NULL;
876 IssmDouble * B2_cum=NULL;
877 IssmDouble* h =NULL;
878 IssmDouble* B1 =NULL;
879 IssmDouble* B2 =NULL;
880 IssmDouble* exp1 = NULL;
881 IssmDouble* exp2 = NULL;
882 IssmDouble* Qs1 = NULL;
883 IssmDouble* Qs2 = NULL;
884
885 // convert effective radius [mm] to grain size [m]
886 gsz=xNew<IssmDouble>(m);
887 for(int i=0;i<m;i++) gsz[i]= (re[i] * 2) / 1000;
888
889 // Spectral fractions [0.3-0.8um 0.8-1.5um 1.5-2.8um]
890 // (Lefebre et al., 2003)
891 IssmDouble sF[3] = {0.606, 0.301, 0.093};
892
893 // initialize variables
894 B1_cum=xNew<IssmDouble>(m+1);
895 B2_cum=xNew<IssmDouble>(m+1);
896 for(int i=0;i<m+1;i++){
897 B1_cum[i]=1;
898 B2_cum[i]=1;
899 }
900
901
902 // spectral albedos:
903 // 0.3 - 0.8um
904 IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz[0],0.5));
905 // 0.8 - 1.5um
906 IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz[0],0.5));
907 // 1.5 - 2.8um
908 IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
909
910 // separate net shortwave radiative flux into spectral ranges
911 IssmDouble swfS[3];
912 swfS[0] = (sF[0] * dsw) * (1 - a0);
913 swfS[1] = (sF[1] * dsw) * (1 - a1);
914 swfS[2] = (sF[2] * dsw) * (1 - a2);
915
916 // absorption coeficient for spectral range:
917 h =xNew<IssmDouble>(m);
918 B1 =xNew<IssmDouble>(m);
919 B2 =xNew<IssmDouble>(m);
920 for(int i=0;i<m;i++) h[i]= d[i] /(pow(gsz[i],0.5));
921 for(int i=0;i<m;i++) B1[i] = .0192 * h[i]; // 0.3 - 0.8um
922 for(int i=0;i<m;i++) B2[i]= .1098 * h[i]; // 0.8 - 1.5um
923 // B3 = +inf // 1.5 - 2.8um
924
925 // cumulative extinction factors
926 exp1 = xNew<IssmDouble>(m);
927 exp2 = xNew<IssmDouble>(m);
928 for(int i=0;i<m;i++) exp1[i]=exp(-B1[i]*dz[i]);
929 for(int i=0;i<m;i++) exp2[i]=exp(-B2[i]*dz[i]);
930
931 for(int i=0;i<m;i++){
932 IssmDouble cum1=exp1[0];
933 IssmDouble cum2=exp2[0];
934 for(int j=1;j<=i;j++){
935 cum1 = cum1*exp1[j];
936 cum2 = cum2*exp2[j];
937 }
938 B1_cum[i+1]=cum1;
939 B2_cum[i+1]=cum2;
940 }
941
942
943 // flux across grid cell boundaries
944 Qs1 = xNew<IssmDouble>(m+1);
945 Qs2 = xNew<IssmDouble>(m+1);
946 for(int i=0;i<m+1;i++){
947 Qs1[i] = swfS[0] * B1_cum[i];
948 Qs2[i] = swfS[1] * B2_cum[i];
949 }
950
951 // net energy flux to each grid cell
952 for(int i=0;i<m;i++) swf[i]= (Qs1[i]-Qs1[i+1]) + (Qs2[i]-Qs2[i+1]);
953
954 // add flux absorbed at surface
955 swf[0] = swf[0]+ swfS[2];
956
957 /*Free ressources: */
958 xDelete<IssmDouble>(gsz);
959 xDelete<IssmDouble>(B1_cum);
960 xDelete<IssmDouble>(B2_cum);
961 xDelete<IssmDouble>(h);
962 xDelete<IssmDouble>(B1);
963 xDelete<IssmDouble>(B2);
964 xDelete<IssmDouble>(exp1);
965 xDelete<IssmDouble>(exp2);
966 xDelete<IssmDouble>(Qs1);
967 xDelete<IssmDouble>(Qs2);
968
969
970 }
971 else{ //function of grid cell density
972
973 /*intermediary: */
974 IssmDouble* B_cum = NULL;
975 IssmDouble* exp_B = NULL;
976 IssmDouble* Qs = NULL;
977 IssmDouble* B = NULL;
978
979 // fraction of sw radiation absorbed in top grid cell (wavelength > 0.8um)
980 IssmDouble SWs = 0.36;
981
982 // SWs and SWss coefficients need to be better constranted. Greuell
983 // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
984 // the // of SW radiation with wavelengths > and < 800 nm
985 // respectively. This, however, may not account for the fact that
986 // the albedo of wavelengths > 800 nm has a much lower albedo.
987
988 // calculate surface shortwave radiation fluxes [W m-2]
989 IssmDouble swf_s = SWs * (1 - as) * dsw;
990
991 // calculate surface shortwave radiation fluxes [W m-2]
992 IssmDouble swf_ss = (1-SWs) * (1 - as) * dsw;
993
994 // SW allowed to penetrate into snowpack
995 IssmDouble Bs = 10; // snow SW extinction coefficient [m-1] (Bassford,2006)
996 IssmDouble Bi = 1.3; // ice SW extinction coefficient [m-1] (Bassford,2006)
997
998 // calculate extinction coefficient B [m-1] vector
999 B=xNew<IssmDouble>(m);
1000 for(int i=0;i<m;i++) B[i] = Bs + (300 - d[i]) * ((Bs - Bi)/(910 - 300));
1001
1002 // cumulative extinction factor
1003 B_cum = xNew<IssmDouble>(m+1);
1004 exp_B = xNew<IssmDouble>(m);
1005 for(int i=0;i<m;i++)exp_B[i]=exp(-B[i]*dz[i]);
1006
1007 B_cum[0]=1;
1008 for(int i=0;i<m;i++){
1009 IssmDouble cum_B=exp_B[0];
1010 for(int j=1;j<=i;j++) cum_B=cum_B*exp_B[j];
1011 B_cum[i+1]= cum_B;
1012 }
1013
1014 // flux across grid cell boundaries
1015 Qs=xNew<IssmDouble>(m+1);
1016 for(int i=0;i<m+1;i++) Qs[i] = swf_ss * B_cum[i];
1017
1018 // net energy flux to each grid cell
1019 for(int i=0;i<m;i++) swf[i] = (Qs[i]-Qs[i+1]);
1020
1021 // add flux absorbed at surface
1022 swf[0] += swf_s;
1023
1024 /*Free ressources:*/
1025 xDelete<IssmDouble>(B_cum);
1026 xDelete<IssmDouble>(exp_B);
1027 xDelete<IssmDouble>(Qs);
1028 xDelete<IssmDouble>(B);
1029 }
1030 }
1031 /*Assign output pointers: */
1032 *pswf=swf;
1033
1034} /*}}}*/
1035void 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, int sid){ /*{{{*/
1036
1037 // Adds precipitation and deposition to the model grid
1038
1039 // Author: Alex Gardner, University of Alberta
1040 // Date last modified: JAN, 2008
1041
1042 /* Description:
1043 adjusts the properties of the top grid cell to account for accumulation
1044 T_air & T = Air and top grid cell temperatures [K]
1045 dz = topgrid cell length [m]
1046 d = density of top grid gell [kg m-3]
1047 P = precipitation [mm w.e.] or [kg m-3]
1048 re = effective grain radius [mm]
1049 gdn = grain dentricity
1050 gsp = grain sphericity*/
1051
1052 // MAIN FUNCTION
1053 // specify constants
1054 const IssmDouble dIce = 910; // density of ice [kg m-3]
1055 const IssmDouble dSnow = 150; // density of snow [kg m-3]
1056 const IssmDouble reNew = 0.1; // new snow grain size [mm]
1057 const IssmDouble gdnNew = 1; // new snow dendricity
1058 const IssmDouble gspNew = 0.5; // new snow sphericity
1059
1060 /*intermediary: */
1061 IssmDouble* mInit=NULL;
1062 bool top=true;
1063 IssmDouble mass, massinit, mass_diff;
1064
1065 /*output: */
1066 IssmDouble* T=NULL;
1067 IssmDouble* dz=NULL;
1068 IssmDouble* d=NULL;
1069 IssmDouble* W=NULL;
1070 IssmDouble* a=NULL;
1071 IssmDouble* re=NULL;
1072 IssmDouble* gdn=NULL;
1073 IssmDouble* gsp=NULL;
1074 int m;
1075
1076 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" accumulation module\n");
1077
1078 /*Recover pointers: */
1079 T=*pT;
1080 dz=*pdz;
1081 d=*pd;
1082 W=*pW;
1083 a=*pa;
1084 re=*pre;
1085 gdn=*pgdn;
1086 gsp=*pgsp;
1087 m=*pm;
1088
1089 // determine initial mass
1090 mInit=xNew<IssmDouble>(m);
1091 for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
1092 massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
1093
1094 if (P > 0){
1095
1096
1097 if (T_air <= 273.15){ // if snow
1098
1099 IssmDouble z_snow = P/dSnow; // depth of snow
1100
1101 // if snow depth is greater than specified min dz, new cell created
1102 if (z_snow > dzMin){
1103
1104 newcell(&T,T_air,top,m); //new cell T
1105 newcell(&dz,z_snow,top,m); //new cell dz
1106 newcell(&d,dSnow,top,m); //new cell d
1107 newcell(&W,0,top,m); //new cell W
1108 newcell(&a,aSnow,top,m); //new cell a
1109 newcell(&re,reNew,top,m); //new cell grain size
1110 newcell(&gdn,gdnNew,top,m); //new cell grain dendricity
1111 newcell(&gsp,gspNew,top,m); //new cell grain sphericity
1112 m=m+1;
1113 }
1114 else { // if snow depth is less than specified minimum dz snow
1115
1116 IssmDouble mass = mInit[0] + P; // grid cell adjust mass
1117
1118 dz[0] = dz[0] + P/dSnow; // adjust grid cell depth
1119 d[0] = mass / dz[0]; // adjust grid cell density
1120
1121 // adjust variables as a linearly weighted function of mass
1122 // adjust temperature (assume P is same temp as air)
1123 T[0] = (T_air * P + T[0] * mInit[0])/mass;
1124
1125 // adjust a, re, gdn & gsp
1126 a[0] = (aSnow * P + a[0] * mInit[0])/mass;
1127 re[0] = (reNew * P + re[0] * mInit[0])/mass;
1128 gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass;
1129 gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass;
1130 }
1131 }
1132 else{ // if rain
1133
1134 /*rain is added by increasing the mass and temperature of the ice
1135 of the top grid cell. Temperatures are set artifically high to
1136 account for the latent heat of fusion. This is the same as
1137 directly adding liquid water to the the snow pack surface but
1138 makes the numerics easier.*/
1139
1140 IssmDouble LF = 0.3345E6; // latent heat of fusion(J kg-1)
1141 IssmDouble CI = 2102; // specific heat capacity of snow/ice (J kg-1 k-1)
1142
1143 // grid cell adjust mass
1144 mass = mInit[0] + P;
1145
1146 // adjust temperature
1147 // liquid: must account for latent heat of fusion
1148 T[0] = (P *(T_air + LF/CI) + T[0] * mInit[0]) / mass;
1149
1150 // adjust grid cell density
1151 d[0] = mass / dz[0];
1152
1153 // if d > the density of ice, d = dIce
1154 if (d[0] > dIce){
1155 d[0] = dIce; // adjust d
1156 dz[0] = mass / d[0]; // dz is adjusted to conserve mass
1157 }
1158 }
1159
1160 // check for conservation of mass
1161 mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i];
1162
1163 mass_diff = mass - massinit - P;
1164
1165 #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode.
1166 mass_diff = round(mass_diff * 100)/100;
1167 if (mass_diff > 0) _error_("mass not conserved in accumulation function");
1168 #endif
1169
1170 }
1171 /*Free ressources:*/
1172 if(mInit)xDelete<IssmDouble>(mInit);
1173
1174 /*Assign output pointers:*/
1175 *pT=T;
1176 *pdz=dz;
1177 *pd=d;
1178 *pW=W;
1179 *pa=a;
1180 *pre=re;
1181 *pgdn=gdn;
1182 *pgsp=gsp;
1183 *pm=m;
1184} /*}}}*/
1185void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, int sid){ /*{{{*/
1186
1187 //// MELT ROUTINE
1188
1189 // Description:
1190 // computes the quantity of meltwater due to snow temperature in excess of
1191 // 0 deg C, determines pore water content and adjusts grid spacing
1192
1193 /*intermediary:*/
1194 IssmDouble* m=NULL;
1195 IssmDouble* maxF=NULL;
1196 IssmDouble* dW=NULL;
1197 IssmDouble* exsW=NULL;
1198 IssmDouble* exsT=NULL;
1199 IssmDouble* surpT=NULL;
1200 IssmDouble* surpE=NULL;
1201 IssmDouble* F=NULL;
1202 IssmDouble* flxDn=NULL;
1203 IssmDouble ER=0;
1204 IssmDouble* EI=NULL;
1205 IssmDouble* EW=NULL;
1206 IssmDouble* M=NULL;
1207 int* D=NULL;
1208
1209 IssmDouble sumM;
1210 IssmDouble sumER;
1211 IssmDouble addE;
1212 IssmDouble mSum0;
1213 IssmDouble sumE0;
1214 IssmDouble mSum1;
1215 IssmDouble sumE1;
1216 IssmDouble dE;
1217 IssmDouble dm;
1218 IssmDouble X;
1219 IssmDouble Wi;
1220 int D_size;
1221 int i;
1222
1223 /*outputs:*/
1224 IssmDouble mAdd;
1225 IssmDouble Rsum;
1226 IssmDouble* T=*pT;
1227 IssmDouble* d=*pd;
1228 IssmDouble* dz=*pdz;
1229 IssmDouble* W=*pW;
1230 IssmDouble* a=*pa;
1231 IssmDouble* re=*pre;
1232 IssmDouble* gdn=*pgdn;
1233 IssmDouble* gsp=*pgsp;
1234 int n=*pn;
1235 IssmDouble* R=0;
1236
1237 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" melt module\n");
1238
1239 //// INITIALIZATION
1240
1241 /*Allocations: */
1242 M=xNewZeroInit<IssmDouble>(n);
1243 maxF=xNew<IssmDouble>(n);
1244 dW=xNew<IssmDouble>(n);
1245
1246 // specify constants
1247 const IssmDouble CtoK = 273.15; // clecius to Kelvin conversion
1248 const IssmDouble CI = 2102; // specific heat capacity of snow/ice (J kg-1 k-1)
1249 const IssmDouble LF = 0.3345E6; // latent heat of fusion(J kg-1)
1250 const IssmDouble dPHC = 830; // pore hole close off density[kg m-3]
1251 const IssmDouble dIce = 910; // density of ice [kg m-3]
1252
1253 // store initial mass [kg] and energy [J]
1254 m=xNew<IssmDouble>(n); for(int i=0;i<n;i++) m[i] = dz[i]* d[i]; // grid cell mass [kg]
1255 EI=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI; // initial enegy of snow/ice
1256 EW=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EW[i]= W[i] * (LF + CtoK * CI); // initial enegy of water
1257
1258 mSum0 = cellsum(W,n) + cellsum(m,n); // total mass [kg]
1259 sumE0 = cellsum(EI,n) + cellsum(EW,n); // total energy [J]
1260
1261 // initialize melt and runoff scalars
1262 Rsum = 0; // runoff [kg]
1263 sumM = 0; // total melt [kg]
1264 mAdd = 0; // mass added/removed to/from base of model [kg]
1265 addE = 0; // energy added/removed to/from base of model [J]
1266
1267 // calculate temperature excess above 0 deg C
1268 exsT=xNewZeroInit<IssmDouble>(n);
1269 for(int i=0;i<n;i++) exsT[i]= fmax(0, T[i] - CtoK); // [K] to [°C]
1270
1271 // new grid point center temperature, T [K]
1272 for(int i=0;i<n;i++) T[i]-=exsT[i];
1273
1274 // specify irreducible water content saturation [fraction]
1275 const IssmDouble Swi = 0.07; // assumed constant after Colbeck, 1974
1276
1277 //// REFREEZE PORE WATER
1278 // check if any pore water
1279 if (cellsum(W,n) > 0){
1280 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" pore water refreeze\n");
1281 // calculate maximum freeze amount, maxF [kg]
1282 for(int i=0;i<n;i++) maxF[i] = fmax(0, -((T[i] - CtoK) * m[i] * CI) / LF);
1283
1284 // freeze pore water and change snow/ice properties
1285 for(int i=0;i<n;i++) dW[i] = fmin(maxF[i], W[i]); // freeze mass [kg]
1286 for(int i=0;i<n;i++) W[i] -= dW[i]; // pore water mass [kg]
1287 for(int i=0;i<n;i++) m[i] += dW[i]; // new mass [kg]
1288 for(int i=0;i<n;i++) d[i] = m[i] / dz[i]; // density [kg m-3]
1289 for(int i=0;i<n;i++) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI)); // temperature [K]
1290
1291 // if pore water froze in ice then adjust d and dz thickness
1292 for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=dIce;
1293 for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
1294 }
1295
1296 // squeeze water from snow pack
1297 exsW=xNew<IssmDouble>(n);
1298 for(int i=0;i<n;i++){
1299 Wi= (910 - d[i]) * Swi * (m[i] / d[i]); // irreducible water content [kg]
1300 exsW[i] = fmax(0, W[i] - Wi); // water "squeezed" from snow [kg]
1301 }
1302
1303 //// MELT, PERCOLATION AND REFREEZE
1304
1305 // run melt algorithm if there is melt water or excess pore water
1306 if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
1307
1308 // _printf_(""MELT OCCURS");
1309 // check to see if thermal energy exceeds energy to melt entire cell
1310 // if so redistribute temperature to lower cells (temperature surplus)
1311 // (maximum T of snow before entire grid cell melts is a constant
1312 // LF/CI = 159.1342)
1313 surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0, exsT [i]- 159.1342);
1314
1315 if (cellsum(surpT,n) > 0 ){
1316 // _printf_("T Surplus");
1317 // calculate surplus energy
1318 surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI / m[i];
1319
1320 int i = 0;
1321 while (cellsum(surpE,n) > 0){
1322 // use surplus energy to increase the temperature of lower cell
1323 T[i+1] = surpE[i] * m[i+1]/CI + T[i+1];
1324 surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342));
1325 surpE[i+1] = surpT[i+1] * CI / m[i+1];
1326
1327 // adjust current cell properties (again 159.1342 is the max T)
1328 T[i] = CtoK + 159.1342;
1329 surpE[i] = 0;
1330 i = i + 1;
1331 }
1332 // recalculate temperature excess above 0 deg C
1333 for(int i=0;i<n;i++) exsT[i] = fmax(0, T[i] - CtoK);
1334 }
1335
1336 // convert temperature excess to melt [kg]
1337 for(int i=0;i<n;i++) M[i] = exsT[i] * d[i] * dz[i] * CI / LF; // melt
1338 sumM = cellsum(M,n); // total melt [kg]
1339
1340 // calculate maximum refreeze amount, maxF [kg]
1341 for(int i=0;i<n;i++)maxF[i] = fmax(0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
1342
1343 // initialize refreeze, runoff, flxDn and dW vectors [kg]
1344 IssmDouble* F = xNewZeroInit<IssmDouble>(n);
1345 IssmDouble* R=xNewZeroInit<IssmDouble>(n);
1346
1347 for(int i=0;i<n;i++)dW[i] = 0;
1348 flxDn=xNewZeroInit<IssmDouble>(n+1); for(int i=0;i<n;i++)flxDn[i+1]=F[i];
1349
1350 // determine the deepest grid cell where melt/pore water is generated
1351 X = 0;
1352 for(int i=n-1;i>=0;i--){
1353 if(M[i]>0 || reCast<int,IssmDouble>(exsW[i])){
1354 X=i;
1355 break;
1356 }
1357 }
1358
1359 //// meltwater percolation
1360 for(int i=0;i<n;i++){
1361 // calculate total melt water entering cell
1362 IssmDouble inM = M[i]+ flxDn[i];
1363
1364 // break loop if there is no meltwater and if depth is > mw_depth
1365 if (inM == 0 && i > X){
1366 break;
1367 }
1368
1369 // if reaches impermeable ice layer all liquid water runs off (R)
1370 else if (d[i] >= dIce){ // dPHC = pore hole close off [kg m-3]
1371 // _printf_("ICE LAYER");
1372 // no water freezes in this cell
1373 // no water percolates to lower cell
1374 // cell ice temperature & density do not change
1375
1376 m[i] = m[i] - M[i]; // mass after melt
1377 Wi = (910-d[i]) * Swi * (m[i]/d[i]); // irreducible water
1378 dW[i] = fmin(inM, Wi - W[i]); // change in pore water
1379 R[i] = fmax(0, inM - dW[i]); // runoff
1380 }
1381 // check if no energy to refreeze meltwater
1382 else if (maxF[i] == 0){
1383 // _printf_("REFREEZE == 0");
1384 // no water freezes in this cell
1385 // cell ice temperature & density do not change
1386
1387 m[i] = m[i] - M[i]; // mass after melt
1388 Wi = (910-d[i]) * Swi * (m[i]/d[i]); // irreducible water
1389 dW[i] = fmin(inM, Wi-W[i]); // change in pore water
1390 flxDn[i+1] = fmax(0, inM-dW[i]); // meltwater out
1391 F[i] = 0; // no freeze
1392 }
1393 // some or all meltwater refreezes
1394 else{
1395 // change in density density and temperature
1396 // _printf_("MELT REFREEZE");
1397 //-----------------------melt water-----------------------------
1398 IssmDouble dz_0 = m[i]/d[i];
1399 IssmDouble dMax = (dIce - d[i])*dz_0; // d max = dIce
1400 IssmDouble F1 = fmin(fmin(inM,dMax),maxF[i]); // maximum refreeze
1401 m[i] = m[i] + F1; // mass after refreeze
1402 d[i] = m[i]/dz_0;
1403
1404 //-----------------------pore water-----------------------------
1405 Wi = (910-d[i])* Swi * dz_0; // irreducible water
1406 dW[i] = fmin(inM - F1, Wi-W[i]); // change in pore water
1407 if (-dW[i]>W[i] ){
1408 dW[i]= W[i];
1409 }
1410 IssmDouble F2 = 0;
1411
1412 if (dW[i] < 0){ // excess pore water
1413 dMax = (dIce - d[i])*dz_0; // maximum refreeze
1414 IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze
1415 F2 = fmin(-dW[i], maxF2); // pore water refreeze
1416 m[i] = m[i] + F2; // mass after refreeze
1417 d[i] = m[i]/dz_0;
1418 }
1419
1420 flxDn[i+1] = inM - F1 - dW[i] - F2; // meltwater out
1421 T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
1422
1423
1424 // check if an ice layer forms
1425 if (d[i] == dIce){
1426 // _printf_("ICE LAYER FORMS");
1427 // excess water runs off
1428 R[i] = flxDn[i+1];
1429 // no water percolates to lower cell
1430 flxDn[i+1] = 0;
1431 }
1432 }
1433 }
1434
1435
1436 //// GRID CELL SPACING AND MODEL DEPTH
1437 for(int i=0;i<n;i++)if (W[i] < 0) _error_("negative pore water generated in melt equations");
1438
1439 // delete all cells with zero mass
1440 // adjust pore water
1441 for(int i=0;i<n;i++)W[i] += dW[i];
1442
1443 // delete all cells with zero mass
1444 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; D=xNew<int>(D_size);
1445 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0){ D[D_size] = i; D_size++;}
1446
1447 celldelete(&m,n,D,D_size);
1448 celldelete(&W,n,D,D_size);
1449 celldelete(&d,n,D,D_size);
1450 celldelete(&T,n,D,D_size);
1451 celldelete(&a,n,D,D_size);
1452 celldelete(&re,n,D,D_size);
1453 celldelete(&gdn,n,D,D_size);
1454 celldelete(&gsp,n,D,D_size);
1455 celldelete(&EI,n,D,D_size);
1456 celldelete(&EW,n,D,D_size);
1457 celldelete(&R,n,D,D_size);
1458 n=D_size;
1459 xDelete<int>(D);
1460
1461 // calculate new grid lengths
1462 for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
1463
1464 //calculate Rsum:
1465 Rsum=cellsum(R,n);
1466
1467 /*Free ressources:*/
1468 xDelete<IssmDouble>(F);
1469 xDelete<IssmDouble>(R);
1470 }
1471
1472 // check if depth is too small
1473 X = 0;
1474 for(int i=n-1;i>=0;i--){
1475 if(dz[i]<dzMin){
1476 X=i;
1477 break;
1478 }
1479 }
1480
1481 for (int i = 0; i<=X;i++){
1482 if (dz [i] < dzMin){ // merge top cells
1483 // _printf_("dz > dzMin * 2')
1484 // adjust variables as a linearly weighted function of mass
1485 IssmDouble m_new = m[i] + m[i+1];
1486 T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
1487 a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
1488 re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
1489 gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
1490 gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
1491
1492 // merge with underlying grid cell and delete old cell
1493 dz [i+1] = dz[i] + dz[i+1]; // combine cell depths
1494 d[i+1] = m_new / dz[i+1]; // combine top densities
1495 W[i+1] = W[i+1] + W[i]; // combine liquid water
1496 m[i+1] = m_new; // combine top masses
1497
1498 // set cell to 99999 for deletion
1499 m[i] = 99999;
1500 }
1501 }
1502
1503 // delete combined cells
1504 D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size);
1505 D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
1506
1507 celldelete(&m,n,D,D_size);
1508 celldelete(&W,n,D,D_size);
1509 celldelete(&dz,n,D,D_size);
1510 celldelete(&d,n,D,D_size);
1511 celldelete(&T,n,D,D_size);
1512 celldelete(&a,n,D,D_size);
1513 celldelete(&re,n,D,D_size);
1514 celldelete(&gdn,n,D,D_size);
1515 celldelete(&gsp,n,D,D_size);
1516 celldelete(&EI,n,D,D_size);
1517 celldelete(&EW,n,D,D_size);
1518 n=D_size;
1519 xDelete<int>(D);
1520
1521 // check if any of the top 10 cell depths are too large
1522 X=0;
1523 for(int i=9;i>=0;i--){
1524 if(dz[i]> 2* dzMin){
1525 X=i;
1526 break;
1527 }
1528 }
1529
1530 i=0;
1531 while(i<=X){
1532 if (dz [i] > dzMin *2){
1533
1534 // _printf_("dz > dzMin * 2");
1535 // split in two
1536 cellsplit(&dz, n, i,.5);
1537 cellsplit(&W, n, i,.5);
1538 cellsplit(&m, n, i,.5);
1539 cellsplit(&T, n, i,1.0);
1540 cellsplit(&d, n, i,1.0);
1541 cellsplit(&a, n, i,1.0);
1542 cellsplit(&EI, n, i,1.0);
1543 cellsplit(&EW, n, i,1.0);
1544 cellsplit(&re, n, i,1.0);
1545 cellsplit(&gdn, n, i,1.0);
1546 cellsplit(&gsp, n, i,1.0);
1547 n++;
1548 X=X+1;
1549 }
1550 else i++;
1551 }
1552
1553 //// CORRECT FOR TOTAL MODEL DEPTH
1554 // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
1555 // INTERPRETATION
1556
1557 // // calculate total model depth
1558 // z = sum(dz);
1559 //
1560 // if (z < zMin){ // check if model is too shallow
1561 // _printf_("z < zMin')
1562 // // mass and energy to be added
1563 // mAdd = m(end) + W(end);
1564 // addE = T(end) * m(end) * CI;
1565 //
1566 // // add a grid cell of the same size and temperature to the bottom
1567 // dz = [dz; dz(end)];
1568 // T = [T; T(end)];
1569 // W = [W; W(end)];
1570 // m = [m; m(end)];
1571 // d = [d; d(end)];
1572 // a = [a; a(end)];
1573 // re = [re; re(end)];
1574 // gdn = [gdn; gdn(end)];
1575 // gsp = [gsp; gsp(end)];
1576 // }
1577 // else (if z > zMax){ // check if model is too deep
1578 // _printf_("z > zMax')
1579 // // mass and energy loss
1580 // mAdd = -(m(end) + W(end));
1581 // addE = -(T(end) * m(end) * CI);
1582 //
1583 // // add a grid cell of the same size and temperature to the bottom
1584 // dz(end) = []; T(end) = []; W(end) = []; m(end) = [];
1585 // d(end) = []; a(end) = []; re(end) = []; gdn(end) = [];
1586 // gsp(end) = [];
1587 // }
1588
1589 //// CHECK FOR MASS AND ENERGY CONSERVATION
1590
1591 // calculate final mass [kg] and energy [J]
1592 sumER = Rsum * (LF + CtoK * CI);
1593 for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
1594 for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
1595
1596 mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
1597 sumE1 = cellsum(EI,n) + cellsum(EW,n);
1598
1599 /*checks: */
1600 for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
1601
1602 /*only in forward mode! avoid round in AD mode as it is not differentiable: */
1603 #ifndef _HAVE_ADOLC_
1604 dm = round(mSum0 - mSum1 + mAdd);
1605 dE = round(sumE0 - sumE1 - sumER + addE);
1606 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
1607 << "dm: " << dm << " dE: " << dE << "\n");
1608 #endif
1609
1610 /*Free ressources:*/
1611 if(m)xDelete<IssmDouble>(m);
1612 if(EI)xDelete<IssmDouble>(EI);
1613 if(EW)xDelete<IssmDouble>(EW);
1614 if(maxF)xDelete<IssmDouble>(maxF);
1615 if(dW)xDelete<IssmDouble>(dW);
1616 if(exsW)xDelete<IssmDouble>(exsW);
1617 if(exsT)xDelete<IssmDouble>(exsT);
1618 if(surpT)xDelete<IssmDouble>(surpT);
1619 if(surpE)xDelete<IssmDouble>(surpE);
1620 if(flxDn)xDelete<IssmDouble>(flxDn);
1621 if(D)xDelete<int>(D);
1622 if(M)xDelete<IssmDouble>(M);
1623
1624 /*Assign output pointers:*/
1625 *pM=sumM;
1626 *pR=Rsum;
1627 *pmAdd=mAdd;
1628
1629 *pT=T;
1630 *pd=d;
1631 *pdz=dz;
1632 *pW=W;
1633 *pa=a;
1634 *pre=re;
1635 *pgdn=gdn;
1636 *pgsp=gsp;
1637 *pn=n;
1638
1639} /*}}}*/
1640void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble dIce, int m, int sid){ /*{{{*/
1641
1642 //// 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???]
1643
1644 //// FUNCTION INFO
1645
1646 // Author: Alex Gardner, University of Alberta
1647 // Date last modified: FEB, 2008
1648
1649 // Description:
1650 // computes the densification of snow/firn using the emperical model of
1651 // Herron and Langway (1980) or the semi-emperical model of Anthern et al.
1652 // (2010)
1653
1654 // Inputs:
1655 // denIdx = densification model to use:
1656 // 1 = emperical model of Herron and Langway (1980)
1657 // 2 = semi-imerical model of Anthern et al. (2010)
1658 // 3 = physical model from Appendix B of Anthern et al. (2010)
1659 // d = initial snow/firn density [kg m-3]
1660 // T = temperature [K]
1661 // dz = grid cell size [m]
1662 // C = average accumulation rate [kg m-2 yr-1]
1663 // dt = time lapsed [s]
1664 // re = effective grain radius [mm];
1665 // Ta = mean annual temperature
1666
1667 // Reference:
1668 // Herron and Langway (1980), Anthern et al. (2010)
1669
1670 //// FOR TESTING
1671 // denIdx = 2;
1672 // d = 800;
1673 // T = 270;
1674 // dz = 0.005;
1675 // C = 200;
1676 // dt = 60*60;
1677 // re = 0.7;
1678 // Tmean = 273.15-18;
1679
1680 //// MAIN FUNCTION
1681 // specify constants
1682 dt = dt / 86400; // convert from [s] to [d]
1683 // R = 8.314 // gas constant [mol-1 K-1]
1684 // Ec = 60 // activation energy for self-diffusion of water
1685 // // molecules through the ice tattice [kJ mol-1]
1686 // Eg = 42.4 // activation energy for grain growth [kJ mol-1]
1687
1688 /*intermediary: */
1689 IssmDouble c0,c1,H;
1690
1691 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" densification module\n");
1692
1693 // initial mass
1694 IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
1695
1696 /*allocations and initialization of overburden pressure and factor H: */
1697 IssmDouble* cumdz = xNew<IssmDouble>(m-1);
1698 cumdz[0]=dz[0];
1699 for(int i=1;i<m-1;i++)cumdz[i]=cumdz[i-1]+dz[i];
1700
1701 IssmDouble* obp = xNew<IssmDouble>(m);
1702 obp[0]=0;
1703 for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
1704
1705 // calculate new snow/firn density for:
1706 // snow with densities <= 550 [kg m-3]
1707 // snow with densities > 550 [kg m-3]
1708
1709
1710 for(int i=0;i<m;i++){
1711 switch (denIdx){
1712 case 1: // Herron and Langway (1980)
1713 c0 = (11 * exp(-10160 / (T[i] * 8.314))) * C/1000;
1714 c1 = (575 * exp(-21400 / (T[i]* 8.314))) * pow(C/1000,.5);
1715 break;
1716 case 2: // Arthern et al. (2010) [semi-emperical]
1717 // common variable
1718 // NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
1719 H = exp((-60000./(T[i] * 8.314)) + (42400./(Tmean * 8.314))) * (C * 9.81);
1720 c0 = 0.07 * H;
1721 c1 = 0.03 * H;
1722 break;
1723
1724 case 3: // Arthern et al. (2010) [physical model eqn. B1]
1725
1726 // common variable
1727 H = exp((-60/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000,2.0);
1728 c0 = 9.2e-9 * H;
1729 c1 = 3.7e-9 * H;
1730 break;
1731
1732 case 4: // Li and Zwally (2004)
1733 c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(273.15 - T[i],-2.061);
1734 c1 = c0;
1735 break;
1736
1737 case 5: // Helsen et al. (2008)
1738 // common variable
1739 c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(273.15 - T[i],-2.061);
1740 c1 = c0;
1741 break;
1742 }
1743
1744 // new snow density
1745 if(d[i] <= 550) d[i] = d[i] + (c0 * (dIce - d[i]) / 365 * dt);
1746 else d[i] = d[i] + (c1 * (dIce - d[i]) / 365 * dt);
1747
1748 //disp((num2str(nanmean(c0 .* (dIce - d(idx)) / 365 * dt))))
1749
1750 // do not allow densities to exceed the density of ice
1751 if(d[i]>dIce)d[i]=dIce;
1752
1753 // calculate new grid cell length
1754 dz[i] = mass_init[i] / d[i];
1755 }
1756
1757 /*Free ressources:*/
1758 xDelete<IssmDouble>(mass_init);
1759 xDelete<IssmDouble>(cumdz);
1760 xDelete<IssmDouble>(obp);
1761
1762} /*}}}*/
1763void 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, int sid){ /*{{{*/
1764
1765 //// TURBULENT HEAT FLUX
1766
1767 // Description:
1768 // computed the surface sensible and latent heat fluxes [W m-2], this
1769 // function also calculated the mass loss/acreation due to
1770 // condensation/evaporation [kg]
1771
1772 // Reference:
1773 // Dingman, 2002.
1774
1775 //// INPUTS:
1776 // Ta: 2m air temperature [K]
1777 // Ts: snow/firn/ice surface temperature [K]
1778 // V: wind speed [m s^-^1]
1779 // eAir: screen level vapor pressure [Pa]
1780 // pAir: surface pressure [Pa]
1781 // ds: surface density [kg/m^3]
1782 // Ws: surface liquid water content [kg/m^2]
1783 // Vz: height above ground at which wind (V) eas sampled [m]
1784 // Tz: height above ground at which temperature (T) was sampled [m]
1785
1786 //// FUNCTION INITILIZATION
1787
1788 // CA = 1005; // heat capacity of air (J kg-1 k-1)
1789 // LF = 0.3345E6; // latent heat of fusion(J kg-1)
1790 // LV = 2.495E6; // latent heat of vaporization(J kg-1)
1791 // dIce = 910; // density of ice [kg m-3]
1792
1793 /*intermediary:*/
1794 IssmDouble d_air;
1795 IssmDouble Ri;
1796 IssmDouble z0;
1797 IssmDouble coef_M,coef_H;
1798 IssmDouble An, C;
1799 IssmDouble L, eS;
1800
1801 /*output: */
1802 IssmDouble shf, lhf, EC;
1803
1804 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" turbulentFlux module\n");
1805
1806 // calculated air density [kg/m3]
1807 d_air = 0.029 * pAir /(8.314 * Ta);
1808
1809 //// Determine Surface Roughness
1810 // Bougamont, 2006
1811 // wind/temperature surface roughness height [m]
1812 if (ds < 910 && Ws == 0) z0 = 0.00012; // 0.12 mm for dry snow
1813 else if (ds >= 910) z0 = 0.0032; // 3.2 mm for ice
1814 else z0 = 0.0013; // 1.3 mm for wet snow
1815
1816 //// Monin–Obukhov Stability Correction
1817 // Reference:
1818 // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
1819 // Journal of Climatology, 2, 65-84.
1820
1821 // if V = 0 goes to infinity therfore if V = 0 change
1822 if(V< .01) V=.01;
1823
1824 // calculate the Bulk Richardson Number (Ri)
1825 Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
1826
1827 // calculate Monin–Obukhov stability factors 'coef_M' and 'coef_H'
1828
1829 // do not allow Ri to exceed 0.19
1830 if(Ri>.19)Ri= 0.19;
1831
1832 // calculate momentum 'coef_M' stability factor
1833 if (Ri > 0) coef_M = pow(1-5.2*Ri,-1); // if stable
1834 else coef_M = pow(1-18*Ri,-0.25);
1835
1836 // calculate heat/wind 'coef_H' stability factor
1837 if (Ri < -0.03) coef_H = 1.3 * coef_M;
1838 else coef_H = coef_M;
1839
1840 //// Bulk-transfer coefficient
1841 An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient
1842 C = An * d_air * V; // shf & lhf common coefficient
1843
1844 //// Sensible Heat
1845 // calculate the sensible heat flux [W m-2](Patterson, 1998)
1846 shf = C * 1005 * (Ta - Ts);
1847
1848 // adjust using Monin–Obukhov stability theory
1849 shf = shf / (coef_M * coef_H);
1850
1851 //// Latent Heat
1852 // determine if snow pack is melting & calcualte surface vapour pressure
1853 // over ice or liquid water
1854 if (Ts >= 273.15){
1855 L = 2.495E6;
1856
1857 // for an ice surface Murphy and Koop, 2005 [Equation 7]
1858 eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts)- 0.00728332 * Ts);
1859 }
1860 else{
1861 L = 2.8295E6; // latent heat of sublimation
1862 // for liquid surface (assume liquid on surface when Ts == 0 deg C)
1863 // Wright (1997), US Meteorological Handbook from Murphy and Koop,
1864 // 2005 Apendix A
1865 eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15));
1866 }
1867
1868 // Latent heat flux [W m-2]
1869 lhf = C * L * (eAir - eS) * 0.622 / pAir;
1870
1871 // adjust using Monin–Obukhov stability theory (if lhf '+' then there is
1872 // energy and mass gained at the surface, if '-' then there is mass and
1873 // energy loss at the surface.
1874 lhf = lhf / (coef_M * coef_H);
1875
1876 // mass loss (-)/acreation(+) due to evaporation/condensation [kg]
1877 EC = lhf * 86400 / L;
1878
1879 /*assign output poitners: */
1880 *pshf=shf;
1881 *plhf=lhf;
1882 *pEC=EC;
1883
1884} /*}}}*/
Note: See TracBrowser for help on using the repository browser.