﻿program flat_plate_convection
    implicit none
    
    ! Variables
    double precision :: L, T_s, T_inf, U_inf, rho, mu, nu, k_f, Pr, Cp
    double precision :: Re_L, Re_x, Re_crit, x_crit
    double precision :: Nu_lam, Nu_turb, Nu_mixed, h_avg
    double precision :: Q_total, A, W
    double precision :: delta_lam, delta_turb, Cf_lam, Cf_turb
    double precision :: delta_t_lam
    character(len=20) :: flow_regime
    integer :: i, n_points
    double precision :: x, dx, local_Re, local_Nu, local_h
    double precision :: local_Cf, local_delta
    
    ! Read inputs
    read(*,*) L        ! Plate length [m]
    read(*,*) W        ! Plate width [m]
    read(*,*) T_s      ! Surface temperature [°C]
    read(*,*) T_inf    ! Free stream temperature [°C]
    read(*,*) U_inf    ! Free stream velocity [m/s]
    read(*,*) rho      ! Fluid density [kg/m³]
    read(*,*) mu       ! Dynamic viscosity [Pa·s]
    read(*,*) k_f      ! Thermal conductivity [W/m·K]
    read(*,*) Pr       ! Prandtl number
    read(*,*) Cp       ! Specific heat [J/kg·K]
    
    ! Calculate properties
    nu = mu / rho
    A = L * W
    Re_crit = 5.0d5
    
    ! Reynolds number at trailing edge
    Re_L = U_inf * L / nu
    
    ! Determine flow regime
    if (Re_L < Re_crit) then
        flow_regime = 'Laminar'
    else
        flow_regime = 'Mixed (Lam+Turb)'
    end if
    
    ! Critical length (transition point)
    x_crit = Re_crit * nu / U_inf
    if (x_crit > L) x_crit = L
    
    ! ============================================
    ! AVERAGE NUSSELT NUMBER CALCULATIONS
    ! ============================================
    
    ! Laminar correlation (entire plate laminar): Nu = 0.664 Re^0.5 Pr^(1/3)
    Nu_lam = 0.664d0 * Re_L**0.5d0 * Pr**(1.0d0/3.0d0)
    
    ! Turbulent correlation (entire plate turbulent): Nu = 0.037 Re^0.8 Pr^(1/3)
    Nu_turb = 0.037d0 * Re_L**0.8d0 * Pr**(1.0d0/3.0d0)
    
    ! Mixed correlation: Nu = (0.037 Re^0.8 - 871) Pr^(1/3)
    if (Re_L > Re_crit) then
        Nu_mixed = (0.037d0 * Re_L**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
    else
        Nu_mixed = Nu_lam
    end if
    
    ! Average heat transfer coefficient
    h_avg = Nu_mixed * k_f / L
    
    ! Total heat transfer
    Q_total = h_avg * A * abs(T_s - T_inf)
    
    ! Boundary layer thicknesses at trailing edge
    if (Re_L < Re_crit) then
        delta_lam = 5.0d0 * L / sqrt(Re_L)
        delta_turb = 0.0d0
        Cf_lam = 1.328d0 / sqrt(Re_L)
        Cf_turb = 0.0d0
    else
        delta_lam = 5.0d0 * x_crit / sqrt(Re_crit)
        delta_turb = 0.37d0 * L / Re_L**0.2d0
        Cf_lam = 1.328d0 / sqrt(Re_crit)
        Cf_turb = 0.074d0 / Re_L**0.2d0
    end if
    
    ! Thermal boundary layer thickness (laminar)
    delta_t_lam = delta_lam / Pr**(1.0d0/3.0d0)
    
    ! ============================================
    ! OUTPUT RESULTS
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   FORCED CONVECTION - FLOW OVER A FLAT PLATE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ----------------------------------------'
    write(*,'(A,F12.4,A)') '  Plate Length (L)        = ', L, ' m'
    write(*,'(A,F12.4,A)') '  Plate Width (W)         = ', W, ' m'
    write(*,'(A,F12.4,A)') '  Surface Area (A)        = ', A, ' m²'
    write(*,'(A,F12.2,A)') '  Surface Temperature     = ', T_s, ' °C'
    write(*,'(A,F12.2,A)') '  Free Stream Temperature = ', T_inf, ' °C'
    write(*,'(A,F12.4,A)') '  Free Stream Velocity    = ', U_inf, ' m/s'
    write(*,*)
    write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Density (ρ)             = ', rho, ' kg/m³'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity (μ)   = ', mu, ' Pa·s'
    write(*,'(A,ES12.4,A)') '  Kinematic Viscosity (ν) = ', nu, ' m²/s'
    write(*,'(A,F12.4,A)')  '  Thermal Conductivity    = ', k_f, ' W/m·K'
    write(*,'(A,F12.4)')    '  Prandtl Number (Pr)     = ', Pr
    write(*,*)
    write(*,'(A)') '--- FLOW ANALYSIS -------------------------------------------'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re_L)  = ', Re_L
    write(*,'(A,ES12.4)')   '  Critical Re             = ', Re_crit
    write(*,'(A,A)')        '  Flow Regime             = ', trim(flow_regime)
    write(*,'(A,F12.4,A)')  '  Transition Point (x_cr) = ', x_crit, ' m'
    write(*,*)
    write(*,'(A)') '--- HEAT TRANSFER RESULTS -----------------------------------'
    write(*,'(A,F12.4)')    '  Average Nusselt Number  = ', Nu_mixed
    write(*,'(A,F12.4,A)')  '  Average h               = ', h_avg, ' W/m²·K'
    write(*,'(A,F12.4,A)')  '  Total Heat Transfer (Q) = ', Q_total, ' W'
    write(*,*)
    write(*,'(A)') '--- BOUNDARY LAYER ------------------------------------------'
    if (Re_L < Re_crit) then
        write(*,'(A,F12.6,A)') '  δ at trailing edge      = ', delta_lam*1000, ' mm'
        write(*,'(A,F12.6,A)') '  δ_thermal at trailing    = ', delta_t_lam*1000, ' mm'
        write(*,'(A,ES12.4)')  '  Average Cf (laminar)    = ', Cf_lam
    else
        write(*,'(A,F12.6,A)') '  δ_lam at transition     = ', delta_lam*1000, ' mm'
        write(*,'(A,F12.6,A)') '  δ_turb at trailing edge = ', delta_turb*1000, ' mm'
        write(*,'(A,ES12.4)')  '  Cf_lam (at transition)  = ', Cf_lam
        write(*,'(A,ES12.4)')  '  Cf_turb (average)       = ', Cf_turb
    end if
    write(*,*)
    
    ! ============================================
    ! LOCAL PROFILE DATA (for chart)
    ! ============================================
    write(*,'(A)') '--- LOCAL PROFILE ALONG PLATE --------------------------------'
    write(*,'(A)') '  x [m]       Re_x          Nu_x         h_x [W/m2.K]  d [mm]'
    write(*,'(A)') '  ------------------------------------------------------------'
    
    n_points = 40
    dx = L / dble(n_points)
    
    do i = 1, n_points
        x = dble(i) * dx
        local_Re = U_inf * x / nu
        
        if (local_Re < Re_crit) then
            ! Laminar
            local_Nu = 0.332d0 * sqrt(local_Re) * Pr**(1.0d0/3.0d0)
            local_delta = 5.0d0 * x / sqrt(local_Re)
            local_Cf = 0.664d0 / sqrt(local_Re)
        else
            ! Turbulent
            local_Nu = 0.0296d0 * local_Re**0.8d0 * Pr**(1.0d0/3.0d0)
            local_delta = 0.37d0 * x / local_Re**0.2d0
            local_Cf = 0.0592d0 / local_Re**0.2d0
        end if
        
        local_h = local_Nu * k_f / x
        
        write(*,'(F8.4,2X,ES12.4,2X,F10.2,4X,F10.4,4X,F10.4)') &
            x, local_Re, local_Nu, local_h, local_delta*1000
    end do
    
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Laminar:    Nu_x = 0.332 Re_x^0.5 Pr^(1/3)'
    write(*,'(A)') '  Turbulent:  Nu_x = 0.0296 Re_x^0.8 Pr^(1/3)'
    write(*,'(A)') '  Average:    Nu_L = (0.037 Re_L^0.8 - 871) Pr^(1/3)'
    write(*,'(A)') '  Valid for:  0.6 < Pr < 60, Re_crit = 5×10⁵'
    
end program flat_plate_convection
