program pipeline_linepack
    implicit none
    integer :: iostat_val
    double precision :: L_km, Di_mm, P1_barA, P2_barA, T_degC, P_min_barA
    double precision :: MW, Z_avg, Q_demand_Sm3h
    double precision :: L_m, Di_m, V_geom_m3, P_avg_barA, T_K
    double precision :: V_std_total_Sm3, V_std_min_Sm3, V_usable_Sm3
    double precision :: Mass_pack_tons, buffer_hours
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: P_std = 1.01325d0, T_std = 288.15d0, Ru = 8314.5d0

    ! Read inputs
    read(*,*,iostat=iostat_val) L_km           ! Pipeline Length [km]
    read(*,*,iostat=iostat_val) Di_mm          ! Pipeline Inside Diameter [mm]
    read(*,*,iostat=iostat_val) P1_barA        ! Inlet Pressure P1 [bar a]
    read(*,*,iostat=iostat_val) P2_barA        ! Outlet Pressure P2 [bar a]
    read(*,*,iostat=iostat_val) T_degC         ! Average Gas Temperature [deg C]
    read(*,*,iostat=iostat_val) P_min_barA     ! Minimum Allowable Delivery P [bar a]
    read(*,*,iostat=iostat_val) MW             ! Gas Molecular Weight [g/mol] (Methane ~16.04)
    read(*,*,iostat=iostat_val) Z_avg          ! Average Compressibility Factor Z
    read(*,*,iostat=iostat_val) Q_demand_Sm3h  ! Consumer Draw Demand Rate [Sm3/h]

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

    if (P1_barA < P2_barA) then
        write(*,*) 'ERROR: Inlet pressure P1 must be >= Outlet pressure P2.'
        stop
    end if

    ! Conversions
    L_m = L_km * 1000.0d0
    Di_m = Di_mm * 1.0d-3
    T_K = T_degC + 273.15d0
    V_geom_m3 = (PI / 4.0d0) * (Di_m**2) * L_m

    ! AGA Average Linepack Pressure
    if (abs(P1_barA - P2_barA) < 1.0d-4) then
        P_avg_barA = P1_barA
    else
        P_avg_barA = (2.0d0 / 3.0d0) * (P1_barA + (P2_barA**2) / (P1_barA + P2_barA))
    end if

    ! Total Standard Gas Inventory [Sm3]
    V_std_total_Sm3 = V_geom_m3 * (P_avg_barA / P_std) * (T_std / T_K) * (1.0d0 / Z_avg)
    Mass_pack_tons = (V_std_total_Sm3 * (P_std * 1.0d5) * (MW * 1.0d-3)) / (Ru * T_std * 1000.0d0)

    ! Minimum inventory at delivery threshold
    V_std_min_Sm3 = V_geom_m3 * (P_min_barA / P_std) * (T_std / T_K) * (1.0d0 / Z_avg)
    V_usable_Sm3 = max(0.0d0, V_std_total_Sm3 - V_std_min_Sm3)

    if (Q_demand_Sm3h > 1.0d-3) then
        buffer_hours = V_usable_Sm3 / Q_demand_Sm3h
    else
        buffer_hours = 0.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — PIPELINE LINEPACK INVENTORY ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Pipeline Length           = ', L_km, ' km'
    write(*,'(A,F10.2,A)')  'Pipeline Diameter Di      = ', Di_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Geometric Volume          = ', V_geom_m3, ' m3'
    write(*,'(A,F10.2,A)')  'Average Linepack Pressure = ', P_avg_barA, ' bar(a)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F12.1,A)')  'Total Standard Linepack   = ', V_std_total_Sm3, ' Sm3'
    write(*,'(A,F10.2,A)')  'Total In-Situ Gas Mass    = ', Mass_pack_tons, ' Metric Tons'
    write(*,'(A,F12.1,A)')  'Usable Draft Buffer Stock = ', V_usable_Sm3, ' Sm3'
    write(*,'(A,F10.2,A)')  'Buffer Autonomy Duration  = ', buffer_hours, ' hours'
    write(*,'(A)') '============================================================'

end program pipeline_linepack
