program check_valve_slam
    implicit none
    integer :: valve_type, iostat_val, i, n_pts
    double precision :: D_mm, D_m, L_m, v0_ms, a_wave_ms, rho_kgm3
    double precision :: P_stat_bar, dv_dt_ms2, dP_allow_bar
    double precision :: tc_close_s, v_reverse_ms, dP_slam_bar, P_peak_bar
    double precision :: Area, Vol_pipe_m3, V_air_m3, V_air_L
    character(len=32) :: valve_name
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) D_mm
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid pipe diameter.'
        stop
    end if
    read(*,*,iostat=iostat_val) L_m
    read(*,*,iostat=iostat_val) v0_ms
    read(*,*,iostat=iostat_val) a_wave_ms
    read(*,*,iostat=iostat_val) rho_kgm3
    read(*,*,iostat=iostat_val) P_stat_bar
    read(*,*,iostat=iostat_val) dv_dt_ms2
    read(*,*,iostat=iostat_val) valve_type      ! 1=Swing Check, 2=Dual Plate, 3=Nozzle Non-Slam
    read(*,*,iostat=iostat_val) dP_allow_bar

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

    ! Conversions
    D_m = D_mm * 1.0d-3
    Area = (PI / 4.0d0) * (D_m**2)
    Vol_pipe_m3 = Area * L_m

    ! Dynamic closure response based on valve type
    if (valve_type == 1) then
        valve_name = 'Standard Swing Check Valve'
        tc_close_s = 0.35d0
    else if (valve_type == 2) then
        valve_name = 'Spring-Assisted Dual Plate'
        tc_close_s = 0.12d0
    else
        valve_name = 'Axial Nozzle Non-Slam Valve'
        tc_close_s = 0.035d0
    end if

    ! Maximum reverse velocity before seating
    v_reverse_ms = dv_dt_ms2 * tc_close_s

    ! Joukowsky Slam Surge Pressure
    dP_slam_bar = (rho_kgm3 * a_wave_ms * v_reverse_ms) * 1.0d-5
    P_peak_bar = P_stat_bar + dP_slam_bar

    ! Air vessel volume sizing (Thorley & Parmakian formulation)
    if (dP_allow_bar > 0.01d0) then
        V_air_m3 = (2.0d0 * Area * L_m * (v0_ms**2) * rho_kgm3) / &
                   (2.0d0 * 1.0d5 * dP_allow_bar * (1.0d0 + (P_stat_bar / max(0.5d0, dP_allow_bar))))
        V_air_m3 = max(0.01d0, V_air_m3)
    else
        V_air_m3 = 1.0d0
    end if
    V_air_L = V_air_m3 * 1000.0d0

    ! Formatted Output
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CHECK VALVE SLAM & AIR CHAMBER SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Pipe Diameter             = ', D_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Pipeline Length           = ', L_m, ' m'
    write(*,'(A,F10.2,A)')  'Initial Flow Velocity     = ', v0_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Acoustic Wave Speed a     = ', a_wave_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Deceleration Rate dv/dt   = ', dv_dt_ms2, ' m/s2'
    write(*,'(A,A)')        'Selected Check Valve      = ', trim(valve_name)
    write(*,'(A,F10.4,A)')  'Valve Dynamic Action Time = ', tc_close_s, ' s'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4,A)')  'Maximum Reverse Velocity  = ', v_reverse_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Slam Surge Pressure Spike = ', dP_slam_bar, ' bar'
    write(*,'(A,F10.2,A)')  'Total Peak Pressure       = ', P_peak_bar, ' bar'
    write(*,'(A,F10.2,A)')  'Recommended Air Chamber   = ', V_air_L, ' Liters'
    write(*,'(A)') '============================================================'

end program check_valve_slam
