program flange_rating
    implicit none
    integer :: mg,cls,gt,iostat_val,i,ns,ci
    double precision :: tempF,P_max,cf,mf,gm,gy
    double precision :: T_tab(11),P_tab(11),cls_fac(6)
    integer :: cls_vals(6)
    double precision :: T_sw,P_sw,T_C
    character(len=40)::mgname,gname,clname
    T_tab=(/ -20.d0,100.d0,200.d0,300.d0,400.d0,500.d0,600.d0,650.d0,700.d0,750.d0,800.d0 /)
    P_tab=(/ 285.d0,260.d0,230.d0,200.d0,170.d0,140.d0,110.d0,95.d0,80.d0,65.d0,50.d0 /)
    cls_fac=(/ 1.d0,2.63d0,5.26d0,7.89d0,13.16d0,21.93d0 /)
    cls_vals=(/ 150,300,600,900,1500,2500 /)
    read(*,*,iostat=iostat_val) mg; read(*,*,iostat=iostat_val) cls
    read(*,*,iostat=iostat_val) tempF; read(*,*,iostat=iostat_val) gt
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    select case(mg)
    case(1);mgname='1.1 - Carbon Steel (A105)';mf=1.0d0
    case(2);mgname='1.2 - C-1/2Mo (A182-F1)';mf=1.05d0
    case(3);mgname='1.3 - Low Alloy (A182-F11)';mf=1.10d0
    case(4);mgname='2.1 - 304 SS (A182-F304)';mf=1.0d0
    case default;mgname='2.2 - 316 SS (A182-F316)';mf=1.0d0;mg=5
    end select
    ci=1
    do i=1,6; if(cls_vals(i)==cls)ci=i; end do
    cf=cls_fac(ci)
    select case(gt)
    case(1);gname='Spiral Wound (316/graphite)';gm=3.0d0;gy=10000.0d0
    case(2);gname='Ring Joint (RTJ)';gm=5.5d0;gy=18000.0d0
    case(3);gname='Flat Sheet (compressed)';gm=2.0d0;gy=3700.0d0
    case default;gname='PTFE Envelope';gm=2.5d0;gy=6500.0d0;gt=4
    end select
    call interp_PT(tempF,T_tab,P_tab,11,P_max)
    P_max=P_max*cf*mf
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   ASME B16.5 FLANGE PRESSURE-TEMPERATURE RATING'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')       '  Material Group            = ',trim(mgname)
    write(*,'(A,I6)')      '  Pressure Class            = ',cls
    write(*,'(A,F10.1,A)') '  Design Temperature        = ',tempF,' F'
    write(*,'(A,F10.1,A)') '  Design Temperature        = ',(tempF-32.0d0)*5.0d0/9.0d0,' C'
    write(*,'(A,A)')       '  Gasket Type               = ',trim(gname)
    write(*,*)
    write(*,'(A)')'--- RATING FACTORS ------------------------------------------'
    write(*,'(A,F10.4)')   '  Class Factor              = ',cf
    write(*,'(A,F10.4)')   '  Material Factor           = ',mf
    write(*,'(A,F10.4)')   '  Gasket Factor m           = ',gm
    write(*,'(A,F10.1,A)') '  Gasket Seating Stress y   = ',gy,' psi'
    write(*,*)
    write(*,'(A)')'--- MAX ALLOWABLE PRESSURE ----------------------------------'
    write(*,'(A,F12.1,A)') '  P_max (at design T)       = ',P_max,' psig'
    write(*,'(A,F12.1,A)') '  P_max (at design T)       = ',P_max*6.895d0,' kPa_g'
    write(*,'(A,F12.2,A)') '  P_max (at design T)       = ',P_max/14.696d0,' atm'
    write(*,*)
    ns=30
    write(*,'(A)')'--- P-T RATING CURVE ----------------------------------------'
    write(*,'(A)')'  T[F]        T[C]        P_max[psig]   P_max[kPa]'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns
      T_sw=-20.0d0+dble(i-1)*(800.0d0-(-20.0d0))/dble(ns-1)
      call interp_PT(T_sw,T_tab,P_tab,11,P_sw)
      P_sw=P_sw*cf*mf
      T_C=(T_sw-32.0d0)*5.0d0/9.0d0
      write(*,'(F8.1,4X,F8.1,4X,F12.1,4X,F12.1)') T_sw,T_C,P_sw,P_sw*6.895d0
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  ASME B16.5 Table 2 (Class 150, Group 1.1) with scaling.'
    write(*,'(A)')'  Class factors: 150x1, 300x2.63, 600x5.26, etc.'
    write(*,'(A)')'  Gasket m and y per ASME B16.20.'
contains
    subroutine interp_PT(T,Tt,Pt,n,Pout)
      integer,intent(in)::n
      double precision,intent(in)::T,Tt(n),Pt(n)
      double precision,intent(out)::Pout
      integer::j
      double precision::frac
      if(T<=Tt(1))then;Pout=Pt(1);return;end if
      if(T>=Tt(n))then;Pout=Pt(n);return;end if
      do j=1,n-1
        if(T>=Tt(j).and.T<=Tt(j+1))then
          frac=(T-Tt(j))/(Tt(j+1)-Tt(j))
          Pout=Pt(j)+frac*(Pt(j+1)-Pt(j))
          return
        end if
      end do
      Pout=Pt(n)
    end subroutine
end program flange_rating
