program control_valve_noise
    implicit none
    integer :: iostat_val, regime_code
    double precision :: P1_barA, P2_barA, T1_degC, mdot_kgh, MW, gamma_k
    double precision :: Pipe_ID_mm, Pipe_wt_mm, Cv_flow
    double precision :: T1_K, mdot_kgs, r_press, r_crit, U_jet, W_mech, eta_acoust
    double precision :: W_acoust, L_w_dB, TL_dB, Lp_1m_dBA
    double precision, parameter :: W0 = 1.0d-12 ! reference sound power 1 pW
    character(len=32) :: regime_str

    ! Read inputs
    read(*,*,iostat=iostat_val) P1_barA
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inlet pressure.'
        stop
    end if
    read(*,*,iostat=iostat_val) P2_barA
    read(*,*,iostat=iostat_val) T1_degC
    read(*,*,iostat=iostat_val) mdot_kgh
    read(*,*,iostat=iostat_val) MW
    read(*,*,iostat=iostat_val) gamma_k
    read(*,*,iostat=iostat_val) Pipe_ID_mm
    read(*,*,iostat=iostat_val) Pipe_wt_mm

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read control valve noise inputs.'
        stop
    end if

    ! Conversions
    T1_K = T1_degC + 273.15d0
    mdot_kgs = mdot_kgh / 3600.0d0
    r_press = P2_barA / max(0.01d0, P1_barA)
    r_crit = (2.0d0 / (gamma_k + 1.0d0))**(gamma_k / (gamma_k - 1.0d0))

    ! Vena Contracta Jet Velocity & Regimes (IEC 60534-8-3)
    if (r_press > r_crit) then
        regime_code = 1
        regime_str = 'Subcritical Subsonic Flow'
        U_jet = sqrt(2.0d0 * (gamma_k / (gamma_k - 1.0d0)) * (8314.5d0 / MW) * T1_K * (1.0d0 - r_press**((gamma_k - 1.0d0) / gamma_k)))
        eta_acoust = 1.0d-4 * (U_jet / 340.0d0)**3.5d0
    else
        regime_code = 2
        regime_str = 'Choked / Shock Turbulent Flow'
        U_jet = sqrt(2.0d0 * (gamma_k / (gamma_k + 1.0d0)) * (8314.5d0 / MW) * T1_K)
        eta_acoust = 2.5d-4 * (1.0d0 + 0.6d0 * (1.0d0 - r_press / r_crit))
    end if

    ! Mechanical Stream Power & Acoustic Power
    W_mech = 0.5d0 * mdot_kgs * (U_jet**2)
    W_acoust = max(1.0d-12, eta_acoust * W_mech)
    L_w_dB = 10.0d0 * log10(W_acoust / W0)

    ! Pipe Wall Acoustic Transmission Loss (TL)
    TL_dB = 15.0d0 + 10.0d0 * log10(max(1.0d0, Pipe_wt_mm)) + 5.0d0 * log10(max(20.0d0, Pipe_ID_mm))

    ! Sound Pressure Level at 1 meter external distance
    Lp_1m_dBA = L_w_dB - TL_dB - 10.0d0 * log10(2.0d0 * 3.14159d0 * 1.0d0) - 5.0d0
    Lp_1m_dBA = max(20.0d0, Lp_1m_dBA)

    ! Formatted Output
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CONTROL VALVE AERODYNAMIC NOISE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Inlet Pressure P1         = ', P1_barA, ' bar(a)'
    write(*,'(A,F10.2,A)')  'Outlet Pressure P2        = ', P2_barA, ' bar(a)'
    write(*,'(A,F10.2,A)')  'Mass Flow Rate W          = ', mdot_kgh, ' kg/h'
    write(*,'(A,A)')        'Expansion Flow Regime     = ', trim(regime_str)
    write(*,'(A,F10.2,A)')  'Vena Contracta Jet Speed  = ', U_jet, ' m/s'
    write(*,'(A,ES12.4)')   'Acoustic Efficiency Factor= ', eta_acoust
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Internal Sound Power Lw   = ', L_w_dB, ' dB'
    write(*,'(A,F10.2,A)')  'Pipe Transmission Loss TL = ', TL_dB, ' dB'
    write(*,'(A,F10.2,A)')  'External SPL at 1 meter   = ', Lp_1m_dBA, ' dBA'
    write(*,'(A)') '============================================================'

end program control_valve_noise