program natural_convection
    implicit none

    ! Inputs
    integer :: geom_type
    double precision :: dim1, dim2
    double precision :: T_s, T_inf
    double precision :: rho, mu, k_f, Pr, Cp, beta

    ! Constants
    double precision, parameter :: g = 9.80665d0
    double precision, parameter :: PI = 3.141592653589793d0

    ! Physical properties
    double precision :: nu, L_c, A_s, T_f, dT
    double precision :: Gr, Ra, Nu_avg, h_avg, Q

    ! Local profile variables
    integer :: i, n_points
    double precision :: x, dx, local_Ra, local_Gr, local_Nu, local_h, local_delta
    character(len=50) :: geom_name
    integer :: iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry type input.'
        stop
    end if

    read(*,*,iostat=iostat_val) dim1
    read(*,*,iostat=iostat_val) dim2
    read(*,*,iostat=iostat_val) T_s
    read(*,*,iostat=iostat_val) T_inf
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) k_f
    read(*,*,iostat=iostat_val) Pr
    read(*,*,iostat=iostat_val) Cp
    read(*,*,iostat=iostat_val) beta

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all convection parameters.'
        stop
    end if

    ! Basic validations
    if (dim1 <= 0.0d0) then
        write(*,*) 'ERROR: Dimension 1 must be positive.'
        stop
    end if
    if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. k_f <= 0.0d0 .or. Pr <= 0.0d0 .or. Cp <= 0.0d0) then
        write(*,*) 'ERROR: Fluid properties must be positive.'
        stop
    end if

    ! Compute film temperature and kinematic viscosity
    T_f = (T_s + T_inf) / 2.0d0
    nu = mu / rho
    dT = abs(T_s - T_inf)

    ! Setup Geometry parameters
    select case (geom_type)
    case (1)
        geom_name = "Vertical Plate"
        L_c = dim1
        A_s = dim1 * dim2
    case (2)
        geom_name = "Horizontal Plate (Upper Hot / Lower Cold)"
        if (dim2 <= 0.0d0) dim2 = dim1 ! Default to square if W is zero
        L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2))
        A_s = dim1 * dim2
    case (3)
        geom_name = "Horizontal Plate (Lower Hot / Upper Cold)"
        if (dim2 <= 0.0d0) dim2 = dim1
        L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2))
        A_s = dim1 * dim2
    case (4)
        geom_name = "Horizontal Cylinder"
        if (dim2 <= 0.0d0) dim2 = 1.0d0 ! Default to 1m length
        L_c = dim1
        A_s = PI * dim1 * dim2
    case (5)
        geom_name = "Sphere"
        L_c = dim1
        A_s = PI * dim1**2
    case default
        write(*,*) 'ERROR: Invalid geometry type selected.'
        stop
    end select

    ! Grashof and Rayleigh numbers
    Gr = (g * beta * dT * L_c**3) / (nu**2)
    Ra = Gr * Pr

    ! Nusselt Number Calculations
    select case (geom_type)
    case (1)
        ! Churchill and Chu (1975) for Vertical Plate (All Ra)
        Nu_avg = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
                 ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
    case (2)
        ! Horizontal Plate (Upper Hot / Lower Cold)
        if (Ra < 1.0d4) then
            Nu_avg = 0.54d0 * Ra**0.25d0 ! Extrapolated/lower-laminar
        else if (Ra <= 1.0d7) then
            Nu_avg = 0.54d0 * Ra**0.25d0
        else if (Ra <= 1.0d11) then
            Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0)
        else
            Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0) ! Extrapolated/upper-turbulent
        end if
    case (3)
        ! Horizontal Plate (Lower Hot / Upper Cold)
        ! Valid for 10^5 to 10^11
        Nu_avg = 0.27d0 * Ra**0.25d0
    case (4)
        ! Churchill and Chu (1975) for Horizontal Cylinder (Ra <= 10^12)
        Nu_avg = (0.60d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
                 ((1.0d0 + (0.559d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
    case (5)
        ! Churchill (1983) for Sphere (Ra <= 10^11 and Pr >= 0.7)
        Nu_avg = 2.0d0 + (0.589d0 * Ra**0.25d0) / &
                 ((1.0d0 + (0.469d0/Pr)**(9.0d0/16.0d0))**(4.0d0/9.0d0))
    end select

    ! Average heat transfer coefficient
    h_avg = Nu_avg * k_f / L_c

    ! Heat transfer rate
    Q = h_avg * A_s * dT

    ! Output results
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   NATURAL CONVECTION HEAT TRANSFER ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Geometry Name           = ', trim(geom_name)
    write(*,'(A,I2)')         '  Geometry Type Code      = ', geom_type
    write(*,'(A,F12.4,A)')  '  Characteristic Length   = ', L_c, ' m'
    write(*,'(A,F12.4,A)')  '  Surface Area (As)       = ', A_s, ' m2'
    write(*,'(A,F12.2,A)')  '  Surface Temperature     = ', T_s, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Ambient Temperature     = ', T_inf, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Film Temperature        = ', T_f, ' deg-C'
    write(*,*)
    write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Density (rho)             = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity (mu)   = ', mu, ' Pa.s'
    write(*,'(A,ES12.4,A)') '  Kinematic Viscosity (nu) = ', nu, ' m2/s'
    write(*,'(A,F12.4,A)')  '  Thermal Conductivity    = ', k_f, ' W/m.K'
    write(*,'(A,F12.4)')    '  Prandtl Number (Pr)     = ', Pr
    write(*,'(A,F12.2,A)')  '  Specific Heat (Cp)      = ', Cp, ' J/kg.K'
    write(*,'(A,ES12.4,A)') '  Expansion Coeff (beta)   = ', beta, ' 1/K'
    write(*,*)
    write(*,'(A)') '--- CONVECTION RESULTS --------------------------------------'
    write(*,'(A,ES12.4)')   '  Grashof Number (Gr)     = ', Gr
    write(*,'(A,ES12.4)')   '  Rayleigh Number (Ra)    = ', Ra
    write(*,'(A,F12.4)')    '  Average Nusselt Number  = ', Nu_avg
    write(*,'(A,F12.4,A)')  '  Average h               = ', h_avg, ' W/m2.K'
    write(*,'(A,F12.4,A)')  '  Heat Transfer Rate (Q)  = ', Q, ' W'
    write(*,*)

    ! ============================================
    ! LOCAL PROFILE DATA (along coordinate x)
    ! ============================================
    write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
    write(*,'(A)') '  x [m]       Ra_x          Nu_x         h_x [W/m2.K]  d [mm]'
    write(*,'(A)') '  ----------------------------------------------------------'

    n_points = 40
    dx = L_c / dble(n_points)

    do i = 1, n_points
        x = dble(i) * dx
        local_Gr = (g * beta * dT * x**3) / (nu**2)
        local_Ra = local_Gr * Pr

        if (dT <= 1.0d-10 .or. local_Gr <= 1.0d-10) then
            local_Nu = 0.0d0
            local_delta = 0.0d0
        else
            if (local_Ra < 1.0d9) then
                ! Laminar boundary layer correlations
                local_Nu = 0.508d0 * (Pr / (0.952d0 + Pr))**0.25d0 * local_Ra**0.25d0
                local_delta = 3.93d0 * x * ((0.952d0 + Pr) / (Pr**2 * local_Gr))**0.25d0
            else
                ! Turbulent boundary layer approximations
                local_Nu = 0.0292d0 * local_Ra**0.4d0
                local_delta = 0.15d0 * x * local_Ra**(-0.1d0)
            end if
        end if

        if (x > 0.0d0) then
            local_h = local_Nu * k_f / x
        else
            local_h = 0.0d0
        end if

        write(*,'(F8.4,2X,ES12.4,2X,F10.2,4X,F10.4,4X,F10.4)') &
            x, local_Ra, local_Nu, local_h, local_delta*1000.0d0
    end do

    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    select case (geom_type)
    case (1)
        write(*,'(A)') '  Vertical Plate: Churchill and Chu (1975) Correlation'
    case (2)
        write(*,'(A)') '  Horizontal Plate (Upper Hot): Nu = 0.54 Ra^0.25 (lam), Nu = 0.15 Ra^(1/3) (turb)'
    case (3)
        write(*,'(A)') '  Horizontal Plate (Lower Hot): Nu = 0.27 Ra^0.25'
    case (4)
        write(*,'(A)') '  Horizontal Cylinder: Churchill and Chu (1975) Correlation'
    case (5)
        write(*,'(A)') '  Sphere: Churchill (1983) Correlation'
    end select
    write(*,'(A)') '  Boundary Layer: Laminar similarity solutions for natural convection.'

end program natural_convection
