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

Last change on this file since 23394 was 23394, checked in by Mathieu Morlighem, 6 years ago

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