c begin file core_ANC.f c c This file contains core routines for the ancillary equations for vapor c pressure, saturated liquid density and saturated vapor density. c c contained here are: c subroutine SETPS (nread,icomp,ierr,herr) c subroutine PSATT (t,x,p,ierr,herr) c subroutine PSATK (icomp,t,p,ierr,herr) c subroutine DPSATK (icomp,t,dpdt,ierr,herr) c subroutine D2PSTK (icomp,t,d2pdt2,ierr,herr) c subroutine TSATP (p,x,t,ierr,herr) c subroutine SETPL (nread,icomp,ierr,herr) c subroutine PLSATT (t,x,p,ierr,herr) c subroutine PLSATK (icomp,t,p,ierr,herr) c subroutine SETDL (nread,icomp,ierr,herr) c subroutine DLSATT (t,x,d,ierr,herr) c subroutine DLSATK (icomp,t,d,ierr,herr) c subroutine SETDV (nread,icomp,ierr,herr) c subroutine DVSATT (t,x,d,ierr,herr) c subroutine DVSATK (icomp,t,d,ierr,herr) c subroutine TSATD (d,x,t,ierr,herr) c subroutine SOLVEA (iflag,pd,x,t,ierr,herr) c c ====================================================================== c ====================================================================== c subroutine SETPS (nread,icomp,ierr,herr) c c set up working arrays for use with the vapor pressure ancillary equation c c inputs: c nread--file to read data from (file should have already been c opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c c outputs: c ierr--error flag: 0 = successful c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (npsk=15) !max number of terms in summation character*1 htab,hnull character*3 hps,hpsk character*255 herr common /HCHAR/ htab,hnull common /PSMOD/ hps,hpsk(n0:nx) common /WLMPS/ pstmin(n0:nx),pstmax(n0:nx) common /WNTPS/ nps1(n0:nx),nps2(n0:nx),nps3(n0:nx), & nps4(n0:nx),nps5(n0:nx),nps6(n0:nx) common /WCFPS/ psk(n0:nx,npsk),psexp(n0:nx,npsk) common /WRDPS/ pstrd(n0:nx),psprd(n0:nx) c c read data from file c write (*,*) ' SETPS--read component',icomp,' from unit',nread read (nread,*) pstmin(icomp) !lower temperature limit read (nread,*) pstmax(icomp) !upper temperature limit c the pressure and density limit are not presently used, c but are contained in the file for consistency and possible future use; c skip over them in reading the file read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) pstrd(icomp),psprd(icomp)!reducing parameters read (nread,*) nps1(icomp),nps2(icomp),nps3(icomp), & nps4(icomp),nps5(icomp),nps6(icomp) if (nps1(icomp).gt.0) then do k=1,nps1(icomp) read (nread,*) psk(icomp,k),psexp(icomp,k) enddo endif ierr=0 herr=' ' c RETURN end !subroutine SETPS c c ====================================================================== c subroutine PSATT (t,x,p,ierr,herr) c c compute mixture or pure fluid vapor pressure with appropriate c ancillary equation c c inputs: c t--temperature (K) c x--composition [array of mol frac] c output: c p--vapor pressure [kPa] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown vapor pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture character*1 htab,hnull character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull c ierr=0 herr=' ' p=0.0d0 if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c special case--pure component call PSATK (icomp,t,p,ierr,herr) else pp=0.0d0 do i=1,nc call PSATK (i,t,p,ierr,herr) !Add Raoult's Law here pp=pp+x(i)*p enddo p=pp end if c write (*,1200) t,p c1200 format (' PSATT--t,p: ',2f11.6) c RETURN end !subroutine PSATT c c ====================================================================== c subroutine PSATK (icomp,t,p,ierr,herr) c c compute pure fluid vapor pressure with the appropriate ancillary equation c c inputs: c icomp--component i c t--temperature (K) c output: c p--vapor pressure [kPa] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown vapor pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (npsk=15) !max number of terms in summation character*1 htab,hnull character*3 hps,hpsk character*255 herr common /HCHAR/ htab,hnull common /PSMOD/ hps,hpsk(n0:nx) common /WLMPS/ pstmin(n0:nx),pstmax(n0:nx) common /WNTPS/ nps1(n0:nx),nps2(n0:nx),nps3(n0:nx), & nps4(n0:nx),nps5(n0:nx),nps6(n0:nx) common /WCFPS/ psk(n0:nx,npsk),psexp(n0:nx,npsk) common /WRDPS/ pstrd(n0:nx),psprd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' p=0.0d0 if (t.lt.pstmin(icomp)) then p=ptp(icomp) ierr=1 write (herr,1097) t,pstmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[PSATK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then p=pcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[PSATK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'PS' in the .fld file i=0 i=ICHAR(hpsk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hpsk(icomp)(1:2).eq.'PS' .and. i.gt.0) then tr=ABS(1.d0-t/pstrd(icomp)) if (MOD(i,2).eq.0) tr=SQRT(tr) !Even values of i pr=0.0d0 if (nps1(icomp).ne.0) then do k=1,nps1(icomp) pr=pr+psk(icomp,k)*tr**psexp(icomp,k) enddo endif if (i.eq.1 .or. i.eq.2) pr=1.d0+pr if (i.eq.3 .or. i.eq.4) pr=EXP(pr) if (i.eq.5 .or. i.eq.6) pr=EXP(pstrd(icomp)/t*pr) p=psprd(icomp)*pr c do not return error message if fluid contains no vapor pressure line elseif (hpsk(icomp).eq.'NBS') then p=pcrit(icomp) ierr=2 herr='[PSATK error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) else p=0.0d0 ierr=100 write (herr,1099) hpsk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[PSATK error 100] ', & 'unknown vapor pressure equation: (',a3,')',a1) end if c write (*,1200) icomp,t,p c1200 format (' PSATK--icomp,t,p: ',i4,2f11.6) c RETURN end !subroutine PSATK c c ====================================================================== c subroutine DPSATK (icomp,t,dpdt,ierr,herr) c c compute the first derivative of the vapor pressure with the c appropriate equation c c inputs: c icomp--component i c t--temperature (K) c output: c dpdt--dp/dt [kPa/K] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown vapor pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 05-30-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (npsk=15) !max number of terms in summation character*1 htab,hnull character*3 hps,hpsk character*255 herr common /HCHAR/ htab,hnull common /PSMOD/ hps,hpsk(n0:nx) common /WLMPS/ pstmin(n0:nx),pstmax(n0:nx) common /WNTPS/ nps1(n0:nx),nps2(n0:nx),nps3(n0:nx), & nps4(n0:nx),nps5(n0:nx),nps6(n0:nx) common /WCFPS/ psk(n0:nx,npsk),psexp(n0:nx,npsk) common /WRDPS/ pstrd(n0:nx),psprd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' dpdt=0.0d0 if (t.lt.pstmin(icomp)) then dpdt=ptp(icomp) ierr=1 write (herr,1097) t,pstmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[DPSATK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then dpdt=pcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[DPSATK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'PS' in the .fld file i=0 i=ICHAR(hpsk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hpsk(icomp)(1:2).eq.'PS' .and. i.gt.0) then tr=ABS(1.d0-t/pstrd(icomp)) if (MOD(i,2).eq.0) tr=SQRT(tr) !Even values of i pr=0.0d0 dpr=0.0d0 if (nps1(icomp).ne.0) then do k=1,nps1(icomp) pe=psexp(icomp,k) pr=pr+psk(icomp,k)*tr**pe if (MOD(i,2).eq.1) then !odd values of i dpr=dpr-pe*psk(icomp,k)*tr**(pe-1.d0)/pstrd(icomp) else !even values of i dpr=dpr-pe/2.d0*psk(icomp,k)*tr**(pe-2.d0)/pstrd(icomp) endif enddo endif tr=pstrd(icomp)/t if (i.eq.1 .or. i.eq.2) pr=dpr if (i.eq.3 .or. i.eq.4) pr=EXP(pr)*dpr if (i.eq.5 .or. i.eq.6) pr=EXP(tr*pr)*tr*(dpr-pr/t) dpdt=psprd(icomp)*pr c do not return error message if fluid contains no vapor pressure line elseif (hpsk(icomp).eq.'NBS') then dpdt=pcrit(icomp) ierr=2 herr='[PSATK error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) else dpdt=0.0d0 ierr=100 write (herr,1099) hpsk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[DPSATK error 100] ', & 'unknown vapor pressure equation: (',a3,')',a1) end if c write (*,1200) icomp,t,dpdt c1200 format (' DPSATK--icomp,t,dp/dt: ',i4,2f11.6) c RETURN end !subroutine DPSATK c c ====================================================================== c subroutine D2PSTK (icomp,t,d2pdt2,ierr,herr) c c compute the second derivative of the vapor pressure with the c appropriate equation c c inputs: c icomp--component i c t--temperature (K) c output: c d2pdt2--d2p/dt2 [kPa/K^2] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown vapor pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 05-30-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (npsk=15) !max number of terms in summation character*1 htab,hnull character*3 hps,hpsk character*255 herr common /HCHAR/ htab,hnull common /PSMOD/ hps,hpsk(n0:nx) common /WLMPS/ pstmin(n0:nx),pstmax(n0:nx) common /WNTPS/ nps1(n0:nx),nps2(n0:nx),nps3(n0:nx), & nps4(n0:nx),nps5(n0:nx),nps6(n0:nx) common /WCFPS/ psk(n0:nx,npsk),psexp(n0:nx,npsk) common /WRDPS/ pstrd(n0:nx),psprd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' d2pdt2=0.0d0 if (t.lt.pstmin(icomp)) then d2pdt2=ptp(icomp) ierr=1 write (herr,1097) t,pstmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[D2PSTK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then d2pdt2=pcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[D2PSTK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'PS' in the .fld file i=0 i=ICHAR(hpsk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hpsk(icomp)(1:2).eq.'PS' .and. i.gt.0) then tr=ABS(1.d0-t/pstrd(icomp)) if (MOD(i,2).eq.0) tr=SQRT(tr) !Even values of i pr=0.0d0 dpr=0.0d0 d2pr=0.0d0 if (nps1(icomp).ne.0) then do k=1,nps1(icomp) pe=psexp(icomp,k) pr=pr+psk(icomp,k)*tr**pe if (MOD(i,2).eq.1) then !odd values of i dpr=dpr-pe*psk(icomp,k)*tr**(pe-1.d0)/pstrd(icomp) d2pr=d2pr+pe*(pe-1.d0)*psk(icomp,k)*tr**(pe-2.d0) & /pstrd(icomp)**2 else !even values of i dpr=dpr-pe/2.d0*psk(icomp,k)*tr**(pe-2.d0)/pstrd(icomp) d2pr=d2pr+pe*(pe-2.d0)/4.d0*psk(icomp,k)*tr**(pe-4.d0) & /pstrd(icomp)**2 endif enddo endif tr=pstrd(icomp)/t if (i.eq.1 .or. i.eq.2) pr=d2pr if (i.eq.3 .or. i.eq.4) pr=EXP(pr)*(dpr**2+d2pr) if (i.eq.5 .or. i.eq.6) pr=EXP(tr*pr)* & ((tr*(dpr-pr/t))**2+2.d0*tr/t**2*pr-2.d0*tr/t*dpr+tr*d2pr) d2pdt2=psprd(icomp)*pr c do not return error message if fluid contains no vapor pressure line elseif (hpsk(icomp).eq.'NBS') then d2pdt2=pcrit(icomp) ierr=2 herr='[D2PSTK error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) else d2pdt2=0.0d0 ierr=100 write (herr,1099) hpsk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[D2PSTK error 100] ', & 'unknown vapor pressure equation: (',a3,')',a1) end if c write (*,1200) icomp,t,d2pdt2 c1200 format (' D2PSTK--icomp,t,dp/dt: ',i4,2f11.6) c RETURN end !subroutine D2PSTK c c ====================================================================== c subroutine TSATP (p,x,t,ierr,herr) c c compute the vapor temperature as a function of pressure. c c inputs: c p--vapor pressure [kPa] c x--composition [array of mol frac] c output: c t--temperature [K] c ierr--error flag: 0 = successful c 1 = error: ppc c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) character*1 htab,hnull character*3 hps,hpsk character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull common /PSMOD/ hps,hpsk(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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) if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c if (hpsk(icomp).eq.'NBS') then !Return if no vapor pressure eq. t=tcrit(icomp) ierr=2 herr='[TSATP error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) RETURN endif if (p.lt.ptp(icomp)) then ierr=1 t=ttp(icomp) write (herr,1148) p/1000.0d0,ptp(icomp)/1000.0d0,hnull call ERRMSG (ierr,herr) 1148 format ('[TSATP error 1] ', & 'pressure less than triple point pressure; P =', & g11.5,' MPa, Ptp =',g11.5,' MPa.',a1) RETURN endif if (p.gt.pcrit(icomp)) then ierr=141 t=tcrit(icomp) write (herr,1149) p/1000.0d0,pcrit(icomp)/1000.0d0,hnull call ERRMSG (ierr,herr) 1149 format ('[TSATP error 141] ', & 'pressure greater than critical pressure; P =', & g11.5,' MPa, Pc =',g11.5,' MPa.',a1) RETURN endif endif c call SOLVEA (0,p,x,t,ierr,herr) c RETURN end !subroutine TSATP c c ====================================================================== c subroutine SETPL (nread,icomp,ierr,herr) c c set up working arrays for use with the liquid pressure ancillary equation c c inputs: c nread--file to read data from (file should have already been c opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c c outputs: c ierr--error flag: 0 = successful c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado c 09-18-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (nplk=15) !max number of terms in summation character*1 htab,hnull character*3 hpl,hplk character*255 herr common /HCHAR/ htab,hnull common /PLMOD/ hpl,hplk(n0:nx) common /WLMPL/ pltmin(n0:nx),pltmax(n0:nx) common /WNTPL/ npl1(n0:nx),npl2(n0:nx),npl3(n0:nx), & npl4(n0:nx),npl5(n0:nx),npl6(n0:nx) common /WCFPL/ plk(n0:nx,nplk),plexp(n0:nx,nplk) common /WRDPL/ pltrd(n0:nx),plprd(n0:nx) c c read data from file c write (*,*) ' SETPL--read component',icomp,' from unit',nread read (nread,*) pltmin(icomp) !lower temperature limit read (nread,*) pltmax(icomp) !upper temperature limit c the pressure and density limit are not presently used, c but are contained in the file for consistency and possible future use; c skip over them in reading the file read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) pltrd(icomp),plprd(icomp)!reducing parameters read (nread,*) npl1(icomp),npl2(icomp),npl3(icomp), & npl4(icomp),npl5(icomp),npl6(icomp) if (npl1(icomp).gt.0) then do k=1,npl1(icomp) read (nread,*) plk(icomp,k),plexp(icomp,k) enddo endif ierr=0 herr=' ' c RETURN end !subroutine SETPL c c ====================================================================== c subroutine PLSATT (t,x,p,ierr,herr) c c compute mixture or pure fluid liquid pressure with appropriate c ancillary equation c (used only for pseudo-pure fluids) c c inputs: c t--temperature (K) c x--composition [array of mol frac] c output: c p--liquid pressure [kPa] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown liquid pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 09-18-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture character*1 htab,hnull character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull c ierr=0 herr=' ' p=0.0d0 if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c special case--pure component call PLSATK (icomp,t,p,ierr,herr) else pp=0.0d0 do i=1,nc call PLSATK (i,t,p,ierr,herr) !Add Raoult's Law here pp=pp+x(i)*p enddo p=pp end if c write (*,1200) t,p c1200 format (' PLSATT--t,p: ',2f11.6) c RETURN end !subroutine PLSATT c c ====================================================================== c subroutine PLSATK (icomp,t,p,ierr,herr) c c compute pure fluid liquid pressure with the appropriate ancillary equation c (used only for pseudo-pure fluids) c c inputs: c icomp--component i c t--temperature (K) c output: c p--liquid pressure [kPa] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown liquid pressure equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (nplk=15) !max number of terms in summation character*1 htab,hnull character*3 hpl,hplk character*255 herr common /HCHAR/ htab,hnull common /PLMOD/ hpl,hplk(n0:nx) common /WLMPL/ pltmin(n0:nx),pltmax(n0:nx) common /WNTPL/ npl1(n0:nx),npl2(n0:nx),npl3(n0:nx), & npl4(n0:nx),npl5(n0:nx),npl6(n0:nx) common /WCFPL/ plk(n0:nx,nplk),plexp(n0:nx,nplk) common /WRDPL/ pltrd(n0:nx),plprd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' p=0.0d0 if (t.lt.pltmin(icomp)) then p=ptp(icomp) ierr=1 write (herr,1097) t,pltmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[PLSATK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then p=pcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[PLSATK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'PL' in the .fld file i=0 i=ICHAR(hplk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hplk(icomp)(1:2).eq.'PL' .and. i.gt.0) then tr=ABS(1.d0-t/pltrd(icomp)) if (MOD(i,2).eq.0) tr=SQRT(tr) !Even values of i pr=0.0d0 if (npl1(icomp).ne.0) then do k=1,npl1(icomp) pr=pr+plk(icomp,k)*tr**plexp(icomp,k) enddo endif if (i.eq.1 .or. i.eq.2) pr=1.d0+pr if (i.eq.3 .or. i.eq.4) pr=EXP(pr) if (i.eq.5 .or. i.eq.6) pr=EXP(pltrd(icomp)/t*pr) p=plprd(icomp)*pr c do not return error message if fluid contains no liquid pressure line elseif (hplk(icomp).eq.'NBS') then p=pcrit(icomp) ierr=2 herr='[PLSATK error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) else p=0.0d0 ierr=100 write (herr,1099) hplk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[PLSATK error 100] ', & 'unknown liquid pressure equation: (',a3,')',a1) end if c write (*,1200) icomp,t,p c1200 format (' PLSATK--icomp,t,p: ',i4,2f11.6) c RETURN end !subroutine PLSATK c c ====================================================================== c subroutine TSATPL (p,x,t,ierr,herr) c c compute the liquid temperature as a function of liquid pressure. c (used only for pseudo-pure fluids) c c inputs: c p--liquid pressure [kPa] c x--composition [array of mol frac] c output: c t--temperature [K] c ierr--error flag: 0 = successful c 1 = error: ppc c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO c 09-19-02 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) character*1 htab,hnull character*3 hpl,hplk character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull common /PLMOD/ hpl,hplk(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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) if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c if (hplk(icomp).eq.'NBS') then !Return if no liquid pressure eq. t=tcrit(icomp) ierr=2 herr='[TSATPL error 2] vapor pressure equation not available' call ERRMSG (ierr,herr) RETURN endif if (p.lt.ptp(icomp)) then ierr=1 t=ttp(icomp) write (herr,1148) p/1000.0d0,ptp(icomp)/1000.0d0,hnull call ERRMSG (ierr,herr) 1148 format ('[TSATPL error 1] ', & 'pressure less than triple point pressure; P =', & g11.5,' MPa, Ptp =',g11.5,' MPa.',a1) RETURN endif if (p.gt.pcrit(icomp)) then ierr=141 t=tcrit(icomp) write (herr,1149) p/1000.0d0,pcrit(icomp)/1000.0d0,hnull call ERRMSG (ierr,herr) 1149 format ('[TSATPL error 141] ', & 'pressure greater than critical pressure; P =', & g11.5,' MPa, Pc =',g11.5,' MPa.',a1) RETURN endif endif c call SOLVEA (4,p,x,t,ierr,herr) c RETURN end !subroutine TSATPL c c ====================================================================== c subroutine SETDL (nread,icomp,ierr,herr) c c set up working arrays for use with the saturated liquid density c ancillary equation c c inputs: c nread--file to read data from (file should have already been c opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c c outputs: c ierr--error flag: 0 = successful c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado c 09-18-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (ndlk=15) !max number of terms in summation character*1 htab,hnull character*3 hdl,hdlk character*255 herr common /HCHAR/ htab,hnull common /DLMOD/ hdl,hdlk(n0:nx) common /WLMDL/ dltmin(n0:nx),dltmax(n0:nx) common /WNTDL/ ndl1(n0:nx),ndl2(n0:nx),ndl3(n0:nx), & ndl4(n0:nx),ndl5(n0:nx),ndl6(n0:nx) common /WCFDL/ dlk(n0:nx,ndlk),dlexp(n0:nx,ndlk) common /WRDDL/ dltrd(n0:nx),dldrd(n0:nx) c c read data from file c write (*,*) ' SETDL--read component',icomp,' from unit',nread read (nread,*) dltmin(icomp) !lower temperature limit read (nread,*) dltmax(icomp) !upper temperature limit c the pressure and density limit are not presently used, c but are contained in the file for consistency and possible future use; c skip over them in reading the file read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) dltrd(icomp),dldrd(icomp)!reducing parameters read (nread,*) ndl1(icomp),ndl2(icomp),ndl3(icomp), & ndl4(icomp),ndl5(icomp),ndl6(icomp) if (ndl1(icomp).gt.0) then do k=1,ndl1(icomp) read (nread,*) dlk(icomp,k),dlexp(icomp,k) enddo endif ierr=0 herr=' ' c RETURN end !subroutine SETDL c c ====================================================================== c subroutine DLSATT (t,x,d,ierr,herr) c c compute mixture or pure fluid saturated liquid density with appropriate c ancillary equation c c inputs: c t--temperature (K) c x--composition [array of mol frac] c output: c d--saturated liquid density [mol/L] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown saturated liquid density equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture character*1 htab,hnull character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull c ierr=0 herr=' ' d=0.0d0 if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c special case--pure component call DLSATK (icomp,t,d,ierr,herr) else dd=0.0d0 do i=1,nc call DLSATK (i,t,d,ierr,herr) !Add Raoult's Law here dd=dd+x(i)*d enddo d=dd end if c write (*,1200) t,d c1200 format (' DLSATT--t,d: ',2f11.6) c RETURN end !subroutine DLSATT c c ====================================================================== c subroutine DLSATK (icomp,t,d,ierr,herr) c c compute pure fluid saturated liquid density with appropriate equation c c inputs: c icomp--component i c t--temperature (K) c output: c d--saturated liquid density [mol/L] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown saturated liquid density equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (ndlk=15) !max number of terms in summation character*1 htab,hnull character*3 hdl,hdlk character*255 herr common /HCHAR/ htab,hnull common /DLMOD/ hdl,hdlk(n0:nx) common /WLMDL/ dltmin(n0:nx),dltmax(n0:nx) common /WNTDL/ ndl1(n0:nx),ndl2(n0:nx),ndl3(n0:nx), & ndl4(n0:nx),ndl5(n0:nx),ndl6(n0:nx) common /WCFDL/ dlk(n0:nx,ndlk),dlexp(n0:nx,ndlk) common /WRDDL/ dltrd(n0:nx),dldrd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' d=0.0d0 if (t.lt.dltmin(icomp)) then d=dtp(icomp) ierr=1 write (herr,1097) t,dltmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[DLSATK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then d=dcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[DLSATK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'DL' in the .fld file i=0 i=ICHAR(hdlk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hdlk(icomp)(1:2).eq.'DL' .and. i.gt.0) then tr=ABS(1.d0-t/dltrd(icomp)) if (MOD(i,2).eq.0) tr=tr**(1.d0/3.d0) !Even values of i dr=0.0d0 if (ndl1(icomp).ne.0) then do k=1,ndl1(icomp) dr=dr+dlk(icomp,k)*tr**dlexp(icomp,k) enddo endif if (i.eq.1 .or. i.eq.2) dr=1.d0+dr if (i.eq.3 .or. i.eq.4) dr=EXP(dr) if (i.eq.5 .or. i.eq.6) dr=EXP(dltrd(icomp)/t*dr) d=dldrd(icomp)*dr c do not return error message if fluid contains no saturated liquid density equation elseif (hdlk(icomp).eq.'NBS') then d=dcrit(icomp) ierr=2 herr='[DLSATK error 2] liquid density equation not available' call ERRMSG (ierr,herr) else d=0.0d0 ierr=100 write (herr,1099) hdlk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[DLSATK error 100] ', & 'unknown saturated liquid density equation: (',a3,')',a1) end if c write (*,1200) icomp,t,d c1200 format (' DLSATK--icomp,t,d: ',i4,2f11.6) c RETURN end !subroutine DLSATK c c ====================================================================== c subroutine SETDV (nread,icomp,ierr,herr) c c set up working arrays for use with the saturated vapor density c ancillary equation c c inputs: c nread--file to read data from (file should have already been c opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c c outputs: c ierr--error flag: 0 = successful c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (ndvk=15) !max number of terms in summation character*1 htab,hnull character*3 hdv,hdvk character*255 herr common /HCHAR/ htab,hnull common /DVMOD/ hdv,hdvk(n0:nx) common /WLMDV/ dvtmin(n0:nx),dvtmax(n0:nx) common /WNTDV/ ndv1(n0:nx),ndv2(n0:nx),ndv3(n0:nx), & ndv4(n0:nx),ndv5(n0:nx),ndv6(n0:nx) common /WCFDV/ dvk(n0:nx,ndvk),dvexp(n0:nx,ndvk) common /WRDDV/ dvtrd(n0:nx),dvdrd(n0:nx) c c read data from file c write (*,*) ' SETDV--read component',icomp,' from unit',nread read (nread,*) dvtmin(icomp) !lower temperature limit read (nread,*) dvtmax(icomp) !upper temperature limit c the pressure and density limit are not presently used, c but are contained in the file for consistency and possible future use; c skip over them in reading the file read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) dvtrd(icomp),dvdrd(icomp)!reducing parameters read (nread,*) ndv1(icomp),ndv2(icomp),ndv3(icomp), & ndv4(icomp),ndv5(icomp),ndv6(icomp) if (ndv1(icomp).gt.0) then do k=1,ndv1(icomp) read (nread,*) dvk(icomp,k),dvexp(icomp,k) enddo endif ierr=0 herr=' ' c RETURN end !subroutine SETDV c c ====================================================================== c subroutine DVSATT (t,x,d,ierr,herr) c c compute mixture or pure fluid saturated vapor density with appropriate c ancillary equation c c inputs: c t--temperature (K) c x--composition [array of mol frac] c output: c d--saturated vapor density [mol/L] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown saturated vapor density equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture character*1 htab,hnull character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull c ierr=0 herr=' ' d=0.0d0 if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c special case--pure component call DVSATK (icomp,t,d,ierr,herr) else dd=0.0d0 do i=1,nc call DVSATK (i,t,d,ierr,herr) !Add Raoult's Law here dd=dd+x(i)*d enddo d=dd end if c write (*,1200) t,d c1200 format (' DVSATT--t,d: ',2f11.6) c RETURN end !subroutine DVSATT c c ====================================================================== c subroutine DVSATK (icomp,t,d,ierr,herr) c c compute pure fluid saturated vapor density with appropriate equation c c inputs: c icomp--component i c t--temperature (K) c output: c d--saturated vapor density [mol/L] c ierr--error flag: 0 = successful c 1 = error: TTc c 100 = error: unknown saturated vapor density equation c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) parameter (ndvk=15) !max number of terms in summation character*1 htab,hnull character*3 hdv,hdvk character*255 herr common /HCHAR/ htab,hnull common /DVMOD/ hdv,hdvk(n0:nx) common /WLMDV/ dvtmin(n0:nx),dvtmax(n0:nx) common /WNTDV/ ndv1(n0:nx),ndv2(n0:nx),ndv3(n0:nx), & ndv4(n0:nx),ndv5(n0:nx),ndv6(n0:nx) common /WCFDV/ dvk(n0:nx,ndvk),dvexp(n0:nx,ndvk) common /WRDDV/ dvtrd(n0:nx),dvdrd(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 ierr=0 herr=' ' d=0.0d0 if (t.lt.dvtmin(icomp)) then d=dtp(icomp) ierr=1 write (herr,1097) t,dvtmin(icomp),hnull call ERRMSG (ierr,herr) 1097 format ('[DVSATK error 1] ', & 'temperature less than triple point temperature; T =', & g11.5,' K, Ttp =',g11.5,' K.',a1) RETURN endif c if (t.gt.tcrit(icomp)) then d=dcrit(icomp) ierr=121 write (herr,1098) t,tcrit(icomp),hnull call ERRMSG (ierr,herr) 1098 format ('[DVSATK error 121] ', & 'temperature greater than critical temperature; T =', & g11.5,' K, Tc =',g11.5,' K.',a1) RETURN endif c c i is the value following 'DV' in the .fld file i=0 i=ICHAR(hdvk(icomp)(3:3))-48 if (i.lt.0) i=0 if (i.gt.9) i=0 if (hdvk(icomp)(1:2).eq.'DV' .and. i.gt.0) then tr=ABS(1.d0-t/dvtrd(icomp)) if (MOD(i,2).eq.0) tr=tr**(1.d0/3.d0) !Even values of i dr=0.0d0 if (ndv1(icomp).ne.0) then do k=1,ndv1(icomp) dr=dr+dvk(icomp,k)*tr**dvexp(icomp,k) enddo endif if (i.eq.1 .or. i.eq.2) dr=1.d0+dr if (i.eq.3 .or. i.eq.4) dr=EXP(dr) if (i.eq.5 .or. i.eq.6) dr=EXP(dvtrd(icomp)/t*dr) d=dvdrd(icomp)*dr c do not return error message if fluid contains no saturated vapor density equation elseif (hdvk(icomp).eq.'NBS') then d=dcrit(icomp) ierr=2 herr='[DLSATK error 2] liquid density equation not available' call ERRMSG (ierr,herr) else d=0.0d0 ierr=100 write (herr,1099) hdvk(icomp),hnull call ERRMSG (ierr,herr) 1099 format ('[DVSATK error 100] ', & 'unknown saturated vapor density equation: (',a3,')',a1) end if c write (*,1200) icomp,t,d c1200 format (' DVSATK--icomp,t,d: ',i4,2f11.6) c RETURN end !subroutine DVSATK c c ====================================================================== c subroutine TSATD (d,x,t,ierr,herr) c c compute the saturated temperature as a function of saturated density. c c inputs: c d--saturated density [mol/L] c x--composition [array of mol frac] c output: c t--temperature [K] c ierr--error flag: 0 = successful c 2 = error: ppc c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chem Properties Div, Boulder, CO c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-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) character*1 htab,hnull character*3 hdl,hdlk,hdv,hdvk character*255 herr dimension x(ncmax) common /NCOMP/ nc,ic common /HCHAR/ htab,hnull common /DLMOD/ hdl,hdlk(n0:nx) common /DVMOD/ hdv,hdvk(n0:nx) common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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 i=1 if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c Return if no saturated ancillary equation t=ttp(icomp) if (hdlk(icomp).eq.'NBS' .and. d.gt.dcrit(icomp)) RETURN if (hdvk(icomp).eq.'NBS' .and. d.lt.dcrit(icomp)) RETURN c if (d.lt.dtpv(1)) then !dtpv has not been set yet! c ierr=2 c write (herr,1148) d,dtpv(1),hnull c call ERRMSG (ierr,herr) c1148 format ('[TSATD error 2] ', c & 'density less than triple point density; d =', c & g11.5,' mol/L, dtp =',g11.5,' mol/L.',a1) c RETURN c endif if (d.gt.dtp(icomp)) then ierr=141 write (herr,1149) d,dtp(icomp),hnull call ERRMSG (ierr,herr) 1149 format ('[TSATD error 141] ', & 'density greater than triple point density; d=', & g11.5,' mol/L, dtp =',g11.5,' mol/L.',a1) RETURN endif if (d.ge.dcrit(icomp)) i=1 if (d.lt.dcrit(icomp)) i=2 endif c call SOLVEA (i,d,x,t,ierr,herr) c RETURN end !subroutine TSATD c c ====================================================================== c subroutine SOLVEA (iflag,pd,x,t,ierr,herr) c c solve for temperature given vapor pressure or saturated liquid c or vapor density c c inputs: c iflag--use: vapor pressure equation when iflag=0 c saturated liquid density equation when iflag=1 c saturated vapor density equation when iflag=2 c liquid pressure equation when iflag=4 c pd--pressure [kPa] or density [mol/L] c x--composition [array of mol frac] c output: c t--temperature (K) c ierr--error flag: 0 = successful c 124 = no convergence c herr--error string (character*255 variable if ierr<>0) c c written by E.W. Lemmon, NIST Physical & Chemical Properties Division, Boulder, Colorado c 03-13-01 EWL, original version c implicit double precision (a-h,o-z) implicit integer (i-n) 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*255 herr common /NCOMP/ nc,ic common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(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) dimension x(ncmax) c tol=0.5d-8 dx=0.001d0 dt=0.d0 icomp=1 if (ic.ne.0) icomp=ic c use midpoint of triple and critical values as starting point ttrp=ttp(icomp) tc=tcrit(icomp) t=(ttp(icomp)+tcrit(icomp))/2.d0 c calculate vapor pressure (or saturated density) and its derivative, c then update temperature iteratively do i=1,50 dy=0 if (t+dx.ge.tc) dy=dx !Prevent numerical deriv. going above Tc if (iflag.eq.0) then call PSATT (t-dy ,x,pp,ierr,herr) call PSATT (t-dy+dx,x,p1,ierr,herr) elseif (iflag.eq.1) then call DLSATT (t-dy ,x,pp,ierr,herr) call DLSATT (t-dy+dx,x,p1,ierr,herr) elseif (iflag.eq.2) then call DVSATT (t-dy ,x,pp,ierr,herr) call DVSATT (t-dy+dx,x,p1,ierr,herr) elseif (iflag.eq.4) then call PLSATT (t-dy ,x,pp,ierr,herr) call PLSATT (t-dy+dx,x,p1,ierr,herr) endif if (ierr.ne.0) return told=t if (ABS(pd).lt.1.d-15 .or. ABS(p1).lt.1.d-15) then t=0.5d0*(t-dt) else if (iflag.eq.1) then dpt=(p1-pp)/dx dt=-(pp-pd)/dpt !1st order Newton's method else dpt=ABS(log(p1)-log(pp))/dx dt=-(log(pp)-log(pd))/dpt !1st order Newton's method endif t=t+dt endif if (ABS(dt).lt.tol) RETURN if (t.lt.ttrp) t=(told+ttrp)/2.d0 if (t.ge.tc) t=(told+tc)/2.d0 enddo ierr=124 t=-9.99999d6 write (herr,1100) call ERRMSG (ierr,herr) 1100 format('[SOLVEA error 124] maximum number of iterations exceeded') RETURN end !subroutine SOLVEA c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_ANC.f c ======================================================================