Changeset 12346

Timestamp:

06/02/12 08:27:24 (13 years ago)

Author:

lemorzad

Message:

debug pdd methode

Location:

issm/trunk/src/c

Files:

• modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (modified) (3 diffs)
• objects/Elements/Tria.cpp (modified) (9 diffs)

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

@@ -64 +64 @@
     tstar = it*DT-siglim;
     tint = tsint;
-    pddt = 0;
+    pddt = 0.;
     for ( jj = 0; jj < 600; jj++){
       if (tint > (tstar+siglim)){break;}
@@ -90 +90 @@
     //    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 +98 @@
     pds[it] = pd*snormfac;  // gaussian integral lookup table for snow fraction
   }
-  pds[itm+1] = 0;
+  pds[itm+1] = 0.;
   //     *******END initialize PDD
   

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

@@ -2088 +2088 @@
    int    i,iqj,imonth;
    double agd[NUMVERTICES];             // surface mass balance
+   double agd2[NUMVERTICES]; 
    double saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    double smelt[NUMVERTICES] = {0};     // yearly melt
@@ -2097 +2098 @@
 
    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 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
@@ -2123 +2124 @@
    double Tsurf[NUMVERTICES] = {0};     // average annual temperature    
    
-   double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
-   double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12];
-   double deltm=1/12;
-   int    ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
+   double h[NUMVERTICES],s[NUMVERTICES]; // ,b[NUMVERTICES]
+   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.
@@ -2136 +2137 @@
    double fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
    double pddtj[NUMVERTICES], hmx2;
-
+int iii;
+
+   /*Recover monthly temperatures and precipitation*/
+   Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
+   Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
+   GaussTria* gauss=new GaussTria();
+   double time,yts;
+   this->parameters->FindParam(&time,TimeEnum);
+   this->parameters->FindParam(&yts,ConstantsYtsEnum);
+   for(int month=0;month<12;month++) {
+     for(int iv=0;iv<NUMVERTICES;iv++) {
+       gauss->GaussVertex(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
+     }
+   }
+// if(this->id==1) {
+//   input2->Echo();}
+   // if(this->id==1){
+   //   printf("-----------------time ------------------------\n");
+   //     printf("%f \n",time);
+   //   printf("-----------------precipitation ------------------------\n");
+   //   for (iii=0; iii<3; iii++){
+   //     printf("%f %f %f %f\n",monthlyprec[iii][1],monthlyprec[iii][3],monthlyprec[iii][5],monthlyprec[iii][11]);
+   //   }
+     //  printf("----------------monthlytemperatures ------------------------\n");
+    //   for (iii=0; iii<3; iii++){
+    //     printf("%f %f %f %f %f %f\n",monthlytemperatures[iii][12],monthlytemperatures[iii][13],monthlytemperatures[iii][15],monthlytemperatures[iii][17],monthlytemperatures[iii][20],monthlytemperatures[iii][24]);
+    //   }
+    // }
    /*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();
-        double time,yts;
-        this->parameters->FindParam(&time,TimeEnum);
-        this->parameters->FindParam(&yts,ConstantsYtsEnum);
-        for(int month=0;month<12;month++){
-                for(int iv=0;iv<NUMVERTICES;iv++){
-                        gauss->GaussVertex(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;i<NUMVERTICES;i++)
-     for(imonth=0;imonth<12;imonth++){
-       t[i][imonth]=ttmp[i];
-       prec[i][imonth]=prectmp[i];
-   }
-
+   //GetInputListOnVertices(&monthlyprec[0],SurfaceforcingsPrecipitationEnum);
+   
    /*Get material parameters :*/
    rho_ice=matpar->GetRhoIce();
@@ -2180 +2191 @@
      imonth =  ismon[iqj];
      for (i = 0; i < NUMVERTICES; i++){
-       st=(s[i]-s0t[i])/1000;
-       tstar = t[i][imonth] - lapser *max(st,sealev);
+       st=(s[i]-s0t[i])/1000.;
+       tstar = monthlytemperatures[i][imonth] - lapser *max(st,sealev);
        Tsurf[i] = tstar*deltm+Tsurf[i];        
        
@@ -2192 +2203 @@
        
        /******exp des/elev precip reduction*******/
-       sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
+       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;           
+       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;
-       
+   // if(time==473040000){
+   // if(this->id==2) {
+   //   printf("----------------qm----------------------\n");
+   //   //printf("%f %f %f \n",qm[0],qm[1],qm[2]);
+   //   printf("%d %f \n",i,qm[i]);
+   //   printf("%f \n",qmpt);
+   //   printf("%f \n",q);
+   //   printf("%f \n",monthlyprec[i][imonth]);
+   //   printf("%f \n",sconv);
+   //   printf("%f \n",pd);
+   //   printf("%f \n",tstar);
+   //   printf("%f \n",Tsurf[i]);
+
+   //   //printf("%f\n",saccu[i]);
+   //   //printf("%f\n",agd[i]);
+   //  }}      
        /*********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 (iqj>5 &&  iqj<9){ Tsum[i]=Tsum[i]+tstar;} 
        if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
        else if (tstar> -siglim){
@@ -2213 +2239 @@
      }
    } // end of seasonal loop 
-   
+ 
    //******************************************************************
    for(i=0;i<NUMVERTICES;i++){
@@ -2224 +2250 @@
        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;
@@ -2306 +2332 @@
      if(Tsurf[i]<0) {
        Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
-     
+ 
      agd[i] = -smelt[i]+saccu[i];
+// if(this->id==2){
+//   printf("----------------melt accu----------------------\n");
+//   printf("%f\n",-smelt[i]);
+//   printf("%f\n",saccu[i]);
+//   printf("%f\n",agd[i]);
+//  }
+     agd[i] = agd[i]/yts;
      pddtj[i]=pddt;
-     
+   }       //end of the for loop over the vertices
+
      /*Update inputs*/    
-     this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
+     this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
      // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
-     this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
-     
-   }       //end of the for loop over the vertices
+     // this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
+
+ 
+   // if(this->id==2){
+   //   for(iii=0;iii<3;iii++){agd2[iii]=agd[iii]*yts;}
+   //   printf("----------------surface mass balance 2------------------------\n");
+   //     printf("%f %f %f\n",agd2[0],agd2[1],agd2[2]);
+   // }
+
 }
 /*}}}*/