GetSatVapPres Function

  function GetSatVapPres(TDryBulb) result(SatVapPres)
    !+ Return saturation vapor pressure given dry-bulb temperature.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1  eqn 5
    !+ Important note: the ASHRAE formulae are defined above and below the freezing point but have
    !+ a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae
    !+ should be defined above and below the triple point of water (not the feezing point) in which case 
    !+ the discontinuity vanishes. It is essential to use the triple point of water otherwise function
    !+ GetTDewPointFromVapPres, which inverts the present function, does not converge properly around
    !+ the freezing point.

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real              ::  SatVapPres
      !+ Vapor pressure of saturated air in Psi [IP] or Pa [SI]
    real              ::  LnPws
      !+ Log of Vapor Pressure of saturated air (dimensionless)
    real              ::  T
      !+ Dry bulb temperature in R [IP] or K [SI]

    if (isIP()) then
      if (TDryBulb < -148.0 .or. TDryBulb > 392.0) then
        error stop "Error: dry bulb temperature must be in range [-148, 392]°F"
      end if

      T = GetTRankineFromTFahrenheit(TDryBulb)

      if (TDryBulb <= TRIPLE_POINT_WATER_IP) then
        LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T**2    &
                + 3.5575832E-10 * T**3 - 9.0344688E-14 * T**4 + 4.1635019 * log(T))
      else
        LnPws = -1.0440397E+04 / T - 1.1294650E+01 - 2.7022355E-02* T + 1.2890360E-05 * T**2  &
                - 2.4780681E-09 * T**3 + 6.5459673 * log(T)
      end if

      else
        if (TDryBulb < -100.0 .or. TDryBulb > 200.0) then
          error stop "Error: dry bulb temperature must be in range [-100, 200]°C"
        end if

        T = GetTKelvinFromTCelsius(TDryBulb)

        if (TDryBulb <= TRIPLE_POINT_WATER_SI) then
          LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T**2    &
                  + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T)
        else
          LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T**2   &
                  - 1.4452093E-08 * T**3 + 6.5459673 * log(T)
        end if
      end if

    SatVapPres = exp(LnPws)
  end function GetSatVapPres