diff --git a/src/Calibrator/Phase.cpp b/src/Calibrator/Phase.cpp index a93a1b57..8e99c64c 100644 --- 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 index a9eff7be..c864a755 100644 --- 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 index b510cd34..1b5d3b21 100644 --- 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 index a1d88bee..2b2616f8 100644 --- 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 index ccfa38cb..caf6186e 100644 --- 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(); }