Updating gfsfile to save the surface pressure height as a HGTS.

radixio/gfsfile.cc

@@ -801,14 +801,6 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     // initialize space for results vector
     results[ll] = std::vector<float>(columns.size(), 0.0f);
 
-    // if "HGTS" has been requested add level as the default HGTS
-    auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
-    if (hIt != columns.end())
-    {
-      results[ll][hIt - columns.begin()] =
-          hypsometric(msle, plevel, temp) - mSurfaceHeight;
-    }
-
     // check for time
     auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
     if (timeIt != columns.end())
@@ -846,14 +838,21 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
       // convert surface specific humidity to specific humidity for relhum
       // conversion
       if (varb.compare("SPH2") == 0) varb = "SPHU";
-      if (varb.compare("PRSS") == 0)
+      // only calculate surface height from surface pressure if geopotential
+      // height was already used
+      if (varb.compare("PRSS") == 0 && mSurfaceHeight == 0.f)
       {
         mSurfaceHeight = hypsometric(msle, vdata[kk][ll], temp);
+        varb           = "HGTS";
+        vdata[kk][ll]  = mSurfaceHeight;
       }
       if (varb.compare("SHGT") == 0)
       {
+        // update varb to indicate a height for saving
+        varb = "HGTS";
         // geopotential surface height
         mSurfaceHeight = vdata[kk][ll];
+        radix_line("Geopotential height(" << mSurfaceHeight << ")");
       }
       // check whether relative humidity has been requested, and convert
       // specific humidity if this is present instead

radixio/gfsfile.hh

@@ -183,7 +183,7 @@ class RADIX_PUBLIC GFSFile
 
   Label mLabel;
   Header mHeader;
-  float mSurfaceHeight;
+  float mSurfaceHeight = 0.f;
   /**
    *  \brief gbl2xy converts global lat&lon to xy coordinates for grid lookup
    */