Changeset 18717 for issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.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/shared/Elements/PddSurfaceMassBalance.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp

@@ -6 +6 @@
 #include "../Numerics/numerics.h"
 
-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water, IssmDouble desfac, IssmDouble s0p){ 
+IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, 
+                                IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
+                                IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, 
+                                IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, 
+                                IssmDouble rho_water, IssmDouble desfac, IssmDouble s0t,
+                                IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, 
+                                IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ 
 
   // output:
@@ -20 +26 @@
   IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
-  IssmDouble lapser=6.5, sealev=0.;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
-  // IssmDouble desfac = 0.5;                 // desert elevation factor
-  // IssmDouble s0p=0.;         // should be set to elevation from precip source
-  IssmDouble s0t=0.;         // should be set to elevation from temperature source
+  //IssmDouble  sealev=0.;         // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
+  //IssmDouble  Pfac=0.5,Tdiff=0.5;
+  IssmDouble  rtlaps;
+  // IssmDouble lapser=6.5         // lapse rate
+  // IssmDouble desfac = 0.5;      // desert elevation factor
+  // IssmDouble s0p=0.;            // should be set to elevation from precip source
+  // IssmDouble s0t=0.;         // should be set to elevation from temperature source
+  IssmDouble tdiffh;  
   IssmDouble st;             // elevation between altitude of the temp record and current altitude
   IssmDouble sp;             // elevation between altitude of the prec record and current altitude
@@ -68 +78 @@
 
   // seasonal loop
-    for (iqj = 0; iqj < 12; iqj++){
-      imonth =  ismon[iqj];
-
-      st=(s-s0t)/1000.;
-      tstar = monthlytemperatures[imonth] - lapser *max(st,sealev);
-      Tsurf = tstar*deltm+Tsurf;        
+  for (iqj = 0; iqj < 12; iqj++){
+    imonth =  ismon[iqj];
+    
+    /*********compute lapse rate ****************/
+    st=(s-s0t)/1000.;
+    rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
+    
+    /*********compute Surface temperature *******/
+    tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
+    tstar = tdiffh + TemperaturesPresentday[imonth] - rtlaps *max(st,sealevTime*0.001);
+    Tsurf = tstar*deltm+Tsurf;        
 
       /*********compute PD ****************/