program hot_tapping_pipeline
    implicit none
    integer :: iostat_val, weld_pass
    double precision :: Do_mm, trun_mm, SMYS_MPa, Pop_bar, Pop_MPa, v_fluid_ms
    double precision :: d_branch_mm, voltage_V, current_A, v_weld_mms, arc_eff
    double precision :: heat_input_kJmm, t_min_safe_mm, T_inner_degC, t85_cooling_s
    character(len=32) :: burn_risk, crack_risk
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Do_mm           ! Pipe Outside Diameter [mm]
    read(*,*,iostat=iostat_val) trun_mm         ! Pipe Wall Thickness [mm]
    read(*,*,iostat=iostat_val) SMYS_MPa        ! Yield Strength SMYS [MPa] (e.g. 358 for X52)
    read(*,*,iostat=iostat_val) Pop_bar         ! Operating Internal Pressure [bar g]
    read(*,*,iostat=iostat_val) v_fluid_ms      ! Internal Fluid Flow Velocity [m/s]
    read(*,*,iostat=iostat_val) d_branch_mm     ! Branch Hot Tap Diameter [mm]
    read(*,*,iostat=iostat_val) voltage_V       ! Welding Arc Voltage [V] (e.g. 22 V)
    read(*,*,iostat=iostat_val) current_A       ! Welding Current [A] (e.g. 110 A)
    read(*,*,iostat=iostat_val) v_weld_mms      ! Travel Speed [mm/s] (e.g. 2.5 mm/s)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs for hot tapping analysis.'
        stop
    end if

    Pop_MPa = Pop_bar * 0.1d0
    arc_eff = 0.80d0 ! SMAW manual stick electrode

    ! Heat Input [kJ/mm]
    heat_input_kJmm = (arc_eff * voltage_V * current_A) / (1000.0d0 * max(0.5d0, v_weld_mms))

    ! Minimum Safe Wall Thickness under ASME B31.8 / API 2201 (Design Factor F=0.72)
    t_min_safe_mm = (Pop_MPa * Do_mm) / (2.0d0 * 0.72d0 * SMYS_MPa)
    t_min_safe_mm = max(4.8d0, t_min_safe_mm) ! API absolute threshold 4.8 mm

    ! Battelle Inner Surface Peak Temperature [deg C]
    T_inner_degC = 20.0d0 + (heat_input_kJmm * 1000.0d0) / &
                   (2.0d0 * PI * 45.0d0 * (trun_mm * 1.0d-3) * (1.0d0 + 0.25d0 * sqrt(v_fluid_ms)))

    if (T_inner_degC > 980.0d0 .or. trun_mm < t_min_safe_mm) then
        burn_risk = 'HIGH BURN-THROUGH HAZARD'
    else if (T_inner_degC > 750.0d0) then
        burn_risk = 'MODERATE BURN-THROUGH RISK'
    else
        burn_risk = 'SAFE AGAINST BURN-THROUGH'
    end if

    ! Cooling time t_8/5 [seconds] (Rykalin 2D/3D heat transfer)
    t85_cooling_s = (6700.0d0 * (heat_input_kJmm**2)) / &
                    ((trun_mm**2) * (1.0d0 + 0.35d0 * v_fluid_ms))
    t85_cooling_s = max(0.5d0, t85_cooling_s)

    if (t85_cooling_s < 3.5d0) then
        crack_risk = 'HIGH HAZ CRACKING RISK (Hard Martensite)'
    else if (t85_cooling_s < 6.0d0) then
        crack_risk = 'MODERATE HAZ HARDNESS'
    else
        crack_risk = 'OPTIMAL DUCTILE MICROSTRUCTURE'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — API RP 2201 HOT TAPPING & WELDING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Pipe Outside Diameter     = ', Do_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Actual Wall Thickness     = ', trun_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Operating Pressure Pop    = ', Pop_bar, ' bar(g)'
    write(*,'(A,F10.2,A)')  'Internal Fluid Velocity   = ', v_fluid_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Welding Heat Input HI     = ', heat_input_kJmm, ' kJ/mm'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Minimum Safe Wall Thk     = ', t_min_safe_mm, ' mm'
    write(*,'(A,F10.1,A)')  'Peak Inner Wall Temp      = ', T_inner_degC, ' C (Limit 980 C)'
    write(*,'(A,A)')        'Burn-Through Safety Check = ', trim(burn_risk)
    write(*,'(A,F10.2,A)')  'Cooling Time t_8/5        = ', t85_cooling_s, ' seconds'
    write(*,'(A,A)')        'HAZ Cold Cracking Risk    = ', trim(crack_risk)
    write(*,'(A)') '============================================================'

end program hot_tapping_pipeline
