program pressurized_vessel_blowdown_transient
    implicit none
    integer :: iostat_val, gas_type
    double precision :: P0_bar, T0_C, V_tank_m3, d_orifice_mm, P_back_bar, Cd_val, P_target_bar
    double precision :: gamma_gas, M_molar, R_univ, T0_K, P0_Pa, P_back_Pa, P_crit_ratio, d_m, A_orifice
    double precision :: m0_gas_kg, P_t, T_t_K, T_min_C, m_dot_0, t_depress_sec, t_depress_min
    character(len=32) :: mdmt_risk
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) gas_type        ! 1=Methane Natural Gas, 2=Nitrogen/Air, 3=Hydrogen, 4=CO2
    read(*,*,iostat=iostat_val) P0_bar          ! Initial Vessel Pressure [bar abs] (e.g. 70.0)
    read(*,*,iostat=iostat_val) T0_C            ! Initial Temperature [deg C] (e.g. 25.0)
    read(*,*,iostat=iostat_val) V_tank_m3       ! Vessel Volume [m3] (e.g. 15.0)
    read(*,*,iostat=iostat_val) d_orifice_mm    ! Blowdown Orifice Diameter [mm] (e.g. 25.0)
    read(*,*,iostat=iostat_val) P_back_bar      ! Backpressure / Flare Header [bar abs] (e.g. 1.2)
    read(*,*,iostat=iostat_val) Cd_val          ! Orifice Discharge Coefficient (e.g. 0.65)
    read(*,*,iostat=iostat_val) P_target_bar    ! Target Depressurization Press [bar abs] (e.g. 6.9 for API 521)

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

    if (P0_bar <= P_back_bar .or. V_tank_m3 <= 0.0d0 .or. d_orifice_mm <= 0.0d0) then
        write(*,*) 'ERROR: Initial pressure must exceed backpressure and volume must be positive.'
        stop
    end if

    R_univ = 8314.46d0 ! J/(kmol.K)
    T0_K = T0_C + 273.15d0
    P0_Pa = P0_bar * 1.0d5
    P_back_Pa = P_back_bar * 1.0d5

    if (gas_type == 1) then ! Methane
        gamma_gas = 1.31d0; M_molar = 16.04d0
    else if (gas_type == 2) then ! Nitrogen/Air
        gamma_gas = 1.40d0; M_molar = 28.96d0
    else if (gas_type == 3) then ! Hydrogen
        gamma_gas = 1.41d0; M_molar = 2.016d0
    else ! CO2
        gamma_gas = 1.28d0; M_molar = 44.01d0
    end if

    d_m = d_orifice_mm / 1000.0d0
    A_orifice = PI * (d_m**2) / 4.0d0

    ! Initial Mass of Gas in Vessel [kg]
    m0_gas_kg = (P0_Pa * V_tank_m3 * M_molar) / (R_univ * T0_K)

    ! Initial Choked Peak Mass Flow [kg/s]
    m_dot_0 = Cd_val * A_orifice * P0_Pa * sqrt((gamma_gas * M_molar) / (R_univ * T0_K)) * &
              ((2.0d0 / (gamma_gas + 1.0d0))**((gamma_gas + 1.0d0) / (2.0d0 * (gamma_gas - 1.0d0))))

    ! Final Temperature at Target Pressure (Isentropic cooling limit)
    T_t_K = T0_K * ((P_target_bar / P0_bar)**((gamma_gas - 1.0d0) / gamma_gas))
    T_min_C = T_t_K - 273.15d0

    ! Blowdown Duration to P_target (API 521 analytical integration)
    t_depress_sec = (2.0d0 * V_tank_m3 / ((gamma_gas - 1.0d0) * Cd_val * A_orifice * sqrt(gamma_gas * R_univ * T0_K / M_molar))) * &
                    (((P0_bar / P_target_bar)**((gamma_gas - 1.0d0) / (2.0d0 * gamma_gas))) - 1.0d0)
    t_depress_min = t_depress_sec / 60.0d0

    ! MDMT Low-Temperature Brittle Fracture Risk
    if (T_min_C < -46.0d0) then
        mdmt_risk = 'HIGH RISK (BELOW -46 C MDMT)'
    else if (T_min_C < -29.0d0) then
        mdmt_risk = 'MODERATE RISK (BELOW -29 C)'
    else
        mdmt_risk = 'LOW RISK (ABOVE -29 C MDMT)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — PRESSURIZED VESSEL TRANSIENT BLOWDOWN'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Initial Pressure / Temp P0, T0= ', P0_bar, ' bar / ', T0_C, ' deg C'
    write(*,'(A,F10.2,A,F10.1,A)') 'Vessel Volume / Orifice Diam  = ', V_tank_m3, ' m3 / ', d_orifice_mm, ' mm'
    write(*,'(A,F10.2,A)')         'Total Gas Inventory Mass      = ', m0_gas_kg, ' kg'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'PEAK INITIAL CHOKED FLOW RATE = ', m_dot_0, ' kg/s'
    write(*,'(A,F10.2,A,F10.1,A)') 'DEPRESSURIZATION TIME (t_blow)= ', t_depress_min, ' min (', t_depress_sec, ' s)'
    write(*,'(A,F10.1,A)')  'MINIMUM FLUID TEMPERATURE     = ', T_min_C, ' deg C'
    write(*,'(A,A)')        'MDMT BRITTLE FRACTURE RISK    = ', trim(mdmt_risk)
    write(*,'(A)') '============================================================'

end program pressurized_vessel_blowdown_transient