c begin file setup.f c c This file contains the routines which initialize the models, fluid- c specific parameters and coefficients, etc. The subroutine SETUP must c be called before any of the other property routines are called. c Call(s) to SETMOD and SETREF are optional and may be used to specify c non-standard models and reference states. c c contained here are: c subroutine SETUP (nc,hfiles,hfmix,hrf,ierr,herr) c subroutine SETFLD (icomp,hfile,ierr,herr) c subroutine SETREF (hrf,ixflag,x0,h0,s0,t0,p0,ierr,herr) c subroutine SETMOD (nc,htype,hmix,hcomp,ierr,herr) c subroutine SETKTV (icomp,jcomp,hmodij,fij,hfmix,ierr,herr) c subroutine GETKTV (icomp,jcomp,hmodij,fij,hfmix,hfij,hbinp,hmxrul) c subroutine RFFILE (hfilei,hflref) c block data BDSET c c these routines set the values in the following common blocks c common /NCOMP/ ncomp c common /CCAS/ hcas(n0:nx) c common /CNAM/ hname(n0:nx) c common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) c common /TRNMOD/ heta,hetak(nrf0:nx),htcx,htcxk(nrf0:nx) c common /STNMOD/ hsten,hstenk(n0:nx) c common /REFST/ hrefst,hrefdf(n0:nx) c common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) c common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), c & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), c & accen(n0:nx),dipole(n0:nx) c c various arrays are dimensioned with parameter statements c parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) c parameter (nrf0=0) !lower limit for transport ref fluid arrays c c ====================================================================== c ====================================================================== c subroutine SETUP (nc,hfiles,hfmix,hrf,ierr,herr) c c define models and initialize arrays c c A call to this routine is required. c c inputs: c nc--number of components (1 for pure fluid) [integer] c hfiles--array of file names specifying fluid/mixture components c [character*80 variable] for each of the nc components; c e.g., :fluids:R134a.fld (Mac) or fluids\R134a.fld (DOS) or c [full_path]/fluids/R134a.fld (UNIX) c hfmix--mixture coefficients [character*80] c file name containing coefficients for mixture model, c if applicable c e.g., :fluids:HMX.bnc c hrf--reference state for thermodynamic calculations [character*3] c 'DEF': default reference state as specified in fluid file c is applied to each pure component c 'NBP': h,s = 0 at pure component normal boiling point(s) c 'ASH': h,s = 0 for sat liquid at -40 C (ASHRAE convention) c 'IIR': h = 200, s = 1.0 for sat liq at 0 C (IIR convention) c other choices are possible, but these require a separate c call to SETREF c outputs: c ierr--error flag: 0 = successful c 101 = error in opening file c 102 = error in file or premature end of file c -103 = unknown model encountered in file c 104 = error in setup of model c 105 = specified model not found c 111 = error in opening mixture file c 112 = mixture file of wrong type c herr--error string (character*255 variable if ierr<>0) c [fluid parameters, etc. returned via various common blocks] c c explanation of parameters (used to dimension arrays) c ncmax: maximum number of mixture components c n0: lower bound on component arrays (-ncmax to accomodate c multiple ECS-thermo reference fluids) c nx: same as ncmax c nrf0: lower bound on arrays associated with transport props c (element 0 stores reference fluid information) c c explanation of commons c /NCOMP/: nc--number of components c c /EOSMOD/ equation of state (thermodynamic) models c hpheq--model for phase equilibria calcs (not currently used) c heos--equation of state model c hmxeos--models for mixture components c hmodcp--ideal heat capacity model to use for component i c (set by call to SETeos routine) c c /TRNMOD/ transport property models c heta--viscosity model for mixture c hetak--viscosity models for mixture components c htcx--thermal conductivity model for mixture c htcxk--conductivity models for mixture components c c /STNMOD/ surface tension models c hsten--surface tension model for mixture c hstenk--surface tension models for mixture components c c /REFST/ c hrefst--reference state (as specified by input argument hrf) c hrefdf--default reference state for each component c c /CCON/: constants for each of the NCMAX components; these values c are taken from the coefficients/array of selected model c wm(i): molecular mass (g/mol) c ttp(i): triple point temperature (K) c tnbp(i): normal boiling point temperature (K) c tc(i): critical temperature (K) c pc(i): critical pressure (kPa) c rhoc(i): critical density (mol/L) c Zcrit(i): critical compressibility factor (Pc/(R*Tc*rhoc)) c accen(i): accentric factor c dipole(i): dipole moment [debye] (eval at Tnbp if t-dependent) c c /CREF/: reference state for each of the NCMAX components c tref(i): reference temperature for enthalpy and entropy c rhoref(i): reference density for enthalpy and entropy c href(i): enthalpy at tref(i), rhoref(i) for component i c sref(i): entropy at tref(i), rhoref(i) for component i c c /CCAS/: CAS numbers for the n0..ncmax components c hcas(i): "i" is mixture component number (n0..ncmax) c c /CNAM/: short-hand names for the nc components c hname(i): "i" is mixture component number c c /CHAR/: characters used to delimit/terminate output c htab: tab (or other character) to delimit output tables c hnull: ASCII null character to terminate error strings c (for compatibility with mixed-language DLLs) c c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-19-94 MM, original version c 07-21-95 MM, restructure to add subordinate (model-specific) c set-up routines c 07-24-95 MM, fluid const (e.g. Tc) set in SETeos rather than here c 09-10-95 MM, add common /CHAR/ and define contents c 09-13-95 MM, add ierr, herr to argument list c 09-25-95 MM, new argument list to SATT (outputs in order p, rho, x) c 10-03-95 MM, read coefficients from files, change argument list c 11-02-95 MM, expand to include mixture Helmholtz (HMX) model c add common block CNAM with fluid names c 11-08-95 MM, initialize hmxeos array, even if heos not HMX c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c kludge to load R134a as ECS reference fluid c 12-08-95 MM, move file read to separate subroutine SETFLD c 12-12-95 MM, remove R134a reference fluid kludge c 01-09-96 MM, add call to SETHMX (read HMX.bnc file) c 01-10-96 MM, add check for inputs same as previous call c 01-11-96 MM, move reference state calculation to SETREF c 01-19-96 MM, reset flag ksetrf to zero on call to SETUP c 01-31-96 MM, implement 'NBS' option for EOS (and transport model) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 02-29-96 MM, cover case of pure-fluid model specified for mixture c 03-13-96 MM, changes to merge with Klein's transport code, c add common /TRNMOD/ c 03-19-96 MM, add dipole moment to /CCON/ c 03-20-96 MM, reduce argument list, move model specification to SETMOD c add mix models to /TRNMOD/ c 03-21-96 MM, replace /MODEL/ with /EOSMOD/, /STNMOD/, /REFST/ c 03-27-96 MM, fix bug when nc = 1 but overall model is a mix model c 05-08-96 MM, add /MXINFO/ load hfmix into hmfile(0) c 05-14-96 MM, add call to SETPH0 model (Helmholtz form) c disable check if inputs same as previous (!temporary) c 11-04-96 MM, change nmxpar from 4 to 6 c 11-19-96 MM, add commons related to critical lines, print parameters c 11-25-96 MM, move debug print of mixture pars to SETHMX c 02-20-97 MM, add hrefdf to /REFST/, new common /CREFDF/ (for default ref st) c 02-24-97 MM, add /CREMOD/ (pointer to transport critical enhancement models) c 03-25-97 MM, bug fix: x0 used but not dimensioned c 03-28-97 MM, bug: models initialized to 'NBS' only for first call to SETUP, c subsequent calls find past settings or 'NUL' if ncomp has increased; c add flag in /RESETM/ to fix c 05-27-97 MM, reset gas constant when lreset=.true. c 07-03-97 MM, use 'DEF' reference state if 'OTH' is specified c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: SETUP c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters character*1 htab,hnull character*3 hrf,hrefst,hrefdf character*3 hpheq,heos,hmxeos,hmodcp character*3 heta,hetak,htcx,htcxk character*3 hsten,hstenk character*3 hph,heqn,heqnk,hvis,hvisk,hcond,hcondk,hsurf,hsurfk character*3 hetacr,htcxcr character*12 hcas,hname character*80 hfiles(ncmax),hfile(n0:nx),hfmix character*80 hfsav,hmxsav character*255 herr,herr1 c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin logical lreset dimension x0(ncmax) c common /NCOMP/ ncomp common /Gcnst/ R common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /TRNMOD/ heta,hetak(nrf0:nx),htcx,htcxk(nrf0:nx) common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax) common /STNMOD/ hsten,hstenk(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /REFST/ hrefst,hrefdf(n0:nx) common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /CREFDF/ tdef(n0:nx),pdef(n0:nx),hdef(n0:nx),sdef(n0:nx) common /CCAS/ hcas(n0:nx) common /CNAM/ hname(n0:nx) common /CHAR/ htab,hnull common /SETSAV/ hfsav(n0:nx),hmxsav common /IRFSAV/ ixfsav,ksetrf common /MODSAV/ hph,heqn,heqnk(nx),hvis,hvisk(nx), & hcond,hcondk(nx),hsurf,hsurfk(nx) c /MXINFO/ contains information on the mixing rules and parameters c hmfile specifies the files from which mixing rules originate c (stored in order 1,2 1,3 2,3; element zero is file called in SETUP) c hbin provides documentation for the current binary parameters c hrule contains descriptions of the currently available mixing rules c hbpar contains descriptions of the binary parameters (e.g. Kt, Kv) c associated with the currently available mixing rules common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) common /MXRULE/ nrule common /RESETM/ lreset !flag indicating need to reset all models c htab=CHAR(9) !tab character for output tables hnull=CHAR(0) !null character to terminate error strings ierr=0 herr=hnull * write (*,1002) heos,(hmxeos(j),j=n0,ncmax), * & heta,(hetak(j),j=nrf0,nx), * & htcx,(htcxk(j),j=nrf0,nx), * & hsten,(hstenk(j),j=n0,ncmax) *1002 format (' SETUP--heos, hmxeos: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)/ * & ' (input) visc models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' tcx models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' surf ten models: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)) * write (*,1003) heqn,(heqnk(j),j=1,ncmax), * & hvis,(hvisk(j),j=1,nx), * & hcond,(hcondk(j),j=1,nx), * & hsurf,(hsurfk(j),j=1,ncmax) *1003 format (/ * & ' SETUP--heqn, heqnk : ',a3,31x,5(a3,1x)/ * & ' (saved) visc models: ',a3,31x,5(a3,1x)/ * & ' tcx models: ',a3,31x,5(a3,1x)/ * & ' surf ten models: ',a3,31x,5(a3,1x)) c c compare inputs to previous values (stored in /xxxMOD/ and /SETSAV/) c if nothing has changed SETUP can be bypassed c lsame=.false. * write (*,*) ' SETUP: starting value of lsame = ',lsame if (nc.eq.ncomp .and. hrf.eq.hrefst) then * write (*,*) ' SETUP: ncomp and ref state same as previous' if ((heqn.eq.heos .or. heqn.eq.'NBS') .and. !EOS model & (hvis.eq.heta .or. hvis.eq.'NBS') .and. !viscosity model & (hcond.eq.htcx .or. hcond.eq.'NBS') .and. !therm cond model & (hsurf.eq.hsten .or. hsurf.eq.'NBS') .and. !surface ten model & hph.eq.hpheq) then !phase equil method, not implemented lsame=.true. * write (*,*) ' SETUP: after check of models lsame = ',lsame do 98 i=1,nc c check that .fld file and component models are the same as last call if (hfiles(i).ne.hfsav(i) .or. & (heqnk(i).ne.hmxeos(i) .and. heqnk(i).ne.'NBS') .or. c different models for surface tension not implemented c & (hsurfk(i).ne.hstenk(i) .and. hsurfk(i).ne.'NBS') .or. & (hvisk(i).ne.hetak(i) .and. hvisk(i).ne.'NBS') .or. & (hcondk(i).ne.htcxk(i) .and. hcondk(i).ne.'NBS')) then lsame=.false. end if * write (*,*) ' SETUP: after component models lsame = ',lsame 98 continue if (hmxsav.ne.hfmix) lsame=.false. !mixture coeff file * write (*,*) ' SETUP: after mixture file lsame = ',lsame end if c c !temporary--disable check of inputs vs. previous inputs lsame=.false. c if (lsame) then * write (*,*) * write (*,*) ' SETUP--inputs are identical, return' * write (*,*) RETURN end if end if c c check state of flag lreset; a value of .true. indicates need to reset c all models to 'NBS'; this occurs when SETUP is called a second time c in a given application if (lreset) then * write (*,*) ' SETUP--lreset=.true., call SETMOD to reset models' c gas constant, CODATA recommended value R=8.314510d0 c gas constant, Moldover's value c R=8.314471d0 call SETMOD (nc,'NBS',heos,hmxeos,ierr,herr) * write (*,1004) heos,(hmxeos(j),j=n0,ncmax), * & heta,(hetak(j),j=nrf0,nx), * & htcx,(htcxk(j),j=nrf0,nx), * & hsten,(hstenk(j),j=n0,ncmax) *1004 format (' SETUP--heos, hmxeos: ',2(a3,2x,5(1x,a3),3x)/ * & ' (reset) visc models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' tcx models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' surf ten models: ',2(a3,2x,5(1x,a3),3x)) end if c c set pointers to models c c special case if equation of state is NIST recommendation c only mixture choice at present is the HMX model if (heos.eq.'NBS' .or. heos.eq.'nbs') then if (nc.eq.1) then heos='NBS' else heos='HMX' end if do 122 i=1,nc hmxeos(i)='NBS' 122 continue else if (heos.eq.'HMX') then c mixture model is HMX, use component models as set in SETMOD else if (nc.gt.1) then c this condition should not be accessed; a mix model other than HMX c has been specified, but only choice at present is the HMX model ierr=-105 herr='[SETUP warning -105] unknown mixture model specified,' & //' the HMX model will be used.'//hnull call ERRMSG (ierr,herr) heos='HMX' do 128 i=1,nc hmxeos(i)='NBS' 128 continue end if end if c c special case if viscosity model is NIST recommendation c only mixture choice at present is the ECS model if (heta.eq.'NBS' .or. heta.eq.'nbs') then if (nc.eq.1) then heta='NBS' else heta='ECS' end if do 142 i=1,nc hetak(i)='NBS' 142 continue else if (heta.eq.'ECS') then c mixture model is ECS, use component models as set in SETMOD else if (nc.gt.1) then c this condition should not be accessed; a mix model other than ECS c has been specified, but only choice at present is the ECS model * write (*,*) ' SETUP--unknown mix viscosity model, ECS used' heta='ECS' do 148 i=1,nc hetak(i)='NBS' 148 continue end if end if c c special case if thermal conductivity model is NIST recommendation c only mixture choice at present is the ECS model if (htcx.eq.'NBS' .or. htcx.eq.'nbs') then if (nc.eq.1) then htcx='NBS' else htcx='ECS' end if do 162 i=1,nc htcxk(i)='NBS' 162 continue else if (htcx.eq.'ECS') then c mixture model is ECS, use component models as set in SETMOD else if (nc.gt.1) then c this condition should not be accessed; a mix model other than ECS c has been specified, but only choice at present is the ECS model * write (*,*) ' SETUP--unknown mix conductivity model, ECS used' htcx='ECS' do 168 i=1,nc htcxk(i)='NBS' 168 continue end if end if c c special case if surface tension model is NIST recommendation if (hsten.eq.'NBS' .or. hsten.eq.'nbs') then if (nc.eq.1) then hsten='NBS' else hsten='STM' end if do 172 i=1,nc hstenk(i)='NBS' 172 continue c else if (hsten.eq.'STM' .or. hsten.eq.'STX') then c mixture model is STM or STX, both are contained in routine STN; c use component models as specified in SETMOD c else c if (nc.gt.1) then c this condition should not be accessed; a mix model other than STM or STX c has been specified, but no other choices are available, at present c write (*,*) ' SETUP--unknown surface tension model, STM used' c hsten='STM' c do 178 i=1,nc c hstenk(i)='NBS' c 178 continue c end if end if c c reset flag ksetrf and copy inputs to common c ksetrf=0 do 182 i=1,nc hfile(i)=hfiles(i) hfsav(i)=hfiles(i) 182 continue hph='NBS' hmxsav=hfmix ncomp=nc c c fill up info arrays for unused/undefined components/reference fluid c do 184 i=n0,0 hfile(i)='NUL' hmxeos(i)='NUL' hcas(i)='not_defined' hname(i)='not_defined' hmxeos(i)='NUL' hstenk(i)='NUL' 184 continue do 186 i=nrf0,0 hetak(0)='NUL' htcxk(0)='NUL' 186 continue if (nc.lt.ncmax) then do 188 i=nc+1,ncmax hcas(i)='not_defined' hname(i)='not_defined' hmxeos(i)='NUL' hetak(i)='NUL' htcxk(i)='NUL' hstenk(i)='NUL' 188 continue end if c c store fluid constants for each of the nc components c do 200 icomp=1,nc call SETFLD (icomp,hfile,ierr,herr) if (ierr.eq.101) then c error in opening file--further processing pointless RETURN end if 200 continue c c store fluid constants for ECS reference fluids (if any) c do 260 icomp=n0,0 if (hfile(icomp).ne.'NUL') then call SETFLD (icomp,hfile,ierr,herr) if (ierr.eq.101) then c error in opening file--further processing pointless RETURN end if end if 260 continue c c set up mixture model, if applicable c if (heos.eq.'HMX') then c store file containing default mixture parameters hmfile(0)=hfmix * write (*,*) ' SETUP--about to call SETHMX' call SETHMX (hfmix,ierr1,herr1) c write (*,*) 'ierr from SETHMX: ',ierr1,': ',herr1 call ERRMSG (ierr1,herr1) end if c c set reference state c ixflag=1 !only 'pure fluid' reference state can be called here if (hrf.eq.'OTH') then ierr=-105 herr='[SETUP warning -105] must use routine SETREF for (OTH) ' & //'reference state choice; will use default (DEF) choice' & //hnull call ERRMSG (ierr,herr) c use the default reference state hrf='DEF' c fill up reference arrays so that calculations do not crash c do 320 i=n0,ncomp c if (tc(i).gt.233.15) then c tref(i)=233.15d0 c else c tref(i)=0.90*tc(i) c end if c sref(i)=0.0d0 c href(i)=0.0d0 c rhoref(i)=1.0d0 c write (*,1005) i,tref(i),href(i),sref(i),rhoref(i) c1005 format (1x,' SETUP--T,h,s,D-ref: ',i3,f7.2,f10.1,2f10.3) c 320 continue c else end if call SETREF (hrf,ixflag,x0,h0,s0,t0,p0,ierr,herr) c write (*,*) ' SETUP--ierr returned from SETREF: ',ierr call ERRMSG (ierr,herr) c any error from SETREF takes precedence over error from SETHMX if (ierr.eq.0) then ierr=ierr1 herr=herr1 end if c c save current models into /MODSAV/--check on subsequent calls to c SETUP, if nothing has changed, can bypass call c hph=hpheq heqn=heos hvis=heta hcond=htcx hsurf=hsten do 440 k=1,nx heqnk(k)=hmxeos(k) hvisk(k)=hetak(k) hcondk(k)=htcxk(k) hsurfk(k)=hstenk(k) 440 continue c * write (*,1442) heos,(hmxeos(j),j=n0,ncmax), * & heta,(hetak(j),j=nrf0,nx), * & (hetacr(j),j=nrf0,ncomp) * write (*,1443) htcx,(htcxk(j),j=nrf0,nx), * & (htcxcr(j),j=nrf0,ncomp) * write (*,1444) hsten,(hstenk(j),j=n0,ncmax) *1442 format (/ * & ' SETUP--heos, hmxeos: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)/ * & ' (exit) visc models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' visc crit mods: ',3x,25x,a3,3x,5(a3,1x)) *1443 format (' tcx models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' t.c. crit mods: ',3x,25x,a3,3x,5(a3,1x)) *1444 format (' surf ten models: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)) c c set flag to reset all models to 'NBS' on subsequent call to SETUP lreset=.true. * write (*,*) 'ierr at SETUP exit: ',ierr,': ',herr c RETURN c end !subroutine SETUP c c ====================================================================== c subroutine SETFLD (icomp,hfile,ierr,herr) c c open a fluid file and read model coefficients (or get from block data) c c inputs: c icomp--pointer specifying component number c zero and negative values are used for ECS reference fluid(s) c hfile--array of file names specifying fluid/mixture components c [character*80 variable] for each of the components; c --or-- c when hf(i) is of the form: c BDATA:nn-nn-nn c use coefficients stored in block data for fluid with CAS c number specified by nn-nn-nn; c e.g. to use stored formulation for R134a, c hf(i) = 'BDATA:811-97-2' c outputs: c ierr--error flag: 0 = successful c 101 = error in opening file c 102 = error in file or premature end of file c -103 = unknown model encountered in file c 104 = error in setup of model c 105 = specified model not found c herr--error string (character*255 variable if ierr<>0) c [fluid parameters, etc. returned via various common blocks] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 12-08-95 MM, original version--extracted from old subroutine SETUP c 12-12-95 MM, add calls to SETECS, SETETA, SETTCX c 12-14-95 MM, add call to RFFILE (add directory to ref fluid file) c 01-12-96 MM, do not reset hmxeos(0) if in file 0 (reference fluid) c 01-23-96 MM, always call SETETA, SETTCX, etc. if found in .fld file c 01-31-96 MM, implement 'NBS' option for EOS (and transport model) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 02-29-96 MM, check only hmxeos(i) [not heos] for call to SETxxx c 03-13-96 MM, changes to merge with Klein's transport code c add common /TRNMOD/; add calls to SETTRN, SETVSi, SETTCi c 03-15-96 MM, always load ECS-transport model (unless htran = 'NUL') c 03-19-96 MM, read dipole moment and add to /CCON/ c 03-20-96 MM, add mixture models to /TRNMOD/ c 03-21-96 MM, replace /MODEL/ with /EOSMOD/, /STNMOD/ c 03-27-96 MM, add calls to set up surface tension c 03-28-96 MM, add flags ltcx,leta to cover case of fluid-specific c model available for only one of the transport props c 06-17-96 MM, check only 'CP' rather than 'CPP' to allow CP1 c 06-18-96 MM, add calls to SETVS4-6 and SETTC2-4 for future use c 10-16-96 MM, add call to SETCI2 (collision integral AUX function) c 10-30-96 MM, add calls to SETTC5, SETTC6 c 11-06-96 MM, check that i.ge.nrf0 before calling collision integral c 01-21-97 MM, add call to SETCI1 (collision integral AUX function) c 02-06-97 MM, add /CNAM80/ to store full chemical name c 02-20-97 MM, read in default reference state from fluid file, c add hrefdf to /REFST/, new common /CREFDF/ c 02-24-97 MM, add calls to transport critical models and /CREMOD/ c 03-26-97 EWL, minor changes for Fortran 90 compatibility c 03-27-97 EWL, if file not found, search in likely directories c 08-22-97 MM, use double backslash on file open to avoid problem with Unix machines c 10-28-97 MM, put fluid file version no. into common c 12-01-97 MM, read in synonyms, add to /CNAM80/ c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*1 h1,hstar character*3 heqn,heta,htcx,hetak,htcxk character*3 hetemp,httemp character*3 hpheq,heos,hmxeos,hmodcp character*3 hsten,hstenk character*3 hmodci,hmaux character*3 hvs,htc character*3 hflag,htype character*3 hrefst,hrefdf character*3 hetacr,htcxcr character*12 hcasn(n0:nx),hcas,hname character*80 hfile(n0:nx),href,hnam80,hfilei,hsyn1,hsyn2 character*255 herr,herr1 logical lbdata,leta,ltcx c common /NCOMP/ ncomp common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /TRNMOD/ heta,hetak(nrf0:nx),htcx,htcxk(nrf0:nx) common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax) common /STNMOD/ hsten,hstenk(n0:nx) common /OMGMOD/ hmodci(nrf0:nx) common /REFST/ hrefst,hrefdf(n0:nx) common /CREFDF/ tdef(n0:nx),pdef(n0:nx),hdef(n0:nx),sdef(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /CCAS/ hcas(n0:nx) common /CNAM/ hname(n0:nx) common /CNAM80/ hnam80(n0:nx),hsyn1(n0:nx),hsyn2(n0:nx) common /CHAR/ htab,hnull common /VERS/ verfl(n0:nx),vermx !fluid & mix file version nos. c i=icomp ltcx=.false. leta=.false. * write (*,1102) i,hfile(i) *1102 format (/1x,'input file for fluid ',i3,': (',a40,')') * write (*,1003) heos,(hmxeos(j),j=n0,ncmax), * & heta,(hetak(j),j=nrf0,nx),htcx,(htcxk(j),j=nrf0,nx) * & ,hsten,(hstenk(j),j=n0,ncmax) *1003 format (' SETFLD--heos, hmxeos: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)/ * & ' (input) visc models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' tcx models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' surf ten models: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)) if (hfile(i)(1:5).eq.'BDATA' .or. hfile(i)(1:5).eq.'bdata') then c get coefficients from common blocks lbdata=.true. nread=0 hcasn(i)=hfile(i)(7:18) hcas(i)=hcasn(i) * write (*,1104) hcasn(i) *1104 format (1x,'coefficients from block data; CAS no.: (',a12,')') else c read coefficients from file c write (*,1105) hfile(i) c1105 format (1x,'reading coefficients from file: (',a40,')') lbdata=.false. nread=12 !logical unit for file reads c open (unit=nread,file=hfile(i),err=998) open (unit=nread,file=hfile(i),status='old',iostat=ioerr) if (ioerr.ne.0) then c if error in opening file, try looking in likely directory c use double backslash on file open to avoid problem with Unix machines hfilei='\refprop\fluids\\'//hfile(i) open (unit=nread,file=hfilei,status='old',iostat=ioerr) if (ioerr.ne.0) then hfilei='\refprop6\fluids\\'//hfile(i) open (unit=nread,file=hfilei,status='old',iostat=ioerr) if (ioerr.ne.0) then hfilei='MPW:code:Refprop613:'//hfile(i) open (unit=nread,file=hfilei,status='old',err=998) if (ioerr.ne.0) then hfilei='MPW:code:Refprop613:fluids:'//hfile(i) open (unit=nread,file=hfilei,status='old',err=998) end if end if end if * write (*,1106) hfile(i),hfilei *1106 format (1x,'initial file specification failed: ',a80/ * & 1x,'reading coefficients from file: ',a80) end if read (nread,2012) hname(i) !short name read (nread,2012) hcasn(i) !CAS number hcas(i)=hcasn(i) read (nread,2080) hnam80(i) !full chemical name c write (*,*) ' SETUP--full name: ',hnam80(i) read (nread,2080) hsyn1(i) !synonym 1 read (nread,2080) hsyn2(i) !synonym 2 c note: read in generic fluid constants; these may be reset by c some SETeos routines to correspond with the values used c in the respective models read (nread,*) wm(i) !molecular weight [g/mol] read (nread,*) ttp(i) !triple point temperature [K] read (nread,*) tnbp(i) !normal boiling point [K] read (nread,*) tc(i) !critical temperature [K] read (nread,*) pc(i) !critical pressure [kPa] read (nread,*) rhoc(i) !critical density [mol/L] read (nread,*) accen(i) !accentric factor [-] read (nread,*) dipole(i) !dipole moment [debye] read (nread,2003) hrefdf(i) !default reference state if (hrefdf(i).eq.'OTH') then c for 'OTHer' reference state read in reference T,p,h,s read (nread,*) tdef(i),pdef(i),hdef(i),sdef(i) end if read (nread,*) verfl(i) !version number * write (*,*) ' SETUP--end of constants; default = ',hrefdf(i) end if c c search for key characters in cols 1 and 2-4 (if reading from file) c and branch to subsidiary setup routines for specified models c 160 if (.not.lbdata) then !read from file read (nread,2013,end=199,err=199) hstar,hflag c write (*,*) ' SETUP--input line: ',hstar,hflag c hstar = '#' or '@' indicates start of model specification c # indicates NIST-recommended model c @ indicates other model(s) c hflag indicates type of model: c 'EOS' = equation of state c 'TRN' = transport property model (i.e. ECS model) c 'ETA' = pure fluid viscosity model c 'TCX' = pure fluid thermal conductivity model c 'AUX' = auxiliary model, such as ideal heat capacity if (hstar.eq.'#' .or. hstar.eq.'@') then read (nread,2003) htype c htype is 3-letter key for particular models, for example: c 'NBS' = NIST-recommended model c 'BWR' = modified Bennedict-Webb-Rubin equation of state c 'FEQ' = fundamental (Helmholtz) equation of state c 'ECS' = extended corresponding states model c 'HMX' = mixture Helmholtz model c 'CPP' = polynomial fit of ideal gas heat capacity c 'VSi' = pure fluid viscosity model #i c 'TCi' = pure fluid thermal conductivity model #i c write (*,*) ' SETUP--hstar,hflag: ',hstar,hflag c write (*,*) ' SETUP--htype: ',htype c do 120 k=1,1000 !skip over source comments read (nread,2001) h1 if (h1.ne.'?') goto 124 !comment block terminated by '!' 120 continue 124 continue end if end if c if (lbdata) then hflag='EOS' htype=heos end if c if (hflag.eq.'EOS' .or. hflag.eq.'eos') then c c set up equation of state c if (hstar.eq.'#') then c special case for NIST-recommended model if (ncomp.eq.1 .and. heos.eq.'NBS') then heos=htype end if if (i.ge.1 .and. hmxeos(i).eq.'NBS') then hmxeos(i)=htype end if end if c if (htype.eq.'BWR') then c modified Bennedict-Webb-Rubin equation of state found in file c if (heos.eq.'BWR' .or. hmxeos(i).eq.'BWR') then if (hmxeos(i).eq.'BWR') then * write (*,*) ' SETUP--about to call SETBWR' call SETBWR (nread,i,hcasn(i),ierr1,herr1) end if else if (htype.eq.'FEQ') then c fundamental (Helmholtz) equation of state found in file c if (heos.eq.'FEQ' .or. hmxeos(i).eq.'FEQ') then if (hmxeos(i).eq.'FEQ') then * write (*,*) ' SETUP--about to call SETFEQ' call SETFEQ (nread,i,hcasn(i),ierr1,herr1) end if else if (htype.eq.'ECS') then c extended corresponding states (ECS) thermo model found in file c if (heos.eq.'ECS' .or. hmxeos(i).eq.'ECS') then if (hmxeos(i).eq.'ECS') then * write (*,*) ' SETUP--about to call SETECS (thermo)' call SETECS (nread,i,hcasn(i),href,heqn,ierr1,herr1) c add directory information to reference fluid file call RFFILE (hfile(i),href) iref=-i hfile(iref)=href hmxeos(iref)=heqn end if else c unidentified model found in file ierr=-103 write (herr,2103) htype,i,hnull call ERRMSG (ierr,herr) end if if (ierr1.ne.0) then ierr=104 write (herr,2104) hmxeos(i),i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if end if c if (lbdata) then hflag='AUX' htype=hmodcp(i) end if c if (hflag.eq.'AUX') then c c set up auxiliary model(s) c (model(s) are specified in call to appropriate SETmod routine) c * write (*,*) ' SETUP--aux model (',htype,') found in fld file ' c write (*,*) ' SETUP--i,htype,hmodcp(i): ',i,htype,hmodcp(i) hmaux=hmodcp(i) !possible use in error message if (htype(1:2).eq.'CP' .and. hmodcp(i).eq.htype) then c ideal gas heat capacity function * write (*,*) ' SETUP--about to call SETCPP' call SETCPP (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'PH0' .and. hmodcp(i).eq.'PH0') then c Helmholtz for ideal gas heat state * write (*,*) ' SETUP--about to call SETPH0' call SETPH0 (nread,i,hcasn(i),ierr1,herr1) else if (i.ge.nrf0) then c do not call collision integral model (or other auxiliary models c associated with the transport properties) if component number < 0 c (there is only one reference fluid for transport props so arrays c are not dimensioned for the negative component numbers set up to c accomodate the multiple reference fluids in the ECS-thermo model) if (htype(1:2).eq.'CI') hmaux=hmodci(i) !possible use in if (htype(1:2).eq.'TK') hmaux=htcxcr(i) !error message if (htype.eq.'CI1' .and. hmodci(i).eq.htype) then c functions for collision integral (used in transport correlations) * write (*,*) ' SETUP--about to call SETCI1' call SETCI1 (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'CI2' .and. hmodci(i).eq.htype) then * write (*,*) ' SETUP--about to call SETCI2' call SETCI2 (nread,i,hcasn(i),ierr1,herr1) c thermal conductivity critical enhancement models else if (htype.eq.'TK1' .and. htcxcr(i).eq.htype) then * write (*,*) ' SETUP--about to call SETTK1' call SETTK1 (nread,i,hcasn(i),ierr1,herr1) c else if (htype.eq.'TK2' .and. htcxcr(i).eq.htype) then c the TK2 model is linked with TC2 and is initialized directly by SETTC2 c call SETTK2 (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'TK3' .and. htcxcr(i).eq.htype) then * write (*,*) ' SETUP--about to call SETTK3' call SETTK3 (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'TK4' .and. htcxcr(i).eq.htype) then * write (*,*) ' SETUP--about to call SETTK4' call SETTK4 (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'TK5' .and. htcxcr(i).eq.htype) then * write (*,*) ' SETUP--about to call SETTK5' call SETTK5 (nread,i,hcasn(i),ierr1,herr1) else if (htype.eq.'TK6' .and. htcxcr(i).eq.htype) then * write (*,*) ' SETUP--about to call SETTK6' call SETTK6 (nread,i,hcasn(i),ierr1,herr1) end if else ierr1=104 end if if (ierr1.ne.0) then ierr=104 write (herr,2104) hmaux,i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if end if c if (lbdata) then hflag='TRN' end if c c transport property models c c do not call transport models if component number less than 0 c (there is only one reference fluid for transport props so arrays c are not dimensioned for the negative component numbers set up to c accomodate the multiple reference fluids in the ECS-thermo model) if (i.ge.nrf0) then c c the 'TRN' models apply to both the viscosity and thermal conductivity c if (hflag.eq.'TRN' .or. hflag.eq.'trn') then if (hstar.eq.'#') then c special case for NIST-recommended model if (ncomp.eq.1) then if (heta.eq.'NBS') heta=htype if (htcx.eq.'NBS') htcx=htype end if if (i.ge.1 .and. hetak(i).eq.'NBS') then hetak(i)=htype end if if (i.ge.1 .and. htcxk(i).eq.'NBS') then htcxk(i)=htype end if c in the case that the TRN modle is the NIST-recommended one, the ECS c method is equivalent to a fluid-specific model for purposes of the c leta and ltcx flags leta=.true. ltcx=.true. else c load the 'TRN' model if a fluid-specific 'ETA' or 'TCX' has not yet c been loaded (e.g. if a fluid-specific correlation is not available); c this is indicated by the flags leta and ltcx; if a fluid-specific c correlation is encountered later in the file this will be overwritten; c this is necessary to cover the case of a fluid-specific (and NIST-rec) c model available for only one of the transport properties if (.not.leta) hetemp=htype if (.not.ltcx) httemp=htype end if c if (htype.eq.'ECS') then c extended corresponding states model found in file c always load the ECS model (if any transport props are asked for) c for possible use if a pure fluid correlation is out of range if (heta.ne.'NUL' .or. htcx.ne.'NUL') then * write (*,*) ' SETUP--about to call SETTRN (ECS-transport)' call SETTRN (nread,i,hcasn(i),href,heqn,hvs,htc,ierr,herr) herr1=herr if (ierr.ne.0) then ierr=104 write (herr,2104) htype,i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if if (i.ge.1) then c store info for ECS reference fluid c (but not if current fluid is the reference fluid) c add directory information to reference fluid file call RFFILE (hfile(i),href) * write (*,1160) heqn,hvs,htc *1160 format (' SETFLD--transport ref fluid mods:',3(1x,a3)) hfile(0)=href !.fld file for reference fluid hmxeos(0)=heqn !model for ref fluid eqn of state hetak(0)=hvs !model for ref fluid viscosity htcxk(0)=htc !model for ref fluid conductivity end if end if else c unidentified model found in file ierr=-103 write (herr,2103) htype,i,hnull call ERRMSG (ierr,herr) end if end if !end of ECS-transport block c c pure fluid viscosity models c if (hflag.eq.'ETA' .or. hflag.eq.'eta') then if (hstar.eq.'#') then c special case for NIST-recommended model if (ncomp.eq.1 .and. heta.eq.'NBS') then heta=htype end if if (i.ge.1 .and. hetak(i).eq.'NBS') then hetak(i)=htype end if end if c if (htype.eq.'VS1') then c viscosity model "1" found in .fld file if (hetak(i).eq.'VS1') then * write (*,*) ' SETUP--about to call SETVS1' call SETVS1 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else if (htype.eq.'VS2') then c viscosity model "2" found in .fld file if (hetak(i).eq.'VS2') then * write (*,*) ' SETUP--about to call SETVS2' call SETVS2 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else if (htype.eq.'VS3') then c viscosity model "3" found in .fld file if (hetak(i).eq.'VS3') then * write (*,*) ' SETUP--about to call SETVS3' call SETVS3 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else if (htype.eq.'VS4') then c viscosity model "4" found in .fld file if (hetak(i).eq.'VS4') then * write (*,*) ' SETUP--about to call SETVS4' call SETVS4 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else if (htype.eq.'VS5') then c viscosity model "5" found in .fld file if (hetak(i).eq.'VS5') then * write (*,*) ' SETUP--about to call SETVS5' call SETVS5 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else if (htype.eq.'VS6') then c viscosity model "6" found in .fld file if (hetak(i).eq.'VS6') then * write (*,*) ' SETUP--about to call SETVS6' call SETVS6 (nread,i,hcasn(i),ierrvs,herr1) leta=.true. !flag that fluid-specific vis model loaded end if else c unidentified model found in file ierr=-103 write (herr,2103) htype,i,hnull call ERRMSG (ierr,herr) end if end if if (ierrvs.ne.0) then ierr=104 write (herr,2104) hetak(i),i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if c end of viscosity block c c pure fluid thermal conductivity models c if (hflag.eq.'TCX' .or. hflag.eq.'tcx') then if (hstar.eq.'#') then c special case for NIST-recommended model if (ncomp.eq.1 .and. htcx.eq.'NBS') then htcx=htype end if if (i.ge.1 .and. htcxk(i).eq.'NBS') then htcxk(i)=htype end if end if c if (htype.eq.'TC1') then c thermal conductivity model "1" found in .fld file if (htcxk(i).eq.'TC1') then * write (*,*) ' SETUP--about to call SETTC1' call SETTC1 (nread,i,hcasn(i),ierrtc,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else if (htype.eq.'TC2') then c thermal conductivity model "2" found in .fld file if (htcxk(i).eq.'TC2') then * write (*,*) ' SETUP--about to call SETTC2' call SETTC2 (nread,i,hcasn(i),ierrtc,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else if (htype.eq.'TC3') then c thermal conductivity model "3" found in .fld file if (htcxk(i).eq.'TC3') then * write (*,*) ' SETUP--about to call SETTC3' call SETTC3 (nread,i,hcasn(i),ierr,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else if (htype.eq.'TC4') then c thermal conductivity model "4" found in .fld file if (htcxk(i).eq.'TC4') then * write (*,*) ' SETUP--about to call SETTC4' call SETTC4 (nread,i,hcasn(i),ierrtc,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else if (htype.eq.'TC5') then c thermal conductivity model "5" found in .fld file if (htcxk(i).eq.'TC5') then * write (*,*) ' SETUP--about to call SETTC5' call SETTC5 (nread,i,hcasn(i),ierrtc,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else if (htype.eq.'TC6') then c thermal conductivity model "6" found in .fld file if (htcxk(i).eq.'TC6') then * write (*,*) ' SETUP--about to call SETTC6' call SETTC6 (nread,i,hcasn(i),ierrtc,herr1) ltcx=.true. !flag that fluid-specific t.c. model loaded end if else c unidentified model found in file ierr=-103 write (herr,2103) htype,i,hnull call ERRMSG (ierr,herr) end if end if if (ierrtc.ne.0) then ierr=104 write (herr,2104) htcxk(i),i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if c end of thermal conductivity block end if !end of overall transport properties block c c pure fluid surface tension models c c do not call surface tension models if component number less than 0 if (i.ge.nrf0) then if (hflag.eq.'STN' .or. hflag.eq.'stn') then if (hstar.eq.'#') then c special case for NIST-recommended model if (ncomp.eq.1 .and. hsten.eq.'NBS') then hsten=htype end if if (i.ge.1 .and. hstenk(i).eq.'NBS') then hstenk(i)=htype end if end if c if (htype.eq.'ST1') then c surface tension model "1" found in .fld file if (hstenk(i).eq.'ST1') then * write (*,*) ' SETUP--about to call SETST1' call SETST1 (nread,i,hcasn(i),ierrst,herr1) end if c else if (htype.eq.'ST2') then c surface tension model "2" found in .fld file c if (hstenk(i).eq.'TC2') then c write (*,*) ' SETUP--about to call SETST2' c call SETST2 (nread,i,hcasn(i),ierrst,herr1) c end if else c unidentified model found in file ierr=-103 write (herr,2103) htype,i,hnull call ERRMSG (ierr,herr) end if end if if (ierrst.ne.0) then ierr=104 write (herr,2104) hstenk(i),i,herr1(1:182),hnull call ERRMSG (ierr,herr) end if end if !end of surface tension block c if (lbdata) then RETURN !next component end if c if (hflag.eq.'END') then c if fluid-specific 'ETA' and/or 'TCX' models have not been loaded, use c the 'TRN' model (e.g. the ECS model) if (i.ge.nrf0) then if (.not.leta) hetak(i)=hetemp if (.not.ltcx) htcxk(i)=httemp end if c close the file, return to setup and process the next component rewind (nread) close (unit=nread,err=998) * write (*,1190) hfile(i) *1190 format (' SETUP--file closed: ',a80) RETURN else goto 160 end if 199 continue ierr=102 write (herr,2102) i,hnull 2102 format ('[SETUP error 102] error in file or premature end of', & ' file for component #',i3,'.',a1) call ERRMSG (ierr,herr) c rewind and close the file rewind (nread) close (unit=nread,err=998) * write (*,1199) hfile(i) *1199 format (1x,'SETFLD--premature end of file: ',a80) c if fluid-specific 'ETA' and/or 'TCX' models have not been loaded, use c the 'TRN' model (e.g. the ECS model) if (i.ge.nrf0) then if (.not.leta) hetak(i)=hetemp if (.not.ltcx) htcxk(i)=httemp end if RETURN c 998 ierr=101 write (herr,2101) i,hfile(i),hnull 2101 format ('[SETUP error 101] error in opening file for component #', & i3,'; filename = (',a80,').',a1) call ERRMSG (ierr,herr) RETURN c 2001 format (a1) 2003 format (a3) 2012 format (a12) 2013 format (a1,a3) 2080 format (a80) 2103 format ('[SETUP warning -103] unknown model (',a3, & ') encountered in file for component #',i3,'.',a1) 2104 format ('[SETUP error 104] error in setup of (',a3, & ') model for component #',i3,': ',a182,a1) c end !subroutine SETFLD c c ====================================================================== c subroutine SETREF (hrf,ixflag,x0,h0,s0,t0,p0,ierr,herr) c c set reference state enthalpy and entropy c c This subroutine must be called after SETUP; it need not be called at c all if the reference state specified in the call to SETUP is to be c used. c c inputs: c hrf--reference state for thermodynamic calculations [character*3] c 'NBP': h,s = 0 at normal boiling point(s) c 'ASH': h,s = 0 for sat liquid at -40 C (ASHRAE convention) c 'IIR': h = 200, s = 1.0 for sat liq at 0 C (IIR convention) c 'DEF': default reference state as specified in fluid file c is applied to each component (ixflag = 1 is used) c 'OTH': other, as specified by h0, s0, t0, p0 c ixflag--composition flag: 1 = ref state applied to pure components c 2 = ref state applied to mixture x0 c following input has meaning only if ixflag = 2 c x0--composition for which h0, s0 apply; array(1:nc) [mol frac] c this is useful for mixtures of a predefined composition, c e.g. refrigerant blends such as R410A c following inputs have meaning only if hrf = 'OTH' c h0--reference state enthalpy at t0,p0 {x0} [J/mol] c s0--reference state entropy at t0,p0 {x0} [J/mol-K] c t0--reference state temperature [K] c t0 = -1 indicates saturated liquid at normal boiling point c (bubble point for a mixture) c p0--reference state pressure [kPa] c p0 = -1 indicates saturated liquid at t0 {and x0} c p0 = -2 indicates saturated vapor at t0 {and x0} c outputs: c ierr--error flag: 0 = successful c 22 = Tmin > Tref for IIR reference state c 23 = Tcrit < Tref for IIR reference state c 24 = Tmin > Tref for ASHRAE reference state c 25 = Tcrit < Tref for ASHRAE reference state c 26 = Tmin > Tnbp for NBP reference state c 27 = Tref, Pref for OTH ref state outside limits c -28 = can't apply 'DEF' to mixture; c will apply to pure components c -29 = unknown reference state specified; c will use 'DEF' c herr--error string (character*255 variable if ierr<>0) c [fluid parameters, etc. returned via various common blocks] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 01-11-96 MM, original version--extracted from old subroutine SETUP c 01-16-96 MM, add check for inputs same as previous call c 01-19-96 MM, bug fixes and set flag ksetrf to 1 when SETREF called c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 03-19-96 MM, add dipole moment to /CCON/ c 03-21-96 MM, replace /MODEL/ with /REFST/ c 06-18-96 MM, check ref state temperature against LIMITX c 01-07-97 MM, check limits for case of NBP ref state (re: CO2 bug) c 01-10-97 MM, pass error flag from limits errors to calling routine c 02-20-97 MM, add 'DEF' option (default ref state read from fld file) c add hrefdf to /REFST/, new common /CREFDF/ c 03-25-97 MM, xliq,xvap used but not dimensioned; change to xl,xv c 03-26-97 EWL, changes for Fortran 90 compatibility c 05-14-97 MM, fix bug in checking for inputs same as last call c 10-01-97 MM, add compiler switches to allow access by DLL c 10-16-97 MM, allow ixflag = 2 for 'DEF' when all fluids are same c 11-24-97 MM, check for errors on calls to SATT, SATP, TPFLSH c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: SETREF c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) character*1 htab,hnull character*3 hrf,hjunk character*3 hrefst,hrefdf,hrefi(n0:nx) character*255 herr,herr2,herrsv dimension x(ncmax),x0(ncmax),xl(ncmax),xv(ncmax),prefi(n0:nx) c common /NCOMP/ ncomp common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /REFST/ hrefst,hrefdf(n0:nx) common /CHAR/ htab,hnull common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /CREFDF/ tdef(n0:nx),pdef(n0:nx),hdef(n0:nx),sdef(n0:nx) common /REFSAV/ x0sav(nx),h0sav,s0sav,t0sav,p0sav common /IRFSAV/ ixfsav,ksetrf c c initialize error code and string ierr=0 herr=hnull * write (*,1000) hrf,ixflag,ncomp *1000 format (/1x,' SETREF--entering reference state = ',a3, * & '; with ixflag, ncomp =',2i4) * write (*,1001) 'input',hrf,ixflag,h0,s0,t0,p0,(x0(i),i=1,ncomp) * write (*,1001) 'saved',hrefst,ixfsav,h0sav,s0sav,t0sav,p0sav, * & (x0sav(i),i=1,ncomp) *1001 format (' SETREF--',a5,' href,xflag,h0,s0,t0,p0: ',a3,i3,4f10.4/ * & ' ',5x,' x0: ',6x,5f10.4) c c compare inputs to previous values (stored in /REFST/ and /REFSAV/) c if nothing has changed SETREF can be bypassed c l1same=.false. !flag for initialization of pure components l2same=.false. !flag for application of ref state to mixture ixflg=ixflag !input value changed if input href not valid if (hrf.eq.hrefst) then l1same=.true. if (hrf.eq.'OTH') then c check h0, s0, t0, p0 only if 'OTHer' reference state is used if (abs(s0-s0sav).gt.1.0d-10 .or. & abs(h0-h0sav).gt.1.0d-10 .or. & abs(t0-t0sav).gt.1.0d-10 .or. & abs(p0-p0sav).gt.1.0d-10) then l1same=.false. end if end if if (ixflag.eq.ixfsav) then l2same=.true. if (ixflag.eq.2) then c check x0 only if ixflag = 2 do 100 i=1,ncomp if (abs(x0(i)-x0sav(i)).gt.1.0d-10) then l2same=.false. end if 100 continue end if end if if (l1same .and. l2same .and. ksetrf.eq.1) then c return if all inputs are same and also SETUP has not been executed c since last call to SETREF (as indicated by ksetrf) * write (*,*) * write (*,*) ' SETREF--inputs are identical, return' RETURN end if end if c * write (*,*) if (l1same .and. ksetrf.eq.1 .and. ixflag.ge.ixfsav) then c can skip initialization of pure components if: c --reference state is the same (as indicated by l1same) c --SETUP has not been called since last call to SETREF (indicated by ksetrf) c --value of ixflag is >= ixfsav, i.e. this call applies ref state to mixture c while last call applied to the pures, or if both applied to pures, c but not if last call was applied to mix, this call to pures * write (*,*) ' SETREF--skip initialization of pure components' else c must initialize pure component reference states do 200 i=n0,ncomp c following is temporary assignment, needed for SATP and/or THERM c (SATP calls ENTRO, and ENTRO needs rhoref, but SATP is used to c calculate rhoref; thus need to put in temporary value here) rhoref(i)=1.0d0 200 continue c if (hrf.eq.'DEF' .and. ixflag.eq.2) then c check if all fluids have same reference state hjunk=hrefdf(1) ldef=.true. do 202 i=1,ncomp if (hjunk.ne.hrefdf(i)) then ldef=.false. end if 202 continue if (.not.ldef) then ierrsv=-28 herrsv='[SETREF warning -28] sorry, cannot apply reference ' & //'state (DEF) at a specified mixture composition ' & //'when components have different default reference ' & //'states; will apply defaults to each pure component' & //hnull call ERRMSG (ierrsv,herrsv) ixflg=1 end if end if c c initialize reference temperature, entropy, and enthalpy c * write (*,*) ' SETREF--240 loop from/to: ',-ncomp,ncomp do 240 i=-ncomp,ncomp if (hrf.eq.'IIR' .or. hrf.eq.'ASH' .or. hrf.eq.'NBP' .or. & hrf.eq.'OTH') then hrefi(i)=hrf ixflg=ixflag else !DEFault or unknown reference state hrefi(i)=hrefdf(i) ixflg=1 !must apply ref state to pures only if (hrf.ne.'DEF') then ierr=-29 herr='[SETREF warning -29] reference state not found, using' & //' default reference state'//hnull call ERRMSG (ierr,herr) end if end if c c write (*,*) ' SETREF--i,hrefi: ',i,' ',hrefi(i) if (hrefi(i).eq.'IIR') then c IIR reference state c write (*,*) ' SETREF--IIR ref state with ixflag =',ixflg tref(i)=273.15d0 prefi(i)=-1.0d0 !code for reference state of sat liquid href(i)=200.0d0*wm(i) !reference enthalpy is 200.0 kJ/kg sref(i)=wm(i) !reference entropy is 1.0 kJ/(kg.K) else if (hrefi(i).eq.'ASH') then c ASHRAE reference state c write (*,*) ' SETREF--ASH ref state with ixflag =',ixflg tref(i)=233.15d0 prefi(i)=-1.0d0 !code for reference state of sat liquid href(i)=0.0d0 sref(i)=0.0d0 else if (hrefi(i).eq.'NBP') then c normal boiling point referene state c write (*,*) ' SETREF--NBP ref state with ixflag =',ixflg tref(i)=tnbp(i) prefi(i)=-1.0d0 !code for ref state of saturated liquid sref(i)=0.0d0 href(i)=0.0d0 else if (hrefi(i).eq.'OTH') then c 'other' reference state--use input values of h0, s0, t0, p0 c write (*,*) ' SETREF--OTHer ref state with ixflag =',ixflg if (hrf.eq.'DEF') then c if input to SETREF is 'DEF' reference t,p,h,s come from fluid file tref(i)=tdef(i) prefi(i)=pdef(i) href(i)=hdef(i) sref(i)=sdef(i) else c otherwise, reference t,p,h,s are inputs to SETREF tref(i)=t0 prefi(i)=p0 href(i)=h0 sref(i)=s0 end if else if (hrefi(i).eq.'NUL') then c this case should be encountered only for i<0 when any model other c than the ECS model is specified; fill arrays with dummy values tref(i)=300.0 prefi(i)=1.0d-6 sref(i)=0.0d0 href(i)=0.0d0 else c unknown reference state (use normal boiling point) ierr=-29 write (herr,1238) hrefi(i),i,hnull 1238 format ('[SETREF warning -29] reference state (',a3,') for ', & 'component ',i3,' not found; will use ', & 'NBP reference state. ',a1) call ERRMSG (ierr,herr) tref(i)=tnbp(i) prefi(i)=-1.0d0 !code for ref state of saturated liquid sref(i)=0.0d0 href(i)=0.0d0 end if 240 continue c c check limits of equation of state against reference state temperature c do 280 i=1,ncomp Dx=0.0d0 px=0.0d0 do 278 j=1,ncomp x(j)=0.0d0 278 continue x(i)=1.0d0 call LIMITX ('EOS',t,Dx,px,x,tmin,tmax,Dmax,pmax,ierr2,herr2) c if (hrefi(i).eq.'IIR') then c IIR reference state if (tmin.gt.273.15d0) then c minimum temperature of EOS is above IIR reference state temperature c issue error message ierr=22 write (herr,1022) i,tmin,hrefdf(i),hnull 1022 format ('[SETUP warning 22] ', & 'minimum temperature for component ',i3,' is greater ', & 'than IIR reference temperature; Tmin =',1pe11.3, & ' K, Tref = 273.15 K.; will use default (',a3,'). ',a1) else if (tc(i).lt.274.15d0) then c critical temperature is too low to use IIR reference state c issue error message ierr=23 write (herr,1023) i,tc(i),hrefdf(i),hnull 1023 format ('[SETUP warning 23] ', & 'critical temperature for component ',i3,' is less ', & 'than IIR reference temperature; Tcrit =',1pe11.3, & ' K, Tref = 273.15 K; will use default (',a3,'). ',a1) end if end if c else if (hrefi(i).eq.'ASH') then c ASHRAE reference state if (tmin.gt.233.15d0) then c minimum temperature of EOS is above reference state temperature ierr=24 write (herr,1024) i,tmin,hrefdf(i),hnull 1024 format ('[SETUP warning 24] ', & 'minimum temperature for component ',i3,' is greater ', & 'than ASHRAE reference temperature; Tmin =',1pe11.3, & ' K, Tref = 233.15 K; will use default (',a3,'). ',a1) else if (tc(i).lt.234.15d0) then c critical temperature is too low to use IIR reference state ierr=25 write (herr,1025) i,tc(i),hrefdf(i),hnull 1025 format ('[SETUP warning 25] ', & 'critical temperature for component ',i3,' is less ', & 'than ASHRAE reference temperature; Tcrit =',1pe11.3, & ' K, Tref = 233.15 K; will use default (',a3,'). ',a1) end if end if else if (hrefi(i).eq.'NBP') then c Normal Boiling Point reference state if (tmin.gt.tnbp(i)) then c minimum temperature of EOS is above boiling point temperature c issue error message (this is a problem with CO2) ierr=26 write (herr,1026) i,tmin,tnbp(i),hrefdf(i),hnull 1026 format ('[SETUP warning 26] ', & 'minimum temperature for component ',i3,' is greater ', & 'than boiling point temperature; Tmin =',1pe11.3, & ' K, Tnbp =',1pe11.3,' K; will use default (',a3,'). ',a1) end if end if c if (ABS(ierr).ge.22) then c if specified reference state is outside limits for component i, c use default state (read in from fluid file); this code assumes c that the specification in the fluid file is valid call ERRMSG (ierr,herr) if (hrefdf(i).eq.'IIR') then * write (*,1282) i,hrefi(i),hrefdf(i) *1282 format (2x,'SETREF--specified reference state for component' * & ,i4,' (',a3,'), is not valid; will apply default ' * & ,'reference state (',a3,') to pures. ') tref(i)=273.15d0 prefi(i)=-1.0d0 href(i)=200.0d0*wm(i) !reference enthalpy is 200.0 kJ/kg sref(i)=wm(i) !reference entropy is 1.0 kJ/(kg.K) else if (hrefdf(i).eq.'ASH') then * write (*,1282) i,hrefi(i),hrefdf(i) tref(i)=233.15d0 prefi(i)=-1.0d0 href(i)=0.0d0 sref(i)=0.0d0 else if (hrefdf(i).eq.'NBP') then * write (*,1282) i,hrefi(i),hrefdf(i) tref(i)=tnbp(i) prefi(i)=-1.0d0 href(i)=0.0d0 sref(i)=0.0d0 else if (hrefdf(i).eq.'OTH') then * write (*,1282) i,hrefi(i),hrefdf(i) tref(i)=tdef(i) prefi(i)=pdef(i) href(i)=hdef(i) sref(i)=sdef(i) else !unknown reference state; set to reasonble values * write (*,*) ' SETREF--default reference state is unknown' tref(i)=0.8d0*tc(i) prefi(i)=-1.0d0 href(i)=0.0d0 sref(i)=0.0d0 end if hrefi(i)=hrefdf(i) ixflg=1 !must apply ref state to pures only end if c save any error flags/messages to pass up the chain; then reset ierrsv=ierr herrsv=herr ierr=0 herr=hnull 280 continue c c apply reference enthalpy and entropy to each component; this for c case of ixflag =1; also used as initialization for ixflag = 2 c do 480 i=1,ncomp do 476 j=1,ncomp x(j)=0.0d0 !initialize all compositions to zero 476 continue x(i)=1.0d0 !if mixture, calculate pure component i c write (*,1002) i,tref(i),href(i),sref(i) c1002 format (1x,' SETREF--i,tref,href,sref: ',i3,f7.2,f10.1,f10.3) if (hrefi(i).eq.'NBP') then c find bubble point at one atmosphere for each pure component pnbp=101.325d0 kbub=1 c write (*,*) ' SETREF--i,x(i): ',i,x(i) call SATP (pnbp,x,kbub,tbub,rhol,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) rhoref(i)=rhol tref(i)=tbub else if (prefi(i).lt.0.0d0) then c find saturated liquid or vapor density (includes ASH and IIR states) if (prefi(i).lt.-1.99d0) then kbub=2 call SATT (tref(i),x,kbub,p,rhol,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) rhoref(i)=rhov else kbub=1 call SATT (tref(i),x,kbub,p,rhol,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) rhoref(i)=rhol end if else if (hrefi(i).eq.'OTH') then c find density at specified t,p (either input to SETREF or default value) t=tref(i) p=prefi(i) call TPFLSH (t,p,x,D,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,herr) rhoref(i)=D end if if (ierr.ge.1) then c error encountered in calculating reference state--set reference c density to critical density? rhoref(i)=rhoc(i) ierr=119 write (herr,1119) i,herr(1:162),hnull 1119 format ('[SETREF error 119] convergence failure in ', & 'calculating reference state for component #',i3, & ': ',a162,a1) end if call THERM (tref(i),rhoref(i),x,ptherm,e,h,s,cv,cp,w,hjt) sref(i)=s !store reference entropy in common href(i)=h !ditto for enthalpy c write (*,1004) i,tref(i),href(i),sref(i),rhoref(i) c1004 format (1x,' SETREF--after SATT, etc.: ',i3,f7.2,f10.1,2f10.3) 480 continue end if !end of initialization for pure components c if (ixflg.eq.2) then c c reference state applied to mixture of specified composition c this option is not available if the 'DEFault' state is chosen c if (hrf.eq.'OTH') then tmix=t0 pref=p0 smix=s0 hmix=h0 else if (hrf.eq.'IIR') then c IIR reference state tmix=273.15d0 pref=-1.0d0 !code for ref state of saturated liquid smix=1.0d0 !reference entropy is 1.0 kJ/(kg.K) hmix=200.0d0 !reference enthalpy is 200. kJ/kg else if (hrf.eq.'ASH') then c ASHRAE reference state tmix=233.15d0 pref=-1.0d0 !code for reference state of sat liquid smix=0.0d0 hmix=0.0d0 else c boiling point or unknown ref state c find bubble point of specified mixture pnbp=101.325d0 call SATP (pnbp,x0,1,tbub,Dmix,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) tmix=tbub pref=-1.0d0 !code for reference state of sat liquid smix=0.0d0 hmix=0.0d0 end if c check mix reference temperature against limits, if outside limits c take average of component Tref's (would be different only if one or c more components were outside limits) call LIMITX ('EOS',tmix,Dx,px,x0,tmin,tmax,Dmx,pmx,ie2,herr2) call ERRMSG (ie2,herr2) if (tmix.lt.tmin .or. tmix.gt.tmax) then tmix=0.0d0 do 670 i=1,ncomp tmix=tmix+x0(i)*tref(i) 670 continue end if c predefined reference states are on mass basis wmm=WMOL(x0) smix=smix*wmm hmix=hmix*wmm end if c if (pref.lt.0.0d0) then c find liquid or vapor density (includes IIR and ASHrae states) if (pref.lt.-1.99d0) then kbub=2 call SATT (tmix,x0,kbub,p,rhol,Dmix,xl,xv,ierr,herr) call ERRMSG (ierr,herr) else kbub=1 call SATT (tmix,x0,kbub,p,Dmix,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) end if else c find density at specified t0, p0 call LIMITX ('EOS',t0,Dx,p0,x0,tmin,tmax,Dmx,pmx,ierr,herr2) call ERRMSG (ierr,herr2) if (ierr.ge.1) then ierr=27 herr='[SETREF error 27] specified reference state for '// & 'mixture is outside equation limits: '//herr2(1:170)// & hnull RETURN end if call TPFLSH (t0,p0,x0,D,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,herr) tmix=t0 Dmix=D end if c call THERM (tmix,Dmix,x0,ptherm,e,h,s,cv,cp,w,hjt) do 680 i=1,ncomp c modify sref, href by amount needed for specified mixture c to agree with s0, h0 c write (*,1005) s,smix,h,hmix c1005 format (1x,' SETREF s,smix,h,hmix: ',2f10.4,2f10.1) sref(i)=sref(i)+s-smix href(i)=href(i)+h-hmix c write (*,1006) i,tref(i),href(i),sref(i),rhoref(i) c1006 format (1x,' SETREF--after mix adjust: ',i3,f7.2,f10.1,2f10.3) 680 continue end if c c variable ksetrf is flag to indicate that SETREF has been called c it is reset to 0 each time SETUP is called (indicating that SETREF c must be recalled) ksetrf=1 c c copy inputs to common block c hrefst=hrf ixfsav=ixflg if (hrf.eq.'OTH') then c save h0, s0, t0, p0 only if 'OTHer' reference state is used h0sav=h0 s0sav=s0 t0sav=t0 p0sav=p0 end if if (ixflg.eq.2) then c save x0 only if ixflg = 2 do 720 i=1,ncomp 720 x0sav(i)=x0(i) end if c c set value of hrf if 'DEF' was input value c c if (hrf.eq.'DEF') then c hrf=hrefi(1) c do 800 i=1,ncomp c if (hrefi(i).ne.hrf) then c hrf='MIX' c end if c 800 continue c end if c c reset error flag/message (the more significant errors are those which c arise from out-of-bounds conditions, not any potential errors in calls c to SATT, etc.) if (ierrsv.ne.0) then ierr=ierrsv herr=herrsv end if c RETURN c end !subroutine SETREF c c ====================================================================== c subroutine SETMOD (nc,htype,hmix,hcomp,ierr,herr) c c set model(s) other than the NIST-recommended ('NBS') ones c c This subroutine must be called before SETUP; it need not be called c at all if the default (NIST-recommended) models are desired. c c inputs: c nc--number of components (1 for pure fluid) [integer] c htype--flag indicating which models are to be set [character*3] c 'EOS': equation of state for thermodynamic properties c 'ETA': viscosity c 'TCX': thermal conductivity c 'STN': surface tension c 'NBS': reset all of the above model types and all c subsidiary component models to 'NBS'; c values of hmix and hcomp are ignored c hmix--mixture model to use for the property specified in the c argument htype [character*3]; c this can a pure fluid model if number of components = 1 c 'NBS': use NIST recommendation for specified fluid/mixture c (valid input for all properties, pure or mix) c c some allowable choices for thermodynamic properties: c pure fluids: c 'FEQ': Helmholtz free energy model c 'BWR': pure fluid modified Bennedict-Webb-Rubin (MBWR) c 'ECS': pure fluid thermo extended corresponding states c mixture: c 'HMX': mixture Helmholtz model for thermodynamic properties c c some allowable choices for viscosity: c pure fluids: c 'ECS': extended corresponding states (all fluids) c 'VS1': the 'composite' model for R134a, R152a, NH3, etc. c 'VS2': Younglove-Ely model for hydrocarbons c mixture: c 'ECS': extended corresponding states c c some allowable choices for thermal conductivity: c pure fluids: c 'ECS': extended corresponding states (all fluids) c 'TC1': the 'composite' model for R134a, R152a, etc. c 'TC2': Younglove-Ely model for hydrocarbons c mixture: c 'ECS': extended corresponding states c c some allowable choices for surface tension: c pure fluids: c 'ST1': surface tension as f(tau); tau = 1 - T/Tc c mixture: c 'STX': surface tension mixture model c c hcomp--component model(s) to use for property specified in htype c [array (1..nc) of character*3] c 'NBS': NIST recommendation for specified fluid/mixture c (valid input for all properties, pure or mix) c other choices are identical to the pure fluid c choices listed above c outputs: c ierr--error flag: 0 = successful c no errors are returned from this routine, c ierr and herr included to maintain parallel c structure and for possible future use c herr--error string (character*255 variable if ierr<>0) c [fluid parameters, etc. returned via various common blocks] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 03-20-96 MM, original version c 03-21-96 MM, replace /MODEL/ with /EOSMOD/, /STNMOD/ c 03-28-96 MM, add option if htype='NBS' reset all models c 03-28-97 MM, add initialization of hnull (repeat of SETUP initialization); c add /RESETM/ to check for initialization of models c 03-31-97 MM, reset all models if lreset true, even for htype<>'NBS' c 04-10-97 MM, add and initialize /CREMOD/ (critical enhancements) c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: SETMOD c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 htype,hmix,hcomp(1:ncmax) character*3 heos,hpheq,hmxeos,hsten,hstenk,hmodcp character*3 heta,hetak,htcx,htcxk character*3 hetacr,htcxcr character*255 herr logical lreset c common /NCOMP/ ncomp common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /TRNMOD/ heta,hetak(nrf0:nx),htcx,htcxk(nrf0:nx) common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax) common /STNMOD/ hsten,hstenk(n0:nx) common /CHAR/ htab,hnull common /RESETM/ lreset !flag indicating need to reset all models c c repeat initialization of following from SETUP in case SETMOD is called c before first call to SETUP hnull=CHAR(0) !null character to terminate error strings ierr=0 herr=hnull c if (htype.eq.'NBS' .or. htype.eq.'nbs' .or. lreset) then c reset all models * write (*,*) ' SETMOD--about to reset all models to NBS' heos='NBS' hpheq='NBS' !phase equilibria model, not currently implemented do 20 i=n0,ncmax hmxeos(i)='NBS' hmodcp(i)='NBS' 20 continue heta='NBS' do 40 i=nrf0,ncmax hetak(i)='NBS' hetacr(i)='NUL' 40 continue htcx='NBS' do 60 i=nrf0,ncmax htcxk(i)='NBS' htcxcr(i)='NUL' 60 continue hsten='NBS' do 80 i=n0,ncmax hstenk(i)='NBS' 80 continue lreset=.false. end if if (htype.eq.'EOS' .or. htype.eq.'eos') then c equation of state specification heos=hmix do 100 i=1,nc hmxeos(i)=hcomp(i) 100 continue else if (htype.eq.'ETA' .or. htype.eq.'eta') then c viscosity specification heta=hmix do 200 i=1,nc hetak(i)=hcomp(i) 200 continue else if (htype.eq.'TCX' .or. htype.eq.'tcx') then c viscosity specification htcx=hmix do 300 i=1,nc htcxk(i)=hcomp(i) 300 continue else if (htype.eq.'STN' .or. htype.eq.'stn') then c surface tension specification hsten=hmix do 400 i=1,nc hstenk(i)=hcomp(i) 400 continue end if * write (*,1003) heos,(hmxeos(j),j=n0,ncmax), * & heta,(hetak(j),j=nrf0,nx), * & htcx,(htcxk(j),j=nrf0,nx), * & hsten,(hstenk(j),j=n0,ncmax) *1003 format (' SETMOD--heos, hmxeos: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)/ * & ' (exit) visc models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' t.c. models: ',a3,25x,a3,3x,5(a3,1x)/ * & ' surf ten models: ',a3,3x,5(a3,1x),a5,3x,5(a3,1x)) c RETURN c end !subroutine SETMOD c c ====================================================================== c subroutine SETKTV (icomp,jcomp,hmodij,fij,hfmix,ierr,herr) c c set mixture model and/or parameters c c This subroutine must be called after SETUP, but before any call to c SETREF; it need not be called at all if the default mixture c parameters (those read in by SETUP) are to be used. c c inputs: c icomp--component i c jcomp--component j c hmodij--mixing rule for the binary pair i,j [character*3] c e.g. 'LJ1' (Lemmon-Jacobsen model) c 'LM1' (modified Lemmon-Jacobsen model) or c 'LIN' (linear mixing rules) c 'RST' indicates reset all pairs to values from c original call to SETUP (i.e. those read from file) c [all other inputs are ignored] c fij--binary mixture parameters [array of dimension nmxpar; c currently nmxpar is set to 6] c the parameters will vary depending on hmodij; c for example, for the modified Lemmon-Jacobsen model (LM1): c fij(1) = Kt c fij(2) = Kv c fij(3) = Fpq c fij(4) not used c fij(5) not used c fij(6) not used c hfmix--file name [character*80] containing generalized parameters c for the binary mixture model; this will usually be the same c as the corresponding input to SETUP (e.g.,':fluids:HMX.bnc') c outputs: c ierr--error flag: 0 = successful c 111 = error in opening mixture file c 112 = mixture file of wrong type c -113 = illegal i,j specification c (i = j or i > nc or j > nc) c herr--error string (character*255 variable if ierr<>0) c [mixture parameters returned via various common blocks] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-07-96 MM, original version c 05-08-96 MM, add /MXINFO/ and load corresponding information c 11-04-96 MM, change nmxpar from 4 to 6 c 11-22-96 MM, 'RST' option now calls SETHMX, i.e. resets all pairs c 11-25-96 MM, reorder such that 'RST' option is tested first c 10-01-97 MM, add compiler switches to allow access by DLL c 02-06-98 MM, initial use of icount unitialized, use nbin instead c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: SETKTV c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing func character*1 htab,hnull character*3 hmodij,hmodmx,htype character*3 hpheq,heos,hmxeos,hmodcp character*12 hname character*80 hfmix character*255 herr c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin c dimension fij(nmxpar) logical lij(ncmax,ncmax) c common /NCOMP/ ncomp common /CHAR/ htab,hnull common /CNAM/ hname(n0:nx) common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary c the namix(i,j) are the number of a(i,j,k,1..3) terms c the iexpmx(i,j,k) are powers of the exponentials in the general term common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) c /MXINFO/ contains information on the mixing rules and parameters c hmfile specifies the files from which mixing rules originate c (stored in order 1,2 1,3 2,3; element zero is file called in SETUP) c hbin provides documentation for the current binary parameters c hrule contains descriptions of the currently available mixing rules c hbpar contains descriptions of the binary parameters (e.g. Kt, Kv) c associated with the currently available mixing rules common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) common /MXRULE/ nrule c common to pass flags associated with reference state common /IRFSAV/ ixfsav,ksetrf c ierr=0 herr=hnull c c set flag indicating that reference state needs to be recalculated c (for use with the ixflag = 2 option of SETREF) ixfsav=0 c special case: reset parameters to original values if (hmodij.eq.'RST' .or. hmodij.eq.'rst') then c retrieve file specified on original call to SETUP hfmix=hmfile(0) call SETHMX (hfmix,ierr,herr) c do 100 i=1,icount do 100 i=1,nbin hmfile(i)=hmfile(0) 100 continue * write (*,*) ' SETKTV--reset all parameters to original values' * write (*,*) ' i j mod Kt Kv Fpq alpha' d do 120 i=i,ncomp d do 110 j=1,ncomp * write (*,1010) i,j,hmodmx(i,j),(fmix(i,j,k),k=1,nmxpar) d 110 continue d 120 continue *1010 format (1x,2i3,2x,a3,4f10.5) RETURN c c replace following with above call to SETHMX; this changes behavior c of routine for case of hmodij = 'RST' to reset all pairs, this c change eliminates the following code (which is nearly duplicate with c corresponding code in SETHMX) at expense of loss of option of c resetting a single pair c c hmfile(ibin)=hmfile(0) c* write (*,*) 'SETKTV--reset parameters for i,j =',icomp,jcomp c* write (*,1100) hfmix c*1100 format (/1x,'reading mixture coefficients from file: (',a40,')') c nread=12 !logical unit for file reads c open (unit=nread,file=hfmix,err=999) c rewind (nread) !rewind file (RDBNC may have left it at the end) c read (nread,2003) htype c if (htype.ne.heos) then c ierr=112 c herr='[SETKTV error 112] mixture file of wrong type'//hnull c call ERRMSG (ierr,herr) cc rewind and close the file c rewind (nread) c close (nread) c RETURN c end if c call RDBNC (nread,icomp,jcomp,hmodij,fij,hbinp,ierr,herr) c if (ierr.ne.0) then cc binary parameter(s) not found in file, must estimate them c call ESTBNC (icomp,jcomp,hmodij,fij,hbinp,ierr,herr) c* write (*,1102) icomp,jcomp,hmodij,(fij(k),k=1,nmxpar) c*1102 format (1x,' SETKTV--outputs from ESTBNC: ',2i3,2x,a3,8f8.4) c end if c hbin(ibin)=hbinp c hmodmx(icomp,jcomp)=hmodij c hmodmx(jcomp,icomp)=hmodij !pair j,i is same as i,j c do 160 k=1,nmxpar c fmix(icomp,jcomp,k)=fij(k) c fmix(jcomp,icomp,k)=fij(k) !pair j,i is same as i,j c 160 continue cc read coefficients for mixing rules from file c call RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr) c if (.not.lij(icomp,jcomp)) then c ierr=-117 c herr=' SETHMX ERROR--mixing term not found for one or '// c & 'more binary pair(s)'//hnull c call ERRMSG (ierr,herr) c end if cc rewind and close the file c rewind (nread) c close (nread) c RETURN end if c c check if icomp or jcomp > ncomp, also do not allow change of c pure fluid or i=j parameters c if (icomp.gt.ncomp .or. jcomp.gt.ncomp) then ierr=-113 herr='[SETKTV warning -113] i > ncomp and/or j > ncomp'//hnull * write (*,*) herr RETURN else if (ncomp.le.1 .or. icomp.eq.jcomp) then ierr=-113 herr='[SETKTV warning -113] binary parameters for i = j '// & 'cannot be changed.'//hnull * write (*,*) herr RETURN end if c c match input icomp,jcomp with binary pair number icount=0 ibin=1 !initialize only do 144 i=1,ncomp-1 do 140 j=i+1,ncomp icount=icount+1 if ((i.eq.icomp .and. j.eq.jcomp) .or. & (i.eq.jcomp .and. j.eq.icomp)) then ibin=icount end if 140 continue 144 continue c do 200 k=1,nmxpar fmix(icomp,jcomp,k)=fij(k) fmix(jcomp,icomp,k)=fij(k) !pair (j,i) is same as (i,j) 200 continue hbin(ibin)='The binary parameter(s) for ('// & hname(icomp)//'+'//hname(jcomp)// & ') have been modified from their original values.' c c check if the mixture model has changed for the binary pair, if so, c must read coefficients for that model c if (hmodij.ne.hmodmx(icomp,jcomp)) then hmodmx(icomp,jcomp)=hmodij hmodmx(jcomp,icomp)=hmodij !pair (j,i) is same as (i,j) * write (*,1105) hfmix *1105 format (/1x,' SETKTV--mixture parameters from file: (',a80,')') nread=12 !logical unit for file reads open (unit=nread,file=hfmix,err=999) rewind (nread) !rewind file (RDBNC may have left it at the end) read (nread,2003) htype if (htype.ne.heos) then ierr=112 herr='[SETKTV error 112] mixture file of wrong type'//hnull * write (*,*) ' SETKTV--heos,htype: (',heos,'),(',htype,')' call ERRMSG (ierr,herr) RETURN end if c store file containing mixture parameters hmfile(ibin)=hfmix * write (*,*) ' SETKTV--new file for ibin =',ibin,' ',hmfile(ibin) call RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr) c rewind and close the file rewind (nread) close (nread) end if c * write (*,*) ' SETKTV: return' * write (*,*) ' i j mod Kt Kv Fpq alpha' d i=icomp d j=jcomp * write (*,1280) i,j,hmodmx(i,j),(fmix(i,j,k),k=1,nmxpar) *1280 format (1x,2i3,2x,a3,6f10.5) RETURN c 999 ierr=111 write (herr,2111) hfmix,hnull 2111 format ('[SETKTV error 111] error in opening mixture file', & '; filename = (',a80,').',a1) call ERRMSG (ierr,herr) RETURN c 2003 format (a3) c end !subroutine SETKTV c c ====================================================================== c subroutine GETKTV (icomp,jcomp,hmodij,fij,hfmix,hfij,hbinp,hmxrul) c c retrieve mixture model and parameter info for a specified binary c c This subroutine should not be called until after a call to SETUP. c c inputs: c icomp--component i c jcomp--component j c outputs: c hmodij--mixing rule for the binary pair i,j (e.g. LJ1 or LIN) c [character*3] c fij--binary mixture parameters [array of dimension nmxpar; c currently nmxpar is set to 6]; the parameters will vary c depending on hmodij; c hfmix--file name [character*80] containing parameters for the c binary mixture model c hfij--description of the binary mixture parameters [character*8 c array of dimension nmxpar] c for example, for the modified Lemmon-Jacobsen model (LM1): c fij(1) = Kt c fij(2) = Kv c fij(3) = Fpq c fij(4) = 'not used' c fij(5) = 'not used' c fij(6) = 'not used' c hbinp--documentation for the binary parameters [character*255] c terminated with ASCII null character c hmxrul--description of the mixing rule [character*80] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-08-96 MM, original version c 08-01-96 MM, bug fix: ibin not defined for icomp=jcomp c 11-04-96 MM, change nmxpar from 4 to 6 c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: GETKTV c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing func character*1 htab,hnull character*3 hmodij,hmodmx character*3 hpheq,heos,hmxeos,hmodcp character*8 hfij(nmxpar) character*80 hfmix,hmxrul character*255 hbinp c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin,herr c dimension fij(nmxpar) c common /NCOMP/ ncomp common /CHAR/ htab,hnull common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary c the namix(i,j) are the number of a(i,j,k,1..3) terms c the iexpmx(i,j,k) are powers of the exponentials in the general term common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) c /MXINFO/ contains information on the mixing rules and parameters c hmfile specifies the files from which mixing rules originate c (stored in order 1,2 1,3 2,3; element zero is file called in SETUP) c hbin provides documentation for the current binary parameters c hrule contains descriptions of the currently available mixing rules c hbpar contains descriptions of the binary parameters (e.g. Kt, Kv) c associated with the currently available mixing rules common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) common /MXRULE/ nrule c hmodij=hmodmx(icomp,jcomp) do 80 k=1,nrule if (hmodij.eq.hrule(k)(1:3)) then irule=k goto 82 end if 80 continue * write (*,*) ' GETKTV ERROR--mixing rule not found: ',hmodij ierr=-117 write (herr,2117) icomp,jcomp,hnull 2117 format ('[GETKTV error 117] mixing rule not found for components', & i3,' +',i3,'.',a1) call ERRMSG (ierr,herr) RETURN 82 continue ibin=0 if (icomp.eq.jcomp .or. ncomp.le.1) then hbinp='trivial case--identical components'//hnull else c match input icomp,jcomp with binary pair number icount=0 do 92 i=1,ncomp-1 do 90 j=i+1,ncomp icount=icount+1 if ((i.eq.icomp .and. j.eq.jcomp) .or. & (i.eq.jcomp .and. j.eq.icomp)) then ibin=icount c GUI combines hbinp with 3-letter code, so can use only 251 char hbinp=hbin(ibin)(1:251)//hnull end if 90 continue 92 continue end if c file containing mixture parameters c write (*,*) ' GETKTV--file for ibin =',ibin,' ',hmfile(ibin) hfmix=hmfile(ibin) hmxrul=hrule(irule)(1:79)//hnull do 100 k=1,nmxpar fij(k)=fmix(icomp,jcomp,k) hfij(k)=hbpar(irule,k) 100 continue c * write (*,1278) hfmix,hbinp,hmxrul *1278 format (/1x,' GETKTV hfile: ',a80/ * & 1x,' hbinp: ',a255/ * & 1x,' hmxrul: ',a80) * write (*,1280) (hfij(k),k=1,nmxpar), * & icomp,jcomp,hmodmx(icomp,jcomp),(fmix(icomp,jcomp,k),k=1,nmxpar) *1280 format (1x,' i j mod ',6(a8,2x)/ * & 1x,2i3,2x,a3,6f10.5) RETURN c end !subroutine GETKTV c c ====================================================================== c subroutine RFFILE (hfilei,hflref) c c provide a full path specification for the '.fld' file for the ECS c reference fluid given a specification for a component '.fld' file c e.g. 'C:\NISTPROP\fluids\R123.fld' and the reference fluid pointer c which is stored in that file, e.g. 'R134a.fld' c c N.B. This routine assumes that the reference fluid '.fld' file c is in the same directory as the component file(s) c c inputs: c hfilei--file name specifying mixture component (character*80) c hflref--file name specifying reference fluid (character*80) c output: c hflref--file name specifying reference fluid (character*80) c the output value of hflref consists of the directory c information gleaned from hfilei concatenated with the c input value of hflref c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 12-14-95 MM, original version c 12-18-95 MM, use htemp to concatenate character strings (cannot c concatenate string with itself on some compilers) c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) character*80 hfilei,hflref,htemp c if (hfilei(1:5).eq.'BDATA' .or. hfilei(1:5).eq.'bdata') then c special case if component data is read from block data; c assume that reference fluid is in same directory as main program c i.e., hflref is simply fluid name + '.fld', e.g. 'R134a.fld' c --or-- c the pointer to the reference fluid in block data is also to c block data, e.g. 'BDATA:811-97-2' for R134a else c usual case of fluid read from '.fld' file ncfile=LEN(hfilei) jpoint=0 do 100 i=1,ncfile j=1+ncfile-i c search for directory delimiters c (\,:,/ for DOS, Macintosh, and Unix, respectively) jchar=ICHAR(hfilei(j:j)) if (jchar.eq.92 .or. jchar.eq.58 .or. jchar .eq.47) then jpoint=j goto 101 end if 100 continue 101 continue !file delimiter character has been found c write (*,1004) ncfile,jpoint c1004 format (1x,' RFFILE--ncfile, jpoint: ',2i4) if (jpoint.le.0) then c no path name provided in specification of component .fld file c assume that reference fluid is in same directory as main program c i.e., hflref is simply fluid name + '.fld', e.g. R134a.fld else htemp=hfilei(1:jpoint)//hflref hflref=htemp end if end if * write (*,1005) hflref *1005 format (1x,' RFFILE--full path for reference fld file: ',a80) c RETURN end !subroutine RFFILE c c ====================================================================== c block data BDSET c c This block data initializes the common blocks which save input c arguments to SETUP and SETREF. This is done to check if arguments c are unchanged since previous call; if so, routine can be bypassed. c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 01-10-96 MM, original version c 01-16-96 MM, add saves for SETREF c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 03-20-96 MM, changes to accomodate SETMOD: initialize /TRNMOD/ c initialize heos, etc to 'NBS' c 03-21-96 MM, replace /MODEL/ with /EOSMOD/, /STNMOD/, /REFST/ c 03-22-96 MM, initialize heta,htcx,hsten to 'NBS', rather than 'NUL' c 03-27-96 MM, ditto for component models, hmxeos,heta,htcxk,hstenk c 02-20-97 MM, initialize default ref states in /REFST/ c 03-28-97 MM, add /RESETM/ and set lreset=.true. so that models are c initialized by call to SETMOD rather than here c 07-15-97 MM, move /FLAGS/ here from ftn_pas.f c 11-25-97 MM, insert 2 extra 9's into the funny numbers in /FLAGS/ c 12-16-97 MM, set last open variable in /FLAGS/ to 'error' c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nx0=nx-n0+1) character*3 hrefst,hrefdf character*3 hph,heqn,heqnk,hvis,hvisk,hcond,hcondk,hsurf,hsurfk character*80 hfsav,hmxsav logical lreset c common /NCOMP/ ncomp common /REFST/ hrefst,hrefdf(n0:nx) common /SETSAV/ hfsav(n0:nx),hmxsav common /REFSAV/ x0sav(nx),h0sav,s0sav,t0sav,p0sav common /IRFSAV/ ixfsav,ksetrf common /MODSAV/ hph,heqn,heqnk(nx),hvis,hvisk(nx), & hcond,hcondk(nx),hsurf,hsurfk(nx) common /RESETM/ lreset !flag indicating need to reset all models c flags indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,xerr,xnotd,xnotc c data ncomp /0/ data lreset /.true./ !reset all models on initial call to SETUP data hfsav /nx0*'NUL'/ data hmxsav /'NUL'/ data hrefst /'NUL'/ c initialize commons which save inputs to SETUP from one call to next data x0sav /nx*0.0d0/ data h0sav,s0sav,t0sav,p0sav /4*0.0d0/ data ixfsav,ksetrf /2*0/ data hph,heqn,heqnk /2*'NUL',nx*'NUL'/ data hvis,hvisk /'NUL',nx*'NUL'/ data hcond,hcondk /'NUL',nx*'NUL'/ data hsurf,hsurfk /'NUL',nx*'NUL'/ data hrefdf /nx0*'NUL'/ data xnota /-9.99991d6/ !flag indicating 'not applicable' data x2ph /-9.99992d6/ !flag indicating '2-phase' data xsubc /-9.99993d6/ !flag indicating 'subcooled' data xsuph /-9.99994d6/ !flag indicating 'superheated' data xsupc /-9.99995d6/ !flag indicating 'supercritical' data xinf /-9.99996d6/ !flag indicating 'infinite' data xerr /-9.99997d6/ !flag indicating 'error' data xnotd /-9.99998d6/ !flag indicating property 'not defined' data xnotc /-9.99999d6/ !flag indicating property 'not calculated' c end !block data BDSET c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file setup.f c ======================================================================