Updating hpaToAtitdue to account for molar mass of dry air.

radixio/gfsfile.cc

@@ -803,12 +803,8 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
     if (hIt != columns.end())
     {
-      float v_temp = radix::virtualTemperature(temp, plevel, relHum);
-      radix_tagged_line("Converting temp(" << temp << ") to vtemp(" << v_temp
-                                           << ")");
-
       results[ll][hIt - columns.begin()] =
-          hpaToAltitude(plevel, v_temp, msle) - surfaceAltitude;
+          hpaToAltitude(plevel, relHum, temp, msle) - surfaceAltitude;
     }
     // check for time
     auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
@@ -846,10 +842,8 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
       if (varb.compare("RH2M") == 0) varb = "RELH";
       if (varb.compare("PRSS") == 0)
       {
-        float v_temp = radix::virtualTemperature(temp, vdata[kk][ll], relHum);
-        radix_tagged_line("Converting temp(" << temp << ") to vtemp(" << v_temp
-                                             << ")");
-        surfaceAltitude = hpaToAltitude(vdata[kk][ll], v_temp, msle) - 2.f;
+        surfaceAltitude =
+            hpaToAltitude(vdata[kk][ll], relHum, temp, msle) - 2.f;
       }
       // check whether relative humidity has been requested, and convert
       // specific humidity if this is present instead

radixmath/constants.hh

@@ -50,6 +50,15 @@ const Real SPECIFIC_GAS_CONSTANT = 8.3144598;
  * @brief GRAVITATIONAL_ACCELERATION m/s^2
  */
 const Real GRAVITATIONAL_ACCELERATION = 9.80616;
+
+/**
+ * @brief MOLAR_MASS_DRY_AIR kg/mol
+ */
+const Real MOLAR_MASS_DRY_AIR = 0.028964;
+/**
+ * @brief MOLAR_MASS_WATER_VAPER kg/mol
+ */
+const Real MOLAR_MASS_WATER_VAPER = 0.018016;
 /**
  * Type def unsigned long to Identifier
  */

radixmath/tests/tstUtil.cc

@@ -76,7 +76,7 @@ TEST(radix, hpaToAltitude)
 
     // calculate virtual temperature
     temperature     = virtualTemperature(temperature, pressure, relhumidity);
-    double altitude = hpaToAltitude(pressure, temperature);
+    double altitude = hpaToAltitude(pressure, relhumidity, temperature);
     EXPECT_DOUBLE_EQ(blessed_altitude, altitude);
   }
 }

radixmath/util.cc

@@ -156,11 +156,15 @@ Real cspanf(Real value, Real begin, Real end)
   }
 }
 
-Real hpaToAltitude(Real hpa, Real temperature, Real msle)
+Real hpaToAltitude(Real hpa, Real relHum, Real temperature, Real msle)
 {
-  //  return ((1 - std::pow(hpa / msle, 0.19022256039)) * temperature) / 0.0065;
+  // convert temperature to virtual temperature
+  temperature = virtualTemperature(temperature, hpa, relHum);
+
   Real result = (SPECIFIC_GAS_CONSTANT * temperature) /
                 GRAVITATIONAL_ACCELERATION * std::log(msle / hpa);
+  // divide out the kg/mol using the molar mass of dry air and water
+  result = result / MOLAR_MASS_DRY_AIR;
   radix_tagged_line(result << " = hpaToAltitude(" << hpa << ", " << temperature
                            << ", " << msle << ")");
   return result;

radixmath/util.hh

@@ -39,11 +39,12 @@ Real virtualTemperature(Real temperature, Real pressure, Real relHum);
 /**
  * @brief hpaToAltitude converts hectopascals or millibars to altitude (meters)
  * @param hpa hectorpascals in millibars
- * @param temperature in kelvin (default=288.15)
+ * @param relHum in %
+ * @param temperature in kelvin
  * @param msle mean sea level pressure (default=1013.25)
  * @return altitude in meters
  */
-Real RADIX_PUBLIC hpaToAltitude(Real hpa, Real temperature = 288.15,
+Real RADIX_PUBLIC hpaToAltitude(Real hpa, Real relHum, Real temperature,
                                 Real msle = 1013.25);
 
 /**