ASME Pressure Vessel (Sec VIII)

program pressure_vessel_asme
    implicit none
    integer :: iostat_val, mat_idx
    double precision :: P_bar, P_MPa, T_degC, D_mm, R_mm, S_allow_MPa, E_joint, CA_mm
    double precision :: t_hoop_mm, t_long_mm, t_shell_mm
    double precision :: t_ellip_mm, t_tori_mm, t_hemi_mm
    double precision :: mawp_shell_bar, t_corroded_mm, t_nom_mm
    character(len=32) :: mat_name

    ! Read inputs
    read(*,*,iostat=iostat_val) P_bar           ! Design Internal Pressure [bar g]
    read(*,*,iostat=iostat_val) T_degC          ! Design Temperature [deg C]
    read(*,*,iostat=iostat_val) D_mm            ! Inside Diameter [mm]
    read(*,*,iostat=iostat_val) S_allow_MPa     ! Material Allowable Stress S [MPa]
    read(*,*,iostat=iostat_val) E_joint         ! Joint Efficiency E (0.7 - 1.0)
    read(*,*,iostat=iostat_val) CA_mm           ! Corrosion Allowance CA [mm]
    read(*,*,iostat=iostat_val) t_nom_mm        ! Nominal Selected Thickness [mm]

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs for ASME pressure vessel sizing.'
        stop
    end if

    P_MPa = P_bar * 0.1d0
    R_mm = D_mm / 2.0d0

    ! UG-27 Shell Thickness Formulas
    t_hoop_mm = (P_MPa * R_mm) / (S_allow_MPa * E_joint - 0.6d0 * P_MPa) + CA_mm
    t_long_mm = (P_MPa * R_mm) / (2.0d0 * S_allow_MPa * E_joint + 0.4d0 * P_MPa) + CA_mm
    t_shell_mm = max(t_hoop_mm, t_long_mm)

    ! UG-32 Formed Heads Formulas
    t_ellip_mm = (P_MPa * D_mm) / (2.0d0 * S_allow_MPa * E_joint - 0.2d0 * P_MPa) + CA_mm
    t_tori_mm  = (0.885d0 * P_MPa * D_mm) / (S_allow_MPa * E_joint - 0.1d0 * P_MPa) + CA_mm
    t_hemi_mm  = (P_MPa * R_mm) / (2.0d0 * S_allow_MPa * E_joint - 0.2d0 * P_MPa) + CA_mm

    ! MAWP of Cylindrical Shell
    t_corroded_mm = max(0.5d0, t_nom_mm - CA_mm)
    mawp_shell_bar = ((S_allow_MPa * E_joint * t_corroded_mm) / (R_mm + 0.6d0 * t_corroded_mm)) * 10.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ASME SECTION VIII DIV 1 SIZING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Design Pressure P         = ', P_bar, ' bar(g)'
    write(*,'(A,F10.1,A)')  'Inside Diameter D         = ', D_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Allowable Stress S        = ', S_allow_MPa, ' MPa'
    write(*,'(A,F10.2)')    'Joint Efficiency Factor E = ', E_joint
    write(*,'(A,F10.2,A)')  'Corrosion Allowance CA    = ', CA_mm, ' mm'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Required Shell Thickness  = ', t_shell_mm, ' mm'
    write(*,'(A,F10.2,A)')  '2:1 Ellipsoidal Head Thk  = ', t_ellip_mm, ' mm'
    write(*,'(A,F10.2,A)')  'ASME Torispherical Head   = ', t_tori_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Hemispherical Head Thk    = ', t_hemi_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Shell MAWP (Nominal Thk)  = ', mawp_shell_bar, ' bar(g)'
    write(*,'(A)') '============================================================'

end program pressure_vessel_asme