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

Last change on this file since 22249 was 22249, checked in by schlegel, 7 years ago

CHG: update Gemb code to use tolerances on all absolute comparisons, make stricter archive tolerances, Smb output is now in m/yr ice so that it can be run with issm mass transport on
File size: 14.2 KB
RevLine 
[17085]1/*!\file SurfaceMassBalancex
2 * \brief: calculates SMB
3 */
4
5#include "./SurfaceMassBalancex.h"
6#include "../../shared/shared.h"
7#include "../../toolkits/toolkits.h"
8
9
10void SmbGradientsx(FemModel* femmodel){/*{{{*/
11
12 // void SurfaceMassBalancex(hd,agd,ni){
13 // INPUT parameters: ni: working size of arrays
14 // INPUT: surface elevation (m): hd(NA)
15 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[18001]16 int v;
17 IssmDouble rho_water; // density of fresh water
18 IssmDouble rho_ice; // density of ice
[21469]19 IssmDouble yts; // conversion factor year to second
[17085]20
[18001]21 /*Loop over all the elements of this partition*/
[17085]22 for(int i=0;i<femmodel->elements->Size();i++){
[18521]23 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18001]24
25 /*Allocate all arrays*/
26 int numvertices = element->GetNumberOfVertices();
27 IssmDouble* Href = xNew<IssmDouble>(numvertices); // reference elevation from which deviations are used to calculate the SMB adjustment
28 IssmDouble* Smbref = xNew<IssmDouble>(numvertices); // reference SMB to which deviations are added
29 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // Hs-SMB relation parameter
30 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // Hs-SMB relation paremeter
31 IssmDouble* s = xNew<IssmDouble>(numvertices); // surface elevation (m)
32 IssmDouble* smb = xNew<IssmDouble>(numvertices);
33
34 /*Recover SmbGradients*/
[19527]35 element->GetInputListOnVertices(Href,SmbHrefEnum);
36 element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
37 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
38 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
[18001]39
[18266]40 /*Recover surface elevation at vertices: */
[18001]41 element->GetInputListOnVertices(s,SurfaceEnum);
42
43 /*Get material parameters :*/
[18946]44 rho_ice=element->matpar->GetMaterialParameter(MaterialsRhoIceEnum);
[18949]45 rho_water=element->matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
[18001]46
[21469]47 /* Get constants */
48 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
49
[18001]50 // loop over all vertices
51 for(v=0;v<numvertices;v++){
52 if(Smbref[v]>0){
53 smb[v]=Smbref[v]+b_pos[v]*(s[v]-Href[v]);
54 }
55 else{
56 smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
57 }
[21469]58
[18266]59 smb[v]=smb[v]/1000*rho_water/rho_ice; // SMB in m/y ice
[18001]60 } //end of the loop over the vertices
61
62 /*Add input to element and Free memory*/
[19527]63 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]64 xDelete<IssmDouble>(Href);
65 xDelete<IssmDouble>(Smbref);
66 xDelete<IssmDouble>(b_pos);
67 xDelete<IssmDouble>(b_neg);
68 xDelete<IssmDouble>(s);
69 xDelete<IssmDouble>(smb);
[17085]70 }
71
72}/*}}}*/
[21469]73void SmbGradientsElax(FemModel* femmodel){/*{{{*/
74
75 // void SurfaceMassBalancex(hd,agd,ni){
76 // INPUT parameters: ni: working size of arrays
77 // INPUT: surface elevation (m): hd(NA)
78 // OUTPUT: surface mass-balance (m/yr ice): agd(NA)
79 int v;
80
81 /*Loop over all the elements of this partition*/
82 for(int i=0;i<femmodel->elements->Size();i++){
83 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
84
85 /*Allocate all arrays*/
86 int numvertices = element->GetNumberOfVertices();
87 IssmDouble* ela = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
88 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
89 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
90 IssmDouble* b_max = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
91 IssmDouble* b_min = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
92 IssmDouble* s = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
93 IssmDouble* smb = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
94
95 /*Recover ELA, SMB gradients, and caps*/
96 element->GetInputListOnVertices(ela,SmbElaEnum);
97 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
98 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
99 element->GetInputListOnVertices(b_max,SmbBMaxEnum);
100 element->GetInputListOnVertices(b_min,SmbBMinEnum);
101
102 /*Recover surface elevation at vertices: */
103 element->GetInputListOnVertices(s,SurfaceEnum);
104
105 /*Loop over all vertices, calculate SMB*/
106 for(v=0;v<numvertices;v++){
107 // if surface is above the ELA
108 if(s[v]>ela[v]){
109 smb[v]=b_pos[v]*(s[v]-ela[v]);
110 }
111 // if surface is below or equal to the ELA
112 else{
113 smb[v]=b_neg[v]*(s[v]-ela[v]);
114 }
115
116 // if SMB is larger than upper cap, set SMB to upper cap
117 if(smb[v]>b_max[v]){
118 smb[v]=b_max[v];
119 }
120 // if SMB is smaller than lower cap, set SMB to lower cap
121 if(smb[v]<b_min[v]){
122 smb[v]=b_min[v];
123 }
124 } //end of the loop over the vertices
125
126 /*Add input to element and Free memory*/
127 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
128 xDelete<IssmDouble>(ela);
129 xDelete<IssmDouble>(b_pos);
130 xDelete<IssmDouble>(b_neg);
131 xDelete<IssmDouble>(b_max);
132 xDelete<IssmDouble>(b_min);
133 xDelete<IssmDouble>(s);
134 xDelete<IssmDouble>(smb);
135
136 }
137
138}/*}}}*/
[17085]139void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
140
141 for(int i=0;i<femmodel->elements->Size();i++){
[18521]142 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17085]143 element->Delta18oParameterization();
144 }
145
146}/*}}}*/
[18968]147void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
148
149 for(int i=0;i<femmodel->elements->Size();i++){
150 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
151 element->MungsmtpParameterization();
152 }
153
154}/*}}}*/
[19172]155void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
156
157 for(int i=0;i<femmodel->elements->Size();i++){
158 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
159 element->Delta18opdParameterization();
160 }
161
162}/*}}}*/
[17085]163void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
164
165 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
166 // note "v" prefix means 12 monthly means, ie time dimension
167 // INPUT parameters: ni: working size of arrays
168 // INPUT: surface elevation (m): hd(NA)
169 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
170 // ,NTIME)
171 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
172 // OUTPUT: mass-balance (m/yr ice): agd(NA)
173 // mean annual surface temperature (degrees C): Tsurf(NA)
174
175 int i, it, jj, itm;
176 IssmDouble DT = 0.02, sigfac, snormfac;
177 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
178 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
179 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
180 IssmDouble siglimc, siglim0, siglim0c;
181 IssmDouble tstep, tsint, tint, tstepc;
182 int NPDMAX = 1504, NPDCMAX = 1454;
183 //IssmDouble pdds[NPDMAX]={0};
184 //IssmDouble pds[NPDCMAX]={0};
185 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
186 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
187 IssmDouble tstar; // monthly mean surface temp
188
[18968]189 bool ismungsm;
190
[17085]191 IssmDouble *pdds = NULL;
192 IssmDouble *pds = NULL;
193 Element *element = NULL;
194
195 pdds=xNew<IssmDouble>(NPDMAX+1);
196 pds=xNew<IssmDouble>(NPDCMAX+1);
197
[18968]198 // Get ismungsm parameter
[19527]199 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
[18968]200
[17085]201 /* initialize PDD (creation of a lookup table)*/
202 tstep = 0.1;
203 tsint = tstep*0.5;
204 sigfac = -1.0/(2.0*pow(signorm,2));
205 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
206 siglim = 2.5*signorm;
207 siglimc = 2.5*signormc;
208 siglim0 = siglim/DT + 0.5;
209 siglim0c = siglimc/DT + 0.5;
210 PDup = siglimc+PDCUT;
211
212 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
213
214 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
215 for(it = 0; it < itm; it++){
216 // tstar = REAL(it)*DT-siglim;
217 tstar = it*DT-siglim;
218 tint = tsint;
219 pddt = 0.;
220 for ( jj = 0; jj < 600; jj++){
221 if (tint > (tstar+siglim)){break;}
222 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
223 tint = tint+tstep;
224 }
225 pdds[it] = pddt*snormfac;
226 }
227 pdds[itm+1] = siglim + DT;
228
229 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
230 tstepc = 0.1;
231 tsint = PDCUT-tstepc*0.5;
232 signormc = signorm - 0.5;
233 sigfac = -1.0/(2.0*pow(signormc,2));
234 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
235 siglimc = 2.5*signormc ;
236 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
237 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
238 for(it = 0; it < itm; it++ ){
239 tstar = it*DT-siglimc;
240 // tstar = REAL(it)*DT-siglimc;
241 tint = tsint; // start against upper bound
242 pd = 0.;
243 for (jj = 0; jj < 600; jj++){
244 if (tint<(tstar-siglimc)) {break;}
245 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
246 tint = tint-tstepc;
247 }
248 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
249 }
250 pds[itm+1] = 0.;
251 // *******END initialize PDD
252
253 for(i=0;i<femmodel->elements->Size();i++){
[18521]254 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18968]255 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm);
[17085]256 }
257 /*free ressouces: */
258 xDelete<IssmDouble>(pdds);
259 xDelete<IssmDouble>(pds);
260}/*}}}*/
261void SmbHenningx(FemModel* femmodel){/*{{{*/
262
[17087]263 /*Intermediaries*/
[17403]264 IssmDouble z_critical = 1675.;
265 IssmDouble dz = 0;
266 IssmDouble a = -15.86;
267 IssmDouble b = 0.00969;
268 IssmDouble c = -0.235;
269 IssmDouble f = 1.;
270 IssmDouble g = -0.0011;
271 IssmDouble h = -1.54e-5;
272 IssmDouble smb,smbref,anomaly,yts,z;
[22249]273
274 /* Get constants */
275 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
276 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
277 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
278 /*Mathieu original*/
279 /*IssmDouble smb,smbref,z;*/
280
[17087]281 /*Loop over all the elements of this partition*/
[17085]282 for(int i=0;i<femmodel->elements->Size();i++){
[18521]283 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17087]284
285 /*Get reference SMB (uncorrected) and allocate all arrays*/
286 int numvertices = element->GetNumberOfVertices();
287 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
288 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
289 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
290 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
[19527]291 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
[17087]292
293 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
294 for(int v=0;v<numvertices;v++){
295
[17403]296 /*Get vertex elevation, anoma smb*/
[17087]297 z = surfacelist[v];
[17403]298 anomaly = smblistref[v];
[17087]299
[22249]300 /* Henning edited acc. to Riannes equations*/
301 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
302 z_critical = z_critical + dz;
303
304 /* Calculate smb acc. to the surface elevation z */
305 if(z<z_critical){
[17403]306 smb = a + b*z + c;
[17087]307 }
308 else{
[22249]309 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
[17087]310 }
[22249]311
[18584]312 /* Compute smb including anomaly,
313 correct for number of seconds in a year [s/yr]*/
314 smb = smb/yts + anomaly;
315
316
[17087]317 /*Update array accordingly*/
318 smblist[v] = smb;
319
320 }
321
322 /*Add input to element and Free memory*/
[19527]323 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
[17087]324 xDelete<IssmDouble>(surfacelist);
325 xDelete<IssmDouble>(smblistref);
326 xDelete<IssmDouble>(smblist);
[17085]327 }
328
329}/*}}}*/
[18001]330void SmbComponentsx(FemModel* femmodel){/*{{{*/
331
332 // void SmbComponentsx(acc,evap,runoff,ni){
333 // INPUT parameters: ni: working size of arrays
334 // INPUT: surface accumulation (m/yr water equivalent): acc
335 // surface evaporation (m/yr water equivalent): evap
336 // surface runoff (m/yr water equivalent): runoff
337 // OUTPUT: mass-balance (m/yr ice): agd(NA)
338 int v;
339
340 /*Loop over all the elements of this partition*/
341 for(int i=0;i<femmodel->elements->Size();i++){
[18521]342 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18001]343
344 /*Allocate all arrays*/
345 int numvertices = element->GetNumberOfVertices();
346 IssmDouble* acc = xNew<IssmDouble>(numvertices);
347 IssmDouble* evap = xNew<IssmDouble>(numvertices);
348 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
349 IssmDouble* smb = xNew<IssmDouble>(numvertices);
350
351 /*Recover Smb Components*/
[19527]352 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
353 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
354 element->GetInputListOnVertices(runoff,SmbRunoffEnum);
[18001]355
356 // loop over all vertices
357 for(v=0;v<numvertices;v++){
358 smb[v]=acc[v]-evap[v]-runoff[v];
359 } //end of the loop over the vertices
360
361 /*Add input to element and Free memory*/
[19527]362 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]363 xDelete<IssmDouble>(acc);
364 xDelete<IssmDouble>(evap);
365 xDelete<IssmDouble>(runoff);
366 xDelete<IssmDouble>(smb);
367 }
368
369}/*}}}*/
370void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
371
372 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
373 // INPUT parameters: ni: working size of arrays
374 // INPUT: surface accumulation (m/yr water equivalent): acc
375 // surface evaporation (m/yr water equivalent): evap
376 // surface melt (m/yr water equivalent): melt
377 // refreeze of surface melt (m/yr water equivalent): refreeze
378 // OUTPUT: mass-balance (m/yr ice): agd(NA)
379 int v;
380
381 /*Loop over all the elements of this partition*/
382 for(int i=0;i<femmodel->elements->Size();i++){
[18521]383 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18001]384
385 /*Allocate all arrays*/
386 int numvertices = element->GetNumberOfVertices();
387 IssmDouble* acc = xNew<IssmDouble>(numvertices);
388 IssmDouble* evap = xNew<IssmDouble>(numvertices);
389 IssmDouble* melt = xNew<IssmDouble>(numvertices);
390 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
391 IssmDouble* smb = xNew<IssmDouble>(numvertices);
392
393 /*Recover Smb Components*/
[19527]394 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
395 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
396 element->GetInputListOnVertices(melt,SmbMeltEnum);
397 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
[18001]398
399 // loop over all vertices
400 for(v=0;v<numvertices;v++){
401 smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
402 } //end of the loop over the vertices
403
404 /*Add input to element and Free memory*/
[19527]405 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]406 xDelete<IssmDouble>(acc);
407 xDelete<IssmDouble>(evap);
408 xDelete<IssmDouble>(melt);
409 xDelete<IssmDouble>(refreeze);
410 xDelete<IssmDouble>(smb);
411 }
412
413}/*}}}*/
Note: See TracBrowser for help on using the repository browser.