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

Last change on this file since 27297 was 27297, checked in by rueckamp, 2 years ago

NEW: SMB debris model based on Mayer & Licuilli 2021
File size: 25.8 KB
RevLine 
[17085]1/*!\file SurfaceMassBalancex
[23366]2 * \brief: calculates SMB
[17085]3 */
4
[26479]5#include <config.h>
[17085]6#include "./SurfaceMassBalancex.h"
7#include "../../shared/shared.h"
8#include "../../toolkits/toolkits.h"
[23814]9#include "../modules.h"
[25379]10#include "../../classes/Inputs/TransientInput.h"
[26477]11#include "../../shared/Random/random.h"
[17085]12
[23814]13void SmbForcingx(FemModel* femmodel){/*{{{*/
14
15 // void SmbForcingx(smb,ni){
16 // INPUT parameters: ni: working size of arrays
17 // OUTPUT: mass-balance (m/yr ice): agd(NA)
18
19}/*}}}*/
[17085]20void SmbGradientsx(FemModel* femmodel){/*{{{*/
21
22 // void SurfaceMassBalancex(hd,agd,ni){
23 // INPUT parameters: ni: working size of arrays
24 // INPUT: surface elevation (m): hd(NA)
25 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[18001]26 int v;
27 IssmDouble rho_water; // density of fresh water
28 IssmDouble rho_ice; // density of ice
[21469]29 IssmDouble yts; // conversion factor year to second
[17085]30
[18001]31 /*Loop over all the elements of this partition*/
[25539]32 for(Object* & object : femmodel->elements->objects){
33 Element* element=xDynamicCast<Element*>(object);
[18001]34
35 /*Allocate all arrays*/
36 int numvertices = element->GetNumberOfVertices();
37 IssmDouble* Href = xNew<IssmDouble>(numvertices); // reference elevation from which deviations are used to calculate the SMB adjustment
38 IssmDouble* Smbref = xNew<IssmDouble>(numvertices); // reference SMB to which deviations are added
39 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // Hs-SMB relation parameter
40 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // Hs-SMB relation paremeter
41 IssmDouble* s = xNew<IssmDouble>(numvertices); // surface elevation (m)
42 IssmDouble* smb = xNew<IssmDouble>(numvertices);
43
44 /*Recover SmbGradients*/
[19527]45 element->GetInputListOnVertices(Href,SmbHrefEnum);
46 element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
47 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
48 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
[18001]49
[18266]50 /*Recover surface elevation at vertices: */
[18001]51 element->GetInputListOnVertices(s,SurfaceEnum);
52
53 /*Get material parameters :*/
[23644]54 rho_ice=element->FindParam(MaterialsRhoIceEnum);
55 rho_water=element->FindParam(MaterialsRhoFreshwaterEnum);
[18001]56
[21469]57 /* Get constants */
58 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
59
[18001]60 // loop over all vertices
61 for(v=0;v<numvertices;v++){
62 if(Smbref[v]>0){
63 smb[v]=Smbref[v]+b_pos[v]*(s[v]-Href[v]);
64 }
65 else{
66 smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
67 }
[21469]68
[18266]69 smb[v]=smb[v]/1000*rho_water/rho_ice; // SMB in m/y ice
[18001]70 } //end of the loop over the vertices
71
72 /*Add input to element and Free memory*/
[25379]73 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]74 xDelete<IssmDouble>(Href);
75 xDelete<IssmDouble>(Smbref);
76 xDelete<IssmDouble>(b_pos);
77 xDelete<IssmDouble>(b_neg);
78 xDelete<IssmDouble>(s);
79 xDelete<IssmDouble>(smb);
[17085]80 }
81
82}/*}}}*/
[21469]83void SmbGradientsElax(FemModel* femmodel){/*{{{*/
84
85 // void SurfaceMassBalancex(hd,agd,ni){
86 // INPUT parameters: ni: working size of arrays
87 // INPUT: surface elevation (m): hd(NA)
88 // OUTPUT: surface mass-balance (m/yr ice): agd(NA)
89 int v;
90
91 /*Loop over all the elements of this partition*/
[25539]92 for(Object* & object : femmodel->elements->objects){
93 Element* element=xDynamicCast<Element*>(object);
[21469]94
95 /*Allocate all arrays*/
96 int numvertices = element->GetNumberOfVertices();
97 IssmDouble* ela = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
98 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
99 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
100 IssmDouble* b_max = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
101 IssmDouble* b_min = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
102 IssmDouble* s = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
103 IssmDouble* smb = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
104
105 /*Recover ELA, SMB gradients, and caps*/
106 element->GetInputListOnVertices(ela,SmbElaEnum);
107 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
108 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
109 element->GetInputListOnVertices(b_max,SmbBMaxEnum);
110 element->GetInputListOnVertices(b_min,SmbBMinEnum);
111
112 /*Recover surface elevation at vertices: */
113 element->GetInputListOnVertices(s,SurfaceEnum);
114
115 /*Loop over all vertices, calculate SMB*/
116 for(v=0;v<numvertices;v++){
117 // if surface is above the ELA
[23366]118 if(s[v]>ela[v]){
[21469]119 smb[v]=b_pos[v]*(s[v]-ela[v]);
120 }
121 // if surface is below or equal to the ELA
122 else{
123 smb[v]=b_neg[v]*(s[v]-ela[v]);
124 }
125
126 // if SMB is larger than upper cap, set SMB to upper cap
127 if(smb[v]>b_max[v]){
128 smb[v]=b_max[v];
129 }
130 // if SMB is smaller than lower cap, set SMB to lower cap
131 if(smb[v]<b_min[v]){
132 smb[v]=b_min[v];
133 }
134 } //end of the loop over the vertices
135
136 /*Add input to element and Free memory*/
[25379]137 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[21469]138 xDelete<IssmDouble>(ela);
139 xDelete<IssmDouble>(b_pos);
140 xDelete<IssmDouble>(b_neg);
141 xDelete<IssmDouble>(b_max);
142 xDelete<IssmDouble>(b_min);
143 xDelete<IssmDouble>(s);
144 xDelete<IssmDouble>(smb);
145
146 }
147
148}/*}}}*/
[27250]149void Smbarmax(FemModel* femmodel){/*{{{*/
[26483]150
[26615]151 /*Get time parameters*/
[27250]152 IssmDouble time,dt,starttime,tstep_arma;
[26615]153 femmodel->parameters->FindParam(&time,TimeEnum);
154 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
155 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
[27250]156 femmodel->parameters->FindParam(&tstep_arma,SmbARMATimestepEnum);
[26482]157
[27250]158 /*Determine if this is a time step for the ARMA model*/
159 bool isstepforarma = false;
[26479]160
[26615]161 #ifndef _HAVE_AD_
[27250]162 if((fmod(time,tstep_arma)<fmod((time-dt),tstep_arma)) || (time<=starttime+dt) || tstep_arma==dt) isstepforarma = true;
[26615]163 #else
164 _error_("not implemented yet");
165 #endif
[26479]166
[26615]167 /*Load parameters*/
168 bool isstochastic;
169 bool issmbstochastic = false;
[27260]170 int M,N,Narlagcoefs,Nmalagcoefs,arorder,maorder,numbasins,numelevbins,my_rank;
[26615]171 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
[27250]172 femmodel->parameters->FindParam(&arorder,SmbARMAarOrderEnum);
[27260]173 femmodel->parameters->FindParam(&maorder,SmbARMAmaOrderEnum);
[26947]174 femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum);
[27250]175 IssmDouble tinit_arma;
[27246]176 IssmDouble* termconstant = NULL;
177 IssmDouble* trend = NULL;
[27250]178 IssmDouble* arlagcoefs = NULL;
[27260]179 IssmDouble* malagcoefs = NULL;
[26947]180 IssmDouble* lapserates = NULL;
181 IssmDouble* elevbins = NULL;
182 IssmDouble* refelevation = NULL;
[26526]183
[27250]184 femmodel->parameters->FindParam(&tinit_arma,SmbARMAInitialTimeEnum);
185 femmodel->parameters->FindParam(&termconstant,&M,SmbARMAconstEnum); _assert_(M==numbasins);
186 femmodel->parameters->FindParam(&trend,&M,SmbARMAtrendEnum); _assert_(M==numbasins);
187 femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,SmbARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(Narlagcoefs==arorder);
[27260]188 femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,SmbARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder);
[27250]189 femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins);
190 femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==numelevbins-1);
191 femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins);
[26481]192
[26615]193 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
[26526]194 if(isstochastic){
[26615]195 int numstochasticfields;
[26526]196 int* stochasticfields;
197 femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
198 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
199 for(int i=0;i<numstochasticfields;i++){
[27250]200 if(stochasticfields[i]==SMBarmaEnum) issmbstochastic = true;
[26615]201 }
202 xDelete<int>(stochasticfields);
203 }
[27250]204 /*Time elapsed with respect to ARMA model initial time*/
205 IssmDouble telapsed_arma = time-tinit_arma;
[26481]206
[26615]207 /*Loop over each element to compute SMB at vertices*/
208 for(Object* &object:femmodel->elements->objects){
209 Element* element = xDynamicCast<Element*>(object);
[27250]210 /*Compute ARMA*/
[27260]211 element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,issmbstochastic,SMBarmaEnum);
[26810]212 /*Compute lapse rate adjustment*/
[26947]213 element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation);
[26810]214 }
[26481]215
[26615]216 /*Cleanup*/
[27246]217 xDelete<IssmDouble>(termconstant);
218 xDelete<IssmDouble>(trend);
[27250]219 xDelete<IssmDouble>(arlagcoefs);
[27260]220 xDelete<IssmDouble>(malagcoefs);
[26947]221 xDelete<IssmDouble>(lapserates);
222 xDelete<IssmDouble>(elevbins);
[26810]223 xDelete<IssmDouble>(refelevation);
[26477]224}/*}}}*/
[17085]225void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
226
[25539]227 for(Object* & object : femmodel->elements->objects){
228 Element* element=xDynamicCast<Element*>(object);
[17085]229 element->Delta18oParameterization();
230 }
231
232}/*}}}*/
[18968]233void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
234
[25539]235 for(Object* & object : femmodel->elements->objects){
236 Element* element=xDynamicCast<Element*>(object);
[18968]237 element->MungsmtpParameterization();
238 }
239
240}/*}}}*/
[19172]241void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
242
[25539]243 for(Object* & object : femmodel->elements->objects){
244 Element* element=xDynamicCast<Element*>(object);
[19172]245 element->Delta18opdParameterization();
246 }
247
248}/*}}}*/
[17085]249void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
250
251 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
252 // note "v" prefix means 12 monthly means, ie time dimension
253 // INPUT parameters: ni: working size of arrays
254 // INPUT: surface elevation (m): hd(NA)
255 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
[23366]256 // ,NTIME)
[17085]257 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
258 // OUTPUT: mass-balance (m/yr ice): agd(NA)
259 // mean annual surface temperature (degrees C): Tsurf(NA)
260
[25539]261 int it, jj, itm;
[17085]262 IssmDouble DT = 0.02, sigfac, snormfac;
[23366]263 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
[17085]264 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
265 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
266 IssmDouble siglimc, siglim0, siglim0c;
267 IssmDouble tstep, tsint, tint, tstepc;
268 int NPDMAX = 1504, NPDCMAX = 1454;
[23366]269 //IssmDouble pdds[NPDMAX]={0};
[17085]270 //IssmDouble pds[NPDCMAX]={0};
271 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
272 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
273 IssmDouble tstar; // monthly mean surface temp
274
[18968]275 bool ismungsm;
[22448]276 bool issetpddfac;
[18968]277
[17085]278 IssmDouble *pdds = NULL;
279 IssmDouble *pds = NULL;
280 Element *element = NULL;
281
[23366]282 pdds=xNew<IssmDouble>(NPDMAX+1);
283 pds=xNew<IssmDouble>(NPDCMAX+1);
[17085]284
[18968]285 // Get ismungsm parameter
[19527]286 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
[18968]287
[22448]288 // Get issetpddfac parameter
289 femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
290
[17085]291 /* initialize PDD (creation of a lookup table)*/
292 tstep = 0.1;
293 tsint = tstep*0.5;
294 sigfac = -1.0/(2.0*pow(signorm,2));
295 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
296 siglim = 2.5*signorm;
297 siglimc = 2.5*signormc;
298 siglim0 = siglim/DT + 0.5;
299 siglim0c = siglimc/DT + 0.5;
300 PDup = siglimc+PDCUT;
301
302 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
303
304 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
[23366]305 for(it = 0; it < itm; it++){
[17085]306 // tstar = REAL(it)*DT-siglim;
307 tstar = it*DT-siglim;
308 tint = tsint;
309 pddt = 0.;
310 for ( jj = 0; jj < 600; jj++){
311 if (tint > (tstar+siglim)){break;}
312 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
313 tint = tint+tstep;
314 }
315 pdds[it] = pddt*snormfac;
316 }
317 pdds[itm+1] = siglim + DT;
318
319 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
320 tstepc = 0.1;
321 tsint = PDCUT-tstepc*0.5;
322 signormc = signorm - 0.5;
323 sigfac = -1.0/(2.0*pow(signormc,2));
324 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
325 siglimc = 2.5*signormc ;
326 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
327 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
328 for(it = 0; it < itm; it++ ){
329 tstar = it*DT-siglimc;
330 // tstar = REAL(it)*DT-siglimc;
331 tint = tsint; // start against upper bound
332 pd = 0.;
333 for (jj = 0; jj < 600; jj++){
334 if (tint<(tstar-siglimc)) {break;}
335 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
336 tint = tint-tstepc;
337 }
338 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
339 }
340 pds[itm+1] = 0.;
341 // *******END initialize PDD
342
[25539]343 for(Object* & object : femmodel->elements->objects){
344 element=xDynamicCast<Element*>(object);
[22448]345 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
[17085]346 }
347 /*free ressouces: */
348 xDelete<IssmDouble>(pdds);
349 xDelete<IssmDouble>(pds);
350}/*}}}*/
[23317]351void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
[23366]352
[23317]353 bool isfirnwarming;
[23328]354 femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
[23366]355
[25539]356 for(Object* & object : femmodel->elements->objects){
357 Element* element=xDynamicCast<Element*>(object);
[23317]358 element->PositiveDegreeDaySicopolis(isfirnwarming);
359 }
360
361}/*}}}*/
[17085]362void SmbHenningx(FemModel* femmodel){/*{{{*/
363
[17087]364 /*Intermediaries*/
[17403]365 IssmDouble z_critical = 1675.;
366 IssmDouble dz = 0;
367 IssmDouble a = -15.86;
368 IssmDouble b = 0.00969;
369 IssmDouble c = -0.235;
370 IssmDouble f = 1.;
371 IssmDouble g = -0.0011;
372 IssmDouble h = -1.54e-5;
373 IssmDouble smb,smbref,anomaly,yts,z;
[22249]374
375 /* Get constants */
376 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
377 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
378 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
379 /*Mathieu original*/
380 /*IssmDouble smb,smbref,z;*/
381
[17087]382 /*Loop over all the elements of this partition*/
[25539]383 for(Object* & object : femmodel->elements->objects){
384 Element* element=xDynamicCast<Element*>(object);
[17087]385
386 /*Get reference SMB (uncorrected) and allocate all arrays*/
387 int numvertices = element->GetNumberOfVertices();
388 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
389 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
390 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
391 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
[19527]392 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
[17087]393
394 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
395 for(int v=0;v<numvertices;v++){
396
[17403]397 /*Get vertex elevation, anoma smb*/
[17087]398 z = surfacelist[v];
[17403]399 anomaly = smblistref[v];
[17087]400
[22249]401 /* Henning edited acc. to Riannes equations*/
402 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
403 z_critical = z_critical + dz;
404
405 /* Calculate smb acc. to the surface elevation z */
406 if(z<z_critical){
[17403]407 smb = a + b*z + c;
[17087]408 }
409 else{
[22249]410 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
[17087]411 }
[22249]412
[18584]413 /* Compute smb including anomaly,
414 correct for number of seconds in a year [s/yr]*/
415 smb = smb/yts + anomaly;
416
[17087]417 /*Update array accordingly*/
418 smblist[v] = smb;
419
420 }
421
422 /*Add input to element and Free memory*/
[25379]423 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
[17087]424 xDelete<IssmDouble>(surfacelist);
425 xDelete<IssmDouble>(smblistref);
426 xDelete<IssmDouble>(smblist);
[17085]427 }
428
429}/*}}}*/
[18001]430void SmbComponentsx(FemModel* femmodel){/*{{{*/
431
432 // void SmbComponentsx(acc,evap,runoff,ni){
433 // INPUT parameters: ni: working size of arrays
434 // INPUT: surface accumulation (m/yr water equivalent): acc
435 // surface evaporation (m/yr water equivalent): evap
436 // surface runoff (m/yr water equivalent): runoff
437 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[23814]438
[18001]439 /*Loop over all the elements of this partition*/
[25539]440 for(Object* & object : femmodel->elements->objects){
441 Element* element=xDynamicCast<Element*>(object);
[18001]442
443 /*Allocate all arrays*/
444 int numvertices = element->GetNumberOfVertices();
[23366]445 IssmDouble* acc = xNew<IssmDouble>(numvertices);
[18001]446 IssmDouble* evap = xNew<IssmDouble>(numvertices);
[23366]447 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
[18001]448 IssmDouble* smb = xNew<IssmDouble>(numvertices);
449
450 /*Recover Smb Components*/
[26208]451 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
452 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
453 element->GetInputListOnVertices(runoff,SmbRunoffEnum);
[18001]454
455 // loop over all vertices
[24335]456 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-runoff[v];
[18001]457
458 /*Add input to element and Free memory*/
[25379]459 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]460 xDelete<IssmDouble>(acc);
461 xDelete<IssmDouble>(evap);
462 xDelete<IssmDouble>(runoff);
463 xDelete<IssmDouble>(smb);
464 }
465
466}/*}}}*/
467void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
468
469 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
470 // INPUT parameters: ni: working size of arrays
471 // INPUT: surface accumulation (m/yr water equivalent): acc
472 // surface evaporation (m/yr water equivalent): evap
473 // surface melt (m/yr water equivalent): melt
474 // refreeze of surface melt (m/yr water equivalent): refreeze
475 // OUTPUT: mass-balance (m/yr ice): agd(NA)
[23814]476
[18001]477 /*Loop over all the elements of this partition*/
[25539]478 for(Object* & object : femmodel->elements->objects){
479 Element* element=xDynamicCast<Element*>(object);
[18001]480
481 /*Allocate all arrays*/
482 int numvertices = element->GetNumberOfVertices();
483 IssmDouble* acc = xNew<IssmDouble>(numvertices);
[23366]484 IssmDouble* evap = xNew<IssmDouble>(numvertices);
[18001]485 IssmDouble* melt = xNew<IssmDouble>(numvertices);
486 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
487 IssmDouble* smb = xNew<IssmDouble>(numvertices);
488
489 /*Recover Smb Components*/
[26208]490 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
491 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
492 element->GetInputListOnVertices(melt,SmbMeltEnum);
493 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
[18001]494
495 // loop over all vertices
[24335]496 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
[18001]497
498 /*Add input to element and Free memory*/
[25379]499 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
[18001]500 xDelete<IssmDouble>(acc);
501 xDelete<IssmDouble>(evap);
502 xDelete<IssmDouble>(melt);
503 xDelete<IssmDouble>(refreeze);
504 xDelete<IssmDouble>(smb);
505 }
506
507}/*}}}*/
[27297]508void SmbDebrisMLx(FemModel* femmodel){/*{{{*/
509
510 // The function is based on:
511 // Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2
512 // Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7
513 // function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90
514
515 /*Intermediaries*/
516 // altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002)
517 IssmDouble LW=2.9; // W/m^2 /100m 2.9
518 IssmDouble SW=1.3; // W/m^2 /100m 1.3
519 IssmDouble HumidityG=0; // % /100m rough estimate
520 IssmDouble AirTemp=0.7; // C /100m
521 IssmDouble WindSpeed=0.02; // m/s /100m rough estimate 0.2
522
523 // accumulation follows a linear increase above the ELA up to a plateau
524 IssmDouble AccG=0.1; // m w.e. /100m
525 IssmDouble AccMax=1.; // m w.e.
526 IssmDouble ReferenceElevation=2200.; // m M&L
527 IssmDouble AblationDays=120.; //
528
529 IssmDouble In=100.; // Wm^-2 incoming long wave
530 IssmDouble Q=500.; // Wm^-2 incoming short wave
531 IssmDouble K=0.585; // Wm^-1K^-1 thermal conductivity 0.585
532 IssmDouble Qm=0.0012; // kg m^-3 measured humiditiy level
533 IssmDouble Qh=0.006 ; // kg m^-3 saturated humidity level
534 IssmDouble Tm=2.; // C air temperature
535 IssmDouble Rhoaa=1.22; // kgm^-3 air densitiy
536 IssmDouble Um=1.5; // ms^-1 measured wind speed
537 IssmDouble Xm=1.5; // ms^-1 measurement height
538 IssmDouble Xr=0.01; // ms^-1 surface roughness 0.01
539 IssmDouble Alphad=0.07; // debris albedo 0.07
540 IssmDouble Alphai=0.4; // ice ablbedo
541 IssmDouble Ustar=0.16; // ms^-1 friction velocity 0.16
542 IssmDouble Ca=1000.; // jkg^-1K^-1 specific heat capacity of air
543 IssmDouble Lm;//=3.34E+05; // jkg^-1K^-1 latent heat of ice melt
544 IssmDouble Lv=2.50E+06; // jkg^-1K^-1 latent heat of evaporation
545 IssmDouble Tf=273.; // K water freeezing temperature
546 IssmDouble Eps=0.95; // thermal emissivity
547 IssmDouble Rhoi=900.; // kgm^-3 ice density
548 IssmDouble Sigma=5.67E-08; // Wm^-2K^-4 Stefan Boltzmann constant
549 IssmDouble Kstar=0.4; // von kármán constant
550 IssmDouble Gamma=180.; // m^-1 wind speed attenuation 234
551 IssmDouble PhiD;//=0.005; // debris packing fraction 0.01
552 IssmDouble Humidity=0.2; // relative humidity
553
554 IssmDouble smb,yts,z,debris;
555 IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
556
557 /*Get material parameters and constants */
558 //femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as benchmark was run with different densities for debris and FS
559 femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum);
560 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
561 PhiD=0.;
562 //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
563
564 /* Loop over all the elements of this partition */
565 for(Object* & object : femmodel->elements->objects){
566 Element* element=xDynamicCast<Element*>(object);
567
568 /* Allocate all arrays */
569 int numvertices=element->GetNumberOfVertices();
570 IssmDouble* surfacelist=xNew<IssmDouble>(numvertices);
571 IssmDouble* smb=xNew<IssmDouble>(numvertices);
572 IssmDouble* debriscover=xNew<IssmDouble>(numvertices);
573 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
574
575 /* Get inputs */
576 element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
577
578 /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
579 for(int v=0;v<numvertices;v++){
580
581 /*Get vertex elevation */
582 z=surfacelist[v];
583
584 /* Get debris cover */
585 debris=debriscover[v];
586
587 /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
588 IssmDouble n=debris/dk;
589 IssmDouble nmax=1000;
590 IssmDouble Alphaeff;
591 if(n>nmax){
592 Alphaeff=Alphad;
593 } else {
594 Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
595 }*/
596
597 // M&L
598 IssmDouble Alphaeff=Alphad;
599
600 /* compute smb */
601 for (int ismb=0;ismb<2;ismb++){
602 if(ismb==0){
603 // calc a reference smb to identify accum and melt region; debris only develops in ablation area
604 debris=0.;
605 }else{
606 // only in the meltregime debris develops
607 if(-MassBalanceCmDayDebris<0) debris=debriscover[v];
608 }
609 MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
610 (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+
611 (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
612 WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z-
613 ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+
614 ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
615 WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
616 K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z-
617 ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)*
618 Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.)
619 -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
620 }
621
622 /* account form ablation days, and convert to m/s */
623 MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts;
624
625 /*Update array accordingly*/
626 smb[v]=MassBalanceMYearDebris;
627 }
628
629 /*Add input to element and Free memory*/
630 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
631 xDelete<IssmDouble>(surfacelist);
632 xDelete<IssmDouble>(smb);
633 xDelete<IssmDouble>(debriscover);
634 }
635}/*}}}*/
[23366]636void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
637
[25539]638 for(Object* & object : femmodel->elements->objects){
639 Element* element=xDynamicCast<Element*>(object);
[23366]640 element->SmbGradCompParameterization();
641 }
642
643}/*}}}*/
[23540]644#ifdef _HAVE_SEMIC_
645void SmbSemicx(FemModel* femmodel){/*{{{*/
646
[25539]647 for(Object* & object : femmodel->elements->objects){
648 Element* element=xDynamicCast<Element*>(object);
[23540]649 element->SmbSemic();
650 }
651
652}/*}}}*/
653#else
654void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");}
655#endif //_HAVE_SEMIC_
Note: See TracBrowser for help on using the repository browser.