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

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

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