program psychrometric
    implicit none
    integer :: proc_type, iostat_val, i, n_sweep
    double precision :: T_in, RH_in, P_atm, V_dot
    double precision :: T_out, RH_out
    double precision :: T_in2, RH_in2, V_dot2
    double precision :: Ps_in, W_in, h_in, v_in, Tdp_in, Twb_in
    double precision :: Ps_out, W_out, h_out, v_out, Tdp_out
    double precision :: mdot_da, mdot_da2, Q_total, Q_sens, Q_lat
    double precision :: m_cond, SHR, BPF, T_coil
    double precision :: W_mix, h_mix, T_mix, RH_mix, Ps_mix
    double precision :: frac, T_sw, Q_sw, Ps_sw, W_sw, h_sw, mdot_sw
    double precision :: Pw_in, Pw_out, Pw_in2
    character(len=50) :: proc_name

    read(*,*,iostat=iostat_val) proc_type
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid process type.'; stop; end if
    read(*,*,iostat=iostat_val) T_in
    read(*,*,iostat=iostat_val) RH_in
    read(*,*,iostat=iostat_val) P_atm
    read(*,*,iostat=iostat_val) V_dot
    read(*,*,iostat=iostat_val) T_out
    read(*,*,iostat=iostat_val) RH_out
    read(*,*,iostat=iostat_val) T_in2
    read(*,*,iostat=iostat_val) RH_in2
    read(*,*,iostat=iostat_val) V_dot2
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if

    if (P_atm <= 0.0d0) P_atm = 101.325d0
    if (V_dot <= 0.0d0) V_dot = 1.0d0
    if (RH_in < 0.0d0) RH_in = 0.0d0
    if (RH_in > 100.0d0) RH_in = 100.0d0
    if (RH_out < 0.0d0) RH_out = 0.0d0
    if (RH_out > 100.0d0) RH_out = 100.0d0
    if (proc_type < 1 .or. proc_type > 5) proc_type = 1

    select case(proc_type)
    case(1); proc_name='Sensible Heating'
    case(2); proc_name='Sensible Cooling (above dew point)'
    case(3); proc_name='Cooling & Dehumidification'
    case(4); proc_name='Adiabatic Humidification (Evap Cooling)'
    case(5); proc_name='Mixing of Two Airstreams'
    end select

    ! ── Psychrometric functions (inline) ──────────────────────
    ! P_sat (Magnus formula, kPa, T in C)
    Ps_in = 0.6105d0 * exp(17.27d0*T_in/(T_in+237.3d0))
    Pw_in = RH_in/100.0d0 * Ps_in
    W_in = 0.622d0 * Pw_in / max(P_atm - Pw_in, 0.01d0)
    h_in = 1.006d0*T_in + W_in*(2501.0d0 + 1.86d0*T_in)
    v_in = 0.2871d0*(T_in+273.15d0)/(P_atm - Pw_in) * (1.0d0+1.6078d0*W_in)
    ! dew point (inverse Magnus)
    if (Pw_in > 0.001d0) then
        Tdp_in = 237.3d0*log(Pw_in/0.6105d0)/(17.27d0 - log(Pw_in/0.6105d0))
    else
        Tdp_in = -40.0d0
    end if
    ! wet bulb (Stull approximation 2011)
    Twb_in = T_in * atan(0.151977d0*sqrt(RH_in+8.313659d0)) + &
             atan(T_in+RH_in) - atan(RH_in-1.676331d0) + &
             0.00391838d0*RH_in**1.5d0*atan(0.023101d0*RH_in) - 4.686035d0

    mdot_da = V_dot / max(v_in, 0.01d0)

    ! ── Process calculations ──────────────────────────────────
    Q_total=0; Q_sens=0; Q_lat=0; m_cond=0; SHR=0; BPF=0; T_coil=0

    if (proc_type == 1 .or. proc_type == 2) then
        ! Sensible heating/cooling: W constant
        W_out = W_in
        Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
        if (Ps_out > 0.0d0) then
            RH_out = W_out * P_atm / (0.622d0 + W_out) / Ps_out * 100.0d0
        else
            RH_out = 0.0d0
        end if
        if (RH_out > 100.0d0) RH_out = 100.0d0
        h_out = 1.006d0*T_out + W_out*(2501.0d0+1.86d0*T_out)
        Q_total = mdot_da * (h_out - h_in)
        Q_sens = Q_total
        Q_lat = 0.0d0
        SHR = 1.0d0

    else if (proc_type == 3) then
        ! Cooling & dehumidification
        Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
        Pw_out = RH_out/100.0d0 * Ps_out
        W_out = 0.622d0*Pw_out/max(P_atm-Pw_out,0.01d0)
        if (W_out > W_in) W_out = W_in  ! can't add moisture in cooling
        h_out = 1.006d0*T_out + W_out*(2501.0d0+1.86d0*T_out)
        Q_total = mdot_da * (h_in - h_out)  ! positive = heat removed
        Q_sens = mdot_da * 1.006d0 * (T_in - T_out)
        Q_lat = Q_total - Q_sens
        if (Q_total > 0.0d0) SHR = Q_sens / Q_total
        m_cond = mdot_da * max(W_in - W_out, 0.0d0)
        ! Coil temperature estimate (ADP)
        if (RH_out > 90.0d0) then
            T_coil = T_out - 2.0d0
        else
            T_coil = T_out - 5.0d0
        end if
        if (abs(T_in - T_coil) > 0.01d0) then
            BPF = (T_out - T_coil) / (T_in - T_coil)
        end if

    else if (proc_type == 4) then
        ! Adiabatic humidification (evaporative cooling)
        ! h ≈ constant, follow wet-bulb line
        ! T_out ≈ T_wb + (T_in - T_wb)*(1-eta), eta~0.85
        T_out = Twb_in + (T_in - Twb_in)*0.15d0  ! 85% effectiveness
        h_out = h_in  ! adiabatic
        W_out = (h_out - 1.006d0*T_out) / (2501.0d0 + 1.86d0*T_out)
        if (W_out < 0.0d0) W_out = 0.0d0
        Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
        Pw_out = W_out * P_atm / (0.622d0 + W_out)
        if (Ps_out > 0.0d0) then
            RH_out = Pw_out / Ps_out * 100.0d0
        else
            RH_out = 100.0d0
        end if
        if (RH_out > 100.0d0) RH_out = 100.0d0
        Q_total = 0.0d0
        Q_sens = mdot_da * 1.006d0 * (T_out - T_in)
        Q_lat = -Q_sens  ! equal and opposite

    else if (proc_type == 5) then
        ! Mixing
        Ps_out = 0.6105d0*exp(17.27d0*T_in2/(T_in2+237.3d0))
        Pw_in2 = RH_in2/100.0d0 * Ps_out
        if (V_dot2 <= 0.0d0) V_dot2 = V_dot * 0.5d0
        ! Stream 2 properties
        W_out = 0.622d0*Pw_in2/max(P_atm - Pw_in2, 0.01d0)  ! W of stream 2
        v_out = 0.2871d0*(T_in2+273.15d0)/(P_atm-Pw_in2)*(1.0d0+1.6078d0*W_out)
        mdot_da2 = V_dot2 / max(v_out, 0.01d0)
        h_out = 1.006d0*T_in2 + W_out*(2501.0d0+1.86d0*T_in2)

        ! Mixed state
        W_mix = (mdot_da*W_in + mdot_da2*W_out) / max(mdot_da+mdot_da2, 1.0d-10)
        h_mix = (mdot_da*h_in + mdot_da2*h_out) / max(mdot_da+mdot_da2, 1.0d-10)
        ! T_mix from h_mix = 1.006*T + W_mix*(2501+1.86*T)
        T_mix = (h_mix - W_mix*2501.0d0) / (1.006d0 + W_mix*1.86d0)
        Ps_mix = 0.6105d0*exp(17.27d0*T_mix/(T_mix+237.3d0))
        Pw_out = W_mix*P_atm/(0.622d0+W_mix)
        if (Ps_mix > 0.0d0) then
            RH_mix = Pw_out/Ps_mix*100.0d0
        else
            RH_mix = 0.0d0
        end if
        if (RH_mix > 100.0d0) RH_mix = 100.0d0

        ! Set "out" = mixed for display
        T_out = T_mix; RH_out = RH_mix; W_out = W_mix; h_out = h_mix
        Q_total = 0.0d0  ! adiabatic mixing
    end if

    ! Outlet dew point
    Pw_out = W_out * P_atm / max(0.622d0 + W_out, 1.0d-10)
    if (Pw_out > 0.001d0) then
        Tdp_out = 237.3d0*log(Pw_out/0.6105d0)/(17.27d0-log(Pw_out/0.6105d0))
    else
        Tdp_out = -40.0d0
    end if
    v_out = 0.2871d0*(T_out+273.15d0)/(P_atm-Pw_out)*(1.0d0+1.6078d0*W_out)

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   PSYCHROMETRIC PROCESS ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Process Type              = ', trim(proc_name)
    write(*,'(A,F10.2,A)')  '  Atmospheric Pressure      = ', P_atm, ' kPa'
    write(*,'(A,F10.4,A)')  '  Airflow Rate              = ', V_dot, ' m3/s'
    write(*,'(A,F10.4,A)')  '  Dry Air Mass Flow         = ', mdot_da, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- INLET AIR STATE -----------------------------------------'
    write(*,'(A,F10.2,A)')  '  T_db (dry bulb)           = ', T_in, ' C'
    write(*,'(A,F10.2,A)')  '  RH                        = ', RH_in, ' percent'
    write(*,'(A,F10.6,A)')  '  W (humidity ratio)        = ', W_in, ' kg/kg_da'
    write(*,'(A,F10.4,A)')  '  W                         = ', W_in*1000.0d0, ' g/kg_da'
    write(*,'(A,F10.4,A)')  '  h (enthalpy)              = ', h_in, ' kJ/kg_da'
    write(*,'(A,F10.4,A)')  '  v (specific volume)       = ', v_in, ' m3/kg_da'
    write(*,'(A,F10.2,A)')  '  T_dp (dew point)          = ', Tdp_in, ' C'
    write(*,'(A,F10.2,A)')  '  T_wb (wet bulb)           = ', Twb_in, ' C'
    write(*,'(A,F10.4,A)')  '  P_sat                     = ', Ps_in, ' kPa'
    write(*,*)
    write(*,'(A)') '--- OUTLET AIR STATE ----------------------------------------'
    write(*,'(A,F10.2,A)')  '  T_db (dry bulb)           = ', T_out, ' C'
    write(*,'(A,F10.2,A)')  '  RH                        = ', RH_out, ' percent'
    write(*,'(A,F10.6,A)')  '  W (humidity ratio)        = ', W_out, ' kg/kg_da'
    write(*,'(A,F10.4,A)')  '  W                         = ', W_out*1000.0d0, ' g/kg_da'
    write(*,'(A,F10.4,A)')  '  h (enthalpy)              = ', h_out, ' kJ/kg_da'
    write(*,'(A,F10.4,A)')  '  v (specific volume)       = ', v_out, ' m3/kg_da'
    write(*,'(A,F10.2,A)')  '  T_dp (dew point)          = ', Tdp_out, ' C'
    write(*,*)
    write(*,'(A)') '--- PROCESS ENERGY ANALYSIS ---------------------------------'
    write(*,'(A,F12.4,A)')  '  Q_total                   = ', Q_total, ' kW'
    write(*,'(A,F12.4,A)')  '  Q_sensible                = ', Q_sens, ' kW'
    write(*,'(A,F12.4,A)')  '  Q_latent                  = ', Q_lat, ' kW'
    write(*,'(A,F10.4)')    '  SHR (sensible heat ratio) = ', SHR
    if (proc_type == 3) then
    write(*,'(A,F12.6,A)')  '  Condensate Rate           = ', m_cond, ' kg/s'
    write(*,'(A,F12.4,A)')  '  Condensate Rate           = ', m_cond*3600.0d0, ' kg/h'
    write(*,'(A,F10.4)')    '  Bypass Factor BPF         = ', BPF
    write(*,'(A,F10.2,A)')  '  Coil Temp (est. ADP)      = ', T_coil, ' C'
    end if
    if (proc_type == 5) then
    write(*,'(A,F10.2,A)')  '  Stream 2 T_db             = ', T_in2, ' C'
    write(*,'(A,F10.2,A)')  '  Stream 2 RH               = ', RH_in2, ' percent'
    write(*,'(A,F10.4,A)')  '  Stream 2 flow             = ', V_dot2, ' m3/s'
    write(*,'(A,F10.4,A)')  '  Stream 2 mdot_da          = ', mdot_da2, ' kg/s'
    end if
    write(*,*)

    ! ── Sweep ─────────────────────────────────────────────────
    n_sweep = 40
    if (proc_type /= 5) then
        write(*,'(A)') '--- SENSITIVITY: Q VS OUTLET TEMPERATURE --------------------'
        write(*,'(A)') '  T_out[C]      Q_total[kW]   W_out[g/kg]'
        write(*,'(A)') '  -----------------------------------------------------------'
        do i = 1, n_sweep
            T_sw = (T_in-30.0d0) + dble(i-1)*60.0d0/dble(n_sweep-1)
            Ps_sw = 0.6105d0*exp(17.27d0*T_sw/(T_sw+237.3d0))
            if (proc_type <= 2) then
                W_sw = W_in
            else if (proc_type == 3) then
                Pw_out = RH_out/100.0d0 * Ps_sw
                W_sw = 0.622d0*Pw_out/max(P_atm-Pw_out,0.01d0)
                if (W_sw > W_in) W_sw = W_in
            else
                h_sw = h_in
                W_sw = (h_sw - 1.006d0*T_sw)/(2501.0d0+1.86d0*T_sw)
                if (W_sw < 0.0d0) W_sw = 0.0d0
            end if
            h_sw = 1.006d0*T_sw + W_sw*(2501.0d0+1.86d0*T_sw)
            if (proc_type == 3) then
                Q_sw = mdot_da*(h_in - h_sw)
            else
                Q_sw = mdot_da*(h_sw - h_in)
            end if
            write(*,'(F10.2,4X,F12.4,4X,F10.4)') T_sw, Q_sw, W_sw*1000.0d0
        end do
    else
        write(*,'(A)') '--- SENSITIVITY: MIXING FRACTION ----------------------------'
        write(*,'(A)') '  Fraction_2    T_mix[C]      RH_mix[%]     W_mix[g/kg]'
        write(*,'(A)') '  -----------------------------------------------------------'
        do i = 1, n_sweep
            frac = dble(i-1)/dble(n_sweep-1)
            W_sw = (1.0d0-frac)*W_in + frac*W_out
            h_sw = (1.0d0-frac)*h_in + frac*h_out
            T_sw = (h_sw - W_sw*2501.0d0)/(1.006d0+W_sw*1.86d0)
            Ps_sw = 0.6105d0*exp(17.27d0*T_sw/(T_sw+237.3d0))
            Pw_out = W_sw*P_atm/(0.622d0+W_sw)
            if (Ps_sw > 0.0d0) then
                mdot_sw = Pw_out/Ps_sw*100.0d0
            else
                mdot_sw = 0.0d0
            end if
            if (mdot_sw > 100.0d0) mdot_sw = 100.0d0
            write(*,'(F10.4,4X,F10.2,4X,F10.2,4X,F10.4)') frac, T_sw, mdot_sw, W_sw*1000.0d0
        end do
    end if
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  P_sat: Magnus formula 0.6105*exp(17.27*T/(T+237.3)) [kPa]'
    write(*,'(A)') '  W = 0.622*Pw/(P_atm - Pw) [kg_w/kg_da]'
    write(*,'(A)') '  h = 1.006*T + W*(2501 + 1.86*T) [kJ/kg_da]'
    write(*,'(A)') '  T_dp: inverse Magnus equation'
    write(*,'(A)') '  T_wb: Stull (2011) approximation'
    write(*,'(A)') '  SHR = Q_sensible / Q_total'

end program psychrometric