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

Last change on this file since 17085 was 17085, checked in by Mathieu Morlighem, 11 years ago

CHG: merging all SMB modules into one for simplicity
File size: 5.4 KB
Line 
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);
39 default:
40 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
41 }
42
43}/*}}}*/
44
45void SmbGradientsx(FemModel* femmodel){/*{{{*/
46
47 // void SurfaceMassBalancex(hd,agd,ni){
48 // INPUT parameters: ni: working size of arrays
49 // INPUT: surface elevation (m): hd(NA)
50 // OUTPUT: mass-balance (m/yr ice): agd(NA)
51
52 for(int i=0;i<femmodel->elements->Size();i++){
53 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
54 element->SmbGradients();
55 }
56
57}/*}}}*/
58void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
59
60 for(int i=0;i<femmodel->elements->Size();i++){
61 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
62 element->Delta18oParameterization();
63 }
64
65}/*}}}*/
66void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
67
68 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
69 // note "v" prefix means 12 monthly means, ie time dimension
70 // INPUT parameters: ni: working size of arrays
71 // INPUT: surface elevation (m): hd(NA)
72 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
73 // ,NTIME)
74 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
75 // OUTPUT: mass-balance (m/yr ice): agd(NA)
76 // mean annual surface temperature (degrees C): Tsurf(NA)
77
78 int i, it, jj, itm;
79 IssmDouble DT = 0.02, sigfac, snormfac;
80 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
81 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
82 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
83 IssmDouble siglimc, siglim0, siglim0c;
84 IssmDouble tstep, tsint, tint, tstepc;
85 int NPDMAX = 1504, NPDCMAX = 1454;
86 //IssmDouble pdds[NPDMAX]={0};
87 //IssmDouble pds[NPDCMAX]={0};
88 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
89 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
90 IssmDouble tstar; // monthly mean surface temp
91
92 IssmDouble *pdds = NULL;
93 IssmDouble *pds = NULL;
94 Element *element = NULL;
95
96 pdds=xNew<IssmDouble>(NPDMAX+1);
97 pds=xNew<IssmDouble>(NPDCMAX+1);
98
99 /* initialize PDD (creation of a lookup table)*/
100 tstep = 0.1;
101 tsint = tstep*0.5;
102 sigfac = -1.0/(2.0*pow(signorm,2));
103 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
104 siglim = 2.5*signorm;
105 siglimc = 2.5*signormc;
106 siglim0 = siglim/DT + 0.5;
107 siglim0c = siglimc/DT + 0.5;
108 PDup = siglimc+PDCUT;
109
110 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
111
112 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
113 for(it = 0; it < itm; it++){
114 // tstar = REAL(it)*DT-siglim;
115 tstar = it*DT-siglim;
116 tint = tsint;
117 pddt = 0.;
118 for ( jj = 0; jj < 600; jj++){
119 if (tint > (tstar+siglim)){break;}
120 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
121 tint = tint+tstep;
122 }
123 pdds[it] = pddt*snormfac;
124 }
125 pdds[itm+1] = siglim + DT;
126
127 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
128 tstepc = 0.1;
129 tsint = PDCUT-tstepc*0.5;
130 signormc = signorm - 0.5;
131 sigfac = -1.0/(2.0*pow(signormc,2));
132 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
133 siglimc = 2.5*signormc ;
134 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
135 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
136 for(it = 0; it < itm; it++ ){
137 tstar = it*DT-siglimc;
138 // tstar = REAL(it)*DT-siglimc;
139 tint = tsint; // start against upper bound
140 pd = 0.;
141 for (jj = 0; jj < 600; jj++){
142 if (tint<(tstar-siglimc)) {break;}
143 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
144 tint = tint-tstepc;
145 }
146 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
147 }
148 pds[itm+1] = 0.;
149 // *******END initialize PDD
150
151 for(i=0;i<femmodel->elements->Size();i++){
152 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
153 element->PositiveDegreeDay(pdds,pds,signorm);
154 }
155
156 /*free ressouces: */
157 xDelete<IssmDouble>(pdds);
158 xDelete<IssmDouble>(pds);
159}/*}}}*/
160void SmbHenningx(FemModel* femmodel){/*{{{*/
161
162 // void SurfaceMassBalancex(hd,agd,ni){
163 // INPUT parameters: ni: working size of arrays
164 // INPUT: surface elevation (m): hd(NA)
165 // OUTPUT: mass-balance (m/yr ice): agd(NA)
166
167 for(int i=0;i<femmodel->elements->Size();i++){
168 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
169 }
170
171}/*}}}*/
Note: See TracBrowser for help on using the repository browser.