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

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago

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