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

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

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