program natural_convection_cavity_enclosure
    implicit none
    integer :: iostat_val
    double precision :: H_height_m, W_width_m, L_depth_m, Th_hot_C, Tc_cold_C
    double precision :: rho_dens, mu_cP, k_cond, cp_spec, beta_exp
    double precision :: mu_Pas, nu_visc, alpha_diff, DeltaT, Ra_W, Pr_num, AspectRatio
    double precision :: Nu_avg, h_avg, Q_dot_W, delta_bl_mm, u_buoy_ms
    double precision, parameter :: G_GRAV = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_height_m    ! Cavity Height [m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) W_width_m     ! Cavity Spacing Width [m] (e.g. 0.05)
    read(*,*,iostat=iostat_val) L_depth_m     ! Cavity Depth [m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) Th_hot_C      ! Hot Wall Temp [deg C] (e.g. 50.0)
    read(*,*,iostat=iostat_val) Tc_cold_C     ! Cold Wall Temp [deg C] (e.g. 10.0)
    read(*,*,iostat=iostat_val) rho_dens      ! Fluid Density [kg/m3] (e.g. 1.164 for air)
    read(*,*,iostat=iostat_val) mu_cP         ! Fluid Viscosity [cP] (e.g. 0.0185 for air)
    read(*,*,iostat=iostat_val) k_cond        ! Thermal Conductivity [W/(m.K)] (e.g. 0.026)
    read(*,*,iostat=iostat_val) cp_spec       ! Specific Heat [J/(kg.K)] (e.g. 1005.0)
    read(*,*,iostat=iostat_val) beta_exp      ! Volumetric Expansion Coeff [1/K] (e.g. 0.00333)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for cavity natural convection calculation.'
        stop
    end if

    if (H_height_m <= 0.0d0 .or. W_width_m <= 0.0d0 .or. L_depth_m <= 0.0d0) then
        write(*,*) 'ERROR: Geometry dimensions must be positive.'
        stop
    end if

    DeltaT = abs(Th_hot_C - Tc_cold_C)
    if (DeltaT <= 0.0d0) DeltaT = 1.0d-3

    mu_Pas = mu_cP * 1.0d-3
    nu_visc = mu_Pas / rho_dens
    alpha_diff = k_cond / (rho_dens * cp_spec)
    Pr_num = nu_visc / alpha_diff
    AspectRatio = H_height_m / W_width_m

    Ra_W = (G_GRAV * beta_exp * DeltaT * (W_width_m**3)) / (nu_visc * alpha_diff)

    ! Catton / Berkovsky-Polevikov Correlation for Vertical Enclosure
    if (Ra_W < 1000.0d0) then
        Nu_avg = 1.0d0 ! Pure conduction
    else if (Ra_W < 1.0d4) then
        Nu_avg = 0.18d0 * ((Ra_W * Pr_num / (0.2d0 + Pr_num))**0.29d0)
    else
        Nu_avg = 0.22d0 * ((Ra_W * Pr_num / (0.2d0 + Pr_num))**0.28d0) * (AspectRatio**(-0.25d0))
    end if
    if (Nu_avg < 1.0d0) Nu_avg = 1.0d0

    h_avg = (Nu_avg * k_cond) / W_width_m
    Q_dot_W = h_avg * (H_height_m * L_depth_m) * DeltaT
    delta_bl_mm = (W_width_m / (max(1.0d0, Ra_W)**0.25d0)) * 1000.0d0
    u_buoy_ms = sqrt(G_GRAV * beta_exp * DeltaT * H_height_m)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CAVITY NATURAL CONVECTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cavity H x W (Aspect A)   = ', H_height_m, ' m x ', W_width_m, ' m (A = ', AspectRatio, ')'
    write(*,'(A,F10.1,A,F10.1,A)') 'Hot Wall Th / Cold Wall Tc= ', Th_hot_C, ' deg C / ', Tc_cold_C, ' deg C'
    write(*,'(A,E12.4,A,F10.2)')   'Cavity Rayleigh Ra_W / Pr = ', Ra_W, ' / ', Pr_num
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2)')           'AVERAGE NUSSELT (Nu_avg)  = ', Nu_avg
    write(*,'(A,F10.2,A)')  'Heat Transfer Coeff h_avg = ', h_avg, ' W/(m2.K)'
    write(*,'(A,F10.2,A,F10.2,A)') 'TOTAL HEAT FLOW (Q_dot)   = ', Q_dot_W, ' W (', Q_dot_W/1000.0d0, ' kW)'
    write(*,'(A,F10.2,A,F10.3,A)') 'BL Thickness / Max Buoy U = ', delta_bl_mm, ' mm / ', u_buoy_ms, ' m/s'
    write(*,'(A)') '============================================================'

end program natural_convection_cavity_enclosure