[issm-svn] r12640 - issm/trunk/src/c
lemorzad at issm.ess.uci.edu
lemorzad at issm.ess.uci.edu
Tue Jul 17 09:11:29 PDT 2012
Author: lemorzad
Date: 2012-07-17 09:11:29 -0700 (Tue, 17 Jul 2012)
New Revision: 12640
Added:
issm/trunk/src/c/shared/Elements/PddSurfaceMassBalance.cpp
Modified:
issm/trunk/src/c/Makefile.am
issm/trunk/src/c/objects/Elements/Penta.cpp
issm/trunk/src/c/objects/Elements/Tria.cpp
issm/trunk/src/c/shared/Elements/elements.h
Log:
adding PddSurfaceMassBalance. Numerical calculation of PositiveDegreeDay
Modified: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am 2012-07-16 22:23:55 UTC (rev 12639)
+++ issm/trunk/src/c/Makefile.am 2012-07-17 16:11:29 UTC (rev 12640)
@@ -204,6 +204,7 @@
./shared/Elements/GetLocalDofList.cpp\
./shared/Elements/GetGlobalDofList.cpp\
./shared/Elements/GetNumberOfDofs.cpp\
+ ./shared/Elements/PddSurfaceMassBalance.cpp\
./shared/String/sharedstring.h\
./shared/Wrapper/wrappershared.h\
./shared/Wrapper/ModuleBoot.cpp\
Modified: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp 2012-07-16 22:23:55 UTC (rev 12639)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp 2012-07-17 16:11:29 UTC (rev 12640)
@@ -2257,60 +2257,12 @@
/*FUNCTION Penta::PositiveDegreeDay{{{1*/
void Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
- int i,iqj,imonth;
- double agd[NUMVERTICES]; // surface mass balance
- double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
- double smelt[NUMVERTICES] = {0}; // yearly melt
- double precrunoff[NUMVERTICES]; // yearly runoff
- double prect; // total precipitation during 1 year taking into account des. ef.
- double water; //water=rain + snowmelt
- 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. 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
- double siglimc, siglim0, siglim0c;
- double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
- double DT = 0.02;
- double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
-
- double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
- double qm[NUMVERTICES] = {0}; // snow part of the precipitation
- double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment
- double qmp[NUMVERTICES] = {0}; // desertification taken into account
- double pdd[NUMVERTICES] = {0};
- double frzndd[NUMVERTICES] = {0};
-
- double tstar; // monthly mean surface temp
- double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature
- double Tsurf[NUMVERTICES] = {0}; // average annual temperature
-
- double h[NUMVERTICES],s[NUMVERTICES]; // ,b[NUMVERTICES]
+ double agd[NUMVERTICES]; // surface mass balance
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 h[NUMVERTICES],s[NUMVERTICES]; // ,b
+ double rho_water,rho_ice;
+ int i;
- double snwm; // snow that could have been melted in a year.
- double snwmf; // ablation factor for snow per positive degree day.
- double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
-
- double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
- double supice,supcap,diffndd;
- double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
- double pddtj[NUMVERTICES], hmx2;
-
- //printf("ok1p\n");
- if(!IsOnBed()) return;
-
/*Recover monthly temperatures and precipitation*/
Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
@@ -2327,165 +2279,24 @@
monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
}
}
- /*Recover info a the vertices: */
- GetInputListOnVertices(&h[0],ThicknessEnum);
- GetInputListOnVertices(&s[0],SurfaceEnum);
- //printf("ok2p\n");
+ /*Recover info at the vertices: */
+ GetInputListOnVertices(&h[0],ThicknessEnum);
+ GetInputListOnVertices(&s[0],SurfaceEnum);
- /*Get material parameters :*/
- rho_ice=matpar->GetRhoIce();
- rho_water=matpar->GetRhoFreshwater();
-
- sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
-
- /*PDD constant*/
- siglim = 2.5*signorm;
- siglimc = 2.5*signormc;
- siglim0 = siglim/DT + 0.5;
- siglim0c = siglimc/DT + 0.5;
- PDup = siglimc+PDCUT;
-
- // seasonal loop
- for (iqj = 0; iqj < 12; iqj++){
- imonth = ismon[iqj];
- for (i = 0; i < NUMVERTICES; i++){
- 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.;
- if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
- else {
- pd = 0.;}
-
- /******exp des/elev precip reduction*******/
- 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] + 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;
+ /*Get material parameters :*/
+ rho_ice=matpar->GetRhoIce();
+ rho_water=matpar->GetRhoFreshwater();
- /*********compute PDD************/
- // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
- // gaussian=T_m, so ndd=-(Tsurf-pdd)
- 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)];
- pdd[i] = pdd[i] + pddsig*deltm;
- frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
- else{frzndd[i] = frzndd[i] - tstar*deltm; }
- }
- } // end of seasonal loop
-
- //******************************************************************
- for(i=0;i<NUMVERTICES;i++){
- saccu[i] = qm[i];
- prect = qmp[i]; // total precipitation during 1 year taking into account des. ef.
- Tsum[i]=Tsum[i]/3;
-
- /***** determine PDD factors *****/
- 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;
- }
- else{
- snwmf=4.3*0.001;
- smf=8.3*0.001;
- }
- snwmf=0.95*snwmf;
- smf=0.95*smf;
-
- /***** compute PDD ablation and refreezing *****/
- pddt = pdd[i] *365;
- snwm = snwmf*pddt; // snow that could have been melted in a year
- hmx2 = min(h[i],dfrz); // refreeze active layer max depth: dfrz
-
- if(snwm < saccu[i]) {
- water=prect-saccu[i] + snwm; //water=rain + snowmelt
- // l 2.2= capillary factor
- // Should refreezing be controlled by frzndd or by mean annual Tsurf?
- // dCovrLm concept is of warming of active layer (thickness =d@=1-
- // >2m)
- // problem with water seepage into ice: should be sealed after
- // refreezing
- // so everything needs to be predicated on 1 year scale, except for
- // thermal
- // conductivity through ice
- // also, need to account that melt season has low accum, so what's
- // going to
- // hold the meltwater around for refreezing? And melt-time will have
- // low seasonal frzndd
-
- // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002
-
- supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
- supcap=min(2.2*(saccu[i]-snwm),water);
- runoff=snwm - supice; //meltwater only, does not include rain
- }
- else { //all snow melted
- supice= min(hmx2*CovrLm*frzndd[i], prect );
- runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
- supcap=0;
- }
- // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
- // except pdd melt heat source is atmosphere, while refreeze is
- // ground/ice stored interim
- // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
- // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
- // <0
- // assume ndd>pdd, little melt => little supice
- // bottom line: compare for Tsurf<0 : supice and no supice case,
- // expect Tsurf difference
- // except some of cooling flux comes from atmosphere//
- // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
- // does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
- // < 0.1
-
- // make more sense to just use residual pdd-ndd except that pdd
- // residual not clear yet
- // frzndd should not be used up by refreezing in snow, so stick in
- // supcap.
- diffndd=0;
- if (frzndd[i]>0) {
- diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
- frzndd[i]=frzndd[i]-diffndd;
- }
- if(runoff<0){
- saccu[i]= saccu[i] -runoff;
- smelt[i] = 0;
- precrunoff[i]=prect-saccu[i];
- //here assume pdd residual is 0, =>
- Tsurf[i]= max(Tsurf[i],-frzndd[i]);
- }
- else {
- smelt[i] = runoff;
- precrunoff[i]=prect-max(0.,supice)-saccu[i];}
- //here really need pdd balance, try 0.5 fudge factor?
- //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
- //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];
- agd[i] = agd[i]/yts;
- pddtj[i]=pddt;
- }//end of the for loop over the vertices
+ /*measure the surface mass balance*/
+ for (i = 0; i < NUMVERTICES; i++){
+ agd[i]=PddSurfaceMassBlance(&monthlytemperatures[0][0], &monthlyprec[0][0], i, pdds, pds, signorm, yts, h[i], s[i], rho_ice, rho_water);
+ }
- /*Update inputs*/
- this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
- //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
- this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
+ /*Update inputs*/
+ this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
+ //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
+ this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
}
/*}}}*/
/*FUNCTION Penta::PotentialSheetUngrounding{{{1*/
Modified: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp 2012-07-16 22:23:55 UTC (rev 12639)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp 2012-07-17 16:11:29 UTC (rev 12640)
@@ -2088,61 +2088,12 @@
/*FUNCTION Tria::PositiveDegreeDay{{{1*/
void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
- int i,iqj,imonth;
double agd[NUMVERTICES]; // surface mass balance
- double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
- double smelt[NUMVERTICES] = {0}; // yearly melt
- double precrunoff[NUMVERTICES]; // yearly runoff
- double prect; // total precipitation during 1 year taking into account des. ef.
- double water; //water=rain + snowmelt
- 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. 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
- double siglimc, siglim0, siglim0c;
- double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
- double DT = 0.02;
- double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
-
- double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
- double qm[NUMVERTICES] = {0}; // snow part of the precipitation
- double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment
- double qmp[NUMVERTICES] = {0}; // desertification taken into account
- double pdd[NUMVERTICES] = {0};
- double frzndd[NUMVERTICES] = {0};
-
- double tstar; // monthly mean surface temp
- double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature
- double Tsurf[NUMVERTICES] = {0}; // average annual temperature
-
- 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 h[NUMVERTICES],s[NUMVERTICES]; // ,b
+ double rho_water,rho_ice;
+ int i;
- double snwm; // snow that could have been melted in a year.
- double snwmf; // ablation factor for snow per positive degree day.
- double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
-
- double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
- double supice,supcap,diffndd;
- double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
- double pddtj[NUMVERTICES], hmx2;
-
- //printf("ok1t\n");
- //if(!IsOnBed()) return;
- //printf("ok2t\n");
-
/*Recover monthly temperatures and precipitation*/
Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
@@ -2159,163 +2110,23 @@
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(&h[0],ThicknessEnum);
+ GetInputListOnVertices(&s[0],SurfaceEnum);
- /*Get material parameters :*/
- rho_ice=matpar->GetRhoIce();
- rho_water=matpar->GetRhoFreshwater();
-
- sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
-
- /*PDD constant*/
- siglim = 2.5*signorm;
- siglimc = 2.5*signormc;
- siglim0 = siglim/DT + 0.5;
- siglim0c = siglimc/DT + 0.5;
- PDup = siglimc+PDCUT;
-
- // seasonal loop
- for (iqj = 0; iqj < 12; iqj++){
- imonth = ismon[iqj];
- for (i = 0; i < NUMVERTICES; i++){
- 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.;
- if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
- else {
- pd = 0.;}
-
- /******exp des/elev precip reduction*******/
- 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] + 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;
+ /*Get material parameters :*/
+ rho_ice=matpar->GetRhoIce();
+ rho_water=matpar->GetRhoFreshwater();
- /*********compute PDD************/
- // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
- // gaussian=T_m, so ndd=-(Tsurf-pdd)
- 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)];
- pdd[i] = pdd[i] + pddsig*deltm;
- frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
- else{frzndd[i] = frzndd[i] - tstar*deltm; }
- }
- } // end of seasonal loop
-
- //******************************************************************
- for(i=0;i<NUMVERTICES;i++){
- saccu[i] = qm[i];
- prect = qmp[i]; // total precipitation during 1 year taking into account des. ef.
- Tsum[i]=Tsum[i]/3;
-
- /***** determine PDD factors *****/
- 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;
- }
- else{
- snwmf=4.3*0.001;
- smf=8.3*0.001;
- }
- snwmf=0.95*snwmf;
- smf=0.95*smf;
-
- /***** compute PDD ablation and refreezing *****/
- pddt = pdd[i] *365;
- snwm = snwmf*pddt; // snow that could have been melted in a year
- hmx2 = min(h[i],dfrz); // refreeze active layer max depth: dfrz
-
- if(snwm < saccu[i]) {
- water=prect-saccu[i] + snwm; //water=rain + snowmelt
- // l 2.2= capillary factor
- // Should refreezing be controlled by frzndd or by mean annual Tsurf?
- // dCovrLm concept is of warming of active layer (thickness =d@=1-
- // >2m)
- // problem with water seepage into ice: should be sealed after
- // refreezing
- // so everything needs to be predicated on 1 year scale, except for
- // thermal
- // conductivity through ice
- // also, need to account that melt season has low accum, so what's
- // going to
- // hold the meltwater around for refreezing? And melt-time will have
- // low seasonal frzndd
-
- // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002
-
- supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
- supcap=min(2.2*(saccu[i]-snwm),water);
- runoff=snwm - supice; //meltwater only, does not include rain
- }
- else { //all snow melted
- supice= min(hmx2*CovrLm*frzndd[i], prect );
- runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
- supcap=0;
- }
- // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
- // except pdd melt heat source is atmosphere, while refreeze is
- // ground/ice stored interim
- // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
- // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
- // <0
- // assume ndd>pdd, little melt => little supice
- // bottom line: compare for Tsurf<0 : supice and no supice case,
- // expect Tsurf difference
- // except some of cooling flux comes from atmosphere//
- // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
- // does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
- // < 0.1
-
- // make more sense to just use residual pdd-ndd except that pdd
- // residual not clear yet
- // frzndd should not be used up by refreezing in snow, so stick in
- // supcap.
- diffndd=0;
- if (frzndd[i]>0) {
- diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
- frzndd[i]=frzndd[i]-diffndd;
- }
- if(runoff<0){
- saccu[i]= saccu[i] -runoff;
- smelt[i] = 0;
- precrunoff[i]=prect-saccu[i];
- //here assume pdd residual is 0, =>
- Tsurf[i]= max(Tsurf[i],-frzndd[i]);
- }
- else {
- smelt[i] = runoff;
- precrunoff[i]=prect-max(0.,supice)-saccu[i];}
- //here really need pdd balance, try 0.5 fudge factor?
- //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
- //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];
- agd[i] = agd[i]/yts;
- pddtj[i]=pddt;
+ /*measure the surface mass balance*/
+ for (i = 0; i < NUMVERTICES; i++){
+ agd[i]=PddSurfaceMassBlance(&monthlytemperatures[0][0], &monthlyprec[0][0], i, pdds, pds, signorm, yts, h[i], s[i], rho_ice, rho_water);
+ }
- } //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]));
+ /*Update inputs*/
+ this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
+ // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
}
/*}}}*/
/*FUNCTION Tria::ProcessResultsUnits{{{1*/
Modified: issm/trunk/src/c/shared/Elements/elements.h
===================================================================
--- issm/trunk/src/c/shared/Elements/elements.h 2012-07-16 22:23:55 UTC (rev 12639)
+++ issm/trunk/src/c/shared/Elements/elements.h 2012-07-17 16:11:29 UTC (rev 12640)
@@ -12,6 +12,7 @@
double Paterson(double temperature);
double Arrhenius(double temperature,double depth,double n);
+double PddSurfaceMassBlance(double* monthlytemperatures, double* monthlyprec, int i, double* pdds, double* pds, double signorm, double yts, double h, double s, double rho_ice, double rho_water);
void GetVerticesCoordinates(double* xyz, Node** nodes, int numvertices);
int GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum);
int* GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum);
More information about the issm-svn
mailing list