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

Last change on this file since 27856 was 27856, checked in by rueckamp, 20 months ago

CHG: Replace SMBdebrisML with SMBdebrisEvatt; adding new options for SMBdebrisEvatt
File size: 19.1 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}/*}}}*/
[27250]149void Smbarmax(FemModel* femmodel){/*{{{*/
[26483]150
[26615]151 /*Get time parameters*/
[27250]152 IssmDouble time,dt,starttime,tstep_arma;
[26615]153 femmodel->parameters->FindParam(&time,TimeEnum);
154 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
155 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
[27250]156 femmodel->parameters->FindParam(&tstep_arma,SmbARMATimestepEnum);
[26482]157
[27250]158 /*Determine if this is a time step for the ARMA model*/
159 bool isstepforarma = false;
[26479]160
[26615]161 #ifndef _HAVE_AD_
[27250]162 if((fmod(time,tstep_arma)<fmod((time-dt),tstep_arma)) || (time<=starttime+dt) || tstep_arma==dt) isstepforarma = true;
[26615]163 #else
164 _error_("not implemented yet");
165 #endif
[26479]166
[26615]167 /*Load parameters*/
168 bool isstochastic;
169 bool issmbstochastic = false;
[27318]170 int M,N,arorder,maorder,numbasins,numparams,numbreaks,numelevbins,my_rank;
[26615]171 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
[27318]172 femmodel->parameters->FindParam(&numparams,SmbNumParamsEnum);
173 femmodel->parameters->FindParam(&numbreaks,SmbNumBreaksEnum);
174 femmodel->parameters->FindParam(&arorder,SmbARMAarOrderEnum);
[27260]175 femmodel->parameters->FindParam(&maorder,SmbARMAmaOrderEnum);
[26947]176 femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum);
[27318]177 IssmDouble* datebreaks = NULL;
178 IssmDouble* arlagcoefs = NULL;
[27260]179 IssmDouble* malagcoefs = NULL;
[27318]180 IssmDouble* polyparams = NULL;
[26947]181 IssmDouble* lapserates = NULL;
182 IssmDouble* elevbins = NULL;
183 IssmDouble* refelevation = NULL;
[26526]184
[27320]185 femmodel->parameters->FindParam(&datebreaks,&M,&N,SmbARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1));
[27319]186 femmodel->parameters->FindParam(&polyparams,&M,&N,SmbARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams);
[27318]187 femmodel->parameters->FindParam(&arlagcoefs,&M,&N,SmbARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder);
188 femmodel->parameters->FindParam(&malagcoefs,&M,&N,SmbARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder);
[27466]189 femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins*12);
190 femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==(numelevbins-1)*12);
[27250]191 femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins);
[26481]192
[26615]193 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
[26526]194 if(isstochastic){
[26615]195 int numstochasticfields;
[26526]196 int* stochasticfields;
197 femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
198 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
199 for(int i=0;i<numstochasticfields;i++){
[27250]200 if(stochasticfields[i]==SMBarmaEnum) issmbstochastic = true;
[26615]201 }
202 xDelete<int>(stochasticfields);
203 }
[26481]204
[26615]205 /*Loop over each element to compute SMB at vertices*/
206 for(Object* &object:femmodel->elements->objects){
207 Element* element = xDynamicCast<Element*>(object);
[27250]208 /*Compute ARMA*/
[27318]209 element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,issmbstochastic,SMBarmaEnum);
[26810]210 /*Compute lapse rate adjustment*/
[26947]211 element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation);
[26810]212 }
[26481]213
[26615]214 /*Cleanup*/
[27250]215 xDelete<IssmDouble>(arlagcoefs);
[27260]216 xDelete<IssmDouble>(malagcoefs);
[27318]217 xDelete<IssmDouble>(polyparams);
218 xDelete<IssmDouble>(datebreaks);
[26947]219 xDelete<IssmDouble>(lapserates);
220 xDelete<IssmDouble>(elevbins);
[26810]221 xDelete<IssmDouble>(refelevation);
[26477]222}/*}}}*/
[17085]223void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
224
[25539]225 for(Object* & object : femmodel->elements->objects){
226 Element* element=xDynamicCast<Element*>(object);
[17085]227 element->Delta18oParameterization();
228 }
229
230}/*}}}*/
[18968]231void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
232
[25539]233 for(Object* & object : femmodel->elements->objects){
234 Element* element=xDynamicCast<Element*>(object);
[18968]235 element->MungsmtpParameterization();
236 }
237
238}/*}}}*/
[19172]239void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
240
[25539]241 for(Object* & object : femmodel->elements->objects){
242 Element* element=xDynamicCast<Element*>(object);
[19172]243 element->Delta18opdParameterization();
244 }
245
246}/*}}}*/
[17085]247void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
248
249 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
250 // note "v" prefix means 12 monthly means, ie time dimension
251 // INPUT parameters: ni: working size of arrays
252 // INPUT: surface elevation (m): hd(NA)
253 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
[23366]254 // ,NTIME)
[17085]255 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
256 // OUTPUT: mass-balance (m/yr ice): agd(NA)
257 // mean annual surface temperature (degrees C): Tsurf(NA)
258
[25539]259 int it, jj, itm;
[17085]260 IssmDouble DT = 0.02, sigfac, snormfac;
[23366]261 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
[17085]262 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
263 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
264 IssmDouble siglimc, siglim0, siglim0c;
265 IssmDouble tstep, tsint, tint, tstepc;
266 int NPDMAX = 1504, NPDCMAX = 1454;
[23366]267 //IssmDouble pdds[NPDMAX]={0};
[17085]268 //IssmDouble pds[NPDCMAX]={0};
269 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
270 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
271 IssmDouble tstar; // monthly mean surface temp
272
[18968]273 bool ismungsm;
[22448]274 bool issetpddfac;
[18968]275
[17085]276 IssmDouble *pdds = NULL;
277 IssmDouble *pds = NULL;
278 Element *element = NULL;
279
[23366]280 pdds=xNew<IssmDouble>(NPDMAX+1);
281 pds=xNew<IssmDouble>(NPDCMAX+1);
[17085]282
[18968]283 // Get ismungsm parameter
[19527]284 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
[18968]285
[22448]286 // Get issetpddfac parameter
287 femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
288
[17085]289 /* initialize PDD (creation of a lookup table)*/
290 tstep = 0.1;
291 tsint = tstep*0.5;
292 sigfac = -1.0/(2.0*pow(signorm,2));
293 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
294 siglim = 2.5*signorm;
295 siglimc = 2.5*signormc;
296 siglim0 = siglim/DT + 0.5;
297 siglim0c = siglimc/DT + 0.5;
298 PDup = siglimc+PDCUT;
299
300 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
301
302 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
[23366]303 for(it = 0; it < itm; it++){
[17085]304 // tstar = REAL(it)*DT-siglim;
305 tstar = it*DT-siglim;
306 tint = tsint;
307 pddt = 0.;
308 for ( jj = 0; jj < 600; jj++){
309 if (tint > (tstar+siglim)){break;}
310 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
311 tint = tint+tstep;
312 }
313 pdds[it] = pddt*snormfac;
314 }
315 pdds[itm+1] = siglim + DT;
316
317 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
318 tstepc = 0.1;
319 tsint = PDCUT-tstepc*0.5;
320 signormc = signorm - 0.5;
321 sigfac = -1.0/(2.0*pow(signormc,2));
322 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
323 siglimc = 2.5*signormc ;
324 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
325 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
326 for(it = 0; it < itm; it++ ){
327 tstar = it*DT-siglimc;
328 // tstar = REAL(it)*DT-siglimc;
329 tint = tsint; // start against upper bound
330 pd = 0.;
331 for (jj = 0; jj < 600; jj++){
332 if (tint<(tstar-siglimc)) {break;}
333 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
334 tint = tint-tstepc;
335 }
336 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
337 }
338 pds[itm+1] = 0.;
339 // *******END initialize PDD
340
[25539]341 for(Object* & object : femmodel->elements->objects){
342 element=xDynamicCast<Element*>(object);
[22448]343 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
[17085]344 }
345 /*free ressouces: */
346 xDelete<IssmDouble>(pdds);
347 xDelete<IssmDouble>(pds);
348}/*}}}*/
[23317]349void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
[23366]350
[23317]351 bool isfirnwarming;
[23328]352 femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
[23366]353
[25539]354 for(Object* & object : femmodel->elements->objects){
355 Element* element=xDynamicCast<Element*>(object);
[23317]356 element->PositiveDegreeDaySicopolis(isfirnwarming);
357 }
358
359}/*}}}*/
[17085]360void SmbHenningx(FemModel* femmodel){/*{{{*/
361
[17087]362 /*Intermediaries*/
[17403]363 IssmDouble z_critical = 1675.;
364 IssmDouble dz = 0;
365 IssmDouble a = -15.86;
366 IssmDouble b = 0.00969;
367 IssmDouble c = -0.235;
368 IssmDouble f = 1.;
369 IssmDouble g = -0.0011;
370 IssmDouble h = -1.54e-5;
371 IssmDouble smb,smbref,anomaly,yts,z;
[22249]372
373 /* Get constants */
374 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
375 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
376 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
377 /*Mathieu original*/
378 /*IssmDouble smb,smbref,z;*/
379
[17087]380 /*Loop over all the elements of this partition*/
[25539]381 for(Object* & object : femmodel->elements->objects){
382 Element* element=xDynamicCast<Element*>(object);
[17087]383
384 /*Get reference SMB (uncorrected) and allocate all arrays*/
385 int numvertices = element->GetNumberOfVertices();
386 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
387 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
388 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
389 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
[19527]390 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
[17087]391
392 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
393 for(int v=0;v<numvertices;v++){
394
[17403]395 /*Get vertex elevation, anoma smb*/
[17087]396 z = surfacelist[v];
[17403]397 anomaly = smblistref[v];
[17087]398
[22249]399 /* Henning edited acc. to Riannes equations*/
400 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
401 z_critical = z_critical + dz;
402
403 /* Calculate smb acc. to the surface elevation z */
404 if(z<z_critical){
[17403]405 smb = a + b*z + c;
[17087]406 }
407 else{
[22249]408 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
[17087]409 }
[22249]410
[18584]411 /* Compute smb including anomaly,
412 correct for number of seconds in a year [s/yr]*/
413 smb = smb/yts + anomaly;
414
[17087]415 /*Update array accordingly*/
416 smblist[v] = smb;
417
418 }
419
420 /*Add input to element and Free memory*/
[25379]421 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
[17087]422 xDelete<IssmDouble>(surfacelist);
423 xDelete<IssmDouble>(smblistref);
424 xDelete<IssmDouble>(smblist);
[17085]425 }
426
427}/*}}}*/
[18001]428void SmbComponentsx(FemModel* femmodel){/*{{{*/
429
430 // void SmbComponentsx(acc,evap,runoff,ni){
431 // INPUT parameters: ni: working size of arrays
432 // INPUT: surface accumulation (m/yr water equivalent): acc
433 // surface evaporation (m/yr water equivalent): evap
434 // surface runoff (m/yr water equivalent): runoff
435 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[23814]436
[18001]437 /*Loop over all the elements of this partition*/
[25539]438 for(Object* & object : femmodel->elements->objects){
439 Element* element=xDynamicCast<Element*>(object);
[18001]440
441 /*Allocate all arrays*/
442 int numvertices = element->GetNumberOfVertices();
[23366]443 IssmDouble* acc = xNew<IssmDouble>(numvertices);
[18001]444 IssmDouble* evap = xNew<IssmDouble>(numvertices);
[23366]445 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
[18001]446 IssmDouble* smb = xNew<IssmDouble>(numvertices);
447
448 /*Recover Smb Components*/
[26208]449 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
450 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
451 element->GetInputListOnVertices(runoff,SmbRunoffEnum);
[18001]452
453 // loop over all vertices
[24335]454 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-runoff[v];
[18001]455
456 /*Add input to element and Free memory*/
[25379]457 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]458 xDelete<IssmDouble>(acc);
459 xDelete<IssmDouble>(evap);
460 xDelete<IssmDouble>(runoff);
461 xDelete<IssmDouble>(smb);
462 }
463
464}/*}}}*/
465void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
466
467 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
468 // INPUT parameters: ni: working size of arrays
469 // INPUT: surface accumulation (m/yr water equivalent): acc
470 // surface evaporation (m/yr water equivalent): evap
471 // surface melt (m/yr water equivalent): melt
472 // refreeze of surface melt (m/yr water equivalent): refreeze
473 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[23814]474
[18001]475 /*Loop over all the elements of this partition*/
[25539]476 for(Object* & object : femmodel->elements->objects){
477 Element* element=xDynamicCast<Element*>(object);
[18001]478
479 /*Allocate all arrays*/
480 int numvertices = element->GetNumberOfVertices();
481 IssmDouble* acc = xNew<IssmDouble>(numvertices);
[23366]482 IssmDouble* evap = xNew<IssmDouble>(numvertices);
[18001]483 IssmDouble* melt = xNew<IssmDouble>(numvertices);
484 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
485 IssmDouble* smb = xNew<IssmDouble>(numvertices);
486
487 /*Recover Smb Components*/
[26208]488 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
489 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
490 element->GetInputListOnVertices(melt,SmbMeltEnum);
491 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
[18001]492
493 // loop over all vertices
[24335]494 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
[18001]495
496 /*Add input to element and Free memory*/
[25379]497 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]498 xDelete<IssmDouble>(acc);
499 xDelete<IssmDouble>(evap);
500 xDelete<IssmDouble>(melt);
501 xDelete<IssmDouble>(refreeze);
502 xDelete<IssmDouble>(smb);
503 }
504
505}/*}}}*/
[27856]506void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
507 for(Object* & object : femmodel->elements->objects){
508 Element* element=xDynamicCast<Element*>(object);
509 element->SmbDebrisEvatt();
510 }
[27297]511}/*}}}*/
[23366]512void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
513
[25539]514 for(Object* & object : femmodel->elements->objects){
515 Element* element=xDynamicCast<Element*>(object);
[23366]516 element->SmbGradCompParameterization();
517 }
518
519}/*}}}*/
[23540]520#ifdef _HAVE_SEMIC_
[27498]521void SmbSemicx(FemModel* femmodel,int ismethod){/*{{{*/
[23540]522
[25539]523 for(Object* & object : femmodel->elements->objects){
524 Element* element=xDynamicCast<Element*>(object);
[27498]525 if (ismethod == 1) element->SmbSemicTransient(); // Inwoo's version.
526 else element->SmbSemic(); // original SmbSEMIC
[23540]527 }
528
529}/*}}}*/
530#else
531void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");}
532#endif //_HAVE_SEMIC_
Note: See TracBrowser for help on using the repository browser.