program turbine_blade_film_cooling
    implicit none
    integer :: iostat_val
    double precision :: D_hole_mm, P_hole_mm, alpha_deg, x_dist_mm, Tw_target_C
    double precision :: Tinf_gas_C, Uinf_gas_ms, rho_gas, Tc_coolant_C, Uc_coolant_ms, rho_coolant, h0_conv
    double precision :: D_m, x_m, Tinf_K, Tc_K, Tw_K, DR_ratio, M_blowing, I_momentum
    double precision :: eta_center, eta_lateral, Taw_C, q_flux_cooled_kWm2, q_flux_uncooled_kWm2, Flux_reduction_pct

    ! Read inputs
    read(*,*,iostat=iostat_val) D_hole_mm       ! Hole Diameter [mm] (e.g. 1.5)
    read(*,*,iostat=iostat_val) P_hole_mm       ! Hole Pitch P [mm] (e.g. 4.5)
    read(*,*,iostat=iostat_val) alpha_deg       ! Hole Inclination Angle [deg] (e.g. 30.0)
    read(*,*,iostat=iostat_val) x_dist_mm       ! Downstream Distance x [mm] (e.g. 30.0)
    read(*,*,iostat=iostat_val) Tw_target_C     ! Blade Internal Metal Temp [deg C] (e.g. 800.0)
    read(*,*,iostat=iostat_val) Tinf_gas_C      ! Hot Mainstream Gas Temp [deg C] (e.g. 1450.0)
    read(*,*,iostat=iostat_val) Uinf_gas_ms     ! Hot Gas Velocity [m/s] (e.g. 250.0)
    read(*,*,iostat=iostat_val) rho_gas         ! Hot Gas Density [kg/m3] (e.g. 2.0)
    read(*,*,iostat=iostat_val) Tc_coolant_C    ! Coolant Air Temp [deg C] (e.g. 550.0)
    read(*,*,iostat=iostat_val) Uc_coolant_ms   ! Coolant Jet Velocity [m/s] (e.g. 180.0)
    read(*,*,iostat=iostat_val) rho_coolant     ! Coolant Density [kg/m3] (e.g. 4.2)
    read(*,*,iostat=iostat_val) h0_conv         ! Baseline HTC h0 [W/(m2.K)] (e.g. 1200.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for turbine film cooling calculation.'
        stop
    end if

    if (D_hole_mm <= 0.0d0 .or. x_dist_mm < 0.0d0 .or. Uinf_gas_ms <= 0.0d0) then
        write(*,*) 'ERROR: Hole diameter, distance, and velocities must be positive.'
        stop
    end if

    D_m = D_hole_mm * 1.0d-3
    x_m = max(1.0d-4, x_dist_mm * 1.0d-3)
    Tinf_K = Tinf_gas_C + 273.15d0
    Tc_K = Tc_coolant_C + 273.15d0
    Tw_K = Tw_target_C + 273.15d0

    ! Dimensionless Parameters
    DR_ratio = rho_coolant / rho_gas
    M_blowing = (rho_coolant * Uc_coolant_ms) / (rho_gas * Uinf_gas_ms)
    I_momentum = (rho_coolant * (Uc_coolant_ms**2)) / (rho_gas * (Uinf_gas_ms**2))

    ! Goldstein & Baldauf Film Cooling Effectiveness Correlation
    if (I_momentum < 0.8d0) then
        ! Attached Jet Regime
        eta_center = 1.0d0 / (1.0d0 + 0.329d0 * ((x_dist_mm / (max(0.1d0, M_blowing) * D_hole_mm))**0.8d0))
    else
        ! Detached / Lifted Jet Regime with reattachment
        eta_center = (0.75d0 / sqrt(I_momentum)) / (1.0d0 + 0.45d0 * ((x_dist_mm / (max(0.1d0, M_blowing) * D_hole_mm))**0.75d0))
    end if
    if (eta_center > 1.0d0) eta_center = 1.0d0

    ! Lateral Average Film Effectiveness
    eta_lateral = eta_center * (D_hole_mm / P_hole_mm) * 1.8d0
    if (eta_lateral > eta_center) eta_lateral = eta_center

    Taw_C = Tinf_gas_C - eta_lateral * (Tinf_gas_C - Tc_coolant_C)
    q_flux_cooled_kWm2 = (h0_conv * (Taw_C - Tw_target_C)) / 1000.0d0
    q_flux_uncooled_kWm2 = (h0_conv * (Tinf_gas_C - Tw_target_C)) / 1000.0d0
    Flux_reduction_pct = (1.0d0 - q_flux_cooled_kWm2 / max(0.01d0, q_flux_uncooled_kWm2)) * 100.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TURBINE BLADE FILM COOLING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Hole Diameter / Pitch     = ', D_hole_mm, ' mm / ', P_hole_mm, ' mm (x = ', x_dist_mm, ' mm)'
    write(*,'(A,F10.2,A,F10.2)')   'Blowing Ratio M / Mom. I  = ', M_blowing, ' / ', I_momentum
    write(*,'(A,F10.1,A,F10.1,A)') 'Gas Temp / Coolant Temp   = ', Tinf_gas_C, ' deg C / ', Tc_coolant_C, ' deg C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4,A,F10.4)')   'FILM EFFECTIVENESS eta    = ', eta_lateral, ' (Centerline eta_0 = ', eta_center, ')'
    write(*,'(A,F10.1,A)')  'ADIABATIC WALL TEMP Taw   = ', Taw_C, ' deg C'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cooled / Uncooled Heat Flux= ', q_flux_cooled_kWm2, ' kW/m2 / ', q_flux_uncooled_kWm2, ' kW/m2'
    write(*,'(A,F10.2,A)')  'Heat Flux Reduction       = ', Flux_reduction_pct, ' %'
    write(*,'(A)') '============================================================'

end program turbine_blade_film_cooling