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

Last change on this file since 26477 was 26477, checked in by vverjans, 3 years ago

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