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

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

NEW: randomflag allows to use fixed random seed in random.cpp, zero-initialization of smbspin for SmbAutoregression
File size: 20.3 KB
Line 
1/*!\file SurfaceMassBalancex
2 * \brief: calculates SMB
3 */
4
5#include <config.h>
6#include "./SurfaceMassBalancex.h"
7#include "../../shared/shared.h"
8#include "../../toolkits/toolkits.h"
9#include "../modules.h"
10#include "../../classes/Inputs/TransientInput.h"
11#include "../../shared/Random/random.h"
12
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}/*}}}*/
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)
26 int v;
27 IssmDouble rho_water; // density of fresh water
28 IssmDouble rho_ice; // density of ice
29 IssmDouble yts; // conversion factor year to second
30
31 /*Loop over all the elements of this partition*/
32 for(Object* & object : femmodel->elements->objects){
33 Element* element=xDynamicCast<Element*>(object);
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*/
45 element->GetInputListOnVertices(Href,SmbHrefEnum);
46 element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
47 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
48 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
49
50 /*Recover surface elevation at vertices: */
51 element->GetInputListOnVertices(s,SurfaceEnum);
52
53 /*Get material parameters :*/
54 rho_ice=element->FindParam(MaterialsRhoIceEnum);
55 rho_water=element->FindParam(MaterialsRhoFreshwaterEnum);
56
57 /* Get constants */
58 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
59
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 }
68
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*/
73 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
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);
80 }
81
82}/*}}}*/
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*/
92 for(Object* & object : femmodel->elements->objects){
93 Element* element=xDynamicCast<Element*>(object);
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
118 if(s[v]>ela[v]){
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*/
137 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
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}/*}}}*/
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 = xNew<IssmDouble>(numbasins);
157 IssmDouble* beta1 = xNew<IssmDouble>(numbasins);
158 IssmDouble* phi = xNew<IssmDouble>(numbasins*arorder);
159 IssmDouble* covmat = xNew<IssmDouble>(numbasins*numbasins);
160 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
161 femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
162 femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
163 femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);
164 femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);
165 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
166 femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
167
168 /*AR model spin-up*/
169 int nspin{2*arorder+5};
170 IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin);
171 my_rank=IssmComm::GetRank();
172 if(my_rank==0){
173 bool randomflag{};
174 femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
175 int fixedseed;
176 for(int i=0;i<nspin;i++){
177 IssmDouble* temparray = xNew<IssmDouble>(numbasins);
178 /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
179 if(randomflag) fixedseed=-1;
180 else fixedseed = i;
181 multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
182 for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
183 xDelete<IssmDouble>(temparray);
184 }
185 }
186 ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
187
188 /*Initialize SmbValuesAutoregressionEnum*/
189 for(Object* &object:femmodel->elements->objects){
190 Element* element = xDynamicCast<Element*>(object); //generate element object
191 element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
192 }
193
194 /*Cleanup*/
195 xDelete<IssmDouble>(beta0);
196 xDelete<IssmDouble>(beta1);
197 xDelete<IssmDouble>(phi);
198 xDelete<IssmDouble>(noisespin);
199 xDelete<IssmDouble>(covmat);
200}/*}}}*/
201void Smbautoregressionx(FemModel* femmodel){/*{{{*/
202
203 /*Get time parameters*/
204 IssmDouble time,dt,starttime,tstep_ar;
205 femmodel->parameters->FindParam(&time,TimeEnum);
206 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
207 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
208 femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
209
210 /*Initialize module at first time step*/
211 if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
212 /*Determine if this is a time step for the AR model*/
213 bool isstepforar = false;
214
215 #ifndef _HAVE_AD_
216 if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
217 #else
218 _error_("not implemented yet");
219 #endif
220
221 /*Load parameters*/
222 int M,N,Nphi,arorder,numbasins,my_rank;
223 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
224 femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
225 IssmDouble tinit_ar;
226 IssmDouble* beta0 = xNew<IssmDouble>(numbasins);
227 IssmDouble* beta1 = xNew<IssmDouble>(numbasins);
228 IssmDouble* phi = xNew<IssmDouble>(numbasins*arorder);
229 IssmDouble* noiseterms = xNew<IssmDouble>(numbasins);
230 IssmDouble* covmat = xNew<IssmDouble>(numbasins*numbasins);
231 femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
232 femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);
233 femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);
234 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
235 femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
236
237 /*Time elapsed with respect to AR model initial time*/
238 IssmDouble telapsed_ar = time-tinit_ar;
239
240 /*Before looping through elements: compute noise term specific to each basin from covmat*/
241 my_rank=IssmComm::GetRank();
242 if(my_rank==0){
243 bool randomflag{};
244 femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
245 _printf_("randomflag: "<<randomflag<<'\n');
246 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
247 int fixedseed;
248 if(randomflag) fixedseed=-1;
249 else fixedseed = static_cast<int>((time-starttime)/dt);
250 multivariateNormal(&noiseterms,numbasins,0.0,covmat,fixedseed);
251 }
252 ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
253
254 /*Loop over each element to compute SMB at vertices*/
255 for(Object* &object:femmodel->elements->objects){
256 Element* element = xDynamicCast<Element*>(object);
257 element->Smbautoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms);
258 }
259
260 /*Cleanup*/
261 xDelete<IssmDouble>(beta0);
262 xDelete<IssmDouble>(beta1);
263 xDelete<IssmDouble>(phi);
264 xDelete<IssmDouble>(noiseterms);
265 xDelete<IssmDouble>(covmat);
266}/*}}}*/
267void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
268
269 for(Object* & object : femmodel->elements->objects){
270 Element* element=xDynamicCast<Element*>(object);
271 element->Delta18oParameterization();
272 }
273
274}/*}}}*/
275void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
276
277 for(Object* & object : femmodel->elements->objects){
278 Element* element=xDynamicCast<Element*>(object);
279 element->MungsmtpParameterization();
280 }
281
282}/*}}}*/
283void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
284
285 for(Object* & object : femmodel->elements->objects){
286 Element* element=xDynamicCast<Element*>(object);
287 element->Delta18opdParameterization();
288 }
289
290}/*}}}*/
291void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
292
293 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
294 // note "v" prefix means 12 monthly means, ie time dimension
295 // INPUT parameters: ni: working size of arrays
296 // INPUT: surface elevation (m): hd(NA)
297 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
298 // ,NTIME)
299 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
300 // OUTPUT: mass-balance (m/yr ice): agd(NA)
301 // mean annual surface temperature (degrees C): Tsurf(NA)
302
303 int it, jj, itm;
304 IssmDouble DT = 0.02, sigfac, snormfac;
305 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
306 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
307 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
308 IssmDouble siglimc, siglim0, siglim0c;
309 IssmDouble tstep, tsint, tint, tstepc;
310 int NPDMAX = 1504, NPDCMAX = 1454;
311 //IssmDouble pdds[NPDMAX]={0};
312 //IssmDouble pds[NPDCMAX]={0};
313 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
314 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
315 IssmDouble tstar; // monthly mean surface temp
316
317 bool ismungsm;
318 bool issetpddfac;
319
320 IssmDouble *pdds = NULL;
321 IssmDouble *pds = NULL;
322 Element *element = NULL;
323
324 pdds=xNew<IssmDouble>(NPDMAX+1);
325 pds=xNew<IssmDouble>(NPDCMAX+1);
326
327 // Get ismungsm parameter
328 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
329
330 // Get issetpddfac parameter
331 femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
332
333 /* initialize PDD (creation of a lookup table)*/
334 tstep = 0.1;
335 tsint = tstep*0.5;
336 sigfac = -1.0/(2.0*pow(signorm,2));
337 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
338 siglim = 2.5*signorm;
339 siglimc = 2.5*signormc;
340 siglim0 = siglim/DT + 0.5;
341 siglim0c = siglimc/DT + 0.5;
342 PDup = siglimc+PDCUT;
343
344 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
345
346 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
347 for(it = 0; it < itm; it++){
348 // tstar = REAL(it)*DT-siglim;
349 tstar = it*DT-siglim;
350 tint = tsint;
351 pddt = 0.;
352 for ( jj = 0; jj < 600; jj++){
353 if (tint > (tstar+siglim)){break;}
354 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
355 tint = tint+tstep;
356 }
357 pdds[it] = pddt*snormfac;
358 }
359 pdds[itm+1] = siglim + DT;
360
361 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
362 tstepc = 0.1;
363 tsint = PDCUT-tstepc*0.5;
364 signormc = signorm - 0.5;
365 sigfac = -1.0/(2.0*pow(signormc,2));
366 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
367 siglimc = 2.5*signormc ;
368 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
369 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
370 for(it = 0; it < itm; it++ ){
371 tstar = it*DT-siglimc;
372 // tstar = REAL(it)*DT-siglimc;
373 tint = tsint; // start against upper bound
374 pd = 0.;
375 for (jj = 0; jj < 600; jj++){
376 if (tint<(tstar-siglimc)) {break;}
377 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
378 tint = tint-tstepc;
379 }
380 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
381 }
382 pds[itm+1] = 0.;
383 // *******END initialize PDD
384
385 for(Object* & object : femmodel->elements->objects){
386 element=xDynamicCast<Element*>(object);
387 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
388 }
389 /*free ressouces: */
390 xDelete<IssmDouble>(pdds);
391 xDelete<IssmDouble>(pds);
392}/*}}}*/
393void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
394
395 bool isfirnwarming;
396 femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
397
398 for(Object* & object : femmodel->elements->objects){
399 Element* element=xDynamicCast<Element*>(object);
400 element->PositiveDegreeDaySicopolis(isfirnwarming);
401 }
402
403}/*}}}*/
404void SmbHenningx(FemModel* femmodel){/*{{{*/
405
406 /*Intermediaries*/
407 IssmDouble z_critical = 1675.;
408 IssmDouble dz = 0;
409 IssmDouble a = -15.86;
410 IssmDouble b = 0.00969;
411 IssmDouble c = -0.235;
412 IssmDouble f = 1.;
413 IssmDouble g = -0.0011;
414 IssmDouble h = -1.54e-5;
415 IssmDouble smb,smbref,anomaly,yts,z;
416
417 /* Get constants */
418 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
419 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
420 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
421 /*Mathieu original*/
422 /*IssmDouble smb,smbref,z;*/
423
424 /*Loop over all the elements of this partition*/
425 for(Object* & object : femmodel->elements->objects){
426 Element* element=xDynamicCast<Element*>(object);
427
428 /*Get reference SMB (uncorrected) and allocate all arrays*/
429 int numvertices = element->GetNumberOfVertices();
430 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
431 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
432 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
433 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
434 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
435
436 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
437 for(int v=0;v<numvertices;v++){
438
439 /*Get vertex elevation, anoma smb*/
440 z = surfacelist[v];
441 anomaly = smblistref[v];
442
443 /* Henning edited acc. to Riannes equations*/
444 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
445 z_critical = z_critical + dz;
446
447 /* Calculate smb acc. to the surface elevation z */
448 if(z<z_critical){
449 smb = a + b*z + c;
450 }
451 else{
452 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
453 }
454
455 /* Compute smb including anomaly,
456 correct for number of seconds in a year [s/yr]*/
457 smb = smb/yts + anomaly;
458
459 /*Update array accordingly*/
460 smblist[v] = smb;
461
462 }
463
464 /*Add input to element and Free memory*/
465 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
466 xDelete<IssmDouble>(surfacelist);
467 xDelete<IssmDouble>(smblistref);
468 xDelete<IssmDouble>(smblist);
469 }
470
471}/*}}}*/
472void SmbComponentsx(FemModel* femmodel){/*{{{*/
473
474 // void SmbComponentsx(acc,evap,runoff,ni){
475 // INPUT parameters: ni: working size of arrays
476 // INPUT: surface accumulation (m/yr water equivalent): acc
477 // surface evaporation (m/yr water equivalent): evap
478 // surface runoff (m/yr water equivalent): runoff
479 // OUTPUT: mass-balance (m/yr ice): agd(NA)
480
481 /*Loop over all the elements of this partition*/
482 for(Object* & object : femmodel->elements->objects){
483 Element* element=xDynamicCast<Element*>(object);
484
485 /*Allocate all arrays*/
486 int numvertices = element->GetNumberOfVertices();
487 IssmDouble* acc = xNew<IssmDouble>(numvertices);
488 IssmDouble* evap = xNew<IssmDouble>(numvertices);
489 IssmDouble* runoff = 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(runoff,SmbRunoffEnum);
496
497 // loop over all vertices
498 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-runoff[v];
499
500 /*Add input to element and Free memory*/
501 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
502 xDelete<IssmDouble>(acc);
503 xDelete<IssmDouble>(evap);
504 xDelete<IssmDouble>(runoff);
505 xDelete<IssmDouble>(smb);
506 }
507
508}/*}}}*/
509void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
510
511 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
512 // INPUT parameters: ni: working size of arrays
513 // INPUT: surface accumulation (m/yr water equivalent): acc
514 // surface evaporation (m/yr water equivalent): evap
515 // surface melt (m/yr water equivalent): melt
516 // refreeze of surface melt (m/yr water equivalent): refreeze
517 // OUTPUT: mass-balance (m/yr ice): agd(NA)
518
519 /*Loop over all the elements of this partition*/
520 for(Object* & object : femmodel->elements->objects){
521 Element* element=xDynamicCast<Element*>(object);
522
523 /*Allocate all arrays*/
524 int numvertices = element->GetNumberOfVertices();
525 IssmDouble* acc = xNew<IssmDouble>(numvertices);
526 IssmDouble* evap = xNew<IssmDouble>(numvertices);
527 IssmDouble* melt = xNew<IssmDouble>(numvertices);
528 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
529 IssmDouble* smb = xNew<IssmDouble>(numvertices);
530
531 /*Recover Smb Components*/
532 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
533 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
534 element->GetInputListOnVertices(melt,SmbMeltEnum);
535 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
536
537 // loop over all vertices
538 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
539
540 /*Add input to element and Free memory*/
541 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
542 xDelete<IssmDouble>(acc);
543 xDelete<IssmDouble>(evap);
544 xDelete<IssmDouble>(melt);
545 xDelete<IssmDouble>(refreeze);
546 xDelete<IssmDouble>(smb);
547 }
548
549}/*}}}*/
550void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
551
552 for(Object* & object : femmodel->elements->objects){
553 Element* element=xDynamicCast<Element*>(object);
554 element->SmbGradCompParameterization();
555 }
556
557}/*}}}*/
558#ifdef _HAVE_SEMIC_
559void SmbSemicx(FemModel* femmodel){/*{{{*/
560
561 for(Object* & object : femmodel->elements->objects){
562 Element* element=xDynamicCast<Element*>(object);
563 element->SmbSemic();
564 }
565
566}/*}}}*/
567#else
568void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");}
569#endif //_HAVE_SEMIC_
Note: See TracBrowser for help on using the repository browser.