Exposing the gfsfile surface height as a member.

radixio/gfsfile.cc

@@ -585,6 +585,9 @@ const std::vector<std::string> &GFSFile::profileTimes() const
 {
   return mProfiles;
 }
+
+float GFSFile::surfaceHeight() const { return mSurfaceHeight; }
+
 std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
                                                int day, int year, int hour,
                                                std::vector<std::string> columns)
@@ -716,13 +719,12 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
   delete rstream;
 
   // SOUND section of Fortran
-  float tpot             = 0.0f;
-  float temp             = 288.15f;  // default to surface temperature
-  float relHum           = 0.0f;
-  bool sfcwnd            = false;
-  float offset           = 0.0f;
-  float plevel           = 0.0f;
-  double surfaceAltitude = 0;
+  float tpot   = 0.0f;
+  float temp   = 288.15f;  // default to surface temperature
+  float relHum = 0.0f;
+  bool sfcwnd  = false;
+  float offset = 0.0f;
+  float plevel = 0.0f;
 
   std::vector<std::vector<float>> results(mHeader.nz);
   for (int ll = 0; ll < mHeader.nz; ++ll)
@@ -804,7 +806,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     if (hIt != columns.end())
     {
       results[ll][hIt - columns.begin()] =
-          hypsometric(msle, plevel, temp) - surfaceAltitude;
+          hypsometric(msle, plevel, temp) - mSurfaceHeight;
     }
 
     // check for time
@@ -846,7 +848,12 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
       if (varb.compare("SPH2") == 0) varb = "SPHU";
       if (varb.compare("PRSS") == 0)
       {
-        surfaceAltitude = hypsometric(msle, vdata[kk][ll], temp);
+        mSurfaceHeight = hypsometric(msle, vdata[kk][ll], temp);
+      }
+      if (varb.compare("SHGT") == 0)
+      {
+        // geopotential surface height
+        mSurfaceHeight = vdata[kk][ll];
       }
       // check whether relative humidity has been requested, and convert
       // specific humidity if this is present instead

radixio/gfsfile.hh

@@ -183,6 +183,7 @@ class RADIX_PUBLIC GFSFile
 
   Label mLabel;
   Header mHeader;
+  float mSurfaceHeight;
   /**
    *  \brief gbl2xy converts global lat&lon to xy coordinates for grid lookup
    */
@@ -282,6 +283,12 @@ class RADIX_PUBLIC GFSFile
    * @return std::vector<std::string>
    */
   const std::vector<std::string> &profileTimes() const;
+
+  /**
+   * @brief surfaceHeight Get the geopotential height 'SHGT' in meters
+   * @return
+   */
+  float surfaceHeight() const;
 };  // class GFSFile
 }  // namespace radix
 #endif /** RADIXIO_GFSFILE_HH_ */