program steam_humidifier_dispersion
    implicit none
    integer :: n_tubes, iostat_val
    double precision :: Q_air_m3h, H_duct_mm, W_duct_mm, Tdb1_C, RH1_pct, RH2_pct
    double precision :: Pvs1_kPa, Pv1_kPa, w1_kgkg, Pvs2_kPa, Pv2_kPa, w2_kgkg
    double precision :: rho_air, mdot_air_kgs, mdot_steam_kgh, v_air_ms
    double precision :: D_abs_m, P_boiler_kW

    ! Read inputs
    read(*,*,iostat=iostat_val) Q_air_m3h       ! Air Volume Flow Rate [m3/h] (e.g. 6000.0)
    read(*,*,iostat=iostat_val) H_duct_mm       ! Duct Height [mm] (e.g. 600.0)
    read(*,*,iostat=iostat_val) W_duct_mm       ! Duct Width [mm] (e.g. 800.0)
    read(*,*,iostat=iostat_val) Tdb1_C          ! Air Temp [deg C] (e.g. 21.0)
    read(*,*,iostat=iostat_val) RH1_pct         ! Entering Relative Humidity [%] (e.g. 20.0)
    read(*,*,iostat=iostat_val) RH2_pct         ! Target Relative Humidity [%] (e.g. 50.0)
    read(*,*,iostat=iostat_val) n_tubes         ! Number of Dispersion Tubes (1=Single, >1=Rapid Grid)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for steam humidifier dispersion.'
        stop
    end if

    if (RH2_pct <= RH1_pct .or. RH2_pct >= 98.0d0) then
        write(*,*) 'ERROR: Target RH must be greater than entering RH and < 98%.'
        stop
    end if

    ! Psychrometrics
    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 / (101.325d0 - Pv1_kPa)

    Pv2_kPa = (RH2_pct / 100.0d0) * Pvs1_kPa
    w2_kgkg = 0.62198d0 * Pv2_kPa / (101.325d0 - Pv2_kPa)

    rho_air = 101325.0d0 / (287.058d0 * (Tdb1_C + 273.15d0))
    mdot_air_kgs = (Q_air_m3h / 3600.0d0) * rho_air
    mdot_steam_kgh = mdot_air_kgs * (w2_kgkg - w1_kgkg) * 3600.0d0

    v_air_ms = (Q_air_m3h / 3600.0d0) / ((H_duct_mm * 1.0d-3) * (W_duct_mm * 1.0d-3))

    ! Non-wetting absorption distance [m]
    D_abs_m = (1.85d0 / dble(max(1, n_tubes)**0.75d0)) * max(1.0d0, v_air_ms / 2.5d0) * &
              (((RH2_pct - RH1_pct) / max(1.0d0, 100.0d0 - RH2_pct))**0.65d0)

    ! Electric steam generator power [kW]
    P_boiler_kW = (mdot_steam_kgh * 2550.0d0) / 3600.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — STEAM HUMIDIFICATION & DISPERSION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F6.2,A)') 'Air Flow / Velocity       = ', Q_air_m3h, ' m3/h / ', v_air_ms, ' m/s'
    write(*,'(A,F10.1,A,F6.1,A)') 'Entering / Target RH      = ', RH1_pct, ' % / ', RH2_pct, ' %'
    write(*,'(A,F10.2,A,F6.2,A)') 'Humidity Ratio (w1 / w2)  = ', w1_kgkg*1000.0d0, ' g/kg / ', w2_kgkg*1000.0d0, ' g/kg'
    write(*,'(A,I2)')       'Dispersion Manifold Tubes = ', n_tubes
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Steam Generation Rate     = ', mdot_steam_kgh, ' kg/h'
    write(*,'(A,F10.2,A)')  'Non-Wetting Absorpt Dist  = ', D_abs_m, ' meters'
    write(*,'(A,F10.2,A)')  'Humidifier Electric Power = ', P_boiler_kW, ' kW'
    write(*,'(A)') '============================================================'

end program steam_humidifier_dispersion
