source: issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp@ 21341

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago

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