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

Last change on this file since 24686 was 24686, checked in by Mathieu Morlighem, 5 years ago

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