Real Gas Compressibility Factor Z

program compressibility_z
    implicit none
    integer::gas,iostat_val,i,ns,j
    double precision::T,P,Tc,Pc,omega,Tc_c,Pc_c,om_c
    double precision::Tr,Pr,Z_pit,Z_vdw,Z_pr,B0,B1
    double precision::a_vdw,b_vdw,Ap,Bp,Z_iter
    double precision::a_pr,b_pr,kappa,alpha_pr,A_pr_coeff,B_pr_coeff
    double precision::P_sw,Z1,Z2,Z3
    double precision,parameter::R=8.314462d0
    character(len=30)::gname
    read(*,*,iostat=iostat_val) gas; read(*,*,iostat=iostat_val) T
    read(*,*,iostat=iostat_val) P; read(*,*,iostat=iostat_val) Tc_c
    read(*,*,iostat=iostat_val) Pc_c; read(*,*,iostat=iostat_val) om_c
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    select case(gas)
    case(1);gname='Nitrogen (N2)';Tc=126.2d0;Pc=3.39d0;omega=0.037d0
    case(2);gname='Oxygen (O2)';Tc=154.6d0;Pc=5.04d0;omega=0.022d0
    case(3);gname='CO2';Tc=304.2d0;Pc=7.38d0;omega=0.224d0
    case(4);gname='Methane (CH4)';Tc=190.6d0;Pc=4.60d0;omega=0.012d0
    case(5);gname='Hydrogen (H2)';Tc=33.2d0;Pc=1.30d0;omega=-0.217d0
    case(6);gname='Air';Tc=132.5d0;Pc=3.77d0;omega=0.035d0
    case default;gname='Custom';Tc=Tc_c;Pc=Pc_c;omega=om_c;gas=7
    end select
    if(Tc<=0)Tc=300.0d0; if(Pc<=0)Pc=5.0d0; if(T<=0)T=300.0d0; if(P<=0)P=1.0d0
    Tr=T/Tc; Pr=P/Pc
    ! Pitzer
    B0=0.083d0-0.422d0/Tr**1.6d0
    B1=0.139d0-0.172d0/Tr**4.2d0
    Z_pit=1.0d0+(B0+omega*B1)*Pr/Tr
    ! Van der Waals
    a_vdw=27.0d0/64.0d0*(R*Tc)**2/(Pc*1.0d6)
    b_vdw=R*Tc/(8.0d0*Pc*1.0d6)
    Ap=a_vdw*P*1.0d6/(R*T)**2
    Bp=b_vdw*P*1.0d6/(R*T)
    Z_vdw=1.0d0
    do j=1,50; Z_vdw=1.0d0+Bp-Ap/(Z_vdw+Bp/Z_vdw); end do
    if(Z_vdw<0.1d0)Z_vdw=0.1d0
    ! Peng-Robinson
    kappa=0.37464d0+1.54226d0*omega-0.26992d0*omega**2
    alpha_pr=(1.0d0+kappa*(1.0d0-sqrt(Tr)))**2
    a_pr=0.45724d0*(R*Tc)**2/(Pc*1.0d6)*alpha_pr
    b_pr=0.07780d0*R*Tc/(Pc*1.0d6)
    A_pr_coeff=a_pr*P*1.0d6/(R*T)**2
    B_pr_coeff=b_pr*P*1.0d6/(R*T)
    Z_pr=1.0d0
    do j=1,80
      Z_pr=1.0d0+B_pr_coeff-(A_pr_coeff-B_pr_coeff-B_pr_coeff**2)/(Z_pr*(Z_pr+B_pr_coeff)+B_pr_coeff*(Z_pr-B_pr_coeff))
      if(Z_pr<0.05d0)Z_pr=0.05d0
    end do
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   COMPRESSIBILITY FACTOR (Z) CALCULATOR'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')       '  Gas                       = ',trim(gname)
    write(*,'(A,F10.2,A)') '  Temperature T             = ',T,' K'
    write(*,'(A,F10.4,A)') '  Pressure P                = ',P,' MPa'
    write(*,'(A,F10.2,A)') '  Tc                        = ',Tc,' K'
    write(*,'(A,F10.4,A)') '  Pc                        = ',Pc,' MPa'
    write(*,'(A,F10.4)')   '  Acentric Factor omega     = ',omega
    write(*,*)
    write(*,'(A)')'--- REDUCED PROPERTIES --------------------------------------'
    write(*,'(A,F10.4)')   '  Tr = T/Tc                 = ',Tr
    write(*,'(A,F10.4)')   '  Pr = P/Pc                 = ',Pr
    write(*,*)
    write(*,'(A)')'--- COMPRESSIBILITY FACTOR Z --------------------------------'
    write(*,'(A,F10.6)')   '  Z (Ideal Gas)             = ',1.0d0
    write(*,'(A,F10.6)')   '  Z (Pitzer Correlation)    = ',Z_pit
    write(*,'(A,F10.6)')   '  Z (Van der Waals)         = ',Z_vdw
    write(*,'(A,F10.6)')   '  Z (Peng-Robinson)         = ',Z_pr
    write(*,*)
    ns=40
    write(*,'(A)')'--- Z VS PRESSURE SWEEP -------------------------------------'
    write(*,'(A)')'  P[MPa]      Z_Pitzer    Z_VdW       Z_PR'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns
      P_sw=0.1d0+dble(i-1)*(50.0d0-0.1d0)/dble(ns-1)
      Pr=P_sw/Pc
      B0=0.083d0-0.422d0/Tr**1.6d0;B1=0.139d0-0.172d0/Tr**4.2d0
      Z1=1.0d0+(B0+omega*B1)*Pr/Tr
      Ap=a_vdw*P_sw*1.0d6/(R*T)**2;Bp=b_vdw*P_sw*1.0d6/(R*T)
      Z2=1.0d0;do j=1,50;Z2=1.0d0+Bp-Ap/(Z2+Bp/Z2);end do;if(Z2<0.05d0)Z2=0.05d0
      A_PR=a_pr*P_sw*1.0d6/(R*T)**2;B_PR=b_pr*P_sw*1.0d6/(R*T)
      Z3=1.0d0;do j=1,80;Z3=1.0d0+B_PR-(A_PR-B_PR-B_PR**2)/(Z3*(Z3+B_PR)+B_PR*(Z3-B_PR));if(Z3<0.05d0)Z3=0.05d0;end do
      write(*,'(F8.3,4X,F10.6,2X,F10.6,2X,F10.6)') P_sw,Z1,Z2,Z3
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  Pitzer: Z=1+(B0+omega*B1)*Pr/Tr'
    write(*,'(A)')'  Van der Waals: iterative cubic solution'
    write(*,'(A)')'  Peng-Robinson: iterative cubic with kappa correlation'
end program compressibility_z