GetTDewPointFromVapPres Function

  function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint)
    !+ Return dew-point temperature given dry-bulb temperature and vapor pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 and 6
    !+ Notes:
    !+ The dew point temperature is solved by inverting the equation giving water vapor pressure
    !+ at saturation from temperature rather than using the regressions provided
    !+ by ASHRAE (eqn. 37 and 38) which are much less accurate and have a
    !+ narrower range of validity.
    !+ The Newton-Raphson (NR) method is used on the logarithm of water vapour
    !+ pressure as a function of temperature, which is a very smooth function
    !+ Convergence is usually achieved in 3 to 5 iterations.
    !+ TDryBulb is not really needed here, just used for convenience.

    real, intent(in)    ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real                ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real                ::  lnVP
      !+ Natural logarithm of partial pressure of water vapor pressure in moist air
    real                ::  d_lnVP
      !+ Derivative of function, calculated numerically
    real                ::  lnVP_iter
      !+ Value of log of vapor water pressure used in NR calculation
    real                ::  TDewPoint_iter
      !+ Value of TDewPoint used in NR calculation
    real, dimension(2)  ::  BOUNDS
      !+ Valid temperature range in °F [IP] or °C [SI]
    integer             :: index
      !+ Index used in the calculation

    ! Bounds and step size as a function of the system of units
    if (isIP()) then
        BOUNDS(1) = -148.0
        BOUNDS(2) =  392.0
    else
        BOUNDS(1) = -100.0
        BOUNDS(2) =  200.0
    end if

    ! Validity check -- bounds outside which a solution cannot be found
    if (VapPres < GetSatVapPres(BOUNDS(1)) .or. VapPres > GetSatVapPres(BOUNDS(2))) then
      error stop "Error: partial pressure of water vapor is outside range of validity of equations"
    end if

    ! We use NR to approximate the solution.
    TDewPoint = TDryBulb
    lnVP = log(VapPres)
    index = 1

    do while (.true.)
      TDewPoint_iter = TDewPoint ! TDewPoint_iter used in NR calculation
      lnVP_iter = log(GetSatVapPres(TDewPoint_iter))

      ! Derivative of function, calculated analytically
      d_lnVP = dLnPws_(TDewPoint_iter)

      ! New estimate, bounded by the search domain defined above
      TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP
      TDewPoint = max(TDewPoint, BOUNDS(1))
      TDewPoint = min(TDewPoint, BOUNDS(2))

      if (abs(TDewPoint - TDewPoint_iter) <= PSYCHROLIB_TOLERANCE) then
        exit
      end if

      if (index > MAX_ITER_COUNT) then
        error stop "Convergence not reached in GetTDewPointFromVapPres. Stopping."
      end if

      index = index + 1
    end do

  TDewPoint = min(TDewPoint, TDryBulb)
  end function GetTDewPointFromVapPres