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

Last change on this file since 19236 was 19236, checked in by schlegel, 11 years ago

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