c begin file core_PH0.f c c This file contains the functions implementing the ideal-gas part of c the reduced Helmholtz free energy form of the pure fluid equation of c state (the so-called "fundamental equation"). c c N.B. The Helmholtz energy consists of ideal and residual (real-gas) c terms; this routine calculates only the ideal part. c c This function computes pure component properties only. c c contained here are: c subroutine SETPH0 (nread,icomp,hcasno,ierr,herr) c function PH0PH0 (icomp,itau,idel,t,rho) c block data SAVPH0 c c these routines use the following common blocks from other files c common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) c common /CHAR/ htab,hnull c common /Gcnst/ R 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 parameter (nph0mx=10) !max number of terms in Cp0 polynomial c c ====================================================================== c ====================================================================== c subroutine SETPH0 (nread,icomp,hcasno,ierr,herr) c c set up working arrays for ideal-gas part of the Helmholtz energy c implements an expression of the form: c phi0 = Sum[ai*log(tau**ti)] + Sum[aj*tau**tj] c + Sum[ak*log(1-EXP(bk*tau))] c c inputs: c nread--file to read data from c <= 0 get data from block data (not currently implemented) c >0 read from logical unit nread (file should have already c been opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid; c zero and negative numbers designate ECS reference fluids c hcasno--CAS number of component icomp c (not req'd if reading from file--included to maintain c parallel structure with other routines) c c outputs: c ierr--error flag: 0 = successful c 1 = error (e.g. fluid not found) c herr--error string (character*255 variable if ierr<>0) c coefficients, etc. returned via arrays in common /xxxPH0/ c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-13-96 MM, original version, based (loosely) on SETCPP c 06-14-96 MM, fix bug in reading/storing coefficients c 11-13-97 MM, (re)initialize contents of /PH0SAV/ when a new fluid is read in c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nph0mx=10) !max number of terms in phi0 function parameter (mxph0=2) !max number of fluids in block data character*1 htab,hnull character*12 hcasno,hcas character*255 herr common /CHAR/ htab,hnull c commons associated with the nc components of current interest c ("working" commons and arrays) common /CASPH0/ hcas(mxph0) common /WNTPH0/ nlog(n0:nx),ntau(n0:nx),nexp(n0:nx) common /WLMPH0/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCFPH0/ ai(n0:nx,nph0mx),ti(n0:nx,nph0mx) common /PH0SAV/ ph0sav(n0:nx),ph1sav(n0:nx),ph2sav(n0:nx), & tsav(0:2,n0:nx),rhosav(0:2,n0:nx) c c (re)initialize contents of /FEQSAV/ when a new fluid is read in do 120 i=n0,nx ph0sav(i)=0.0d0 ph1sav(i)=0.0d0 ph2sav(i)=0.0d0 do 100 j=0,2 tsav(i,j)=0.0d0 rhosav(i,j)=0.0d0 100 continue 120 continue c if (nread.le.0) then c get coefficients from block data c identify specified fluid with entries in database via match of CAS no do 200 k=1,mxph0 if (hcasno.eq.hcas(k)) then c write (*,*) ' SETPH0--ERROR, block data read not implemented' ierr=1 herr=' SETPH0 ERROR--block data option not implemented'//hnull RETURN end if 200 continue ierr=1 herr=' SETPH0 ERROR--input fluid (block data) not found'//hnull RETURN else c read data from file c write (*,*) ' SETPH0--read component',icomp,' from unit',nread read (nread,*) tmin(icomp) !lower temperature limit read (nread,*) tmax(icomp) !upper temperature limit read (nread,*) pmax(icomp) !upper pressure limit read (nread,*) rhomax(icomp) !upper density limit c read number of terms for each of the various types read (nread,*) nlog(icomp),ntau(icomp),nexp(icomp) iterm=0 if (nlog(icomp).ge.1) then c read coefficients for terms of the form [ai*log(tau**ti)] do 240 i=1,nlog(icomp) iterm=iterm+1 read (nread,*) ai(icomp,iterm),ti(icomp,iterm) c write (*,*) ' SETPH0 log term: i,ai: ',iterm,ai(icomp,iterm) 240 continue end if if (ntau(icomp).ge.1) then c read coefficients for terms of the form [ai*tau**ti] do 260 j=1,ntau(icomp) iterm=iterm+1 read (nread,*) ai(icomp,iterm),ti(icomp,iterm) c write (*,*) ' SETPH0 log term: i,ai: ',iterm,ai(icomp,iterm) 260 continue end if if (nexp(icomp).ge.1) then c read coefficients for terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) do 280 k=1,nexp(icomp) iterm=iterm+1 read (nread,*) ai(icomp,iterm),ti(icomp,iterm) c write (*,*) ' SETPH0 log term: i,ai: ',iterm,ai(icomp,iterm) 280 continue end if c write (*,*) ' SETPH0--final coeff: ',ai(icomp,iterm) ierr=0 herr=hnull end if c RETURN 2003 format (a3) end !subroutine SETPH0 c c ====================================================================== c function PH0PH0 (icomp,itau,idel,t,rho) c c compute the ideal gas part of the reduced Helmholtz energy or a c derivative as functions of temperature and pressure; for c use with a Helmholtz-explicit equation of state c c inputs: c icomp--pointer specifying component (1..nc) c itau--flag specifying order of temperature derivative to calc c idel--flag specifying order of density derivative to calculate c (the density derivatives are not used in the calculation c of any property, and are not implemented) c when itau = 0 and idel = 0, compute A0/RT c when itau = 1 and idel = 0, 1st temperature derivative c when itau = 2 and idel = 0, 2nd temperature derivative c t--temperature (K) c rho--density (mol/L) c output (as function value): c ph0ph0--ideal-gas part of the Helmholtz energy in reduced form (A/RT) c derivatives (as specified by itau and idel) are multiplied c by the corresponding power of tau; i.e. when itau = 1, the c quantity returned is tau*d(ph0ph0)/d(tau) and when itau = 2, c the quantity returned is tau*tau*d2(ph0ph0)/d(tau)**2 c c N.B. While the real-gas part of the Helmholtz energy is calculated c in terms of dimensionless temperature and density, the ideal- c gas part is calculated in terms of absolute temperature and c density. (This distinction is necessary for mixtures.) c c The Helmholtz energy consists of ideal-gas and residual c (real-gas) terms; this routine calculates only the ideal part. c c This function computes pure component properties only. c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-13-96 MM, original version, based (loosely) on PH0CPP c 06-14-96 MM, add rhosav to /PH0SAV/ c 07-08-96 MM, change derivative outputs: tau*d(phi)/d(tau), etc c 08-20-97 MM, call ERRMSG if itau out of range; drop idel=idel 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) parameter (nph0mx=10) !max number of terms in phi0 function character*1 htab,hnull character*255 herr common /CHAR/ htab,hnull common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /Gcnst/ R c commons associated with the nc components of current interest c ("working" commons and arrays) common /WNTPH0/ nlog(n0:nx),ntau(n0:nx),nexp(n0:nx) common /WLMPH0/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCFPH0/ ai(n0:nx,nph0mx),ti(n0:nx,nph0mx) c saved values from previous calls to this routine common /PH0SAV/ ph0sav(n0:nx),ph1sav(n0:nx),ph2sav(n0:nx), & tsav(0:2,n0:nx),rhosav(0:2,n0:nx) c c compute reduced Helmholtz; first check if t same as previous call c PH0PH0=0.0d0 !initialize in case of error if (itau.eq.0) then if (abs(t-tsav(0,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(0,icomp)).lt.1.0d-10) then PH0PH0=ph0sav(icomp) c write (*,*) ' PH0PH0--using saved phi for itau = 0, t =',t else call REDK (icomp,t0,D0) tau=t0/t del=rho/D0 c write (*,1030) icomp,t,rho,t0,D0,tau,del c1030 format (1x,'PH0PH0--icomp,t,rho,t0,D0,tau,del: ',i4,6e14.6) iterm=0 phisum=log(del) c & +href(icomp)/R/t-sref(icomp)/R !see U. Idaho class notes if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] taulog=log(tau) do 40 i=1,nlog(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm)*taulog c write (*,1040) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1040 format (' PH0PH0--log term iterm,ai,ti,phisum: ',i3,3f12.6) 40 continue end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] do 60 j=1,ntau(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*tau**ti(icomp,iterm) c write (*,1060) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1060 format (' PH0PH0--tau term iterm,ai,ti,phisum: ',i3,3f12.6) 60 continue end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) do 80 k=1,nexp(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm) & *log(1.0d0-EXP(ti(icomp,iterm)*tau)) c write (*,1080) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1080 format (' PH0PH0--exp term iterm,ai,ti,phisum: ',i3,3f12.6) 80 continue end if PH0PH0=phisum c save information for possible use on next call to function tsav(0,icomp)=t rhosav(0,icomp)=rho ph0sav(icomp)=phisum end if c c compute derivative w.r.t. tau (dimensionless temperature) c else if (itau.eq.1) then if (abs(t-tsav(1,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(1,icomp)).lt.1.0d-10) then PH0PH0=ph1sav(icomp) else call REDK (icomp,t0,D0) tau=t0/t iterm=0 phisum=0.0d0 c phisum=href(icomp)/(R*t0) !see U. Idaho class notes if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] tauinv=1.0d0/tau do 140 i=1,nlog(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm)*tauinv 140 continue end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] iterm=iterm+1 do 160 j=1,ntau(icomp) phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *tau**(ti(icomp,iterm)-1.0d0) 160 continue end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) iterm=iterm+1 do 180 k=1,nexp(icomp) exptau=EXP(ti(icomp,iterm)*tau) phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm) & *exptau/(1.0d0-exptau) 180 continue end if c save information for possible use on next call to function tsav(1,icomp)=t rhosav(1,icomp)=rho ph1sav(icomp)=phisum*tau !return tau*d(ph0ph0)/d(tau) PH0PH0=ph1sav(icomp) end if c c compute 2nd derivative w.r.t. tau (dimensionless temperature) c else if (itau.eq.2) then if (abs(t-tsav(2,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(2,icomp)).lt.1.0d-10) then PH0PH0=ph2sav(icomp) else call REDK (icomp,t0,D0) tau=t0/t iterm=0 phisum=0.0d0 if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] tauin2=1.0d0/tau**2 do 240 i=1,nlog(icomp) iterm=iterm+1 phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)*tauin2 240 continue end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] iterm=iterm+1 do 260 j=1,ntau(icomp) phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *(ti(icomp,iterm)-1.0d0)*tau**(ti(icomp,iterm)-2.0d0) 260 continue end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) iterm=iterm+1 do 280 k=1,nexp(icomp) exptau=EXP(ti(icomp,iterm)*tau) phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)**2 & *exptau/(1.0d0-exptau)**2 280 continue end if c save information for possible use on next call to function tsav(2,icomp)=t rhosav(2,icomp)=rho c return tau**2*d2(ph0ph0)/d(tau**2) ph2sav(icomp)=phisum*tau*tau PH0PH0=ph2sav(icomp) end if else c c invalid value of itau c ierr=99 write (herr,1099) itau,idel,hnull 1099 format ('[PH0PH0 warning] invalid input; itau =',i4,'; idel =', & i4,a1) call ERRMSG (ierr,herr) PH0PH0=0.0d0 end if c c write (*,*) ' PH0PH0: output phi: ',ph0ph0 c RETURN end !function PH0PH0 c c ====================================================================== c block data SAVPH0 c c initialize "saved" variables used in function PH0PH0 c c explanation of commons and constituent arrays c /PH0SAV/ used to save information between calls to PH0PH0 c ph0sav(i): phi0 (A0/RT) for component i on last call c ph1sav(i): d(phi0)/d(tau) for component i on last call c ph2sav(i): d2(phi0)/d(tau)2 for component i on last call c tsav(itau,i): previous temperature for given value of itau c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-13-96 MM, original version, based (loosely) on SAVCPP c 06-14-96 MM, add rhosav to /PH0SAV/ c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nxsav=nx-n0+1) parameter (mxph0=2) !max number of fluids in block data character*12 hcas c common /PH0SAV/ ph0sav(n0:nx),ph1sav(n0:nx),ph2sav(n0:nx), & tsav(0:2,n0:nx),rhosav(0:2,n0:nx) common /CASPH0/ hcas(mxph0) c data ph0sav /nxsav*0.0d0/ data ph1sav /nxsav*0.0d0/ data ph2sav /nxsav*0.0d0/ data tsav /nxsav*0.0d0,nxsav*0.0d0,nxsav*0.0d0/ data rhosav /nxsav*0.0d0,nxsav*0.0d0,nxsav*0.0d0/ c initialize common associated with coefficients stored in block data c (not implemented) data hcas /mxph0*'not defined'/ c end !block data SAVPH0 c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_PH0.f c ======================================================================