program stack_effect_ventilation
    implicit none
    integer :: iostat_val
    double precision :: H_stack_m, A1_m2, A2_m2, Ti_degC, To_degC
    double precision :: v_wind_ms, dCp_val, V_bldg_m3, Cd_val
    double precision :: Ti_K, To_K, rho_o, dP_stack_Pa, dP_wind_Pa, dP_total_Pa
    double precision :: A_eff_m2, Q_vent_m3s, Q_vent_m3h, ACH_val
    double precision, parameter :: g = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_stack_m       ! Stack / Atrium Height [m] (e.g. 15.0)
    read(*,*,iostat=iostat_val) A1_m2           ! Inlet Opening Area [m2] (e.g. 3.0)
    read(*,*,iostat=iostat_val) A2_m2           ! Outlet Opening Area [m2] (e.g. 3.5)
    read(*,*,iostat=iostat_val) Ti_degC         ! Indoor Air Temp [deg C] (e.g. 24.0)
    read(*,*,iostat=iostat_val) To_degC         ! Outdoor Air Temp [deg C] (e.g. 12.0)
    read(*,*,iostat=iostat_val) v_wind_ms       ! Wind Velocity [m/s] (e.g. 3.5)
    read(*,*,iostat=iostat_val) dCp_val         ! Wind Pressure Coeff Delta Cp (e.g. 0.70)
    read(*,*,iostat=iostat_val) V_bldg_m3       ! Building Volume [m3] (e.g. 4500.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for stack effect ventilation.'
        stop
    end if

    Ti_K = Ti_degC + 273.15d0
    To_K = To_degC + 273.15d0
    rho_o = 101325.0d0 / (287.058d0 * To_K)
    Cd_val = 0.62d0

    ! Stack Pressure Difference [Pa]
    dP_stack_Pa = rho_o * g * H_stack_m * ((Ti_K - To_K) / Ti_K)

    ! Wind Pressure Difference [Pa]
    dP_wind_Pa = dCp_val * 0.5d0 * rho_o * (v_wind_ms**2)

    ! Combined driving head (quadrature sum per British Standard BS 5925 / CIBSE)
    dP_total_Pa = sqrt((dP_stack_Pa**2) + (dP_wind_Pa**2))

    ! Effective Opening Area [m2]
    A_eff_m2 = (A1_m2 * A2_m2) / sqrt((A1_m2**2) + (A2_m2**2))

    ! Ventilation Volumetric Flow Rate
    Q_vent_m3s = Cd_val * A_eff_m2 * sqrt((2.0d0 * max(0.01d0, dP_total_Pa)) / rho_o)
    Q_vent_m3h = Q_vent_m3s * 3600.0d0
    ACH_val = Q_vent_m3h / max(1.0d0, V_bldg_m3)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — STACK EFFECT & NATURAL VENTILATION'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Stack / Atrium Height     = ', H_stack_m, ' m'
    write(*,'(A,F10.2,A,F10.2,A)') 'Opening Areas (In/Out)    = ', A1_m2, ' m2 / ', A2_m2, ' m2'
    write(*,'(A,F10.2,A)')  'Effective Opening Area    = ', A_eff_m2, ' m2'
    write(*,'(A,F10.1,A,F10.1,A)') 'Indoor / Outdoor Temp     = ', Ti_degC, ' C / ', To_degC, ' C'
    write(*,'(A,F10.2,A)')  'Wind Speed                = ', v_wind_ms, ' m/s'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Thermal Stack Pressure    = ', dP_stack_Pa, ' Pa'
    write(*,'(A,F10.2,A)')  'Wind Surface Pressure     = ', dP_wind_Pa, ' Pa'
    write(*,'(A,F10.2,A)')  'Total Driving Pressure    = ', dP_total_Pa, ' Pa'
    write(*,'(A,F10.2,A)')  'Ventilation Airflow Rate  = ', Q_vent_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  'Air Changes per Hour (ACH)= ', ACH_val, ' h-1'
    write(*,'(A)') '============================================================'

end program stack_effect_ventilation
