Changeset 12640 for issm/trunk/src/c/objects/Elements/Penta.cpp

Timestamp:

07/17/12 09:11:29 (13 years ago)

Author:

lemorzad

Message:

adding PddSurfaceMassBalance. Numerical calculation of PositiveDegreeDay

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (2 diffs)

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

@@ -2258 +2258 @@
 void  Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
 
-   int    i,iqj,imonth;
-   double agd[NUMVERTICES];  // surface mass balance
-   double saccu[NUMVERTICES] = {0};     // yearly surface accumulation
-   double smelt[NUMVERTICES] = {0};     // yearly melt
-   double precrunoff[NUMVERTICES];      // yearly runoff
-   double prect; // total precipitation during 1 year taking into account des. ef.
-   double water; //water=rain + snowmelt 
-   double runoff; //meltwater only, does not include rain 
-   double sconv; //rhow_rain/rhoi / 12 months
-
-   double rho_water,rho_ice,density;
-   double lapser=6.5/1000., sealev=0.;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
-   double desfac = 0.5;                 //desert elevation factor
-   double s0p[NUMVERTICES]={0};         //should be set to elevation from precip source
-   double s0t[NUMVERTICES]={0};         //should be set to elevation from temperature source
-   double st;             // elevation between altitude of the temp record and current altitude
-   double sp;             // elevation between altitude of the prec record and current altitude
-
-   // PDD and PD constants and variables
-   double siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
-   double signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
-   double siglimc, siglim0, siglim0c;
-   double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
-   double DT = 0.02;
-   double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
-   
-   double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
-   double qm[NUMVERTICES] = {0};        // snow part of the precipitation 
-   double qmt[NUMVERTICES] = {0};       // precipitation without desertification effect adjustment
-   double qmp[NUMVERTICES] = {0};       // desertification taken into account
-   double pdd[NUMVERTICES] = {0};     
-   double frzndd[NUMVERTICES] = {0};  
-
-   double tstar;                        // monthly mean surface temp
-   double Tsum[NUMVERTICES]= {0};       // average summer (JJA) temperature
-   double Tsurf[NUMVERTICES] = {0};     // average annual temperature    
-   
-   double h[NUMVERTICES],s[NUMVERTICES]; // ,b[NUMVERTICES]
+   double agd[NUMVERTICES];             // surface mass balance
    double monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
-   double deltm=1./12.;
-   int    ismon[12]={11,0,1,2,3,4,5,6,7,8,9,10};
-
-   double snwm;  // snow that could have been melted in a year.
-   double snwmf; //  ablation factor for snow per positive degree day.
-   double smf;   //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
-
-   double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
-   double supice,supcap,diffndd;
-   double fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
-   double pddtj[NUMVERTICES], hmx2;
-
-   //printf("ok1p\n");
-   if(!IsOnBed()) return;
+   double h[NUMVERTICES],s[NUMVERTICES]; // ,b
+   double rho_water,rho_ice;
+   int    i;
 
    /*Recover monthly temperatures and precipitation*/
@@ -2328 +2280 @@
      }
    } 
-   /*Recover info a the vertices: */
-   GetInputListOnVertices(&h[0],ThicknessEnum);
-   GetInputListOnVertices(&s[0],SurfaceEnum);
-
-   //printf("ok2p\n");
-
-   /*Get material parameters :*/
-   rho_ice=matpar->GetRhoIce();
-   rho_water=matpar->GetRhoFreshwater();
-   
-   sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
-     
-     /*PDD constant*/
-   siglim = 2.5*signorm; 
-   siglimc = 2.5*signormc;
-   siglim0 =  siglim/DT + 0.5;
-   siglim0c =  siglimc/DT + 0.5;
-   PDup = siglimc+PDCUT;
-   
-   // seasonal loop
-   for (iqj = 0; iqj < 12; iqj++){
-     imonth =  ismon[iqj];
-     for (i = 0; i < NUMVERTICES; i++){
-       st=(s[i]-s0t[i])/1000.;
-       tstar = monthlytemperatures[i][imonth] - lapser *max(st,sealev);
-       Tsurf[i] = tstar*deltm+Tsurf[i];        
-       
-       /*********compute PD ****************/
-       if (tstar < PDup){
-         pd = 1.;
-         if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
-       else { 
-         pd = 0.;}
-       
-       /******exp des/elev precip reduction*******/
-       sp=(s[i]-s0p[i])/1000.; // deselev effect is wrt chng in topo
-       if (sp>0.0){q = exp(-desfac*sp);}
-       else {q = 1.0;}
-       
-       qmt[i]= qmt[i] + monthlyprec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent per month
-       qmpt= q*monthlyprec[i][imonth]*sconv;           
-       qmp[i]= qmp[i] + qmpt;
-       qm[i]= qm[i] + qmpt*pd;
-
-       /*********compute PDD************/
-       // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
-       // gaussian=T_m, so ndd=-(Tsurf-pdd)
-       if (iqj>5 &&  iqj<9){ Tsum[i]=Tsum[i]+tstar;} 
-       if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
-       else if (tstar> -siglim){
-         pddsig=pdds[int(tstar/DT + siglim0)];
-         pdd[i] = pdd[i] + pddsig*deltm;
-         frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
-       else{frzndd[i] = frzndd[i] - tstar*deltm; }
-     }
-   } // end of seasonal loop 
-   
-   //******************************************************************
-   for(i=0;i<NUMVERTICES;i++){
-     saccu[i] = qm[i];
-     prect = qmp[i];     // total precipitation during 1 year taking into account des. ef.
-     Tsum[i]=Tsum[i]/3;
-     
-     /***** determine PDD factors *****/
-     if(Tsum[i]< -1.) {
-       snwmf=2.65*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
-       smf=17.22*0.001;    //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
-     }
-     else if(Tsum[i]< 10){
-       snwmf = (0.15*Tsum[i] + 2.8)*0.001;
-       smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
-     }
-     else{
-       snwmf=4.3*0.001;
-       smf=8.3*0.001;
-     }
-     snwmf=0.95*snwmf;
-     smf=0.95*smf;
-     
-     /*****  compute PDD ablation and refreezing *****/
-     pddt = pdd[i] *365;
-     snwm = snwmf*pddt;       // snow that could have been melted in a year
-     hmx2 = min(h[i],dfrz);   // refreeze active layer max depth: dfrz
-     
-     if(snwm < saccu[i]) {
-       water=prect-saccu[i] + snwm; //water=rain + snowmelt
-       //     l 2.2= capillary factor
-       //     Should refreezing be controlled by frzndd or by mean annual Tsurf?
-       //     dCovrLm concept is of warming of active layer (thickness =d@=1-
-       //     >2m)
-       //     problem with water seepage into ice: should be sealed after
-       //     refreezing
-       //     so everything needs to be predicated on 1 year scale, except for
-       //     thermal
-       //     conductivity through ice
-       //     also, need to account that melt season has low accum, so what's
-       //     going to
-       //     hold the meltwater around for refreezing? And melt-time will have
-       //     low seasonal frzndd
-       
-       //      Superimposed ice :  Pfeffer et al. 1991, Tarasov 2002
-       
-       supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
-       supcap=min(2.2*(saccu[i]-snwm),water);
-       runoff=snwm - supice;  //meltwater only, does not include rain
-     }
-     else {  //all snow melted
-       supice= min(hmx2*CovrLm*frzndd[i], prect );
-       runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
-       supcap=0;
-     }
-     //     pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
-     //     except pdd melt heat source is atmosphere, while refreeze is
-     //     ground/ice stored interim
-     //     assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
-     //     assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
-     //     <0
-     //     assume ndd>pdd, little melt => little supice 
-     //     bottom line: compare for Tsurf<0 : supice and no supice case,
-     //     expect Tsurf difference
-     //     except some of cooling flux comes from atmosphere//
-     //     1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
-     //     does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
-     //     < 0.1 
-     
-     //     make more sense to just use residual pdd-ndd except that pdd
-     //     residual not clear yet
-     //     frzndd should not be used up by refreezing in snow, so stick in
-     //     supcap.
-     diffndd=0;
-     if (frzndd[i]>0) {
-       diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
-       frzndd[i]=frzndd[i]-diffndd;
-     }
-     if(runoff<0){
-       saccu[i]= saccu[i] -runoff;
-       smelt[i] = 0;
-       precrunoff[i]=prect-saccu[i];
-       //here assume pdd residual is 0, => 
-       Tsurf[i]= max(Tsurf[i],-frzndd[i]);
-     }
-     else {
-       smelt[i] = runoff;
-       precrunoff[i]=prect-max(0.,supice)-saccu[i];}
-     //here really need pdd balance, try 0.5 fudge factor?
-     //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
-     //yet from site plots, can be ice free with Tsurf=-5.5C
-     if(Tsurf[i]<0) {
-       Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
-     
-     agd[i] = -smelt[i]+saccu[i];
-     agd[i] = agd[i]/yts;
-     pddtj[i]=pddt;
-   }//end of the for loop over the vertices
-
-        /*Update inputs*/    
-        this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
-        //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
-        this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
+
+  /*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 (i = 0; i < NUMVERTICES; i++){
+     agd[i]=PddSurfaceMassBlance(&monthlytemperatures[0][0], &monthlyprec[0][0], i, pdds, pds, signorm, yts, h[i], s[i], rho_ice, rho_water);
+   }
+
+   /*Update inputs*/    
+   this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
+   //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
+   this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
 }
 /*}}}*/