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

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

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