[issm-svn] r12346 - issm/trunk/src/c
lemorzad at issm.ess.uci.edu
lemorzad at issm.ess.uci.edu
Sat Jun 2 08:27:24 PDT 2012
Author: lemorzad
Date: 2012-06-02 08:27:24 -0700 (Sat, 02 Jun 2012)
New Revision: 12346
Modified:
issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
issm/trunk/src/c/objects/Elements/Tria.cpp
Log:
debug pdd methode
Modified: issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
===================================================================
--- issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-06-02 03:48:42 UTC (rev 12345)
+++ issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-06-02 15:27:24 UTC (rev 12346)
@@ -63,7 +63,7 @@
// tstar = REAL(it)*DT-siglim;
tstar = it*DT-siglim;
tint = tsint;
- pddt = 0;
+ pddt = 0.;
for ( jj = 0; jj < 600; jj++){
if (tint > (tstar+siglim)){break;}
pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
@@ -89,7 +89,7 @@
tstar = it*DT-siglimc;
// tstar = REAL(it)*DT-siglimc;
tint = tsint; // start against upper bound
- pd = 0;
+ pd = 0.;
for (jj = 0; jj < 600; jj++){
if (tint<(tstar-siglimc)) {break;}
pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
@@ -97,7 +97,7 @@
}
pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
}
- pds[itm+1] = 0;
+ pds[itm+1] = 0.;
// *******END initialize PDD
for(i=0;i<elements->Size();i++){
Modified: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp 2012-06-02 03:48:42 UTC (rev 12345)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp 2012-06-02 15:27:24 UTC (rev 12346)
@@ -2087,6 +2087,7 @@
int i,iqj,imonth;
double agd[NUMVERTICES]; // surface mass balance
+ double agd2[NUMVERTICES];
double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
double smelt[NUMVERTICES] = {0}; // yearly melt
double precrunoff[NUMVERTICES]; // yearly runoff
@@ -2096,7 +2097,7 @@
double sconv; //rhow_rain/rhoi / 12 months
double rho_water,rho_ice,density;
- double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
+ double lapser=6.5/1000., sealev=0.; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
double desfac = 0.5; // desert elevation factor
double s0p[NUMVERTICES]={0}; // should be set to elevation from precip source
double s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source
@@ -2122,10 +2123,10 @@
double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature
double Tsurf[NUMVERTICES] = {0}; // average annual temperature
- double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
- double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12];
- double deltm=1/12;
- int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
+ double h[NUMVERTICES],s[NUMVERTICES]; // ,b[NUMVERTICES]
+ double monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+ double deltm=1./12.;
+ int ismon[12]={11,0,1,2,3,4,5,6,7,8,9,10};
double snwm; // snow that could have been melted in a year.
double snwmf; // ablation factor for snow per positive degree day.
@@ -2135,33 +2136,43 @@
double supice,supcap,diffndd;
double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
double pddtj[NUMVERTICES], hmx2;
+int iii;
+ /*Recover monthly temperatures and precipitation*/
+ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
+ Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
+ GaussTria* gauss=new GaussTria();
+ double time,yts;
+ this->parameters->FindParam(&time,TimeEnum);
+ this->parameters->FindParam(&yts,ConstantsYtsEnum);
+ for(int month=0;month<12;month++) {
+ for(int iv=0;iv<NUMVERTICES;iv++) {
+ gauss->GaussVertex(iv);
+ input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
+ monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
+ input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
+ monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
+ }
+ }
+// if(this->id==1) {
+// input2->Echo();}
+ // if(this->id==1){
+ // printf("-----------------time ------------------------\n");
+ // printf("%f \n",time);
+ // printf("-----------------precipitation ------------------------\n");
+ // for (iii=0; iii<3; iii++){
+ // printf("%f %f %f %f\n",monthlyprec[iii][1],monthlyprec[iii][3],monthlyprec[iii][5],monthlyprec[iii][11]);
+ // }
+ // printf("----------------monthlytemperatures ------------------------\n");
+ // for (iii=0; iii<3; iii++){
+ // printf("%f %f %f %f %f %f\n",monthlytemperatures[iii][12],monthlytemperatures[iii][13],monthlytemperatures[iii][15],monthlytemperatures[iii][17],monthlytemperatures[iii][20],monthlytemperatures[iii][24]);
+ // }
+ // }
/*Recover info at the vertices: */
GetInputListOnVertices(&h[0],ThicknessEnum);
GetInputListOnVertices(&s[0],SurfaceEnum);
- GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
- GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
-
- /*Recover monthly temperatures*/
- Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
- GaussTria* gauss=new GaussTria();
- double time,yts;
- this->parameters->FindParam(&time,TimeEnum);
- this->parameters->FindParam(&yts,ConstantsYtsEnum);
- for(int month=0;month<12;month++){
- for(int iv=0;iv<NUMVERTICES;iv++){
- gauss->GaussVertex(iv);
- input->GetInputValue(&t[iv][month],gauss,time+month/12*yts);
- t[iv][month]=t[iv][month]-273.15; // conversion from Kelvin to celcius
- }
- }
-
- for(i=0;i<NUMVERTICES;i++)
- for(imonth=0;imonth<12;imonth++){
- t[i][imonth]=ttmp[i];
- prec[i][imonth]=prectmp[i];
- }
-
+ //GetInputListOnVertices(&monthlyprec[0],SurfaceforcingsPrecipitationEnum);
+
/*Get material parameters :*/
rho_ice=matpar->GetRhoIce();
rho_water=matpar->GetRhoFreshwater();
@@ -2179,8 +2190,8 @@
for (iqj = 0; iqj < 12; iqj++){
imonth = ismon[iqj];
for (i = 0; i < NUMVERTICES; i++){
- st=(s[i]-s0t[i])/1000;
- tstar = t[i][imonth] - lapser *max(st,sealev);
+ st=(s[i]-s0t[i])/1000.;
+ tstar = monthlytemperatures[i][imonth] - lapser *max(st,sealev);
Tsurf[i] = tstar*deltm+Tsurf[i];
/*********compute PD ****************/
@@ -2191,19 +2202,34 @@
pd = 0;}
/******exp des/elev precip reduction*******/
- sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
+ sp=(s[i]-s0p[i])/1000.; // deselev effect is wrt chng in topo
if (sp>0.0){q = exp(-desfac*sp);}
else {q = 1.0;}
- qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent per month
- qmpt= q*prec[i][imonth]*sconv;
+ qmt[i]= qmt[i] + monthlyprec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent per month
+ qmpt= q*monthlyprec[i][imonth]*sconv;
qmp[i]= qmp[i] + qmpt;
qm[i]= qm[i] + qmpt*pd;
-
+ // if(time==473040000){
+ // if(this->id==2) {
+ // printf("----------------qm----------------------\n");
+ // //printf("%f %f %f \n",qm[0],qm[1],qm[2]);
+ // printf("%d %f \n",i,qm[i]);
+ // printf("%f \n",qmpt);
+ // printf("%f \n",q);
+ // printf("%f \n",monthlyprec[i][imonth]);
+ // printf("%f \n",sconv);
+ // printf("%f \n",pd);
+ // printf("%f \n",tstar);
+ // printf("%f \n",Tsurf[i]);
+
+ // //printf("%f\n",saccu[i]);
+ // //printf("%f\n",agd[i]);
+ // }}
/*********compute PDD************/
// ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
// gaussian=T_m, so ndd=-(Tsurf-pdd)
- if (iqj>6 && iqj<10){ Tsum[i]=Tsum[i]+tstar;}
+ if (iqj>5 && iqj<9){ Tsum[i]=Tsum[i]+tstar;}
if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
else if (tstar> -siglim){
pddsig=pdds[int(tstar/DT + siglim0)];
@@ -2212,7 +2238,7 @@
else{frzndd[i] = frzndd[i] - tstar*deltm; }
}
} // end of seasonal loop
-
+
//******************************************************************
for(i=0;i<NUMVERTICES;i++){
saccu[i] = qm[i];
@@ -2223,7 +2249,7 @@
if(Tsum[i]< -1.) {
snwmf=2.65*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
- }
+ }
else if(Tsum[i]< 10){
snwmf = (0.15*Tsum[i] + 2.8)*0.001;
smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
@@ -2305,16 +2331,30 @@
//yet from site plots, can be ice free with Tsurf=-5.5C
if(Tsurf[i]<0) {
Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
-
+
agd[i] = -smelt[i]+saccu[i];
+// if(this->id==2){
+// printf("----------------melt accu----------------------\n");
+// printf("%f\n",-smelt[i]);
+// printf("%f\n",saccu[i]);
+// printf("%f\n",agd[i]);
+// }
+ agd[i] = agd[i]/yts;
pddtj[i]=pddt;
-
+ } //end of the for loop over the vertices
+
/*Update inputs*/
- this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
+ this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
// this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
- this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
-
- } //end of the for loop over the vertices
+ // this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
+
+
+ // if(this->id==2){
+ // for(iii=0;iii<3;iii++){agd2[iii]=agd[iii]*yts;}
+ // printf("----------------surface mass balance 2------------------------\n");
+ // printf("%f %f %f\n",agd2[0],agd2[1],agd2[2]);
+ // }
+
}
/*}}}*/
/*FUNCTION Tria::ProcessResultsUnits{{{1*/
More information about the issm-svn
mailing list