Index: /issm/trunk-jpl/externalpackages/vim/addons/vim/syntax/c.vim =================================================================== --- /issm/trunk-jpl/externalpackages/vim/addons/vim/syntax/c.vim (revision 12703) +++ /issm/trunk-jpl/externalpackages/vim/addons/vim/syntax/c.vim (revision 12704) @@ -706,14 +706,4 @@ syn keyword cConstant SurfaceforcingsPrecipitationEnum syn keyword cConstant SurfaceforcingsMassBalanceEnum -syn keyword cConstant SurfaceforcingsIspddEnum -syn keyword cConstant SurfaceforcingsIssmbgradientsEnum -syn keyword cConstant SurfaceforcingsMonthlytemperaturesEnum -syn keyword cConstant SurfaceforcingsHcEnum -syn keyword cConstant SurfaceforcingsSmbPosMaxEnum -syn keyword cConstant SurfaceforcingsSmbPosMinEnum -syn keyword cConstant SurfaceforcingsAPosEnum -syn keyword cConstant SurfaceforcingsBPosEnum -syn keyword cConstant SurfaceforcingsANegEnum -syn keyword cConstant SurfaceforcingsBNegEnum syn keyword cConstant ThermalMaxiterEnum syn keyword cConstant ThermalPenaltyFactorEnum Index: /issm/trunk-jpl/src/c/Container/Elements.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Elements.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/Container/Elements.cpp (revision 12704) @@ -237,5 +237,5 @@ for(int j=0;jSize();j++){ Element* element=(Element*)this->GetObjectByOffset(j); - element->GetVectorFromResults(vector,i,resultssizes[i]); + element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]); } vector->Assemble(); Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12703) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12704) @@ -205,4 +205,5 @@ ./shared/Elements/GetGlobalDofList.cpp\ ./shared/Elements/GetNumberOfDofs.cpp\ + ./shared/Elements/PddSurfaceMassBalance.cpp\ ./shared/String/sharedstring.h\ ./shared/Wrapper/wrappershared.h\ Index: /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp (revision 12704) @@ -8,10 +8,10 @@ #endif +#include +#include #include "../../shared/shared.h" #include "../../io/io.h" #include "../../include/include.h" #include "./matlabio.h" - -#include /*FUNCTION OptionDoubleParse {{{*/ Index: /issm/trunk-jpl/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (revision 12704) @@ -64,5 +64,5 @@ tstar = it*DT-siglim; tint = tsint; - pddt = 0; + pddt = 0.; for ( jj = 0; jj < 600; jj++){ if (tint > (tstar+siglim)){break;} @@ -90,5 +90,5 @@ // 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;} @@ -98,5 +98,5 @@ pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction } - pds[itm+1] = 0; + pds[itm+1] = 0.; // *******END initialize PDD Index: /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.cpp (revision 12704) @@ -126,2 +126,12 @@ } /*}}}*/ +/*FUNCTION DoubleElementResult::GetVectorFromResults{{{1*/ +void DoubleElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){ + + _error_("cannot return vector on vertices"); +} /*}}}*/ +/*FUNCTION DoubleElementResult::GetElementVectorFromResults{{{1*/ +void DoubleElementResult::GetElementVectorFromResults(Vector* vector,int dof){ + + vector->SetValue(dof,value,INS_VAL); +} /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.h (revision 12704) @@ -48,6 +48,6 @@ /*DoubleElementResult management: {{{*/ int InstanceEnum(); - void GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){_error2_("not implemented");}; - void GetElementVectorFromResults(Vector* vector,int dof){_error2_("not implemented");}; + void GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs); + void GetElementVectorFromResults(Vector* vector,int dof); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12704) @@ -63,5 +63,5 @@ virtual void InputScale(int enum_type,IssmDouble scale_factor)=0; virtual void GetVectorFromInputs(Vector* vector, int name_enum)=0; - virtual void GetVectorFromResults(Vector* vector,int id,int interp)=0; + virtual void GetVectorFromResults(Vector* vector,int id,int enum_in,int interp)=0; virtual void InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max)=0; virtual bool InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums)=0; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12704) @@ -1080,8 +1080,11 @@ /*}}}*/ /*FUNCTION Penta::GetVectorFromResults{{{*/ -void Penta::GetVectorFromResults(Vector* vector,int offset,int interp){ +void Penta::GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp){ /*Get result*/ ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset); + if(elementresult->InstanceEnum()!=enum_in){ + _error_("Results of offset %i is %s, when %s was expected",offset,EnumToStringx(elementresult->InstanceEnum()),EnumToStringx(enum_in)); + } if(interp==P1Enum){ int doflist1[NUMVERTICES]; @@ -2255,224 +2258,43 @@ void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){ - - int i,iqj,imonth; - IssmDouble agd[NUMVERTICES]; // surface and basal - IssmDouble saccu[NUMVERTICES] = {0}; // yearly surface accumulation - IssmDouble smelt[NUMVERTICES] = {0}; // yearly melt - IssmDouble precrunoff[NUMVERTICES]; // yearly runoff - IssmDouble prect; // total precipitation during 1 year taking into account des. ef. - IssmDouble water; //water=rain + snowmelt - IssmDouble runoff; //meltwater only, does not include rain - IssmDouble sconv; //rhow_rain/rhoi / 12 months - - IssmDouble rho_water,rho_ice,density; - IssmDouble lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. - IssmDouble desfac = 0.5; //desert elevation factor - IssmDouble s0p[NUMVERTICES]={0}; //should be set to elevation from precip source - IssmDouble s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source - IssmDouble st; // elevation between altitude of the temp record and current altitude - IssmDouble sp; // elevation between altitude of the prec record and current altitude - - - // PDD and PD constants and variables - IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day - IssmDouble siglimc, siglim0, siglim0c; - IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) - IssmDouble DT = 0.02; - IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - - IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. - IssmDouble qm[NUMVERTICES] = {0}; // snow part of the precipitation - IssmDouble qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment - IssmDouble qmp[NUMVERTICES] = {0}; // desertification taken into account - IssmDouble pdd[NUMVERTICES] = {0}; - IssmDouble frzndd[NUMVERTICES] = {0}; - - IssmDouble tstar; // monthly mean surface temp - IssmDouble Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - IssmDouble Tsurf[NUMVERTICES] = {0}; // average annual temperature - - IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] - IssmDouble t[NUMVERTICES][12],prec[NUMVERTICES][12]; - IssmDouble deltm=1/12; - int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11}; - - IssmDouble snwm; // snow that could have been melted in a year. - IssmDouble snwmf; // ablation factor for snow per positive degree day. - IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). - - IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr - IssmDouble supice,supcap,diffndd; - IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice - IssmDouble pddtj[NUMVERTICES], hmx2; - - /*Recover info at the vertices: */ - GetInputListOnVertices(&h[0],ThicknessEnum); - GetInputListOnVertices(&s[0],SurfaceEnum); - GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum); - GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum); - - for(i=0;iGetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); + GaussPenta* gauss=new GaussPenta(); + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(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); + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + + /*measure the surface mass balance*/ + for (int iv = 0; iv < NUMVERTICES; iv++){ + agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water); } - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - //rho_freshwater=matpar->GetRhoWater(); - - sconv=(1000/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 = t[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] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent - qmpt= q*prec[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 (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;i2m) - // 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]; - 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])); - - } //end of the for loop over the vertices + /*Update inputs*/ + this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); + //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); + this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum); } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12704) @@ -89,5 +89,5 @@ IssmDouble GetZcoord(GaussPenta* gauss); void GetVectorFromInputs(Vector* vector,int name_enum); - void GetVectorFromResults(Vector* vector,int offset,int interp); + void GetVectorFromResults(Vector* vector,int offset,int name_enum,int interp); int Sid(); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12704) @@ -1192,8 +1192,11 @@ /*}}}*/ /*FUNCTION Tria::GetVectorFromResults{{{*/ -void Tria::GetVectorFromResults(Vector* vector,int offset,int interp){ +void Tria::GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp){ /*Get result*/ ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset); + if(elementresult->InstanceEnum()!=enum_in){ + _error_("Results of offset %i is %s, when %s was expected",offset,EnumToStringx(elementresult->InstanceEnum()),EnumToStringx(enum_in)); + } if(interp==P1Enum){ int doflist1[NUMVERTICES]; @@ -2086,234 +2089,42 @@ void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){ - int i,iqj,imonth; IssmDouble agd[NUMVERTICES]; // surface mass balance - IssmDouble saccu[NUMVERTICES] = {0}; // yearly surface accumulation - IssmDouble smelt[NUMVERTICES] = {0}; // yearly melt - IssmDouble precrunoff[NUMVERTICES]; // yearly runoff - IssmDouble prect; // total precipitation during 1 year taking into account des. ef. - IssmDouble water; //water=rain + snowmelt - IssmDouble runoff; //meltwater only, does not include rain - IssmDouble sconv; //rhow_rain/rhoi / 12 months - - IssmDouble rho_water,rho_ice,density; - IssmDouble lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper - IssmDouble desfac = 0.5; // desert elevation factor - IssmDouble s0p[NUMVERTICES]={0}; // should be set to elevation from precip source - IssmDouble s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source - IssmDouble st; // elevation between altitude of the temp record and current altitude - IssmDouble sp; // elevation between altitude of the prec record and current altitude - - // PDD and PD constants and variables - IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm - IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day - IssmDouble siglimc, siglim0, siglim0c; - IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) - IssmDouble DT = 0.02; - IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - - IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. - IssmDouble qm[NUMVERTICES] = {0}; // snow part of the precipitation - IssmDouble qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment - IssmDouble qmp[NUMVERTICES] = {0}; // desertification taken into account - IssmDouble pdd[NUMVERTICES] = {0}; - IssmDouble frzndd[NUMVERTICES] = {0}; - - IssmDouble tstar; // monthly mean surface temp - IssmDouble Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature - IssmDouble Tsurf[NUMVERTICES] = {0}; // average annual temperature - - IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES] - IssmDouble t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12]; - IssmDouble deltm=1/12; - int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11}; - - IssmDouble snwm; // snow that could have been melted in a year. - IssmDouble snwmf; // ablation factor for snow per positive degree day. - IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). - - IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr - IssmDouble supice,supcap,diffndd; - IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice - IssmDouble pddtj[NUMVERTICES], hmx2; - - /*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(); - IssmDouble time,yts; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - for(int month=0;month<12;month++){ - for(int iv=0;ivGaussVertex(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;iGetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); + GaussTria* gauss=new GaussTria(); + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(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 + } } - /*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 = t[i][imonth] - lapser *max(st,sealev); - Tsurf[i] = tstar*deltm+Tsurf[i]; - - /*********compute PD ****************/ - if (tstar < PDup){ - pd = 1; - if (tstar >= -siglimc){ pd = pds[reCast(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] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent per month - qmpt= q*prec[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 (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;} - else if (tstar> -siglim){ - pddsig=pdds[reCast(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;i2m) - // 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]; - pddtj[i]=pddt; - - /*Update inputs*/ - this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom - // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); - this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0])); - - } //end of the for loop over the vertices + /*Recover info at the vertices: */ + GetInputListOnVertices(&h[0],ThicknessEnum); + GetInputListOnVertices(&s[0],SurfaceEnum); + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + + /*measure the surface mass balance*/ + for (int iv = 0; ivinputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); + // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12704) @@ -89,5 +89,5 @@ void GetSolutionFromInputs(Vector* solution); void GetVectorFromInputs(Vector* vector, int name_enum); - void GetVectorFromResults(Vector* vector,int offset,int interp); + void GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp); void InputArtificialNoise(int enum_type,IssmDouble min, IssmDouble max); bool InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums); Index: /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.h (revision 12704) @@ -49,6 +49,7 @@ void GetInputValue(IssmDouble* pvalue); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h (revision 12704) @@ -54,6 +54,7 @@ void GetInputValue(IssmDouble* pvalue); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/DatasetInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/DatasetInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/DatasetInput.h (revision 12704) @@ -44,11 +44,12 @@ void Configure(Parameters* parameters); /*}}}*/ - /*numeriics: {{{*/ + /*numerics: {{{*/ void GetInputValue(bool* pvalue){_error2_("not implemented yet");}; void GetInputValue(int* pvalue){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error2_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index); void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.h (revision 12704) @@ -48,6 +48,7 @@ void GetInputValue(IssmDouble* pvalue); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/Input.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/Input.h (revision 12704) @@ -21,5 +21,5 @@ virtual ~Input(){}; - /*Virtual functions:{{{*/ + virtual int InstanceEnum()=0; virtual void GetInputValue(bool* pvalue)=0; @@ -27,6 +27,7 @@ virtual void GetInputValue(IssmDouble* pvalue)=0; virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss)=0; + virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0; virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time)=0; - virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0; + virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time)=0; virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index)=0; virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index)=0; @@ -65,7 +66,4 @@ virtual Input* PointwiseMin(Input* inputmin)=0; virtual ElementResult* SpawnResult(int step, IssmDouble time)=0; - - /*}}}*/ - }; #endif Index: /issm/trunk-jpl/src/c/objects/Inputs/IntInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/IntInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/IntInput.h (revision 12704) @@ -49,6 +49,7 @@ void GetInputValue(IssmDouble* pvalue); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.h (revision 12704) @@ -49,6 +49,7 @@ void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error2_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp (revision 12704) @@ -172,6 +172,34 @@ } /*}}}*/ +/*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/ +void TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){ + IssmDouble time; + + /*First, recover current time from parameters: */ + this->parameters->FindParam(&time,TimeEnum); + + /*Retrieve interpolated values for this time step: */ + Input* input=GetTimeInput(time); + + /*Call input function*/ + input->GetInputValue(pvalue,gauss); + + delete input; +} +/*}}}*/ /*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){{{*/ void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){ + + /*Retrieve interpolated values for this time step: */ + Input* input=GetTimeInput(time); + + /*Call input function*/ + input->GetInputValue(pvalue,gauss); + + delete input; +} +/*}}}*/ +/*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){{{*/ +void TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){ /*Retrieve interpolated values for this time step: */ Index: /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/TransientInput.h (revision 12704) @@ -51,6 +51,7 @@ void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time); - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.h (revision 12703) +++ /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.h (revision 12704) @@ -49,6 +49,7 @@ void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");} void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){_error2_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 12704) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 12704) @@ -0,0 +1,205 @@ +/* file: PddSurfaceMassBlance.cpp + Calculating the surface mass balance using the positive degree day method. + */ + +#include "./elements.h" + +double PddSurfaceMassBlance(double* monthlytemperatures, double* monthlyprec, double* pdds, double* pds, double signorm, double yts, double h, double s, double rho_ice, double rho_water){ + + // output: + double B; // surface mass balance, melt+accumulation + + int iqj,imonth, j; + + double saccu; // yearly surface accumulation + double smelt; // yearly melt + double precrunoff; // 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 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=0.; // should be set to elevation from precip source + double s0t=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 = 0.; // snow part of the precipitation + double qmt = 0.; // precipitation without desertification effect adjustment + double qmp = 0.; // desertification taken into account + double pdd = 0.; + double frzndd = 0.; + + double tstar; // monthly mean surface temp + double Tsum= 0.; // average summer (JJA) temperature + double Tsurf = 0.; // average annual temperature + + + 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. + 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, hmx2; + + 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]; + + st=(s-s0t)/1000.; + tstar = monthlytemperatures[imonth] - lapser *max(st,sealev); + Tsurf = tstar*deltm+Tsurf; + + /*********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-s0p)/1000.; // deselev effect is wrt chng in topo + if (sp>0.0){q = exp(-desfac*sp);} + else {q = 1.0;} + + qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month + qmpt= q*monthlyprec[imonth]*sconv; + qmp= qmp + qmpt; + qm= qm + qmpt*pd; + + /*********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=Tsum+tstar;} + if (tstar >= siglim) {pdd = pdd + tstar*deltm;} + else if (tstar> -siglim){ + pddsig=pdds[int(tstar/DT + siglim0)]; + pdd = pdd + pddsig*deltm; + frzndd = frzndd - (tstar-pddsig)*deltm;} + else{frzndd = frzndd - tstar*deltm; } + } // end of seasonal loop + + //****************************************************************** + saccu = qm; + prect = qmp; // total precipitation during 1 year taking into account des. ef. + Tsum=Tsum/3; + + /***** determine PDD factors *****/ + if(Tsum< -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< 10){ + snwmf = (0.15*Tsum + 2.8)*0.001; + smf = (0.0067*pow((10.-Tsum),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 *365; + snwm = snwmf*pddt; // snow that could have been melted in a year + hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz + + if(snwm < saccu) { + water=prect-saccu + 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+2.2*(saccu-snwm), water); // superimposed ice + supcap=min(2.2*(saccu-snwm),water); + runoff=snwm - supice; //meltwater only, does not include rain + } + else { //all snow melted + supice= min(hmx2*CovrLm*frzndd, prect ); + runoff= saccu + smf*(pddt-saccu/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>0) { + diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd); + frzndd=frzndd-diffndd; + } + if(runoff<0){ + saccu= saccu -runoff; + smelt = 0; + precrunoff=prect-saccu; + //here assume pdd residual is 0, => + Tsurf= max(Tsurf,-frzndd); + } + else { + smelt = runoff; + precrunoff=prect-max(0.,supice)-saccu;} + //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<0) { + Tsurf= min(Tsurf+fsupT*diffndd , 0.);} + + B = -smelt+saccu; + B = B/yts; + pddtj=pddt; + + return B; +} Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12703) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12704) @@ -13,4 +13,5 @@ IssmDouble Paterson(IssmDouble temperature); IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); +IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water); void GetVerticesCoordinates(IssmDouble* xyz, Node** nodes, int numvertices); int GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum); Index: /issm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp (revision 12703) +++ /issm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp (revision 12704) @@ -7,4 +7,5 @@ #include +#include #include "../Alloc/alloc.h" #include "../../include/include.h" Index: /issm/trunk-jpl/src/m/model/extrude.m =================================================================== --- /issm/trunk-jpl/src/m/model/extrude.m (revision 12703) +++ /issm/trunk-jpl/src/m/model/extrude.m (revision 12704) @@ -143,4 +143,5 @@ 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