Wind Turbine BEM & Betz Limit Sizer

program wind_turbine_blade_element_momentum_bem
    implicit none
    integer :: iostat_val, num_blades
    double precision :: D_rotor_m, V_wind_ms, N_rpm, pitch_deg, rho_air_kgm3, Weibull_c_ms
    double precision :: R_blade_m, A_swept_m2, TSR_lambda, U_tip_ms, P_wind_MW
    double precision :: a_induction, a_prime, Cp_aero, Ct_thrust, P_aero_MW, P_elec_MW
    double precision :: Thrust_kN, AEP_GWh, capacity_factor_pct
    character(len=32) :: aero_status
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: CP_BETZ = 16.0d0 / 27.0d0

    ! Read inputs
    read(*,*,iostat=iostat_val) D_rotor_m        ! Rotor Diameter [m] (e.g. 126.0)
    read(*,*,iostat=iostat_val) V_wind_ms        ! Rated Wind Speed [m/s] (e.g. 11.5)
    read(*,*,iostat=iostat_val) N_rpm            ! Rotor Rotational Speed [rpm] (e.g. 12.0)
    read(*,*,iostat=iostat_val) pitch_deg        ! Blade Pitch Angle [deg] (e.g. 0.0)
    read(*,*,iostat=iostat_val) num_blades       ! Number of Blades (e.g. 3)
    read(*,*,iostat=iostat_val) rho_air_kgm3     ! Air Density [kg/m3] (e.g. 1.225)
    read(*,*,iostat=iostat_val) Weibull_c_ms     ! Site Annual Mean Wind Weibull c [m/s] (e.g. 8.5)

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

    if (D_rotor_m <= 0.0d0 .or. V_wind_ms <= 0.0d0 .or. N_rpm <= 0.0d0 .or. num_blades <= 0) then
        write(*,*) 'ERROR: Diameter, wind speed, RPM, and blades count must be positive.'
        stop
    end if

    R_blade_m = D_rotor_m / 2.0d0
    A_swept_m2 = PI * (R_blade_m**2)
    U_tip_ms = 2.0d0 * PI * (N_rpm / 60.0d0) * R_blade_m
    TSR_lambda = U_tip_ms / V_wind_ms

    ! Wind Kinetic Power in Swept Area [MW]
    P_wind_MW = (0.50d0 * rho_air_kgm3 * A_swept_m2 * (V_wind_ms**3)) / 1.0d6

    ! Empirical Heuristic BEM formulation for Modern 3-Blade HAWT
    ! Cp as function of Tip Speed Ratio (TSR) and pitch
    Cp_aero = 0.5176d0 * (116.0d0 / (TSR_lambda + 0.08d0 * pitch_deg - 0.035d0 / (pitch_deg**3 + 1.0d0)) - &
              0.40d0 * pitch_deg - 5.0d0) * exp(-21.0d0 / (TSR_lambda + 0.08d0 * pitch_deg - 0.035d0 / (pitch_deg**3 + 1.0d0))) + &
              0.0068d0 * TSR_lambda

    Cp_aero = max(0.05d0, min(0.495d0, Cp_aero))

    ! Axial Induction Factor estimate from Cp
    a_induction = (1.0d0 - sqrt(max(0.0d0, 1.0d0 - Cp_aero))) / 2.0d0
    a_induction = min(0.40d0, max(0.05d0, a_induction))

    ! Thrust Coefficient Ct (Glauert formulation)
    if (a_induction <= 0.33d0) then
        Ct_thrust = 4.0d0 * a_induction * (1.0d0 - a_induction)
    else
        Ct_thrust = 4.0d0 * (0.33d0**2 + (1.0d0 - 2.0d0 * 0.33d0) * a_induction)
    end if

    P_aero_MW = P_wind_MW * Cp_aero
    P_elec_MW = P_aero_MW * 0.94d0 ! 94% gearbox & generator efficiency
    Thrust_kN = (0.50d0 * rho_air_kgm3 * A_swept_m2 * (V_wind_ms**2) * Ct_thrust) / 1000.0d0

    ! Annual Energy Production (AEP) with Rayleigh / Weibull k=2
    AEP_GWh = P_elec_MW * 8760.0d0 * 0.42d0 * ((Weibull_c_ms / 8.5d0)**2.2d0) / 1000.0d0
    capacity_factor_pct = (AEP_GWh * 1000.0d0) / (P_elec_MW * 8760.0d0) * 100.0d0

    if (Cp_aero > 0.46d0) then
        aero_status = 'NEAR-OPTIMAL BEM EFFICIENCY'
    else if (TSR_lambda > 9.0d0) then
        aero_status = 'HIGH TSR / TIP LOSS PROFILE'
    else
        aero_status = 'STALL / OFF-DESIGN OPERATION'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — WIND TURBINE BEM & BETZ SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Tip Speed Ratio (TSR λ) / Tip= ', TSR_lambda, ' / ', U_tip_ms, ' m/s'
    write(*,'(A,F10.3,A,F10.3,A)') 'Power Coeff Cp / Betz Limit  = ', Cp_aero, ' / ', CP_BETZ, ''
    write(*,'(A,A)')               'Aerodynamic Regime Status    = ', trim(aero_status)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'ELECTRICAL POWER GENERATED   = ', P_elec_MW, ' MW (Aero: ', P_aero_MW, ' MW)'
    write(*,'(A,F10.1,A,F10.3,A)') 'Total Rotor Thrust Load      = ', Thrust_kN, ' kN (Ct = ', Ct_thrust, ')'
    write(*,'(A,F10.2,A,F10.1,A)') 'Estimated AEP / Capacity Fact= ', AEP_GWh, ' GWh/yr / ', capacity_factor_pct, ' %'
    write(*,'(A,F10.1,A,F10.1,A)') 'Rotor Swept Area / Diam      = ', A_swept_m2, ' m2 / ', D_rotor_m, ' m'
    write(*,'(A)') '============================================================'

end program wind_turbine_blade_element_momentum_bem