Use standard atmosphere for pressure to estimate wetbulb

* To save memory, use a standard atmosphere in the calculation of wetbulb temperature when deriving precipitation Phase. * This saves memory since the pressure forecast field is no longer needed * The impact on the computed wetbulb temperature is likely small * In the future add an option to compute the "true" pressure
metno · Feb 11, 2015 · 42fba0c · 42fba0c
1 parent 7647a7e
commit 42fba0c
Show file tree

Hide file tree

Showing 5 changed files with 38 additions and 5 deletions.
diff --git a/src/Calibrator/Phase.cpp b/src/Calibrator/Phase.cpp
@@ -3,10 +3,12 @@
 #include "../Util.h"
 #include "../File/File.h"
 #include "../ParameterFile.h"
+#include "../Downscaler/Pressure.h"
 CalibratorPhase::CalibratorPhase(const ParameterFile* iParameterFile) :
       Calibrator(),
       mParameterFile(iParameterFile),
       mMinPrecip(0.2),
+      mEstimatePressure(true),
       mUseWetbulb(1) {
    if(iParameterFile->getNumParameters() != 2) {
       Util::error("Parameter file '" + iParameterFile->getFilename() + "' does not have two datacolumns");
@@ -19,6 +21,7 @@ bool CalibratorPhase::calibrateCore(File& iFile) const {
    int nEns = iFile.getNumEns();
    int nTime = iFile.getNumTime();
    iFile.initNewVariable(Variable::Phase);
+   vec2 elevs = iFile.getElevs();
 
 
    // Loop over offsets
@@ -33,19 +36,25 @@ bool CalibratorPhase::calibrateCore(File& iFile) const {
       FieldPtr rh;
       if(mUseWetbulb) {
          // Only load these fields if they are to be used, to save memory
-         pressure = iFile.getField(Variable::P, t);
-         rh       = iFile.getField(Variable::RH, t);
+         rh = iFile.getField(Variable::RH, t);
+         if(!mEstimatePressure)
+            pressure = iFile.getField(Variable::P, t);
       }
 
       #pragma omp parallel for
       for(int i = 0; i < nLat; i++) {
          for(int j = 0; j < nLon; j++) {
+            float currElev = elevs[i][j];
             for(int e = 0; e < nEns; e++) {
                float currDryTemp  = (*temp)(i,j,e);
                float currTemp     = currDryTemp;
                float currPrecip   = (*precip)(i,j,e);
                if(mUseWetbulb) {
-                  float currPressure = (*pressure)(i,j,e);
+                  float currPressure;
+                  if(mEstimatePressure)
+                     currPressure = DownscalerPressure::calcPressure(0, 101325, currElev);
+                  else
+                     currPressure = (*pressure)(i,j,e);
                   float currRh       = (*rh)(i,j,e);
                   float currWetbulb  = getWetbulb(currDryTemp, currPressure, currRh);
                   if(Util::isValid(snowSleetThreshold) && Util::isValid(sleetRainThreshold)

diff --git a/src/Calibrator/Phase.h b/src/Calibrator/Phase.h
@@ -25,5 +25,8 @@ class CalibratorPhase : public Calibrator {
       const ParameterFile* mParameterFile;
       float mMinPrecip;
       bool mUseWetbulb;
+      //! If true compute pressure based on standard atmosphere (instead of using forecasted data)
+      //! This is likely a good enough approximation when computing wetbulb temperature and saves memory.
+      bool mEstimatePressure;
 };
 #endif
diff --git a/src/Downscaler/Pressure.cpp b/src/Downscaler/Pressure.cpp
@@ -44,8 +44,9 @@ void DownscalerPressure::downscaleCore(const File& iInput, File& iOutput) const
                   ofield(i,j,e) = ifield(Icenter,Jcenter,e);
                }
                else {
-                  float dElev = currElev - nearestElev;
-                  ofield(i,j,e) = ifield(Icenter,Jcenter,e) * exp(mConstant * (dElev));
+                  float nearestPressure = ifield(Icenter,Jcenter,e);
+                  float currPressure = calcPressure(nearestElev, nearestPressure, currElev);
+                  ofield(i,j,e) = currPressure;
                }
             }
          }
@@ -58,3 +59,13 @@ std::string DownscalerPressure::description() {
       << "                                elevation difference and a standard atmosphere." << std::endl;
    return ss.str();
 }
+
+float DownscalerPressure::calcPressure(float iElev0, float iPressure0, float iElev1) {
+   if(Util::isValid(iElev0) && Util::isValid(iPressure0) && Util::isValid(iElev1)) {
+      float dElev = iElev1 - iElev0;
+      return iPressure0 * exp(mConstant * (dElev));
+   }
+   else {
+      return Util::MV;
+   }
+}
diff --git a/src/Downscaler/Pressure.h b/src/Downscaler/Pressure.h
@@ -13,6 +13,7 @@ class DownscalerPressure : public Downscaler {
       DownscalerPressure(Variable::Type iVariable);
       static std::string description();
       std::string name() const {return "pressure";};
+      static float calcPressure(float iElev0, float iPressure0, float iElev1);
    private:
       void downscaleCore(const File& iInput, File& iOutput) const;
       static const float mConstant;

diff --git a/src/Testing/DownscalerPressure.cpp b/src/Testing/DownscalerPressure.cpp
@@ -61,6 +61,15 @@ namespace {
       EXPECT_FLOAT_EQ(294.83279, toT(0,2,0)); // 304 347m->600m
       EXPECT_FLOAT_EQ(320.89322, toT(0,3,0)); // 304 347m->-100m
    }
+   TEST_F(TestDownscalerPressure, calcPressure) {
+      EXPECT_FLOAT_EQ(101325, DownscalerPressure::calcPressure(0, 101325, 0));
+      EXPECT_FLOAT_EQ(89777.391, DownscalerPressure::calcPressure(0, 101325, 1000));
+      EXPECT_FLOAT_EQ(67504.562, DownscalerPressure::calcPressure(300, 70000, 600));
+
+      EXPECT_FLOAT_EQ(Util::MV, DownscalerPressure::calcPressure(Util::MV, 101325, 89874.6));
+      EXPECT_FLOAT_EQ(Util::MV, DownscalerPressure::calcPressure(0, Util::MV, 89874.6));
+      EXPECT_FLOAT_EQ(Util::MV, DownscalerPressure::calcPressure(0, 101325, Util::MV));
+   }
    TEST_F(TestDownscalerPressure, description) {
       DownscalerPressure::description();
    }