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

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

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