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

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

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