c  begin file prop_sub.f
c
c  This file contains the basic (non-iterative) routines to calculate
c  various properties of fluids and mixtures.  These routines must first
c  be initialized by a call to the subroutine SETUP.
c
c  contained here are:
c     subroutine CRITP (x,tcrit,pcrit,Dcrit,ierr,herr)
c     subroutine THERM (t,rho,x,p,e,h,s,cv,cp,w,hjt)
c     subroutine THERM2 (t,rho,x,p,e,h,s,cv,cp,w,Z,hjt,A,G,xkappa,beta,
c    &                   dPdrho,d2PdD2,dPT,drhodT,drhodP,
c    &                   spare1,spare2,spare3,spare4)
c     subroutine THERM0 (t,rho,x,p0,e0,h0,s0,cv0,cp0,w0,A0,G0)
c     subroutine ENTRO (t,rho,x,s)
c     subroutine ENTHAL (t,rho,x,h)
c     subroutine ENERGY (t,rho,x,e)
c     subroutine CVCP (t,rho,x,cv,cp)
c     subroutine CVCPK (icomp,t,rho,cv,cp)
c     subroutine GIBBS (t,rho,x,Ar,Gr)
c     subroutine AG (t,rho,x,a,g)
c     subroutine PRESS (t,rho,x,p)
c     subroutine DPDD (t,rho,x,dpdrho)
c     subroutine DPDDK (icomp,t,rho,dpdrho)
c     subroutine DPDD2 (t,rho,x,d2PdD2)
c     subroutine DPDT (t,rho,x,dpt)
c     subroutine DPDTK (icomp,t,rho,dpt)
c     subroutine DDDP (t,rho,x,drhodp)
c     subroutine DDDT (t,rho,x,drhodt)
c     subroutine DHD1(t,rho,x,dhdt_d,dhdt_p,dhdd_t,dhdd_p,dhdp_t,dhdp_d)
c     subroutine FGCTY (t,rho,x,f)
c     subroutine ACTVY (t,rho,x,gamma)
c     subroutine VIRB (t,x,b)
c     subroutine DBDT (t,x,dbt)
c     subroutine VIRC (t,x,c)
c     subroutine VIRD (t,x,d)
c     subroutine B12 (t,x,b)
c     subroutine EXCESS (t,p,x,vE,eE,hE,sE)
c     subroutine FPV (t,rho,p,x,f)
c     subroutine RMIX (x)
c
c  these routines use the following common blocks from other files
c     common /MODEL/ hrefst,heos,hpheq,h2eos(n0:nx),hmixp,htran,hsten
c     common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
c     common /NCOMP/ nc,ic
c     common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
c  various arrays are dimensioned with parameter statements
c     parameter (ncmax=20)        !max number of components in mixture
c     parameter (nrefmx=10)       !max number of fluids for transport ECS
c     parameter (n0=-ncmax-nrefmx,nx=ncmax)
c
c ======================================================================
c ======================================================================
c
      subroutine CRITP (x,tcrit,pcrit,Dcrit,ierr,herr)
c
c  critical parameters as a function of composition
c
c  input:
c        x--composition [array of mol frac]
c  outputs:
c    tcrit--critical temperature [K]
c    pcrit--critical pressure [kPa]
c    Dcrit--critical density [mol/L]
c     ierr--error flag:  0 = successful
c                        1 = did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-20-94  MM, original version
c  07-21-95  MM, call CRTBWR instead of accessing arrays directly
c  08-07-95  MM, add call to Fundamental (Helmholtz) EOS
c  09-13-95  MM, add ierr, herr to argument list
c  10-03-95  MM, change /MODEL/--models specified by strings
c  11-02-95  MM, add call mixture Helmholtz model (HMX)
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c                add call to ECS model
c  03-19-19  MM, add dipole moment to /CCON/
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  03-07-07 EWL, add check for tc, pc, dc less than 0
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: CRITP
c     dll_export CRITP
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*1 htab,hnull
      character*3 hpheq,heos,hmxeos,hmodcp
      character*255 herr
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      dimension x(ncmax)
c
      ierr=0
      herr=' '
      if (heos.eq.'FEQ') then
c  pure fluid Fundamental (Helmholtz) EOS
        icomp=1
        call CRTFEQ (icomp,tcrit,pcrit,Dcrit)
      else if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state
        icomp=1
        call CRTBWR (icomp,tcrit,pcrit,Dcrit)
      else if (heos.eq.'ECS') then
c  pure fluid ECS-thermo model
        icomp=1
        call CRTECS (icomp,tcrit,pcrit,Dcrit)
      else if (heos.eq.'HMX') then
c  mixture Helmholtz model
c       write (*,1022) (x(i),i=1,nc)
c1022   format (1x,' CRITP--about to call CRTHMX w/ x = ',5f12.8)
        call CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr)
      else if (heos.eq.'AGA') then
        call CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr)
      else if (heos.eq.'PR') then
        call CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr)
      else
        ierr=1
        herr='[CRITP error] Specified model not found'//hnull
        call ERRMSG (ierr,herr)
      end if
      if (tcrit.le.0) tcrit=300
      if (pcrit.le.0) pcrit=1000
      if (dcrit.le.0) dcrit=10
c
      RETURN
      end                                              !subroutine CRITP
c
c ======================================================================
c
      subroutine THERM (t,rho,x,p,e,h,s,cv,cp,w,hjt)
c
c  compute thermal quantities as a function of temperature, density,
c  and compositions using core functions (Helmholtz free energy, ideal
c  gas heat capacity and various derivatives and integrals)
c
c  Based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994,
c  Appendix A for pressure-explicit equations (e.g. MBWR) and
c  Baehr & Tillner-Roth, Thermodynamic Properties of Environmentally
c  Acceptable Refrigerants, Berlin:  Springer-Verlag (1995) for
c  Helmholtz-explicit equations (e.g. FEQ).
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c        p--pressure [kPa]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [J/mol-K]
c       Cv--isochoric heat capacity [J/mol-K]
c       Cp--isobaric heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c      hjt--isenthalpic Joule-Thompson coefficient [K/kPa]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-11-94  MM, original version
c  08-04-95  MM, add calls to Fundamental (Helmholtz) EOS
c  10-03-95  MM, change /MODEL/--models specified by strings
c  10-10-95  MM, compute ideal gas pressure and pass to PHI0
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-06-95  MM, add calls to mixture ideal gas function
c  11-08-95  MM, convert calls to PHI0, CP0, CPI, CPT to mixture form
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  12-13-95  MM, compute entropy using Cp0, etc rather than PHI0
c  01-18-96  MM, fix s and h ref state for HMX: s = s - sum[x(i)*sref(i)]
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-19-96  MM, add dipole moment to /CCON/
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  04-19-96  MM, change call to PHI0:  pass rho instead of pideal
c                calculate e,h,s using PHI0 rather than Cp0
c  07-05-96  MM, change e, Cv:  PHI0 returns tau*d(phi0)/d(tau), etc.
c  04-22-97  MM, lower bound on rho for s calc set to 1.0d-20
c  10-01-97  MM, add compiler switch to allow access by DLL
c  03-30-98 EWL, add Joule-Thompson coeff to MBWR case
c  03-31-98  MM, special case for Joule-Thompson coeff for rho = 0
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  12-02-98 EWL, remove compositional dependence for pure fluids
c  12-02-98 EWL, restructure to closely mimic THERM2
c  12-22-98 EWL, recalculate R for mixtures based on values for pure fluids
c  05-25-00 EWL, moved calculation of Z to bottom AFTER p is calculated!
c  09-00-00 EWL, removed the del from del*phi01, etc.  The del's and tau's
c                are now included in the core routines.  Put the reducing
c                variables tz and rhoz directly in the common blocks.
c  09-05-02 EWL, change check on rho from 1.d-10 to 1.0d-8 to avoid
c                zero's coming back from core_feq at low rhos.
c  10-04-06 EWL, change remaining checks on rho from 1.d-10 to 1.0d-8
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: THERM
c     dll_export THERM
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      dimension x(ncmax)
c
      call RMIX(x)
      RT=R*t
      w2=0.d0
      if (rho.lt.1.0d-20) then
c  entropy calc will crash if rho = 0
        rhos=1.0d-20
      else
        rhos=rho
      end if
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        a=ABWR(icomp,t,rho)
        dadt=DABWR(icomp,t,rho)
        cpiint=CPI(t,x)
        cptint=CPT(t,x)
        dPdrho=DPDBWR(icomp,t,rho)
        dPT=DPTBWR(icomp,t,rho)
        e=a-t*dadt
c    &    +cpiint-R*(t-tref(icomp))  ! R*tref is const, merge w/ href
     &    +cpiint-RT
     &    -href(icomp)
c
        s=-dadt+R*log(rhoref(icomp)/rhos)+cptint-R*log(t/tref(icomp))
     &    -sref(icomp)
        cv=-t*D2ABWR(icomp,t,rho)+CP0(t,x)-R
        if (rho.gt.1.0d-8) then
          cp=cv+t/rho**2*dPT**2/dPdrho
        else
          cp=CP0(t,x)
        end if
        if (rho.gt.1.0d-8) then
          h=e+p/rho
          hjt=(t/rho*dPT/dPdrho-1.0d0)/rho/cp
        else
          h=e+RT
          call VIRB (t,x,b)
          call DBDT (t,x,dbt)
          hjt=(dbt*t-b)/cp
        end if
c  if any of the factors in speed of sound are negative (e.g. in two-
c  phase region) return 0.0
c       w=SQRT(1.0d3/wm*cp/cv*dPdrho)
        w2=cp/cv*dPdrho
        if (w2.gt.0.0d0) then
          w=SQRT(1.0d3/wm(icomp)*w2)
        else
          w=0.0d0
        end if
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)  !real-gas terms
          phi01=PHIK(icomp,0,1,tau,del)
          phi10=PHIK(icomp,1,0,tau,del)
          phi11=PHIK(icomp,1,1,tau,del)
          phi02=PHIK(icomp,0,2,tau,del)
          phi20=PHIK(icomp,2,0,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi00=PHIX(0,0,tau,del,x)   !real-gas terms
          phi01=PHIX(0,1,tau,del,x)
          phi10=PHIX(1,0,tau,del,x)
          phi11=PHIX(1,1,tau,del,x)
          phi02=PHIX(0,2,tau,del,x)
          phi20=PHIX(2,0,tau,del,x)
        end if
c       write (*,1003) t,rho,phi00,phi01,phi10,phi11,phi02,phi20
c1003   format (1x,' THERM--t,rho,PHIs:   ',f8.2,f12.6,6d16.6)
c
        phig00=PHI0(0,0,t,rhos,x)      !ideal-gas terms
        phig10=PHI0(1,0,t,rho,x)
        phig20=PHI0(2,0,t,rho,x)
c       write (*,1005) (x(j),j=1,ncmax)
c1005   format (1x,' THERM--output x(i): ',5f14.8)
c       write (*,1024) phig00,phig10,phig20
c1024   format (1x,' THERM--phig-00/01/02:',20x,3d16.6)
        p=RT*rho*(1.0d0+phi01)
        e=RT*(phig10+phi10)
        if (nc.eq.1 .or. ic.ne.0) then
          icomp=1
          if (ic.ne.0) icomp=ic
          e=e-href(icomp)
        else
          do i=1,nc
            e=e-x(i)*href(i)
          enddo
        endif
        h=e+RT*(1.0d0+phi01)
        s=R*(phig10+phi10-phig00-phi00)
c       write (*,*) ' THERM--t,rho,x,s,sref:  ',t,rho,x(1),s,sref(1)
        if (nc.eq.1 .or. ic.ne.0) then
          icomp=1
          if (ic.ne.0) icomp=ic
          s=s-sref(icomp)
        else
          do i=1,nc
            s=s-x(i)*sref(i)
          enddo
        endif
        cv=-R*(phi20+phig20)
c       write (*,*) ' THERM--tau,del,Cv_resid:  ',tau,del,phi20
        cp=cv+R*(1.0d0+phi01-phi11)**2/(1.0d0+2.0d0*phi01+phi02)
        if (cv.ne.0.d0) w2=RT*cp/cv*(1.0d0+2.0d0*phi01+phi02)
c  if any of the factors in speed of sound are negative (e.g. in two-
c  phase region) return 0.0
        if (w2.gt.0.0d0) then
          w=SQRT(w2*1.0d3/WMOL(x))  !convert from molar to mass units
        else
          w=0.0d0
        end if
        if (rho.gt.1.0d-8) then
          hjt=-1.0d0/(cp*rho)*(phi01+phi02+phi11)/
     &        (1.0d0+2.0d0*phi01+phi02)
        else
          call VIRB (t,x,b)
          call DBDT (t,x,dbt)
          hjt=(dbt*t-b)/cp
        end if
      end if
c
      RETURN
      end                                              !subroutine THERM
c
c ======================================================================
c
      subroutine THERM2 (t,rho,x,p,e,h,s,cv,cp,w,Z,hjt,A,G,
     &                   xkappa,beta,dPdrho,d2PdD2,dPT,drhodT,drhodP,
     &                   spare1,spare2,spare3,spare4)
c
c  compute thermal quantities as a function of temperature, density,
c  and compositions using core functions (Helmholtz free energy, ideal
c  gas heat capacity and various derivatives and integrals)
c
c  this routine is the same as THERM, except that additional properties
c  are calculated
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c        p--pressure [kPa]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [J/mol-K]
c       Cv--isochoric heat capacity [J/mol-K]
c       Cp--isobaric heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c        Z--compressibility factor (= PV/RT) [dimensionless]
c      hjt--isenthalpic Joule-Thompson coefficient [K/kPa]
c        A--Helmholtz energy [J/mol]
c        G--Gibbs free energy [J/mol]
c   xkappa--isothermal compressibility (= -1/V dV/dP = 1/rho dD/dP) [1/kPa]
c     beta--volume expansivity (= 1/V dV/dT = -1/rho dD/dT) [1/K]
c   dPdrho--derivative dP/drho [kPa-L/mol]
c   d2PdD2--derivative d^2P/drho^2 [kPa-L^2/mol^2]
c      dPT--derivative dP/dT [kPa/K]
c   drhodT--derivative drho/dT [mol/(L-K)]
c   drhodP--derivative drho/dP [mol/(L-kPa)]
c   sparei--4 space holders for possible future properties
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  03-16-98  MM, original version; based on THERM
c  03-30-98 EWL, add Joule-Thompson coeff to MBWR case
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  12-02-98 EWL, remove compositional dependence for pure fluids
c  12-02-98 EWL, restructure to closely mimic THERM
c  12-02-98 EWL, add the reference states to the calculation of A and G
c  12-22-98 EWL, recalculate R for mixtures based on values for pure fluids
c  02-11-99 EWL, skip calculation of d2PdD2 if rho=0
c  04-27-01 EWL, change order of calculation of a and g in BWR section
c  04-27-01 DGF, change sign before sref in the calculation of a and g
c  09-05-02 EWL, add ideal gas isothermal compressibility and d2pdD2
c  09-05-02 EWL, change check on rho from 1.d-10 to 1.0d-8 to avoid
c                zero's coming back from core_feq at low rhos.
c  06-14-06 EWL, change dPdD to dPdrho, etc, to avoid compiler problems with subroutine dPdD, etc.
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: THERM2
c     dll_export THERM2
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c  common block containing flags to GUI
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      dimension x(ncmax)
c

      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
      spare1=xnotc   !flag indicating not calculated
      spare2=xnotc
      spare3=xnotc
      spare4=xnotc
c
      call RMIX(x)
      RT=R*t
      call DBDT (t,x,dbt)
      call VIRB (t,x,b)
      rhos=rho
c  entropy calc will crash if rho = 0
      if (rho.lt.1.0d-20) rhos=1.0d-20
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        ar=ABWR(icomp,t,rho)
        dadt=DABWR(icomp,t,rho)
        cpiint=CPI(t,x)
        cptint=CPT(t,x)
        dPdrho=DPDBWR(icomp,t,rho)
        dPT=DPTBWR(icomp,t,rho)
        e=ar-t*dadt
c    &    +cpiint-R*(t-tref(icomp))  ! R*tref is const, merge w/ href
     &    +cpiint-RT
     &    -href(icomp)
c  additional properties added to THERM2 not in THERM
        d2PdD2=D2PBWR(icomp,t,rho)
        drhodT=-dPT/dPdrho
        drhodP=1.0d0/dPdrho
c
        s=-dadt+R*log(rhoref(icomp)/rhos)+cptint-R*log(t/tref(icomp))
     &    -sref(icomp)
        cv=-t*D2ABWR(icomp,t,rho)+CP0(t,x)-R
        A=e-t*s
        if (rho.gt.1.0d-8) then
          G=A+p/rho
          beta=-drhodT/rho
          xkappa=drhodP/rho
        else
          G=A+R*t
          beta=1.0d0/t  !if rho = 0, then ideal-gas behavior
          xkappa=xnotc
          if (p.gt.0.d0) xkappa=1.d0/p
        end if
        if (rho.gt.1.0d-8) then
          cp=cv+t/rho**2*dPT**2/dPdrho
        else
          cp=CP0(t,x)
        end if
        if (rho.gt.1.0d-8) then
          h=e+p/rho
          hjt=(t/rho*DPTBWR(icomp,t,rho)/dPdrho-1.0d0)/rho/cp
        else
          h=e+RT
          hjt=(dbt*t-b)/cp
        end if
c  if any of the factors in speed of sound are negative (e.g. in two-
c  phase region) return 0.0
c       w=SQRT(1.0d3/wm*cp/cv*dPdrho)
        w2=cp/cv*dPdrho
        if (w2.gt.0.0d0) then
          w=SQRT(1.0d3/wm(icomp)*w2)
        else
          w=0.0d0
        end if
c
      else
c  call general PHIK or PHIX routines for all other models
        if (icomp.ne.0) then
c  pure fluid
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)  !real-gas terms
          phi01=PHIK(icomp,0,1,tau,del)
          phi10=PHIK(icomp,1,0,tau,del)
          phi11=PHIK(icomp,1,1,tau,del)
          phi02=PHIK(icomp,0,2,tau,del)
          phi20=PHIK(icomp,2,0,tau,del)
          phi03=PHIK(icomp,0,3,tau,del)

        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi00=PHIX(0,0,tau,del,x)   !real-gas terms
          phi01=PHIX(0,1,tau,del,x)
          phi10=PHIX(1,0,tau,del,x)
          phi11=PHIX(1,1,tau,del,x)
          phi02=PHIX(0,2,tau,del,x)
          phi20=PHIX(2,0,tau,del,x)
          phi03=PHIX(0,3,tau,del,x)
        end if
c       write (*,1003) t,rho,phi00,phi01,phi10,phi11,phi02,phi20
c1003   format (1x,' THERM--t,rho,PHIs:   ',f8.2,f12.6,6d16.6)
c
        phig00=PHI0(0,0,t,rhos,x)      !ideal-gas terms
        phig10=PHI0(1,0,t,rho,x)
        phig20=PHI0(2,0,t,rho,x)
c       write (*,1005) (x(j),j=1,ncmax)
c1005   format (1x,' THERM--output x(i): ',5f14.8)
c       write (*,1024) phig00,phig10,phig20
c1024   format (1x,' THERM--phig-00/01/02:',20x,3d16.6)
        p=RT*rho*(1.0d0+phi01)
        e=RT*(phig10+phi10)
        if (icomp.ne.0) then
          e=e-href(icomp)
        else
          do i=1,nc
            e=e-x(i)*href(i)
          enddo
        endif
        h=e+RT*(1.0d0+phi01)
        s=R*(phig10+phi10-phig00-phi00)
c       write (*,*) ' THERM--t,rho,x,s,sref:  ',t,rho,x(1),s,sref(1)
        if (icomp.ne.0) then
          s=s-sref(icomp)
        else
          do i=1,nc
            s=s-x(i)*sref(i)
          enddo
        endif
        cv=-R*(phi20+phig20)
c       write (*,*) ' THERM--tau,del,Cv_resid:  ',tau,del,phi20
        cp=cv+R*(1.0d0+phi01-phi11)**2/
     &     (1.0d0+2.0d0*phi01+phi02)
        w2=RT*cp/cv*(1.0d0+2.0d0*phi01+phi02)
c  if any of the factors in speed of sound are negative (e.g. in two-
c  phase region) return 0.0
        if (w2.gt.0.0d0) then
          w=SQRT(w2*1.0d3/WMOL(x))  !convert from molar to mass units
        else
          w=0.0d0
        end if
        if (rho.gt.1.0d-8) then
          hjt=-1.0d0/(cp*rho)*(phi01+phi02+phi11)/
     &        (1.0d0+2.0d0*phi01+phi02)
        else
          hjt=(dbt*t-b)/cp
        end if
c  additional properties added to THERM2 not in THERM
        A=RT*(phi00+phig00)
        G=A+RT*(1.0d0+phi01)
        if (icomp.ne.0) then
          a=a-href(icomp)+sref(icomp)*t
          g=g-href(icomp)+sref(icomp)*t
        else
          do i=1,nc
            a=a-x(i)*(href(i)-sref(i)*t)
            g=g-x(i)*(href(i)-sref(i)*t)
          enddo
        endif
        dPdrho=RT*(1.0d0+2.0d0*phi01+phi02)
        dPT=R*rho*(1.0d0+phi01-phi11)
        drhodP=1.0d0/(RT*(1.0d0+2.0d0*phi01+phi02))
        drhodT=-rho*(1.0d0+phi01-phi11)/(t*(1.0d0+2.0d0*phi01+phi02))
        if (rho.gt.1.0d-8) then
          d2PdD2=RT/rho*(2.0d0*phi01+4.0d0*phi02+phi03)
          beta=-drhodT/rho
          xkappa=drhodP/rho
        else
          d2PdD2=2.d0*b*R*t
          beta=1.0d0/t  !if rho = 0, then ideal-gas behavior
          xkappa=xnotc
          if (p.gt.0.d0) xkappa=1.d0/p
        end if
      end if
      if (rho.lt.1.0d-20) then
        Z=1.0d0       !if rho = 0, then ideal-gas behavior
      else
        Z=p/(RT*rho)
      end if
c
      RETURN
      end                                             !subroutine THERM2
c
c ======================================================================
c
      subroutine THERM0 (t,rho,x,p0,e0,h0,s0,cv0,cp00,w0,A0,G0)
c
c  compute ideal gas thermal quantities as a function of temperature, density,
c  and compositions using core functions
c
c  this routine is the same as THERM, except it only calculates ideal gas
c  properties (Z=1) at any temperature and density
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c       p0--pressure [kPa]
c       e0--internal energy [J/mol]
c       h0--enthalpy [J/mol]
c       s0--entropy [J/mol-K]
c      Cv0--isochoric heat capacity [J/mol-K]
c     Cp00--isobaric heat capacity [J/mol-K]
c       w0--speed of sound [m/s]
c       A0--Helmholtz energy [J/mol]
c       G0--Gibbs free energy [J/mol]
c
c  11-26-02 EWL, original version; based on THERM
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: THERM0
c     dll_export THERM0
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      RT=R*t
      rhos=rho
c  entropy calc will crash if rho = 0
      if (rho.lt.1.0d-20) rhos=1.0d-20
      phig00=PHI0(0,0,t,rhos,x)      !ideal-gas terms
      phig10=PHI0(1,0,t,rho,x)
      phig20=PHI0(2,0,t,rho,x)
      p0=RT*rho
      e0=RT*phig10
      s0=R*(phig10-phig00)
      A0=RT*phig00
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
        e0=e0-href(icomp)
        s0=s0-sref(icomp)
        a0=a0-href(icomp)+sref(icomp)*t
      else
        do i=1,nc
          e0=e0-x(i)*href(i)
          s0=s0-x(i)*sref(i)
          a0=a0-x(i)*(href(i)-sref(i)*t)
        enddo
      endif
      cv0=-R*phig20
      cp00=cv0+R
      h0=e0+RT
      G0=A0+RT
      w2=RT*cp00/cv0
      if (w2.gt.0.0d0) then
        w0=SQRT(w2*1.0d3/WMOL(x))
      else
        w0=0.0d0
      end if
c
      RETURN
      end                                             !subroutine THERM0
c
c ======================================================================
c
      subroutine ENTRO (t,rho,x,s)
c
c  compute entropy as a function of temperature, density and composition
c  using core functions (temperature derivative of Helmholtz free energy
c  and ideal gas integrals)
c
c  based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994,
c  equations A5, A19 - A26
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c        s--entropy [J/mol-K]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-05-94  MM, original version
c  10-03-95  MM, change /MODEL/--models specified by strings
c  10-10-95  MM, compute ideal gas pressure and pass to PHI0
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-08-95  MM, convert calls to PHI0, CP0, CPI, CPT to mixture form
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  01-19-96  MM, fix ref state for HMX: s = s - sum[x(i)*sref(i)]
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  04-19-96  MM, change call to PHI0:  pass rho instead of pideal
c                calculate s using PHI0 rather than Cp0
c  07-19-96  MM, change general calls to PHI0 (same as THERM)
c  04-22-97  MM, lower bound on rho for s calc set to 1.0d-20
c  10-01-97  MM, add compiler switch to allow access by DLL
c  12-02-98 EWL, remove compositional dependence for pure fluids
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: ENTRO
c     dll_export ENTRO
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (rho.lt.1.0d-20) then
c  entropy calc will crash if rho = 0
        rhos=1.0d-20
      else
        rhos=rho
      end if
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        s=-DABWR(icomp,t,rho)+R*log(rhoref(icomp)/rhos)+CPT(t,x)
     &    -R*log(t/tref(icomp))-sref(icomp)
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)  !real-gas terms
          phi10=PHIK(icomp,1,0,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi00=PHIX(0,0,tau,del,x)  !real-gas terms
          phi10=PHIX(1,0,tau,del,x)
        end if
c
        phig00=PHI0(0,0,t,rhos,x)      !ideal-gas terms
        phig10=PHI0(1,0,t,rho,x)
        s=R*(phig10+phi10-phig00-phi00)
        if (nc.eq.1 .or. ic.ne.0) then
          icomp=1
          if (ic.ne.0) icomp=ic
          s=s-sref(icomp)
        else
          do i=1,nc
            s=s-x(i)*sref(i)
          enddo
        endif
      end if
c
      RETURN
      end                                              !subroutine ENTRO
c
c ======================================================================
c
      subroutine ENTHAL (t,rho,x,h)
c
c  compute enthalpy as a function of temperature, density, and
c  composition using core functions (Helmholtz free energy and ideal
c  gas integrals)
c
c  based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994,
c  equations A7, A18, A19
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c        h--enthalpy [J/mol]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-06-94  MM, original version
c  10-03-95  MM, change /MODEL/--models specified by strings
c  10-10-95  MM, compute ideal gas pressure and pass to PHI0
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-08-95  MM, convert calls to PHI0, CP0, CPI, CPT to mixture form
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  01-19-96  MM, fix ref state for HMX: h = h - sum[x(i)*href(i)]
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  07-19-96  MM, change general calls to PHI0 (same as THERM)
c  10-01-97  MM, add compiler switch to allow access by DLL
c  12-02-98 EWL, remove compositional dependence for pure fluids
c  04-05-00 EWL, check for rho=0 and avoid division by zero
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: ENTHAL
c     dll_export ENTHAL
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        if (rho.lt.1.d-8) then
          h=ABWR(icomp,t,rho)-t*DABWR(icomp,t,rho)+CPI(t,x)-href(icomp)
        else
          h=ABWR(icomp,t,rho)-t*DABWR(icomp,t,rho)+
     &     PBWR(icomp,t,rho)/rho-R*t+CPI(t,x)-href(icomp)
        endif
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)
          phi10=PHIK(icomp,1,0,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi10=PHIX(1,0,tau,del,x)
        end if
c
        RT=R*t
        phig10=PHI0(1,0,t,rho,x)
        e=RT*(phig10+phi10)
        if (nc.eq.1 .or. ic.ne.0) then
          icomp=1
          if (ic.ne.0) icomp=ic
          e=e-href(icomp)
        else
          do i=1,nc
            e=e-x(i)*href(i)
          enddo
        endif
        h=e+RT*(1.0d0+phi01)
      end if
c
      RETURN
      end                                             !subroutine ENTHAL
c
c ======================================================================
c
      subroutine ENERGY (t,rho,x,e)
c
c  compute energy as a function of temperature, density, and
c  composition using core functions (Helmholtz free energy and ideal
c  gas integrals)
c
c  based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c        e--energy [J/mol]
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  12-13-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: ENERGY
c     dll_export ENERGY
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        e=ABWR(icomp,t,rho)-t*DABWR(icomp,t,rho)-R*t+CPI(t,x)
     &   -href(icomp)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi10=PHIK(icomp,1,0,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi10=PHIX(1,0,tau,del,x)
        end if
c
        phig10=PHI0(1,0,t,rho,x)
        e=R*t*(phig10+phi10)
        if (nc.eq.1 .or. ic.ne.0) then
          icomp=1
          if (ic.ne.0) icomp=ic
          e=e-href(icomp)
        else
          do i=1,nc
            e=e-x(i)*href(i)
          enddo
        endif
      end if
c
      RETURN
      end                                             !subroutine ENERGY
c
c ======================================================================
c
      subroutine CVCP (t,rho,x,cv,cp)
c
c  compute isochoric (constant volume) and isobaric (constant pressure)
c  heat capacity as functions of temperature, density, and composition
c  using core functions
c
c  based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994,
c  equation A15, A16
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c       cv--isochoric heat capacity [J/mol-K]
c       cp--isobaric heat capacity [J/mol-K]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-06-94  MM, original version
c  10-03-95  MM, change /MODEL/--models specified by strings
c  10-10-95  MM, compute ideal gas pressure and pass to PHI0
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-08-95  MM, convert calls to PHI0, CP0, CPI, CPT to mixture form
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  07-19-96  MM, change general calls to PHI0 (same as THERM)
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: CVCP
c     dll_export CVCP
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        cv=-t*D2ABWR(icomp,t,rho)+CP0(t,x)-R
        if (rho.gt.1.0d-8) then
          cp=cv+t/rho**2*DPTBWR(icomp,t,rho)**2/DPDBWR(icomp,t,rho)
        else
          cp=cv+R
        end if
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi02=PHIK(icomp,0,2,tau,del)
          phi11=PHIK(icomp,1,1,tau,del)
          phi20=PHIK(icomp,2,0,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi02=PHIX(0,2,tau,del,x)
          phi11=PHIX(1,1,tau,del,x)
          phi20=PHIX(2,0,tau,del,x)
        end if
c
        phig20=PHI0(2,0,t,rho,x)         !ideal-gas term
        cv=-R*(phi20+phig20)
        cp=cv+R*(1.0d0+phi01-phi11)**2/
     &     (1.0d0+2.0d0*phi01+phi02)
      end if
c
      RETURN
      end                                               !subroutine CVCP
c
c ======================================================================
c
      subroutine CVCPK (icomp,t,rho,cv,cp)
c
c  compute isochoric (constant volume) and isobaric (constant pressure)
c  heat capacity as functions of temperature for a given component
c
c  analogous to CVCP, except for component icomp, this is used by transport
c  routines to calculate Cv & Cp for the reference fluid (component zero)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  outputs:
c       cv--isochoric heat capacity [J/mol-K]
c       cp--isobaric heat capacity [J/mol-K]
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  06-16-97  MM, original version; based on CVCP
c  10-01-97  MM, add compiler switch to allow access by DLL
c  03-06-98  MM, check hmxeos, not heos, for 'BWR' (crash for icomp = 0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: CVCPK
c     dll_export CVCPK
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      if (hmxeos(icomp).eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        if (rho.gt.1.0d-8) then
          cv=-t*D2ABWR(icomp,t,rho)+CP0K(icomp,t)-R
          cp=cv+t/rho**2*DPTBWR(icomp,t,rho)**2/DPDBWR(icomp,t,rho)
        else
          cp=CP0K(icomp,t)
          cv=cp-R
        end if
c
      else
c  call general PHIK or PHIX routines for all other models
        tau=tz(icomp)/t
        del=rho/rhoz(icomp)
        phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
        phi02=PHIK(icomp,0,2,tau,del)
        phi11=PHIK(icomp,1,1,tau,del)
        phi20=PHIK(icomp,2,0,tau,del)
c
        phig20=PHI0K(icomp,2,0,t,rho)  !ideal-gas term
        cv=-R*(phi20+phig20)
        cp=cv+R*(1.0d0+phi01-phi11)**2/(1.0d0+2.0d0*phi01+phi02)
      end if
c
      RETURN
      end                                              !subroutine CVCPK
c
c ======================================================================
c
      subroutine GIBBS (t,rho,x,Ar,Gr)
c
c  compute residual Helmholtz and Gibbs free energy as a function of
c  temperature, density, and composition using core functions
c
c  N.B.  The quantity calculated is
c
c        G(T,rho) - G0(T,P*) = G(T,rho) - G0(T,rho) + RTln(RTrho/P*)
c
c        where G0 is the ideal gas state and P* is a reference pressure
c        which is equal to the current pressure of interest.  Since Gr
c        is used only as a difference in phase equilibria calculations
c        where the temperature and pressure of the phases are equal, the
c        (RT/P*) part of the log term will cancel and is omitted.
c
c        "normal" (not residual) A and G are computed by subroutine AG
c
c  based on derivations in Younglove & McLinden, JPCRD 23 #5, 1994,
c  equations A8 - A12
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c       Ar--residual Helmholtz free energy [J/mol]
c       Gr--residual Gibbs free energy [J/mol]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-07-94  MM, original version
c  08-07-95  MM, add calls to Fundamental (Helmholtz) EOS
c  10-03-95  MM, change /MODEL/--models specified by strings
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: GIBBS
c     dll_export GIBBS
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        Ar=ABWR(icomp,t,rho)
        Gr=Ar+PBWR(icomp,t,rho)/rho+R*t*(-1.0d0+log(rho))
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)  !real-gas terms
          phi01=PHIK(icomp,0,1,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi00=PHIX(0,0,tau,del,x)  !real-gas terms
          phi01=PHIX(0,1,tau,del,x)
        end if
c
        RT=R*t
        Ar=RT*phi00
        Gr=Ar+RT*(1.0d0+phi01)+RT*(-1.0d0+log(rho))
      end if
c
      RETURN
      end                                              !subroutine GIBBS
c
c ======================================================================
c
      subroutine AG (t,rho,x,a,g)
c
c  compute Helmholtz and Gibbs energies as a function of temperature,
c  density, and composition.
c
c  N.B.  These are not residual values (those are calculated by GIBBS).
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c        a--Helmholtz energy [J/mol]
c        g--Gibbs free energy [J/mol]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c   3-27-98 EWL, original version
c  12-02-98 EWL, reorganize code so to eliminate x(i) in pure fluid calculation
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: AG
c     dll_export AG
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (rho.lt.1.0d-20) then
c  calc will crash if rho = 0
        rhos=1.0d-20
      else
        rhos=rho
      end if
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        a=ABWR(icomp,t,rho)
        cpiint=CPI(t,x)
        cptint=CPT(t,x)
        a=a+cpiint-R*t-href(icomp)
     &   -t*(R*log(rhoref(icomp)/rhos)+cptint-R*log(t/tref(icomp))
     &    -sref(icomp))
        if (rho.ne.0.d0) then
          g=a+p/rho
        else
          g=a+R*t
        endif
c
      else
        RT=R*t
        phig00=PHI0(0,0,t,rhos,x)      !ideal-gas terms
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)  !real-gas terms
          phi01=PHIK(icomp,0,1,tau,del)
          a=RT*(phig00+phi00)-(href(icomp)-t*sref(icomp))
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi00=PHIX(0,0,tau,del,x)   !real-gas terms
          phi01=PHIX(0,1,tau,del,x)
          a=RT*(phig00+phi00)
          do i=1,nc
            a=a-x(i)*(href(i)-t*sref(i))
          enddo
        end if
        g=a+RT*(1.0d0+phi01)
      end if
c
      RETURN
      end                                                 !subroutine AG
c
c ======================================================================
c
      subroutine PRESS (t,rho,x,p)
c
c  compute pressure as a function of temperature,
c  density, and composition using core functions
c
c  direct implementation of core function of corresponding model
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c        p--pressure [kPa]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-19-94  MM, original version
c  08-07-95  MM, add calls to Fundamental (Helmholtz) EOS
c  10-03-95  MM, change /MODEL/--models specified by strings
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PRESS
c     dll_export PRESS
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          p=R*t*rho*(1.0d0+PHIK(icomp,0,1,tau,del))
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          p=R*t*rho*(1.0d0+PHIX(0,1,tau,del,x))
        end if
      end if
c
      RETURN
      end                                              !subroutine PRESS
c
c ======================================================================
c
      subroutine DPDD (t,rho,x,dpdrho)
c
c  compute partial derivative of pressure w.r.t. density at constant
c  temperature as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c   dpdrho--dP/drho [kPa-L/mol]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  04-23-95  MM, original version
c  08-07-95  MM, add calls to Fundamental (Helmholtz) EOS
c  10-03-95  MM, change /MODEL/--models specified by strings
c  11-03-95  MM, add calls to mixture Helmholtz (HMX) model
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  10-16-96  MM, change name from DPRHO to DPDD
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DPDD
c     dll_export DPDD
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        dpdrho=DPDBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi02=PHIK(icomp,0,2,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi02=PHIX(0,2,tau,del,x)
        end if
        dpdrho=R*t*(1.0d0+2.0d0*phi01+phi02)
      end if
c
      RETURN
      end                                               !subroutine DPDD
c
c ======================================================================
c
      subroutine DPDDK (icomp,t,rho,dPdrho)
c
c  compute partial derivative of pressure w.r.t. density at constant
c  temperature as a function of temperature and density for a specified
c  component
c
c  analogous to DPDD, except for component icomp, this is used by transport
c  routines to calculate dP/dD for the reference fluid (component zero)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output:
c   dPdrho--dP/drho [kPa-L/mol]
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  06-16-97  MM, original version; based on DPDD
c  09-29-97  MM, if component uses MBWR, call DPDBWR
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DPDDK
c     dll_export DPDDK
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      if (hmxeos(icomp).eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        dpdrho=DPDBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        tau=tz(icomp)/t
        del=rho/rhoz(icomp)
        phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
        phi02=PHIK(icomp,0,2,tau,del)
        dpdrho=R*t*(1.0d0+2.0d0*phi01+phi02)
      end if
c
      RETURN
      end                                              !subroutine DPDDK
c
c ======================================================================
c
      subroutine DPDD2 (t,rho,x,d2PdD2)
c
c  compute second partial derivative of pressure w.r.t. density at const
c  temperature as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c   d2pdD2--d^2P/drho^2 [kPa-L^2/mol^2]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  06-03-97 EWL, original version
c  10-01-97  MM, add compiler switch to allow access by DLL
c  02-11-99 EWL, skip calculation of d2PdD2 if rho=0
c  09-05-02 EWL, add ideal gas d2PdD2
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DPDD2
c     dll_export DPDD2
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        d2PdD2=D2PBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi02=PHIK(icomp,0,2,tau,del)
          phi03=PHIK(icomp,0,3,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi02=PHIX(0,2,tau,del,x)
          phi03=PHIX(0,3,tau,del,x)
        end if
        if (rho.gt.1.0d-8) then
          d2PdD2=R*t/rho*(2.0d0*phi01+4.0d0*phi02+phi03)
        else
          call VIRB (t,x,b)
          d2PdD2=2.d0*b*R*t
        end if
      end if
c
      RETURN
      end                                              !subroutine DPDD2
c
c ======================================================================
c
      subroutine DPDT (t,rho,x,dpt)
c
c  compute partial derivative of pressure w.r.t. temperature at constant
c  density as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c      dpt--dP/dT [kPa/K]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  10-16-96  MM, original version, based on DPDD
c  10-28-96  MM, insert missing rho into form using PHI's
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DPDT
c     dll_export DPDT
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        dpt=DPTBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi11=PHIK(icomp,1,1,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi11=PHIX(1,1,tau,del,x)
        end if
        dpt=R*rho*(1.0d0+phi01-phi11)
      end if
c
      RETURN
      end                                               !subroutine DPDT
c
c ======================================================================
c
      subroutine DPDTK (icomp,t,rho,dpt)
c
c  compute partial derivative of pressure w.r.t. temperature at constant
c  density as a function of temperature and density for a specified component
c
c  analogous to DPDT, except for component icomp, this is used by transport
c  routines to calculate dP/dT
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output:
c      dpt--dP/dT [kPa/K]
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  07-07-98 EWL, original version, based on DPDT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DPDTK
c     dll_export DPDTK
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      if (hmxeos(icomp).eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        dpt=DPTBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        tau=tz(icomp)/t
        del=rho/rhoz(icomp)
        phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
        phi11=PHIK(icomp,1,1,tau,del)
        dpt=R*rho*(1.0d0+phi01-phi11)
      end if
c
      RETURN
      end                                              !subroutine DPDTK
c
c ======================================================================
c
      subroutine DDDP (t,rho,x,drhodp)
c
c  compute partial derivative of density w.r.t. pressure at constant
c  temperature as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c   drhodp--drho/dP [mol/(L-kPa)]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  08-29-97  MM, original version, based on DPDD (just the inverse)
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DDDP
c     dll_export DDDP
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        drhodp=1.0d0/DPDBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi02=PHIK(icomp,0,2,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi02=PHIX(0,2,tau,del,x)
        end if
        drhodp=1.0d0/(R*t*(1.0d0+2.0d0*phi01+phi02))
      end if
c
      RETURN
      end                                               !subroutine DDDP
c
c ======================================================================
c
      subroutine DDDT (t,rho,x,drhodt)
c
c  compute partial derivative of density w.r.t. temperature at constant
c  pressure as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c   drhodt--drho/dT [mol/(L-K)]
c
c   d(rho)/d(T) = -d(rho)/dP x dP/dT = -dP/dT / (dP/d(rho))
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  08-29-97  MM, original version, based on DPDD and DPDT
c  10-01-97  MM, add compiler switch to allow access by DLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DDDT
c     dll_export DDDT
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        drhodt=-DPTBWR(icomp,t,rho)/DPDBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi11=PHIK(icomp,1,1,tau,del)
          phi02=PHIK(icomp,0,2,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi11=PHIX(1,1,tau,del,x)
          phi02=PHIX(0,2,tau,del,x)
        end if
        drhodt=-rho*(1.0d0+phi01-phi11)/(t*(1.0d0+2.0d0*phi01+phi02))
      end if
c
      RETURN
      end                                               !subroutine DDDT
c
c ======================================================================
c
      subroutine DHD1(t,rho,x,dhdt_d,dhdt_p,dhdd_t,dhdd_p,dhdp_t,dhdp_d)
c
c  compute partial derivatives of enthalpy w.r.t. t, p, or rho at constant
c  t, p, or rho as a function of temperature, density, and composition
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c      dht--dH/dT [J/mol-K]
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  10-02-00 EWL, original version
c  05-30-06 EWL, change subroutine name from DHDT to DHD1, and add other
c                derivates of h with respect to t, p, or rho
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DHD1
c     dll_export DHD1
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      dimension x(ncmax)
c
      call RMIX(x)
      rhos=rho
      if (rho.lt.1.0d-10) rhos=1.0d-10
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        tau=tz(icomp)/t
        del=rho/rhoz(icomp)
        phi01=PHIBWR(icomp,0,1,tau,del)  !real-gas terms
        phi11=PHIBWR(icomp,1,1,tau,del)
        phi20=PHIBWR(icomp,2,0,tau,del)
        phi02=PHIBWR(icomp,0,2,tau,del)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
          phi11=PHIK(icomp,1,1,tau,del)
          phi20=PHIK(icomp,2,0,tau,del)
          phi02=PHIK(icomp,0,2,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)  !real-gas terms
          phi11=PHIX(1,1,tau,del,x)
          phi20=PHIX(2,0,tau,del,x)
          phi02=PHIX(0,2,tau,del,x)
        end if
      end if
      phig20=PHI0(2,0,t,rho,x)
      phig11=PHI0(1,1,t,rho,x)
      call THERM2 (t,rhos,x,p,e,h,s,cv,cp,w,Z,hjt,A,G,xkappa,beta,
     &      dPdrho,d2PdD2,dPT,drhodT,drhodP,spare1,spare2,spare3,spare4)
      dhdt_p=cp
      dhdt_d=R*(-phig20-phi20+phi01-phi11+1.d0)
      if (rho.gt.1.0d-8) then
        dhdp_t=1.d0/rho+t*drhodT/rho**2
        dhdd_t=R*T/rho*(phig11+phi11+phi01+phi02)
        dhdp_d=dhdp_t+dhdt_p/dPT
        dhdd_p=dhdd_t+dhdt_d/drhodT
      else
        call VIRB (t,x,b)
        call DBDT (t,x,dbt)
        dhdp_t=1.d0/rhos+t*drhodT/rhos**2
        dhdd_t=-r*t**2*dbt+r*t*b
        dhdp_d=xinf
        dhdd_p=xinf
      endif
c
      RETURN
      end                                               !subroutine DHD1
c
c ======================================================================
c
      subroutine FGCTY (t,rho,x,f)
c
c  compute fugacity for each of the nc components of a mixture by
c  numerical differentiation (using central differences) of the
c  dimensionless residual Helmholtz energy
c
c  based on derivations in E.W. Lemmon, MS Thesis, University of Idaho
c  (1991); section 3.2
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  output:
c        f--array (1..nc) of fugacities [kPa]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  12-15-95  MM, original version
c  12-18-95  MM, add pure component fugacity as a special case
c  01-08-96  MM, bug on call to PHIFEQ (wrong arguments)
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
c  03-19-19  MM, add dipole moment to /CCON/
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  12-16-97  MM, add check for rho = 0; overflow on exponent (set to xerr)
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  12-22-98 EWL, recalculate R for mixtures based on values for pure fluids
c  05-08-06 EWL, modify how delp and deln are calculated for x>0.9999
c  01-25-07 EWL, change default f(i) from 0 to 1.  Skip calculation if x(i)=0
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: FGCTY
c     dll_export FGCTY
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      dimension x(ncmax),f(ncmax)
      dimension xplus(ncmax),xminus(ncmax)
      character*3 hpheq,heos,hmxeos,hmodcp
      common /NCOMP/ nc,ic
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
c  flags indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,xerr,xnotd,xnotc
      delmol=1.0d-4
      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
c
c  fill output fugacity array with zeros (final value for undefined
c  components and insurance against problems for others)
      do i=1,nc
        f(i)=1.0d0
      enddo
c
      call RMIX(x)
c  check for zero input density
      if (rho.lt.1.0d-20) RETURN
c
      RTrho=R*t*rho
      if (icomp.eq.0) then
        do i=1,nc
          if (ABS(x(i)-1.0d0).lt.1.d-14) icomp=i
        enddo
      endif
      if (icomp.ne.0) then
c  pure component
        if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--use BWR-specific routines
          Ar=ABWR(icomp,t,rho)
          p=PBWR(icomp,t,rho)
          f(icomp)=RTrho*exp(Ar/(R*t)+p/RTrho-1.0d0)
        else
c  for other models, use general PHIK routines
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi00=PHIK(icomp,0,0,tau,del)
          phi01=PHIK(icomp,0,1,tau,del)
c  check for potential under- or over-flow (can happen in 2-phase, but
c  the fugacity is meaningless there anyway)
          arg=phi00+phi01
          if (ABS(arg).lt.500.0d0) then
            f(icomp)=RTrho*exp(arg)
          else
            f(icomp)=xerr
          end if
        end if
      else
c
c  mixture
        deln=delmol       !initialize only
        do i=1,nc
c  compute positive and negative increments to number of moles
          if (x(i).gt.0.d0) then
          if (x(i).lt.delmol) then
c  special case--composition of component i is nearly zero
ccc         deln=-x(i)
ccc         delp=2.0d0*delmol-deln
            deln=-x(i)/2.d0
            delp=-deln
          else if (x(i).gt.1.0d0-deln) then
c  special case--composition of component i is nearly one (pure fluid)
ccc         delp=1.0d0-x(i)
ccc         deln=-2.0d0*delmol+delp
            delp=(1.0d0-x(i))/2.d0
            deln=-delp
          else
c  general case:  deln < x(i) < 1 - deln
            delp=delmol
            deln=-delmol
          end if
c         write (*,*) ' FGCTY--delp,deln: ',delp,deln
          delp1=1.0d0/(1.0d0+delp)
          deln1=1.0d0/(1.0d0+deln)
c  since total number of moles is now 1 + (delp or deln), all of the
c  compositions have changed
          do j=1,nc
            xplus(j)=x(j)*delp1
            xminus(j)=x(j)*deln1
          enddo
          xplus(i)=(x(i)+delp)*delp1
          xminus(i)=(x(i)+deln)*deln1
c  derivative is at constant volume, so must adjust density
          Dplus=rho*(1.0d0+delp)
          Dminus=rho*(1.0d0+deln)
c  compute residual Helmholtz at 'plus' and 'minus' density and composition
c  could call subroutine GIBBS here, but more efficient to directly call
c  the core routines (via PHIX)
          call REDX (xplus,t0,rho0)
          tau=t0/t
          del=Dplus/rho0
          Aplus=PHIX(0,0,tau,del,xplus)          !real-gas terms
          call REDX (xminus,t0,rho0)
          tau=t0/t
          del=Dminus/rho0
          Aminus=PHIX(0,0,tau,del,xminus)        !real-gas terms
          dadn=((1.0d0+delp)*Aplus-(1.0d0+deln)*Aminus)/(delp-deln)
c         write (*,*) ' FGCTY--delp,deln:  ',delp,deln
c         write (*,*) ' FGCTY--i,A+, A-, dAdN: ',i,Aplus,Aminus,dadn
c  check for potential under- or over-flow (can happen in 2-phase, but
c  the fugacity is meaningless there anyway)
          if (ABS(dadn).lt.100.0d0) then
            f(i)=x(i)*RTrho*exp(dadn)             !A is dimensionless
          else
            f(i)=xerr
          end if
          end if
        enddo
      end if
c
      RETURN
      end                                              !subroutine FGCTY
cc
cc ======================================================================
cc
c      subroutine ACTVY (t,rho,x,gamma)
cc
cc  !temp--work in progress:  this routine returns gamma = 1.0
cc
cc  compute activity coefficients for each of the nc components of a
cc  mixture by numerical differentiation
cc
cc  inputs:
cc        t--temperature [K]
cc      rho--molar density [mol/L]
cc        x--composition array [mol frac]
cc  output:
cc    gamma--array (1..nc) of activity coefficients
cc
cc  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
cc  12-18-95  MM, original version
cc  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
cc                add Zcrit to common /CCON/
cc                replace calls to PHIHMX, PHIFEQ with general PHIX, PHIK
cc  03-19-19  MM, add dipole moment to /CCON/
cc  03-22-96  MM, replace /MODEL/ with /EOSMOD/
cc  10-01-97  MM, add compiler switch to allow access by DLL
cc  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
cc
c      implicit double precision (a-h,o-z)
c      implicit integer (i-k,m,n)
c
ccDEC$ ATTRIBUTES DLLEXPORT :: ACTVY
c      dll_export ACTVY
cc
c      parameter (ncmax=20)        !max number of components in mixture
c      parameter (nrefmx=10)       !max number of fluids for transport ECS
c      parameter (n0=-ncmax-nrefmx,nx=ncmax)
c      dimension x(ncmax),gamma(ncmax)
c      dimension xplus(ncmax),xminus(ncmax)
c      character*3 hpheq,heos,hmxeos,hmodcp
c      common /NCOMP/ nc,ic
c      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
c      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
c     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
c     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
c     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
c      data delmol/1.0d-4/
cc
c      call RMIX(x)
cc  fill output array with zeros for undefined components
c      do i=nc,nc
c        gamma(i)=0.0d0
c      enddo
cc
c      if (nc.eq.1 .or. ic.ne.0) then
cc  pure component
c        gamma(1)=1.0d0
c      else
c        deln=delmol       !initialize only
cc  mixture
c        do i=1,nc
cc  compute positive and negative increments to number of moles
c        if (x(i).lt.delmol) then
cc  special case--composition of component i is nearly zero
c          deln=-x(i)
c          delp=2.0d0*delmol-deln
c        else if (x(i).gt.1.0d0-deln) then
cc  special case--composition of component i is nearly one (pure fluid)
c          delp=1.0d0-x(i)
c          deln=-2.0d0*delmol+delp
c        else
cc  general case:  deln < x(i) < 1 - deln
c          delp=delmol
c          deln=-delmol
c        end if
cc       write (*,*) ' ACTVY--delp,deln: ',delp,deln
c        delp1=1.0d0/(1.0d0+delp)
c        deln1=1.0d0/(1.0d0+deln)
cc  since total number of moles is now 1 + (delp or deln), all of the
cc  compositions have changed
c        do j=1,nc
c          xplus(j)=x(j)*delp1
c          xminus(j)=x(j)*deln1
c        enddo
c        xplus(i)=(x(i)+delp)*delp1
c        xminus(i)=(x(i)+deln)*deln1
cc  derivative is at constant volume, so must adjust density
c        Dplus=rho*(1.0d0+delp)
c        Dminus=rho*(1.0d0+deln)
cc  compute at 'plus' and 'minus' density and composition
cc  could call general subroutines here, but more efficient to directly
cc  call the core routines (via PHIX)
c        call REDX (xplus,t0,rho0)
c        tau=t0/t
c        del=Dplus/rho0
c        Aplus=PHIX(0,0,tau,del,xplus)          !real-gas terms
c        call REDX (xminus,t0,rho0)
c        tau=t0/t
c        del=Dminus/rho0
c        Aminus=PHIX(0,0,tau,del,xminus)        !real-gas terms
c        gamma(i)=1.0d0
c        enddo
c      end if
cc
c      RETURN
c      end                                              !subroutine ACTVY
c
c ======================================================================
c
      subroutine VIRB (t,x,b)
c
c  compute second virial coefficient as a function of temperature
c  and composition.
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac]
c  outputs:
c        b--second virial coefficient [L/mol]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  03-27-98 EWL, original version
c  08-30-04 EWL, change rho to 0.00000001
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: VIRB
c     dll_export VIRB
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      rho=0.00000001d0
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        b=(p/rho/R/t-1.0d0)/rho
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,0,1,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(0,1,tau,del,x)
        end if
        b=phi01/rho
      end if
c
      RETURN
      end                                               !subroutine VIRB
c
c ======================================================================
c
      subroutine DBDT (t,x,dbt)
c
c  compute the 2nd derivative of B (B is the second virial coefficient) with
c  respect to T as a function of temperature and composition.
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac]
c  outputs:
c      dbt--2nd derivative of B with respect to T [L/mol-K]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  07-30-01 EWL, original version
c  08-30-04 EWL, change rho to 0.00000001
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: DBDT
c     dll_export DBDT
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      rho=0.00000001d0
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        dpt=DPTBWR(icomp,t,rho)
        dbt=(dpt - p/t)/rho**2/t/R
c
      else
c  call general PHIK or PHIX routines for all other models
       if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi01=PHIK(icomp,1,1,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi01=PHIX(1,1,tau,del,x)
        end if
        dbt=-phi01/rho/t
      end if
c
      RETURN
      end                                               !subroutine DBDT
c
c ======================================================================
c
      subroutine VIRC (t,x,c)
c
c  compute the third virial coefficient as a function of temperature
c  and composition.
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac]
c  outputs:
c        c--third virial coefficient [(L/mol)^2]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c   3-27-98 EWL, original version
c  12-02-98 EWL, change rho to 0.0001 to avoid numerical problems in BWR calc.
c  08-30-04 EWL, change rho to 0.00000001
c  05-24-06 EWL, change rho to 0.000001 for the BWR
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: VIRC
c     dll_export VIRC
c
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
c     rho=0.00000001d0
      rho=0.000001d0
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        p=PBWR(icomp,t,rho)
        dpd=DPDBWR(icomp,t,rho)
        c=((dpd-2.0d0*p/rho)/R/t+1.0d0)/rho**2
c
      else
c  call general PHIK or PHIX routines for all other models
       if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi02=PHIK(icomp,0,2,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi02=PHIX(0,2,tau,del,x)
        end if
        c=phi02/rho**2
      end if
c
      RETURN
      end                                               !subroutine VIRC
c
c ======================================================================
c
      subroutine VIRD (t,x,d)
c
c  compute the fourth virial coefficient as a function of temperature
c  and composition.
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac]
c  outputs:
c        d--fourth virial coefficient [(L/mol)^3]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  11-26-01 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      call RMIX(x)
      rho=0.00000001d0
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
c       icomp=1
c       p=PBWR(icomp,t,rho)
c       dpd=DPDBWR(icomp,t,rho)
c  need to update with correct formula:
c       c=((dpd-2.0d0*p/rho)/R/t+1.0d0)/rho**2
        d=0
c
      else
c  call general PHIK or PHIX routines for all other models
       if (nc.eq.1 .or. ic.ne.0) then
c  pure fluid
          icomp=1
          if (ic.ne.0) icomp=ic
          tau=tz(icomp)/t
          del=rho/rhoz(icomp)
          phi03=PHIK(icomp,0,3,tau,del)
        else
c  mixture
          call REDX (x,t0,rho0)
          tau=t0/t
          del=rho/rho0
          phi03=PHIX(0,3,tau,del,x)
        end if
        d=phi03/rho**3
      end if
c
      RETURN
      end                                               !subroutine VIRD
c
c ======================================================================
c
      subroutine B12 (t,x,b)
c
c  compute b12 as a function of temperature and composition.
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac]
c  outputs:
c        b--b12 [(L/mol)^2]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  04-19-01 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax),xx(ncmax)
c
      if (nc.ne.2) then   !Only calculate b12 for a binary
        b=0
        RETURN
      endif
      call VIRB(t,x,bx)
      xx(1)=1.d0
      xx(2)=0.d0
      call VIRB(t,xx,b1)
      xx(1)=0.d0
      xx(2)=1.d0
      call VIRB(t,xx,b2)
      b=(bx-x(1)**2*b1-x(2)**2*b2)/2.d0/x(1)/x(2)
      RETURN
      end                                                !subroutine B12
c
c ======================================================================
c
      subroutine EXCESS (t,p,x,vE,eE,hE,sE)
c
c  compute excess properties as a function of temperature, pressure,
c  and composition.
c
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        x--composition [array of mol frac]
c  outputs:
c        vE--excess volume [L/mol]
c        eE--excess energy [J/mol]
c        hE--excess enthalpy [J/mol]
c        sE--excess entropy [J/mol-K]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  04-25-02 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax),xx(ncmax),xliq(ncmax),xvap(ncmax)
      character*255 herr
c
      call TPFLSH (t,p,x,d,dl,dv,xliq,xvap,q,eE,hE,sE,cv,cp,w,ierr,herr)
      vE=0.d0
      if (d.ne.0.d0) vE=1.d0/d

      do i=1,nc
        do j=1,nc
          xx(j)=0
        enddo
        xx(i)=1
        call TPFLSH (t,p,xx,d,dl,dv,xliq,xvap,q,e,h,s,cv,cp,w,ierr,herr)
        if (d.ne.0.d0) vE=vE-x(i)/d
        eE=eE-x(i)*e
        hE=hE-x(i)*h
        sE=sE-x(i)*s
      enddo

      RETURN
      end                                             !subroutine EXCESS
c
c ======================================================================
c
      subroutine FPV (t,rho,p,x,f)
c
c  Compute the supercompressibility factor, Fpv.
c
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c  outputs:
c        f--Fpv = sqrt[Z(60 F, 14.73 psia)/Z(T,P)]
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  11-07-02 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: FPV
c     dll_export FPV
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
      character*255 herr
c
      tfpv = 288.705555555556d0    !60 F
      pfpv = 101.55977492837d0     !14.73 psia
      call TPRHO (tfpv,pfpv,x,2,0,dfpv,ierr,herr)
      if (p.gt.0.d0) then
        f=SQRT(pfpv/tfpv/dfpv*rho*t/p)
      else
        f=SQRT(pfpv/tfpv/dfpv/R)
      endif

      RETURN
      end                                                !subroutine FPV
Cc
Cc ======================================================================
Cc
C      subroutine SPECGR (t,rho,p,gr)
Cc
Cc  Compute the specific gravity (relative density).
Cc
Cc  inputs:
Cc        t--temperature [K]
Cc        p--pressure [kPa]
Cc      rho--molar density [mol/L]
Cc  outputs:
Cc       gr--specific gravity [dimensionless]
Cc
Cc  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
Cc  11-07-02 EWL, original version
Cc
C      implicit double precision (a-h,o-z)
C      implicit integer (i-k,m,n)
C      parameter (ncmax=20)        !max number of components in mixture
C      parameter (nrefmx=10)       !max number of fluids for transport ECS
C      parameter (n0=-ncmax-nrefmx,nx=ncmax)
C      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
Cc
C      rhoair=1.d0      !Need to add formulation for air here
C      if (rhoair.gt.0.d0) then
C        gr=rho/rhoair
C      else
C        gr=1.d0
C      endif
C
C      RETURN
C      end                                             !subroutine SPECGR
c
c ======================================================================
c
      subroutine RMIX (x)
c
c  inputs:
c        x--composition [array of mol frac]
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  01-19-01 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /NCOMP/ nc,ic
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax)
c
      if (nc.gt.1 .and. ic.gt.0) then
        R=Reos(ic)
      elseif (nc.gt.1) then
        R=0.0d0
        do i=1,nc
          R=R+x(i)*Reos(i)
        enddo
      endif
      if (R.lt.1.d-10) R=8.314472d0  !Check for bad x(i)
      RETURN
      end                                               !subroutine RMIX
c
c ======================================================================
c
      subroutine THERM3 (t,rho,x,
     &           xkappa,beta,xisenk,xkt,betas,bs,xkkt,thrott,pi,spht)
c
c  Compute miscellaneous thermodynamic properties
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c
c  outputs:
c   xkappa--Isothermal compressibility
c     beta--Volume expansivity
c   xisenk--Isentropic expansion coefficient
c      xkt--Isothermal expansion coefficient
c    betas--Adiabatic compressibility
c       bs--Adiabatic bulk modulus
c     xkkt--Isothermal bulk modulus
c   thrott--Isothermal throttling coefficient
c       pi--Internal pressure
c     spht--Specific heat input
c
c  written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO
c  06-16-06 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: THERM3
c     dll_export THERM3
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      dimension x(ncmax)
c
      call THERM2 (t,rho,x,p,e,h,s,cv,cp,w,Z,hjt,A,G,xkappa,beta,
     &             dPdrho,d2PdD2,dPT,drhodT,drhodP,
     &             spare1,spare2,spare3,spare4)
      wm=WMOL(x)
      if (p.le.0.d0) then
       xkt=1
       xisenk=w**2/R/T*wm*0.001d0
      else
       xkt=rho/p*dPdrho                !Isothermal expansion coefficient
       xisenk=w**2*rho/p*wm*0.001d0    !Isentropic expansion coefficient
      endif
      betas=xnotc
      if (rho.gt.0.d0 .and. w.gt.0.d0)
     &  betas=1.d0/rho/w**2/wm*1000.d0 !Adiabatic compressibility
      bs=xisenk*p                      !Adiabatic bulk modulus
      xkkt=xkt*p                       !Isothermal bulk modulus
      thrott=-hjt*cp            !Isothermal throttling coefficient
      pi=t*dpt-p                       !Internal pressure
      if (dpt.ne.0.d0) then
        spht=rho*cp*dpdrho/dpt         !Specific heat input
      else
        spht=cp*t
      end If
      RETURN
      end                                             !subroutine THERM3
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file prop_sub.f
c ======================================================================
