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

Last change on this file since 19172 was 19172, checked in by lemorzad, 10 years ago

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