From d6c3304f87725b04e2e664b37d1875cc43aa7e97 Mon Sep 17 00:00:00 2001
From: Tim Moore <timoore@redhat.com>
Date: Mon, 26 Oct 2009 06:54:39 +0100
Subject: [PATCH] Generate sea level pressure from metar

This needs to account for the current temperature.
---
 src/Environment/atmosphere.cxx       | 15 +++++++++++
 src/Environment/atmosphere.hxx       |  9 +++++++
 src/Environment/environment_ctrl.cxx | 40 ++++++++++++++++++++++------
 3 files changed, 56 insertions(+), 8 deletions(-)

diff --git a/src/Environment/atmosphere.cxx b/src/Environment/atmosphere.cxx
index a15c943ed..099669586 100644
--- a/src/Environment/atmosphere.cxx
+++ b/src/Environment/atmosphere.cxx
@@ -260,6 +260,21 @@ double FGAtmo::QNH(const double field_elev, const double field_press) {
     return rslt;
 }
 
+// Invert the QNH calculation to get the field pressure from a metar
+// report.
+// field pressure _in pascals_ 
+//  ... caller gets to convert to inHg or millibars
+// Field elevation in m
+// Altimeter setting (QNH) in pascals
+// Valid for fields within the troposphere only.
+double FGAtmo::fieldPressure(const double field_elev, const double qnh)
+{
+    using namespace atmodel;
+    static const double nn = ISA::lam0 * Rgas / g / mm;
+    const double pratio = pow(qnh / ISA::P0, nn);
+    return ISA::P0 * pow(pratio - field_elev * ISA::lam0 / ISA::T0, 1.0 / nn);
+}
+
 void FGAltimeter::dump_stack1(const double Tref) {
     using namespace atmodel;
     const int bs(200);
diff --git a/src/Environment/atmosphere.hxx b/src/Environment/atmosphere.hxx
index 55ccfb3f9..1ee4e1b6b 100644
--- a/src/Environment/atmosphere.hxx
+++ b/src/Environment/atmosphere.hxx
@@ -110,6 +110,15 @@ public:
 // Field pressure in pascals
 // Valid for fields within the troposphere only.
     double QNH(const double field_elev, const double field_press);
+/**
+ * Invert the QNH calculation to get the field pressure from a metar
+ * report. Valid for fields within the troposphere only.
+ * @param field_elev field elevation in m
+ * @param qnh altimeter setting in pascals
+ * @return field pressure _in pascals_. Caller gets to convert to inHg
+ * or millibars 
+ */
+    static double fieldPressure(const double field_elev, const double qnh);
 };
 
 
diff --git a/src/Environment/environment_ctrl.cxx b/src/Environment/environment_ctrl.cxx
index 5b9cfcfb9..330f8e010 100644
--- a/src/Environment/environment_ctrl.cxx
+++ b/src/Environment/environment_ctrl.cxx
@@ -34,6 +34,7 @@
 #include <Main/fg_props.hxx>
 #include <Main/util.hxx>
 
+#include "atmosphere.hxx"
 #include "fgmetar.hxx"
 #include "environment_ctrl.hxx"
 
@@ -433,6 +434,25 @@ static inline double convert_to_180( double d )
 	return d > 180.0 ? d - 360.0 : d;
 }
 
+// Return the sea level pressure for a metar observation, in inHg.
+// This is different from QNH because it accounts for the current
+// temperature at the observation point.
+// metarPressure in inHg
+// fieldHt in ft
+// fieldTemp in C
+
+static double reducePressureSl(double metarPressure, double fieldHt,
+                               double fieldTemp)
+{
+    double elev = fieldHt * SG_FEET_TO_METER;
+    double fieldPressure
+        = FGAtmo::fieldPressure(elev, metarPressure * atmodel::inHg);
+    double slPressure = P_layer(0, elev, fieldPressure,
+                                fieldTemp + atmodel::freezing,
+                                atmodel::ISA::lam0);
+    return slPressure / atmodel::inHg;
+}
+
 void
 FGMetarCtrl::update(double dt)
 {
@@ -444,6 +464,7 @@ FGMetarCtrl::update(double dt)
 
 	bool reinit_required = false;
 	bool layer_rebuild_required = false;
+        double station_elevation_ft = station_elevation_n->getDoubleValue();
 
 	if (first_update) {
 		double dir = base_wind_dir_n->getDoubleValue()+magnetic_variation_n->getDoubleValue();
@@ -455,7 +476,10 @@ FGMetarCtrl::update(double dt)
 		fgDefaultWeatherValue("visibility-m", metarvis);
 
 		double metarpressure = pressure_n->getDoubleValue();
-		fgDefaultWeatherValue("pressure-sea-level-inhg", metarpressure);
+		fgDefaultWeatherValue("pressure-sea-level-inhg",
+                                      reducePressureSl(metarpressure,
+                                                       station_elevation_ft,
+                                                       temperature_n->getDoubleValue()));
 
 		// We haven't already loaded a METAR, so apply it immediately.
 		vector<SGPropertyNode_ptr> layers = clouds_n->getChildren("layer");
@@ -597,8 +621,11 @@ FGMetarCtrl::update(double dt)
 
 		double pressure = boundary_sea_level_pressure_n->getDoubleValue();
 		double metarpressure = pressure_n->getDoubleValue();
-		if( pressure != metarpressure ) {
-			pressure = interpolate_val( pressure, metarpressure, MaxPressureChangeInHgSec );
+                double newpressure = reducePressureSl(metarpressure,
+                                                      station_elevation_ft,
+                                                      temperature_n->getDoubleValue());
+		if( pressure != newpressure ) {
+			pressure = interpolate_val( pressure, newpressure, MaxPressureChangeInHgSec );
 			fgDefaultWeatherValue("pressure-sea-level-inhg", pressure);
 			reinit_required = true;
 		}
@@ -666,11 +693,8 @@ FGMetarCtrl::update(double dt)
 			}
 		}
 	}
-	{
-		double station_elevation_ft = station_elevation_n->getDoubleValue();
-		set_temp_at_altitude(temperature_n->getDoubleValue(), station_elevation_ft);
-		set_dewpoint_at_altitude(dewpoint_n->getDoubleValue(), station_elevation_ft);
-	}
+        set_temp_at_altitude(temperature_n->getDoubleValue(), station_elevation_ft);
+        set_dewpoint_at_altitude(dewpoint_n->getDoubleValue(), station_elevation_ft);
 	//TODO: check if temperature/dewpoint have changed. This requires reinit.
 
 	// Force an update of the 3D clouds