Ice Sheet System Model  4.18
Code documentation
Functions | Variables
Gembx.cpp File Reference
#include "./SurfaceMassBalancex.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../modules.h"
#include "../../classes/Inputs2/TransientInput2.h"

Go to the source code of this file.

Functions

void Gembx (FemModel *femmodel)
 
void GembgridInitialize (IssmDouble **pdz, int *psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY)
 
IssmDouble Marbouty (IssmDouble T, IssmDouble d, IssmDouble dT)
 
void grainGrowth (IssmDouble **pre, IssmDouble **pgdn, IssmDouble **pgsp, IssmDouble *T, IssmDouble *dz, IssmDouble *d, IssmDouble *W, IssmDouble smb_dt, int m, int aIdx, int sid)
 
void 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)
 
void 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)
 
void shortwave (IssmDouble **pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble *d, IssmDouble *dz, IssmDouble *re, IssmDouble dIce, int m, int sid)
 
void 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)
 
void 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)
 
void densification (IssmDouble **pd, IssmDouble **pdz, IssmDouble *T, IssmDouble *re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid)
 
void 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)
 

Variables

const double Pi = 3.141592653589793
 
const double CtoK = 273.15
 
const double dts = 86400.0
 
const double Ttol = 1e-10
 
const double Dtol = 1e-11
 
const double Gdntol = 1e-10
 
const double Wtol = 1e-13
 
const double Ptol = 1e-6
 
const double CI = 2102.0
 
const double LF = 0.3345E6
 
const double LV = 2.495E6
 
const double LS = 2.8295E6
 
const double SB = 5.67E-8
 
const double CA = 1005.0
 
const double R = 8.314
 
const double Delflag =-99999
 

Function Documentation

◆ Gembx()

void Gembx ( FemModel *  femmodel)

Definition at line 34 of file Gembx.cpp.

34  { /*{{{*/
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(int i=0;i<femmodel->elements->Size();i++){
62  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
63 
64  timeclim=time;
65  if (isclimatology){
66  //If this is a climatology, we need to repeat the forcing after the final time
67  TransientInput2* Ta_input_tr = element->inputs2->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 } /*}}}*/

◆ GembgridInitialize()

void GembgridInitialize ( IssmDouble **  pdz,
int *  psize,
IssmDouble  zTop,
IssmDouble  dzTop,
IssmDouble  zMax,
IssmDouble  zY 
)

Definition at line 87 of file Gembx.cpp.

87  { /*{{{*/
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 } /*}}}*/

◆ Marbouty()

IssmDouble Marbouty ( IssmDouble  T,
IssmDouble  d,
IssmDouble  dT 
)

Definition at line 166 of file Gembx.cpp.

166  { /*{{{*/
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 } /*}}}*/

◆ grainGrowth()

void grainGrowth ( IssmDouble **  pre,
IssmDouble **  pgdn,
IssmDouble **  pgsp,
IssmDouble *  T,
IssmDouble *  dz,
IssmDouble *  d,
IssmDouble *  W,
IssmDouble  smb_dt,
int  m,
int  aIdx,
int  sid 
)

Definition at line 199 of file Gembx.cpp.

199  { /*{{{*/
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 } /*}}}*/

◆ albedo()

void 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 
)

Definition at line 404 of file Gembx.cpp.

404  { /*{{{*/
405 
407  // 0 : direct input from aValue parameter
408  // 1 : effective grain radius (Gardner & Sharp, 2009)
409  // 2 : effective grain radius (Brun et al., 2009)
410  // 3 : density and cloud amount (Greuell & Konzelmann, 1994)
411  // 4 : exponential time decay & wetness (Bougamont & Bamber, 2005)
412 
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 
447  // a = grid cell albedo
448 
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  if (d[0]<dPHC-Dtol){
490 
491  IssmDouble sF[3] = {0.606, 0.301, 0.093};
492 
493  // convert effective radius to grain size in meters
494  IssmDouble gsz = (re[0] * 2.0) / 1000.0;
495 
496  // spectral range:
497  // 0.3 - 0.8um
498  IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5));
499  // 0.8 - 1.5um
500  IssmDouble a1 = max(0., 0.95 - 15.4 *pow(gsz,0.5));
501  // 1.5 - 2.8um
502  IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));
503 
504  // broadband surface albedo
505  a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;
506 
507  // In a layer < 10cm, account for mix of ice and snow,
508  // after P. Alexander et al., 2014
509  IssmDouble depthsnow=0.0;
510  IssmDouble aice=0.0;
511  int lice=0;
512  for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
513  depthsnow=depthsnow+dz[l];
514  lice=l+1;
515  }
516  if (depthsnow<=0.1+Dtol & lice<m & d[lice]>=dPHC-Dtol){
517  aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce);
518  a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1);
519  }
520  }
521  else if (d[0]<dIce-Dtol){ //For continuity of albedo in firn i.e. P. Alexander et al., 2014
522 
523  //ai=ai_max + (as_min - ai_max)*(dI-dIce)/(dC-dIce)
524  //dC is pore close off (830 kg m^-3)
525  //dI is density of the upper firn layer
526 
527  a[0] = ai_max + (as_min - ai_max)*(d[0]-dIce)/(dPHC-dIce);
528 
529  }
530  else{ //surface layer is density of ice
531 
532  //When density is > dIce (typically 910 kg m^-3, 920 is used by Alexander in MAR),
533  //ai=ai_min + (ai_max - ai_min)*e^(-1*(Msw(t)/K))
534  //K is a scale factor (set to 200 kg m^-2)
535  //Msw(t) is the time-dependent accumulated amount of excessive surface meltwater
536  // before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly)
537  IssmDouble M = Msurf+W[0];
538  a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min);
539 
540  }
541  }
542  else if(aIdx==3){
543 
544  // a as a function of density
545 
546  // calculate albedo
547  a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
548  }
549  else if(aIdx==4){
550 
551  // exponential time decay & wetness
552 
553  // change in albedo with time:
554  // (d_a) = (a - a_old)/(t0)
555  // where: t0 = timescale for albedo decay
556 
557  dt = dt / dts; // convert from [s] to [d]
558 
559  // initialize variables
560  IssmDouble* t0=xNew<IssmDouble>(m);
561  IssmDouble* T=xNew<IssmDouble>(m);
562  IssmDouble* t0warm=xNew<IssmDouble>(m);
563  IssmDouble* d_a=xNew<IssmDouble>(m);
564 
565  // specify constants
566  // a_wet = 0.15; // water albedo (0.15)
567  // a_new = aSnow // new snow albedo (0.64 - 0.89)
568  // a_old = aIce; // old snow/ice albedo (0.27-0.53)
569  // t0_wet = t0wet; // time scale for wet snow (15-21.9) [d]
570  // t0_dry = t0dry; // warm snow timescale [15] [d]
571  // K = 7 // time scale temperature coef. (7) [d]
572  // W0 = 300; // 200 - 600 [mm]
573  const IssmDouble z_snow = 15.0; // 16 - 32 [mm]
574 
575  // determine timescale for albedo decay
576  for(int i=0;i<m;i++)if(W[i]>0.0+Wtol)t0[i]=t0wet; // wet snow timescale
577  for(int i=0;i<m;i++)T[i]=TK[i] - CtoK; // change T from K to degC
578  for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry;
579  for(int i=0;i<m;i++)if(fabs(W[i])<Wtol && T[i]>=-10.0-Ttol)t0[i]= t0warm[i];
580  for(int i=0;i<m;i++)if(T[i]<-10.0-Ttol) t0[i] = 10.0 * K + t0dry; // 'cold' snow timescale
581 
582  // calculate new albedo
583  for(int i=0;i<m;i++)d_a[i] = (a[i] - aIce) / t0[i] * dt; // change in albedo
584  for(int i=0;i<m;i++)a[i] -= d_a[i]; // new albedo
585 
586  // modification of albedo due to thin layer of snow or solid
587  // condensation (deposition) at the surface
588 
589  // check if condensation occurs & if it is deposited in solid phase
590  if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0; // add cond to precip [mm]
591 
592  a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);
593 
594  //----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------
595  // modification of albedo due to thin layer of water on the surface
596  // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);
597 
598  /*Free resources:*/
599  xDelete<IssmDouble>(t0);
600  xDelete<IssmDouble>(T);
601  xDelete<IssmDouble>(t0warm);
602  xDelete<IssmDouble>(d_a);
603 
604  }
605  else _error_("albedo method switch should range from 0 to 4!");
606  }
607 
608  // Check for erroneous values
609  if (a[0] > 1) _printf_("albedo > 1.0\n");
610  else if (a[0] < 0) _printf_("albedo is negative\n");
611  else if (xIsNan(a[0])) _error_("albedo == NAN\n");
612 
613  /*Assign output pointers:*/
614  *pa=a;
615 
616 } /*}}}*/

◆ thermo()

void 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 
)

Definition at line 617 of file Gembx.cpp.

617  { /*{{{*/
618 
619  /* ENGLACIAL THERMODYNAMICS*/
620 
621  /* Description:
622  computes new temperature profile accounting for energy absorption and
623  thermal diffusion.*/
624 
625  // INPUTS
626  // T: grid cell temperature [k]
627  // dz: grid cell depth [m]
628  // d: grid cell density [kg m-3]
629  // swf: shortwave radiation fluxes [W m-2]
630  // dlwrf: downward longwave radiation fluxes [W m-2]
631  // Ta: 2 m air temperature
632  // V: wind velocity [m s-1]
633  // eAir: screen level vapor pressure [Pa]
634  // Ws: surface water content [kg]
635  // dt0: time step of input data [s]
636  // elev: surface elevation [m a.s.l.]
637  // Vz: air temperature height above surface [m]
638  // Tz: wind height above surface [m]
639  // thermo_scaling: scaling factor to multiply the thermal diffusion timestep (delta t)
640 
641  // OUTPUTS
642  // T: grid cell temperature [k]
643  // EC: evaporation/condensation [kg]
644  // ulwrf: upward longwave radiation flux [W m-2]
645 
646  /*intermediary: */
647  IssmDouble* K = NULL;
648  IssmDouble* KU = NULL;
649  IssmDouble* KD = NULL;
650  IssmDouble* KP = NULL;
651  IssmDouble* Au = NULL;
652  IssmDouble* Ad = NULL;
653  IssmDouble* Ap = NULL;
654  IssmDouble* Nu = NULL;
655  IssmDouble* Nd = NULL;
656  IssmDouble* Np = NULL;
657  IssmDouble* dzU = NULL;
658  IssmDouble* dzD = NULL;
659  IssmDouble* sw = NULL;
660  IssmDouble* dT_sw = NULL;
661  IssmDouble* T0 = NULL;
662  IssmDouble* Tu = NULL;
663  IssmDouble* Td = NULL;
664 
665  IssmDouble z0=0.0;
666  IssmDouble dt=0.0;
667  IssmDouble max_fdt=0.0;
668  IssmDouble Ts=0.0;
669  IssmDouble L=0.0;
670  IssmDouble eS=0.0;
671  IssmDouble Ri=0.0;
672  IssmDouble coefM=0.0;
673  IssmDouble coefH=0.0;
674  IssmDouble An=0.0;
675  IssmDouble C=0.0;
676  IssmDouble shf=0.0;
677  IssmDouble ds=0.0;
678  IssmDouble dAir=0.0;
679  IssmDouble TCs=0.0;
680  IssmDouble lhf=0.0;
681  IssmDouble EC_day=0.0;
682  IssmDouble dT_turb=0.0;
683  IssmDouble turb=0.0;
684  IssmDouble ulw=0.0;
685  IssmDouble dT_ulw=0.0;
686  IssmDouble dlw=0.0;
687  IssmDouble dT_dlw=0.0;
688 
689  /*outputs:*/
690  IssmDouble EC=0.0;
691  IssmDouble* T=*pT;
692  IssmDouble ulwrf=0.0;
693 
694  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" thermal module\n");
695 
696  ds = d[0]; // density of top grid cell
697 
698  // calculated air density [kg/m3]
699  dAir = 0.029 * pAir /(R * Ta);
700 
701  // thermal capacity of top grid cell [J/k]
702  TCs = d[0]*dz[0]*CI;
703 
704  //initialize Evaporation - Condenstation
705  EC = 0.0;
706 
707  // check if all SW applied to surface or distributed throught subsurface
708  // swIdx = length(swf) > 1
709 
710  // SURFACE ROUGHNESS (Bougamont, 2006)
711  // wind/temperature surface roughness height [m]
712  if (ds < dIce-Dtol && fabs(Ws) < Wtol) z0 = 0.00012; // 0.12 mm for dry snow
713  else if (ds >= dIce-Dtol) z0 = 0.0032; // 3.2 mm for ice
714  else z0 = 0.0013; // 1.3 mm for wet snow
715 
716  // if V = 0 goes to infinity therfore if V = 0 change
717  if(V<0.01-Dtol)V=0.01;
718 
719  // Bulk-transfer coefficient for turbulent fluxes
720  //An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient
721  An = pow(0.4,2) / (log(Tz/z0)*log(Vz/z0)); // Bulk-transfer coefficient
722  C = An * dAir * V; // shf & lhf common coefficient
723 
724  // THERMAL CONDUCTIVITY (Sturm, 1997: J. Glaciology)
725  // calculate new K profile [W m-1 K-1]
726 
727  // initialize conductivity
728  K= xNewZeroInit<IssmDouble>(m);
729 
730  // for snow and firn (density < 910 kg m-3) (Sturn et al, 1997)
731  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));
732 
733  // for ice (density >= 910 kg m-3)
734  for(int i=0;i<m;i++) if(d[i]>=dIce-Dtol) K[i] = 9.828 * exp(-5.7E-3*T[i]);
735 
736  // THERMAL DIFFUSION COEFFICIENTS
737 
738  // A discretization scheme which truncates the Taylor-Series expansion
739  // after the 3rd term is used. See Patankar 1980, Ch. 3&4
740 
741  // discretized heat equation:
742 
743  // Tp = (Au*Tuo+ Ad*Tdo+ (Ap-Au-Ad)Tpo+ S) / Ap
744 
745  // where neighbor coefficients Au, Ap, & Ad are
746 
747  // Au = [dz_u/2KP + dz_p/2KE]^-1
748  // Ad = [dz_d/2KP + dz_d/2KD]^-1
749  // Ap = d*CI*dz/Dt
750 
751  // and u & d represent grid points up and down from the center grid point
752  // point p and o identifies previous time step values. S is a source term.
753 
754  // u, d, and p conductivities
755  KU = xNew<IssmDouble>(m);
756  KD = xNew<IssmDouble>(m);
757  KP = xNew<IssmDouble>(m);
758 
759  KU[0] = Delflag;
760  KD[m-1] = Delflag;
761  for(int i=1;i<m;i++) KU[i]= K[i-1];
762  for(int i=0;i<m-1;i++) KD[i] = K[i+1];
763  for(int i=0;i<m;i++) KP[i] = K[i];
764 
765  // determine u, d & p cell widths
766  dzU = xNew<IssmDouble>(m);
767  dzD = xNew<IssmDouble>(m);
768  dzU[0]=Delflag;
769  dzD[m-1]=Delflag;
770 
771  for(int i=1;i<m;i++) dzU[i]= dz[i-1];
772  for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
773 
774  // determine minimum acceptable delta t (diffusion number > 1/2) [s]
775  // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability
776  dt=1e12;
777  for(int i=0;i<m;i++)dt = min(dt,CI * pow(dz[i],2) * d[i] / (3 * K[i]) * thermo_scaling);
778 
779  // smallest possible even integer of 60 min where diffusion number > 1/2
780  // must go evenly into one hour or the data frequency if it is smaller
781 
782  // all integer factors of the number of second in a day (86400 [s])
783  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,
784  72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600};
785 
786  // return the min integer factor that is < dt
787  max_fdt=f[0];
788  for(int i=0;i<45;i++){
789  if (f[i]<dt-Dtol)if(f[i]>=max_fdt-Dtol)max_fdt=f[i];
790  }
791  dt=max_fdt;
792 
793  // determine mean (harmonic mean) of K/dz for u, d, & p
794  Au = xNew<IssmDouble>(m);
795  Ad = xNew<IssmDouble>(m);
796  Ap = xNew<IssmDouble>(m);
797  for(int i=0;i<m;i++){
798  Au[i] = pow((dzU[i]/2.0/KP[i] + dz[i]/2.0/KU[i]),-1.0);
799  Ad[i] = pow((dzD[i]/2.0/KP[i] + dz[i]/2.0/KD[i]),-1.0);
800  Ap[i] = (d[i]*dz[i]*CI)/dt;
801  }
802 
803  // create "neighbor" coefficient matrix
804  Nu = xNew<IssmDouble>(m);
805  Nd = xNew<IssmDouble>(m);
806  Np = xNew<IssmDouble>(m);
807  for(int i=0;i<m;i++){
808  Nu[i] = Au[i] / Ap[i];
809  Nd[i] = Ad[i] / Ap[i];
810  Np[i]= 1.0 - Nu[i] - Nd[i];
811  }
812 
813  // specify boundary conditions: constant flux at bottom
814  Nu[m-1] = 0.0;
815  Np[m-1] = 1.0;
816 
817  // zero flux at surface
818  Np[0] = 1.0 - Nd[0];
819 
820  // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
821  Nu[0] = 0.0;
822  Nd[m-1] = 0.0;
823 
824  /* RADIATIVE FLUXES*/
825 
826  // energy supplied by shortwave radiation [J]
827  sw = xNew<IssmDouble>(m);
828  for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
829 
830  // temperature change due to SW
831  dT_sw = xNew<IssmDouble>(m);
832  for(int i=0;i<m;i++) dT_sw[i]= sw[i] / (CI * d[i] * dz[i]);
833 
834  // Upward longwave radiation flux is calculated from the snow surface
835  // temperature which is set equal to the average temperature of the
836  // top grid cells.
837 
838  // energy supplied by downward longwave radiation to the top grid cell [J]
839  dlw = dlwrf * dt;
840 
841  // temperature change due to dlw_surf
842  dT_dlw = dlw / TCs;
843 
844  // PREALLOCATE ARRAYS BEFORE LOOP FOR IMPROVED PERFORMANCE
845  T0 = xNewZeroInit<IssmDouble>(m+2);
846  Tu=xNew<IssmDouble>(m);
847  Td=xNew<IssmDouble>(m);
848 
849  /* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
850  for (IssmDouble i=1;i<=dt0;i+=dt){
851 
852  // PART OF ENERGY CONSERVATION CHECK
853  // store initial temperature
854  //T_init = T;
855 
856  // calculate temperature of snow surface (Ts)
857  // when incoming SW radition is allowed to penetrate the surface,
858  // the modeled energy balance becomes very sensitive to how Ts is
859  // calculated. The estimated enegy balance & melt are significanly
860  // less when Ts is taken as the mean of the x top grid cells.
861  if(m>1) Ts = (T[0] + T[1])/2.0;
862  else Ts = T[0];
863  Ts = min(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC)
864 
865  //TURBULENT HEAT FLUX
866 
867  // Monin-Obukhov Stability Correction
868  // Reference:
869  // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
870  // Journal of Climatology, 2, 65-84.
871 
872  // calculate the Bulk Richardson Number (Ri)
873  Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
874 
875  // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
876 
877  // do not allow Ri to exceed 0.19
878  Ri = min(Ri, 0.19);
879 
880  // calculate momentum 'coefM' stability factor
881  if (Ri > 0.0+Ttol){
882  // if stable
883  coefM = 1.0/(1.0-5.2*Ri);
884  }
885  else {
886  coefM =pow (1.0-18*Ri,-0.25);
887  }
888 
889  // calculate heat/wind 'coef_H' stability factor
890  if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
891  else coefH = coefM;
892 
894  // calculate the sensible heat flux [W m-2](Patterson, 1998)
895  shf = C * CA * (Ta - Ts);
896 
897  // adjust using Monin-Obukhov stability theory
898  shf = shf / (coefM * coefH);
899 
901  // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
902  if (Ts >= CtoK-Ttol){
903  L = LV; //for liquid water at 273.15 k to vapor
904 
905  //for liquid surface (assume liquid on surface when Ts == 0 deg C)
906  // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
907  eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
908  }
909  else{
910  L = LS; // latent heat of sublimation
911 
912  // for an ice surface Murphy and Koop, 2005 [Equation 7]
913  eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
914  }
915 
916  // Latent heat flux [W m-2]
917  lhf = C * L * (eAir - eS) * 0.622 / pAir;
918 
919  // adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface,
920  // if '-' then there is mass and energy loss at the surface.
921  lhf = lhf / (coefM * coefH);
922 
923  //mass loss (-)/acreation(+) due to evaporation/condensation [kg]
924  EC_day = lhf * dts / L;
925 
926  // temperature change due turbulent fluxes
927  turb = (shf + lhf)* dt;
928  dT_turb = turb / TCs;
929 
930  // upward longwave contribution
931  ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here
932  ulwrf = ulwrf - ulw/dt0;
933 
934  dT_ulw = ulw / TCs;
935 
936  // new grid point temperature
937 
938  //SW penetrates surface
939  if(!isconstrainsurfaceT) for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
940  if(!isconstrainsurfaceT) T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
941 
942  // temperature diffusion
943  for(int j=0;j<m;j++) T0[1+j]=T[j];
944  for(int j=0;j<m;j++) Tu[j] = T0[j];
945  for(int j=0;j<m;j++) Td[j] = T0[2+j];
946  for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
947 
948  // calculate cumulative evaporation (+)/condensation(-)
949  if(!isconstrainsurfaceT) EC = EC + (EC_day/dts)*dt;
950 
951  /* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
952  //energy flux across lower boundary (energy supplied by underling ice)
953  base_flux = Ad(-1)*(T_init()-T_init(-1)) * dt;
954 
955  E_used = sum((T - T_init) * (d*dz*CI));
956  E_sup = ((sum(swf) * dt) + dlw + ulw + turb + base_flux);
957 
958  E_diff = E_used - E_sup;
959 
960  if abs(E_diff) > 1E-6 || isnan(E_diff)
961  disp(T(1))
962  _error_("energy not conserved in thermodynamics equations");
963  */
964  }
965 
966  /*Free resources:*/
967  xDelete<IssmDouble>(K);
968  xDelete<IssmDouble>(KU);
969  xDelete<IssmDouble>(KD);
970  xDelete<IssmDouble>(KP);
971  xDelete<IssmDouble>(Au);
972  xDelete<IssmDouble>(Ad);
973  xDelete<IssmDouble>(Ap);
974  xDelete<IssmDouble>(Nu);
975  xDelete<IssmDouble>(Nd);
976  xDelete<IssmDouble>(Np);
977  xDelete<IssmDouble>(dzU);
978  xDelete<IssmDouble>(dzD);
979  xDelete<IssmDouble>(sw);
980  xDelete<IssmDouble>(dT_sw);
981  xDelete<IssmDouble>(T0);
982  xDelete<IssmDouble>(Tu);
983  xDelete<IssmDouble>(Td);
984 
985  /*Assign output pointers:*/
986  *pEC=EC;
987  *pT=T;
988  *pulwrf=ulwrf;
989 
990 } /*}}}*/

◆ shortwave()

void shortwave ( IssmDouble **  pswf,
int  swIdx,
int  aIdx,
IssmDouble  dsw,
IssmDouble  as,
IssmDouble *  d,
IssmDouble *  dz,
IssmDouble *  re,
IssmDouble  dIce,
int  m,
int  sid 
)

Definition at line 991 of file Gembx.cpp.

991  { /*{{{*/
992 
993  // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
994 
995  // swIdx = 0 : all absorbed SW energy is assigned to the top grid cell
996 
997  // swIdx = 1 : absorbed SW is distributed with depth as a function of:
998  // 1 : snow density (taken from Bassford, 2004)
999  // 2 : grain size in 3 spectral bands (Brun et al., 1992)
1000 
1001  // Inputs
1002  // swIdx = shortwave allowed to penetrate surface (0 = No, 1 = Yes)
1003  // aIdx = method for calculating albedo (1-4)
1004  // dsw = downward shortwave radiative flux [w m-2]
1005  // as = surface albedo
1006  // d = grid cell density [kg m-3]
1007  // dz = grid cell depth [m]
1008  // re = grid cell effective grain radius [mm]
1009 
1010  // Outputs
1011  // swf = absorbed shortwave radiation [W m-2]
1012  //
1013 
1014  /*outputs: */
1015  IssmDouble* swf=NULL;
1016 
1017  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" shortwave module\n");
1018 
1019  /*Initialize and allocate: */
1020  swf=xNewZeroInit<IssmDouble>(m);
1021 
1022  // SHORTWAVE FUNCTION
1023  if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
1024 
1025  // calculate surface shortwave radiation fluxes [W m-2]
1026  swf[0] = (1.0 - as) * dsw;
1027  }
1028  else{ // sw radation is absorbed at depth within the glacier
1029 
1030  if (aIdx == 2){ // function of effective radius (3 spectral bands)
1031 
1032  IssmDouble * gsz=NULL;
1033  IssmDouble * B1_cum=NULL;
1034  IssmDouble * B2_cum=NULL;
1035  IssmDouble* h =NULL;
1036  IssmDouble* B1 =NULL;
1037  IssmDouble* B2 =NULL;
1038  IssmDouble* exp1 = NULL;
1039  IssmDouble* exp2 = NULL;
1040  IssmDouble* Qs1 = NULL;
1041  IssmDouble* Qs2 = NULL;
1042 
1043  // convert effective radius [mm] to grain size [m]
1044  gsz=xNew<IssmDouble>(m);
1045  for(int i=0;i<m;i++) gsz[i]= (re[i] * 2.0) / 1000.0;
1046 
1047  // Spectral fractions [0.3-0.8um 0.8-1.5um 1.5-2.8um]
1048  // (Lefebre et al., 2003)
1049  IssmDouble sF[3] = {0.606, 0.301, 0.093};
1050 
1051  // initialize variables
1052  B1_cum=xNew<IssmDouble>(m+1);
1053  B2_cum=xNew<IssmDouble>(m+1);
1054  for(int i=0;i<m+1;i++){
1055  B1_cum[i]=1.0;
1056  B2_cum[i]=1.0;
1057  }
1058 
1059  // spectral albedos:
1060  // 0.3 - 0.8um
1061  IssmDouble a0 = min(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
1062  // 0.8 - 1.5um
1063  IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
1064  // 1.5 - 2.8um
1065  IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
1066 
1067  // separate net shortwave radiative flux into spectral ranges
1068  IssmDouble swfS[3];
1069  swfS[0] = (sF[0] * dsw) * (1.0 - a0);
1070  swfS[1] = (sF[1] * dsw) * (1.0 - a1);
1071  swfS[2] = (sF[2] * dsw) * (1.0 - a2);
1072 
1073  // absorption coeficient for spectral range:
1074  h =xNew<IssmDouble>(m);
1075  B1 =xNew<IssmDouble>(m);
1076  B2 =xNew<IssmDouble>(m);
1077  for(int i=0;i<m;i++) h[i]= d[i] /(pow(gsz[i],0.5));
1078  for(int i=0;i<m;i++) B1[i] = 0.0192 * h[i]; // 0.3 - 0.8um
1079  for(int i=0;i<m;i++) B2[i]= 0.1098 * h[i]; // 0.8 - 1.5um
1080  // B3 = +inf // 1.5 - 2.8um
1081 
1082  // cumulative extinction factors
1083  exp1 = xNew<IssmDouble>(m);
1084  exp2 = xNew<IssmDouble>(m);
1085  for(int i=0;i<m;i++) exp1[i]=exp(-B1[i]*dz[i]);
1086  for(int i=0;i<m;i++) exp2[i]=exp(-B2[i]*dz[i]);
1087 
1088  for(int i=0;i<m;i++){
1089  IssmDouble cum1=exp1[0];
1090  IssmDouble cum2=exp2[0];
1091  for(int j=1;j<=i;j++){
1092  cum1 = cum1*exp1[j];
1093  cum2 = cum2*exp2[j];
1094  }
1095  B1_cum[i+1]=cum1;
1096  B2_cum[i+1]=cum2;
1097  }
1098 
1099  // flux across grid cell boundaries
1100  Qs1 = xNew<IssmDouble>(m+1);
1101  Qs2 = xNew<IssmDouble>(m+1);
1102  for(int i=0;i<m+1;i++){
1103  Qs1[i] = swfS[0] * B1_cum[i];
1104  Qs2[i] = swfS[1] * B2_cum[i];
1105  }
1106 
1107  // net energy flux to each grid cell
1108  for(int i=0;i<m;i++) swf[i]= (Qs1[i]-Qs1[i+1]) + (Qs2[i]-Qs2[i+1]);
1109 
1110  // add flux absorbed at surface
1111  swf[0] = swf[0]+ swfS[2];
1112 
1113  /*Free resources: */
1114  xDelete<IssmDouble>(gsz);
1115  xDelete<IssmDouble>(B1_cum);
1116  xDelete<IssmDouble>(B2_cum);
1117  xDelete<IssmDouble>(h);
1118  xDelete<IssmDouble>(B1);
1119  xDelete<IssmDouble>(B2);
1120  xDelete<IssmDouble>(exp1);
1121  xDelete<IssmDouble>(exp2);
1122  xDelete<IssmDouble>(Qs1);
1123  xDelete<IssmDouble>(Qs2);
1124 
1125  }
1126  else{ //function of grid cell density
1127 
1128  /*intermediary: */
1129  IssmDouble* B_cum = NULL;
1130  IssmDouble* exp_B = NULL;
1131  IssmDouble* Qs = NULL;
1132  IssmDouble* B = NULL;
1133 
1134  // fraction of sw radiation absorbed in top grid cell (wavelength > 0.8um)
1135  IssmDouble SWs = 0.36;
1136 
1137  // SWs and SWss coefficients need to be better constranted. Greuell
1138  // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
1139  // the % of SW radiation with wavelengths > and < 800 nm
1140  // respectively. This, however, may not account for the fact that
1141  // the albedo of wavelengths > 800 nm has a much lower albedo.
1142 
1143  // calculate surface shortwave radiation fluxes [W m-2]
1144  IssmDouble swf_s = SWs * (1.0 - as) * dsw;
1145 
1146  // calculate surface shortwave radiation fluxes [W m-2]
1147  IssmDouble swf_ss = (1.0-SWs) * (1.0 - as) * dsw;
1148 
1149  // SW allowed to penetrate into snowpack
1150  IssmDouble Bs = 10.0; // snow SW extinction coefficient [m-1] (Bassford,2006)
1151  IssmDouble Bi = 1.3; // ice SW extinction coefficient [m-1] (Bassford,2006)
1152 
1153  // calculate extinction coefficient B [m-1] vector
1154  B=xNew<IssmDouble>(m);
1155  for(int i=0;i<m;i++) B[i] = Bs + (300.0 - d[i]) * ((Bs - Bi)/(dIce - 300.0));
1156 
1157  // cumulative extinction factor
1158  B_cum = xNew<IssmDouble>(m+1);
1159  exp_B = xNew<IssmDouble>(m);
1160  for(int i=0;i<m;i++)exp_B[i]=exp(-B[i]*dz[i]);
1161 
1162  B_cum[0]=1.0;
1163  for(int i=0;i<m;i++){
1164  IssmDouble cum_B=exp_B[0];
1165  for(int j=1;j<=i;j++) cum_B=cum_B*exp_B[j];
1166  B_cum[i+1]= cum_B;
1167  }
1168 
1169  // flux across grid cell boundaries
1170  Qs=xNew<IssmDouble>(m+1);
1171  for(int i=0;i<m+1;i++) Qs[i] = swf_ss * B_cum[i];
1172 
1173  // net energy flux to each grid cell
1174  for(int i=0;i<m;i++) swf[i] = (Qs[i]-Qs[i+1]);
1175 
1176  // add flux absorbed at surface
1177  swf[0] += swf_s;
1178 
1179  /*Free resources:*/
1180  xDelete<IssmDouble>(B_cum);
1181  xDelete<IssmDouble>(exp_B);
1182  xDelete<IssmDouble>(Qs);
1183  xDelete<IssmDouble>(B);
1184  }
1185  }
1186  /*Assign output pointers: */
1187  *pswf=swf;
1188 
1189 } /*}}}*/

◆ accumulation()

void 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 
)

Definition at line 1190 of file Gembx.cpp.

1190  { /*{{{*/
1191 
1192  // Adds precipitation and deposition to the model grid
1193 
1194  // Author: Alex Gardner, University of Alberta
1195  // Date last modified: JAN, 2008
1196 
1197  /* Description:
1198  adjusts the properties of the top grid cell to account for accumulation
1199  T_air & T = Air and top grid cell temperatures [K]
1200  Tmean = average surface temperature [K]
1201  Vmean = average wind velocity [m s-1]
1202  V = wind velocity [m s-1]
1203  C = average accumulation rate [kg m-2 yr-1]
1204  dz = topgrid cell length [m]
1205  d = density of top grid gell [kg m-3]
1206  P = precipitation [mm w.e.] or [kg m-3]
1207  re = effective grain radius [mm]
1208  gdn = grain dentricity
1209  gsp = grain sphericity*/
1210 
1211  // MAIN FUNCTION
1212  // specify constants
1213  IssmDouble dSnow = 150; // density of snow [kg m-3]
1214  const IssmDouble reNew = 0.05; // new snow grain size [mm]
1215  const IssmDouble gdnNew = 1.0; // new snow dendricity
1216  const IssmDouble gspNew = 0.5; // new snow sphericity
1217 
1218  /*intermediary: */
1219  IssmDouble* mInit=NULL;
1220  bool top=true;
1221  IssmDouble mass=0;
1222  IssmDouble massinit=0;
1223  IssmDouble mass_diff=0;
1224 
1225  /*output: */
1226  IssmDouble* T=NULL;
1227  IssmDouble* dz=NULL;
1228  IssmDouble* d=NULL;
1229  IssmDouble* W=NULL;
1230  IssmDouble* a=NULL;
1231  IssmDouble* re=NULL;
1232  IssmDouble* gdn=NULL;
1233  IssmDouble* gsp=NULL;
1234  int m=0;
1235 
1236  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" accumulation module\n");
1237 
1238  /*Recover pointers: */
1239  T=*pT;
1240  dz=*pdz;
1241  d=*pd;
1242  W=*pW;
1243  a=*pa;
1244  re=*pre;
1245  gdn=*pgdn;
1246  gsp=*pgsp;
1247  m=*pm;
1248 
1249  //Density of fresh snow [kg m-3]
1250  switch (dsnowIdx){
1251  case 0: // Default value defined above
1252  break;
1253 
1254  case 1: // Density of Antarctica snow
1255  dSnow = 350;
1256  break;
1257 
1258  case 2: // Density of Greenland snow, Fausto et al., 2008
1259  dSnow = 315;
1260  break;
1261 
1262  case 3: //Surface snow accumulation density from Kaspers et al., 2004, Antarctica
1263  //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W
1264  // 7.36x10-2 1.06x10-3 6.69x10-2 4.77x10-3
1265  dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
1266  break;
1267 
1268  case 4: // Kuipers Munneke and others (2015), Greenland
1269  dSnow = 481.0 + 4.834*(Tmean-CtoK);
1270  break;
1271  }
1272 
1273  // determine initial mass
1274  mInit=xNew<IssmDouble>(m);
1275  for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
1276  massinit=0.0;
1277  for(int i=0;i<m;i++)massinit+=mInit[i];
1278 
1279  if (P > 0.0+Ptol){
1280 
1281  if (T_air <= CtoK+Ttol){ // if snow
1282 
1283  IssmDouble z_snow = P/dSnow; // depth of snow
1284  IssmDouble dfall = gdnNew;
1285  IssmDouble sfall = gspNew;
1286  IssmDouble refall = reNew;
1287 
1288  //From Vionnet et al., 2012 (Crocus)
1289  dfall = min(max(1.29 - 0.17*V,0.20),1.0);
1290  sfall = min(max(0.08*V + 0.38,0.5),0.9);
1291  refall = max((0.1 + (1.0-dfall)*0.25 + (0.5-sfall)*0.1)/2.0,Gdntol);
1292 
1293  // if snow depth is greater than specified min dz, new cell created
1294  if (z_snow > dzMin+Dtol){
1295 
1296  newcell(&T,T_air,top,m); //new cell T
1297  newcell(&dz,z_snow,top,m); //new cell dz
1298  newcell(&d,dSnow,top,m); //new cell d
1299  newcell(&W,0.0,top,m); //new cell W
1300  newcell(&a,aSnow,top,m); //new cell a
1301  newcell(&re,refall,top,m); //new cell grain size
1302  newcell(&gdn,dfall,top,m); //new cell grain dendricity
1303  newcell(&gsp,sfall,top,m); //new cell grain sphericity
1304  m=m+1;
1305  }
1306  else { // if snow depth is less than specified minimum dz snow
1307 
1308  mass = mInit[0] + P; // grid cell adjust mass
1309 
1310  dz[0] = dz[0] + P/dSnow; // adjust grid cell depth
1311  d[0] = mass / dz[0]; // adjust grid cell density
1312 
1313  // adjust variables as a linearly weighted function of mass
1314  // adjust temperature (assume P is same temp as air)
1315  T[0] = (T_air * P + T[0] * mInit[0])/mass;
1316 
1317  // adjust a, re, gdn & gsp
1318  if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
1319  gdn[0] = (dfall * P + gdn[0] * mInit[0])/mass;
1320  gsp[0] = (sfall * P + gsp[0] * mInit[0])/mass;
1321  re[0] = max((0.1 + (1.0-gdn[0])*0.25 + (0.5-gsp[0])*0.1)/2.0,Gdntol);
1322  }
1323  }
1324  else{ // if rain
1325 
1326  /*rain is added by increasing the mass and temperature of the ice
1327  of the top grid cell. Temperatures are set artifically high to
1328  account for the latent heat of fusion. This is the same as
1329  directly adding liquid water to the the snow pack surface but
1330  makes the numerics easier.*/
1331 
1332  // grid cell adjust mass
1333  mass = mInit[0] + P;
1334 
1335  // adjust temperature
1336  // liquid: must account for latent heat of fusion
1337  T[0] = (P *(T_air + LF/CI) + T[0] * mInit[0]) / mass;
1338 
1339  // adjust grid cell density
1340  d[0] = mass / dz[0];
1341 
1342  // if d > the density of ice, d = dIce
1343  if (d[0] > dIce-Dtol){
1344  d[0] = dIce; // adjust d
1345  dz[0] = mass / d[0]; // dz is adjusted to conserve mass
1346  }
1347  }
1348 
1349  // check for conservation of mass
1350  mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i];
1351 
1352  mass_diff = mass - massinit - P;
1353 
1354  #ifndef _HAVE_AD_ //avoid round operation. only check in forward mode.
1355  mass_diff = round(mass_diff * 100.0)/100.0;
1356  if (mass_diff > 0) _error_("mass not conserved in accumulation function");
1357  #endif
1358 
1359  }
1360  /*Free resources:*/
1361  if(mInit)xDelete<IssmDouble>(mInit);
1362 
1363  /*Assign output pointers:*/
1364  *pT=T;
1365  *pdz=dz;
1366  *pd=d;
1367  *pW=W;
1368  *pa=a;
1369  *pre=re;
1370  *pgdn=gdn;
1371  *pgsp=gsp;
1372  *pm=m;
1373 } /*}}}*/

◆ melt()

void 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 
)

Definition at line 1374 of file Gembx.cpp.

1374  { /*{{{*/
1375 
1377 
1378  // Description:
1379  // computes the quantity of meltwater due to snow temperature in excess of
1380  // 0 deg C, determines pore water content and adjusts grid spacing
1381 
1382  /*intermediary:*/
1383  IssmDouble* m=NULL;
1384  IssmDouble* maxF=NULL;
1385  IssmDouble* dW=NULL;
1386  IssmDouble* exsW=NULL;
1387  IssmDouble* exsT=NULL;
1388  IssmDouble* surpT=NULL;
1389  IssmDouble* surpE=NULL;
1390  IssmDouble* flxDn=NULL;
1391  IssmDouble ER=0.0;
1392  IssmDouble* EI=NULL;
1393  IssmDouble* EW=NULL;
1394  IssmDouble* M=NULL;
1395  int* D=NULL;
1396 
1397  IssmDouble sumM=0.0;
1398  IssmDouble sumER=0.0;
1399  IssmDouble addE=0.0;
1400  IssmDouble mSum0=0.0;
1401  IssmDouble sumE0=0.0;
1402  IssmDouble mSum1=0.0;
1403  IssmDouble sumE1=0.0;
1404  IssmDouble dE=0.0;
1405  IssmDouble dm=0.0;
1406  IssmDouble X=0.0;
1407  IssmDouble Wi=0.0;
1408 
1409  IssmDouble Ztot=0.0;
1410  IssmDouble T_bot=0.0;
1411  IssmDouble m_bot=0.0;
1412  IssmDouble dz_bot=0.0;
1413  IssmDouble d_bot=0.0;
1414  IssmDouble W_bot=0.0;
1415  IssmDouble a_bot=0.0;
1416  IssmDouble re_bot=0.0;
1417  IssmDouble gdn_bot=0.0;
1418  IssmDouble gsp_bot=0.0;
1419  IssmDouble EI_bot=0.0;
1420  IssmDouble EW_bot=0.0;
1421  bool top=false;
1422 
1423  IssmDouble* Zcum=NULL;
1424  IssmDouble* dzMin2=NULL;
1425  IssmDouble* dzMax2=NULL;
1426  IssmDouble zY2=zY;
1427  bool lastCellFlag = false;
1428  int X1=0;
1429  int X2=0;
1430 
1431  int D_size = 0;
1432 
1433  /*outputs:*/
1434  IssmDouble Msurf = 0.0;
1435  IssmDouble mAdd = 0.0;
1436  IssmDouble surplusE = 0.0;
1437  IssmDouble surplusT = 0.0;
1438  IssmDouble dz_add = 0.0;
1439  IssmDouble Rsum = 0.0;
1440  IssmDouble* T=*pT;
1441  IssmDouble* d=*pd;
1442  IssmDouble* dz=*pdz;
1443  IssmDouble* W=*pW;
1444  IssmDouble* a=*pa;
1445  IssmDouble* re=*pre;
1446  IssmDouble* gdn=*pgdn;
1447  IssmDouble* gsp=*pgsp;
1448  int n=*pn;
1449  IssmDouble* R=NULL;
1450 
1451  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" melt module\n");
1452 
1454 
1455  /*Allocations: */
1456  M=xNewZeroInit<IssmDouble>(n);
1457  maxF=xNewZeroInit<IssmDouble>(n);
1458  dW=xNewZeroInit<IssmDouble>(n);
1459 
1460  // store initial mass [kg] and energy [J]
1461  m=xNew<IssmDouble>(n); for(int i=0;i<n;i++) m[i] = dz[i]* d[i]; // grid cell mass [kg]
1462  EI=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI; // initial enegy of snow/ice
1463  EW=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EW[i]= W[i] * (LF + CtoK * CI); // initial enegy of water
1464 
1465  mSum0 = cellsum(W,n) + cellsum(m,n); // total mass [kg]
1466  sumE0 = cellsum(EI,n) + cellsum(EW,n); // total energy [J]
1467 
1468  // initialize melt and runoff scalars
1469  Rsum = 0.0; // runoff [kg]
1470  sumM = 0.0; // total melt [kg]
1471  mAdd = 0.0; // mass added/removed to/from base of model [kg]
1472  addE = 0.0; // energy added/removed to/from base of model [J]
1473  dz_add=0.0; // thickness of the layer added/removed to/from base of model [m]
1474 
1475  // calculate temperature excess above 0 deg C
1476  exsT=xNewZeroInit<IssmDouble>(n);
1477  for(int i=0;i<n;i++) exsT[i]= max(0.0, T[i] - CtoK); // [K] to [degC]
1478 
1479  // new grid point center temperature, T [K]
1480  // for(int i=0;i<n;i++) T[i]-=exsT[i];
1481  for(int i=0;i<n;i++) T[i]=min(T[i],CtoK);
1482 
1483  // specify irreducible water content saturation [fraction]
1484  const IssmDouble Swi = 0.07; // assumed constant after Colbeck, 1974
1485  const IssmDouble dPHC = 830; //Pore closeoff density
1486 
1488  // check if any pore water
1489  if (cellsum(W,n) >0.0+Wtol){
1490  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" pore water refreeze\n");
1491  // calculate maximum freeze amount, maxF [kg]
1492  for(int i=0;i<n;i++) maxF[i] = max(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
1493 
1494  // freeze pore water and change snow/ice properties
1495  for(int i=0;i<n;i++) dW[i] = min(maxF[i], W[i]); // freeze mass [kg]
1496  for(int i=0;i<n;i++) W[i] -= dW[i]; // pore water mass [kg]
1497  for(int i=0;i<n;i++) m[i] += dW[i]; // new mass [kg]
1498  for(int i=0;i<n;i++) d[i] = m[i] / dz[i]; // density [kg m-3]
1499  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]
1500 
1501  // if pore water froze in ice then adjust d and dz thickness
1502  for(int i=0;i<n;i++)if(d[i]> dIce-Dtol)d[i]=dIce;
1503  for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
1504 
1505  }
1506 
1507  // squeeze water from snow pack
1508  exsW=xNew<IssmDouble>(n);
1509  for(int i=0;i<n;i++){
1510  Wi= (dIce - d[i]) * Swi * (m[i] / d[i]); // irreducible water content [kg]
1511  exsW[i] = max(0.0, W[i] - Wi); // water "squeezed" from snow [kg]
1512  }
1513 
1515 
1516  // run melt algorithm if there is melt water or excess pore water
1517  if ((cellsum(exsT,n) > 0.0+Ttol) || (cellsum(exsW,n) > 0.0+Wtol)){
1518  // _printf_(""MELT OCCURS");
1519  // check to see if thermal energy exceeds energy to melt entire cell
1520  // if so redistribute temperature to lower cells (temperature surplus)
1521  // (maximum T of snow before entire grid cell melts is a constant
1522  // LF/CI = 159.1342)
1523  surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = max(0.0, exsT[i]- LF/CI);
1524 
1525  if (cellsum(surpT,n) > 0.0+Ttol ){
1526  // _printf_("T Surplus");
1527  // calculate surplus energy
1528  surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];
1529 
1530  int i = 0;
1531  while (cellsum(surpE,n) > 0.0+Ttol && i<n){
1532 
1533  if (i<n-1){
1534  // use surplus energy to increase the temperature of lower cell
1535  T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
1536 
1537  exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
1538  T[i+1] = min(CtoK, T[i+1]);
1539 
1540  surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
1541  surpE[i+1] = surpT[i+1] * CI * m[i+1];
1542  }
1543  else{
1544  surplusT=max(0.0, exsT[i] - LF/CI);
1545  surplusE=surpE[i];
1546  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0){
1547  _printf0_(" WARNING: surplus energy at the base of GEMB column\n");
1548  }
1549  }
1550 
1551  // adjust current cell properties (again 159.1342 is the max T)
1552  exsT[i] = LF/CI;
1553  surpE[i] = 0.0;
1554  i = i + 1;
1555  }
1556  }
1557 
1558  // convert temperature excess to melt [kg]
1559  IssmDouble Mmax=0.0;
1560  for(int i=0;i<n;i++){
1561  Mmax=exsT[i] * d[i] * dz[i] * CI / LF;
1562  M[i] = min(Mmax, m[i]); // melt
1563  }
1564  Msurf = M[0];
1565  sumM = cellsum(M,n); // total melt [kg]
1566 
1567  // calculate maximum refreeze amount, maxF [kg]
1568  for(int i=0;i<n;i++)maxF[i] = max(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
1569 
1570  // initialize refreeze, runoff, flxDn and dW vectors [kg]
1571  IssmDouble* F = xNewZeroInit<IssmDouble>(n);
1572  IssmDouble* R = xNewZeroInit<IssmDouble>(n);
1573 
1574  for(int i=0;i<n;i++)dW[i] = 0.0;
1575  flxDn=xNewZeroInit<IssmDouble>(n+1);
1576 
1577  // determine the deepest grid cell where melt/pore water is generated
1578  X = 0;
1579  for(int i=n-1;i>=0;i--){
1580  if(M[i]> 0.0+Wtol || exsW[i]> 0.0+Wtol){
1581  X=i;
1582  break;
1583  }
1584  }
1585 
1586  IssmDouble depthice=0.0;
1588  for(int i=0;i<n;i++){
1589  // calculate total melt water entering cell
1590  IssmDouble inM = M[i]+ flxDn[i];
1591 
1592  depthice=0.0;
1593  if (d[i] >= dPHC-Dtol){
1594  for(int l=i;(l<n && d[l]>=dPHC-Dtol);l++) depthice=depthice+dz[l];
1595  }
1596 
1597  // break loop if there is no meltwater and if depth is > mw_depth
1598  if (fabs(inM) < Wtol && i > X){
1599  break;
1600  }
1601  // if reaches impermeable ice layer all liquid water runs off (R)
1602  else if (d[i] >= dIce-Dtol || (d[i] >= dPHC-Dtol && depthice>0.1)){ // dPHC = pore hole close off [kg m-3]
1603  // _printf_("ICE LAYER");
1604  // no water freezes in this cell
1605  // no water percolates to lower cell
1606  // cell ice temperature & density do not change
1607 
1608  m[i] = m[i] - M[i]; // mass after melt
1609  Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water
1610  dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water
1611  R[i] = max(0.0, inM - dW[i]); // runoff
1612  F[i] = 0.0;
1613  }
1614  // check if no energy to refreeze meltwater
1615  else if (fabs(maxF[i]) < Dtol){
1616  // _printf_("REFREEZE == 0");
1617  // no water freezes in this cell
1618  // cell ice temperature & density do not change
1619 
1620  m[i] = m[i] - M[i]; // mass after melt
1621  Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water
1622  dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water
1623  flxDn[i+1] = max(0.0, inM - dW[i]); // meltwater out
1624  R[i] = 0.0;
1625  F[i] = 0.0; // no freeze
1626  }
1627  // some or all meltwater refreezes
1628  else{
1629  // change in density and temperature
1630  // _printf_("MELT REFREEZE");
1631  //-----------------------melt water-----------------------------
1632  m[i] = m[i] - M[i];
1633  IssmDouble dz_0 = m[i]/d[i];
1634  IssmDouble dMax = (dIce - d[i])*dz_0; // d max = dIce
1635  IssmDouble F1 = min(min(inM,dMax),maxF[i]); // maximum refreeze
1636  m[i] = m[i] + F1; // mass after refreeze
1637  d[i] = m[i]/dz_0;
1638 
1639  //-----------------------pore water-----------------------------
1640  Wi = (dIce-d[i])* Swi * dz_0; // irreducible water
1641  dW[i] = min(inM - F1, Wi-W[i]); // change in pore water
1642  if (dW[i] < 0.0-Wtol && -1*dW[i]> W[i]-Wtol ){
1643  dW[i]= -1*W[i];
1644  }
1645  IssmDouble F2 = 0.0;
1646 
1647  if (dW[i] < 0.0-Wtol){ // excess pore water
1648  dMax = (dIce - d[i])*dz_0; // maximum refreeze
1649  IssmDouble maxF2 = min(dMax, maxF[i]-F1); // maximum refreeze
1650  F2 = min(-1*dW[i], maxF2); // pore water refreeze
1651  m[i] = m[i] + F2; // mass after refreeze
1652  d[i] = m[i]/dz_0;
1653  dW[i] = dW[i] - F2;
1654  }
1655 
1656  F[i] = F1 + F2;
1657 
1658  flxDn[i+1] = inM - F1 - dW[i]; // meltwater out
1659  if (m[i]>Wtol){
1660  T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
1661  }
1662 
1663  // check if an ice layer forms
1664  if (fabs(d[i] - dIce) < Dtol){
1665  // _printf_("ICE LAYER FORMS");
1666  // excess water runs off
1667  R[i] = flxDn[i+1];
1668  // no water percolates to lower cell
1669  flxDn[i+1] = 0.0;
1670  }
1671  }
1672  }
1673 
1675  for(int i=0;i<n;i++)if (W[i] < 0.0-Wtol) _error_("negative pore water generated in melt equations");
1676 
1677  // delete all cells with zero mass
1678  // adjust pore water
1679  for(int i=0;i<n;i++)W[i] += dW[i];
1680 
1681  //calculate Rsum:
1682  Rsum=cellsum(R,n) + flxDn[n];
1683 
1684  // delete all cells with zero mass
1685  D_size=0; for(int i=0;i<n;i++)if(m[i]> (0.0+Wtol))D_size++;
1686  D=xNew<int>(D_size);
1687  D_size=0; for(int i=0;i<n;i++)if(m[i]> (0.0+Wtol)){ D[D_size] = i; D_size++;}
1688 
1689  celldelete(&m,n,D,D_size);
1690  celldelete(&W,n,D,D_size);
1691  celldelete(&d,n,D,D_size);
1692  celldelete(&dz,n,D,D_size);
1693  celldelete(&T,n,D,D_size);
1694  celldelete(&a,n,D,D_size);
1695  celldelete(&re,n,D,D_size);
1696  celldelete(&gdn,n,D,D_size);
1697  celldelete(&gsp,n,D,D_size);
1698  celldelete(&EI,n,D,D_size);
1699  celldelete(&EW,n,D,D_size);
1700  n=D_size;
1701  xDelete<int>(D);
1702 
1703  // calculate new grid lengths
1704  for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
1705 
1706  /*Free resources:*/
1707  xDelete<IssmDouble>(F);
1708  xDelete<IssmDouble>(R);
1709  }
1710 
1711  //Merging of cells as they are burried under snow.
1712  Zcum=xNew<IssmDouble>(n);
1713  dzMin2=xNew<IssmDouble>(n);
1714  dzMax2=xNew<IssmDouble>(n);
1715 
1716  X=0;
1717  Zcum[0]=dz[0]; // Compute a cumulative depth vector
1718  for (int i=0;i<n;i++){
1719  if (i==0){
1720  dzMin2[i]=dzMin;
1721  }
1722  else{
1723  Zcum[i]=Zcum[i-1]+dz[i];
1724  if (Zcum[i]<=zTop+Dtol){
1725  dzMin2[i]=dzMin;
1726  X=i;
1727  }
1728  else{
1729  dzMin2[i]=zY2*dzMin2[i-1];
1730  }
1731  }
1732  }
1733 
1734  // Check to see if any cells are too small and need to be merged
1735  for (int i=0; i<n; i++){
1736  if ( (i<=X && dz[i]<dzMin-Dtol) || (i>X && dz[i]<dzMin2[i]-Dtol) ) {
1737 
1738  if (i==n-1){
1739  X2=i;
1740  //find closest cell to merge with
1741  for(int j=n-2;j>=0;j--){
1742  if(m[j]!=Delflag){
1743  X1=j;
1744  break;
1745  }
1746  }
1747  }
1748  else{
1749  X1=i+1;
1750  X2=i;
1751  }
1752 
1753  // adjust variables as a linearly weighted function of mass
1754  IssmDouble m_new = m[X2] + m[X1];
1755  T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
1756  a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
1757  re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
1758  gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
1759  gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
1760 
1761  // merge with underlying grid cell and delete old cell
1762  dz[X1] = dz[X2] + dz[X1]; // combine cell depths
1763  d[X1] = m_new / dz[X1]; // combine top densities
1764  W[X1] = W[X1] + W[X2]; // combine liquid water
1765  m[X1] = m_new; // combine top masses
1766 
1767  // set cell to -99999 for deletion
1768  m[X2] = Delflag;
1769  }
1770  }
1771 
1772  // delete combined cells
1773  D_size=0; for(int i=0;i<n;i++)if(m[i]> Delflag+Wtol)D_size++;
1774  D=xNew<int>(D_size);
1775  D_size=0; for(int i=0;i<n;i++)if(m[i]> Delflag+Wtol){ D[D_size] = i; D_size++;}
1776 
1777  celldelete(&m,n,D,D_size);
1778  celldelete(&W,n,D,D_size);
1779  celldelete(&dz,n,D,D_size);
1780  celldelete(&d,n,D,D_size);
1781  celldelete(&T,n,D,D_size);
1782  celldelete(&a,n,D,D_size);
1783  celldelete(&re,n,D,D_size);
1784  celldelete(&gdn,n,D,D_size);
1785  celldelete(&gsp,n,D,D_size);
1786  celldelete(&EI,n,D,D_size);
1787  celldelete(&EW,n,D,D_size);
1788  n=D_size;
1789  xDelete<int>(D);
1790 
1791  // check if any of the cell depths are too large
1792  X=0;
1793  Zcum[0]=dz[0]; // Compute a cumulative depth vector
1794  for (int i=0;i<n;i++){
1795  if (i==0){
1796  dzMax2[i]=dzMin*2;
1797  }
1798  else{
1799  Zcum[i]=Zcum[i-1]+dz[i];
1800  if (Zcum[i]<=zTop+Dtol){
1801  dzMax2[i]=dzMin*2;
1802  X=i;
1803  }
1804  else{
1805  dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
1806  }
1807  }
1808  }
1809 
1810  for (int j=n-1;j>=0;j--){
1811  if ((j<X && dz[j] > dzMax2[j]+Dtol) || (dz[j] > dzMax2[j]*zY2+Dtol)){
1812  // _printf_("dz > dzMin * 2");
1813  // split in two
1814  cellsplit(&dz, n, j,.5);
1815  cellsplit(&W, n, j,.5);
1816  cellsplit(&m, n, j,.5);
1817  cellsplit(&T, n, j,1.0);
1818  cellsplit(&d, n, j,1.0);
1819  cellsplit(&a, n, j,1.0);
1820  cellsplit(&EI, n, j,.5);
1821  cellsplit(&EW, n, j,.5);
1822  cellsplit(&re, n, j,1.0);
1823  cellsplit(&gdn, n, j,1.0);
1824  cellsplit(&gsp, n, j,1.0);
1825  n++;
1826  }
1827  }
1828 
1830  // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
1831  // INTERPRETATION
1832  // calculate total model depth
1833  Ztot = cellsum(dz,n);
1834 
1835  if (Ztot < zMin-Dtol){
1836  // printf("Total depth < zMin %f \n", Ztot);
1837  // mass and energy to be added
1838  mAdd = m[n-1]+W[n-1];
1839  addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);
1840 
1841  // add a grid cell of the same size and temperature to the bottom
1842  dz_bot=dz[n-1];
1843  T_bot=T[n-1];
1844  W_bot=W[n-1];
1845  m_bot=m[n-1];
1846  d_bot=d[n-1];
1847  a_bot=a[n-1];
1848  re_bot=re[n-1];
1849  gdn_bot=gdn[n-1];
1850  gsp_bot=gsp[n-1];
1851  EI_bot=EI[n-1];
1852  EW_bot=EW[n-1];
1853 
1854  dz_add=dz_bot;
1855 
1856  newcell(&dz,dz_bot,top,n);
1857  newcell(&T,T_bot,top,n);
1858  newcell(&W,W_bot,top,n);
1859  newcell(&m,m_bot,top,n);
1860  newcell(&d,d_bot,top,n);
1861  newcell(&a,a_bot,top,n);
1862  newcell(&re,re_bot,top,n);
1863  newcell(&gdn,gdn_bot,top,n);
1864  newcell(&gsp,gsp_bot,top,n);
1865  newcell(&EI,EI_bot,top,n);
1866  newcell(&EW,EW_bot,top,n);
1867  n=n+1;
1868  }
1869  else if (Ztot > zMax+Dtol){
1870  // printf("Total depth > zMax %f \n", Ztot);
1871  // mass and energy loss
1872  mAdd = -(m[n-1]+W[n-1]);
1873  addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
1874  dz_add=-(dz[n-1]);
1875 
1876  // remove a grid cell from the bottom
1877  D_size=n-1;
1878  D=xNew<int>(D_size);
1879 
1880  for(int i=0;i<n-1;i++) D[i]=i;
1881  celldelete(&dz,n,D,D_size);
1882  celldelete(&T,n,D,D_size);
1883  celldelete(&W,n,D,D_size);
1884  celldelete(&m,n,D,D_size);
1885  celldelete(&d,n,D,D_size);
1886  celldelete(&a,n,D,D_size);
1887  celldelete(&re,n,D,D_size);
1888  celldelete(&gdn,n,D,D_size);
1889  celldelete(&gsp,n,D,D_size);
1890  celldelete(&EI,n,D,D_size);
1891  celldelete(&EW,n,D,D_size);
1892  n=D_size;
1893  xDelete<int>(D);
1894  }
1895 
1897 
1898  // calculate final mass [kg] and energy [J]
1899  sumER = Rsum * (LF + CtoK * CI);
1900  for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
1901  for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
1902 
1903  mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
1904  sumE1 = cellsum(EI,n) + cellsum(EW,n);
1905 
1906  /*checks: */
1907  for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
1908 
1909  /*only in forward mode! avoid round in AD mode as it is not differentiable: */
1910  #ifndef _HAVE_AD_
1911  dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
1912  dE = round(sumE0 - sumE1 - sumER + addE - surplusE);
1913  if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
1914  << "dm: " << dm << " dE: " << dE << "\n");
1915  #endif
1916 
1917  /*Free resources:*/
1918  if(m)xDelete<IssmDouble>(m);
1919  if(EI)xDelete<IssmDouble>(EI);
1920  if(EW)xDelete<IssmDouble>(EW);
1921  if(maxF)xDelete<IssmDouble>(maxF);
1922  if(dW)xDelete<IssmDouble>(dW);
1923  if(exsW)xDelete<IssmDouble>(exsW);
1924  if(exsT)xDelete<IssmDouble>(exsT);
1925  if(surpT)xDelete<IssmDouble>(surpT);
1926  if(surpE)xDelete<IssmDouble>(surpE);
1927  if(flxDn)xDelete<IssmDouble>(flxDn);
1928  if(D)xDelete<int>(D);
1929  if(M)xDelete<IssmDouble>(M);
1930  xDelete<IssmDouble>(Zcum);
1931  xDelete<IssmDouble>(dzMin2);
1932 
1933  /*Assign output pointers:*/
1934  *pMs=Msurf;
1935  *pM=sumM;
1936  *pR=Rsum;
1937  *pmAdd=mAdd;
1938  *pdz_add=dz_add;
1939 
1940  *pT=T;
1941  *pd=d;
1942  *pdz=dz;
1943  *pW=W;
1944  *pa=a;
1945  *pre=re;
1946  *pgdn=gdn;
1947  *pgsp=gsp;
1948  *pn=n;
1949 
1950 } /*}}}*/

◆ densification()

void densification ( IssmDouble **  pd,
IssmDouble **  pdz,
IssmDouble *  T,
IssmDouble *  re,
int  denIdx,
IssmDouble  C,
IssmDouble  dt,
IssmDouble  Tmean,
IssmDouble  dIce,
int  m,
int  sid 
)

Definition at line 1951 of file Gembx.cpp.

1951  { /*{{{*/
1952 
1954 
1956 
1957  // Author: Alex Gardner, University of Alberta
1958  // Date last modified: FEB, 2008
1959 
1960  // Description:
1961  // computes the densification of snow/firn using the emperical model of
1962  // Herron and Langway (1980) or the semi-emperical model of Anthern et al.
1963  // (2010)
1964 
1965  // Inputs:
1966  // denIdx = densification model to use:
1967  // 1 = emperical model of Herron and Langway (1980)
1968  // 2 = semi-imerical model of Anthern et al. (2010)
1969  // 3 = physical model from Appendix B of Anthern et al. (2010)
1970  // d = initial snow/firn density [kg m-3]
1971  // T = temperature [K]
1972  // dz = grid cell size [m]
1973  // C = average accumulation rate [kg m-2 yr-1]
1974  // dt = time lapsed [s]
1975  // re = effective grain radius [mm];
1976  // Ta = mean annual temperature
1977 
1978  // Reference:
1979  // Herron and Langway (1980), Anthern et al. (2010)
1980 
1982  // denIdx = 2;
1983  // d = 800;
1984  // T = 270;
1985  // dz = 0.005;
1986  // C = 200;
1987  // dt = 60*60;
1988  // re = 0.7;
1989  // Tmean = 273.15-18;
1990 
1992  // specify constants
1993  dt = dt / dts; // convert from [s] to [d]
1994  // R = 8.314 // gas constant [mol-1 K-1]
1995  // Ec = 60 // activation energy for self-diffusion of water
1996  // // molecules through the ice tattice [kJ mol-1]
1997  // Eg = 42.4 // activation energy for grain growth [kJ mol-1]
1998 
1999  /*intermediary: */
2000  IssmDouble c0=0.0;
2001  IssmDouble c1=0.0;
2002  IssmDouble H=0.0;
2003  IssmDouble M0=0.0;
2004  IssmDouble M1=0.0;
2005  IssmDouble c0arth=0.0;
2006  IssmDouble c1arth=0.0;
2007 
2008  /*output: */
2009  IssmDouble* dz=NULL;
2010  IssmDouble* d=NULL;
2011 
2012  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" densification module\n");
2013 
2014  /*Recover pointers: */
2015  dz=*pdz;
2016  d=*pd;
2017 
2018  // initial mass
2019  IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
2020 
2021  /*allocations and initialization of overburden pressure and factor H: */
2022  IssmDouble* cumdz = xNew<IssmDouble>(m-1);
2023  cumdz[0]=dz[0];
2024  for(int i=1;i<m-1;i++)cumdz[i]=cumdz[i-1]+dz[i];
2025 
2026  IssmDouble* obp = xNew<IssmDouble>(m);
2027  obp[0]=0.0;
2028  for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
2029 
2030  // calculate new snow/firn density for:
2031  // snow with densities <= 550 [kg m-3]
2032  // snow with densities > 550 [kg m-3]
2033 
2034  for(int i=0;i<m;i++){
2035  switch (denIdx){
2036  case 1: // Herron and Langway (1980)
2037  c0 = (11.0 * exp(-10160.0 / (T[i] * R))) * C/1000.0;
2038  c1 = (575.0 * exp(-21400.0 / (T[i]* R))) * pow(C/1000.0,.5);
2039  break;
2040 
2041  case 2: // Arthern et al. (2010) [semi-emperical]
2042  // common variable
2043  // NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
2044  H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
2045  c0 = 0.07 * H;
2046  c1 = 0.03 * H;
2047  break;
2048 
2049  case 3: // Arthern et al. (2010) [physical model eqn. B1]
2050 
2051  // common variable
2052  H = exp((-60000.0/(T[i] * R))) * obp[i] / pow(re[i]/1000.0,2.0);
2053  c0 = 9.2e-9 * H;
2054  c1 = 3.7e-9 * H;
2055  break;
2056 
2057  case 4: // Li and Zwally (2004)
2058  c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(CtoK - T[i],-2.061);
2059  c1 = c0;
2060  break;
2061 
2062  case 5: // Helsen et al. (2008)
2063  // common variable
2064  c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
2065  c1 = c0;
2066  break;
2067 
2068  case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica
2069  // common variable
2070  // From literature: H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
2071  H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
2072  c0arth = 0.07 * H;
2073  c1arth = 0.03 * H;
2074  //ERA-5
2075  //M0 = max(2.3128 - (0.2480 * log(C)),0.25);
2076  //M1 = max(2.7950 - (0.3318 * log(C)),0.25);
2077  //RACMO
2078  M0 = max(1.6599 - (0.1724 * log(C)),0.25);
2079  M1 = max(2.0102 - (0.2458 * log(C)),0.25);
2080  //From Ligtenberg
2081  //H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
2082  //M0 = max(1.435 - (0.151 * log(C)),0.25);
2083  //M1 = max(2.366 - (0.293 * log(C)),0.25);
2084  c0 = M0*c0arth;
2085  c1 = M1*c1arth;
2086  break;
2087 
2088  case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
2089  // common variable
2090  // From literature: H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
2091  H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
2092  c0arth = 0.07 * H;
2093  c1arth = 0.03 * H;
2094  // ERA5
2095  //M0 = max(1.8554 - (0.1316 * log(C)),0.25);
2096  //M1 = max(2.8901 - (0.3014 * log(C)),0.25);
2097  // RACMO
2098  M0 = max(1.6201 - (0.1450 * log(C)),0.25);
2099  M1 = max(2.5577 - (0.2899 * log(C)),0.25);
2100  // From Kuipers Munneke
2101  //M0 = max(1.042 - (0.0916 * log(C)),0.25);
2102  //M1 = max(1.734 - (0.2039 * log(C)),0.25);
2103  c0 = M0*c0arth;
2104  c1 = M1*c1arth;
2105  break;
2106  }
2107 
2108  // new snow density
2109  if(d[i] <= 550.0+Dtol) d[i] = d[i] + (c0 * (dIce - d[i]) / 365.0 * dt);
2110  else d[i] = d[i] + (c1 * (dIce - d[i]) / 365.0 * dt);
2111 
2112  //disp((num2str(nanmean(c0 .* (dIce - d(idx)) / 365 * dt))))
2113 
2114  // do not allow densities to exceed the density of ice
2115  if(d[i] > dIce-Ptol) d[i]=dIce;
2116 
2117  // calculate new grid cell length
2118  dz[i] = mass_init[i] / d[i];
2119  }
2120 
2121  /*Free resources:*/
2122  xDelete<IssmDouble>(mass_init);
2123  xDelete<IssmDouble>(cumdz);
2124  xDelete<IssmDouble>(obp);
2125 
2126  /*Assign output pointers:*/
2127  *pdz=dz;
2128  *pd=d;
2129 
2130 } /*}}}*/

◆ turbulentFlux()

void 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 
)

Definition at line 2131 of file Gembx.cpp.

2131  { /*{{{*/
2132 
2134 
2135  // Description:
2136  // computed the surface sensible and latent heat fluxes [W m-2], this
2137  // function also calculated the mass loss/acreation due to
2138  // condensation/evaporation [kg]
2139 
2140  // Reference:
2141  // Dingman, 2002.
2142 
2144  // Ta: 2m air temperature [K]
2145  // Ts: snow/firn/ice surface temperature [K]
2146  // V: wind speed [m s^-^1]
2147  // eAir: screen level vapor pressure [Pa]
2148  // pAir: surface pressure [Pa]
2149  // ds: surface density [kg/m^3]
2150  // Ws: surface liquid water content [kg/m^2]
2151  // Vz: height above ground at which wind (V) eas sampled [m]
2152  // Tz: height above ground at which temperature (T) was sampled [m]
2153 
2154  /*intermediary:*/
2155  IssmDouble d_air=0.0;
2156  IssmDouble Ri=0.0;
2157  IssmDouble z0=0.0;
2158  IssmDouble coef_M=0.0;
2159  IssmDouble coef_H=0.0;
2160  IssmDouble An=0.0;
2161  IssmDouble C=0.0;
2162  IssmDouble L=0.0;
2163  IssmDouble eS=0.0;
2164 
2165  /*output: */
2166  IssmDouble shf=0.0;
2167  IssmDouble lhf=0.0;
2168  IssmDouble EC=0.0;
2169 
2170  if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" turbulentFlux module\n");
2171 
2172  // calculated air density [kg/m3]
2173  d_air = 0.029 * pAir /(R * Ta);
2174 
2176  // Bougamont, 2006
2177  // wind/temperature surface roughness height [m]
2178  if (ds < dIce-Dtol && fabs(Ws) < Wtol) z0 = 0.00012; // 0.12 mm for dry snow
2179  else if (ds >= dIce-Dtol) z0 = 0.0032; // 3.2 mm for ice
2180  else z0 = 0.0013; // 1.3 mm for wet snow
2181 
2183  // Reference:
2184  // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
2185  // Journal of Climatology, 2, 65-84.
2186 
2187  // if V = 0 goes to infinity therfore if V = 0 change
2188  if(V< 0.01-Dtol)V=0.01;
2189 
2190  // calculate the Bulk Richardson Number (Ri)
2191  Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
2192 
2193  // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
2194 
2195  // do not allow Ri to exceed 0.19
2196  Ri = min(Ri, 0.19);
2197 
2198  // calculate momentum 'coefM' stability factor
2199  if (Ri > 0.0+Ttol){
2200  // if stable
2201  coef_M = 1.0/(1.0-5.2*Ri);
2202  }
2203  else {
2204  coef_M =pow (1.0-18*Ri,-0.25);
2205  }
2206 
2207  // calculate heat/wind 'coef_H' stability factor
2208  if (Ri <= -0.03+Ttol) coef_H = coef_M/1.3;
2209  else coef_H = coef_M;
2210 
2212  //An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient
2213  An = pow(0.4,2) / (log(Tz/z0)*log(Vz/z0)); // Bulk-transfer coefficient
2214  C = An * d_air * V; // shf & lhf common coefficient
2215 
2217  // calculate the sensible heat flux [W m-2](Patterson, 1998)
2218  shf = C * CA * (Ta - Ts);
2219 
2220  // adjust using Monin-Obukhov stability theory
2221  shf = shf / (coef_M * coef_H);
2222 
2223  // Latent Heat
2224  // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
2225  if (Ts >= CtoK-Ttol){
2226  L = LV; //for liquid water at 273.15 k to vapor
2227 
2228  //for liquid surface (assume liquid on surface when Ts == 0 deg C)
2229  // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
2230  eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
2231  }
2232  else{
2233  L = LS; // latent heat of sublimation
2234 
2235  // for an ice surface Murphy and Koop, 2005 [Equation 7]
2236  eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
2237  }
2238 
2239  // Latent heat flux [W m-2]
2240  lhf = C * L * (eAir - eS) * 0.622 / pAir;
2241 
2242  // adjust using Monin-Obukhov stability theory (if lhf '+' then there is
2243  // energy and mass gained at the surface, if '-' then there is mass and
2244  // energy loss at the surface.
2245  lhf = lhf / (coef_M * coef_H);
2246 
2247  // mass loss (-)/acreation(+) due to evaporation/condensation [kg]
2248  EC = lhf * dts / L;
2249 
2250  /*assign output poitners: */
2251  *pshf=shf;
2252  *plhf=lhf;
2253  *pEC=EC;
2254 
2255 } /*}}}*/

Variable Documentation

◆ Pi

const double Pi = 3.141592653589793

Definition at line 11 of file Gembx.cpp.

◆ CtoK

const double CtoK = 273.15

Definition at line 12 of file Gembx.cpp.

◆ dts

const double dts = 86400.0

Definition at line 13 of file Gembx.cpp.

◆ Ttol

const double Ttol = 1e-10

Definition at line 18 of file Gembx.cpp.

◆ Dtol

const double Dtol = 1e-11

Definition at line 19 of file Gembx.cpp.

◆ Gdntol

const double Gdntol = 1e-10

Definition at line 20 of file Gembx.cpp.

◆ Wtol

const double Wtol = 1e-13

Definition at line 21 of file Gembx.cpp.

◆ Ptol

const double Ptol = 1e-6

Definition at line 22 of file Gembx.cpp.

◆ CI

const double CI = 2102.0

Definition at line 24 of file Gembx.cpp.

◆ LF

const double LF = 0.3345E6

Definition at line 25 of file Gembx.cpp.

◆ LV

const double LV = 2.495E6

Definition at line 26 of file Gembx.cpp.

◆ LS

const double LS = 2.8295E6

Definition at line 27 of file Gembx.cpp.

◆ SB

const double SB = 5.67E-8

Definition at line 28 of file Gembx.cpp.

◆ CA

const double CA = 1005.0

Definition at line 29 of file Gembx.cpp.

◆ R

const double R = 8.314

Definition at line 30 of file Gembx.cpp.

◆ Delflag

const double Delflag =-99999

Definition at line 32 of file Gembx.cpp.

DataSet::Size
int Size()
Definition: DataSet.cpp:399
TimesteppingFinalTimeEnum
@ TimesteppingFinalTimeEnum
Definition: EnumDefinitions.h:430
Delflag
const double Delflag
Definition: Gembx.cpp:32
newcell
void newcell(IssmDouble **pcell, IssmDouble newvalue, bool top, int m)
Definition: MatrixUtils.cpp:577
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
SmbDtEnum
@ SmbDtEnum
Definition: EnumDefinitions.h:357
IssmDouble
double IssmDouble
Definition: types.h:37
xIsNan
int xIsNan(const T &X)
Definition: isnan.h:16
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
TransientInput2::GetTimeInputOffset
int GetTimeInputOffset(IssmDouble time)
Definition: TransientInput2.cpp:545
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Gdntol
const double Gdntol
Definition: Gembx.cpp:20
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
TransientInput2
Definition: TransientInput2.h:13
CtoK
const double CtoK
Definition: Gembx.cpp:12
SmbTaEnum
@ SmbTaEnum
Definition: EnumDefinitions.h:782
CI
const double CI
Definition: Gembx.cpp:24
Element
Definition: Element.h:41
Element::SmbGemb
void SmbGemb(IssmDouble timeinputs, int count)
Definition: Element.cpp:3584
Element::inputs2
Inputs2 * inputs2
Definition: Element.h:47
TransientInput2::GetTimeByOffset
IssmDouble GetTimeByOffset(int offset)
Definition: TransientInput2.cpp:539
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
VerboseSmb
bool VerboseSmb(void)
Definition: Verbosity.cpp:30
SmbIsclimatologyEnum
@ SmbIsclimatologyEnum
Definition: EnumDefinitions.h:363
Wtol
const double Wtol
Definition: Gembx.cpp:21
cellsum
IssmDouble cellsum(IssmDouble *cell, int m)
Definition: MatrixUtils.cpp:604
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Ptol
const double Ptol
Definition: Gembx.cpp:22
TimesteppingInterpForcingsEnum
@ TimesteppingInterpForcingsEnum
Definition: EnumDefinitions.h:431
Inputs2::GetTransientInput
TransientInput2 * GetTransientInput(int enum_type)
Definition: Inputs2.cpp:406
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
Ttol
const double Ttol
Definition: Gembx.cpp:18
R
const double R
Definition: Gembx.cpp:30
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
cellsplit
void cellsplit(IssmDouble **pcell, int m, int i, IssmDouble scale)
Definition: MatrixUtils.cpp:630
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
SB
const double SB
Definition: Gembx.cpp:28
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
dts
const double dts
Definition: Gembx.cpp:13
Dtol
const double Dtol
Definition: Gembx.cpp:19
LS
const double LS
Definition: Gembx.cpp:27
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
Marbouty
IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT)
Definition: Gembx.cpp:166
CA
const double CA
Definition: Gembx.cpp:29
celldelete
void celldelete(IssmDouble **pcell, int m, int *indices, int nind)
Definition: MatrixUtils.cpp:612
Pi
const double Pi
Definition: Gembx.cpp:11
LF
const double LF
Definition: Gembx.cpp:25
LV
const double LV
Definition: Gembx.cpp:26
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:24 for Ice Sheet System Model by   doxygen 1.8.19