source: issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp@ 23644

Last change on this file since 23644 was 23644, checked in by Mathieu Morlighem, 6 years ago

NEW: removing matpar from femmodel
File size: 15.4 KB
Line 
1/*!\file SurfaceMassBalancex
2 * \brief: calculates SMB
3 */
4
5#include "./SurfaceMassBalancex.h"
6#include "../../shared/shared.h"
7#include "../../toolkits/toolkits.h"
8
9void SmbGradientsx(FemModel* femmodel){/*{{{*/
10
11 // void SurfaceMassBalancex(hd,agd,ni){
12 // INPUT parameters: ni: working size of arrays
13 // INPUT: surface elevation (m): hd(NA)
14 // OUTPUT: mass-balance (m/yr ice): agd(NA)
15 int v;
16 IssmDouble rho_water; // density of fresh water
17 IssmDouble rho_ice; // density of ice
18 IssmDouble yts; // conversion factor year to second
19
20 /*Loop over all the elements of this partition*/
21 for(int i=0;i<femmodel->elements->Size();i++){
22 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
23
24 /*Allocate all arrays*/
25 int numvertices = element->GetNumberOfVertices();
26 IssmDouble* Href = xNew<IssmDouble>(numvertices); // reference elevation from which deviations are used to calculate the SMB adjustment
27 IssmDouble* Smbref = xNew<IssmDouble>(numvertices); // reference SMB to which deviations are added
28 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // Hs-SMB relation parameter
29 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // Hs-SMB relation paremeter
30 IssmDouble* s = xNew<IssmDouble>(numvertices); // surface elevation (m)
31 IssmDouble* smb = xNew<IssmDouble>(numvertices);
32
33 /*Recover SmbGradients*/
34 element->GetInputListOnVertices(Href,SmbHrefEnum);
35 element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
36 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
37 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
38
39 /*Recover surface elevation at vertices: */
40 element->GetInputListOnVertices(s,SurfaceEnum);
41
42 /*Get material parameters :*/
43 rho_ice=element->FindParam(MaterialsRhoIceEnum);
44 rho_water=element->FindParam(MaterialsRhoFreshwaterEnum);
45
46 /* Get constants */
47 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
48
49 // loop over all vertices
50 for(v=0;v<numvertices;v++){
51 if(Smbref[v]>0){
52 smb[v]=Smbref[v]+b_pos[v]*(s[v]-Href[v]);
53 }
54 else{
55 smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
56 }
57
58 smb[v]=smb[v]/1000*rho_water/rho_ice; // SMB in m/y ice
59 } //end of the loop over the vertices
60
61 /*Add input to element and Free memory*/
62 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
63 xDelete<IssmDouble>(Href);
64 xDelete<IssmDouble>(Smbref);
65 xDelete<IssmDouble>(b_pos);
66 xDelete<IssmDouble>(b_neg);
67 xDelete<IssmDouble>(s);
68 xDelete<IssmDouble>(smb);
69 }
70
71}/*}}}*/
72void SmbGradientsElax(FemModel* femmodel){/*{{{*/
73
74 // void SurfaceMassBalancex(hd,agd,ni){
75 // INPUT parameters: ni: working size of arrays
76 // INPUT: surface elevation (m): hd(NA)
77 // OUTPUT: surface mass-balance (m/yr ice): agd(NA)
78 int v;
79
80 /*Loop over all the elements of this partition*/
81 for(int i=0;i<femmodel->elements->Size();i++){
82 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
83
84 /*Allocate all arrays*/
85 int numvertices = element->GetNumberOfVertices();
86 IssmDouble* ela = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
87 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
88 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
89 IssmDouble* b_max = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
90 IssmDouble* b_min = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
91 IssmDouble* s = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
92 IssmDouble* smb = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
93
94 /*Recover ELA, SMB gradients, and caps*/
95 element->GetInputListOnVertices(ela,SmbElaEnum);
96 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
97 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
98 element->GetInputListOnVertices(b_max,SmbBMaxEnum);
99 element->GetInputListOnVertices(b_min,SmbBMinEnum);
100
101 /*Recover surface elevation at vertices: */
102 element->GetInputListOnVertices(s,SurfaceEnum);
103
104 /*Loop over all vertices, calculate SMB*/
105 for(v=0;v<numvertices;v++){
106 // if surface is above the ELA
107 if(s[v]>ela[v]){
108 smb[v]=b_pos[v]*(s[v]-ela[v]);
109 }
110 // if surface is below or equal to the ELA
111 else{
112 smb[v]=b_neg[v]*(s[v]-ela[v]);
113 }
114
115 // if SMB is larger than upper cap, set SMB to upper cap
116 if(smb[v]>b_max[v]){
117 smb[v]=b_max[v];
118 }
119 // if SMB is smaller than lower cap, set SMB to lower cap
120 if(smb[v]<b_min[v]){
121 smb[v]=b_min[v];
122 }
123 } //end of the loop over the vertices
124
125 /*Add input to element and Free memory*/
126 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
127 xDelete<IssmDouble>(ela);
128 xDelete<IssmDouble>(b_pos);
129 xDelete<IssmDouble>(b_neg);
130 xDelete<IssmDouble>(b_max);
131 xDelete<IssmDouble>(b_min);
132 xDelete<IssmDouble>(s);
133 xDelete<IssmDouble>(smb);
134
135 }
136
137}/*}}}*/
138void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
139
140 for(int i=0;i<femmodel->elements->Size();i++){
141 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
142 element->Delta18oParameterization();
143 }
144
145}/*}}}*/
146void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
147
148 for(int i=0;i<femmodel->elements->Size();i++){
149 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
150 element->MungsmtpParameterization();
151 }
152
153}/*}}}*/
154void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
155
156 for(int i=0;i<femmodel->elements->Size();i++){
157 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
158 element->Delta18opdParameterization();
159 }
160
161}/*}}}*/
162void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
163
164 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
165 // note "v" prefix means 12 monthly means, ie time dimension
166 // INPUT parameters: ni: working size of arrays
167 // INPUT: surface elevation (m): hd(NA)
168 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
169 // ,NTIME)
170 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
171 // OUTPUT: mass-balance (m/yr ice): agd(NA)
172 // mean annual surface temperature (degrees C): Tsurf(NA)
173
174 int i, it, jj, itm;
175 IssmDouble DT = 0.02, sigfac, snormfac;
176 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
177 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
178 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
179 IssmDouble siglimc, siglim0, siglim0c;
180 IssmDouble tstep, tsint, tint, tstepc;
181 int NPDMAX = 1504, NPDCMAX = 1454;
182 //IssmDouble pdds[NPDMAX]={0};
183 //IssmDouble pds[NPDCMAX]={0};
184 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
185 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
186 IssmDouble tstar; // monthly mean surface temp
187
188 bool ismungsm;
189 bool issetpddfac;
190
191 IssmDouble *pdds = NULL;
192 IssmDouble *pds = NULL;
193 Element *element = NULL;
194
195 pdds=xNew<IssmDouble>(NPDMAX+1);
196 pds=xNew<IssmDouble>(NPDCMAX+1);
197
198 // Get ismungsm parameter
199 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
200
201 // Get issetpddfac parameter
202 femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
203
204 /* initialize PDD (creation of a lookup table)*/
205 tstep = 0.1;
206 tsint = tstep*0.5;
207 sigfac = -1.0/(2.0*pow(signorm,2));
208 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
209 siglim = 2.5*signorm;
210 siglimc = 2.5*signormc;
211 siglim0 = siglim/DT + 0.5;
212 siglim0c = siglimc/DT + 0.5;
213 PDup = siglimc+PDCUT;
214
215 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
216
217 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
218 for(it = 0; it < itm; it++){
219 // tstar = REAL(it)*DT-siglim;
220 tstar = it*DT-siglim;
221 tint = tsint;
222 pddt = 0.;
223 for ( jj = 0; jj < 600; jj++){
224 if (tint > (tstar+siglim)){break;}
225 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
226 tint = tint+tstep;
227 }
228 pdds[it] = pddt*snormfac;
229 }
230 pdds[itm+1] = siglim + DT;
231
232 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
233 tstepc = 0.1;
234 tsint = PDCUT-tstepc*0.5;
235 signormc = signorm - 0.5;
236 sigfac = -1.0/(2.0*pow(signormc,2));
237 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
238 siglimc = 2.5*signormc ;
239 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
240 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
241 for(it = 0; it < itm; it++ ){
242 tstar = it*DT-siglimc;
243 // tstar = REAL(it)*DT-siglimc;
244 tint = tsint; // start against upper bound
245 pd = 0.;
246 for (jj = 0; jj < 600; jj++){
247 if (tint<(tstar-siglimc)) {break;}
248 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
249 tint = tint-tstepc;
250 }
251 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
252 }
253 pds[itm+1] = 0.;
254 // *******END initialize PDD
255
256 for(i=0;i<femmodel->elements->Size();i++){
257 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
258 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
259 }
260 /*free ressouces: */
261 xDelete<IssmDouble>(pdds);
262 xDelete<IssmDouble>(pds);
263}/*}}}*/
264void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
265
266 bool isfirnwarming;
267 femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
268
269 for(int i=0;i<femmodel->elements->Size();i++){
270 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
271 element->PositiveDegreeDaySicopolis(isfirnwarming);
272 }
273
274}/*}}}*/
275void SmbHenningx(FemModel* femmodel){/*{{{*/
276
277 /*Intermediaries*/
278 IssmDouble z_critical = 1675.;
279 IssmDouble dz = 0;
280 IssmDouble a = -15.86;
281 IssmDouble b = 0.00969;
282 IssmDouble c = -0.235;
283 IssmDouble f = 1.;
284 IssmDouble g = -0.0011;
285 IssmDouble h = -1.54e-5;
286 IssmDouble smb,smbref,anomaly,yts,z;
287
288 /* Get constants */
289 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
290 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
291 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
292 /*Mathieu original*/
293 /*IssmDouble smb,smbref,z;*/
294
295 /*Loop over all the elements of this partition*/
296 for(int i=0;i<femmodel->elements->Size();i++){
297 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
298
299 /*Get reference SMB (uncorrected) and allocate all arrays*/
300 int numvertices = element->GetNumberOfVertices();
301 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
302 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
303 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
304 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
305 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
306
307 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
308 for(int v=0;v<numvertices;v++){
309
310 /*Get vertex elevation, anoma smb*/
311 z = surfacelist[v];
312 anomaly = smblistref[v];
313
314 /* Henning edited acc. to Riannes equations*/
315 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
316 z_critical = z_critical + dz;
317
318 /* Calculate smb acc. to the surface elevation z */
319 if(z<z_critical){
320 smb = a + b*z + c;
321 }
322 else{
323 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
324 }
325
326 /* Compute smb including anomaly,
327 correct for number of seconds in a year [s/yr]*/
328 smb = smb/yts + anomaly;
329
330 /*Update array accordingly*/
331 smblist[v] = smb;
332
333 }
334
335 /*Add input to element and Free memory*/
336 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
337 xDelete<IssmDouble>(surfacelist);
338 xDelete<IssmDouble>(smblistref);
339 xDelete<IssmDouble>(smblist);
340 }
341
342}/*}}}*/
343void SmbComponentsx(FemModel* femmodel){/*{{{*/
344
345 // void SmbComponentsx(acc,evap,runoff,ni){
346 // INPUT parameters: ni: working size of arrays
347 // INPUT: surface accumulation (m/yr water equivalent): acc
348 // surface evaporation (m/yr water equivalent): evap
349 // surface runoff (m/yr water equivalent): runoff
350 // OUTPUT: mass-balance (m/yr ice): agd(NA)
351 int v;
352
353 /*Loop over all the elements of this partition*/
354 for(int i=0;i<femmodel->elements->Size();i++){
355 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
356
357 /*Allocate all arrays*/
358 int numvertices = element->GetNumberOfVertices();
359 IssmDouble* acc = xNew<IssmDouble>(numvertices);
360 IssmDouble* evap = xNew<IssmDouble>(numvertices);
361 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
362 IssmDouble* smb = xNew<IssmDouble>(numvertices);
363
364 /*Recover Smb Components*/
365 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
366 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
367 element->GetInputListOnVertices(runoff,SmbRunoffEnum);
368
369 // loop over all vertices
370 for(v=0;v<numvertices;v++){
371 smb[v]=acc[v]-evap[v]-runoff[v];
372 } //end of the loop over the vertices
373
374 /*Add input to element and Free memory*/
375 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
376 xDelete<IssmDouble>(acc);
377 xDelete<IssmDouble>(evap);
378 xDelete<IssmDouble>(runoff);
379 xDelete<IssmDouble>(smb);
380 }
381
382}/*}}}*/
383void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
384
385 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
386 // INPUT parameters: ni: working size of arrays
387 // INPUT: surface accumulation (m/yr water equivalent): acc
388 // surface evaporation (m/yr water equivalent): evap
389 // surface melt (m/yr water equivalent): melt
390 // refreeze of surface melt (m/yr water equivalent): refreeze
391 // OUTPUT: mass-balance (m/yr ice): agd(NA)
392 int v;
393
394 /*Loop over all the elements of this partition*/
395 for(int i=0;i<femmodel->elements->Size();i++){
396 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
397
398 /*Allocate all arrays*/
399 int numvertices = element->GetNumberOfVertices();
400 IssmDouble* acc = xNew<IssmDouble>(numvertices);
401 IssmDouble* evap = xNew<IssmDouble>(numvertices);
402 IssmDouble* melt = xNew<IssmDouble>(numvertices);
403 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
404 IssmDouble* smb = xNew<IssmDouble>(numvertices);
405
406 /*Recover Smb Components*/
407 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
408 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
409 element->GetInputListOnVertices(melt,SmbMeltEnum);
410 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
411
412 // loop over all vertices
413 for(v=0;v<numvertices;v++){
414 smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
415 } //end of the loop over the vertices
416
417 /*Add input to element and Free memory*/
418 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
419 xDelete<IssmDouble>(acc);
420 xDelete<IssmDouble>(evap);
421 xDelete<IssmDouble>(melt);
422 xDelete<IssmDouble>(refreeze);
423 xDelete<IssmDouble>(smb);
424 }
425
426}/*}}}*/
427void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
428
429 // void SurfaceMassBalancex(hd,agd,ni){
430 // INPUT parameters: ni: working size of arrays
431 // INPUT: surface elevation (m): hd(NA)
432 // OUTPUT: mass-balance (m/yr ice): agd(NA)
433
434 for(int i=0;i<femmodel->elements->Size();i++){
435 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
436 element->SmbGradCompParameterization();
437 }
438
439}/*}}}*/
440#ifdef _HAVE_SEMIC_
441void SmbSemicx(FemModel* femmodel){/*{{{*/
442
443 for(int i=0;i<femmodel->elements->Size();i++){
444 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
445 element->SmbSemic();
446 }
447
448}/*}}}*/
449#else
450void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");}
451#endif //_HAVE_SEMIC_
Note: See TracBrowser for help on using the repository browser.