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 ENTRO (t,rho,x,s)
c     subroutine ENTHAL (t,rho,x,h)
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 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,dp2dD2)
c     subroutine FGCTY (t,rho,x,f)
c     subroutine ACTVY (t,rho,x,gamma)
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/ ncomp
c     common /Gcnst/ R
c     common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
c    &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
c    &              accen(n0:nx),dipole(n0:nx)
c
c  various arrays are dimensioned with parameter statements
c     parameter (ncmax=5)        !max number of components in mixture
c     parameter (n0=-ncmax,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 accomodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: CRITP
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      character*1 htab,hnull
      character*3 hpheq,heos,hmxeos,hmodcp
      character*255 herr
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
      dimension x(ncmax)
c
      ierr=0
      herr=hnull
      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
        ierr=1
        herr=' CRITP--ERROR--specified model not found'//herr
      end if
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 accomodate 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 accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: THERM
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      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/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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)
        e=a-t*dadt
c    &    +cpiint-R*(t-tref(icomp))  ! R*tref is const, merge w/ href
     &    +cpiint-R*t
     &    -href(icomp)
        if (abs(rho).gt.1.0d-10) then
          h=e+p/rho
        else
          h=e+R*t
        end if
        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 (abs(rho).gt.1.0d-10) then
          cp=cv+t/rho**2*DPTBWR(icomp,t,rho)**2/DPDBWR(icomp,t,rho)
        else
          cp=CP0(t,x)
        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*DPDBWR(icomp,t,rho))
        w2=cp/cv*DPDBWR(icomp,t,rho)
        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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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
        RT=R*t
c       pideal=RT*rho                 !compute ideal-gas pressure
        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+del*phi01)
c       e=RT*tau*(phig10+phi10)
c  move tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        e=RT*(phig10+tau*phi10)
        do 240 i=1,nc
        e=e-x(i)*href(i)
  240   continue
        h=e+RT*(1.0d0+del*phi01)
c       s=R*(tau*(phig10+phi10)-phig00-phi00)
c  move tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        s=R*(phig10+tau*phi10-phig00-phi00)
c       write (*,*) ' THERM--t,rho,x,s,sref:  ',t,rho,x(1),s,sref(1)
        do 244 i=1,nc
        s=s-x(i)*sref(i)
  244   continue
c       cv=-R*tau*tau*(phi20+phig20)
c  move tau*tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        cv=-R*(tau*tau*phi20+phig20)
c       write (*,*) ' THERM--tau,del,Cv_resid:  ',tau,del,tau*tau*phi20
        delp01=del*phi01
        cp=cv+R*(1.0d0+delp01-del*tau*phi11)**2/
     &     (1.0d0+2.0d0*delp01+del*del*phi02)
        w2=RT*cp/cv*(1.0d0+2.0d0*delp01+del**2*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
        hjt=-1.0d0/(cp*rho0)*(phi01+del*phi02+tau*phi11)/
     &      (1.0d0+2.0d0*delp01+del*del*phi02)
      end if
c
      RETURN
      end                                              !subroutine THERM
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 accomodate 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 accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: ENTRO
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      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
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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)
c  move tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        s=R*(phig10+tau*phi10-phig00-phi00)
        do 244 i=1,nc
        s=s-x(i)*sref(i)
  244   continue
      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 accomodate 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 accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: ENTHAL
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      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
      dimension x(ncmax)
c
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        ic=1
        h=ABWR(ic,t,rho)-t*DABWR(ic,t,rho)+PBWR(ic,t,rho)/rho-R*t
     &    +CPI(t,x)-href(ic)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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)
c  move tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        e=RT*(phig10+tau*phi10)
        do 240 i=1,nc
        e=e-x(i)*href(i)
  240   continue
        h=e+RT*(1.0d0+del*phi01)
      end if
c
      RETURN
      end                                             !subroutine ENTHAL
c
c ======================================================================
c
      subroutine CVCP (t,rho,x,cv,cp)
c
c  compute isochoric (constant volume) and isochoric (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 accomodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: CVCP
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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 (abs(rho).gt.1.0d-10) 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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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
c  move tau*tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        cv=-R*(tau*tau*phi20+phig20)
        delp01=del*phi01
        cp=cv+R*(1.0d0+delp01-del*tau*phi11)**2/
     &     (1.0d0+2.0d0*delp01+del*del*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 isochoric (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 switches to allow access by DLL
c  03-06-98  MM, check hmxeos, not heos, for 'BWR' (crash for icomp = 0)
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: CVCPK
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
c
      if (hmxeos(icomp).eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        if (abs(rho).gt.1.0d-10) 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
        call REDK (icomp,t0,rho0)
        tau=t0/t
        del=rho/rho0
        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
c  move tau*tau*phi_ideal to the PH0xxx routine (e.g. in core_CPP)
        cv=-R*(tau*tau*phi20+phig20)
        delp01=del*phi01
        cp=cv+R*(1.0d0+delp01-del*tau*phi11)**2/
     &     (1.0d0+2.0d0*delp01+del*del*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  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 accomodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: GIBBS
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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+del*phi01)+RT*(-1.0d0+log(rho))
      end if
c
      RETURN
      end                                              !subroutine GIBBS
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 accomodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: PRESS
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          p=R*t*rho*(1.0d0+del*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+del*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 accomodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DPDD
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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*del*phi01+del*del*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     dPdD--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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DPDDK
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
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
        call REDK (icomp,t0,rho0)
        tau=t0/t
        del=rho/rho0
        phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
        phi02=PHIK(icomp,0,2,tau,del)
        dpdrho=R*t*(1.0d0+2.0d0*del*phi01+del*del*phi02)
      end if
c
      RETURN
      end                                              !subroutine DPDDK
c
c ======================================================================
c
      subroutine DPDD2 (t,rho,x,dp2dD2)
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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DPDD2
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=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
c
      if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--call BWR-specific routines
        icomp=1
        dp2dD2=D2PBWR(icomp,t,rho)
c
      else
c  call general PHIK or PHIX routines for all other models
        if (nc.eq.1) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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
        dp2dD2=R*t/rho*(2.0d0*del*phi01+4.0d0*del*del*phi02
     &        +del**3*phi03)
      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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DPDT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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+del*phi01-del*tau*phi11)
      end if
c
      RETURN
      end                                               !subroutine DPDT
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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DDDP
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /Gcnst/ R
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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*del*phi01+del*del*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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: DDDT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      character*3 hpheq,heos,hmxeos,hmodcp
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      common /NCOMP/ nc
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      dimension x(ncmax)
c
      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) then
c  pure fluid
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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+del*phi01-del*tau*phi11)
     &        /(t*(1.0d0+2.0d0*del*phi01+del*del*phi02))
      end if
c
      RETURN
      end                                               !subroutine DDDT
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 accomodate 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 switches to allow access by DLL
c  12-16-97  MM, add check for rho = 0; overflow on exponent (set to xerr)
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: FGCTY
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      dimension x(ncmax),f(ncmax)
      dimension xplus(ncmax),xminus(ncmax)
      character*3 hpheq,heos,hmxeos,hmodcp
      common /NCOMP/ nc
      common /Gcnst/ R
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
c  flags indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,xerr,xnotd,xnotc
      data delmol/1.0d-4/
c
c  fill output fugacity array with zeros (final value for undefined
c  components and insurance against problems for others)
      do 100 i=1,ncmax
      f(i)=0.0d0
  100 continue
c
c  check for zero input density
      if (rho.lt.1.0d-20) then
        do 102 i=1,nc
        f(i)=0.0d0
  102   continue
        RETURN
      end if
c
      RTrho=R*t*rho
      if (nc.eq.1) then
c  pure component
        if (heos.eq.'BWR') then
c  pure fluid MBWR equation of state--use BWR-specific routines
          icomp=1
          Ar=ABWR(icomp,t,rho)
          p=PBWR(icomp,t,rho)
          f(1)=RTrho*exp(Ar/(R*t)+p/RTrho-1.0d0)
        else
c  for other models, use general REDK and PHIK routines
          icomp=1
          call REDK (icomp,t0,rho0)
          tau=t0/t
          del=rho/rho0
          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+del*phi01
          if (ABS(arg).lt.500.0d0) then
            f(1)=RTrho*exp(arg)
          else
            f(1)=xerr
          end if
        end if
      else
c
c  mixture
        deln=delmol       !initialize only
        do 200 i=1,nc
c  compute positive and negative increments to number of moles
        if (x(i).lt.delmol) then
c  special case--composition of component i is nearly zero
          deln=-x(i)
          delp=2.0d0*delmol-deln
        else if (x(i).gt.1.0d0-deln) then
c  special case--composition of component i is neary one (pure fluid)
          delp=1.0d0-x(i)
          deln=-2.0d0*delmol+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 120 j=1,nc
        xplus(j)=x(j)*delp1
        xminus(j)=x(j)*deln1
  120   continue
        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 compsition
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.500.0d0) then
          f(i)=x(i)*RTrho*exp(dadn)             !A is dimensionless
        else
          f(i)=xerr
        end if
  200   continue
      end if
c
      RETURN
      end                                              !subroutine FGCTY
c
c ======================================================================
c
      subroutine ACTVY (t,rho,x,gamma)
c
c  !temp--work in progress:  this routine returns gamma = 1.0
c
c  compute activity coefficients for each of the nc components of a
c  mixture by numerical differentiation
c
c  inputs:
c        t--temperature [K]
c      rho--molar density [mol/L]
c        x--composition array [mol frac]
c  output:
c    gamma--array (1..nc) of activity coefficients
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  12-18-95  MM, original version
c  02-27-96  MM, parameter n0=-ncmax to accomodate 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 switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: ACTVY
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      dimension x(ncmax),gamma(ncmax)
      dimension xplus(ncmax),xminus(ncmax)
      character*3 hpheq,heos,hmxeos,hmodcp
      common /NCOMP/ nc
      common /Gcnst/ R
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
      data delmol/1.0d-4/
c
c  fill output array with zeros for undefined components
      do 100 i=nc,ncmax
      gamma(i)=0.0d0
  100 continue
c
      if (nc.eq.1) then
c  pure component
        gamma(1)=1.0d0
      else
        deln=delmol       !initialize only
c  mixture
        do 200 i=1,nc
c  compute positive and negative increments to number of moles
        if (x(i).lt.delmol) then
c  special case--composition of component i is nearly zero
          deln=-x(i)
          delp=2.0d0*delmol-deln
        else if (x(i).gt.1.0d0-deln) then
c  special case--composition of component i is neary one (pure fluid)
          delp=1.0d0-x(i)
          deln=-2.0d0*delmol+delp
        else
c  general case:  deln < x(i) < 1 - deln
          delp=delmol
          deln=-delmol
        end if
c       write (*,*) ' ACTVY--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 120 j=1,nc
        xplus(j)=x(j)*delp1
        xminus(j)=x(j)*deln1
  120   continue
        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 at 'plus' and 'minus' density and compsition
c  could call general subroutines here, but more efficient to directly
c  call 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
        gamma(i)=1.0d0
  200   continue
      end if
c
      RETURN
      end                                              !subroutine ACTVY
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file prop_sub.f
c ======================================================================