[issm-svn] r11778 - issm/trunk/src
lemorzad at issm.ess.uci.edu
lemorzad at issm.ess.uci.edu
Thu Mar 22 14:45:29 PDT 2012
Author: lemorzad
Date: 2012-03-22 14:45:29 -0700 (Thu, 22 Mar 2012)
New Revision: 11778
Modified:
issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
issm/trunk/src/c/objects/Elements/Element.h
issm/trunk/src/c/objects/Elements/Penta.cpp
issm/trunk/src/c/objects/Elements/Penta.h
issm/trunk/src/c/objects/Elements/Tria.cpp
issm/trunk/src/c/objects/Elements/Tria.h
issm/trunk/src/m/model/parameterization/parameterize.m
Log:
ton message
Modified: issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
===================================================================
--- issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-03-22 21:45:29 UTC (rev 11778)
@@ -10,12 +10,101 @@
#include "../../EnumDefinitions/EnumDefinitions.h"
void PositiveDegreeDayx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
-
- Element* element = NULL;
- for(int i=0;i<elements->Size();i++){
- element=(Element*)elements->GetObjectByOffset(i);
- element->PositiveDegreeDay();
- }
+// void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
+// note "v" prefix means 12 monthly means, ie time dimension
+// INPUT parameters: ni: working size of arrays
+// INPUT: surface elevation (m): hd(NA)
+// monthly mean surface sealevel temperature (degrees C): vTempsea(NA
+// ,NTIME)
+// monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
+// OUTPUT: mass-balance (m/yr ice): agd(NA)
+// mean annual surface temperature (degrees C): Tsurf(NA)
+ int i, it, jj, itm;
+ double DT = 0.02, sigfac, snormfac;
+ double signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
+ double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
+ double signormc; // sigma of the temperature distribution for cloudy day
+ double siglimc=0, siglim0, siglim0c;
+ double tstep, tsint, tint, tstepc;
+ int NPDMAX = 1504, NPDCMAX = 1454;
+ //double pdds[NPDMAX]={0};
+ //double pds[NPDCMAX]={0};
+ double pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
+ double PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
+ double tstar; // monthly mean surface temp
+
+ double* pdds=NULL;
+ double* pds=NULL;
+ Element* element = NULL;
+
+ pdds=(double*)xmalloc(NPDMAX*sizeof(double)+1);
+ pds=(double*)xmalloc(NPDCMAX*sizeof(double)+1);
+
+
+ // PDD constant
+ siglim = 2.5*signorm;
+
+ // initialize PDD (creation of a lookup table)
+ tstep = 0.1;
+ tsint = tstep*0.5;
+ sigfac = -1.0/(2.0*pow(signorm,2));
+ snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
+ siglim = 2.5*signorm;
+ itm = (int)(2*siglim/DT + 1.5);
+
+ if (itm >= NPDMAX){
+ printf("increase NPDMAX in massBalance.cpp\n");
+ exit (1);
+ }
+ for (it = 0; it < itm; it++){
+ // tstar = REAL(it)*DT-siglim;
+ tstar = it*DT-siglim;
+ tint = tsint;
+ pddt = 0;
+ for ( jj = 0; jj < 600; jj++){
+ if (tint > (tstar+siglim)){break;}
+ pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
+ tint = tint+tstep;
+ }
+ pdds[it] = pddt*snormfac;
+ }
+ pdds[itm+1] = siglim + DT;
+
+ //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
+ tstepc = 0.1;
+ tsint = PDCUT-tstepc*0.5;
+ signormc = signorm - 0.5;
+ sigfac = -1.0/(2.0*pow(signormc,2));
+ snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
+ siglimc = 2.5*signormc ;
+ itm = (int)((PDCUT+2.*siglimc)/DT + 1.5);
+ if (itm >= NPDCMAX){
+ printf("'increase NPDCMAX in p35com'\n");
+ exit (1);
+ }
+ for (it = 0; it < itm; it++ ){
+ tstar = it*DT-siglimc;
+ // tstar = REAL(it)*DT-siglimc;
+ tint = tsint; // start against upper bound
+ pd = 0;
+ for (jj = 0; jj < 600; jj++){
+ if (tint<(tstar-siglimc)) {break;}
+ pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
+ tint = tint-tstepc;
+ }
+ pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
+ }
+ pds[itm+1] = 0;
+ // *******END initialize PDD
+
+ for(i=0;i<elements->Size();i++){
+ element=(Element*)elements->GetObjectByOffset(i);
+ element->PositiveDegreeDay(pdds,pds,signorm);
+ }
+ /*free ressouces: */
+ xfree((void**)&pdds);
+ xfree((void**)&pds);
+
}
Modified: issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Element.h 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/objects/Elements/Element.h 2012-03-22 21:45:29 UTC (rev 11778)
@@ -67,7 +67,7 @@
virtual double TimeAdapt()=0;
virtual void MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding)=0;
virtual void PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0;
- virtual void PositiveDegreeDay(void)=0;
+ virtual void PositiveDegreeDay(double* pdds,double* pds,double signorm)=0;
virtual int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0;
virtual void ResetCoordinateSystem()=0;
virtual void SmearFunction(Vec smearedvector,double (*WeightFunction)(double distance,double radius),double radius)=0;
Modified: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp 2012-03-22 21:45:29 UTC (rev 11778)
@@ -2378,9 +2378,227 @@
}
/*}}}*/
/*FUNCTION Penta::PositiveDegreeDay{{{1*/
-void Penta::PositiveDegreeDay(){
+void Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
- _error_("Not implemented yet");
+
+ int i,iqj,imonth;
+ double agd[NUMVERTICES]; // surface and basal
+ 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.
+ 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 siglimc=0, 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],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 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;
+
+ /*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_water=matpar->GetRhoWater();
+ density=rho_ice/rho_water;
+
+ sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
+
+ /*PDD constant*/
+ siglim = 2.5*signorm;
+ 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;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];
+ 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
}
/*}}}*/
/*FUNCTION Penta::PotentialSheetUngrounding{{{1*/
Modified: issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.h 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/objects/Elements/Penta.h 2012-03-22 21:45:29 UTC (rev 11778)
@@ -111,7 +111,7 @@
void ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results);
void PatchFill(int* pcount, Patch* patch);
void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
- void PositiveDegreeDay(void);
+ void PositiveDegreeDay(double* pdds,double* pds,double signorm);
void ProcessResultsUnits(void);
void ResetCoordinateSystem(void);
double SurfaceArea(void);
Modified: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp 2012-03-22 21:45:29 UTC (rev 11778)
@@ -2203,9 +2203,226 @@
}
/*}}}*/
/*FUNCTION Tria::PositiveDegreeDay{{{1*/
-void Tria::PositiveDegreeDay(){
+void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
- _error_("Not implemented yet");
+ int i,iqj,imonth;
+ double agd[NUMVERTICES]; // surface and basal
+ 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.
+ 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 siglimc=0, 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],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 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;
+
+ /*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_water=matpar->GetRhoWater();
+ density=rho_ice/rho_water;
+
+ sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
+
+ /*PDD constant*/
+ siglim = 2.5*signorm;
+ 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;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];
+ 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
}
/*}}}*/
/*FUNCTION Tria::ProcessResultsUnits{{{1*/
Modified: issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.h 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/c/objects/Elements/Tria.h 2012-03-22 21:45:29 UTC (rev 11778)
@@ -106,7 +106,7 @@
void MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
void MigrateGroundingLine(double* oldfloating,double* sheet_ungrounding);
void PotentialSheetUngrounding(Vec potential_sheet_ungrounding);
- void PositiveDegreeDay(void);
+ void PositiveDegreeDay(double* pdds,double* pds,double signorm);
void RequestedOutput(int output_enum,int step,double time);
void ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results);
void PatchFill(int* pcount, Patch* patch);
Modified: issm/trunk/src/m/model/parameterization/parameterize.m
===================================================================
--- issm/trunk/src/m/model/parameterization/parameterize.m 2012-03-22 15:08:37 UTC (rev 11777)
+++ issm/trunk/src/m/model/parameterization/parameterize.m 2012-03-22 21:45:29 UTC (rev 11778)
@@ -28,7 +28,7 @@
try,
eval(temporaryname);
delete([temporaryname '.m']);
-catch me,
+catch,
delete([temporaryname '.m']);
%copy error message
More information about the issm-svn
mailing list