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

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

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