GetTWetBulbFromHumRatio Function

  function GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) result(TWetBulb)
    !+ Return wet-bulb temperature given dry-bulb temperature, humidity ratio, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35 solved for Tstar

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real              ::  TDewPoint
      !+ TDewPoint : Dew-point temperature in °F [IP] or °C [SI]
    real              ::  TWetBulbSup
      !+ Upper value of wet bulb temperature in bissection method (initial guess is from dry bulb temperature) in °F [IP] or °C [SI]
    real              ::  TWetBulbInf
      !+ Lower value of wet bulb temperature in bissection method (initial guess is from dew point temperature) in °F [IP] or °C [SI]
    real              ::  Wstar
      !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  BoundedHumRatio
      !+ Humidity ratio bounded to MIN_HUM_RATIO
    integer           ::  index
      !+ index used in iteration

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio cannot be negative"
    end if
    BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO)

    TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure)

    ! Initial guesses
    TWetBulbSup = TDryBulb
    TWetBulbInf = TDewPoint
    TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.0

    index = 1
    ! Bisection loop
    do while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE)

    ! Compute humidity ratio at temperature Tstar
    Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)

    ! Get new bounds
    if (Wstar > BoundedHumRatio) then
      TWetBulbSup = TWetBulb
    else
      TWetBulbInf = TWetBulb
    end if

    ! New guess of wet bulb temperature
    TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.0

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

    index = index + 1
    end do
  end function GetTWetBulbFromHumRatio