program air_curtain_aerodynamics
    implicit none
    integer :: iostat_val
    double precision :: H_door_m, W_door_m, b0_mm, v0_ms, alpha_deg
    double precision :: Ti_degC, To_degC, v_wind_ms, Cp_wind
    double precision :: Ti_K, To_K, rho_i, rho_o, b0_m, alpha_rad
    double precision :: jet_momentum, thermal_head, wind_head, total_force, Dn_val
    double precision :: E_seal_pct, Q_free_m3h, Q_saved_kW, P_fan_kW
    character(len=32) :: seal_status
    double precision, parameter :: g = 9.80665d0
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_door_m       ! Door Height [m] (e.g. 3.0)
    read(*,*,iostat=iostat_val) W_door_m       ! Door Width [m] (e.g. 2.4)
    read(*,*,iostat=iostat_val) b0_mm          ! Nozzle Slot Width [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) v0_ms          ! Discharge Velocity [m/s] (e.g. 12.0)
    read(*,*,iostat=iostat_val) alpha_deg      ! Outward Angle [deg] (e.g. 15.0)
    read(*,*,iostat=iostat_val) Ti_degC        ! Indoor Temp [deg C] (e.g. 21.0)
    read(*,*,iostat=iostat_val) To_degC        ! Outdoor Temp [deg C] (e.g. 2.0)
    read(*,*,iostat=iostat_val) v_wind_ms      ! Cross Wind Speed [m/s] (e.g. 2.5)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for air curtain analysis.'
        stop
    end if

    Ti_K = Ti_degC + 273.15d0
    To_K = To_degC + 273.15d0
    rho_i = 101325.0d0 / (287.058d0 * Ti_K)
    rho_o = 101325.0d0 / (287.058d0 * To_K)
    b0_m = b0_mm * 1.0d-3
    alpha_rad = alpha_deg * (PI / 180.0d0)
    Cp_wind = 0.65d0

    ! Jet momentum per unit width [N/m]
    jet_momentum = rho_i * b0_m * (v0_ms**2) * cos(alpha_rad)

    ! Opposing forces (Thermal stack + Wind head)
    thermal_head = g * (H_door_m**2) * abs(rho_o - rho_i)
    wind_head = 0.5d0 * Cp_wind * rho_o * (v_wind_ms**2) * H_door_m
    total_force = thermal_head + wind_head

    ! Foster Deflection Number Dn
    Dn_val = jet_momentum / max(0.1d0, total_force)

    if (Dn_val < 0.12d0) then
        seal_status = 'POOR (Jet Broken by Outside Draft)'
        E_seal_pct = max(10.0d0, Dn_val * 350.0d0)
    else if (Dn_val <= 0.40d0) then
        seal_status = 'OPTIMAL (Stable Aerodynamic Barrier)'
        E_seal_pct = 70.0d0 + (Dn_val - 0.12d0) * 50.0d0
        if (E_seal_pct > 85.0d0) E_seal_pct = 85.0d0
    else
        seal_status = 'OVERPOWERED (Excessive Jet Spill)'
        E_seal_pct = 75.0d0
    end if

    ! Unprotected doorway infiltration rate [m3/h]
    Q_free_m3h = (0.35d0 * W_door_m * H_door_m * sqrt(g * H_door_m * abs(Ti_K - To_K) / Ti_K)) * 3600.0d0
    Q_saved_kW = (E_seal_pct / 100.0d0) * ((Q_free_m3h / 3600.0d0) * rho_o * 1.006d0 * abs(Ti_degC - To_degC))
    P_fan_kW = (0.5d0 * rho_i * (v0_ms**3) * b0_m * W_door_m) / (1000.0d0 * 0.60d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — AIR CURTAIN AERODYNAMICS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Doorway Height / Width    = ', H_door_m, ' m / ', W_door_m, ' m'
    write(*,'(A,F10.1,A,F6.1,A)')  'Nozzle Slot / Velocity    = ', b0_mm, ' mm / ', v0_ms, ' m/s'
    write(*,'(A,F10.1,A)')  'Discharge Angle           = ', alpha_deg, ' deg'
    write(*,'(A,F10.1,A,F6.1,A)')  'Indoor / Outdoor Temp     = ', Ti_degC, ' C / ', To_degC, ' C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.3)')    'Foster Deflection No (Dn) = ', Dn_val
    write(*,'(A,A)')        'Aerodynamic Barrier State = ', trim(seal_status)
    write(*,'(A,F10.1,A)')  'Sealing Efficiency E_seal = ', E_seal_pct, ' %'
    write(*,'(A,F10.1,A)')  'Unprotected Infiltration  = ', Q_free_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  'Thermal Energy Saved      = ', Q_saved_kW, ' kW'
    write(*,'(A,F10.2,A)')  'Air Curtain Fan Power     = ', P_fan_kW, ' kW'
    write(*,'(A)') '============================================================'

end program air_curtain_aerodynamics