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

Last change on this file since 23814 was 23814, checked in by schlegel, 6 years ago

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