diff --git a/src/FDM/JSBSim/CMakeLists.txt b/src/FDM/JSBSim/CMakeLists.txt
index f0a62f5ee..04f60458f 100644
--- a/src/FDM/JSBSim/CMakeLists.txt
+++ b/src/FDM/JSBSim/CMakeLists.txt
@@ -23,7 +23,6 @@ set(HEADERS
     input_output/FGOutputFG.h
     input_output/FGOutputFile.h
     input_output/FGOutputSocket.h
-    input_output/FGUDPOutputSocket.h
     input_output/FGOutputTextFile.h
     input_output/FGOutputType.h
     input_output/FGModelLoader.h
@@ -40,6 +39,8 @@ set(HEADERS
     math/FGRealValue.h
     math/FGRungeKutta.h
     math/FGTable.h
+    math/FGFunctionValue.h
+    math/FGTemplateFunc.h
     models/FGAccelerations.h
     models/FGAerodynamics.h
     models/FGAircraft.h
@@ -117,7 +118,6 @@ set(SOURCES
     input_output/FGOutputFG.cpp
     input_output/FGOutputFile.cpp
     input_output/FGOutputSocket.cpp
-    input_output/FGUDPOutputSocket.cpp
     input_output/FGOutputTextFile.cpp
     input_output/FGOutputType.cpp
     input_output/FGModelLoader.cpp
@@ -191,6 +191,8 @@ set(SOURCES
 
 
 add_library(JSBSim STATIC ${SOURCES} ${HEADERS})
+set(VERSION_MESSAGE "compiled from FlightGear ${FLIGHTGEAR_VERSION}")
+add_definitions("-DJSBSIM_VERSION=\"${VERSION_MESSAGE}\"")
 
 target_link_libraries(JSBSim SimGearCore)
 target_include_directories(JSBSim PRIVATE ${CMAKE_SOURCE_DIR}/src/FDM/JSBSim)
diff --git a/src/FDM/JSBSim/FGFDMExec.cpp b/src/FDM/JSBSim/FGFDMExec.cpp
index 26161d94f..10e7ba80b 100644
--- a/src/FDM/JSBSim/FGFDMExec.cpp
+++ b/src/FDM/JSBSim/FGFDMExec.cpp
@@ -169,6 +169,7 @@ FGFDMExec::FGFDMExec(FGPropertyManager* root, unsigned int* fdmctr) : Root(root)
   instance->Tie("simulation/disperse", this, &FGFDMExec::GetDisperse);
   instance->Tie("simulation/randomseed", this, (iPMF)&FGFDMExec::SRand, &FGFDMExec::SRand, false);
   instance->Tie("simulation/terminate", (int *)&Terminate);
+  instance->Tie("simulation/pause", (int *)&holding);
   instance->Tie("simulation/sim-time-sec", this, &FGFDMExec::GetSimTime);
   instance->Tie("simulation/dt", this, &FGFDMExec::GetDeltaT);
   instance->Tie("simulation/jsbsim-debug", this, &FGFDMExec::GetDebugLevel, &FGFDMExec::SetDebugLevel);
@@ -540,7 +541,6 @@ void FGFDMExec::LoadModelConstants(void)
   Aerodynamics->in.Wingspan      = Aircraft->GetWingSpan();
   Auxiliary->in.Wingspan         = Aircraft->GetWingSpan();
   Auxiliary->in.Wingchord        = Aircraft->Getcbar();
-  Auxiliary->in.PitotAngle       = Aircraft->GetPitotAngle();
   GroundReactions->in.vXYZcg     = MassBalance->GetXYZcg();
 
   LoadPlanetConstants();
diff --git a/src/FDM/JSBSim/FGJSBBase.cpp b/src/FDM/JSBSim/FGJSBBase.cpp
index fa587ca53..a57d52898 100644
--- a/src/FDM/JSBSim/FGJSBBase.cpp
+++ b/src/FDM/JSBSim/FGJSBBase.cpp
@@ -86,15 +86,11 @@ const double FGJSBBase::psftoinhg = 0.014138;
 const double FGJSBBase::psftopa = 47.88;
 const double FGJSBBase::ktstofps = 1.68781;
 const double FGJSBBase::fpstokts = 1.0/ktstofps;
-const double FGJSBBase::inchtoft = 0.08333333;
+const double FGJSBBase::inchtoft = 1.0/12;
 const double FGJSBBase::in3tom3 = 1.638706E-5;
 const double FGJSBBase::m3toft3 = 1.0/(fttom*fttom*fttom);
 const double FGJSBBase::inhgtopa = 3386.38;
 const double FGJSBBase::fttom = 0.3048;
-double FGJSBBase::Reng = 1716.56;   // Gas constant for Air (ft-lb/slug-R)
-double FGJSBBase::Rstar = 1545.348; // Universal gas constant
-double FGJSBBase::Mair = 28.9645;   //
-const double FGJSBBase::SHRatio = 1.40;
 
 // Note that definition of lbtoslug by the inverse of slugtolb and not
 // to a different constant you can also get from some tables will make
@@ -108,7 +104,7 @@ const double FGJSBBase::kgtolb = 2.20462;
 const double FGJSBBase::kgtoslug = 0.06852168;
 
 const string FGJSBBase::needed_cfg_version = "2.0";
-const string FGJSBBase::JSBSim_version = "1.0 " __DATE__ " " __TIME__ ;
+const string FGJSBBase::JSBSim_version = JSBSIM_VERSION " " __DATE__ " " __TIME__ ;
 
 queue <FGJSBBase::Message> FGJSBBase::Messages;
 FGJSBBase::Message FGJSBBase::localMsg;
diff --git a/src/FDM/JSBSim/FGJSBBase.h b/src/FDM/JSBSim/FGJSBBase.h
index 8fec460c4..0e7929e41 100644
--- a/src/FDM/JSBSim/FGJSBBase.h
+++ b/src/FDM/JSBSim/FGJSBBase.h
@@ -358,10 +358,6 @@ protected:
   static const double m3toft3;
   static const double inhgtopa;
   static const double fttom;
-  static double Reng;         // Specific Gas Constant,ft^2/(sec^2*R)
-  static double Rstar;
-  static double Mair;
-  static const double SHRatio;
   static const double lbtoslug;
   static const double slugtolb;
   static const double kgtolb;
diff --git a/src/FDM/JSBSim/initialization/FGInitialCondition.cpp b/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
index 10ba4d960..ae0b84559 100644
--- a/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
+++ b/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
@@ -185,9 +185,8 @@ void FGInitialCondition::SetVcalibratedKtsIC(double vcas)
   double rhoSL = Atmosphere->GetDensitySL();
   double mach = MachFromVcalibrated(fabs(vcas)*ktstofps, pressure, pressureSL, rhoSL);
   double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
-  double PitotAngle = Aircraft->GetPitotAngle();
 
-  SetVtrueFpsIC(mach * soundSpeed / (cos(alpha+PitotAngle) * cos(beta)));
+  SetVtrueFpsIC(mach * soundSpeed);
   lastSpeedSet = setvc;
 }
 
@@ -703,7 +702,6 @@ void FGInitialCondition::SetAltitudeASLFtIC(double alt)
   double mach0 = vt / soundSpeed;
   double vc0 = VcalibratedFromMach(mach0, pressure, pressureSL, rhoSL);
   double ve0 = vt * sqrt(rho/rhoSL);
-  double PitotAngle = Aircraft->GetPitotAngle();
 
   double geodLatitude = position.GetGeodLatitudeRad();
   altitudeASL=alt;
@@ -720,9 +718,8 @@ void FGInitialCondition::SetAltitudeASLFtIC(double alt)
 
   switch(lastSpeedSet) {
     case setvc:
-      mach0 = MachFromVcalibrated(vc0 * cos(alpha+PitotAngle) * cos(beta),
-                                  pressure, pressureSL, rhoSL);
-      SetVtrueFpsIC(mach0 * soundSpeed / (cos(alpha+PitotAngle) * cos(beta)));
+      mach0 = MachFromVcalibrated(vc0, pressure, pressureSL, rhoSL);
+      SetVtrueFpsIC(mach0 * soundSpeed);
       break;
     case setmach:
       SetVtrueFpsIC(mach0 * soundSpeed);
@@ -843,8 +840,7 @@ double FGInitialCondition::GetVcalibratedKtsIC(void) const
   double pressureSL = Atmosphere->GetPressureSL();
   double rhoSL = Atmosphere->GetDensitySL();
   double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
-  double PitotAngle = Aircraft->GetPitotAngle();
-  double mach = vt * cos(alpha+PitotAngle) * cos(beta) / soundSpeed;
+  double mach = vt / soundSpeed;
 
   return fpstokts * VcalibratedFromMach(mach, pressure, pressureSL, rhoSL);
 }
diff --git a/src/FDM/JSBSim/input_output/FGInputSocket.cpp b/src/FDM/JSBSim/input_output/FGInputSocket.cpp
index 2f1f9ee5d..7ebc6b970 100644
--- a/src/FDM/JSBSim/input_output/FGInputSocket.cpp
+++ b/src/FDM/JSBSim/input_output/FGInputSocket.cpp
@@ -61,7 +61,8 @@ CLASS IMPLEMENTATION
 
 FGInputSocket::FGInputSocket(FGFDMExec* fdmex) :
   FGInputType(fdmex),
-  socket(0)
+  socket(0),
+  SockProtocol(FGfdmSocket::ptTCP)
 {
 }
 
@@ -95,7 +96,7 @@ bool FGInputSocket::InitModel(void)
 {
   if (FGInputType::InitModel()) {
     delete socket;
-    socket = new FGfdmSocket(SockPort);
+    socket = new FGfdmSocket(SockPort, SockProtocol);
 
     if (socket == 0) return false;
     if (!socket->GetConnectStatus()) return false;
diff --git a/src/FDM/JSBSim/input_output/FGInputSocket.h b/src/FDM/JSBSim/input_output/FGInputSocket.h
index d07951320..ba886ee71 100644
--- a/src/FDM/JSBSim/input_output/FGInputSocket.h
+++ b/src/FDM/JSBSim/input_output/FGInputSocket.h
@@ -92,6 +92,7 @@ protected:
 
   unsigned int SockPort;
   FGfdmSocket* socket;
+  FGfdmSocket::ProtocolType SockProtocol;
   std::string data;
 };
 }
diff --git a/src/FDM/JSBSim/input_output/FGInputType.cpp b/src/FDM/JSBSim/input_output/FGInputType.cpp
index 2a7cb7ed7..4aceed6f9 100644
--- a/src/FDM/JSBSim/input_output/FGInputType.cpp
+++ b/src/FDM/JSBSim/input_output/FGInputType.cpp
@@ -81,7 +81,7 @@ void FGInputType::SetIdx(unsigned int idx)
 bool FGInputType::Load(Element* element)
 {
   // Perform base class Load.
-  if(!FGModel::Load(element))
+  if(!FGModel::Load(element, true))
     return false;
 
   // no common attributes yet (see FGOutputType for example
diff --git a/src/FDM/JSBSim/input_output/FGModelLoader.cpp b/src/FDM/JSBSim/input_output/FGModelLoader.cpp
index 68799d58f..25df074f1 100644
--- a/src/FDM/JSBSim/input_output/FGModelLoader.cpp
+++ b/src/FDM/JSBSim/input_output/FGModelLoader.cpp
@@ -46,7 +46,7 @@ using namespace std;
 
 namespace JSBSim {
 
-IDENT(IdSrc, "$Id: FGModelLoader.cpp,v 1.4 2017/02/25 14:23:18 bcoconni Exp $");
+IDENT(IdSrc, "$Id: FGModelLoader.cpp,v 1.5 2017/03/18 16:17:42 bcoconni Exp $");
 IDENT(IdHdr, ID_MODELLOADER);
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -71,7 +71,7 @@ Element_ptr FGModelLoader::Open(Element *el)
       document = XMLFileRead.LoadXMLDocument(path);
       if (document == 0L) {
         cerr << endl << el->ReadFrom()
-             << "Could not open file: " << path << endl;
+             << "Could not open file: " << fname << endl;
         return NULL;
       }
       CachedFiles[path.utf8Str()] = document;
diff --git a/src/FDM/JSBSim/input_output/FGOutputSocket.cpp b/src/FDM/JSBSim/input_output/FGOutputSocket.cpp
index b7ff26873..f81932a3d 100644
--- a/src/FDM/JSBSim/input_output/FGOutputSocket.cpp
+++ b/src/FDM/JSBSim/input_output/FGOutputSocket.cpp
@@ -55,6 +55,7 @@ INCLUDES
 #include "models/FGFCS.h"
 #include "models/atmosphere/FGWinds.h"
 #include "input_output/FGXMLElement.h"
+#include "math/FGPropertyValue.h"
 
 using namespace std;
 
@@ -260,13 +261,11 @@ void FGOutputSocket::PrintHeaders(void)
   if (SubSystems & ssPropulsion && Propulsion->GetNumEngines() > 0)
     socket->Append(Propulsion->GetPropulsionStrings(","));
 
-  if (OutputProperties.size() > 0) {
-    for (unsigned int i=0;i<OutputProperties.size();i++)
-      if (OutputCaptions[i].size() > 0) {
-        socket->Append(OutputCaptions[i]);
-      } else {
-        socket->Append(OutputProperties[i]->GetPrintableName());
-      }
+  for (unsigned int i=0;i<OutputParameters.size();++i) {
+    if (!OutputCaptions[i].empty())
+      socket->Append(OutputCaptions[i]);
+    else
+      socket->Append(OutputParameters[i]->GetPrintableName());
   }
 
   socket->Send();
@@ -378,8 +377,8 @@ void FGOutputSocket::Print(void)
     socket->Append(Propulsion->GetPropulsionValues(","));
   }
 
-  for (unsigned int i=0;i<OutputProperties.size();i++) {
-    socket->Append(OutputProperties[i]->getDoubleValue());
+  for (unsigned int i=0;i<OutputParameters.size();++i) {
+    socket->Append(OutputParameters[i]->GetValue());
   }
 
   socket->Send();
diff --git a/src/FDM/JSBSim/input_output/FGOutputTextFile.cpp b/src/FDM/JSBSim/input_output/FGOutputTextFile.cpp
index 40e81f716..b000f417c 100644
--- a/src/FDM/JSBSim/input_output/FGOutputTextFile.cpp
+++ b/src/FDM/JSBSim/input_output/FGOutputTextFile.cpp
@@ -59,6 +59,7 @@ INCLUDES
 #include "models/FGFCS.h"
 #include "models/atmosphere/FGWinds.h"
 #include "input_output/FGXMLElement.h"
+#include "math/FGPropertyValue.h"
 
 using namespace std;
 
@@ -76,8 +77,6 @@ bool FGOutputTextFile::Load(Element* el)
   if(!FGOutputFile::Load(el))
     return false;
 
-//  PreLoad(el, PropertyManager);
-
   string type = el->GetAttributeValue("type");
   string delim;
   if (type == "TABULAR") {
@@ -234,14 +233,12 @@ bool FGOutputTextFile::OpenFile(void)
     outstream << delimeter;
     outstream << Propulsion->GetPropulsionStrings(delimeter);
   }
-  if (OutputProperties.size() > 0) {
-    for (unsigned int i=0;i<OutputProperties.size();i++) {
-      if (OutputCaptions[i].size() > 0) {
-        outstream << delimeter << OutputCaptions[i];
-      } else {
-        outstream << delimeter << OutputProperties[i]->GetFullyQualifiedName();
-      }
-    }
+
+  for (unsigned int i=0;i<OutputParameters.size();++i) {
+    if (!OutputCaptions[i].empty())
+      outstream << delimeter << OutputCaptions[i];
+    else
+      outstream << delimeter << OutputParameters[i]->GetFullyQualifiedName();
   }
 
   if (PreFunctions.size() > 0) {
@@ -394,8 +391,8 @@ void FGOutputTextFile::Print(void)
   }
 
   outstream.precision(18);
-  for (unsigned int i=0;i<OutputProperties.size();i++) {
-    outstream << delimeter << OutputProperties[i]->getDoubleValue();
+  for (unsigned int i=0;i<OutputParameters.size();++i) {
+    outstream << delimeter << OutputParameters[i]->GetValue();
   }
   for (unsigned int i=0;i<PreFunctions.size();i++) {
     outstream << delimeter << PreFunctions[i]->getDoubleValue();
diff --git a/src/FDM/JSBSim/input_output/FGOutputType.cpp b/src/FDM/JSBSim/input_output/FGOutputType.cpp
index 25d2edd2d..bd3f727f9 100644
--- a/src/FDM/JSBSim/input_output/FGOutputType.cpp
+++ b/src/FDM/JSBSim/input_output/FGOutputType.cpp
@@ -43,6 +43,8 @@ INCLUDES
 #include "FGOutputType.h"
 #include "input_output/FGXMLElement.h"
 #include "input_output/FGPropertyManager.h"
+#include "math/FGTemplateFunc.h"
+#include "math/FGFunctionValue.h"
 
 namespace JSBSim {
 
@@ -81,7 +83,10 @@ FGOutputType::FGOutputType(FGFDMExec* fdmex) :
 
 FGOutputType::~FGOutputType()
 {
-  OutputProperties.clear();
+  vector<FGPropertyValue*>::iterator it;
+  for (it=OutputParameters.begin(); it != OutputParameters.end(); ++it)
+    delete *it;
+
   Debug(1);
 }
 
@@ -133,25 +138,41 @@ bool FGOutputType::Load(Element* element)
     string property_str = property_element->GetDataLine();
     FGPropertyNode* node = PropertyManager->GetNode(property_str);
     if (!node) {
-      cerr << fgred << highint << endl << "  No property by the name "
+      cerr << property_element->ReadFrom()
+           << fgred << highint << endl << "  No property by the name "
            << property_str << " has been defined. This property will " << endl
            << "  not be logged. You should check your configuration file."
            << reset << endl;
     } else {
-      OutputProperties.push_back(node);
-      if (property_element->HasAttribute("caption")) {
-        OutputCaptions.push_back(property_element->GetAttributeValue("caption"));
-      } else {
-        OutputCaptions.push_back("");
+      if (property_element->HasAttribute("apply")) {
+        string function_str = property_element->GetAttributeValue("apply");
+        FGOutput* Output = FDMExec->GetOutput();
+        FGTemplateFunc* f = Output->GetTemplateFunc(function_str);
+        if (f)
+          OutputParameters.push_back(new FGFunctionValue(node, f));
+        else {
+          cerr << property_element->ReadFrom()
+               << fgred << highint << "  No function by the name "
+               << function_str << " has been defined. This property will "
+               << "not be logged. You should check your configuration file."
+               << reset << endl;
+        }
       }
+      else
+        OutputParameters.push_back(new FGPropertyValue(node));
+
+      if (property_element->HasAttribute("caption"))
+        OutputCaptions.push_back(property_element->GetAttributeValue("caption"));
+      else
+        OutputCaptions.push_back("");
     }
     property_element = element->FindNextElement("property");
   }
 
   double outRate = 1.0;
-  if (!element->GetAttributeValue("rate").empty()) {
+  if (element->HasAttribute("rate"))
     outRate = element->GetAttributeValueAsNumber("rate");
-  }
+
   SetRateHz(outRate);
 
   return true;
@@ -204,6 +225,16 @@ double FGOutputType::GetRateHz(void) const
   return 1.0 / (rate * FDMExec->GetDeltaT());
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGOutputType::SetOutputProperties(vector<FGPropertyNode_ptr> & outputProperties)
+{
+  vector<FGPropertyNode_ptr>::iterator it;
+  for (it = outputProperties.begin(); it != outputProperties.end(); ++it)
+    OutputParameters.push_back(new FGPropertyValue(*it));
+}
+
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 //    The bitmasked value choices are as follows:
 //    unset: In this case (the default) JSBSim would only print
@@ -245,9 +276,9 @@ void FGOutputType::Debug(int from)
       if (SubSystems & ssGroundReactions) cout << "    Ground parameters logged" << endl;
       if (SubSystems & ssFCS)             cout << "    FCS parameters logged" << endl;
       if (SubSystems & ssPropulsion)      cout << "    Propulsion parameters logged" << endl;
-      if (OutputProperties.size() > 0)    cout << "    Properties logged:" << endl;
-      for (unsigned int i=0;i<OutputProperties.size();i++) {
-        cout << "      - " << OutputProperties[i]->GetName() << endl;
+      if (!OutputParameters.empty())      cout << "    Properties logged:" << endl;
+      for (unsigned int i=0;i<OutputParameters.size();i++) {
+        cout << "      - " << OutputParameters[i]->GetName() << endl;
       }
     }
   }
diff --git a/src/FDM/JSBSim/input_output/FGOutputType.h b/src/FDM/JSBSim/input_output/FGOutputType.h
index 1f3ccceb1..07599fca5 100644
--- a/src/FDM/JSBSim/input_output/FGOutputType.h
+++ b/src/FDM/JSBSim/input_output/FGOutputType.h
@@ -67,6 +67,7 @@ class FGFCS;
 class FGGroundReactions;
 class FGExternalReactions;
 class FGBuoyantForces;
+class FGPropertyValue;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS DOCUMENTATION
@@ -124,10 +125,7 @@ public:
   /** Set the list of properties that should be output for this output instance.
       @param outputProperties list of properties that should be output
   */
-  void SetOutputProperties(std::vector<FGPropertyNode_ptr> & outputProperties)
-  {
-    OutputProperties = outputProperties;
-  }
+  void SetOutputProperties(std::vector<FGPropertyNode_ptr> & outputProperties);
 
   /** Overwrites the name identifier under which the output will be logged.
       This method is taken into account if it is called before
@@ -199,7 +197,7 @@ public:
 protected:
   unsigned int OutputIdx;
   int SubSystems;
-  std::vector <FGPropertyNode_ptr> OutputProperties;
+  std::vector <FGPropertyValue*> OutputParameters;
   std::vector <std::string> OutputCaptions;
   bool enabled;
 
diff --git a/src/FDM/JSBSim/input_output/FGScript.cpp b/src/FDM/JSBSim/input_output/FGScript.cpp
index 38639a465..58c99fcc3 100644
--- a/src/FDM/JSBSim/input_output/FGScript.cpp
+++ b/src/FDM/JSBSim/input_output/FGScript.cpp
@@ -53,6 +53,7 @@ INCLUDES
 #include "models/FGInput.h"
 #include "math/FGCondition.h"
 #include "math/FGFunction.h"
+#include "math/FGFunctionValue.h"
 
 using namespace std;
 
@@ -88,6 +89,8 @@ FGScript::~FGScript()
     delete Events[i].Condition;
     for (j=0; j<Events[i].Functions.size(); j++)
       delete Events[i].Functions[j];
+    for (j=0; j<Events[i].NotifyProperties.size(); j++)
+      delete Events[i].NotifyProperties[j];
   }
   Events.clear();
 
@@ -206,8 +209,12 @@ bool FGScript::LoadScript(const SGPath& script, double default_dT,
 
   // Now, read output spec if given.
   element = document->FindElement("output");
+  SGPath scriptDir = SGPath(script.dir());
+  if (scriptDir.isNull())
+    scriptDir = SGPath(".");
+
   while (element) {
-    if (!FDMExec->GetOutput()->Load(element))
+    if (!FDMExec->GetOutput()->Load(element, scriptDir))
       return false;
  
     element = document->FindNextElement("output");
@@ -282,14 +289,29 @@ bool FGScript::LoadScript(const SGPath& script, double default_dT,
       while (notify_property_element) {
         notifyPropertyName = notify_property_element->GetDataLine();
 
-        newEvent->NotifyPropertyNames.push_back(notifyPropertyName);
-        newEvent->NotifyProperties.push_back(0);
-          string caption_attribute = notify_property_element->GetAttributeValue("caption");
-          if (caption_attribute.empty()) {
-            newEvent->DisplayString.push_back(notifyPropertyName);
-          } else {
-            newEvent->DisplayString.push_back(caption_attribute);
+        if (notify_property_element->HasAttribute("apply")) {
+          string function_str = notify_property_element->GetAttributeValue("apply");
+          FGOutput* Output = FDMExec->GetOutput();
+          FGTemplateFunc* f = Output->GetTemplateFunc(function_str);
+          if (f)
+            newEvent->NotifyProperties.push_back(new FGFunctionValue(notifyPropertyName, PropertyManager, f));
+          else {
+            cerr << notify_property_element->ReadFrom()
+              << fgred << highint << "  No function by the name "
+              << function_str << " has been defined. This property will "
+              << "not be logged. You should check your configuration file."
+              << reset << endl;
           }
+        }
+        else
+          newEvent->NotifyProperties.push_back(new FGPropertyValue(notifyPropertyName, PropertyManager));
+        
+        string caption_attribute = notify_property_element->GetAttributeValue("caption");
+        if (caption_attribute.empty()) {
+          newEvent->DisplayString.push_back(notifyPropertyName);
+        } else {
+          newEvent->DisplayString.push_back(caption_attribute);
+        }
 
         notify_property_element = notify_element->FindNextElement("property");
       }
@@ -487,13 +509,6 @@ bool FGScript::RunScript(void)
           cout << "    " << thisEvent.Description << endl;
         }
         for (j=0; j<thisEvent.NotifyProperties.size();j++) {
-          if (thisEvent.NotifyProperties[j] == 0) {
-            if (PropertyManager->HasNode(thisEvent.NotifyPropertyNames[j])) {
-              thisEvent.NotifyProperties[j] = PropertyManager->GetNode(thisEvent.NotifyPropertyNames[j]);
-            } else {
-              throw("Could not find property named "+thisEvent.NotifyPropertyNames[j]+" in script.");
-            }
-          }
           cout << "    " << thisEvent.DisplayString[j] << " = " << thisEvent.NotifyProperties[j]->getDoubleValue();
           if (thisEvent.NotifyKML) cout << " <br/>";
           cout << endl;
@@ -657,10 +672,10 @@ void FGScript::Debug(int from)
             } else {
               cout << "  Notifications:" << endl << "    {" << endl;
             }
-            for (unsigned j=0; j<Events[i].NotifyPropertyNames.size();j++) {
-              cout << "      "
-                   << Events[i].NotifyPropertyNames[j]
-                   << endl;
+            for (unsigned j=0; j<Events[i].NotifyProperties.size();j++) {
+					cout << "      "
+						  << Events[i].NotifyProperties[j]->GetPrintableName()
+                    << endl;
             }
             cout << "    }" << endl;
           }
diff --git a/src/FDM/JSBSim/input_output/FGScript.h b/src/FDM/JSBSim/input_output/FGScript.h
index 611648565..95b4a811a 100644
--- a/src/FDM/JSBSim/input_output/FGScript.h
+++ b/src/FDM/JSBSim/input_output/FGScript.h
@@ -60,6 +60,7 @@ namespace JSBSim {
 class FGFDMExec;
 class FGCondition;
 class FGFunction;
+class FGPropertyValue;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS DOCUMENTATION
@@ -227,8 +228,7 @@ private:
     std::string           Description;
     std::vector <FGPropertyNode_ptr>  SetParam;
     std::vector <std::string>  SetParamName;
-    std::vector <FGPropertyNode_ptr>  NotifyProperties;
-    std::vector <std::string>              NotifyPropertyNames;
+    std::vector <FGPropertyValue*>  NotifyProperties;
     std::vector <std::string>              DisplayString;
     std::vector <eAction> Action;
     std::vector <eType>   Type;
diff --git a/src/FDM/JSBSim/input_output/FGUDPInputSocket.cpp b/src/FDM/JSBSim/input_output/FGUDPInputSocket.cpp
index 679083477..e26ee95ce 100644
--- a/src/FDM/JSBSim/input_output/FGUDPInputSocket.cpp
+++ b/src/FDM/JSBSim/input_output/FGUDPInputSocket.cpp
@@ -57,36 +57,21 @@ CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 FGUDPInputSocket::FGUDPInputSocket(FGFDMExec* fdmex) :
-  FGInputType(fdmex),
-  socket(0)
+  FGInputSocket(fdmex), rate(20), oldTimeStamp(0.0)
 {
-  rate = 20;
   SockPort = 5139;
-  oldTimeStamp = 0.0;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-FGUDPInputSocket::~FGUDPInputSocket()
-{
-  delete socket;
+  SockProtocol = FGfdmSocket::ptUDP;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 bool FGUDPInputSocket::Load(Element* el)
 {
-  if (!FGInputType::Load(el))
+  if (!FGInputSocket::Load(el))
     return false;
    
   rate = atoi(el->GetAttributeValue("rate").c_str());
   SetRate(0.5 + 1.0/(FDMExec->GetDeltaT()*rate));
-   
-  SockPort = atoi(el->GetAttributeValue("port").c_str());
-  if (SockPort == 0) {
-    cerr << endl << "No port assigned in input element" << endl;
-    return false;
-  }
   
   Element *property_element = el->FindElement("property");
 
@@ -107,25 +92,8 @@ bool FGUDPInputSocket::Load(Element* el)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-bool FGUDPInputSocket::InitModel(void)
-{
-  if (FGInputType::InitModel()) {
-    delete socket;
-    socket = new FGfdmSocket(SockPort, FGfdmSocket::ptUDP, FGfdmSocket::dIN);
-
-    if (socket == 0) return false;
-    cout << "UDP input socket established on port " << SockPort << endl;
-    return true;
-  }
-
-  return false;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 void FGUDPInputSocket::Read(bool Holding)
 {
-
   if (socket == 0) return;
     
   data = socket->Receive();
@@ -160,9 +128,7 @@ void FGUDPInputSocket::Read(bool Holding)
     for (unsigned int i=1; i<values.size(); i++) {
       InputProperties[i-1]->setDoubleValue(values[i]);
     }
-    
   }
-  
 }
 
 }
diff --git a/src/FDM/JSBSim/input_output/FGUDPInputSocket.h b/src/FDM/JSBSim/input_output/FGUDPInputSocket.h
index 4cc324678..e58d3cd04 100644
--- a/src/FDM/JSBSim/input_output/FGUDPInputSocket.h
+++ b/src/FDM/JSBSim/input_output/FGUDPInputSocket.h
@@ -38,8 +38,7 @@ SENTRY
 INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#include "FGInputType.h"
-#include "input_output/FGfdmSocket.h"
+#include "FGInputSocket.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
@@ -64,26 +63,17 @@ CLASS DOCUMENTATION
 CLASS DECLARATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-class FGUDPInputSocket : public FGInputType
+class FGUDPInputSocket : public FGInputSocket
 {
 public:
   /** Constructor. */
   FGUDPInputSocket(FGFDMExec* fdmex);
 
-  /** Destructor. */
-  ~FGUDPInputSocket();
-
   /** Reads the property names from an XML file.
       @param element The root XML Element of the input file.
   */
   bool Load(Element* el);
 
-  /** Initializes the instance. This method basically opens the socket to which
-      inputs will be directed.
-      @result true if the execution succeeded.
-   */
-  bool InitModel(void);
-
   /// Reads the socket and updates properties accordingly.
   void Read(bool Holding);
 
@@ -92,9 +82,6 @@ protected:
   int rate;
   double oldTimeStamp;
   std::vector<FGPropertyNode_ptr> InputProperties;
-  unsigned int SockPort;
-  FGfdmSocket* socket;
-  std::string data;
 };
 }
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/input_output/FGUDPOutputSocket.cpp b/src/FDM/JSBSim/input_output/FGUDPOutputSocket.cpp
deleted file mode 100644
index 8ec0232af..000000000
--- a/src/FDM/JSBSim/input_output/FGUDPOutputSocket.cpp
+++ /dev/null
@@ -1,146 +0,0 @@
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- Module:       FGUDPOutputSocket.cpp
- Author:       David Culp
- Date started: 03/31/15
- Purpose:      Manage output of property values to a UDP socket
- Called by:    FGOutput
-
- ------------- Copyright (C) 2015 David Culp ----------------
-
- This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU Lesser General Public License as published by the Free Software
- Foundation; either version 2 of the License, or (at your option) any later
- version.
-
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
- details.
-
- You should have received a copy of the GNU Lesser General Public License along with
- this program; if not, write to the Free Software Foundation, Inc., 59 Temple
- Place - Suite 330, Boston, MA  02111-1307, USA.
-
- Further information about the GNU Lesser General Public License can also be found on
- the world wide web at http://www.gnu.org.
-
-FUNCTIONAL DESCRIPTION
---------------------------------------------------------------------------------
-This class sends comma-separated strings over a UDP socket.  The format is that required
-by the QtJSBSim application.
-
-HISTORY
---------------------------------------------------------------------------------
-03/31/15   DC    Created
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-INCLUDES
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-#include <cstring>
-#include <cstdlib>
-
-#include "FGUDPOutputSocket.h"
-#include "FGFDMExec.h"
-#include "input_output/FGPropertyManager.h"
-#include "input_output/FGXMLElement.h"
-
-using namespace std;
-
-namespace JSBSim {
-
-IDENT(IdSrc,"$Id: FGUDPOutputSocket.cpp,v 1.2 2015/08/16 13:19:52 bcoconni Exp $");
-IDENT(IdHdr,ID_UDPOUTPUTSOCKET);
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-CLASS IMPLEMENTATION
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-FGUDPOutputSocket::FGUDPOutputSocket(FGFDMExec* fdmex) :
-  FGOutputType(fdmex),
-  socket(0)
-{
-  FDMExec = fdmex;
-  PropertyManager = fdmex->GetPropertyManager();
-  root = PropertyManager->GetNode(); 
-  root->SetDouble("simulation/null", 0.0);
-  SockName = "localhost";
-  SockPort = 5138;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-FGUDPOutputSocket::~FGUDPOutputSocket()
-{
-  delete socket;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-bool FGUDPOutputSocket::Load(Element* el)
-{
-
-  Element *property_element = el->FindElement("property");
-
-  while (property_element) {
-    string property_str = property_element->GetDataLine();
-    FGPropertyNode* node = PropertyManager->GetNode(property_str);
-    if (!node) {
-      node = PropertyManager->GetNode("simulation/null");
-    }
-    OutputProperties.push_back(node);
-    property_element = el->FindNextElement("property");
-  }
-
-  double outRate = 1.0;
-  if (!el->GetAttributeValue("rate").empty()) {
-    outRate = el->GetAttributeValueAsNumber("rate");
-  }
-  SetRateHz(outRate);
-  
-  SockPort = atoi(el->GetAttributeValue("port").c_str());
-  if (SockPort == 0) {
-    cerr << endl << "No port assigned for output." << endl;
-    return false;
-  }
-
-  return true;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-bool FGUDPOutputSocket::InitModel(void)
-{
-  if (FGOutputType::InitModel()) {
-    delete socket;
-    socket = new FGfdmSocket(SockName, SockPort, FGfdmSocket::ptUDP);
-
-    if (socket == 0) return false;
-    if (!socket->GetConnectStatus()) return false;
-
-    return true;
-  }
-
-  return false;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-void FGUDPOutputSocket::Print(void)
-{
-
-  if (socket == 0) return;
-  if (!socket->GetConnectStatus()) return;
-
-  socket->Clear();
-  socket->Append(FDMExec->GetSimTime());
-
-  for (unsigned int i=0;i<OutputProperties.size();i++) {
-    socket->Append(OutputProperties[i]->getDoubleValue());
-  }
-
-  socket->Send();
-}
-
-}
diff --git a/src/FDM/JSBSim/input_output/FGUDPOutputSocket.h b/src/FDM/JSBSim/input_output/FGUDPOutputSocket.h
deleted file mode 100644
index 9be951e71..000000000
--- a/src/FDM/JSBSim/input_output/FGUDPOutputSocket.h
+++ /dev/null
@@ -1,106 +0,0 @@
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- Header:       FGOutputSocket.h
- Author:       David Culp
- Date started: 03/31/15
-
- ------------- Copyright (C) 2015 David Culp ---------------
-
- This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU Lesser General Public License as published by the Free Software
- Foundation; either version 2 of the License, or (at your option) any later
- version.
-
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
- details.
-
- You should have received a copy of the GNU Lesser General Public License along with
- this program; if not, write to the Free Software Foundation, Inc., 59 Temple
- Place - Suite 330, Boston, MA  02111-1307, USA.
-
- Further information about the GNU Lesser General Public License can also be found on
- the world wide web at http://www.gnu.org.
-
-HISTORY
---------------------------------------------------------------------------------
-03/31/15   DC    Created
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SENTRY
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-#ifndef FGUDPOUTPUTSOCKET_H
-#define FGUDPOUTPUTSOCKET_H
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-INCLUDES
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-#include "FGOutputType.h"
-#include "input_output/FGPropertyManager.h"
-#include "input_output/FGfdmSocket.h"
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-DEFINITIONS
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-#define ID_UDPOUTPUTSOCKET "$Id: FGUDPOutputSocket.h,v 1.1 2015/04/02 02:23:33 dpculp Exp $"
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FORWARD DECLARATIONS
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-namespace JSBSim {
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-CLASS DOCUMENTATION
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-/** Implements the output to a UDP socket. This class outputs data to a socket
-    in a comma-separated strings format. The first string represents a time
-    stamp, and each string thereafter is the value of a property.  If a
-    specified property does not exist, then a zero is sent, so that the
-    number of values sent will always equal the number of properties requested. 
- */
-
-/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-CLASS DECLARATION
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
-class FGUDPOutputSocket : public FGOutputType
-{
-public:
-  /** Constructor. */
-  FGUDPOutputSocket(FGFDMExec* fdmex);
-
-  /** Destructor. */
-  ~FGUDPOutputSocket();
-
-  /** Init the output directives from an XML file.
-      @param element XML Element that is pointing to the output directives
-  */
-  virtual bool Load(Element* el);
-
-  /** Initializes the instance. This method opens the ouput socket.
-      @result true if the execution succeeded.
-   */
-  bool InitModel(void);
-  
-  /// Generates and sends the output datagram.
-  void Print(void);
-
-protected:
-
-  std::string SockName;
-  unsigned int SockPort;
-  FGfdmSocket* socket;
-  FGPropertyManager* PropertyManager;
-  FGPropertyNode* root;
-  FGFDMExec* FDMExec;
-};
-}
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-#endif
-
diff --git a/src/FDM/JSBSim/input_output/FGfdmSocket.cpp b/src/FDM/JSBSim/input_output/FGfdmSocket.cpp
index f6eeb1289..433078481 100644
--- a/src/FDM/JSBSim/input_output/FGfdmSocket.cpp
+++ b/src/FDM/JSBSim/input_output/FGfdmSocket.cpp
@@ -131,46 +131,69 @@ FGfdmSocket::FGfdmSocket(const string& address, int port, int protocol)
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-// assumes UDP socket on localhost, for inbound datagrams
-FGfdmSocket::FGfdmSocket(int port, int protocol, int direction) // assumes UDP
+// assumes TCP or UDP socket on localhost, for inbound datagrams
+FGfdmSocket::FGfdmSocket(int port, int protocol)
 {
   sckt = -1;
   connected = false;
   Protocol = (ProtocolType)protocol;
-  Direction = (DirectionType) direction;
+  string ProtocolName;
  
 #if defined(_MSC_VER) || defined(__MINGW32__)
   if (!LoadWinSockDLL()) return;
 #endif
 
   if (Protocol == ptUDP) {  //use udp protocol
+    ProtocolName = "UDP";
     sckt = socket(AF_INET, SOCK_DGRAM, IPPROTO_UDP);
 #if defined(_MSC_VER) || defined(__MINGW32__)
-	u_long NonBlock = 1; // True
-	ioctlsocket(sckt, FIONBIO, &NonBlock);
+    u_long NonBlock = 1; // True
+    ioctlsocket(sckt, FIONBIO, &NonBlock);
 #else
-	fcntl(sckt, F_SETFL, O_NONBLOCK);
+    fcntl(sckt, F_SETFL, O_NONBLOCK);
 #endif
-    cout << "Creating UDP input socket on port " << port << endl;
   }
+  else {
+    ProtocolName = "TCP";
+    sckt = socket(AF_INET, SOCK_STREAM, 0);
+  }
+  cout << "Creating input " << ProtocolName << " socket on port " << port << endl;
   
-    if (sckt != -1) { 
-      memset(&scktName, 0, sizeof(struct sockaddr_in));
-      scktName.sin_family = AF_INET;
-      scktName.sin_port = htons(port);
+  if (sckt != -1) {
+    memset(&scktName, 0, sizeof(struct sockaddr_in));
+    scktName.sin_family = AF_INET;
+    scktName.sin_port = htons(port);
+    if (Protocol == ptUDP)
       scktName.sin_addr.s_addr = htonl(INADDR_ANY);
-      int len = sizeof(struct sockaddr_in);
-      if (bind(sckt, (struct sockaddr*)&scktName, len) != -1) { 
-        cout << "Successfully bound to UDP input socket on port " << port << endl <<endl;
+    int len = sizeof(struct sockaddr_in);
+    if (bind(sckt, (struct sockaddr*)&scktName, len) != -1) {
+      cout << "Successfully bound to " << ProtocolName << " input socket on port "
+           << port << endl <<endl;
+      if (Protocol == ptTCP) {
+        unsigned long NoBlock = true;
+        if (listen(sckt, 5) >= 0) { // successful listen()
+#if defined(_MSC_VER) || defined(__MINGW32__)
+          ioctlsocket(sckt, FIONBIO, &NoBlock);
+          sckt_in = accept(sckt, (struct sockaddr*)&scktName, &len);
+#else
+          ioctl(sckt, FIONBIO, &NoBlock);
+          sckt_in = accept(sckt, (struct sockaddr*)&scktName, (socklen_t*)&len);
+#endif
+          connected = true;
+        } else {
+          cerr << "Could not listen ..." << endl;
+        }
+      } else
         connected = true;
-      } else {                // unsuccessful
-        cout << "Could not bind to UDP input socket, error = " << errno << endl;
-      }
-    } else {          // unsuccessful
-      cout << "Could not create socket for UDP input, error = " << errno << endl;
+    } else {                // unsuccessful
+        cout << "Could not bind to " << ProtocolName << " input socket, error = "
+             << errno << endl;
     }
-    
-  
+  } else {          // unsuccessful
+      cout << "Could not create " << ProtocolName << " socket for input, error = "
+           << errno << endl;
+  }
+
   Debug(0);
 }
 
@@ -225,49 +248,6 @@ FGfdmSocket::FGfdmSocket(const string& address, int port) // assumes TCP
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGfdmSocket::FGfdmSocket(int port) // assumes TCP
-{
-  connected = false;
-  unsigned long NoBlock = true;
-  Protocol = ptTCP;
-
-  #if defined(_MSC_VER) || defined(__MINGW32__)
-  if (!LoadWinSockDLL()) return;
-  #endif
-
-  sckt = socket(AF_INET, SOCK_STREAM, 0);
-
-  if (sckt >= 0) {  // successful
-    memset(&scktName, 0, sizeof(struct sockaddr_in));
-    scktName.sin_family = AF_INET;
-    scktName.sin_port = htons(port);
-    int len = sizeof(struct sockaddr_in);
-    if (bind(sckt, (struct sockaddr*)&scktName, len) == 0) {   // successful
-      cout << "Successfully bound to socket for input on port " << port << endl;
-      if (listen(sckt, 5) >= 0) { // successful listen()
-        #if defined(_MSC_VER) || defined(__MINGW32__)
-          ioctlsocket(sckt, FIONBIO, &NoBlock);
-          sckt_in = accept(sckt, (struct sockaddr*)&scktName, &len);
-        #else
-          ioctl(sckt, FIONBIO, &NoBlock);
-          sckt_in = accept(sckt, (struct sockaddr*)&scktName, (socklen_t*)&len);
-        #endif
-      } else {
-        cerr << "Could not listen ..." << endl;
-      }
-      connected = true;
-    } else {                // unsuccessful
-      cerr << "Could not bind to socket for input ..." << endl;
-    }
-  } else {          // unsuccessful
-    cerr << "Could not create socket for FDM input, error = " << errno << endl;
-  }
-
-  Debug(0);
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 FGfdmSocket::~FGfdmSocket()
 {
   if (sckt) shutdown(sckt,2);
diff --git a/src/FDM/JSBSim/input_output/FGfdmSocket.h b/src/FDM/JSBSim/input_output/FGfdmSocket.h
index 65151ccf6..8f933f42f 100644
--- a/src/FDM/JSBSim/input_output/FGfdmSocket.h
+++ b/src/FDM/JSBSim/input_output/FGfdmSocket.h
@@ -91,8 +91,7 @@ class FGfdmSocket : public FGJSBBase
 public:
   FGfdmSocket(const std::string&, int);
   FGfdmSocket(const std::string&, int, int);
-  FGfdmSocket(int, int, int);  
-  FGfdmSocket(int);
+  FGfdmSocket(int, int);
   ~FGfdmSocket();
   void Send(void);
   void Send(const char *data, int length);
@@ -109,12 +108,10 @@ public:
   bool GetConnectStatus(void) {return connected;}
 
   enum ProtocolType {ptUDP, ptTCP};
-  enum DirectionType {dIN, dOUT};
  
 private:
   int sckt;
   int sckt_in;
-  DirectionType Direction;
   ProtocolType Protocol;
   struct sockaddr_in scktName;
   struct hostent *host;
diff --git a/src/FDM/JSBSim/math/FGFunction.cpp b/src/FDM/JSBSim/math/FGFunction.cpp
index 4272dc883..65726e922 100644
--- a/src/FDM/JSBSim/math/FGFunction.cpp
+++ b/src/FDM/JSBSim/math/FGFunction.cpp
@@ -33,6 +33,7 @@ INCLUDES
 #include <cstdlib>
 #include <cmath>
 
+#include "simgear/misc/strutils.hxx"
 #include "FGFunction.h"
 #include "FGTable.h"
 #include "FGPropertyValue.h"
@@ -50,6 +51,8 @@ IDENT(IdHdr,ID_FUNCTION);
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
+const double FGFunction::invlog2val = 1.0/log10(2.0);
+
 const std::string FGFunction::property_string = "property";
 const std::string FGFunction::value_string = "value";
 const std::string FGFunction::table_string = "table";
@@ -108,21 +111,21 @@ const std::string FGFunction::ifthen_string = "ifthen";
 const std::string FGFunction::switch_string = "switch";
 const std::string FGFunction::interpolate1d_string = "interpolate1d";
 
-FGFunction::FGFunction(FGPropertyManager* propMan, Element* el, const string& prefix)
-                                      : PropertyManager(propMan), Prefix(prefix)
+FGFunction::FGFunction(FGPropertyManager* PropertyManager, Element* el,
+                       const string& prefix, FGPropertyValue* var)
+  : Prefix(prefix), cached(false), cachedValue(-HUGE_VAL), pCopyTo(0L)
 {
-  Element* element;
-  string operation, property_name;
-  cached = false;
-  cachedValue = -HUGE_VAL;
-  invlog2val = 1.0/log10(2.0);
-  pCopyTo = 0L;
+  Load(PropertyManager, el, var);
+}
 
+void FGFunction::Load(FGPropertyManager* PropertyManager, Element* el,
+                      FGPropertyValue* var)
+{
   Name = el->GetAttributeValue("name");
-  operation = el->GetName();
+  string operation = el->GetName();
 
   if (operation == function_string) {
-    sCopyTo = el->GetAttributeValue("copyto");
+    string sCopyTo = el->GetAttributeValue("copyto");
     if (!sCopyTo.empty()) {
 
       if (sCopyTo.find("#") != string::npos) {
@@ -232,7 +235,7 @@ FGFunction::FGFunction(FGPropertyManager* propMan, Element* el, const string& pr
     cerr << "Bad operation " << operation << " detected in configuration file" << endl;
   }
 
-  element = el->GetElement();
+  Element* element = el->GetElement();
   if (!element && Type != eRandom && Type != eUrandom && Type != ePi) {
     cerr << fgred << highint << endl;
     cerr << "  No element was specified as an argument to the \"" << operation << "\" operation" << endl;
@@ -248,26 +251,36 @@ FGFunction::FGFunction(FGPropertyManager* propMan, Element* el, const string& pr
 
     // data types
     if (operation == property_string || operation == p_string) {
-      property_name = element->GetDataLine();
-      if (property_name.find("#") != string::npos) {
-        if (is_number(Prefix)) {
-          property_name = replace(property_name,"#",Prefix);
+      string property_name = element->GetDataLine();
+
+      if (var && simgear::strutils::strip(property_name) == "#")
+        Parameters.push_back(var);
+      else {
+        if (property_name.find("#") != string::npos) {
+          if (is_number(Prefix)) {
+            property_name = replace(property_name,"#",Prefix);
+          }
+          else
+            cerr << el->ReadFrom()
+                 << fgred << "Illegal use of the special character '#'"
+                 << reset << endl;
+        }
+
+        if (PropertyManager->HasNode(property_name)) {
+          FGPropertyNode* newNode = PropertyManager->GetNode(property_name);
+          Parameters.push_back(new FGPropertyValue( newNode ));
+        } else {
+          // cerr << fgcyan << "Warning: The property " + property_name + " is initially undefined."
+          //      << reset << endl;
+          Parameters.push_back(new FGPropertyValue( property_name,
+                                                    PropertyManager ));
         }
-      }
-      if (PropertyManager->HasNode(property_name)) {
-        FGPropertyNode* newNode = PropertyManager->GetNode(property_name);
-        Parameters.push_back(new FGPropertyValue( newNode ));
-      } else {
-        // cerr << fgcyan << "Warning: The property " + property_name + " is initially undefined."
-        //      << reset << endl;
-        Parameters.push_back(new FGPropertyValue( property_name,
-                                                  PropertyManager ));
       }
     } else if (operation == value_string || operation == v_string) {
       Parameters.push_back(new FGRealValue(element->GetDataAsNumber()));
     } else if (operation == table_string || operation == t_string) {
       Parameters.push_back(new FGTable(PropertyManager, element, Prefix));
-    // operations
+      // operations
     } else if (operation == product_string ||
                operation == difference_string ||
                operation == sum_string ||
@@ -315,28 +328,21 @@ FGFunction::FGFunction(FGPropertyManager* propMan, Element* el, const string& pr
                operation == ifthen_string ||
                operation == switch_string ||
                operation == interpolate1d_string)
-    {
-      Parameters.push_back(new FGFunction(PropertyManager, element, Prefix));
-    } else if (operation != description_string) {
+      {
+        Parameters.push_back(new FGFunction(PropertyManager, element, Prefix, var));
+      } else if (operation != description_string) {
       cerr << "Bad operation " << operation << " detected in configuration file" << endl;
     }
     element = el->GetNextElement();
   }
 
-  bind(el); // Allow any function to save its value
+  bind(el, PropertyManager); // Allow any function to save its value
 
   Debug(0);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGFunction::~FGFunction(void)
-{
-  for (unsigned int i=0; i<Parameters.size(); i++) delete Parameters[i];
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 void FGFunction::cacheValue(bool cache)
 {
   cached = false; // Must set cached to false prior to calling GetValue(), else
@@ -784,7 +790,7 @@ string FGFunction::GetValueAsString(void) const
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void FGFunction::bind(Element* el)
+void FGFunction::bind(Element* el, FGPropertyManager* PropertyManager)
 {
   if ( !Name.empty() ) {
     string tmp;
diff --git a/src/FDM/JSBSim/math/FGFunction.h b/src/FDM/JSBSim/math/FGFunction.h
index 73bab38d2..b704eab44 100644
--- a/src/FDM/JSBSim/math/FGFunction.h
+++ b/src/FDM/JSBSim/math/FGFunction.h
@@ -52,6 +52,7 @@ FORWARD DECLARATIONS
 namespace JSBSim {
 
 class Element;
+class FGPropertyValue;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS DOCUMENTATION
@@ -696,28 +697,32 @@ DECLARATION: FGFunction
 
 // Todo: Does this class need a copy constructor, like FGLGear?
 
-class FGFunction : public FGParameter
+class FGFunction : public FGParameter, public FGJSBBase
 {
 public:
+  /// Default constructor.
+  FGFunction()
+    : cached(false), cachedValue(-HUGE_VAL), pCopyTo(0L) {}
 
-/** Constructor.
+  /** Constructor.
     When this constructor is called, the XML element pointed to in memory by the
     element argument is traversed. If other FGParameter-derived objects (values,
     functions, properties, or tables) are encountered, this instance of the
-    FGFunction object will store a pointer to the found object and pass the relevant
-    Element pointer to the constructor for the new object. In other words, each
-    FGFunction object maintains a list of "child" FGParameter-derived objects which
-    in turn may each contain its own list, and so on. At runtime, each object
-    evaluates its child parameters, which each may have its own child parameters to
-    evaluate.
+    FGFunction object will store a pointer to the found object and pass the
+    relevant Element pointer to the constructor for the new object. In other
+    words, each FGFunction object maintains a list of "child"
+    FGParameter-derived objects which in turn may each contain its own list, and
+    so on. At runtime, each object evaluates its child parameters, which each
+    may have its own child parameters to evaluate.
+
     @param PropertyManager a pointer to the property manager instance.
-    @param element a pointer to the Element object containing the function definition.
-    @param prefix an optional prefix to prepend to the name given to the property
-           that represents this function (if given).
+    @param element a pointer to the Element object containing the function
+                   definition.
+    @param prefix an optional prefix to prepend to the name given to the
+                  property that represents this function (if given).
 */
-  FGFunction(FGPropertyManager* PropertyManager, Element* element, const std::string& prefix="");
-  /// Destructor.
-  virtual ~FGFunction();
+  FGFunction(FGPropertyManager* PropertyManager, Element* element,
+             const std::string& prefix="", FGPropertyValue* var=0L);
 
 /** Retrieves the value of the function object.
     @return the total value of the function. */
@@ -738,12 +743,13 @@ public:
     @param shouldCache specifies whether the function should cache the computed value. */
   void cacheValue(bool shouldCache);
 
+protected:
+  void Load(FGPropertyManager* PropertyManager, Element* element,
+            FGPropertyValue* var);
+  virtual void bind(Element*, FGPropertyManager*);
+
 private:
-  std::vector <FGParameter*> Parameters;
-  FGPropertyManager* const PropertyManager;
-  bool cached;
-  double invlog2val;
-  std::string Prefix;
+  static const double invlog2val;
   static const std::string description_string;
   static const std::string property_string;
   static const std::string value_string;
@@ -799,7 +805,7 @@ private:
   static const std::string ifthen_string;
   static const std::string switch_string;
   static const std::string interpolate1d_string;
-  double cachedValue;
+
   enum functionType {eTopLevel=0, eProduct, eDifference, eSum, eQuotient, ePow, eSqrt, eToRadians,
                      eToDegrees, eExp, eAbs, eSign, eSin, eCos, eTan, eASin, eACos, eATan, eATan2,
                      eMin, eMax, eAvg, eFrac, eInteger, eMod, eRandom, eUrandom, ePi,
@@ -807,12 +813,14 @@ private:
                      eIfThen, eSwitch, eInterpolate1D, eRotation_alpha_local,
                      eRotation_beta_local, eRotation_gamma_local, eRotation_bf_to_wf,
                      eRotation_wf_to_bf} Type;
+  std::string Prefix;
+  bool cached;
+  double cachedValue;
   std::string Name;
-  std::string sCopyTo;        // Property name to copy function value to
+  std::vector <FGParameter_ptr> Parameters;
   FGPropertyNode_ptr pCopyTo; // Property node for CopyTo property string
 
   unsigned int GetBinary(double) const;
-  void bind(Element*);
   void Debug(int from);
 };
 
diff --git a/src/FDM/JSBSim/math/FGFunctionValue.h b/src/FDM/JSBSim/math/FGFunctionValue.h
new file mode 100644
index 000000000..7d4660ed8
--- /dev/null
+++ b/src/FDM/JSBSim/math/FGFunctionValue.h
@@ -0,0 +1,85 @@
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  Header: FGFunctionValue.h
+  Author: Bertrand Coconnier
+  Date started: March 10 2018
+
+  --------- Copyright (C) 2018  B. Coconnier (bcoconni@users.sf.net) -----------
+
+  This program is free software; you can redistribute it and/or modify it under
+  the terms of the GNU Lesser General Public License as published by the Free
+  Software Foundation; either version 2 of the License, or (at your option) any
+  later version.
+
+  This program is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+  details.
+
+  You should have received a copy of the GNU Lesser General Public License along
+  with this program; if not, write to the Free Software Foundation, Inc., 59
+  Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+  Further information about the GNU Lesser General Public License can also be
+  found on the world wide web at http://www.gnu.org.
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  SENTRY
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#ifndef FGFUNCTIONVALUE_H
+#define FGFUNCTIONVALUE_H
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  INCLUDES
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#include "math/FGPropertyValue.h"
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  FORWARD DECLARATIONS
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+namespace JSBSim {
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  CLASS DOCUMENTATION
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+/** Represents a property value on which a function is applied
+    @author Bertrand Coconnier
+*/
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  DECLARATION: FGFunctionValue
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+class FGFunctionValue : public FGPropertyValue
+{
+public:
+
+  FGFunctionValue(FGPropertyNode* propNode, FGTemplateFunc* f)
+    :FGPropertyValue(propNode), function(f) {}
+  FGFunctionValue(std::string propName, FGPropertyManager* propertyManager,
+                  FGTemplateFunc* f)
+    :FGPropertyValue(propName, propertyManager), function(f) {}
+
+  double GetValue(void) const { return function->GetValue(GetNode()); }
+
+  std::string GetName(void) const {
+    return function->GetName() + "(" + FGPropertyValue::GetName() + ")";
+  }
+  std::string GetPrintableName(void) const {
+    return function->GetName() + "(" + FGPropertyValue::GetPrintableName() + ")";
+  }
+  std::string GetFullyQualifiedName(void) const {
+    return function->GetName() + "(" + FGPropertyValue::GetFullyQualifiedName() + ")";
+  }
+
+private:
+  FGTemplateFunc_ptr function;
+};
+
+} // namespace JSBSim
+
+#endif
diff --git a/src/FDM/JSBSim/math/FGParameter.h b/src/FDM/JSBSim/math/FGParameter.h
index c1d66fc20..3b16e5bc8 100644
--- a/src/FDM/JSBSim/math/FGParameter.h
+++ b/src/FDM/JSBSim/math/FGParameter.h
@@ -34,7 +34,8 @@ SENTRY
 INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#include "FGJSBBase.h"
+#include <string>
+#include "simgear/structure/SGSharedPtr.hxx"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
@@ -60,7 +61,7 @@ CLASS DOCUMENTATION
 DECLARATION: FGParameter
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-class FGParameter : public FGJSBBase
+class FGParameter : public SGReferenced
 {
 public:
   virtual ~FGParameter(void) {};
@@ -73,6 +74,8 @@ public:
 protected:
 };
 
+typedef SGSharedPtr<FGParameter> FGParameter_ptr;
+
 } // namespace JSBSim
 
 #endif
diff --git a/src/FDM/JSBSim/math/FGPropertyValue.cpp b/src/FDM/JSBSim/math/FGPropertyValue.cpp
index ff25cb379..dd420aad2 100644
--- a/src/FDM/JSBSim/math/FGPropertyValue.cpp
+++ b/src/FDM/JSBSim/math/FGPropertyValue.cpp
@@ -62,7 +62,7 @@ FGPropertyValue::FGPropertyValue(std::string propName, FGPropertyManager* proper
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-double FGPropertyValue::GetValue(void) const
+FGPropertyNode* FGPropertyValue::GetNode(void) const
 {
   FGPropertyNode* node = PropertyNode;
 
@@ -76,18 +76,45 @@ double FGPropertyValue::GetValue(void) const
     }
   }
 
-  return node->getDoubleValue()*Sign;
+  return node;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGPropertyValue::GetValue(void) const
+{
+
+  return GetNode()->getDoubleValue()*Sign;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 std::string FGPropertyValue::GetName(void) const
 {
-  if (PropertyNode) {
+  if (PropertyNode)
     return PropertyNode->GetName();
-  } else {
+  else
     return PropertyName;
-  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+std::string FGPropertyValue::GetFullyQualifiedName(void) const
+{
+  if (PropertyNode)
+    return PropertyNode->GetFullyQualifiedName();
+  else
+    return PropertyName;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+std::string FGPropertyValue::GetPrintableName(void) const
+{
+  if (PropertyNode)
+    return PropertyNode->GetPrintableName();
+ else
+   return PropertyName;
 }
 
 }
diff --git a/src/FDM/JSBSim/math/FGPropertyValue.h b/src/FDM/JSBSim/math/FGPropertyValue.h
index 265b2fd07..27adb18a6 100644
--- a/src/FDM/JSBSim/math/FGPropertyValue.h
+++ b/src/FDM/JSBSim/math/FGPropertyValue.h
@@ -68,12 +68,16 @@ public:
 
   FGPropertyValue(FGPropertyNode* propNode);
   FGPropertyValue(std::string propName, FGPropertyManager* propertyManager);
-  ~FGPropertyValue() {};
 
-  double GetValue(void) const;
+  virtual double GetValue(void) const;
   void SetNode(FGPropertyNode* node) {PropertyNode = node;}
 
-  std::string GetName(void) const;
+  virtual std::string GetName(void) const;
+  virtual std::string GetFullyQualifiedName(void) const;
+  virtual std::string GetPrintableName(void) const;
+
+protected:
+  FGPropertyNode* GetNode(void) const;
 
 private:
   FGPropertyManager* PropertyManager; // Property root used to do late binding.
diff --git a/src/FDM/JSBSim/math/FGRealValue.cpp b/src/FDM/JSBSim/math/FGRealValue.cpp
index 47693f887..13faa5c2d 100644
--- a/src/FDM/JSBSim/math/FGRealValue.cpp
+++ b/src/FDM/JSBSim/math/FGRealValue.cpp
@@ -29,7 +29,7 @@ INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include "FGRealValue.h"
-
+#include "FGJSBBase.h"
 #include "input_output/string_utilities.h"
 
 using namespace std;
diff --git a/src/FDM/JSBSim/math/FGTable.h b/src/FDM/JSBSim/math/FGTable.h
index d3f04bfae..ff9427066 100644
--- a/src/FDM/JSBSim/math/FGTable.h
+++ b/src/FDM/JSBSim/math/FGTable.h
@@ -240,7 +240,7 @@ combustion_efficiency = Lookup_Combustion_Efficiency->GetValue(equivalence_ratio
 CLASS DECLARATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-class FGTable : public FGParameter
+class FGTable : public FGParameter, public FGJSBBase
 {
 public:
   /// Destructor
diff --git a/src/FDM/JSBSim/math/FGTemplateFunc.h b/src/FDM/JSBSim/math/FGTemplateFunc.h
new file mode 100644
index 000000000..c105b2eed
--- /dev/null
+++ b/src/FDM/JSBSim/math/FGTemplateFunc.h
@@ -0,0 +1,86 @@
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  Header: FGTemplateFunc.h
+  Author: Bertrand Coconnier
+  Date started: March 10 2018
+
+  --------- Copyright (C) 2018  B. Coconnier (bcoconni@users.sf.net) -----------
+
+  This program is free software; you can redistribute it and/or modify it under
+  the terms of the GNU Lesser General Public License as published by the Free
+  Software Foundation; either version 2 of the License, or (at your option) any
+  later version.
+
+  This program is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+  details.
+
+  You should have received a copy of the GNU Lesser General Public License along
+  with this program; if not, write to the Free Software Foundation, Inc., 59
+  Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+  Further information about the GNU Lesser General Public License can also be
+  found on the world wide web at http://www.gnu.org.
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  SENTRY
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#ifndef FGTEMPLATEFUNC_H
+#define FGTEMPLATEFUNC_H
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  INCLUDES
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#include "math/FGFunction.h"
+#include "math/FGPropertyValue.h"
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  FORWARD DECLARATIONS
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+namespace JSBSim {
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  CLASS DOCUMENTATION
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  DECLARATION: FGTemplateFunc
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+class FGTemplateFunc : public FGFunction
+{
+public:
+
+  FGTemplateFunc(FGPropertyManager* PropertyManager, Element* element)
+    : var(0L)
+  {
+    Load(PropertyManager, element, &var);
+    // Since 'var' is a member of FGTemplateFunc, we don't want SGSharedPtr to
+    // destroy 'var' when it would no longer be referenced by any shared
+    // pointers. In order to avoid this, the reference counter is increased an
+    // extra time to make sure that it never reaches 0 when all shared pointers
+    // referencing 'var' are destroyed.
+    get(&var);
+  }
+
+  double GetValue(FGPropertyNode* node) {
+    var.SetNode(node);
+    return FGFunction::GetValue();
+  }
+
+private:
+  /** FGTemplateFunc must not be bound to the property manager. The bind method
+      is therefore overloaded as a no-op */
+  virtual void bind(Element*, FGPropertyManager*) {}
+  FGPropertyValue var;
+};
+
+typedef SGSharedPtr<FGTemplateFunc> FGTemplateFunc_ptr;
+
+} // namespace JSBSim
+
+#endif
diff --git a/src/FDM/JSBSim/models/FGAerodynamics.cpp b/src/FDM/JSBSim/models/FGAerodynamics.cpp
index 8e718fe44..1e0ae7172 100644
--- a/src/FDM/JSBSim/models/FGAerodynamics.cpp
+++ b/src/FDM/JSBSim/models/FGAerodynamics.cpp
@@ -76,7 +76,8 @@ FGAerodynamics::FGAerodynamics(FGFDMExec* FDMExec) : FGModel(FDMExec)
   AxisIdx["Y"] = 1;
   AxisIdx["Z"] = 2;
 
-  axisType = atNone;
+  forceAxisType = atNone;
+  momentAxisType = atNone;
 
   AeroFunctions = new AeroFunctionArray[6];
   AeroFunctionsAtCG = new AeroFunctionArray[6];
@@ -154,7 +155,7 @@ bool FGAerodynamics::Run(bool Holding)
     // Skip the computation if qbar is close to zero to avoid huge values for
     // aero/cl-squared when a non-null lift coincides with a very small aero
     // velocity (i.e. when qbar is close to zero).
-    clsq = (vFw(eLift) + vFwAtCG(eLift))/ (in.Wingarea*in.Qbar);
+    clsq = vFw(eLift) / (in.Wingarea*in.Qbar);
     clsq *= clsq;
   }
 
@@ -179,17 +180,18 @@ bool FGAerodynamics::Run(bool Holding)
 
   if (alphahystmax != 0.0 && alphahystmin != 0.0) {
     if (in.Alpha > alphahystmax) {
-       stall_hyst = 1;
+      stall_hyst = 1;
     } else if (in.Alpha < alphahystmin) {
-       stall_hyst = 0;
+      stall_hyst = 0;
     }
   }
 
   vFw.InitMatrix();
-  vFwAtCG.InitMatrix();
   vFnative.InitMatrix();
   vFnativeAtCG.InitMatrix();
 
+  BuildStabilityTransformMatrices();
+
   for (axis_ctr = 0; axis_ctr < 3; ++axis_ctr) {
     AeroFunctionArray::iterator f;
 
@@ -210,68 +212,46 @@ bool FGAerodynamics::Run(bool Holding)
     }
   }
 
-  // Note that we still need to convert to wind axes here, because it is
-  // used in the L/D calculation, and we still may want to look at Lift
-  // and Drag.
+  switch (forceAxisType) {
+  case atBodyXYZ:       // Forces already in body axes; no manipulation needed
+    vForces = vFnative;
+    vForcesAtCG = vFnativeAtCG;
+    break;
+  case atWind:      // Copy forces into wind axes
+    vFnative(eDrag)*=-1; vFnative(eLift)*=-1;
+    vForces = in.Tw2b*vFnative;
 
-  // JSB 4/27/12 - After use, convert wind axes to produce normal lift
-  // and drag values - not negative ones!
+    vFnativeAtCG(eDrag)*=-1; vFnativeAtCG(eLift)*=-1;
+    vForcesAtCG = in.Tw2b*vFnativeAtCG;
+    break;
+  case atBodyAxialNormal:   // Convert native forces into Axial|Normal|Side system
+    vFnative(eX)*=-1; vFnative(eZ)*=-1;
+    vForces = vFnative;
 
-  // As a clarification, JSBSim assumes that drag and lift values are defined 
-  // in wind axes - BUT with a 180 rotation about the Y axis. That is, lift and
-  // drag will be positive up and aft, respectively, so that they are reported
-  // as positive numbers. However, the wind axes themselves assume that the X
-  // and Z forces are positive forward and down.
+    vFnativeAtCG(eX)*=-1; vFnativeAtCG(eZ)*=-1;
+    vForcesAtCG = vFnativeAtCG;
+    break;
+  case atStability:   // Convert from stability axes to both body and wind axes
+    vFnative(eDrag) *= -1; vFnative(eLift) *= -1;
+    vForces = Ts2b*vFnative;
 
-  switch (axisType) {
-    case atBodyXYZ:       // Forces already in body axes; no manipulation needed
-      vFw = in.Tb2w*vFnative;
-      vForces = vFnative;
-      vFw(eDrag)*=-1; vFw(eLift)*=-1;
-
-      vFwAtCG = in.Tb2w*vFnativeAtCG;
-      vForcesAtCG = vFnativeAtCG;
-      vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
-      break;
-    case atLiftDrag:      // Copy forces into wind axes
-      vFw = vFnative;
-      vFw(eDrag)*=-1; vFw(eLift)*=-1;
-      vForces = in.Tw2b*vFw;
-      vFw(eDrag)*=-1; vFw(eLift)*=-1;
-
-      vFwAtCG = vFnativeAtCG;
-      vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
-      vForcesAtCG = in.Tw2b*vFwAtCG;
-      vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
-      break;
-    case atAxialNormal:   // Convert native forces into Axial|Normal|Side system
-      vFw = in.Tb2w*vFnative;
-      vFnative(eX)*=-1; vFnative(eZ)*=-1;
-      vForces = vFnative;
-
-      vFwAtCG = in.Tb2w*vFnativeAtCG;
-      vFnativeAtCG(eX)*=-1; vFnativeAtCG(eZ)*=-1;
-      vForcesAtCG = vFnativeAtCG;
-      break;
-    default:
-      cerr << endl << "  A proper axis type has NOT been selected. Check "
-                   << "your aerodynamics definition." << endl;
-      exit(-1);
+    vFnativeAtCG(eDrag) *= -1; vFnativeAtCG(eLift) *= -1;
+    vForcesAtCG = Ts2b*vFnativeAtCG;
+    break;
+  default:
+    cerr << endl << "  A proper axis type has NOT been selected. Check "
+         << "your aerodynamics definition." << endl;
+    exit(-1);
   }
-
-  // Calculate lift Lift over Drag
-  if ( fabs(vFw(eDrag) + vFwAtCG(eDrag)) > 0.0)
-    lod = fabs( (vFw(eLift) + vFwAtCG(eLift))/ (vFw(eDrag) + vFwAtCG(eDrag)));
-
-  // Calculate aerodynamic reference point shift, if any. The shift
-  // takes place in the structual axis. That is, if the shift is positive,
-  // it is towards the back (tail) of the vehicle. The AeroRPShift
-  // function should be non-dimensionalized by the wing chord. The
-  // calculated vDeltaRP will be in feet.
+  // Calculate aerodynamic reference point shift, if any. The shift takes place
+  // in the structual axis. That is, if the shift is positive, it is towards the
+  // back (tail) of the vehicle. The AeroRPShift function should be
+  // non-dimensionalized by the wing chord. The calculated vDeltaRP will be in
+  // feet.
   if (AeroRPShift) vDeltaRP(eX) = AeroRPShift->GetValue()*in.Wingchord;
 
-  vDXYZcg(eX) = in.RPBody(eX) - vDeltaRP(eX); // vDeltaRP is given in the structural frame
-  vDXYZcg(eY) = in.RPBody(eY) + vDeltaRP(eY);
+  vDXYZcg(eX) = in.RPBody(eX) - vDeltaRP(eX); // vDeltaRP is given in the
+  vDXYZcg(eY) = in.RPBody(eY) + vDeltaRP(eY); // structural frame.
   vDXYZcg(eZ) = in.RPBody(eZ) - vDeltaRP(eZ);
 
   vMomentsMRC.InitMatrix();
@@ -287,11 +267,50 @@ bool FGAerodynamics::Run(bool Holding)
       vMomentsMRC(axis_ctr+1) += (*f)->GetValue();
     }
   }
-  vMoments = vMomentsMRC + vDXYZcg*vForces; // M = r X F
-  // Now add the "at CG" values to base forces - after the moments have been transferred
+  
+  // Transform moments to bodyXYZ if the moments are specified in stability or
+  // wind axes
+  vMomentsMRCBodyXYZ.InitMatrix();
+  switch (momentAxisType) {
+  case atBodyXYZ:
+    vMomentsMRCBodyXYZ = vMomentsMRC;
+    break;
+  case atStability:
+    vMomentsMRCBodyXYZ = Ts2b*vMomentsMRC;
+    break;
+  case atWind:
+    vMomentsMRCBodyXYZ = in.Tw2b*vMomentsMRC;
+    break;
+  default:
+    cerr << endl << "  A proper axis type has NOT been selected. Check "
+         << "your aerodynamics definition." << endl;
+    exit(-1);
+  }
+
+  vMoments = vMomentsMRCBodyXYZ + vDXYZcg*vForces; // M = r X F
+  
+  // Now add the "at CG" values to base forces - after the moments have been
+  // transferred.
   vForces += vForcesAtCG;
-  vFnative += vFnativeAtCG;
-  vFw += vFwAtCG;
+
+  // Note that we still need to convert to wind axes here, because it is used in
+  // the L/D calculation, and we still may want to look at Lift and Drag.
+  //
+  // JSB 4/27/12 - After use, convert wind axes to produce normal lift and drag
+  // values - not negative ones!
+  //
+  // As a clarification, JSBSim assumes that drag and lift values are defined in
+  // wind axes - BUT with a 180 rotation about the Y axis. That is, lift and
+  // drag will be positive up and aft, respectively, so that they are reported
+  // as positive numbers. However, the wind axes themselves assume that the X
+  // and Z forces are positive forward and down. Same applies to the stability
+  // axes.
+  vFw = in.Tb2w * vForces;
+  vFw(eDrag) *= -1; vFw(eLift) *= -1;
+
+  // Calculate Lift over Drag
+  if ( fabs(vFw(eDrag)) > 0.0)
+    lod = fabs( vFw(eLift) / vFw(eDrag));
 
   RunPostFunctions();
 
@@ -300,6 +319,17 @@ bool FGAerodynamics::Run(bool Holding)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+FGColumnVector3 FGAerodynamics::GetForcesInStabilityAxes(void) const
+{ 
+  FGColumnVector3 vFs = Tb2s*vForces;
+  // Need sign flips since drag is positive and lift is positive in stability axes
+  vFs(eDrag) *= -1; vFs(eLift) *= -1;
+
+  return vFs;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 bool FGAerodynamics::Load(Element *document)
 {
   string axis;
@@ -309,7 +339,7 @@ bool FGAerodynamics::Load(Element *document)
   Name = "Aerodynamics Model: " + document->GetAttributeValue("name");
 
   // Perform base class Pre-Load
-  if (!FGModel::Load(document))
+  if (!FGModel::Load(document, true))
     return false;
 
   DetermineAxisSystem(document); // Determine if Lift/Side/Drag, etc. is used.
@@ -353,7 +383,8 @@ bool FGAerodynamics::Load(Element *document)
       try {
         ca.push_back( new FGFunction(PropertyManager, function_element) );
       } catch (const string& str) {
-        cerr << endl << fgred << "Error loading aerodynamic function in " 
+        cerr << endl << axis_element->ReadFrom()
+             << endl << fgred << "Error loading aerodynamic function in "
              << current_func_name << ":" << str << " Aborting." << reset << endl;
         return false;
       }
@@ -361,7 +392,8 @@ bool FGAerodynamics::Load(Element *document)
         try {
           ca_atCG.push_back( new FGFunction(PropertyManager, function_element) );
         } catch (const string& str) {
-          cerr << endl << fgred << "Error loading aerodynamic function in " 
+          cerr << endl << axis_element->ReadFrom()
+               << endl << fgred << "Error loading aerodynamic function in "
                << current_func_name << ":" << str << " Aborting." << reset << endl;
           return false;
         }
@@ -387,6 +419,9 @@ bool FGAerodynamics::Load(Element *document)
 // a warning message will be given IF the AXIAL|NORMAL specifiers are also given.
 // This is OK, and the warning is due to the SIDE specifier used for both
 // the Lift/Drag and Axial/Normal axis systems.
+// Alternatively the axis name 'X|Y|Z or ROLL|PITCH|YAW' can be specified in 
+// conjunction with a frame 'BODY|STABILITY|WIND', for example:
+// <axis name="X" frame="STABILITY"/> 
 
 void FGAerodynamics::DetermineAxisSystem(Element* document)
 {
@@ -394,41 +429,85 @@ void FGAerodynamics::DetermineAxisSystem(Element* document)
   string axis;
   while (axis_element) {
     axis = axis_element->GetAttributeValue("name");
-    if (axis == "LIFT" || axis == "DRAG") {
-      if (axisType == atNone) axisType = atLiftDrag;
-      else if (axisType != atLiftDrag) {
-        cerr << endl << "  Mixed aerodynamic axis systems have been used in the"
-                     << " aircraft config file. (LIFT DRAG)" << endl;
+    string frame = axis_element->GetAttributeValue("frame");
+    if (axis == "X" || axis == "Y" || axis == "Z") {
+      ProcessAxesNameAndFrame(forceAxisType, axis, frame, axis_element,
+                              "(X Y Z)");
+    } else if (axis == "ROLL" || axis == "PITCH" || axis == "YAW") {
+      ProcessAxesNameAndFrame(momentAxisType, axis, frame, axis_element,
+                              "(ROLL PITCH YAW)");
+    } else if (axis == "LIFT" || axis == "DRAG") {
+      if (forceAxisType == atNone) forceAxisType = atWind;
+      else if (forceAxisType != atWind) {
+        cerr << endl << axis_element->ReadFrom()
+             << endl << "  Mixed aerodynamic axis systems have been used in the"
+             << " aircraft config file. (LIFT DRAG)" << endl;
       }
     } else if (axis == "SIDE") {
-      if (axisType != atNone && axisType != atLiftDrag && axisType != atAxialNormal) {
-        cerr << endl << "  Mixed aerodynamic axis systems have been used in the"
-                     << " aircraft config file. (SIDE)" << endl;
+      if (forceAxisType != atNone && forceAxisType != atWind && forceAxisType != atBodyAxialNormal) {
+        cerr << endl << axis_element->ReadFrom()
+             << endl << "  Mixed aerodynamic axis systems have been used in the"
+             << " aircraft config file. (SIDE)" << endl;
       }
     } else if (axis == "AXIAL" || axis == "NORMAL") {
-      if (axisType == atNone) axisType = atAxialNormal;
-      else if (axisType != atAxialNormal) {
-        cerr << endl << "  Mixed aerodynamic axis systems have been used in the"
-                     << " aircraft config file. (NORMAL AXIAL)" << endl;
+      if (forceAxisType == atNone) forceAxisType = atBodyAxialNormal;
+      else if (forceAxisType != atBodyAxialNormal) {
+        cerr << endl << axis_element->ReadFrom()
+             << endl << "  Mixed aerodynamic axis systems have been used in the"
+             << " aircraft config file. (NORMAL AXIAL)" << endl;
       }
-    } else if (axis == "X" || axis == "Y" || axis == "Z") {
-      if (axisType == atNone) axisType = atBodyXYZ;
-      else if (axisType != atBodyXYZ) {
-        cerr << endl << "  Mixed aerodynamic axis systems have been used in the"
-                     << " aircraft config file. (XYZ)" << endl;
-      }
-    } else if (axis != "ROLL" && axis != "PITCH" && axis != "YAW") { // error
-      cerr << endl << "  An unknown axis type, " << axis << " has been specified"
-                   << " in the aircraft configuration file." << endl;
+    } else { // error
+      cerr << endl << axis_element->ReadFrom()
+           << endl << "  An unknown axis type, " << axis << " has been specified"
+           << " in the aircraft configuration file." << endl;
       exit(-1);
     }
     axis_element = document->FindNextElement("axis");
   }
-  if (axisType == atNone) {
-    axisType = atLiftDrag;
+
+  if (forceAxisType == atNone) {
+    forceAxisType = atWind;
     cerr << endl << "  The aerodynamic axis system has been set by default"
                  << " to the Lift/Side/Drag system." << endl;
   }
+  if (momentAxisType == atNone) {
+    momentAxisType = atBodyXYZ;
+    cerr << endl << "  The aerodynamic moment axis system has been set by default"
+      << " to the bodyXYZ system." << endl;
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGAerodynamics::ProcessAxesNameAndFrame(eAxisType& axisType, const string& name,
+                                             const string& frame, Element* el,
+                                             const string& validNames)
+{
+  if (frame == "BODY" || frame.empty()) {
+    if (axisType == atNone) axisType = atBodyXYZ;
+    else if (axisType != atBodyXYZ)
+      cerr << endl << el->ReadFrom()
+           << endl << " Mixed aerodynamic axis systems have been used in the "
+                   << " aircraft config file." << validNames << " - BODY" << endl;
+  }
+  else if (frame == "STABILITY") {
+    if (axisType == atNone) axisType = atStability;
+    else if (axisType != atStability)
+      cerr << endl << el->ReadFrom()
+           << endl << " Mixed aerodynamic axis systems have been used in the "
+                   << " aircraft config file." << validNames << " - STABILITY" << endl;
+  }
+  else if (frame == "WIND") {
+    if (axisType == atNone) axisType = atWind;
+    else if (axisType != atWind)
+      cerr << endl << el->ReadFrom()
+           << endl << " Mixed aerodynamic axis systems have been used in the "
+                   << " aircraft config file." << validNames << " - WIND" << endl;
+  }
+  else {
+    cerr << endl << " Unknown axis frame type of - " << frame << endl;
+    exit(-1);
+  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -495,15 +574,24 @@ void FGAerodynamics::bind(void)
 {
   typedef double (FGAerodynamics::*PMF)(int) const;
 
-  PropertyManager->Tie("forces/fbx-aero-lbs",  this, 1, (PMF)&FGAerodynamics::GetForces);
-  PropertyManager->Tie("forces/fby-aero-lbs",  this, 2, (PMF)&FGAerodynamics::GetForces);
-  PropertyManager->Tie("forces/fbz-aero-lbs",  this, 3, (PMF)&FGAerodynamics::GetForces);
-  PropertyManager->Tie("moments/l-aero-lbsft", this, 1, (PMF)&FGAerodynamics::GetMoments);
-  PropertyManager->Tie("moments/m-aero-lbsft", this, 2, (PMF)&FGAerodynamics::GetMoments);
-  PropertyManager->Tie("moments/n-aero-lbsft", this, 3, (PMF)&FGAerodynamics::GetMoments);
-  PropertyManager->Tie("forces/fwx-aero-lbs",  this, 1, (PMF)&FGAerodynamics::GetvFw);
-  PropertyManager->Tie("forces/fwy-aero-lbs",  this, 2, (PMF)&FGAerodynamics::GetvFw);
-  PropertyManager->Tie("forces/fwz-aero-lbs",  this, 3, (PMF)&FGAerodynamics::GetvFw);
+  PropertyManager->Tie("forces/fbx-aero-lbs",  this, eX, (PMF)&FGAerodynamics::GetForces);
+  PropertyManager->Tie("forces/fby-aero-lbs",  this, eY, (PMF)&FGAerodynamics::GetForces);
+  PropertyManager->Tie("forces/fbz-aero-lbs",  this, eZ, (PMF)&FGAerodynamics::GetForces);
+  PropertyManager->Tie("moments/l-aero-lbsft", this, eL, (PMF)&FGAerodynamics::GetMoments);
+  PropertyManager->Tie("moments/m-aero-lbsft", this, eM, (PMF)&FGAerodynamics::GetMoments);
+  PropertyManager->Tie("moments/n-aero-lbsft", this, eN, (PMF)&FGAerodynamics::GetMoments);
+  PropertyManager->Tie("forces/fwx-aero-lbs",  this, eDrag, (PMF)&FGAerodynamics::GetvFw);
+  PropertyManager->Tie("forces/fwy-aero-lbs",  this, eSide, (PMF)&FGAerodynamics::GetvFw);
+  PropertyManager->Tie("forces/fwz-aero-lbs",  this, eLift, (PMF)&FGAerodynamics::GetvFw);
+  PropertyManager->Tie("forces/fsx-aero-lbs",  this, eX, (PMF)&FGAerodynamics::GetForcesInStabilityAxes);
+  PropertyManager->Tie("forces/fsy-aero-lbs",  this, eY, (PMF)&FGAerodynamics::GetForcesInStabilityAxes);
+  PropertyManager->Tie("forces/fsz-aero-lbs",  this, eZ, (PMF)&FGAerodynamics::GetForcesInStabilityAxes);
+  PropertyManager->Tie("moments/roll-stab-aero-lbsft", this, eRoll, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes);
+  PropertyManager->Tie("moments/pitch-stab-aero-lbsft", this, ePitch, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes);
+  PropertyManager->Tie("moments/yaw-stab-aero-lbsft", this, eYaw, (PMF)&FGAerodynamics::GetMomentsInStabilityAxes);
+  PropertyManager->Tie("moments/roll-wind-aero-lbsft", this, eRoll, (PMF)&FGAerodynamics::GetMomentsInWindAxes);
+  PropertyManager->Tie("moments/pitch-wind-aero-lbsft", this, ePitch, (PMF)&FGAerodynamics::GetMomentsInWindAxes);
+  PropertyManager->Tie("moments/yaw-wind-aero-lbsft", this, eYaw, (PMF)&FGAerodynamics::GetMomentsInWindAxes);
   PropertyManager->Tie("forces/lod-norm",      this, &FGAerodynamics::GetLoD);
   PropertyManager->Tie("aero/cl-squared",      this, &FGAerodynamics::GetClSquared);
   PropertyManager->Tie("aero/qbar-area", &qbar_area);
@@ -516,6 +604,44 @@ void FGAerodynamics::bind(void)
   PropertyManager->Tie("aero/stall-hyst-norm", this, &FGAerodynamics::GetHysteresisParm);
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//
+// Build transformation matrices for transforming from stability axes to 
+// body axes and to wind axes. Where "a" is alpha and "B" is beta:
+//
+// The transform from body to stability axes is:
+//
+//   cos(a)     0     sin(a)
+//   0          1     0
+//   -sin(a)    0     cos(a)
+//
+// The transform from stability to body axes is:
+//
+//   cos(a)     0     -sin(a)
+//   0          1     0
+//   sin(a)     0     cos(a)
+//
+//
+
+void FGAerodynamics::BuildStabilityTransformMatrices(void)
+{
+  double ca = cos(in.Alpha);
+  double sa = sin(in.Alpha);
+
+  // Stability-to-body
+  Ts2b(1, 1) = ca;
+  Ts2b(1, 2) = 0.0;
+  Ts2b(1, 3) = -sa;
+  Ts2b(2, 1) = 0.0;
+  Ts2b(2, 2) = 1.0;
+  Ts2b(2, 3) = 0.0;
+  Ts2b(3, 1) = sa;
+  Ts2b(3, 2) = 0.0;
+  Ts2b(3, 3) = ca;
+
+  Tb2s = Ts2b.Transposed();
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 //    The bitmasked value choices are as follows:
 //    unset: In this case (the default) JSBSim would only print
@@ -541,19 +667,22 @@ void FGAerodynamics::Debug(int from)
 
   if (debug_lvl & 1) { // Standard console startup message output
     if (from == 2) { // Loader
-      switch (axisType) {
-        case (atLiftDrag):
-          cout << endl << "  Aerodynamics (Lift|Side|Drag axes):" << endl << endl;
-          break;
-        case (atAxialNormal):
-          cout << endl << "  Aerodynamics (Axial|Side|Normal axes):" << endl << endl;
-          break;
-        case (atBodyXYZ):
-          cout << endl << "  Aerodynamics (X|Y|Z axes):" << endl << endl;
-          break;
+      switch (forceAxisType) {
+      case (atWind):
+        cout << endl << "  Aerodynamics (Lift|Side|Drag axes):" << endl << endl;
+        break;
+      case (atBodyAxialNormal):
+        cout << endl << "  Aerodynamics (Axial|Side|Normal axes):" << endl << endl;
+        break;
+      case (atBodyXYZ):
+        cout << endl << "  Aerodynamics (Body X|Y|Z axes):" << endl << endl;
+        break;
+      case (atStability):
+        cout << endl << "  Aerodynamics (Stability X|Y|Z axes):" << endl << endl;
+        break;
       case (atNone):
-          cout << endl << "  Aerodynamics (undefined axes):" << endl << endl;
-          break;
+        cout << endl << "  Aerodynamics (undefined axes):" << endl << endl;
+        break;
       }
     }
   }
diff --git a/src/FDM/JSBSim/models/FGAerodynamics.h b/src/FDM/JSBSim/models/FGAerodynamics.h
index 9b7ec1dad..f0adf6632 100644
--- a/src/FDM/JSBSim/models/FGAerodynamics.h
+++ b/src/FDM/JSBSim/models/FGAerodynamics.h
@@ -141,7 +141,7 @@ public:
       have found the aerodynamics keyword in the configuration file.
       @param element pointer to the current XML element for aerodynamics parameters.
       @return true if successful */
-  bool Load(Element* element);
+  virtual bool Load(Element* element);
 
   /** Gets the total aerodynamic force vector.
       @return a force vector reference. */
@@ -181,6 +181,34 @@ public:
       axis force. */
   double GetvFw(int axis) const { return vFw(axis); }
 
+  /** Retrieves the aerodynamic forces in the stability axes.
+  @return a reference to a column vector containing the stability axis forces. */
+  FGColumnVector3 GetForcesInStabilityAxes(void) const;
+
+  /** Retrieves the aerodynamic forces in the stability axes, given an axis.
+  @param axis the axis to return the force for (eX, eY, eZ).
+  @return a reference to a column vector containing the requested stability
+  axis force. */
+  double GetForcesInStabilityAxes(int n) const { return GetForcesInStabilityAxes()(n); }
+
+  /** Gets the total aerodynamic moment vector about the CG in the stability axes.
+  @return a moment vector reference. */
+  FGColumnVector3 GetMomentsInStabilityAxes(void) const { return Tb2s*vMoments; }
+
+  /** Gets the aerodynamic moment about the CG for an axis.
+  @return the moment about a single axis (as described also in the
+  similar call to GetForces(int n).*/
+  double GetMomentsInStabilityAxes(int n) const { return GetMomentsInStabilityAxes()(n); }
+
+  /** Gets the total aerodynamic moment vector about the CG in the wind axes.
+  @return a moment vector reference. */
+  FGColumnVector3 GetMomentsInWindAxes(void) const { return in.Tb2w*vMoments; }
+
+  /** Gets the aerodynamic moment about the CG for an axis.
+  @return the moment about a single axis (as described also in the
+  similar call to GetForces(int n).*/
+  double GetMomentsInWindAxes(int n) const { return GetMomentsInWindAxes()(n); }
+
   /** Retrieves the lift over drag ratio */
   double GetLoD(void) const { return lod; }
 
@@ -227,21 +255,22 @@ public:
   } in;
 
 private:
-  enum eAxisType {atNone, atLiftDrag, atAxialNormal, atBodyXYZ} axisType;
+  enum eAxisType {atNone, atWind, atBodyAxialNormal, atBodyXYZ, atStability} forceAxisType, momentAxisType;
   typedef std::map<std::string,int> AxisIndex;
   AxisIndex AxisIdx;
   FGFunction* AeroRPShift;
   typedef std::vector <FGFunction*> AeroFunctionArray;
   AeroFunctionArray* AeroFunctions;
+  FGMatrix33 Ts2b, Tb2s;
   FGColumnVector3 vFnative;
   FGColumnVector3 vFw;
   FGColumnVector3 vForces;
   AeroFunctionArray* AeroFunctionsAtCG;
-  FGColumnVector3 vFwAtCG;
   FGColumnVector3 vFnativeAtCG;
   FGColumnVector3 vForcesAtCG;
   FGColumnVector3 vMoments;
   FGColumnVector3 vMomentsMRC;
+  FGColumnVector3 vMomentsMRCBodyXYZ;
   FGColumnVector3 vDXYZcg;
   FGColumnVector3 vDeltaRP;
   double alphaclmax, alphaclmin;
@@ -253,7 +282,11 @@ private:
 
   typedef double (FGAerodynamics::*PMF)(int) const;
   void DetermineAxisSystem(Element* document);
+  void ProcessAxesNameAndFrame(FGAerodynamics::eAxisType& axisType,
+                               const string& name, const string& frame,
+                               Element* el, const string& validNames);
   void bind(void);
+  void BuildStabilityTransformMatrices(void);
 
   void Debug(int from);
 };
diff --git a/src/FDM/JSBSim/models/FGAircraft.cpp b/src/FDM/JSBSim/models/FGAircraft.cpp
index 66dcbfc40..36dff9737 100644
--- a/src/FDM/JSBSim/models/FGAircraft.cpp
+++ b/src/FDM/JSBSim/models/FGAircraft.cpp
@@ -78,7 +78,6 @@ FGAircraft::FGAircraft(FGFDMExec* fdmex) : FGModel(fdmex)
   lbarh = lbarv = 0.0;
   vbarh = vbarv = 0.0;
   WingIncidence = 0.0;
-  PitotAngle = 0.0;
 
   bind();
 
@@ -137,7 +136,7 @@ bool FGAircraft::Load(Element* el)
   string element_name;
   Element* element;
 
-  if (!FGModel::Load(el)) return false;
+  if (!FGModel::Load(el, true)) return false;
 
   if (el->FindElement("wingarea"))
     WingArea = el->FindElementValueAsNumberConvertTo("wingarea", "FT2");
@@ -155,8 +154,6 @@ bool FGAircraft::Load(Element* el)
     VTailArea = el->FindElementValueAsNumberConvertTo("vtailarea", "FT2");
   if (el->FindElement("vtailarm"))
     VTailArm = el->FindElementValueAsNumberConvertTo("vtailarm", "FT");
-  if (el->FindElement("pitot_angle"))
-    PitotAngle = el->FindElementValueAsNumberConvertTo("pitot_angle", "RAD");
 
   // Find all LOCATION elements that descend from this METRICS branch of the
   // config file. This would be CG location, eyepoint, etc.
diff --git a/src/FDM/JSBSim/models/FGAircraft.h b/src/FDM/JSBSim/models/FGAircraft.h
index 0726166ec..fa7f149b6 100644
--- a/src/FDM/JSBSim/models/FGAircraft.h
+++ b/src/FDM/JSBSim/models/FGAircraft.h
@@ -78,7 +78,6 @@ CLASS DOCUMENTATION
         <vtailarea unit="{FT2 | M}"> {number} </vtailarea>
         <vtailarm unit="{FT | M}"> {number} </vtailarm>
         <wing_incidence unit="{RAD | DEG}"> {number} </wing_incidence>
-        <pitot_angle unit="{RAD | DEG}"> {number} </pitot_angle>
         <location name="{AERORP | EYEPOINT | VRP}" unit="{IN | M}">
             <x> {number} </x>
             <y> {number} </y>
@@ -132,7 +131,7 @@ public:
       The executive calls this method to load the aircraft into JSBSim.
       @param el a pointer to the element tree
       @return true if successful */
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
 
   /** Gets the aircraft name
       @return the name of the aircraft as a string type */
@@ -144,7 +143,6 @@ public:
   double GetWingSpan(void) const { return WingSpan; }
   /// Gets the average wing chord
   double Getcbar(void) const { return cbar; }
-  double GetPitotAngle(void) const { return PitotAngle; }
   double GetWingIncidence(void) const { return WingIncidence; }
   double GetWingIncidenceDeg(void) const { return WingIncidence*radtodeg; }
   double GetHTailArea(void) const { return HTailArea; }
@@ -197,7 +195,7 @@ private:
 
   double WingArea, WingSpan, cbar, WingIncidence;
   double HTailArea, VTailArea, HTailArm, VTailArm;
-  double lbarh,lbarv,vbarh,vbarv,PitotAngle;
+  double lbarh,lbarv,vbarh,vbarv;
   std::string AircraftName;
 
   void Debug(int from);
diff --git a/src/FDM/JSBSim/models/FGAtmosphere.cpp b/src/FDM/JSBSim/models/FGAtmosphere.cpp
index b339057fa..bbac2744b 100644
--- a/src/FDM/JSBSim/models/FGAtmosphere.cpp
+++ b/src/FDM/JSBSim/models/FGAtmosphere.cpp
@@ -50,13 +50,24 @@ INCLUDES
 
 namespace JSBSim {
 
-IDENT(IdSrc,"$Id: FGAtmosphere.cpp,v 1.61 2016/01/10 19:22:12 bcoconni Exp $");
+IDENT(IdSrc,"$Id: FGAtmosphere.cpp,v 1.62 2016/01/16 12:05:47 bcoconni Exp $");
 IDENT(IdHdr,ID_ATMOSPHERE);
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
+// Atmosphere constants in British units converted from the SI values specified in the 
+// ISA document - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf
+
+const double KtoDegR = 1.8; // Kelvin to degree Rankine
+const double FGAtmosphere::Rstar = 8.31432 * (FGJSBBase::kgtoslug / (KtoDegR * FGJSBBase::fttom * FGJSBBase::fttom)); // ft*lbf/R/mol
+const double FGAtmosphere::Mair = 28.9645 * FGJSBBase::kgtoslug / 1000.0; // slug/mol
+const double FGAtmosphere::g0 = 9.80665 / FGJSBBase::fttom; // ft/s^2
+double FGAtmosphere::Reng = Rstar / Mair;
+
+const double FGAtmosphere::SHRatio = 1.40;
+
 FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex),
                                                PressureAltitude(0.0),      // ft
                                                DensityAltitude(0.0),       // ft
@@ -84,7 +95,7 @@ bool FGAtmosphere::InitModel(void)
 
   Calculate(0.0);
   SLtemperature = Temperature = 518.67;
-  SLpressure = Pressure = 2116.22;
+  SLpressure = Pressure = 2116.228;
   SLdensity = Density = Pressure/(Reng*Temperature);
   SLsoundspeed = Soundspeed = sqrt(SHRatio*Reng*(Temperature));
 
@@ -130,8 +141,8 @@ void FGAtmosphere::Calculate(double altitude)
     Density = node->GetDouble("atmosphere/override/density");
 
   Soundspeed  = sqrt(SHRatio*Reng*(Temperature));
-  PressureAltitude = altitude;
-  DensityAltitude = altitude;
+  PressureAltitude = CalculatePressureAltitude(Pressure, altitude);
+  DensityAltitude = CalculateDensityAltitude(Density, altitude);
 
   Viscosity = Beta * pow(Temperature, 1.5) / (SutherlandConstant + Temperature);
   KinematicViscosity = Viscosity / Density;
diff --git a/src/FDM/JSBSim/models/FGAtmosphere.h b/src/FDM/JSBSim/models/FGAtmosphere.h
index 4a603bdc7..632473aa2 100644
--- a/src/FDM/JSBSim/models/FGAtmosphere.h
+++ b/src/FDM/JSBSim/models/FGAtmosphere.h
@@ -232,6 +232,20 @@ protected:
   /// Calculate the atmosphere for the given altitude.
   void Calculate(double altitude);
 
+  /// Calculates the density altitude given any temperature or pressure bias.
+  /// Calculated density for the specified geometric altitude given any temperature
+  /// or pressure biases is passed in.
+  /// @param density
+  /// @param geometricAlt
+  virtual double CalculateDensityAltitude(double density, double geometricAlt) { return geometricAlt; }
+
+  /// Calculates the pressure altitude given any temperature or pressure bias.
+  /// Calculated pressure for the specified geometric altitude given any temperature
+  /// or pressure biases is passed in.
+  /// @param pressure
+  /// @param geometricAlt
+  virtual double CalculatePressureAltitude(double pressure, double geometricAlt) { return geometricAlt; }
+
   // Converts to Rankine from one of several unit systems.
   virtual double ConvertToRankine(double t, eTemperature unit) const;
   
@@ -241,6 +255,14 @@ protected:
   // Converts from PSF (pounds per square foot) to one of several unit systems.
   virtual double ConvertFromPSF(double t, ePressure unit=ePSF) const;
 
+  static const double Rstar;  // Universal gas constant - ft*lbf/R/mol
+  static const double Mair;   // Mean molecular weight - slug/mol
+  static const double g0;     // Sea-level acceleration of gravity - ft/s^2
+  static double Reng;         // Specific gas constant - ft*lbf/slug/R  
+
+  static const double SHRatio;
+  
+
   virtual void bind(void);
   void Debug(int from);
 };
diff --git a/src/FDM/JSBSim/models/FGAuxiliary.cpp b/src/FDM/JSBSim/models/FGAuxiliary.cpp
index b4eef4f31..ee5f2a6ca 100644
--- a/src/FDM/JSBSim/models/FGAuxiliary.cpp
+++ b/src/FDM/JSBSim/models/FGAuxiliary.cpp
@@ -68,10 +68,10 @@ FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex)
 
   vcas = veas = 0.0;
   qbar = qbarUW = qbarUV = 0.0;
-  Mach = MachU = MachPitot = 0.0;
+  Mach = MachU = 0.0;
   alpha = beta = 0.0;
   adot = bdot = 0.0;
-  gamma = Vt = Vground = Vpitot = 0.0;
+  gamma = Vt = Vground = 0.0;
   psigt = 0.0;
   day_of_year = 1;
   seconds_in_day = 0.0;
@@ -84,8 +84,6 @@ FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex)
   vAeroUVW.InitMatrix();
   vAeroPQR.InitMatrix();
   vMachUVW.InitMatrix();
-  vWindUVW.InitMatrix();
-  vPitotUVW.InitMatrix();
   vEuler.InitMatrix();
   vEulerRates.InitMatrix();
 
@@ -106,10 +104,10 @@ bool FGAuxiliary::InitModel(void)
 
   vcas = veas = 0.0;
   qbar = qbarUW = qbarUV = 0.0;
-  Mach = MachU = MachPitot = 0.0;
+  Mach = MachU = 0.0;
   alpha = beta = 0.0;
   adot = bdot = 0.0;
-  gamma = Vt = Vground = Vpitot = 0.0;
+  gamma = Vt = Vground = 0.0;
   psigt = 0.0;
   day_of_year = 1;
   seconds_in_day = 0.0;
@@ -206,20 +204,14 @@ bool FGAuxiliary::Run(bool Holding)
   tat = in.Temperature*(1 + 0.2*Mach*Mach); // Total Temperature, isentropic flow
   tatc = RankineToCelsius(tat);
 
-  // Pitot
+  pt = PitotTotalPressure(Mach, in.Pressure);
 
-  vWindUVW(eU) = Vt;
-  vPitotUVW = mTw2p * vWindUVW;
-  Vpitot = vPitotUVW(eU);
-  if (Vpitot < 0.0) Vpitot = 0.0;
-  MachPitot = Vpitot / in.SoundSpeed;
-  pt = PitotTotalPressure(MachPitot, in.Pressure);
-
-  if (abs(MachPitot) > 0.0) {
-    vcas = VcalibratedFromMach(MachPitot, in.Pressure, in.PressureSL, in.DensitySL);
+  if (abs(Mach) > 0.0) {
+    vcas = VcalibratedFromMach(Mach, in.Pressure, in.PressureSL, in.DensitySL);
     veas = sqrt(2 * qbar / in.DensitySL);
     vtrue = 1116.43559 * Mach * sqrt(in.Temperature / 518.67);
-  } else {
+  }
+  else {
     vcas = veas = vtrue = 0.0;
   }
 
@@ -284,23 +276,6 @@ void FGAuxiliary::UpdateWindMatrices(void)
   mTw2b(3,3) =  ca;
 
   mTb2w = mTw2b.Transposed();
-
-  // The pitot frame is the same as the body frame except rotated about the
-  // Y axis by the pitot attachment angle.
-
-  ca = cos(alpha + in.PitotAngle);
-  sa = sin(alpha + in.PitotAngle);
-
-  mTw2p(1,1) =  ca*cb;
-  mTw2p(1,2) = -ca*sb;
-  mTw2p(1,3) = -sa;
-  mTw2p(2,1) =  sb;
-  mTw2p(2,2) =  cb;
-  mTw2p(2,3) =  0.0;
-  mTw2p(3,1) =  sa*cb;
-  mTw2p(3,2) = -sa*sb;
-  mTw2p(3,3) =  ca;
-
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -424,6 +399,7 @@ void FGAuxiliary::bind(void)
   PropertyManager->Tie("aero/h_b-cg-ft", this, &FGAuxiliary::GetHOverBCG);
   PropertyManager->Tie("aero/h_b-mac-ft", this, &FGAuxiliary::GetHOverBMAC);
   PropertyManager->Tie("flight-path/gamma-rad", this, &FGAuxiliary::GetGamma);
+  PropertyManager->Tie("flight-path/gamma-deg", this, inDegrees, (PMF)&FGAuxiliary::GetGamma);
   PropertyManager->Tie("flight-path/psi-gt-rad", this, &FGAuxiliary::GetGroundTrack);
 
   PropertyManager->Tie("position/distance-from-start-lon-mt", this, &FGAuxiliary::GetLongitudeRelativePosition);
diff --git a/src/FDM/JSBSim/models/FGAuxiliary.h b/src/FDM/JSBSim/models/FGAuxiliary.h
index bb8c9f23f..c7f198bdf 100644
--- a/src/FDM/JSBSim/models/FGAuxiliary.h
+++ b/src/FDM/JSBSim/models/FGAuxiliary.h
@@ -236,6 +236,11 @@ public:
   double GetGamma(void)              const { return gamma;         }
   double GetGroundTrack(void)        const { return psigt;         }
 
+  double GetGamma(int unit) const {
+    if (unit == inDegrees) return gamma*radtodeg;
+    else return BadUnits();
+  }
+
   double GetHeadWind(void) const;
   double GetCrossWind(void) const;
 
@@ -289,7 +294,6 @@ public:
     FGColumnVector3 TurbPQR;
     double WindPsi;
     double Vwind;
-    double PitotAngle;
   } in;
 
 private:
@@ -298,7 +302,6 @@ private:
 
   FGMatrix33 mTw2b;
   FGMatrix33 mTb2w;
-  FGMatrix33 mTw2p;
 
   FGColumnVector3 vPilotAccel;
   FGColumnVector3 vPilotAccelN;
@@ -309,12 +312,10 @@ private:
   FGColumnVector3 vEuler;
   FGColumnVector3 vEulerRates;
   FGColumnVector3 vMachUVW;
-  FGColumnVector3 vWindUVW;
-  FGColumnVector3 vPitotUVW;
   FGLocation vLocationVRP;
 
-  double Vt, Vground, Vpitot;
-  double Mach, MachU, MachPitot;
+  double Vt, Vground;
+  double Mach, MachU;
   double qbar, qbarUW, qbarUV;
   double Re; // Reynolds Number = V*c/mu
   double alpha, beta;
diff --git a/src/FDM/JSBSim/models/FGBuoyantForces.cpp b/src/FDM/JSBSim/models/FGBuoyantForces.cpp
index 51ae27f01..ed3356aaa 100644
--- a/src/FDM/JSBSim/models/FGBuoyantForces.cpp
+++ b/src/FDM/JSBSim/models/FGBuoyantForces.cpp
@@ -123,7 +123,7 @@ bool FGBuoyantForces::Load(Element *document)
   Debug(2);
 
   // Perform base class Pre-Load
-  if (!FGModel::Load(document))
+  if (!FGModel::Load(document, true))
     return false;
 
   gas_cell_element = document->FindElement("gas_cell");
diff --git a/src/FDM/JSBSim/models/FGBuoyantForces.h b/src/FDM/JSBSim/models/FGBuoyantForces.h
index b63d90a84..e4a0ffee1 100644
--- a/src/FDM/JSBSim/models/FGBuoyantForces.h
+++ b/src/FDM/JSBSim/models/FGBuoyantForces.h
@@ -128,7 +128,7 @@ public:
       have found the Buoyant_forces keyword in the configuration file.
       @param element pointer to the current XML element for Buoyant forces parameters.
       @return true if successful */
-  bool Load(Element* element);
+  virtual bool Load(Element* element);
 
   /** Gets the total Buoyant force vector.
       @return a force vector in lbs. */
diff --git a/src/FDM/JSBSim/models/FGExternalReactions.cpp b/src/FDM/JSBSim/models/FGExternalReactions.cpp
index 50304959a..46574e412 100644
--- a/src/FDM/JSBSim/models/FGExternalReactions.cpp
+++ b/src/FDM/JSBSim/models/FGExternalReactions.cpp
@@ -72,7 +72,7 @@ FGExternalReactions::FGExternalReactions(FGFDMExec* fdmex) : FGModel(fdmex)
 bool FGExternalReactions::Load(Element* el)
 {
   // Call the base class Load() function to load interface properties.
-  if (!FGModel::Load(el))
+  if (!FGModel::Load(el, true))
     return false;
 
   Debug(2);
diff --git a/src/FDM/JSBSim/models/FGExternalReactions.h b/src/FDM/JSBSim/models/FGExternalReactions.h
index 50766967f..9e6d1160f 100644
--- a/src/FDM/JSBSim/models/FGExternalReactions.h
+++ b/src/FDM/JSBSim/models/FGExternalReactions.h
@@ -155,7 +155,7 @@ public:
       a FGExternalForce object will be instantiated for each force definition.
       @param el a pointer to the XML element holding the external reactions definition.
   */
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
 
   /** Retrieves the total forces defined in the external reactions.
       @return the total force in pounds.
diff --git a/src/FDM/JSBSim/models/FGFCS.cpp b/src/FDM/JSBSim/models/FGFCS.cpp
index cd8c28fe5..c9276c8f5 100644
--- a/src/FDM/JSBSim/models/FGFCS.cpp
+++ b/src/FDM/JSBSim/models/FGFCS.cpp
@@ -490,7 +490,7 @@ bool FGFCS::Load(Element* document)
   }
 
   // Load interface properties from document
-  if (!FGModel::Load(document))
+  if (!FGModel::Load(document, true))
     return false;
 
   Name += document->GetAttributeValue("name");
diff --git a/src/FDM/JSBSim/models/FGFCS.h b/src/FDM/JSBSim/models/FGFCS.h
index cca285672..941475b1d 100644
--- a/src/FDM/JSBSim/models/FGFCS.h
+++ b/src/FDM/JSBSim/models/FGFCS.h
@@ -542,7 +542,7 @@ public:
       Load() is called from FGFDMExec.
       @param el pointer to the Element instance
       @return true if succesful */
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
 
   SGPath FindFullPathName(const SGPath& path) const;
 
diff --git a/src/FDM/JSBSim/models/FGGroundReactions.cpp b/src/FDM/JSBSim/models/FGGroundReactions.cpp
index 1e8ae271c..0c5f48a5b 100644
--- a/src/FDM/JSBSim/models/FGGroundReactions.cpp
+++ b/src/FDM/JSBSim/models/FGGroundReactions.cpp
@@ -158,7 +158,7 @@ bool FGGroundReactions::Load(Element* document)
   Debug(2);
 
   // Perform base class Pre-Load
-  if (!FGModel::Load(document))
+  if (!FGModel::Load(document, true))
     return false;
 
   unsigned int numContacts = document->GetNumElements("contact");
diff --git a/src/FDM/JSBSim/models/FGGroundReactions.h b/src/FDM/JSBSim/models/FGGroundReactions.h
index fbb89e3a4..a7ceed967 100644
--- a/src/FDM/JSBSim/models/FGGroundReactions.h
+++ b/src/FDM/JSBSim/models/FGGroundReactions.h
@@ -93,7 +93,7 @@ public:
                      "Resume" command to be given.
       @return false if no error */
   bool Run(bool Holding);
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
   const FGColumnVector3& GetForces(void) const {return vForces;}
   double GetForces(int idx) const {return vForces(idx);}
   const FGColumnVector3& GetMoments(void) const {return vMoments;}
diff --git a/src/FDM/JSBSim/models/FGMassBalance.cpp b/src/FDM/JSBSim/models/FGMassBalance.cpp
index 33eb377f7..58fb622a8 100644
--- a/src/FDM/JSBSim/models/FGMassBalance.cpp
+++ b/src/FDM/JSBSim/models/FGMassBalance.cpp
@@ -134,7 +134,7 @@ bool FGMassBalance::Load(Element* document)
   Name = "Mass Properties Model: " + document->GetAttributeValue("name");
 
   // Perform base class Pre-Load
-  if (!FGModel::Load(document))
+  if (!FGModel::Load(document, true))
     return false;
 
   SetAircraftBaseInertias(ReadInertiaMatrix(document));
diff --git a/src/FDM/JSBSim/models/FGMassBalance.h b/src/FDM/JSBSim/models/FGMassBalance.h
index fba04d9df..b7c31176c 100644
--- a/src/FDM/JSBSim/models/FGMassBalance.h
+++ b/src/FDM/JSBSim/models/FGMassBalance.h
@@ -111,7 +111,7 @@ public:
   FGMassBalance(FGFDMExec*);
   ~FGMassBalance();
 
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
   bool InitModel(void);
   /** Runs the Mass Balance model; called by the Executive
       Can pass in a value indicating if the executive is directing the simulation to Hold.
diff --git a/src/FDM/JSBSim/models/FGModel.cpp b/src/FDM/JSBSim/models/FGModel.cpp
index c73bf7011..85f1a5c51 100644
--- a/src/FDM/JSBSim/models/FGModel.cpp
+++ b/src/FDM/JSBSim/models/FGModel.cpp
@@ -110,7 +110,7 @@ SGPath FGModel::FindFullPathName(const SGPath& path) const
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-bool FGModel::Load(Element* el)
+bool FGModel::Load(Element* el, bool preLoad)
 {
   FGModelLoader ModelLoader(this);
   Element* document = ModelLoader.Open(el);
@@ -124,17 +124,21 @@ bool FGModel::Load(Element* el)
     return false;
   }
 
-  bool result = FGModelFunctions::Load(document, PropertyManager);
+  bool result = true;
+
+  if (preLoad)
+    result = FGModelFunctions::Load(document, PropertyManager);
 
   if (document != el) {
     el->MergeAttributes(document);
 
-    // After reading interface properties in a file, read properties in the
-    // local model element. This allows general-purpose models to be defined in
-    // a file, with overrides or initial loaded constants supplied in the
-    // relevant element of the aircraft configuration file.
-
-    LocalProperties.Load(el, PropertyManager, true);
+    if (preLoad) {
+      // After reading interface properties in a file, read properties in the
+      // local model element. This allows general-purpose models to be defined
+      // in a file, with overrides or initial loaded constants supplied in the
+      // relevant element of the aircraft configuration file.
+      LocalProperties.Load(el, PropertyManager, true);
+    }
 
     Element* element = document->FindElement();
     while (element) {
diff --git a/src/FDM/JSBSim/models/FGModel.h b/src/FDM/JSBSim/models/FGModel.h
index a436e41c3..aa254ab6e 100644
--- a/src/FDM/JSBSim/models/FGModel.h
+++ b/src/FDM/JSBSim/models/FGModel.h
@@ -108,9 +108,11 @@ protected:
   unsigned int rate;
 
   /** Loads this model.
-      @param el a pointer to the element
+      @param el      a pointer to the element
+      @param preLoad true if model functions and local properties must be
+                     preloaded.
       @return true if model is successfully loaded*/
-  virtual bool Load(Element* el);
+  virtual bool Load(Element* el, bool preLoad);
 
   virtual void Debug(int from);
 
diff --git a/src/FDM/JSBSim/models/FGOutput.cpp b/src/FDM/JSBSim/models/FGOutput.cpp
index d7e15e8b2..88a200a44 100644
--- a/src/FDM/JSBSim/models/FGOutput.cpp
+++ b/src/FDM/JSBSim/models/FGOutput.cpp
@@ -44,10 +44,10 @@ INCLUDES
 #include "input_output/FGOutputSocket.h"
 #include "input_output/FGOutputTextFile.h"
 #include "input_output/FGOutputFG.h"
-#include "input_output/FGUDPOutputSocket.h"
 #include "input_output/FGXMLFileRead.h"
 #include "input_output/FGXMLElement.h"
 #include "input_output/FGModelLoader.h"
+#include "math/FGTemplateFunc.h"
 
 using namespace std;
 
@@ -76,9 +76,9 @@ FGOutput::FGOutput(FGFDMExec* fdmex) : FGModel(fdmex)
 
 FGOutput::~FGOutput()
 {
-  vector<FGOutputType*>::iterator it;
-  for (it = OutputTypes.begin(); it != OutputTypes.end(); ++it)
-    delete (*it);
+  vector<FGOutputType*>::iterator itv;
+  for (itv = OutputTypes.begin(); itv != OutputTypes.end(); ++itv)
+    delete (*itv);
 
   Debug(1);
 }
@@ -221,9 +221,6 @@ bool FGOutput::Load(int subSystems, std::string protocol, std::string type,
   } else if (type == "FLIGHTGEAR") {
     Output = new FGOutputFG(FDMExec);
     name += ":" + port + "/" + protocol;
-  } else if (type == "QTJSBSIM") {
-    Output = new FGUDPOutputSocket(FDMExec);
-    name += ":" + port + "/" + protocol;
   } else if (type == "TERMINAL") {
     // Not done yet
   } else if (type != string("NONE")) {
@@ -246,22 +243,30 @@ bool FGOutput::Load(int subSystems, std::string protocol, std::string type,
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-bool FGOutput::Load(Element* el)
+bool FGOutput::Load(Element* document, const SGPath& dir)
 {
-  // Unlike the other FGModel classes, properties listed in the <output> section
-  // are not intended to create new properties. For that reason, FGOutput
-  // cannot load its XML directives with FGModel::Load().
-  // Instead FGModelLoader::Open() and FGModel::PreLoad() must be explicitely
-  // called.
-  FGModelLoader ModelLoader(this);
-  Element* element = ModelLoader.Open(el);
+  // Optional path to use for included files
+  includePath = dir;
 
-  if (!element) return false;
+  // Perform base class Pre-Load
+  if (!FGModel::Load(document, false))
+    return false;
 
-  FGModel::PreLoad(element, PropertyManager);
+  Element *function = document->FindElement("function");
+
+  while (function) {
+    string fType = function->GetAttributeValue("type");
+
+    if (fType == "template") {
+      string name = function->GetAttributeValue("name");
+      TemplateFunctions[name] = new FGTemplateFunc(PropertyManager, function);
+    }
+
+    function = document->FindNextElement("function");
+  }
 
   size_t idx = OutputTypes.size();
-  string type = element->GetAttributeValue("type");
+  string type = document->GetAttributeValue("type");
   FGOutputType* Output = 0;
 
   if (debug_lvl > 0) cout << endl << "  Output data set: " << idx << "  " << endl;
@@ -276,8 +281,6 @@ bool FGOutput::Load(Element* el)
     Output = new FGOutputSocket(FDMExec);
   } else if (type == "FLIGHTGEAR") {
     Output = new FGOutputFG(FDMExec);
-  } else if (type == "QTJSBSIM") {
-    Output = new FGUDPOutputSocket(FDMExec);
   } else if (type == "TERMINAL") {
     // Not done yet
   } else if (type != string("NONE")) {
@@ -287,8 +290,9 @@ bool FGOutput::Load(Element* el)
   if (!Output) return false;
 
   Output->SetIdx(idx);
-  Output->Load(element);
-  PostLoad(element, PropertyManager);
+  Output->PreLoad(document, PropertyManager);
+  Output->Load(document);
+  Output->PostLoad(document, PropertyManager);
 
   OutputTypes.push_back(Output);
 
@@ -296,6 +300,19 @@ bool FGOutput::Load(Element* el)
   return true;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+SGPath FGOutput::FindFullPathName(const SGPath& path) const
+{
+  // Check optional include path if set
+  if (!includePath.isNull()) {
+    SGPath name = CheckPathName(includePath, path);
+    if (!name.isNull()) return name;
+  }
+
+  return FGModel::FindFullPathName(path);
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 //    The bitmasked value choices are as follows:
 //    unset: In this case (the default) JSBSim would only print
diff --git a/src/FDM/JSBSim/models/FGOutput.h b/src/FDM/JSBSim/models/FGOutput.h
index 9c7062247..b5b7396b7 100644
--- a/src/FDM/JSBSim/models/FGOutput.h
+++ b/src/FDM/JSBSim/models/FGOutput.h
@@ -42,6 +42,7 @@ INCLUDES
 
 #include "FGModel.h"
 #include "input_output/FGOutputType.h"
+#include "math/FGTemplateFunc.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
@@ -199,8 +200,9 @@ public:
   /** Load the output directives and adds a new output instance to the Output
       Manager list.
       @param el XMLElement that is pointing to the output directives
+      @param dir optional directory path to load included files from
       @result true if the execution succeeded. */
-  bool Load(Element* el);
+  virtual bool Load(Element* el, const SGPath& dir = SGPath());
   /** Load the output directives and adds a new output instance to the Output
       Manager list. Unlike the Load() method, the new output instance is not
       generated from output directives read in a XML file but from a list of
@@ -222,9 +224,20 @@ public:
       @result the name identifier.*/
   std::string GetOutputName(unsigned int idx) const;
 
+  SGPath FindFullPathName(const SGPath& path) const;
+
+  FGTemplateFunc* GetTemplateFunc(const std::string& name) {
+    if (TemplateFunctions.count(name))
+      return TemplateFunctions[name];
+    else
+      return NULL;
+  }
+
 private:
   std::vector<FGOutputType*> OutputTypes;
+  std::map<std::string, FGTemplateFunc_ptr> TemplateFunctions;
   bool enabled;
+  SGPath includePath;
 
   void Debug(int from);
 };
diff --git a/src/FDM/JSBSim/models/FGPropulsion.cpp b/src/FDM/JSBSim/models/FGPropulsion.cpp
index 561182ce3..ef9227b83 100644
--- a/src/FDM/JSBSim/models/FGPropulsion.cpp
+++ b/src/FDM/JSBSim/models/FGPropulsion.cpp
@@ -373,7 +373,7 @@ bool FGPropulsion::Load(Element* el)
   Name = "Propulsion Model: " + el->GetAttributeValue("name");
 
   // Perform base class Pre-Load
-  if (!FGModel::Load(el))
+  if (!FGModel::Load(el, true))
     return false;
 
   // Process tank definitions first to establish the number of fuel tanks
diff --git a/src/FDM/JSBSim/models/FGPropulsion.h b/src/FDM/JSBSim/models/FGPropulsion.h
index 728b7c916..b842ffdd4 100644
--- a/src/FDM/JSBSim/models/FGPropulsion.h
+++ b/src/FDM/JSBSim/models/FGPropulsion.h
@@ -128,7 +128,7 @@ public:
       Characteristics of the propulsion system are read in from the config file.
       @param el pointer to an XML element that contains the engine information.
       @return true if successfully loaded, otherwise false */
-  bool Load(Element* el);
+  virtual bool Load(Element* el);
 
   /// Retrieves the number of engines defined for the aircraft.
   unsigned int GetNumEngines(void) const {return (unsigned int)Engines.size();}
diff --git a/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp b/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp
index a5f6e8e54..0cc46a259 100644
--- a/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp
+++ b/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp
@@ -57,14 +57,16 @@ IDENT(IdHdr,ID_STANDARDATMOSPHERE);
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex) : FGAtmosphere(fdmex),
-                                                               TemperatureBias(0.0),
-                                                               TemperatureDeltaGradient(0.0)
+// Effective radius of the earth at a specific latitude per ISA 1976 (converted to ft)
+// r0 = 6356766 m
+const double FGStandardAtmosphere::EarthRadius = 6356766.0/FGJSBBase::fttom;
+
+FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex)
+  : FGAtmosphere(fdmex), TemperatureBias(0.0), TemperatureDeltaGradient(0.0),
+    StdAtmosTemperatureTable(9)
 {
   Name = "FGStandardAtmosphere";
 
-  StdAtmosTemperatureTable = new FGTable(9);
-
   // This is the U.S. Standard Atmosphere table for temperature in degrees
   // Rankine, based on geometric altitude. The table values are often given
   // in literature relative to geopotential altitude. 
@@ -72,36 +74,35 @@ FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex) : FGAtmosphere(fdme
   //                        GeoMet Alt    Temp      GeoPot Alt  GeoMet Alt
   //                           (ft)      (deg R)      (km)        (km)
   //                         --------   --------    ----------  ----------
-  //*StdAtmosTemperatureTable <<      0.00 << 518.67  //    0.000       0.000
-  //                          <<  36151.80 << 389.97  //   11.000      11.019
-  //                          <<  65823.90 << 389.97  //   20.000      20.063
-  //                          << 105518.06 << 411.60  //   32.000      32.162
-  //                          << 155348.07 << 487.20  //   47.000      47.350
-  //                          << 168676.12 << 487.20  //   51.000      51.413
-  //                          << 235570.77 << 386.40  //   71.000      71.802
-  //                          << 282152.08 << 336.50  //   84.852      86.000
-  //                          << 298556.40 << 336.50; //               91.000 - First layer in high altitude regime 
+  //StdAtmosTemperatureTable <<      0.00 << 518.67  //    0.000       0.000
+  //                         <<  36151.80 << 389.97  //   11.000      11.019
+  //                         <<  65823.90 << 389.97  //   20.000      20.063
+  //                         << 105518.06 << 411.60  //   32.000      32.162
+  //                         << 155348.07 << 487.20  //   47.000      47.350
+  //                         << 168676.12 << 487.20  //   51.000      51.413
+  //                         << 235570.77 << 386.40  //   71.000      71.802
+  //                         << 282152.08 << 336.50  //   84.852      86.000
+  //                         << 298556.40 << 336.50; //               91.000 - First layer in high altitude regime 
 
   //                            GeoPot Alt    Temp       GeoPot Alt  GeoMet Alt
   //                               (ft)      (deg R)        (km)        (km)
   //                           -----------   --------     ----------  ----------
-  *StdAtmosTemperatureTable <<      0.0000 << 518.67  //    0.000       0.000
-                            <<  36089.2388 << 389.97  //   11.000      11.019
-                            <<  65616.7979 << 389.97  //   20.000      20.063
-                            << 104986.8766 << 411.57  //   32.000      32.162
-                            << 154199.4751 << 487.17  //   47.000      47.350
-                            << 167322.8346 << 487.17  //   51.000      51.413
-                            << 232939.6325 << 386.37  //   71.000      71.802
-                            << 278385.8268 << 336.50 //   84.852      86.000
-                            << 298556.40   << 336.50; //               91.000 - First layer in high altitude regime 
+  StdAtmosTemperatureTable <<      0.0000 << 518.67  //    0.000       0.000
+                           <<  36089.2388 << 389.97  //   11.000      11.019
+                           <<  65616.7979 << 389.97  //   20.000      20.063
+                           << 104986.8766 << 411.57  //   32.000      32.162
+                           << 154199.4751 << 487.17  //   47.000      47.350
+                           << 167322.8346 << 487.17  //   51.000      51.413
+                           << 232939.6325 << 386.37  //   71.000      71.802
+                           << 278385.8268 << 336.5028  // 84.852      86.000
+                           << 298556.4304 << 336.5028; //             91.000 - First layer in high altitude regime 
 
-  LapseRateVector.resize(StdAtmosTemperatureTable->GetNumRows()-1);
-  PressureBreakpointVector.resize(StdAtmosTemperatureTable->GetNumRows());
+  PressureBreakpointVector.resize(StdAtmosTemperatureTable.GetNumRows());
 
   // Assume the altitude to fade out the gradient at is at the highest
   // altitude in the table. Above that, other functions are used to
   // calculate temperature.
-  GradientFadeoutAltitude = (*StdAtmosTemperatureTable)(StdAtmosTemperatureTable->GetNumRows(),0);
+  GradientFadeoutAltitude = StdAtmosTemperatureTable(StdAtmosTemperatureTable.GetNumRows(),0);
 
   bind();
   Debug(0);
@@ -111,8 +112,6 @@ FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex) : FGAtmosphere(fdme
 
 FGStandardAtmosphere::~FGStandardAtmosphere()
 {
-  delete StdAtmosTemperatureTable;
-  LapseRateVector.clear();
   Debug(1);
 }
 
@@ -120,16 +119,22 @@ FGStandardAtmosphere::~FGStandardAtmosphere()
 
 bool FGStandardAtmosphere::InitModel(void)
 {
-  PressureBreakpointVector[0] = StdSLpressure = 2116.22; // psf
+  PressureBreakpointVector[0] = StdSLpressure = SLpressure = Pressure = 2116.228; // psf
   TemperatureDeltaGradient = 0.0;
   TemperatureBias = 0.0;
   CalculateLapseRates();
   CalculatePressureBreakpoints();
+
+  StdSLtemperature = SLtemperature = StdAtmosTemperatureTable(1, 1);
+  StdSLdensity     = SLdensity = StdSLpressure / (Reng * StdSLtemperature);
+
+  StdPressureBreakpointVector = PressureBreakpointVector;
+
+  CalculateStdDensityBreakpoints();
+
   Calculate(0.0);
-  StdSLtemperature = SLtemperature = Temperature;
-  SLpressure = Pressure;
-  StdSLdensity     = SLdensity = Density;
-  StdSLsoundspeed  = SLsoundspeed = Soundspeed;
+
+  StdSLsoundspeed = SLsoundspeed = Soundspeed;
 
   rSLtemperature = 1/SLtemperature ;
   rSLpressure    = 1/SLpressure    ;
@@ -143,43 +148,39 @@ bool FGStandardAtmosphere::InitModel(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Get the actual pressure as modeled at a specified altitude
-// These calculations are from equations 33a and 33b in the U.S. Standard Atmosphere
-// document referenced in the documentation for this code.
+// These calculations are from equations 33a and 33b in the U.S. Standard
+// Atmosphere document referenced in the documentation for this code.
 
 double FGStandardAtmosphere::GetPressure(double altitude) const
 {
-  unsigned int b=0;
-  double pressure = 0.0;
-  double Lmb, Exp, Tmb, deltaH, factor;
-  double numRows = StdAtmosTemperatureTable->GetNumRows();
+  double GeoPotAlt = GeopotentialAltitude(altitude);
 
   // Iterate through the altitudes to find the current Base Altitude
   // in the table. That is, if the current altitude (the argument passed in)
   // is 20000 ft, then the base altitude from the table is 0.0. If the
   // passed-in altitude is 40000 ft, the base altitude is 36089.2388 ft (and
   // the index "b" is 2 - the second entry in the table).
-  double testAlt = (*StdAtmosTemperatureTable)(b+1,0);
-  double GeoPotAlt = (altitude*20855531.5)/(20855531.5+altitude);
-  while ((GeoPotAlt >= testAlt) && (b <= numRows-2)) {
-    b++;
-    testAlt = (*StdAtmosTemperatureTable)(b+1,0);
-  }
-  if (b>0) b--;
+  double BaseAlt = StdAtmosTemperatureTable(1,0);
+  unsigned int numRows = StdAtmosTemperatureTable.GetNumRows();
+  unsigned int b;
 
-  double BaseAlt = (*StdAtmosTemperatureTable)(b+1,0);
-  Tmb = GetTemperature(BaseAlt);
-  deltaH = GeoPotAlt - BaseAlt;
-
-  if (LapseRateVector[b] != 0.00) {
-    Lmb = LapseRateVector[b];
-    Exp = Mair/(Rstar*Lmb);
-    factor = Tmb/(Tmb + Lmb*deltaH);
-    pressure = PressureBreakpointVector[b]*pow(factor, Exp);
-  } else {
-    pressure = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
+  for (b=0; b < numRows-2; ++b) {
+    double testAlt = StdAtmosTemperatureTable(b+2,0);
+    if (GeoPotAlt < testAlt)
+      break;
+    BaseAlt = testAlt;
   }
 
-  return pressure;
+  double Tmb = GetTemperature(GeometricAltitude(BaseAlt));
+  double deltaH = GeoPotAlt - BaseAlt;
+  double Lmb = LapseRateVector[b];
+
+  if (Lmb != 0.0) {
+    double Exp = g0*Mair / (Rstar*Lmb);
+    double factor = Tmb/(Tmb + Lmb*deltaH);
+    return PressureBreakpointVector[b]*pow(factor, Exp);
+  } else
+    return PressureBreakpointVector[b]*exp(-g0*Mair*deltaH/(Rstar*Tmb));
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -198,11 +199,23 @@ void FGStandardAtmosphere::SetPressureSL(ePressure unit, double pressure)
 
 double FGStandardAtmosphere::GetTemperature(double altitude) const
 {
-  double GeoPotAlt = (altitude*20855531.5)/(20855531.5+altitude);
+  double GeoPotAlt = GeopotentialAltitude(altitude);
 
-  double T = StdAtmosTemperatureTable->GetValue(GeoPotAlt) + TemperatureBias;
-  if (altitude <= GradientFadeoutAltitude)
-    T += TemperatureDeltaGradient * (GradientFadeoutAltitude - altitude);
+  double T;
+
+  if (GeoPotAlt >= 0.0) {
+    T = StdAtmosTemperatureTable.GetValue(GeoPotAlt);
+
+    if (GeoPotAlt <= GradientFadeoutAltitude)
+      T -= TemperatureDeltaGradient * GeoPotAlt;
+  }
+  else {
+    // We don't need to add TemperatureDeltaGradient*GeoPotAlt here because
+    // the lapse rate vector already accounts for the temperature gradient.
+    T = StdAtmosTemperatureTable.GetValue(0.0) + GeoPotAlt*LapseRateVector[0];
+  }
+
+  T += TemperatureBias + TemperatureDeltaGradient * GradientFadeoutAltitude;
 
   return T;
 }
@@ -218,8 +231,12 @@ double FGStandardAtmosphere::GetStdTemperature(double altitude) const
 
   if (altitude < 298556.4) {                // 91 km - station 8
 
-    double GeoPotAlt = (altitude*20855531.5)/(20855531.5+altitude);
-    temp = StdAtmosTemperatureTable->GetValue(GeoPotAlt);
+    double GeoPotAlt = GeopotentialAltitude(altitude);
+
+    if (GeoPotAlt >= 0.0)
+      temp = StdAtmosTemperatureTable.GetValue(GeoPotAlt);
+    else
+      temp = StdAtmosTemperatureTable.GetValue(0.0) + GeoPotAlt*LapseRateVector[0];
 
   } else if (altitude < 360892.4) {        // 110 km - station 9
 
@@ -244,42 +261,34 @@ double FGStandardAtmosphere::GetStdTemperature(double altitude) const
 
 double FGStandardAtmosphere::GetStdPressure(double altitude) const
 {
-  double press=0;
-  if (TemperatureBias == 0.0 && TemperatureDeltaGradient == 0.0 && PressureBreakpointVector[0] == StdSLpressure) {
-    press = GetPressure(altitude);
-  } else if (altitude <= 100000.0) {
-    GetStdPressure100K(altitude);
-  } else {
-    // Cannot currently retrieve the standard pressure
+  double GeoPotAlt = GeopotentialAltitude(altitude);
+
+  // Iterate through the altitudes to find the current Base Altitude
+  // in the table. That is, if the current altitude (the argument passed in)
+  // is 20000 ft, then the base altitude from the table is 0.0. If the
+  // passed-in altitude is 40000 ft, the base altitude is 36089.2388 ft (and
+  // the index "b" is 2 - the second entry in the table).
+  double BaseAlt = StdAtmosTemperatureTable(1,0);
+  unsigned int numRows = StdAtmosTemperatureTable.GetNumRows();
+  unsigned int b;
+
+  for (b=0; b < numRows-2; ++b) {
+    double testAlt = StdAtmosTemperatureTable(b+2,0);
+    if (GeoPotAlt < testAlt)
+      break;
+    BaseAlt = testAlt;
   }
-  return press;
-}
 
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-// This function calculates an approximation of the standard atmospheric pressure
-// up to an altitude of about 100,000 ft. If the temperature and pressure are not
-// altered for local conditions, the GetPressure(h) function should be used,
-// as that is valid to a much higher altitude. This function is accurate to within
-// a couple of psf up to 100K ft. This polynomial fit was determined using Excel.
+  double Tmb = GetStdTemperature(GeometricAltitude(BaseAlt));
+  double deltaH = GeoPotAlt - BaseAlt;
+  double Lmb = LapseRateVector[b];
 
-double FGStandardAtmosphere::GetStdPressure100K(double altitude) const
-{
-  // Limit this equation to input altitudes of 100000 ft.
-  if (altitude > 100000.0) altitude = 100000.0;
-
-  double alt[5];
-  const double coef[5] = {  2116.217,
-                          -7.648932746E-2,
-                           1.0925498604E-6,
-                          -7.1135726027E-12,
-                           1.7470331356E-17 };
-
-  alt[0] = 1;
-  for (int pwr=1; pwr<=4; pwr++) alt[pwr] = alt[pwr-1]*altitude;
-
-  double press = 0.0;
-  for (int ctr=0; ctr<=4; ctr++) press += coef[ctr]*alt[ctr];
-  return press;
+  if (Lmb != 0.0) {
+    double Exp = g0*Mair / (Rstar*Lmb);
+    double factor = Tmb/(Tmb + Lmb*deltaH);
+    return StdPressureBreakpointVector[b]*pow(factor, Exp);
+  } else
+    return StdPressureBreakpointVector[b]*exp(-g0*Mair*deltaH/(Rstar*Tmb));
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -343,7 +352,7 @@ void FGStandardAtmosphere::SetTemperatureGradedDelta(double deltemp, double h, e
   if (unit == eCelsius || unit == eKelvin)
     deltemp *= 1.80; // If temp delta "t" is given in metric, scale up to English
 
-  TemperatureDeltaGradient = deltemp/(GradientFadeoutAltitude - h);
+  TemperatureDeltaGradient = deltemp/(GradientFadeoutAltitude - GeopotentialAltitude(h));
   CalculateLapseRates();
   CalculatePressureBreakpoints();
 }
@@ -376,13 +385,16 @@ void FGStandardAtmosphere::SetTemperatureGradedDelta(double deltemp, double h, e
 
 void FGStandardAtmosphere::CalculateLapseRates()
 {
-  for (unsigned int bh=0; bh<LapseRateVector.size(); bh++)
+  unsigned int numRows = StdAtmosTemperatureTable.GetNumRows();
+  LapseRateVector.clear();
+
+  for (unsigned int bh=0; bh < numRows-1; bh++)
   {
-    double t0 = (*StdAtmosTemperatureTable)(bh+1,1);
-    double t1 = (*StdAtmosTemperatureTable)(bh+2,1);
-    double h0 = (*StdAtmosTemperatureTable)(bh+1,0);
-    double h1 = (*StdAtmosTemperatureTable)(bh+2,0);
-    LapseRateVector[bh] = (t1 - t0) / (h1 - h0) + TemperatureDeltaGradient;
+    double t0 = StdAtmosTemperatureTable(bh+1,1);
+    double t1 = StdAtmosTemperatureTable(bh+2,1);
+    double h0 = StdAtmosTemperatureTable(bh+1,0);
+    double h1 = StdAtmosTemperatureTable(bh+2,0);
+    LapseRateVector.push_back((t1 - t0) / (h1 - h0) - TemperatureDeltaGradient);
   }
 }
 
@@ -391,20 +403,20 @@ void FGStandardAtmosphere::CalculateLapseRates()
 void FGStandardAtmosphere::CalculatePressureBreakpoints()
 {
   for (unsigned int b=0; b<PressureBreakpointVector.size()-1; b++) {
-    double BaseTemp = (*StdAtmosTemperatureTable)(b+1,1);
-    double BaseAlt = (*StdAtmosTemperatureTable)(b+1,0);
-    double UpperAlt = (*StdAtmosTemperatureTable)(b+2,0);
+    double BaseTemp = StdAtmosTemperatureTable(b+1,1);
+    double BaseAlt = StdAtmosTemperatureTable(b+1,0);
+    double UpperAlt = StdAtmosTemperatureTable(b+2,0);
     double deltaH = UpperAlt - BaseAlt;
     double Tmb = BaseTemp
                  + TemperatureBias 
                  + (GradientFadeoutAltitude - BaseAlt)*TemperatureDeltaGradient;
     if (LapseRateVector[b] != 0.00) {
       double Lmb = LapseRateVector[b];
-      double Exp = Mair/(Rstar*Lmb);
+      double Exp = g0*Mair / (Rstar*Lmb);
       double factor = Tmb/(Tmb + Lmb*deltaH);
       PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*pow(factor, Exp);
     } else {
-      PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
+      PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*exp(-g0*Mair*deltaH/(Rstar*Tmb));
     }
   }
 }
@@ -428,6 +440,85 @@ void FGStandardAtmosphere::ResetSLPressure()
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+void FGStandardAtmosphere::CalculateStdDensityBreakpoints()
+{
+  StdDensityBreakpointVector.clear();
+  for (unsigned int i = 0; i < StdPressureBreakpointVector.size(); i++) {
+    StdDensityBreakpointVector.push_back(StdPressureBreakpointVector[i] / (Reng * StdAtmosTemperatureTable(i + 1, 1)));
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGStandardAtmosphere::CalculateDensityAltitude(double density, double geometricAlt)
+{
+  // Work out which layer we're dealing with
+  unsigned int b = 0;
+  for (; b < StdDensityBreakpointVector.size() - 2; b++) {
+    if (density >= StdDensityBreakpointVector[b + 1])
+      break;
+  }
+
+  // Get layer properties
+  double Tmb = StdAtmosTemperatureTable(b + 1, 1);
+  double Hb = StdAtmosTemperatureTable(b + 1, 0);
+  double UpperTemp = StdAtmosTemperatureTable(b + 2, 1);
+  double UpperAlt = StdAtmosTemperatureTable(b + 2, 0);
+  double deltaH = UpperAlt - Hb;
+  double Lmb = (UpperTemp - Tmb) / deltaH;
+  double pb = StdDensityBreakpointVector[b];
+
+  double density_altitude = 0.0;
+
+  // https://en.wikipedia.org/wiki/Barometric_formula for density solved for H
+  if (Lmb != 0.0) {
+    double Exp = -1.0 / (1.0 + (g0*Mair)/(Rstar*Lmb));
+    density_altitude = Hb + (Tmb / Lmb) * (pow(density / pb, Exp) - 1);
+  } else {
+    double Factor = -(Rstar*Tmb) / (g0*Mair);
+    density_altitude = Hb + Factor * log(density / pb);
+  }
+
+  return GeometricAltitude(density_altitude);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGStandardAtmosphere::CalculatePressureAltitude(double pressure, double geometricAlt)
+{
+  // Work out which layer we're dealing with
+  unsigned int b = 0;
+  for (; b < StdPressureBreakpointVector.size() - 2; b++) {
+    if (pressure >= StdPressureBreakpointVector[b + 1])
+      break;
+  }
+
+  // Get layer properties
+  double Tmb = StdAtmosTemperatureTable(b + 1, 1);
+  double Hb = StdAtmosTemperatureTable(b + 1, 0);
+  double UpperTemp = StdAtmosTemperatureTable(b + 2, 1);
+  double UpperAlt = StdAtmosTemperatureTable(b + 2, 0);
+  double deltaH = UpperAlt - Hb;
+  double Lmb = (UpperTemp - Tmb) / deltaH;
+  double Pb = StdPressureBreakpointVector[b];
+
+  double pressure_altitude = 0.0;
+
+  if (Lmb != 0.00) {
+    // Equation 33(a) from ISA document solved for H
+    double Exp = -(Rstar*Lmb) / (g0*Mair);
+    pressure_altitude = Hb + (Tmb / Lmb) * (pow(pressure / Pb, Exp) - 1);
+  } else {
+    // Equation 33(b) from ISA document solved for H
+    double Factor = -(Rstar*Tmb) / (g0*Mair);
+    pressure_altitude = Hb + Factor * log(pressure / Pb);
+  }
+
+  return GeometricAltitude(pressure_altitude);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 void FGStandardAtmosphere::bind(void)
 {
   typedef double (FGStandardAtmosphere::*PMFi)(int) const;
diff --git a/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.h b/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.h
index ec2e24aba..0c5fef10a 100644
--- a/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.h
+++ b/src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.h
@@ -217,9 +217,6 @@ public:
   /// Returns the pressure at a specified altitude in psf.
   virtual double GetPressure(double altitude) const;
 
-  /// Returns the standard pressure at a specified altitude in psf
-  virtual double GetStdPressure100K(double altitude) const;
-
   /// Returns the standard pressure at the specified altitude.
   virtual double GetStdPressure(double altitude) const;
 
@@ -253,9 +250,11 @@ protected:
   double TemperatureDeltaGradient;
   double GradientFadeoutAltitude;
 
-  FGTable* StdAtmosTemperatureTable;
+  FGTable StdAtmosTemperatureTable;
   std::vector<double> LapseRateVector;
   std::vector<double> PressureBreakpointVector;
+  std::vector<double> StdPressureBreakpointVector;
+  std::vector<double> StdDensityBreakpointVector;
 
   /// Recalculate the lapse rate vectors when the temperature profile is altered
   /// in a way that would change the lapse rates, such as when a gradient is applied.
@@ -266,8 +265,42 @@ protected:
   /// altitudes in the standard temperature table.
   void CalculatePressureBreakpoints();
 
+  /// Calculate the atmospheric density breakpoints at the 
+  /// altitudes in the standard temperature table.
+  void CalculateStdDensityBreakpoints();
+
+  /// Convert a geometric altitude to a geopotential altitude
+  double GeopotentialAltitude(double geometalt) const { return (geometalt * EarthRadius) / (EarthRadius + geometalt); }
+
+  /// Convert a geopotential altitude to a geometric altitude
+  double GeometricAltitude(double geopotalt) const { return (geopotalt * EarthRadius) / (EarthRadius - geopotalt); }
+
+  /** Calculates the density altitude given any temperature or pressure bias.
+  Calculated density for the specified geometric altitude given any temperature
+  or pressure biases is passed in.
+  @param density
+  @param geometricAlt
+  @see
+  https://en.wikipedia.org/wiki/Density_altitude
+  https://wahiduddin.net/calc/density_altitude.htm
+  */
+  virtual double CalculateDensityAltitude(double density, double geometricAlt);
+
+  /** Calculates the pressure altitude given any temperature or pressure bias.
+  Calculated density for the specified geometric altitude given any temperature
+  or pressure biases is passed in.
+  @param pressure
+  @param geometricAlt
+  @see
+  https://en.wikipedia.org/wiki/Pressure_altitude
+  */
+  virtual double CalculatePressureAltitude(double pressure, double geometricAlt);
+
   virtual void bind(void);
   void Debug(int from);
+
+  /// Earth radius in ft as defined for ISA 1976
+  static const double EarthRadius;
 };
 
 } // namespace JSBSim
diff --git a/src/FDM/JSBSim/models/flight_control/FGDeadBand.h b/src/FDM/JSBSim/models/flight_control/FGDeadBand.h
index d1aab83e1..5162f880c 100644
--- a/src/FDM/JSBSim/models/flight_control/FGDeadBand.h
+++ b/src/FDM/JSBSim/models/flight_control/FGDeadBand.h
@@ -43,7 +43,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_DEADBAND "$Id: FGDeadBand.h,v 1.10 2013/01/26 17:06:50 bcoconni Exp $"
+#define ID_DEADBAND "$Id$"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -79,7 +79,7 @@ CLASS DOCUMENTATION
     produce no output. For example, say that the width value is 2.0. If the
     input is between -1.0 and +1.0, the output will be zero.
     @author Jon S. Berndt
-    @version $Id: FGDeadBand.h,v 1.10 2013/01/26 17:06:50 bcoconni Exp $
+    @version $Id$
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/flight_control/FGGain.h b/src/FDM/JSBSim/models/flight_control/FGGain.h
index 2105ff2a1..9726fb0ad 100644
--- a/src/FDM/JSBSim/models/flight_control/FGGain.h
+++ b/src/FDM/JSBSim/models/flight_control/FGGain.h
@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_GAIN "$Id: FGGain.h,v 1.15 2013/01/26 17:06:50 bcoconni Exp $"
+#define ID_GAIN "$Id$"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -208,7 +208,7 @@ CLASS DOCUMENTATION
     @endcode
 
     @author Jon S. Berndt
-    @version $Revision: 1.15 $
+    @version $Revision$
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/propulsion/FGPropeller.h b/src/FDM/JSBSim/models/propulsion/FGPropeller.h
index 57c442c50..720e05058 100644
--- a/src/FDM/JSBSim/models/propulsion/FGPropeller.h
+++ b/src/FDM/JSBSim/models/propulsion/FGPropeller.h
@@ -64,6 +64,7 @@ CLASS DOCUMENTATION
 
 ~~~{.xml}
 <sense> {1 | -1} </sense> 
+<p_factor> {number} </p_factor>
 <propeller name="{string}" version="{string}">
   <ixx> {number} </ixx>
   <diameter unit="IN"> {number} </diameter>
@@ -75,7 +76,6 @@ CLASS DOCUMENTATION
   <maxrpm> {number} </maxrpm>
   <constspeed> {number} </constspeed>
   <reversepitch> {number} </reversepitch>
-  <p_factor> {number} </p_factor>
   <ct_factor> {number} </ct_factor>
   <cp_factor> {number} </cp_factor>
 
@@ -122,7 +122,8 @@ CLASS DOCUMENTATION
     \<sense>         - Direction of rotation (1=clockwise as viewed from cockpit,
                         -1=anti-clockwise as viewed from cockpit). Sense is
                        specified in the parent tag of the propeller.
-    \<p_factor>      - P factor.
+    \<p_factor>      - P factor. It is specified in the parent tag of
+                       the propeller.
     \<ct_factor>     - A multiplier for the coefficients of thrust.
     \<cp_factor>     - A multiplier for the coefficients of power.
 </pre>
@@ -133,6 +134,11 @@ coefficient of power (Cp).
 Two tables are optional. They apply a factor to Ct and Cp based on the
 helical tip Mach.
 
+The parameters <sense> and <p_factor> must be specified at the parent level i.e.
+in the <thruster> element. This allows to specify different sense and P factor
+values for each propeller of the model while using the same definition file for
+all the propellers.
+
 In addition to thrust, the propeller applies two moments to the aircraft:
 - The torque that tends to roll the aircraft in the direction opposite to the
 propeller rotation,