Changeset 18717 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

10/31/14 10:58:45 (10 years ago)

Author:

lemorzad

Message:

NEW: modified lapse rate to PDD scheme

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (3 diffs)

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

@@ -1916 +1916 @@
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
-   IssmDouble h[NUMVERTICES],s[NUMVERTICES]; // ,b
-   IssmDouble rho_water,rho_ice,desfac,s0p;
+   IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+   IssmDouble h[NUMVERTICES],s[NUMVERTICES];
+   IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
+   IssmDouble PfacTime,TdiffTime,sealevTime;
 
    /*Recover monthly temperatures and precipitation*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
+   Input*     input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);  _assert_(input3);
+   Input*     input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);   _assert_(input4);
    GaussPenta* gauss=new GaussPenta();
    IssmDouble time,yts;
@@ -1933 +1937 @@
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
        monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
+       input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
+       TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
+       input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
+       TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
      }
    } 
+
+  /*Recover Pfac, Tdiff and sealev at time t:*/
+  this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
+  this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
+  this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
 
   /*Recover info at the vertices: */
@@ -1944 +1957 @@
   rho_water=matpar->GetRhoFreshwater();
 
-  /*Get desertification effect parameters*/
+  /*Get other pdd parameters*/
   desfac=matpar->GetDesFac();
   s0p=matpar->GetS0p();
+  s0t=matpar->GetS0t();
+  rlaps=matpar->GetRlaps();
+  rlapslgm=matpar->GetRlapslgm();
 
    /*measure the surface mass balance*/
    for (int iv = 0; iv < NUMVERTICES; iv++){
-     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, 
-                                  signorm, yts, h[iv], s[iv], rho_ice, rho_water, desfac, s0p);
+     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],  
+                                  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
+                                  pdds,pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water,
+                                  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
    }