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

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

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