program npsh_cavitation_pump_safety
    implicit none
    integer :: iostat_val, fluid_type
    double precision :: P_vessel_bar, T_fluid_C, z_suction_m, Q_flow_m3h, D_pipe_mm, L_suction_m
    double precision :: sum_K_fittings, NPSH_r_m, H_pump_m, roughness_mm
    double precision :: rho, mu, Pv_Pa, g, D_m, A_pipe, V_suction, Re_D, f_darcy, h_friction_m, h_minor_m, h_loss_tot
    double precision :: P_vessel_Pa, NPSH_a_m, margin_m, margin_ratio, thoma_sigma
    character(len=32) :: cav_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type       ! 1=Water, 2=Hydrocarbon Gasoline, 3=Ethanol, 4=Hot Boiler Feed
    read(*,*,iostat=iostat_val) P_vessel_bar     ! Vessel / Atmospheric Pressure [bar abs] (e.g. 1.013)
    read(*,*,iostat=iostat_val) T_fluid_C        ! Fluid Operating Temperature [deg C] (e.g. 25.0)
    read(*,*,iostat=iostat_val) z_suction_m      ! Static Suction Head (+ flooded, - lift) [m] (e.g. 1.5)
    read(*,*,iostat=iostat_val) Q_flow_m3h       ! Pump Flow Rate [m3/h] (e.g. 80.0)
    read(*,*,iostat=iostat_val) D_pipe_mm        ! Suction Pipe Inner Diameter [mm] (e.g. 100.0)
    read(*,*,iostat=iostat_val) L_suction_m      ! Suction Pipe Length [m] (e.g. 6.0)
    read(*,*,iostat=iostat_val) sum_K_fittings   ! Minor Loss Coefficient Sum K (e.g. 2.2)
    read(*,*,iostat=iostat_val) NPSH_r_m         ! Pump Required NPSHr [m] (e.g. 2.8)
    read(*,*,iostat=iostat_val) H_pump_m         ! Pump Total Head [m] (e.g. 35.0)

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

    if (D_pipe_mm <= 0.0d0 .or. Q_flow_m3h <= 0.0d0 .or. P_vessel_bar <= 0.0d0) then
        write(*,*) 'ERROR: Pipe diameter, flow rate and pressure must be positive.'
        stop
    end if

    g = 9.80665d0
    roughness_mm = 0.045d0 ! Commercial steel

    ! Fluid properties & Vapor Pressure
    if (fluid_type == 1 .or. fluid_type == 4) then ! Water / Boiler Feed
        rho = 1000.0d0 - 0.003d0 * (T_fluid_C - 4.0d0)**2
        mu = 0.00179d0 / (1.0d0 + 0.03368d0 * T_fluid_C + 0.000221d0 * T_fluid_C**2)
        Pv_Pa = 10.0d0**(8.07131d0 - (1730.63d0 / (233.426d0 + T_fluid_C))) * 133.322d0
    else if (fluid_type == 2) then ! Hydrocarbon Gasoline
        rho = 740.0d0; mu = 0.0006d0
        Pv_Pa = 10.0d0**(6.95464d0 - (1170.966d0 / (224.366d0 + T_fluid_C))) * 100000.0d0 / 760.0d0
    else ! Ethanol
        rho = 789.0d0; mu = 0.0011d0
        Pv_Pa = 10.0d0**(8.04494d0 - (1554.3d0 / (222.65d0 + T_fluid_C))) * 133.322d0
    end if

    D_m = D_pipe_mm / 1000.0d0
    A_pipe = PI * (D_m**2) / 4.0d0
    V_suction = (Q_flow_m3h / 3600.0d0) / A_pipe
    Re_D = rho * V_suction * D_m / mu

    ! Swamee-Jain friction factor
    if (Re_D < 2300.0d0) then
        f_darcy = 64.0d0 / max(1.0d0, Re_D)
    else
        f_darcy = 0.25d0 / (log10((roughness_mm / (3.7d0 * D_pipe_mm)) + (5.74d0 / (Re_D**0.9d0)))**2)
    end if

    h_friction_m = f_darcy * (L_suction_m / D_m) * (V_suction**2) / (2.0d0 * g)
    h_minor_m = sum_K_fittings * (V_suction**2) / (2.0d0 * g)
    h_loss_tot = h_friction_m + h_minor_m

    P_vessel_Pa = P_vessel_bar * 1.0d5
    NPSH_a_m = ((P_vessel_Pa - Pv_Pa) / (rho * g)) + z_suction_m - h_loss_tot

    margin_m = NPSH_a_m - NPSH_r_m
    margin_ratio = NPSH_a_m / max(0.1d0, NPSH_r_m)
    thoma_sigma = NPSH_a_m / max(1.0d0, H_pump_m)

    if (margin_m >= 1.0d0 .and. margin_ratio >= 1.3d0) then
        cav_status = 'SAFE (NO CAVITATION)'
    else if (margin_m >= 0.5d0) then
        cav_status = 'ACCEPTABLE (HI MARGIN)'
    else if (margin_m >= 0.0d0) then
        cav_status = 'INCIPIENT CAVITATION'
    else
        cav_status = 'SEVERE CAVITATION (DAMAGING)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CENTRIFUGAL PUMP NPSHa & CAVITATION'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Flow Rate / Suction Velocity = ', Q_flow_m3h, ' m3/h / ', V_suction, ' m/s'
    write(*,'(A,F10.3,A,F10.1,A)') 'Vapor Pressure / Fluid Temp  = ', Pv_Pa/1000.0d0, ' kPa / ', T_fluid_C, ' deg C'
    write(*,'(A,F10.2,A,F10.2,A)') 'Static Head / Suction Loss   = ', z_suction_m, ' m / ', h_loss_tot, ' m'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'NPSH AVAILABLE (NPSHa)     = ', NPSH_a_m, ' m'
    write(*,'(A,F10.2,A)')  'NPSH Required (NPSHr)      = ', NPSH_r_m, ' m'
    write(*,'(A,F10.2,A)')  'CAVITATION MARGIN (NPSHa - NPSHr) = ', margin_m, ' m'
    write(*,'(A,F10.2)')    'Margin Ratio (NPSHa/NPSHr) = ', margin_ratio
    write(*,'(A,F10.3)')    'Thoma Cavitation Number    = ', thoma_sigma
    write(*,'(A,A)')        'CAVITATION SAFETY STATUS   = ', trim(cav_status)
    write(*,'(A)') '============================================================'

end program npsh_cavitation_pump_safety