[issm-svn] r12362 - issm/trunk/src
lemorzad at issm.ess.uci.edu
lemorzad at issm.ess.uci.edu
Mon Jun 4 14:33:56 PDT 2012
Author: lemorzad
Date: 2012-06-04 14:33:55 -0700 (Mon, 04 Jun 2012)
New Revision: 12362
Modified:
issm/trunk/src/c/objects/Elements/Penta.cpp
issm/trunk/src/c/objects/Elements/Tria.cpp
issm/trunk/src/m/model/extrude.m
Log:
debuggind pdd methode
Modified: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp 2012-06-04 21:16:54 UTC (rev 12361)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp 2012-06-04 21:33:55 UTC (rev 12362)
@@ -2254,9 +2254,8 @@
/*FUNCTION Penta::PositiveDegreeDay{{{1*/
void Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
-
int i,iqj,imonth;
- double agd[NUMVERTICES]; // surface and basal
+ double agd[NUMVERTICES]; // surface mass balance
double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
double smelt[NUMVERTICES] = {0}; // yearly melt
double precrunoff[NUMVERTICES]; // yearly runoff
@@ -2265,15 +2264,14 @@
double runoff; //meltwater only, does not include rain
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.
+ 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 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
double st; // elevation between altitude of the temp record and current altitude
double sp; // elevation between altitude of the prec record and current altitude
-
// PDD and PD constants and variables
double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
@@ -2293,10 +2291,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][12],prec[NUMVERTICES][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.
@@ -2307,25 +2305,31 @@
double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
double pddtj[NUMVERTICES], hmx2;
+ /*Recover monthly temperatures and precipitation*/
+ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
+ Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
+ GaussPenta* gauss=new GaussPenta();
+ 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
+ }
+ }
/*Recover info at the vertices: */
GetInputListOnVertices(&h[0],ThicknessEnum);
GetInputListOnVertices(&s[0],SurfaceEnum);
- GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
- GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
- for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion 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];
- }
-
/*Get material parameters :*/
rho_ice=matpar->GetRhoIce();
- //rho_freshwater=matpar->GetRhoWater();
+ rho_water=matpar->GetRhoFreshwater();
- sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
+ sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
/*PDD constant*/
siglim = 2.5*signorm;
@@ -2338,31 +2342,31 @@
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 ****************/
if (tstar < PDup){
- pd = 1;
+ pd = 1.;
if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
else {
- pd = 0;}
+ 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
- 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;
-
+
/*********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)];
@@ -2466,13 +2470,11 @@
Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
agd[i] = -smelt[i]+saccu[i];
+ agd[i] = agd[i]/yts;
pddtj[i]=pddt;
-
/*Update inputs*/
- this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
- // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
- this->inputs->AddInput(new PentaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
-
+ this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
+ // this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
} //end of the for loop over the vertices
}
/*}}}*/
Modified: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp 2012-06-04 21:16:54 UTC (rev 12361)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp 2012-06-04 21:33:55 UTC (rev 12362)
@@ -2087,7 +2087,6 @@
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
@@ -2136,8 +2135,9 @@
double supice,supcap,diffndd;
double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
double pddtj[NUMVERTICES], hmx2;
-int iii;
+ if(!IsOnBed()) return;
+
/*Recover monthly temperatures and precipitation*/
Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
@@ -2154,25 +2154,10 @@
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: */
+ /*Recover info at the vertices: */
GetInputListOnVertices(&h[0],ThicknessEnum);
GetInputListOnVertices(&s[0],SurfaceEnum);
- //GetInputListOnVertices(&monthlyprec[0],SurfaceforcingsPrecipitationEnum);
-
+
/*Get material parameters :*/
rho_ice=matpar->GetRhoIce();
rho_water=matpar->GetRhoFreshwater();
@@ -2181,7 +2166,7 @@
/*PDD constant*/
siglim = 2.5*signorm;
- siglimc = 2.5*signormc;
+ siglimc = 2.5*signormc;
siglim0 = siglim/DT + 0.5;
siglim0c = siglimc/DT + 0.5;
PDup = siglimc+PDCUT;
@@ -2196,10 +2181,10 @@
/*********compute PD ****************/
if (tstar < PDup){
- pd = 1;
+ pd = 1.;
if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
else {
- pd = 0;}
+ pd = 0.;}
/******exp des/elev precip reduction*******/
sp=(s[i]-s0p[i])/1000.; // deselev effect is wrt chng in topo
@@ -2210,22 +2195,7 @@
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)
@@ -2333,28 +2303,12 @@
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]));
// this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
- // 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]);
- // }
-
+ } //end of the for loop over the vertices
}
/*}}}*/
/*FUNCTION Tria::ProcessResultsUnits{{{1*/
Modified: issm/trunk/src/m/model/extrude.m
===================================================================
--- issm/trunk/src/m/model/extrude.m 2012-06-04 21:16:54 UTC (rev 12361)
+++ issm/trunk/src/m/model/extrude.m 2012-06-04 21:33:55 UTC (rev 12362)
@@ -142,6 +142,7 @@
md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node');
md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node');
md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
+md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node');
%results
if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end;
More information about the issm-svn
mailing list