Cooling Coil Psychrometrics

program cooling_coil_psychrometrics
    implicit none
    integer :: iostat_val, iter
    double precision :: Tdb1_C, RH1_pct, Q_air_m3h, BF_val, ADP_target_C, dTw_C
    double precision :: P_atm_kPa, Pvs1_kPa, Pv1_kPa, w1_kgkg, h1_kJkg, rho_air
    double precision :: mdot_air_kgs, Pvs_adp_kPa, w_adp_kgkg, h_adp_kJkg
    double precision :: Tdb2_C, w2_kgkg, h2_kJkg, RH2_pct, Pvs2_kPa, Pv2_kPa
    double precision :: Q_total_kW, Q_sens_kW, Q_lat_kW, SHR_calc, condensate_Lh, V_chw_m3h
    double precision, parameter :: P_STD = 101.325d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Tdb1_C        ! Entering Dry Bulb [deg C] (e.g. 28.0)
    read(*,*,iostat=iostat_val) RH1_pct       ! Entering Relative Humidity [%] (e.g. 55.0)
    read(*,*,iostat=iostat_val) Q_air_m3h     ! Air Volume Flow Rate [m3/h] (e.g. 5000.0)
    read(*,*,iostat=iostat_val) BF_val        ! Coil Bypass Factor BF (e.g. 0.10)
    read(*,*,iostat=iostat_val) ADP_target_C  ! Apparatus Dew Point ADP [deg C] (e.g. 9.5)
    read(*,*,iostat=iostat_val) dTw_C         ! Chilled Water Delta T [deg C] (e.g. 6.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for cooling coil psychrometrics.'
        stop
    end if

    P_atm_kPa = P_STD
    BF_val = max(0.01d0, min(0.40d0, BF_val))

    ! Saturation pressure Antoine / Magnus formula
    Pvs1_kPa = 0.61078d0 * exp((17.27d0 * Tdb1_C) / (Tdb1_C + 237.3d0))
    Pv1_kPa = (RH1_pct / 100.0d0) * Pvs1_kPa
    w1_kgkg = 0.62198d0 * Pv1_kPa / (P_atm_kPa - Pv1_kPa)
    h1_kJkg = 1.006d0 * Tdb1_C + w1_kgkg * (2501.0d0 + 1.86d0 * Tdb1_C)
    rho_air = (P_atm_kPa * 1000.0d0) / (287.058d0 * (Tdb1_C + 273.15d0) * (1.0d0 + 1.6078d0 * w1_kgkg))
    mdot_air_kgs = (Q_air_m3h / 3600.0d0) * rho_air

    ! ADP saturation state
    Pvs_adp_kPa = 0.61078d0 * exp((17.27d0 * ADP_target_C) / (ADP_target_C + 237.3d0))
    w_adp_kgkg = 0.62198d0 * Pvs_adp_kPa / (P_atm_kPa - Pvs_adp_kPa)
    h_adp_kJkg = 1.006d0 * ADP_target_C + w_adp_kgkg * (2501.0d0 + 1.86d0 * ADP_target_C)

    ! Leaving air state
    Tdb2_C = ADP_target_C + BF_val * (Tdb1_C - ADP_target_C)
    w2_kgkg = w_adp_kgkg + BF_val * (w1_kgkg - w_adp_kgkg)
    h2_kJkg = 1.006d0 * Tdb2_C + w2_kgkg * (2501.0d0 + 1.86d0 * Tdb2_C)

    Pvs2_kPa = 0.61078d0 * exp((17.27d0 * Tdb2_C) / (Tdb2_C + 237.3d0))
    Pv2_kPa = (w2_kgkg * P_atm_kPa) / (0.62198d0 + w2_kgkg)
    RH2_pct = min(100.0d0, max(0.0d0, (Pv2_kPa / Pvs2_kPa) * 100.0d0))

    ! Capacities
    Q_total_kW = mdot_air_kgs * (h1_kJkg - h2_kJkg)
    Q_sens_kW = mdot_air_kgs * 1.006d0 * (Tdb1_C - Tdb2_C)
    Q_lat_kW = max(0.0d0, Q_total_kW - Q_sens_kW)
    SHR_calc = Q_sens_kW / max(0.01d0, Q_total_kW)

    condensate_Lh = max(0.0d0, mdot_air_kgs * (w1_kgkg - w2_kgkg) * 3600.0d0)
    V_chw_m3h = (Q_total_kW * 3600.0d0) / (998.0d0 * 4.186d0 * max(1.0d0, dTw_C))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — COOLING COIL PSYCHROMETRIC ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F6.1,A)') 'Entering Air State        = ', Tdb1_C, ' C DB, ', RH1_pct, ' % RH'
    write(*,'(A,F10.2,A)')  'Entering Enthalpy h1      = ', h1_kJkg, ' kJ/kg'
    write(*,'(A,F10.4,A)')  'Entering Humidity Ratio w1= ', w1_kgkg * 1000.0d0, ' g/kg'
    write(*,'(A,F10.1,A)')  'Air Flow Rate Q           = ', Q_air_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  'Apparatus Dew Point ADP   = ', ADP_target_C, ' C'
    write(*,'(A,F10.3)')    'Coil Bypass Factor BF     = ', BF_val
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A,F6.1,A)') 'Leaving Air State         = ', Tdb2_C, ' C DB, ', RH2_pct, ' % RH'
    write(*,'(A,F10.2,A)')  'Total Cooling Capacity Qt = ', Q_total_kW, ' kW'
    write(*,'(A,F10.2,A)')  'Sensible Capacity Qs      = ', Q_sens_kW, ' kW'
    write(*,'(A,F10.2,A)')  'Latent Capacity Ql        = ', Q_lat_kW, ' kW'
    write(*,'(A,F10.3)')    'Sensible Heat Ratio SHR   = ', SHR_calc
    write(*,'(A,F10.2,A)')  'Condensate Water Rate     = ', condensate_Lh, ' L/h'
    write(*,'(A,F10.2,A)')  'Chilled Water Flow (6C dT)= ', V_chw_m3h, ' m3/h'
    write(*,'(A)') '============================================================'

end program cooling_coil_psychrometrics