Psychrometrics

program psychrometrics
    implicit none

    ! Constants
    double precision, parameter :: R_a = 0.28705d0     ! kJ/(kg.K) gas constant for dry air
    double precision, parameter :: Cp_da = 1.006d0     ! kJ/(kg.K) specific heat of dry air
    double precision, parameter :: Cp_wv = 1.86d0      ! kJ/(kg.K) specific heat of water vapor
    double precision, parameter :: h_fg = 2501.0d0     ! kJ/kg latent heat at 0 C

    ! Inputs
    integer :: mode  ! 1 = T & RH, 2 = T & Twb, 3 = T & Tdp, 4 = T & w
    double precision :: T, P, val

    ! Solved properties
    double precision :: RH, T_wb, T_dp, w, h, v, Pv, Psat_T
    
    ! Solver bounds
    double precision :: a, b, mid, w_calc, eps
    integer :: i, max_iter
    logical :: has_error

    ! Read inputs from stdin
    read(*,*,iostat=i) mode
    if (i /= 0) then
        write(*,*) 'ERROR: Invalid input mode.'
        stop
    end if

    read(*,*,iostat=i) T
    if (i /= 0) then
        write(*,*) 'ERROR: Invalid dry-bulb temperature.'
        stop
    end if

    read(*,*,iostat=i) P
    if (i /= 0) then
        write(*,*) 'ERROR: Invalid barometric pressure.'
        stop
    end if

    read(*,*,iostat=i) val
    if (i /= 0) then
        write(*,*) 'ERROR: Invalid second parameter value.'
        stop
    end if

    ! Basic validations
    if (P <= 0.0d0) then
        write(*,*) 'ERROR: Barometric pressure must be positive.'
        stop
    end if

    has_error = .false.
    Psat_T = psat(T)

    ! Compute based on combination mode
    select case (mode)
    case (1)
        ! Mode 1: Dry-Bulb T & Relative Humidity RH (%)
        RH = val
        if (RH < 0.0d0) RH = 0.0d0
        if (RH > 100.0d0) RH = 100.0d0
        
        Pv = (RH / 100.0d0) * Psat_T
        if (Pv >= P) then
            Pv = P - 1.0d-4
            RH = (Pv / Psat_T) * 100.0d0
        end if
        
        w = 0.62198d0 * Pv / (P - Pv)
        T_dp = T_dew(Pv)
        call solve_twb(T, P, w, T_dp, T_wb)

    case (2)
        ! Mode 2: Dry-Bulb T & Wet-Bulb T_wb (C)
        T_wb = val
        if (T_wb > T) then
            write(*,*) 'ERROR: Wet-bulb temperature cannot exceed dry-bulb temperature.'
            stop
        end if
        
        w = calculate_w_from_twb(T, P, T_wb)
        if (w < 0.0d0) then
            w = 0.0d0
            Pv = 0.0d0
            RH = 0.0d0
            T_dp = -100.0d0
        else
            Pv = P * w / (0.62198d0 + w)
            if (Pv > Psat_T) then
                ! Clamping to saturation if solver yields physically impossible state
                Pv = Psat_T
                w = 0.62198d0 * Pv / (P - Pv)
                RH = 100.0d0
                T_dp = T
                T_wb = T
            else
                RH = (Pv / Psat_T) * 100.0d0
                T_dp = T_dew(Pv)
            end if
        end if

    case (3)
        ! Mode 3: Dry-Bulb T & Dew-Point T_dp (C)
        T_dp = val
        if (T_dp > T) then
            write(*,*) 'ERROR: Dew-point temperature cannot exceed dry-bulb temperature.'
            stop
        end if
        
        Pv = psat(T_dp)
        if (Pv >= P) then
            write(*,*) 'ERROR: Dew-point saturation pressure exceeds barometric pressure.'
            stop
        end if
        
        w = 0.62198d0 * Pv / (P - Pv)
        RH = (Pv / Psat_T) * 100.0d0
        call solve_twb(T, P, w, T_dp, T_wb)

    case (4)
        ! Mode 4: Dry-Bulb T & Humidity Ratio w (kg/kg)
        w = val
        if (w < 0.0d0) then
            write(*,*) 'ERROR: Humidity ratio cannot be negative.'
            stop
        end if
        
        Pv = P * w / (0.62198d0 + w)
        if (Pv > Psat_T) then
            write(*,*) 'ERROR: Specified humidity ratio exceeds saturation state at T.'
            stop
        end if
        
        RH = (Pv / Psat_T) * 100.0d0
        T_dp = T_dew(Pv)
        call solve_twb(T, P, w, T_dp, T_wb)

    case default
        write(*,*) 'ERROR: Unknown input combination mode.'
        stop
    end select

    ! Calculate Enthalpy (h in kJ/kg dry air)
    h = Cp_da * T + w * (h_fg + Cp_wv * T)

    ! Calculate Specific Volume (v in m3/kg dry air)
    v = R_a * (T + 273.15d0) / (P - Pv)

    ! Print outputs to stdout in clear key-value pairs
    write(*,'(A,I2)') 'MODE=', mode
    write(*,'(A,F14.4)') 'T=', T
    write(*,'(A,F14.4)') 'P=', P
    write(*,'(A,F14.4)') 'T_WB=', T_wb
    write(*,'(A,F14.4)') 'T_DP=', T_dp
    write(*,'(A,F14.4)') 'RH=', RH
    write(*,'(A,F14.8)') 'W=', w
    write(*,'(A,F14.4)') 'H=', h
    write(*,'(A,F14.6)') 'V=', v
    write(*,'(A,F14.6)') 'PV=', Pv
    write(*,'(A,F14.6)') 'PSAT=', Psat_T

contains

    ! Saturation pressure of liquid water or ice (Buck equation) in kPa
    double precision function psat(T_c)
        double precision, intent(in) :: T_c
        if (T_c >= 0.0d0) then
            psat = 0.61121d0 * exp((17.67d0 * T_c) / (T_c + 243.5d0))
        else
            psat = 0.61115d0 * exp((22.452d0 * T_c) / (T_c + 272.62d0))
        end if
    end function psat

    ! Dew point calculation from partial vapor pressure (Buck equation inverse) in C
    double precision function T_dew(Pv_val)
        double precision, intent(in) :: Pv_val
        double precision :: y
        if (Pv_val <= 1.0d-8) then
            T_dew = -100.0d0
            return
        end if
        if (Pv_val >= 0.61121d0) then
            y = log(Pv_val / 0.61121d0)
            T_dew = (243.5d0 * y) / (17.67d0 - y)
        else
            y = log(Pv_val / 0.61115d0)
            T_dew = (272.62d0 * y) / (22.452d0 - y)
        end if
    end function T_dew

    ! adiabatic saturation relation for humidity ratio w given dry-bulb T and wet-bulb Twb
    double precision function calculate_w_from_twb(T_db, P_bar, T_wb_val)
        double precision, intent(in) :: T_db, P_bar, T_wb_val
        double precision :: psat_twb, w_sat_twb
        
        psat_twb = psat(T_wb_val)
        if (P_bar > psat_twb) then
            w_sat_twb = 0.62198d0 * psat_twb / (P_bar - psat_twb)
        else
            w_sat_twb = 9999.0d0
        end if
        
        calculate_w_from_twb = ((2501.0d0 - 2.381d0 * T_wb_val) * w_sat_twb - Cp_da * (T_db - T_wb_val)) / &
                                (2501.0d0 + Cp_wv * T_db - 4.186d0 * T_wb_val)
    end function calculate_w_from_twb

    ! Numerical Bisection solver for wet-bulb temperature
    subroutine solve_twb(T_db, P_bar, w_act, T_dp_val, T_wb_res)
        double precision, intent(in) :: T_db, P_bar, w_act, T_dp_val
        double precision, intent(out) :: T_wb_res
        double precision :: t_low, t_high, t_m, w_m
        integer :: iter

        t_low = T_dp_val
        t_high = T_db
        
        ! Check if already saturated or bounds close
        if (abs(t_high - t_low) < 1.0d-4) then
            T_wb_res = T_db
            return
        end if

        do iter = 1, 80
            t_m = 0.5d0 * (t_low + t_high)
            w_m = calculate_w_from_twb(T_db, P_bar, t_m)
            if (w_m > w_act) then
                t_high = t_m
            else
                t_low = t_m
            end if
            if (abs(t_high - t_low) < 1.0d-5) exit
        end do
        T_wb_res = 0.5d0 * (t_low + t_high)
    end subroutine solve_twb

end program psychrometrics