program relief_valve
    implicit none
    integer::svc,iostat_val,i,ns,best_j
    double precision::Ps,Pb,T,MW,k,W,rho_l,mu_l,op,Kd_in,Kb_in
    double precision::P1,C_coeff,Kd,Kb,A_req,A_mm2,op_sw,A_sw
    double precision::orifice_A(14)
    character(len=4)::orifice_L(14)
    character(len=30)::sname,sel_letter
    double precision,parameter::pi=3.14159265d0
    orifice_L=(/'D   ','E   ','F   ','G   ','H   ','J   ','K   ','L   ','M   ','N   ','P   ','Q   ','R   ','T   '/)
    orifice_A=(/71.d0,126.d0,198.d0,325.d0,506.d0,830.d0,1186.d0,1841.d0,2323.d0,2800.d0,4116.d0,7126.d0,10323.d0,16774.d0/)
    read(*,*,iostat=iostat_val) svc; read(*,*,iostat=iostat_val) Ps
    read(*,*,iostat=iostat_val) Pb; read(*,*,iostat=iostat_val) T
    read(*,*,iostat=iostat_val) MW; read(*,*,iostat=iostat_val) k
    read(*,*,iostat=iostat_val) W; read(*,*,iostat=iostat_val) rho_l
    read(*,*,iostat=iostat_val) mu_l; read(*,*,iostat=iostat_val) op
    read(*,*,iostat=iostat_val) Kd_in; read(*,*,iostat=iostat_val) Kb_in
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    if(Ps<=0)Ps=700.0d0; if(T<=0)T=400.0d0; if(MW<=0)MW=28.0d0
    if(k<=0.or.k>2)k=1.4d0; if(W<=0)W=5000.0d0
    if(op<=0)op=10.0d0; if(rho_l<=0)rho_l=998.0d0
    select case(svc)
    case(1);sname='Gas / Vapor'
    case(2);sname='Steam'
    case default;sname='Liquid';svc=3
    end select
    ! Absolute relief pressure
    P1=(Ps*(1.0d0+op/100.0d0)+101.325d0)  ! kPa abs
    Kd=Kd_in; Kb=Kb_in
    if(svc==1)then
      if(Kd<=0)Kd=0.975d0; if(Kb<=0)Kb=1.0d0
      C_coeff=0.03948d0*sqrt(k*(2.0d0/(k+1.0d0))**((k+1.0d0)/(k-1.0d0)))
      ! A = W*sqrt(T*Z/(M*k)) / (C*Kd*P1*Kb) — W in kg/h, P in kPa
      A_req=W*sqrt(T/(MW))/(C_coeff*Kd*P1*Kb)
      A_mm2=A_req*1.0d6  ! approx scaling
      ! Better: API formula in mm2
      A_mm2=W/(13160.0d0*Kd*P1/1000.0d0*Kb)*sqrt(T/MW)*1000.0d0
      if(A_mm2<1.0d0)A_mm2=1.0d0
    else if(svc==2)then
      if(Kd<=0)Kd=0.975d0
      A_mm2=W/(52.49d0*Kd*P1/1000.0d0)*10.0d0
      if(A_mm2<1.0d0)A_mm2=1.0d0
    else
      if(Kd<=0)Kd=0.65d0
      A_mm2=W/(38.0d0*Kd)*sqrt(rho_l/max(P1-Pb-101.325d0,1.0d0))*100.0d0
      if(A_mm2<1.0d0)A_mm2=1.0d0
    end if
    ! Select orifice
    sel_letter='T'; best_j=14
    do i=1,14
      if(orifice_A(i)>=A_mm2)then;sel_letter=orifice_L(i);best_j=i;exit;end if
    end do
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   API 520 SAFETY / RELIEF VALVE SIZING'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')       '  Service Type              = ',trim(sname)
    write(*,'(A,F10.2,A)') '  Set Pressure P_set        = ',Ps,' kPa_g'
    write(*,'(A,F10.2,A)') '  Back Pressure P_back      = ',Pb,' kPa_g'
    write(*,'(A,F10.2,A)') '  Temperature T             = ',T,' K'
    write(*,'(A,F10.2,A)') '  Molecular Weight MW       = ',MW,' g/mol'
    write(*,'(A,F10.4)')   '  Isentropic Exponent k     = ',k
    write(*,'(A,F10.2,A)') '  Required Flow W           = ',W,' kg/h'
    write(*,'(A,F10.2,A)') '  Overpressure              = ',op,' percent'
    write(*,'(A,F10.4)')   '  Discharge Coeff Kd        = ',Kd
    write(*,'(A,F10.4)')   '  Back Pressure Factor Kb   = ',Kb
    write(*,*)
    write(*,'(A)')'--- RESULTS -------------------------------------------------'
    write(*,'(A,F12.2,A)') '  Relief Pressure P1 (abs)  = ',P1,' kPa'
    write(*,'(A,F12.2,A)') '  Relief Pressure P1 (abs)  = ',P1/6.895d0,' psia'
    write(*,'(A,F12.2,A)') '  Required Orifice Area     = ',A_mm2,' mm2'
    write(*,'(A,F12.4,A)') '  Required Orifice Area     = ',A_mm2/645.16d0,' in2'
    write(*,'(A,A)')       '  Selected Orifice Letter   = ',trim(sel_letter)
    write(*,'(A,F12.2,A)') '  Selected Orifice Area     = ',orifice_A(best_j),' mm2'
    write(*,*)
    write(*,'(A)')'--- STANDARD API ORIFICE TABLE ------------------------------'
    write(*,'(A)')'  Letter    Area[mm2]     Area[in2]'
    write(*,'(A)')'  -------------------------------------------'
    do i=1,14
      write(*,'(A6,4X,F10.1,4X,F10.4)') orifice_L(i),orifice_A(i),orifice_A(i)/645.16d0
    end do
    write(*,*)
    ns=12
    write(*,'(A)')'--- SENSITIVITY: A VS OVERPRESSURE --------------------------'
    write(*,'(A)')'  OP[%]       P1[kPa]     A_req[mm2]    Orifice'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns
      op_sw=10.0d0+dble(i-1)*11.0d0/dble(ns-1)
      P1=(Ps*(1.0d0+op_sw/100.0d0)+101.325d0)
      if(svc==1)then
        A_sw=W/(13160.0d0*Kd*P1/1000.0d0*Kb)*sqrt(T/MW)*1000.0d0
      else if(svc==2)then
        A_sw=W/(52.49d0*Kd*P1/1000.0d0)*10.0d0
      else
        A_sw=W/(38.0d0*Kd)*sqrt(rho_l/max(P1-Pb-101.325d0,1.0d0))*100.0d0
      end if
      if(A_sw<1.0d0)A_sw=1.0d0
      sel_letter='T'
      do best_j=1,14;if(orifice_A(best_j)>=A_sw)then;sel_letter=orifice_L(best_j);exit;end if;end do
      write(*,'(F6.1,4X,F10.2,4X,F12.2,4X,A)') op_sw,P1,A_sw,trim(sel_letter)
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  API 520 Part I for gas/vapor, steam, and liquid.'
    write(*,'(A)')'  Orifice designations per API 526.'
end program relief_valve
