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

Last change on this file since 18584 was 18584, checked in by Mathieu Morlighem, 10 years ago

BUG: fixed m/yr m/sec conversion on SMB
File size: 12.6 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:
[18584]33 if(VerboseSolution())_printf0_(" call smb gradients module\n");
[17085]34 SmbGradientsx(femmodel);
35 break;
36 case SMBhenningEnum:
[18584]37 if(VerboseSolution())_printf0_(" call smb Henning module\n");
[17085]38 SmbHenningx(femmodel);
[17087]39 break;
[18001]40 case SMBcomponentsEnum:
[18584]41 if(VerboseSolution())_printf0_(" call smb Components module\n");
[18001]42 SmbComponentsx(femmodel);
43 break;
44 case SMBmeltcomponentsEnum:
[18584]45 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
[18001]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++){
[18521]66 Element* element=xDynamicCast<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* 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
[18266]83 /*Recover surface elevation at vertices: */
[18001]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 }
[18266]98 smb[v]=smb[v]/1000*rho_water/rho_ice; // SMB in m/y ice
[18001]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);
[17085]109 }
110
111}/*}}}*/
112void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
113
114 for(int i=0;i<femmodel->elements->Size();i++){
[18521]115 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17085]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++){
[18521]206 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17085]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
[17087]216 /*Intermediaries*/
[17403]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
[17087]234 /*Loop over all the elements of this partition*/
[17085]235 for(int i=0;i<femmodel->elements->Size();i++){
[18521]236 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17087]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
[17403]249 /*Get vertex elevation, anoma smb*/
[17087]250 z = surfacelist[v];
[17403]251 anomaly = smblistref[v];
[17087]252
[17403]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;
[17087]260 }
261 else{
[17403]262 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
[17087]263 }
[17403]264
[18584]265 /* Compute smb including anomaly,
266 correct for number of seconds in a year [s/yr]*/
267 smb = smb/yts + anomaly;
268
269
[17087]270 /*Update array accordingly*/
271 smblist[v] = smb;
272
273 }
274
275 /*Add input to element and Free memory*/
276 element->AddInput(SurfaceforcingsMassBalanceEnum,smblist,P1Enum);
277 xDelete<IssmDouble>(surfacelist);
278 xDelete<IssmDouble>(smblistref);
279 xDelete<IssmDouble>(smblist);
[17085]280 }
281
282}/*}}}*/
[18001]283void SmbComponentsx(FemModel* femmodel){/*{{{*/
284
285 // void SmbComponentsx(acc,evap,runoff,ni){
286 // INPUT parameters: ni: working size of arrays
287 // INPUT: surface accumulation (m/yr water equivalent): acc
288 // surface evaporation (m/yr water equivalent): evap
289 // surface runoff (m/yr water equivalent): runoff
290 // OUTPUT: mass-balance (m/yr ice): agd(NA)
291 int v;
292
293 /*Loop over all the elements of this partition*/
294 for(int i=0;i<femmodel->elements->Size();i++){
[18521]295 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18001]296
297 /*Allocate all arrays*/
298 int numvertices = element->GetNumberOfVertices();
299 IssmDouble* acc = xNew<IssmDouble>(numvertices);
300 IssmDouble* evap = xNew<IssmDouble>(numvertices);
301 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
302 IssmDouble* smb = xNew<IssmDouble>(numvertices);
303
304 /*Recover Smb Components*/
305 element->GetInputListOnVertices(acc,SurfaceforcingsAccumulationEnum);
306 element->GetInputListOnVertices(evap,SurfaceforcingsEvaporationEnum);
307 element->GetInputListOnVertices(runoff,SurfaceforcingsRunoffEnum);
308
309 // loop over all vertices
310 for(v=0;v<numvertices;v++){
311 smb[v]=acc[v]-evap[v]-runoff[v];
312 } //end of the loop over the vertices
313
314 /*Add input to element and Free memory*/
315 element->AddInput(SurfaceforcingsMassBalanceEnum,smb,P1Enum);
316 xDelete<IssmDouble>(acc);
317 xDelete<IssmDouble>(evap);
318 xDelete<IssmDouble>(runoff);
319 xDelete<IssmDouble>(smb);
320 }
321
322}/*}}}*/
323void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
324
325 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
326 // INPUT parameters: ni: working size of arrays
327 // INPUT: surface accumulation (m/yr water equivalent): acc
328 // surface evaporation (m/yr water equivalent): evap
329 // surface melt (m/yr water equivalent): melt
330 // refreeze of surface melt (m/yr water equivalent): refreeze
331 // OUTPUT: mass-balance (m/yr ice): agd(NA)
332 int v;
333
334 /*Loop over all the elements of this partition*/
335 for(int i=0;i<femmodel->elements->Size();i++){
[18521]336 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[18001]337
338 /*Allocate all arrays*/
339 int numvertices = element->GetNumberOfVertices();
340 IssmDouble* acc = xNew<IssmDouble>(numvertices);
341 IssmDouble* evap = xNew<IssmDouble>(numvertices);
342 IssmDouble* melt = xNew<IssmDouble>(numvertices);
343 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
344 IssmDouble* smb = xNew<IssmDouble>(numvertices);
345
346 /*Recover Smb Components*/
347 element->GetInputListOnVertices(acc,SurfaceforcingsAccumulationEnum);
348 element->GetInputListOnVertices(evap,SurfaceforcingsEvaporationEnum);
349 element->GetInputListOnVertices(melt,SurfaceforcingsMeltEnum);
350 element->GetInputListOnVertices(refreeze,SurfaceforcingsRefreezeEnum);
351
352 // loop over all vertices
353 for(v=0;v<numvertices;v++){
354 smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
355 } //end of the loop over the vertices
356
357 /*Add input to element and Free memory*/
358 element->AddInput(SurfaceforcingsMassBalanceEnum,smb,P1Enum);
359 xDelete<IssmDouble>(acc);
360 xDelete<IssmDouble>(evap);
361 xDelete<IssmDouble>(melt);
362 xDelete<IssmDouble>(refreeze);
363 xDelete<IssmDouble>(smb);
364 }
365
366}/*}}}*/
Note: See TracBrowser for help on using the repository browser.