Changeset 11994 for issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

Timestamp:

04/16/12 14:48:34 (13 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk and trunk-jpl

Location:

issm/trunk-jpl

Files:

• . (modified) (2 props)
• src/c/objects/Elements/Penta.cpp (modified) (1 diff)

issm/trunk-jpl ¶

@@ -7 +7 @@
 config.status
 configure
-doxygen
 ISSM.paf
 ISSM.ppf
 ISSM.ppf_cache
 libtool
-list
 Makefile
 Makefile.in
 stamp-h1
 svn-commit*
-nightlylog

@@ -7 +7 @@
 config.status
 configure
-doxygen
 ISSM.paf
 ISSM.ppf
 ISSM.ppf_cache
 libtool
-list
 Makefile
 Makefile.in
 stamp-h1
 svn-commit*
-nightlylog

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp ¶

@@ -2410 +2410 @@
 /*}}}*/
 /*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
 }
 /*}}}*/

@@ -7 +7 @@
 config.status
 configure
-doxygen
 ISSM.paf
 ISSM.ppf
 ISSM.ppf_cache
 libtool
-list
 Makefile
 Makefile.in
 stamp-h1
 svn-commit*
-nightlylog