source: issm/oecreview/Archive/19101-20495/ISSM-19136-19137.diff@ 20498

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -2573 +2573 @@
    int        i;
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+   IssmDouble yearlytemperatures[NUMVERTICES]={0.};
    IssmDouble tmp[NUMVERTICES];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
+   IssmDouble mavg=1./12.; //factor for monthly average
    
    /*Get material parameters :*/
    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -2589 +2591 @@
   rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
   rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
 
-   /*Recover monthly temperatures and precipitation and Present day and LGm ones*/
+   /*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
    GaussTria* gauss=new GaussTria();
@@ -2597 +2599 @@
    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
+       yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin 
+       monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
      }
-   }
+   }     
 
   /*Recover Pfac, Tdiff and sealev at time t:
     This parameters are used to interpolate the temperature 
@@ -2625 +2627 @@
      
    /*measure the surface mass balance*/
   for (int iv = 0; iv<NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
+    agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
                                   pdds, pds, signorm, yts, h[iv], s[iv],
                                   desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
                                   rho_water,rho_ice);
-   }
+  }
 
    /*Update inputs*/    
    // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
@@ -2647 +2649 @@
    // NewPrecipitationInput->Configure(this->parameters);
 
 
-   this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
+  this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
+  this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    // this->inputs->AddInput(NewTemperatureInput);
    // this->inputs->AddInput(NewPrecipitationInput);
    // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));

../trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -2165 +2165 @@
    int        i;
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+   IssmDouble yearlytemperatures[NUMVERTICES]={0.};
    IssmDouble tmp[NUMVERTICES];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
+   IssmDouble mavg=1./12.; //factor for monthly average
 
    /*Get material parameters :*/
    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -2193 +2195 @@
      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
+       yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin 
+       monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
      }
    } 
@@ -2240 +2243 @@
 
 
 
+   this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
    this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    // this->inputs->AddInput(NewTemperatureInput);
    // this->inputs->AddInput(NewPrecipitationInput);
    // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
     this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
+    this->InputExtrude(TemperatureEnum,-1);
    // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);