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 ======================================================================