program psv_sizing_api
    implicit none
    integer :: svc_type, iostat_val, i, selected_idx
    double precision :: P_set_bar, P_back_bar, overpressure_pct, T_degC
    double precision :: W_kgh, MW, k_ratio, Z_factor, rho_kgm3, mu_cP
    double precision :: P1_barA, Pb_barA, T_K, P_crit_ratio, r_crit, r_actual
    double precision :: C_factor, Kb, Kw, Kd, Kv, A_req_mm2, W_cap_kgh
    double precision :: std_areas(14)
    character(len=2) :: std_letters(14)
    character(len=2) :: sel_letter
    double precision :: sel_area

    ! API 526 Standard Orifices (Area in mm2)
    data std_letters /'D ', 'E ', 'F ', 'G ', 'H ', 'J ', 'K ', 'L ', 'M ', 'N ', 'P ', 'Q ', 'R ', 'T '/
    data std_areas   /71.0d0, 126.0d0, 198.0d0, 325.0d0, 506.0d0, 830.0d0, 1186.0d0, &
                      1841.0d0, 2323.0d0, 2800.0d0, 4116.0d0, 7129.0d0, 10323.0d0, 16774.0d0/

    ! Read inputs
    read(*,*,iostat=iostat_val) svc_type         ! 1=Gas/Vapor, 2=Steam, 3=Liquid
    read(*,*,iostat=iostat_val) P_set_bar        ! Set Pressure [bar gauge]
    read(*,*,iostat=iostat_val) P_back_bar       ! Backpressure [bar gauge]
    read(*,*,iostat=iostat_val) overpressure_pct ! Overpressure [%] (e.g. 10.0, 16.0, 21.0)
    read(*,*,iostat=iostat_val) T_degC           ! Relieving Temperature [deg C]
    read(*,*,iostat=iostat_val) W_kgh            ! Required Mass Relief Rate [kg/h]
    read(*,*,iostat=iostat_val) MW               ! Molecular Weight [g/mol]
    read(*,*,iostat=iostat_val) k_ratio          ! Isentropic exp k = Cp/Cv
    read(*,*,iostat=iostat_val) Z_factor         ! Compressibility factor Z
    read(*,*,iostat=iostat_val) rho_kgm3         ! Liquid density [kg/m3] (for svc=3)
    read(*,*,iostat=iostat_val) mu_cP            ! Liquid viscosity [cP] (for svc=3)

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

    ! Conversions
    P1_barA = (P_set_bar * (1.0d0 + overpressure_pct / 100.0d0)) + 1.01325d0
    Pb_barA = P_back_bar + 1.01325d0
    T_K = T_degC + 273.15d0

    if (svc_type == 1) then
        ! Gas / Vapor Sizing (API 520 Part I)
        Kd = 0.975d0
        Kb = 1.0d0
        if (Pb_barA / P1_barA > 0.5d0) then
            Kb = max(0.65d0, 1.0d0 - 0.7d0 * (Pb_barA / P1_barA - 0.5d0))
        end if

        ! Gas expansion coefficient C
        C_factor = 0.03948d0 * sqrt(k_ratio * ((2.0d0 / (k_ratio + 1.0d0))**((k_ratio + 1.0d0) / (k_ratio - 1.0d0))))
        r_crit = (2.0d0 / (k_ratio + 1.0d0))**(k_ratio / (k_ratio - 1.0d0))
        r_actual = Pb_barA / P1_barA

        if (r_actual <= r_crit) then
            ! Choked / Critical Flow
            A_req_mm2 = (W_kgh / (C_factor * Kd * P1_barA * Kb * 100.0d0)) * sqrt((T_K * Z_factor) / MW)
        else
            ! Subcritical Flow
            A_req_mm2 = (W_kgh / (C_factor * Kd * P1_barA * Kb * 100.0d0)) * sqrt((T_K * Z_factor) / MW) / &
                        max(0.2d0, sqrt(1.0d0 - ((r_actual - r_crit) / (1.0d0 - r_crit))**2))
        end if

    else if (svc_type == 2) then
        ! Steam Sizing (API 520 / ASME Sec VIII)
        Kd = 0.975d0
        Kb = 1.0d0
        A_req_mm2 = W_kgh / (51.5d0 * Kd * P1_barA * Kb)

    else
        ! Liquid Sizing (API 520 Liquid Formula)
        Kd = 0.65d0
        Kw = 1.0d0
        Kv = 1.0d0
        if (P1_barA <= Pb_barA) then
            write(*,*) 'ERROR: Relieving pressure must exceed backpressure.'
            stop
        end if
        A_req_mm2 = (W_kgh / (11.78d0 * Kd * Kw * Kv * sqrt(max(0.01d0, rho_kgm3 * (P1_barA - Pb_barA)))))
    end if

    ! Select API 526 Orifice Letter
    selected_idx = 14
    do i = 1, 14
        if (std_areas(i) >= A_req_mm2) then
            selected_idx = i
            exit
        end if
    end do
    sel_letter = std_letters(selected_idx)
    sel_area = std_areas(selected_idx)
    W_cap_kgh = W_kgh * (sel_area / max(1.0d-3, A_req_mm2))

    ! Output Formatted Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — API 520 / 526 PSV SIZING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,I2)')       'Service Fluid Type        = ', svc_type
    write(*,'(A,F10.2,A)')  'Set Pressure              = ', P_set_bar, ' bar(g)'
    write(*,'(A,F10.2,A)')  'Relieving Pressure P1     = ', P1_barA, ' bar(a)'
    write(*,'(A,F10.2,A)')  'Backpressure Pb           = ', Pb_barA, ' bar(a)'
    write(*,'(A,F10.2,A)')  'Relieving Temperature     = ', T_degC, ' C'
    write(*,'(A,F10.2,A)')  'Required Relief Rate W    = ', W_kgh, ' kg/h'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Calculated Required Area  = ', A_req_mm2, ' mm2'
    write(*,'(A,F10.4,A)')  'Required Area in Sq In    = ', A_req_mm2 / 645.16d0, ' sq.in'
    write(*,'(A,A,A,F10.2,A)') 'Selected API 526 Orifice  = Letter [', sel_letter, '] Area = ', sel_area, ' mm2'
    write(*,'(A,F10.2,A)')  'Rated Orifice Capacity    = ', W_cap_kgh, ' kg/h'
    write(*,'(A,F10.2,A)')  'Overcapacity Margin       = ', ((sel_area - A_req_mm2) / A_req_mm2) * 100.0d0, ' %'
    write(*,'(A)') '============================================================'

end program psv_sizing_api