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
     real ::  T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH

    ! Bounds and step size as a function of the system of units
    if (isIP()) then
        BOUNDS(1) = -148.0
        BOUNDS(2) =  392.0
        T_WATER_FREEZE = 32.
    else
        BOUNDS(1) = -100.0
        BOUNDS(2) =  200.0
        T_WATER_FREEZE = 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

    ! Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing
    T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.        ! Temperature just below freezing
    T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.        ! Temperature just below freezing
    T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE             ! Temperature just above freezing
    PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW)
    PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH)

    ! Restrict iteration to either left or right part of the saturation vapor pressure curve
    ! to avoid iterating back and forth across the discontinuity of the curve at the freezing point
    ! When the partial pressure of water vapor is within the discontinuity of GetSatVapPres,
    ! simply return the freezing point of water.
    if (VapPres < PWS_FREEZE_LOW) then
        BOUNDS(2) = T_WATER_FREEZE_LOW
    else if (VapPres > PWS_FREEZE_HIGH) then
        BOUNDS(1) = T_WATER_FREEZE_HIGH
    else
        TDewPoint = T_WATER_FREEZE
        return
    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