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

Last change on this file since 17403 was 17403, checked in by hakesson, 11 years ago

CHG: updating Henning's surface mass balance scheme
File size: 7.3 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:
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);
[17087]39 break;
[17085]40 default:
41 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
42 }
43
44}/*}}}*/
45
46void SmbGradientsx(FemModel* femmodel){/*{{{*/
47
48 // void SurfaceMassBalancex(hd,agd,ni){
49 // INPUT parameters: ni: working size of arrays
50 // INPUT: surface elevation (m): hd(NA)
51 // OUTPUT: mass-balance (m/yr ice): agd(NA)
52
53 for(int i=0;i<femmodel->elements->Size();i++){
54 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
55 element->SmbGradients();
56 }
57
58}/*}}}*/
59void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
60
61 for(int i=0;i<femmodel->elements->Size();i++){
62 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
63 element->Delta18oParameterization();
64 }
65
66}/*}}}*/
67void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
68
69 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
70 // note "v" prefix means 12 monthly means, ie time dimension
71 // INPUT parameters: ni: working size of arrays
72 // INPUT: surface elevation (m): hd(NA)
73 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
74 // ,NTIME)
75 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
76 // OUTPUT: mass-balance (m/yr ice): agd(NA)
77 // mean annual surface temperature (degrees C): Tsurf(NA)
78
79 int i, it, jj, itm;
80 IssmDouble DT = 0.02, sigfac, snormfac;
81 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
82 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
83 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
84 IssmDouble siglimc, siglim0, siglim0c;
85 IssmDouble tstep, tsint, tint, tstepc;
86 int NPDMAX = 1504, NPDCMAX = 1454;
87 //IssmDouble pdds[NPDMAX]={0};
88 //IssmDouble pds[NPDCMAX]={0};
89 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
90 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
91 IssmDouble tstar; // monthly mean surface temp
92
93 IssmDouble *pdds = NULL;
94 IssmDouble *pds = NULL;
95 Element *element = NULL;
96
97 pdds=xNew<IssmDouble>(NPDMAX+1);
98 pds=xNew<IssmDouble>(NPDCMAX+1);
99
100 /* initialize PDD (creation of a lookup table)*/
101 tstep = 0.1;
102 tsint = tstep*0.5;
103 sigfac = -1.0/(2.0*pow(signorm,2));
104 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
105 siglim = 2.5*signorm;
106 siglimc = 2.5*signormc;
107 siglim0 = siglim/DT + 0.5;
108 siglim0c = siglimc/DT + 0.5;
109 PDup = siglimc+PDCUT;
110
111 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
112
113 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
114 for(it = 0; it < itm; it++){
115 // tstar = REAL(it)*DT-siglim;
116 tstar = it*DT-siglim;
117 tint = tsint;
118 pddt = 0.;
119 for ( jj = 0; jj < 600; jj++){
120 if (tint > (tstar+siglim)){break;}
121 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
122 tint = tint+tstep;
123 }
124 pdds[it] = pddt*snormfac;
125 }
126 pdds[itm+1] = siglim + DT;
127
128 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
129 tstepc = 0.1;
130 tsint = PDCUT-tstepc*0.5;
131 signormc = signorm - 0.5;
132 sigfac = -1.0/(2.0*pow(signormc,2));
133 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
134 siglimc = 2.5*signormc ;
135 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
136 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
137 for(it = 0; it < itm; it++ ){
138 tstar = it*DT-siglimc;
139 // tstar = REAL(it)*DT-siglimc;
140 tint = tsint; // start against upper bound
141 pd = 0.;
142 for (jj = 0; jj < 600; jj++){
143 if (tint<(tstar-siglimc)) {break;}
144 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
145 tint = tint-tstepc;
146 }
147 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
148 }
149 pds[itm+1] = 0.;
150 // *******END initialize PDD
151
152 for(i=0;i<femmodel->elements->Size();i++){
153 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
154 element->PositiveDegreeDay(pdds,pds,signorm);
155 }
156
157 /*free ressouces: */
158 xDelete<IssmDouble>(pdds);
159 xDelete<IssmDouble>(pds);
160}/*}}}*/
161void SmbHenningx(FemModel* femmodel){/*{{{*/
162
[17087]163 /*Intermediaries*/
[17403]164 IssmDouble z_critical = 1675.;
165 IssmDouble dz = 0;
166 IssmDouble a = -15.86;
167 IssmDouble b = 0.00969;
168 IssmDouble c = -0.235;
169 IssmDouble f = 1.;
170 IssmDouble g = -0.0011;
171 IssmDouble h = -1.54e-5;
172 IssmDouble smb,smbref,anomaly,yts,z;
173
174 /* Get constants */
175 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
176 /*iomodel->Constant(&yts,ConstantsYtsEnum);*/
177 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
178 /*Mathieu original*/
179 /*IssmDouble smb,smbref,z;*/
180
[17087]181 /*Loop over all the elements of this partition*/
[17085]182 for(int i=0;i<femmodel->elements->Size();i++){
183 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
[17087]184
185 /*Get reference SMB (uncorrected) and allocate all arrays*/
186 int numvertices = element->GetNumberOfVertices();
187 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
188 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
189 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
190 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
191 element->GetInputListOnVertices(smblistref,SurfaceforcingsSmbrefEnum);
192
193 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
194 for(int v=0;v<numvertices;v++){
195
[17403]196 /*Get vertex elevation, anoma smb*/
[17087]197 z = surfacelist[v];
[17403]198 anomaly = smblistref[v];
[17087]199
[17403]200 /* Henning edited acc. to Riannes equations*/
201 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
202 z_critical = z_critical + dz;
203
204 /* Calculate smb acc. to the surface elevation z */
205 if(z<z_critical){
206 smb = a + b*z + c;
[17087]207 }
208 else{
[17403]209 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
[17087]210 }
[17403]211
212 /* Compute smb including anomaly,
213 correct for number of seconds in a year [s/yr]*/
214 smb = smb + anomaly*yts;
215
[17087]216 /*Update array accordingly*/
217 smblist[v] = smb;
218
219 }
220
221 /*Add input to element and Free memory*/
222 element->AddInput(SurfaceforcingsMassBalanceEnum,smblist,P1Enum);
223 xDelete<IssmDouble>(surfacelist);
224 xDelete<IssmDouble>(smblistref);
225 xDelete<IssmDouble>(smblist);
[17085]226 }
227
228}/*}}}*/
Note: See TracBrowser for help on using the repository browser.