program water_hammer_allievi_surge
    implicit none
    integer :: iostat_val, pipe_material, fluid_type
    double precision :: L_pipe_m, D_pipe_mm, wall_thick_mm, Q_flow_m3h, P_static_bar, Tc_valve_sec
    double precision :: rho, K_bulk, E_pipe, nu_poisson, D_m, e_m, A_pipe, V0_ms
    double precision :: a_wave_ms, Tp_pipe_sec, dP_joukowsky_bar, dH_joukowsky_m
    double precision :: dP_actual_bar, P_peak_bar, hoop_stress_MPa, allievi_theta
    character(len=32) :: closure_type, safety_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) pipe_material   ! 1=Carbon Steel (E=205 GPa), 2=Ductile Iron (E=170 GPa), 3=HDPE (E=1.0 GPa), 4=PVC (E=3.0 GPa)
    read(*,*,iostat=iostat_val) fluid_type      ! 1=Water (K=2.19 GPa), 2=Crude Oil (K=1.5 GPa), 3=Ethanol (K=1.1 GPa)
    read(*,*,iostat=iostat_val) L_pipe_m        ! Pipeline Length L [m] (e.g. 1500.0)
    read(*,*,iostat=iostat_val) D_pipe_mm       ! Pipe Inner Diameter [mm] (e.g. 300.0)
    read(*,*,iostat=iostat_val) wall_thick_mm   ! Wall Thickness [mm] (e.g. 8.0)
    read(*,*,iostat=iostat_val) Q_flow_m3h      ! Flow Rate [m3/h] (e.g. 500.0)
    read(*,*,iostat=iostat_val) P_static_bar    ! Static Operating Pressure [bar] (e.g. 6.0)
    read(*,*,iostat=iostat_val) Tc_valve_sec    ! Valve Closure Time [s] (e.g. 2.5)

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

    if (L_pipe_m <= 0.0d0 .or. D_pipe_mm <= 0.0d0 .or. wall_thick_mm <= 0.0d0 .or. Q_flow_m3h <= 0.0d0) then
        write(*,*) 'ERROR: Pipe geometry, length, and flow rate must be positive.'
        stop
    end if

    ! Fluid properties
    if (fluid_type == 1) then ! Water
        rho = 998.0d0; K_bulk = 2.19d9
    else if (fluid_type == 2) then ! Oil
        rho = 850.0d0; K_bulk = 1.50d9
    else ! Ethanol
        rho = 789.0d0; K_bulk = 1.10d9
    end if

    ! Pipe material elasticity
    if (pipe_material == 1) then ! Carbon Steel
        E_pipe = 2.05d11; nu_poisson = 0.30d0
    else if (pipe_material == 2) then ! Ductile Iron
        E_pipe = 1.70d11; nu_poisson = 0.28d0
    else if (pipe_material == 3) then ! HDPE
        E_pipe = 1.0d9;   nu_poisson = 0.45d0
    else ! PVC
        E_pipe = 3.0d9;   nu_poisson = 0.38d0
    end if

    D_m = D_pipe_mm / 1000.0d0
    e_m = wall_thick_mm / 1000.0d0
    A_pipe = PI * (D_m**2) / 4.0d0
    V0_ms = (Q_flow_m3h / 3600.0d0) / A_pipe

    ! Acoustic Wave Speed [m/s]
    a_wave_ms = sqrt((K_bulk / rho) / (1.0d0 + (K_bulk / E_pipe) * (D_m / e_m) * (1.0d0 - nu_poisson**2)))

    ! Pipeline Characteristic Period [s]
    Tp_pipe_sec = 2.0d0 * L_pipe_m / a_wave_ms

    ! Joukowsky Maximum Potential Surge (Rapid Closure)
    dP_joukowsky_bar = (rho * a_wave_ms * V0_ms) / 1.0d5
    dH_joukowsky_m = (a_wave_ms * V0_ms) / 9.80665d0

    ! Actual Surge based on Closure Time
    if (Tc_valve_sec <= Tp_pipe_sec) then
        closure_type = 'RAPID CLOSURE (Tc <= 2L/a)'
        dP_actual_bar = dP_joukowsky_bar
    else
        closure_type = 'SLOW CLOSURE (Tc > 2L/a)'
        dP_actual_bar = dP_joukowsky_bar * (Tp_pipe_sec / Tc_valve_sec)
    end if

    P_peak_bar = P_static_bar + dP_actual_bar
    allievi_theta = a_wave_ms * Tc_valve_sec / (2.0d0 * L_pipe_m)

    ! Pipe Hoop Stress [MPa]
    hoop_stress_MPa = ((P_peak_bar * 1.0d5) * D_m) / (2.0d0 * e_m * 1.0d6)

    if (pipe_material == 1 .and. hoop_stress_MPa < 140.0d0) then
        safety_status = 'SAFE (BELOW STEEL ALLOWABLE)'
    else if (pipe_material == 3 .and. hoop_stress_MPa < 10.0d0) then
        safety_status = 'SAFE (BELOW HDPE ALLOWABLE)'
    else if (pipe_material == 4 .and. hoop_stress_MPa < 14.0d0) then
        safety_status = 'SAFE (BELOW PVC ALLOWABLE)'
    else
        safety_status = 'WARNING: HIGH SURGE STRESS'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — WATER HAMMER & SURGE PRESSURE ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2,A)') 'Wave Speed (a) / Flow Velocity= ', a_wave_ms, ' m/s / ', V0_ms, ' m/s'
    write(*,'(A,F10.2,A,F10.2,A)') 'Pipe Period (Tp) / Closure Tc = ', Tp_pipe_sec, ' s / ', Tc_valve_sec, ' s'
    write(*,'(A,A)')               'Closure Regime Classification = ', trim(closure_type)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'MAX JOUKOWSKY SURGE (ΔP_max)  = ', dP_joukowsky_bar, ' bar (', dH_joukowsky_m, ' m)'
    write(*,'(A,F10.2,A)')         'ACTUAL SURGE PRESSURE (ΔP)    = ', dP_actual_bar, ' bar'
    write(*,'(A,F10.2,A)')         'PEAK PIPELINE PRESSURE (Pmax) = ', P_peak_bar, ' bar'
    write(*,'(A,F10.1,A)')         'Pipe Wall Hoop Stress         = ', hoop_stress_MPa, ' MPa'
    write(*,'(A,A)')               'SAFETY INTEGRITY STATUS       = ', trim(safety_status)
    write(*,'(A)') '============================================================'

end program water_hammer_allievi_surge