Ice Sheet System Model  4.18
Code documentation
Functions
SurfaceMassBalancex.h File Reference

header file for SMB More...

#include "../../classes/classes.h"

Go to the source code of this file.

Functions

void SurfaceMassBalancex (FemModel *femmodel)
 
void SmbForcingx (FemModel *femmodel)
 
void SmbGradientsx (FemModel *femmodel)
 
void SmbGradientsElax (FemModel *femmodel)
 
void Delta18oParameterizationx (FemModel *femmodel)
 
void MungsmtpParameterizationx (FemModel *femmodel)
 
void Delta18opdParameterizationx (FemModel *femmodel)
 
void PositiveDegreeDayx (FemModel *femmodel)
 
void PositiveDegreeDaySicopolisx (FemModel *femmodel)
 
void SmbHenningx (FemModel *femmodel)
 
void SmbComponentsx (FemModel *femmodel)
 
void SmbMeltComponentsx (FemModel *femmodel)
 
void SmbGradientsComponentsx (FemModel *femmodel)
 
void SmbSemicx (FemModel *femmodel)
 
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 **a, int aIdx, IssmDouble *re, IssmDouble *dz, IssmDouble *d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble *T, IssmDouble *W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid)
 
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 thermo (IssmDouble *pEC, IssmDouble **T, IssmDouble *pulwrf, IssmDouble *dz, IssmDouble *d, IssmDouble *swf, IssmDouble dlw, 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 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 Ta, 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)
 

Detailed Description

header file for SMB

Definition in file SurfaceMassBalancex.h.

Function Documentation

◆ SurfaceMassBalancex()

void SurfaceMassBalancex ( FemModel *  femmodel)

◆ SmbForcingx()

void SmbForcingx ( FemModel *  femmodel)

Definition at line 11 of file SurfaceMassBalancex.cpp.

11  {/*{{{*/
12 
13  // void SmbForcingx(smb,ni){
14  // INPUT parameters: ni: working size of arrays
15  // OUTPUT: mass-balance (m/yr ice): agd(NA)
16  bool isclimatology;
17  femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
18 
19  if (isclimatology){
20 
21  /*Get time parameters*/
22  IssmDouble time,dt,starttime,finaltime;
23  femmodel->parameters->FindParam(&time,TimeEnum);
24  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
25  femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
26  femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
27 
28  if(time<=starttime+dt){
29  /*FIXME: this is wrong, should be done at the ElementUpdate step of analysis, not here!*/
30  InputDuplicatex(femmodel,SmbMassBalanceEnum,SmbMassBalanceClimateEnum);
31  femmodel->inputs2->DeleteInput(SmbMassBalanceEnum);
32  }
33 
34  /*If this is a climatology, we need to repeat the forcing after the final time*/
35  TransientInput2* smb_input=femmodel->inputs2->GetTransientInput(SmbMassBalanceClimateEnum); _assert_(smb_input);
36 
37  /*Get accumulation climatology value*/
38  int offsetend = smb_input->GetTimeInputOffset(finaltime);
39  IssmDouble time0 = smb_input->GetTimeByOffset(-1);
40  IssmDouble timeend = smb_input->GetTimeByOffset(offsetend);
41 
42  _assert_(timeend>time0);
43  IssmDouble timeclim = time;
44 
45  if(time>time0 && timeend>time0){
46  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
47  if(delta==0.){
48  timeclim=timeend;
49  }
50  else{
51  timeclim=time0+delta;
52  }
53  }
54 
55  /*Loop over all the elements of this partition*/
56  for(int i=0;i<femmodel->elements->Size();i++){
57  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
58 
59  int numvertices = element->GetNumberOfVertices();
60  IssmDouble* smb = xNew<IssmDouble>(numvertices);
61  element->GetInputListOnVerticesAtTime(smb,SmbMassBalanceClimateEnum,timeclim);
62 
63  /*Add input to element and Free memory*/
64  element->AddInput2(SmbMassBalanceEnum,smb,P1Enum);
65  xDelete<IssmDouble>(smb);
66  }
67  }
68 
69 }/*}}}*/

◆ SmbGradientsx()

void SmbGradientsx ( FemModel *  femmodel)

Definition at line 70 of file SurfaceMassBalancex.cpp.

70  {/*{{{*/
71 
72  // void SurfaceMassBalancex(hd,agd,ni){
73  // INPUT parameters: ni: working size of arrays
74  // INPUT: surface elevation (m): hd(NA)
75  // OUTPUT: mass-balance (m/yr ice): agd(NA)
76  int v;
77  IssmDouble rho_water; // density of fresh water
78  IssmDouble rho_ice; // density of ice
79  IssmDouble yts; // conversion factor year to second
80 
81  /*Loop over all the elements of this partition*/
82  for(int i=0;i<femmodel->elements->Size();i++){
83  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
84 
85  /*Allocate all arrays*/
86  int numvertices = element->GetNumberOfVertices();
87  IssmDouble* Href = xNew<IssmDouble>(numvertices); // reference elevation from which deviations are used to calculate the SMB adjustment
88  IssmDouble* Smbref = xNew<IssmDouble>(numvertices); // reference SMB to which deviations are added
89  IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // Hs-SMB relation parameter
90  IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // Hs-SMB relation paremeter
91  IssmDouble* s = xNew<IssmDouble>(numvertices); // surface elevation (m)
92  IssmDouble* smb = xNew<IssmDouble>(numvertices);
93 
94  /*Recover SmbGradients*/
95  element->GetInputListOnVertices(Href,SmbHrefEnum);
96  element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
97  element->GetInputListOnVertices(b_pos,SmbBPosEnum);
98  element->GetInputListOnVertices(b_neg,SmbBNegEnum);
99 
100  /*Recover surface elevation at vertices: */
101  element->GetInputListOnVertices(s,SurfaceEnum);
102 
103  /*Get material parameters :*/
104  rho_ice=element->FindParam(MaterialsRhoIceEnum);
105  rho_water=element->FindParam(MaterialsRhoFreshwaterEnum);
106 
107  /* Get constants */
108  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
109 
110  // loop over all vertices
111  for(v=0;v<numvertices;v++){
112  if(Smbref[v]>0){
113  smb[v]=Smbref[v]+b_pos[v]*(s[v]-Href[v]);
114  }
115  else{
116  smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
117  }
118 
119  smb[v]=smb[v]/1000*rho_water/rho_ice; // SMB in m/y ice
120  } //end of the loop over the vertices
121 
122  /*Add input to element and Free memory*/
123  element->AddInput2(SmbMassBalanceEnum,smb,P1Enum);
124  xDelete<IssmDouble>(Href);
125  xDelete<IssmDouble>(Smbref);
126  xDelete<IssmDouble>(b_pos);
127  xDelete<IssmDouble>(b_neg);
128  xDelete<IssmDouble>(s);
129  xDelete<IssmDouble>(smb);
130  }
131 
132 }/*}}}*/

◆ SmbGradientsElax()

void SmbGradientsElax ( FemModel *  femmodel)

Definition at line 133 of file SurfaceMassBalancex.cpp.

133  {/*{{{*/
134 
135  // void SurfaceMassBalancex(hd,agd,ni){
136  // INPUT parameters: ni: working size of arrays
137  // INPUT: surface elevation (m): hd(NA)
138  // OUTPUT: surface mass-balance (m/yr ice): agd(NA)
139  int v;
140 
141  /*Loop over all the elements of this partition*/
142  for(int i=0;i<femmodel->elements->Size();i++){
143  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
144 
145  /*Allocate all arrays*/
146  int numvertices = element->GetNumberOfVertices();
147  IssmDouble* ela = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
148  IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
149  IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
150  IssmDouble* b_max = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
151  IssmDouble* b_min = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
152  IssmDouble* s = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
153  IssmDouble* smb = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
154 
155  /*Recover ELA, SMB gradients, and caps*/
156  element->GetInputListOnVertices(ela,SmbElaEnum);
157  element->GetInputListOnVertices(b_pos,SmbBPosEnum);
158  element->GetInputListOnVertices(b_neg,SmbBNegEnum);
159  element->GetInputListOnVertices(b_max,SmbBMaxEnum);
160  element->GetInputListOnVertices(b_min,SmbBMinEnum);
161 
162  /*Recover surface elevation at vertices: */
163  element->GetInputListOnVertices(s,SurfaceEnum);
164 
165  /*Loop over all vertices, calculate SMB*/
166  for(v=0;v<numvertices;v++){
167  // if surface is above the ELA
168  if(s[v]>ela[v]){
169  smb[v]=b_pos[v]*(s[v]-ela[v]);
170  }
171  // if surface is below or equal to the ELA
172  else{
173  smb[v]=b_neg[v]*(s[v]-ela[v]);
174  }
175 
176  // if SMB is larger than upper cap, set SMB to upper cap
177  if(smb[v]>b_max[v]){
178  smb[v]=b_max[v];
179  }
180  // if SMB is smaller than lower cap, set SMB to lower cap
181  if(smb[v]<b_min[v]){
182  smb[v]=b_min[v];
183  }
184  } //end of the loop over the vertices
185 
186  /*Add input to element and Free memory*/
187  element->AddInput2(SmbMassBalanceEnum,smb,P1Enum);
188  xDelete<IssmDouble>(ela);
189  xDelete<IssmDouble>(b_pos);
190  xDelete<IssmDouble>(b_neg);
191  xDelete<IssmDouble>(b_max);
192  xDelete<IssmDouble>(b_min);
193  xDelete<IssmDouble>(s);
194  xDelete<IssmDouble>(smb);
195 
196  }
197 
198 }/*}}}*/

◆ Delta18oParameterizationx()

void Delta18oParameterizationx ( FemModel *  femmodel)

Definition at line 199 of file SurfaceMassBalancex.cpp.

199  {/*{{{*/
200 
201  for(int i=0;i<femmodel->elements->Size();i++){
202  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
203  element->Delta18oParameterization();
204  }
205 
206 }/*}}}*/

◆ MungsmtpParameterizationx()

void MungsmtpParameterizationx ( FemModel *  femmodel)

Definition at line 207 of file SurfaceMassBalancex.cpp.

207  {/*{{{*/
208 
209  for(int i=0;i<femmodel->elements->Size();i++){
210  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
211  element->MungsmtpParameterization();
212  }
213 
214 }/*}}}*/

◆ Delta18opdParameterizationx()

void Delta18opdParameterizationx ( FemModel *  femmodel)

Definition at line 215 of file SurfaceMassBalancex.cpp.

215  {/*{{{*/
216 
217  for(int i=0;i<femmodel->elements->Size();i++){
218  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
219  element->Delta18opdParameterization();
220  }
221 
222 }/*}}}*/

◆ PositiveDegreeDayx()

void PositiveDegreeDayx ( FemModel *  femmodel)

Definition at line 223 of file SurfaceMassBalancex.cpp.

223  {/*{{{*/
224 
225  // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
226  // note "v" prefix means 12 monthly means, ie time dimension
227  // INPUT parameters: ni: working size of arrays
228  // INPUT: surface elevation (m): hd(NA)
229  // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
230  // ,NTIME)
231  // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
232  // OUTPUT: mass-balance (m/yr ice): agd(NA)
233  // mean annual surface temperature (degrees C): Tsurf(NA)
234 
235  int i, it, jj, itm;
236  IssmDouble DT = 0.02, sigfac, snormfac;
237  IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
238  IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
239  IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
240  IssmDouble siglimc, siglim0, siglim0c;
241  IssmDouble tstep, tsint, tint, tstepc;
242  int NPDMAX = 1504, NPDCMAX = 1454;
243  //IssmDouble pdds[NPDMAX]={0};
244  //IssmDouble pds[NPDCMAX]={0};
245  IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
246  IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
247  IssmDouble tstar; // monthly mean surface temp
248 
249  bool ismungsm;
250  bool issetpddfac;
251 
252  IssmDouble *pdds = NULL;
253  IssmDouble *pds = NULL;
254  Element *element = NULL;
255 
256  pdds=xNew<IssmDouble>(NPDMAX+1);
257  pds=xNew<IssmDouble>(NPDCMAX+1);
258 
259  // Get ismungsm parameter
260  femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
261 
262  // Get issetpddfac parameter
263  femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
264 
265  /* initialize PDD (creation of a lookup table)*/
266  tstep = 0.1;
267  tsint = tstep*0.5;
268  sigfac = -1.0/(2.0*pow(signorm,2));
269  snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
270  siglim = 2.5*signorm;
271  siglimc = 2.5*signormc;
272  siglim0 = siglim/DT + 0.5;
273  siglim0c = siglimc/DT + 0.5;
274  PDup = siglimc+PDCUT;
275 
276  itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
277 
278  if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
279  for(it = 0; it < itm; it++){
280  // tstar = REAL(it)*DT-siglim;
281  tstar = it*DT-siglim;
282  tint = tsint;
283  pddt = 0.;
284  for ( jj = 0; jj < 600; jj++){
285  if (tint > (tstar+siglim)){break;}
286  pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
287  tint = tint+tstep;
288  }
289  pdds[it] = pddt*snormfac;
290  }
291  pdds[itm+1] = siglim + DT;
292 
293  //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
294  tstepc = 0.1;
295  tsint = PDCUT-tstepc*0.5;
296  signormc = signorm - 0.5;
297  sigfac = -1.0/(2.0*pow(signormc,2));
298  snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
299  siglimc = 2.5*signormc ;
300  itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
301  if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
302  for(it = 0; it < itm; it++ ){
303  tstar = it*DT-siglimc;
304  // tstar = REAL(it)*DT-siglimc;
305  tint = tsint; // start against upper bound
306  pd = 0.;
307  for (jj = 0; jj < 600; jj++){
308  if (tint<(tstar-siglimc)) {break;}
309  pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
310  tint = tint-tstepc;
311  }
312  pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
313  }
314  pds[itm+1] = 0.;
315  // *******END initialize PDD
316 
317  for(i=0;i<femmodel->elements->Size();i++){
318  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
319  element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
320  }
321  /*free ressouces: */
322  xDelete<IssmDouble>(pdds);
323  xDelete<IssmDouble>(pds);
324 }/*}}}*/

◆ PositiveDegreeDaySicopolisx()

void PositiveDegreeDaySicopolisx ( FemModel *  femmodel)

Definition at line 325 of file SurfaceMassBalancex.cpp.

325  {/*{{{*/
326 
327  bool isfirnwarming;
328  femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
329 
330  for(int i=0;i<femmodel->elements->Size();i++){
331  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
332  element->PositiveDegreeDaySicopolis(isfirnwarming);
333  }
334 
335 }/*}}}*/

◆ SmbHenningx()

void SmbHenningx ( FemModel *  femmodel)

Definition at line 336 of file SurfaceMassBalancex.cpp.

336  {/*{{{*/
337 
338  /*Intermediaries*/
339  IssmDouble z_critical = 1675.;
340  IssmDouble dz = 0;
341  IssmDouble a = -15.86;
342  IssmDouble b = 0.00969;
343  IssmDouble c = -0.235;
344  IssmDouble f = 1.;
345  IssmDouble g = -0.0011;
346  IssmDouble h = -1.54e-5;
347  IssmDouble smb,smbref,anomaly,yts,z;
348 
349  /* Get constants */
350  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
351  /*iomodel->FindConstant(&yts,"md.constants.yts");*/
352  /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
353  /*Mathieu original*/
354  /*IssmDouble smb,smbref,z;*/
355 
356  /*Loop over all the elements of this partition*/
357  for(int i=0;i<femmodel->elements->Size();i++){
358  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
359 
360  /*Get reference SMB (uncorrected) and allocate all arrays*/
361  int numvertices = element->GetNumberOfVertices();
362  IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
363  IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
364  IssmDouble* smblist = xNew<IssmDouble>(numvertices);
365  element->GetInputListOnVertices(surfacelist,SurfaceEnum);
366  element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
367 
368  /*Loop over all vertices of element and correct SMB as a function of altitude z*/
369  for(int v=0;v<numvertices;v++){
370 
371  /*Get vertex elevation, anoma smb*/
372  z = surfacelist[v];
373  anomaly = smblistref[v];
374 
375  /* Henning edited acc. to Riannes equations*/
376  /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
377  z_critical = z_critical + dz;
378 
379  /* Calculate smb acc. to the surface elevation z */
380  if(z<z_critical){
381  smb = a + b*z + c;
382  }
383  else{
384  smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
385  }
386 
387  /* Compute smb including anomaly,
388  correct for number of seconds in a year [s/yr]*/
389  smb = smb/yts + anomaly;
390 
391  /*Update array accordingly*/
392  smblist[v] = smb;
393 
394  }
395 
396  /*Add input to element and Free memory*/
397  element->AddInput2(SmbMassBalanceEnum,smblist,P1Enum);
398  xDelete<IssmDouble>(surfacelist);
399  xDelete<IssmDouble>(smblistref);
400  xDelete<IssmDouble>(smblist);
401  }
402 
403 }/*}}}*/

◆ SmbComponentsx()

void SmbComponentsx ( FemModel *  femmodel)

Definition at line 404 of file SurfaceMassBalancex.cpp.

404  {/*{{{*/
405 
406  // void SmbComponentsx(acc,evap,runoff,ni){
407  // INPUT parameters: ni: working size of arrays
408  // INPUT: surface accumulation (m/yr water equivalent): acc
409  // surface evaporation (m/yr water equivalent): evap
410  // surface runoff (m/yr water equivalent): runoff
411  // OUTPUT: mass-balance (m/yr ice): agd(NA)
412  bool isclimatology;
413  femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
414 
415  /*Loop over all the elements of this partition*/
416  for(int i=0;i<femmodel->elements->Size();i++){
417  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
418 
419  /*Allocate all arrays*/
420  int numvertices = element->GetNumberOfVertices();
421  IssmDouble* acc = xNew<IssmDouble>(numvertices);
422  IssmDouble* evap = xNew<IssmDouble>(numvertices);
423  IssmDouble* runoff = xNew<IssmDouble>(numvertices);
424  IssmDouble* smb = xNew<IssmDouble>(numvertices);
425 
426  /*Recover Smb Components*/
427  if(isclimatology){
428 
429  int offsetend;
430  IssmDouble time0,timeend,timeclim;
431  IssmDouble time,starttime,finaltime;
432 
433  /*Get time parameters*/
434  femmodel->parameters->FindParam(&time,TimeEnum);
435  femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
436  femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
437 
438  //If this is a climatology, we need to repeat the forcing after the final time
439  TransientInput2* acc_input = element->inputs2->GetTransientInput(SmbAccumulationEnum); _assert_(acc_input);
440  TransientInput2* evap_input = element->inputs2->GetTransientInput(SmbEvaporationEnum); _assert_(evap_input);
441  TransientInput2* runoff_input = element->inputs2->GetTransientInput(SmbRunoffEnum); _assert_(runoff_input);
442 
443  //Get accumulation climatology value
444  offsetend = acc_input->GetTimeInputOffset(finaltime);
445  time0 = acc_input->GetTimeByOffset(-1);
446  timeend = acc_input->GetTimeByOffset(offsetend);
447  timeclim = time;
448  if(time>time0 & timeend>time0){
449  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
450  if(delta==0.)
451  timeclim=timeend;
452  else
453  timeclim=time0+delta;
454  }
455  element->GetInputListOnVerticesAtTime(acc,SmbAccumulationEnum,timeclim);
456 
457  //Get evaporation climatology value
458  offsetend = evap_input->GetTimeInputOffset(finaltime);
459  time0 = evap_input->GetTimeByOffset(-1);
460  timeend = evap_input->GetTimeByOffset(offsetend);
461  timeclim = time;
462  if(time>time0 & timeend>time0){
463  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
464  if(delta==0.)
465  timeclim=timeend;
466  else
467  timeclim=time0+delta;
468  }
469  element->GetInputListOnVerticesAtTime(evap,SmbEvaporationEnum,timeclim);
470 
471  //Get runoff climatology value
472  offsetend = runoff_input->GetTimeInputOffset(finaltime);
473  time0 = runoff_input->GetTimeByOffset(-1);
474  timeend = runoff_input->GetTimeByOffset(offsetend);
475  timeclim = time;
476  if(time>time0 & timeend>time0){
477  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
478  if(delta==0.)
479  timeclim=timeend;
480  else
481  timeclim=time0+delta;
482  }
483  element->GetInputListOnVerticesAtTime(runoff,SmbRunoffEnum,timeclim);
484  }
485  else{
486  element->GetInputListOnVertices(acc,SmbAccumulationEnum);
487  element->GetInputListOnVertices(evap,SmbEvaporationEnum);
488  element->GetInputListOnVertices(runoff,SmbRunoffEnum);
489  }
490 
491  // loop over all vertices
492  for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-runoff[v];
493 
494  /*Add input to element and Free memory*/
495  element->AddInput2(SmbMassBalanceEnum,smb,P1Enum);
496  xDelete<IssmDouble>(acc);
497  xDelete<IssmDouble>(evap);
498  xDelete<IssmDouble>(runoff);
499  xDelete<IssmDouble>(smb);
500  }
501 
502 }/*}}}*/

◆ SmbMeltComponentsx()

void SmbMeltComponentsx ( FemModel *  femmodel)

Definition at line 503 of file SurfaceMassBalancex.cpp.

503  {/*{{{*/
504 
505  // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
506  // INPUT parameters: ni: working size of arrays
507  // INPUT: surface accumulation (m/yr water equivalent): acc
508  // surface evaporation (m/yr water equivalent): evap
509  // surface melt (m/yr water equivalent): melt
510  // refreeze of surface melt (m/yr water equivalent): refreeze
511  // OUTPUT: mass-balance (m/yr ice): agd(NA)
512  bool isclimatology;
513  femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
514 
515  /*Loop over all the elements of this partition*/
516  for(int i=0;i<femmodel->elements->Size();i++){
517  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
518 
519  /*Allocate all arrays*/
520  int numvertices = element->GetNumberOfVertices();
521  IssmDouble* acc = xNew<IssmDouble>(numvertices);
522  IssmDouble* evap = xNew<IssmDouble>(numvertices);
523  IssmDouble* melt = xNew<IssmDouble>(numvertices);
524  IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
525  IssmDouble* smb = xNew<IssmDouble>(numvertices);
526 
527  /*Recover Smb Components*/
528  if (isclimatology){
529 
530  int offsetend;
531  IssmDouble time0,timeend,timeclim;
532  IssmDouble time,starttime,finaltime;
533  femmodel->parameters->FindParam(&time,TimeEnum);
534  femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
535  femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
536 
537 
538  //If this is a climatology, we need to repeat the forcing after the final time
539  TransientInput2* acc_input = element->inputs2->GetTransientInput(SmbAccumulationEnum); _assert_(acc_input);
540  TransientInput2* evap_input = element->inputs2->GetTransientInput(SmbEvaporationEnum); _assert_(evap_input);
541  TransientInput2* melt_input = element->inputs2->GetTransientInput(SmbMeltEnum); _assert_(melt_input);
542  TransientInput2* refreeze_input = element->inputs2->GetTransientInput(SmbRefreezeEnum); _assert_(refreeze_input);
543 
544  //Get accumulation climatology value
545  offsetend = acc_input->GetTimeInputOffset(finaltime);
546  time0 = acc_input->GetTimeByOffset(-1);
547  timeend = acc_input->GetTimeByOffset(offsetend);
548  timeclim = time;
549  if(time>time0 & timeend>time0){
550  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
551  if(delta==0.)
552  timeclim=timeend;
553  else
554  timeclim=time0+delta;
555  }
556  element->GetInputListOnVerticesAtTime(acc,SmbAccumulationEnum,timeclim);
557 
558  //Get evaporation climatology value
559  offsetend = evap_input->GetTimeInputOffset(finaltime);
560  time0 = evap_input->GetTimeByOffset(-1);
561  timeend = evap_input->GetTimeByOffset(offsetend);
562  timeclim = time;
563  if(time>time0 & timeend>time0){
564  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
565  if(delta==0.)
566  timeclim=timeend;
567  else
568  timeclim=time0+delta;
569  }
570  element->GetInputListOnVerticesAtTime(evap,SmbEvaporationEnum,timeclim);
571 
572  //Get melt climatology value
573  offsetend = melt_input->GetTimeInputOffset(finaltime);
574  time0 = melt_input->GetTimeByOffset(-1);
575  timeend = melt_input->GetTimeByOffset(offsetend);
576  timeclim = time;
577  if(time>time0 & timeend>time0){
578  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
579  if(delta==0.)
580  timeclim=timeend;
581  else
582  timeclim=time0+delta;
583  }
584  element->GetInputListOnVerticesAtTime(melt,SmbMeltEnum,timeclim);
585 
586  //Get refreeze climatology value
587  offsetend = refreeze_input->GetTimeInputOffset(finaltime);
588  time0 = refreeze_input->GetTimeByOffset(-1);
589  timeend = refreeze_input->GetTimeByOffset(offsetend);
590  timeclim = time;
591  if(time>time0 & timeend>time0){
592  IssmDouble delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
593  if(delta==0.)
594  timeclim=timeend;
595  else
596  timeclim=time0+delta;
597  }
598  element->GetInputListOnVerticesAtTime(refreeze,SmbRefreezeEnum,timeclim);
599  }
600  else{
601  element->GetInputListOnVertices(acc,SmbAccumulationEnum);
602  element->GetInputListOnVertices(evap,SmbEvaporationEnum);
603  element->GetInputListOnVertices(melt,SmbMeltEnum);
604  element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
605  }
606 
607  // loop over all vertices
608  for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
609 
610  /*Add input to element and Free memory*/
611  element->AddInput2(SmbMassBalanceEnum,smb,P1Enum);
612  xDelete<IssmDouble>(acc);
613  xDelete<IssmDouble>(evap);
614  xDelete<IssmDouble>(melt);
615  xDelete<IssmDouble>(refreeze);
616  xDelete<IssmDouble>(smb);
617  }
618 
619 }/*}}}*/

◆ SmbGradientsComponentsx()

void SmbGradientsComponentsx ( FemModel *  femmodel)

Definition at line 620 of file SurfaceMassBalancex.cpp.

620  {/*{{{*/
621 
622  for(int i=0;i<femmodel->elements->Size();i++){
623  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
624  element->SmbGradCompParameterization();
625  }
626 
627 }/*}}}*/

◆ SmbSemicx()

void SmbSemicx ( FemModel *  femmodel)

Definition at line 638 of file SurfaceMassBalancex.cpp.

638 {_error_("SEMIC not installed");}

◆ 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 **  a,
int  aIdx,
IssmDouble *  re,
IssmDouble *  dz,
IssmDouble *  d,
IssmDouble  cldFrac,
IssmDouble  aIce,
IssmDouble  aSnow,
IssmDouble  aValue,
IssmDouble  adThresh,
IssmDouble *  T,
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 } /*}}}*/

◆ 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 } /*}}}*/

◆ thermo()

void thermo ( IssmDouble *  pEC,
IssmDouble **  T,
IssmDouble *  pulwrf,
IssmDouble *  dz,
IssmDouble *  d,
IssmDouble *  swf,
IssmDouble  dlw,
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 } /*}}}*/

◆ 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  Ta,
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 } /*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
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: Gembx.cpp:1374
TimesteppingFinalTimeEnum
@ TimesteppingFinalTimeEnum
Definition: EnumDefinitions.h:430
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
Delflag
const double Delflag
Definition: Gembx.cpp:32
newcell
void newcell(IssmDouble **pcell, IssmDouble newvalue, bool top, int m)
Definition: MatrixUtils.cpp:577
SmbBPosEnum
@ SmbBPosEnum
Definition: EnumDefinitions.h:713
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
SmbDtEnum
@ SmbDtEnum
Definition: EnumDefinitions.h:357
IssmDouble
double IssmDouble
Definition: types.h:37
SmbAccumulationEnum
@ SmbAccumulationEnum
Definition: EnumDefinitions.h:708
xIsNan
int xIsNan(const T &X)
Definition: isnan.h:16
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
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
Element::SmbGradCompParameterization
void SmbGradCompParameterization(void)
Definition: Element.cpp:657
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
SmbSmbrefEnum
@ SmbSmbrefEnum
Definition: EnumDefinitions.h:780
TransientInput2
Definition: TransientInput2.h:13
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
SmbBMaxEnum
@ SmbBMaxEnum
Definition: EnumDefinitions.h:710
CtoK
const double CtoK
Definition: Gembx.cpp:12
Element::GetInputListOnVerticesAtTime
void GetInputListOnVerticesAtTime(IssmDouble *pvalue, int enumtype, IssmDouble time)
Definition: Element.cpp:1139
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
SmbTaEnum
@ SmbTaEnum
Definition: EnumDefinitions.h:782
CI
const double CI
Definition: Gembx.cpp:24
Element
Definition: Element.h:41
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Element::SmbGemb
void SmbGemb(IssmDouble timeinputs, int count)
Definition: Element.cpp:3584
Element::PositiveDegreeDay
void PositiveDegreeDay(IssmDouble *pdds, IssmDouble *pds, IssmDouble signorm, bool ismungsm, bool issetpddfac)
Definition: Element.cpp:2788
Element::inputs2
Inputs2 * inputs2
Definition: Element.h:47
TransientInput2::GetTimeByOffset
IssmDouble GetTimeByOffset(int offset)
Definition: TransientInput2.cpp:539
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
SmbIsfirnwarmingEnum
@ SmbIsfirnwarmingEnum
Definition: EnumDefinitions.h:368
Element::MungsmtpParameterization
void MungsmtpParameterization(void)
Definition: Element.cpp:2400
SmbEvaporationEnum
@ SmbEvaporationEnum
Definition: EnumDefinitions.h:738
VerboseSmb
bool VerboseSmb(void)
Definition: Verbosity.cpp:30
Element::PositiveDegreeDaySicopolis
void PositiveDegreeDaySicopolis(bool isfirnwarming)
Definition: Element.cpp:2975
SmbIssetpddfacEnum
@ SmbIssetpddfacEnum
Definition: EnumDefinitions.h:373
SmbIsclimatologyEnum
@ SmbIsclimatologyEnum
Definition: EnumDefinitions.h:363
Element::Delta18opdParameterization
void Delta18opdParameterization(void)
Definition: Element.cpp:531
FemModel::inputs2
Inputs2 * inputs2
Definition: FemModel.h:47
SmbRefreezeEnum
@ SmbRefreezeEnum
Definition: EnumDefinitions.h:770
InputDuplicatex
void InputDuplicatex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDuplicatex.cpp:10
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
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
TimesteppingStartTimeEnum
@ TimesteppingStartTimeEnum
Definition: EnumDefinitions.h:432
SmbElaEnum
@ SmbElaEnum
Definition: EnumDefinitions.h:737
SmbBMinEnum
@ SmbBMinEnum
Definition: EnumDefinitions.h:711
TimesteppingInterpForcingsEnum
@ TimesteppingInterpForcingsEnum
Definition: EnumDefinitions.h:431
Inputs2::GetTransientInput
TransientInput2 * GetTransientInput(int enum_type)
Definition: Inputs2.cpp:406
SmbRunoffEnum
@ SmbRunoffEnum
Definition: EnumDefinitions.h:772
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
SmbHrefEnum
@ SmbHrefEnum
Definition: EnumDefinitions.h:744
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
SmbBNegEnum
@ SmbBNegEnum
Definition: EnumDefinitions.h:712
SB
const double SB
Definition: Gembx.cpp:28
Inputs2::DeleteInput
int DeleteInput(int enum_type)
Definition: Inputs2.cpp:195
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
SmbIsmungsmEnum
@ SmbIsmungsmEnum
Definition: EnumDefinitions.h:371
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
SmbMassBalanceClimateEnum
@ SmbMassBalanceClimateEnum
Definition: EnumDefinitions.h:747
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
Marbouty
IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT)
Definition: Gembx.cpp:166
SmbMeltEnum
@ SmbMeltEnum
Definition: EnumDefinitions.h:754
Element::GetInputListOnVertices
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1131
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
Element::Delta18oParameterization
void Delta18oParameterization(void)
Definition: Element.cpp:432
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