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

Last change on this file since 18521 was 18521, checked in by Mathieu Morlighem, 11 years ago

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