Changeset 11778

Timestamp:

03/22/12 14:45:29 (13 years ago)

Author:

lemorzad

Message:

ton message

Location:

issm/trunk/src

Files:

• c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (modified) (1 diff)
• c/objects/Elements/Element.h (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (1 diff)
• c/objects/Elements/Penta.h (modified) (1 diff)
• c/objects/Elements/Tria.cpp (modified) (1 diff)
• c/objects/Elements/Tria.h (modified) (1 diff)
• m/model/parameterization/parameterize.m (modified) (1 diff)

issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp

@@ -11 +11 @@
 
 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);
+  
 }

issm/trunk/src/c/objects/Elements/Element.h

@@ -68 +68 @@
                 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;

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -2379 +2379 @@
 /*}}}*/
 /*FUNCTION Penta::PositiveDegreeDay{{{1*/
-void  Penta::PositiveDegreeDay(){
-
-        _error_("Not implemented yet");
+void  Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
+
+
+   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
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Penta.h

@@ -112 +112 @@
                 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);

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -2204 +2204 @@
 /*}}}*/
 /*FUNCTION Tria::PositiveDegreeDay{{{1*/
-void  Tria::PositiveDegreeDay(){
-
-        _error_("Not implemented yet");
+void  Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
+
+   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
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Tria.h

@@ -107 +107 @@
                 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);

issm/trunk/src/m/model/parameterization/parameterize.m

@@ -29 +29 @@
         eval(temporaryname);
         delete([temporaryname '.m']);
-catch me,
+catch,
         delete([temporaryname '.m']);