Changeset 12362


Ignore:
Timestamp:
06/04/12 14:33:55 (13 years ago)
Author:
lemorzad
Message:

debuggind pdd methode

Location:
issm/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r12330 r12362  
    22552255void  Penta::PositiveDegreeDay(double* pdds,double* pds,double signorm){
    22562256
    2257 
    22582257   int    i,iqj,imonth;
    2259    double agd[NUMVERTICES];  // surface and basal
     2258   double agd[NUMVERTICES];  // surface mass balance
    22602259   double saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    22612260   double smelt[NUMVERTICES] = {0};     // yearly melt
    … …  
    22662265   double sconv; //rhow_rain/rhoi / 12 months
    22672266
    2268    double  rho_water,rho_ice,density;
    2269    double lapser=6.5/1000, sealev=0;    // lapse rate. degrees per meter.
     2267   double rho_water,rho_ice,density;
     2268   double lapser=6.5/1000., sealev=0.;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
    22702269   double desfac = 0.5;                 //desert elevation factor
    22712270   double s0p[NUMVERTICES]={0};         //should be set to elevation from precip source
    … …  
    22732272   double st;             // elevation between altitude of the temp record and current altitude
    22742273   double sp;             // elevation between altitude of the prec record and current altitude
    2275 
    22762274
    22772275   // PDD and PD constants and variables
    … …  
    22942292   double Tsurf[NUMVERTICES] = {0};     // average annual temperature   
    22952293   
    2296    double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
    2297    double t[NUMVERTICES][12],prec[NUMVERTICES][12];
    2298    double deltm=1/12;
    2299    int    ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
     2294   double h[NUMVERTICES],s[NUMVERTICES]; // ,b[NUMVERTICES]
     2295   double monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
     2296   double deltm=1./12.;
     2297   int    ismon[12]={11,0,1,2,3,4,5,6,7,8,9,10};
    23002298
    23012299   double snwm;  // snow that could have been melted in a year.
    … …  
    23082306   double pddtj[NUMVERTICES], hmx2;
    23092307
     2308   /*Recover monthly temperatures and precipitation*/
     2309   Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
     2310   Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
     2311   GaussPenta* gauss=new GaussPenta();
     2312   double time,yts;
     2313   this->parameters->FindParam(&time,TimeEnum);
     2314   this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2315   for(int month=0;month<12;month++) {
     2316     for(int iv=0;iv<NUMVERTICES;iv++) {
     2317       gauss->GaussVertex(iv);
     2318       input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
     2319       monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
     2320       input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
     2321       monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
     2322     }
     2323   }
    23102324   /*Recover info at the vertices: */
    23112325   GetInputListOnVertices(&h[0],ThicknessEnum);
    23122326   GetInputListOnVertices(&s[0],SurfaceEnum);
    2313    GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
    2314    GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
    2315 
    2316    for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion from Kelvin to celcius
    2317 
    2318    for(i=0;i<NUMVERTICES;i++)
    2319      for(imonth=0;imonth<12;imonth++){
    2320        t[i][imonth]=ttmp[i];
    2321        prec[i][imonth]=prectmp[i];
    2322    }
    23232327
    23242328   /*Get material parameters :*/
    23252329   rho_ice=matpar->GetRhoIce();
    2326    //rho_freshwater=matpar->GetRhoWater();
     2330   rho_water=matpar->GetRhoFreshwater();
    23272331   
    2328    sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
     2332   sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
    23292333     
    23302334     /*PDD constant*/
    … …  
    23392343     imonth =  ismon[iqj];
    23402344     for (i = 0; i < NUMVERTICES; i++){
    2341        st=(s[i]-s0t[i])/1000;
    2342        tstar = t[i][imonth] - lapser *max(st,sealev);
     2345       st=(s[i]-s0t[i])/1000.;
     2346       tstar = monthlytemperatures[i][imonth] - lapser *max(st,sealev);
    23432347       Tsurf[i] = tstar*deltm+Tsurf[i];       
    23442348       
    23452349       /*********compute PD ****************/
    23462350       if (tstar < PDup){
    2347          pd = 1;
     2351         pd = 1.;
    23482352         if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
    23492353       else {
    2350          pd = 0;}
     2354         pd = 0.;}
    23512355       
    23522356       /******exp des/elev precip reduction*******/
    2353        sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
     2357       sp=(s[i]-s0p[i])/1000.; // deselev effect is wrt chng in topo
    23542358       if (sp>0.0){q = exp(-desfac*sp);}
    23552359       else {q = 1.0;}
    23562360       
    2357        qmt[i]= qmt[i] + prec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent
    2358        qmpt= q*prec[i][imonth]*sconv;           
     2361       qmt[i]= qmt[i] + monthlyprec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent per month
     2362       qmpt= q*monthlyprec[i][imonth]*sconv;           
    23592363       qmp[i]= qmp[i] + qmpt;
    23602364       qm[i]= qm[i] + qmpt*pd;
    2361        
     2365
    23622366       /*********compute PDD************/
    23632367       // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
    23642368       // gaussian=T_m, so ndd=-(Tsurf-pdd)
    2365        if (iqj>6 &&  iqj<10){ Tsum[i]=Tsum[i]+tstar;}
     2369       if (iqj>5 &&  iqj<9){ Tsum[i]=Tsum[i]+tstar;}
    23662370       if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
    23672371       else if (tstar> -siglim){
    … …  
    24672471     
    24682472     agd[i] = -smelt[i]+saccu[i];
     2473     agd[i] = agd[i]/yts;
    24692474     pddtj[i]=pddt;
    2470      
    24712475     /*Update inputs*/   
    2472      this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
    2473      // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2474      this->inputs->AddInput(new PentaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
    2475      
     2476     this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
     2477     // this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    24762478   }       //end of the for loop over the vertices
    24772479}

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

    r12346 r12362  
    20882088   int    i,iqj,imonth;
    20892089   double agd[NUMVERTICES];             // surface mass balance
    2090    double agd2[NUMVERTICES];
    20912090   double saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    20922091   double smelt[NUMVERTICES] = {0};     // yearly melt
    … …  
    21372136   double fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
    21382137   double pddtj[NUMVERTICES], hmx2;
    2139 int iii;
     2138
     2139   if(!IsOnBed()) return;
    21402140
    21412141   /*Recover monthly temperatures and precipitation*/
    … …  
    21552155     }
    21562156   }
    2157 // if(this->id==1) {
    2158 //   input2->Echo();}
    2159    // if(this->id==1){
    2160    //   printf("-----------------time ------------------------\n");
    2161    //     printf("%f \n",time);
    2162    //   printf("-----------------precipitation ------------------------\n");
    2163    //   for (iii=0; iii<3; iii++){
    2164    //     printf("%f %f %f %f\n",monthlyprec[iii][1],monthlyprec[iii][3],monthlyprec[iii][5],monthlyprec[iii][11]);
    2165    //   }
    2166      //  printf("----------------monthlytemperatures ------------------------\n");
    2167     //   for (iii=0; iii<3; iii++){
    2168     //     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]);
    2169     //   }
    2170     // }
    2171    /*Recover info at the vertices: */
     2157  /*Recover info at the vertices: */
    21722158   GetInputListOnVertices(&h[0],ThicknessEnum);
    21732159   GetInputListOnVertices(&s[0],SurfaceEnum);
    2174    //GetInputListOnVertices(&monthlyprec[0],SurfaceforcingsPrecipitationEnum);
    2175    
     2160
    21762161   /*Get material parameters :*/
    21772162   rho_ice=matpar->GetRhoIce();
    … …  
    21822167     /*PDD constant*/
    21832168   siglim = 2.5*signorm;
    2184    siglimc = 2.5*signormc; 
     2169   siglimc = 2.5*signormc;
    21852170   siglim0 =  siglim/DT + 0.5;
    21862171   siglim0c =  siglimc/DT + 0.5;
    … …  
    21972182       /*********compute PD ****************/
    21982183       if (tstar < PDup){
    2199          pd = 1;
     2184         pd = 1.;
    22002185         if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
    22012186       else {
    2202          pd = 0;}
     2187         pd = 0.;}
    22032188       
    22042189       /******exp des/elev precip reduction*******/
    … …  
    22112196       qmp[i]= qmp[i] + qmpt;
    22122197       qm[i]= qm[i] + qmpt*pd;
    2213    // if(time==473040000){
    2214    // if(this->id==2) {
    2215    //   printf("----------------qm----------------------\n");
    2216    //   //printf("%f %f %f \n",qm[0],qm[1],qm[2]);
    2217    //   printf("%d %f \n",i,qm[i]);
    2218    //   printf("%f \n",qmpt);
    2219    //   printf("%f \n",q);
    2220    //   printf("%f \n",monthlyprec[i][imonth]);
    2221    //   printf("%f \n",sconv);
    2222    //   printf("%f \n",pd);
    2223    //   printf("%f \n",tstar);
    2224    //   printf("%f \n",Tsurf[i]);
    2225 
    2226    //   //printf("%f\n",saccu[i]);
    2227    //   //printf("%f\n",agd[i]);
    2228    //  }}     
     2198
    22292199       /*********compute PDD************/
    22302200       // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
    … …  
    23342304 
    23352305     agd[i] = -smelt[i]+saccu[i];
    2336 // if(this->id==2){
    2337 //   printf("----------------melt accu----------------------\n");
    2338 //   printf("%f\n",-smelt[i]);
    2339 //   printf("%f\n",saccu[i]);
    2340 //   printf("%f\n",agd[i]);
    2341 //  }
    23422306     agd[i] = agd[i]/yts;
    23432307     pddtj[i]=pddt;
    2344    }       //end of the for loop over the vertices
    2345 
    23462308     /*Update inputs*/   
    23472309     this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
    23482310     // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2349      // this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
    2350 
    2351  
    2352    // if(this->id==2){
    2353    //   for(iii=0;iii<3;iii++){agd2[iii]=agd[iii]*yts;}
    2354    //   printf("----------------surface mass balance 2------------------------\n");
    2355    //     printf("%f %f %f\n",agd2[0],agd2[1],agd2[2]);
    2356    // }
    2357 
     2311   }       //end of the for loop over the vertices
    23582312}
    23592313/*}}}*/

  • issm/trunk/src/m/model/extrude.m

    r11995 r12362  
    143143md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node');
    144144md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
     145md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node');
    145146
    146147%results
Note: See TracChangeset for help on using the changeset viewer.