Pelton Hydro Turbine Runner & Bucket Sizer

program pelton_turbine_runner_bucket_sizing
    implicit none
    integer :: iostat_val, num_nozzles, num_buckets
    double precision :: H_net_m, Q_tot_m3s, N_rpm, Cv_nozzle, ku_speed_ratio, beta2_deg
    double precision :: Vj_ms, Q_jet_m3s, dj_m, dj_mm, U_runner_ms, D_pitch_m, jet_ratio_m
    double precision :: bucket_width_mm, bucket_length_mm, bucket_depth_mm
    double precision :: W_euler_Jkg, eta_hydraulic, eta_overall, P_hyd_MW, P_shaft_MW
    character(len=32) :: pelton_regime
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: G_ACC = 9.80665d0
    double precision, parameter :: RHO_W = 1000.0d0

    ! Read inputs
    read(*,*,iostat=iostat_val) H_net_m         ! Net Available Water Head [m] (e.g. 650.0)
    read(*,*,iostat=iostat_val) Q_tot_m3s       ! Total Water Discharge Q [m3/s] (e.g. 4.50)
    read(*,*,iostat=iostat_val) N_rpm           ! Runner Rotational Speed [rpm] (e.g. 500.0)
    read(*,*,iostat=iostat_val) num_nozzles     ! Number of Injector Nozzles (e.g. 4)
    read(*,*,iostat=iostat_val) Cv_nozzle       ! Nozzle Velocity Coeff Cv (e.g. 0.98)
    read(*,*,iostat=iostat_val) ku_speed_ratio  ! Speed Ratio ku = U/Vj (e.g. 0.465)
    read(*,*,iostat=iostat_val) beta2_deg       ! Bucket Flow Deflection Angle [deg] (e.g. 165.0)

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

    if (H_net_m <= 0.0d0 .or. Q_tot_m3s <= 0.0d0 .or. N_rpm <= 0.0d0 .or. num_nozzles <= 0) then
        write(*,*) 'ERROR: Head, flow rate, speed, and nozzles must be positive.'
        stop
    end if

    ! Water Jet Velocity & Sizing
    Vj_ms = Cv_nozzle * sqrt(2.0d0 * G_ACC * H_net_m)
    Q_jet_m3s = Q_tot_m3s / dble(num_nozzles)
    dj_m = sqrt((4.0d0 * Q_jet_m3s) / (PI * Vj_ms))
    dj_mm = dj_m * 1000.0d0

    ! Runner Wheel Geometry
    U_runner_ms = ku_speed_ratio * Vj_ms
    D_pitch_m = (60.0d0 * U_runner_ms) / (PI * N_rpm)
    jet_ratio_m = D_pitch_m / dj_m

    if (jet_ratio_m < 8.0d0) then
        pelton_regime = 'HEAVILY LOADED / COMPACT'
    else if (jet_ratio_m <= 16.0d0) then
        pelton_regime = 'OPTIMAL HIGH-EFFICIENCY RANGE'
    else
        pelton_regime = 'LIGHTLY LOADED / HIGH HEAD'
    end if

    ! Tygun Formula for Number of Buckets
    num_buckets = int(0.5d0 * jet_ratio_m + 15.0d0 + 0.5d0)

    ! Bucket Dimensions
    bucket_width_mm = 3.0d0 * dj_mm
    bucket_length_mm = 2.6d0 * dj_mm
    bucket_depth_mm = 1.0d0 * dj_mm

    ! Euler Work & Efficiency (95% bucket blade velocity friction)
    W_euler_Jkg = U_runner_ms * (Vj_ms - U_runner_ms) * (1.0d0 + 0.95d0 * cos((180.0d0 - beta2_deg) * PI / 180.0d0))
    eta_hydraulic = W_euler_Jkg / (G_ACC * H_net_m)
    eta_overall = eta_hydraulic * 0.98d0 ! 98% mechanical & windage efficiency

    P_hyd_MW = (RHO_W * G_ACC * Q_tot_m3s * H_net_m) / 1.0d6
    P_shaft_MW = P_hyd_MW * eta_overall

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC -- PELTON RUNNER & BUCKET SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Water Jet Velocity / Jet Diam= ', Vj_ms, ' m/s / ', dj_mm, ' mm'
    write(*,'(A,F10.2,A,F10.1,A)') 'Runner PCD / Jet Ratio (D/dj)= ', D_pitch_m, ' m / ', jet_ratio_m, ''
    write(*,'(A,I4,A,I3,A)')       'Buckets Count / Nozzles Count= ', num_buckets, ' buckets / ', num_nozzles, ' jets'
    write(*,'(A,A)')               'Runner Operating Regime      = ', trim(pelton_regime)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'SHAFT POWER GENERATED        = ', P_shaft_MW, ' MW (Hyd: ', P_hyd_MW, ' MW)'
    write(*,'(A,F10.2,A)')         'Overall Efficiency (eta_o)   = ', eta_overall*100.0d0, ' %'
    write(*,'(A,F10.2,A)')         'Hydraulic Efficiency (eta_h) = ', eta_hydraulic*100.0d0, ' %'
    write(*,'(A,F10.1,A,F10.1,A)') 'Bucket Width / Length / Depth= ', bucket_width_mm, ' mm / ', bucket_length_mm, ' mm'
    write(*,'(A,F10.1,A,F10.2,A)') 'Bucket Depth / Speed Ratio ku= ', bucket_depth_mm, ' mm / ', ku_speed_ratio, ''
    write(*,'(A)') '============================================================'

end program pelton_turbine_runner_bucket_sizing