c  begin file setup2.f
c
c  This file contains routine which initialize the models, fluid-
c  specific parameters and coefficients, etc.  The subroutine SETUP (in
c  file setup.f) must be called before any of the other property routines
c  are called.  Call(s) to SETMOD and SETREF are optional and may be used
c  to specify non-standard models and reference states.
c
c  If the transport models, dielectric models, and surface tension eqs.
c  (vis, tcx, st, de) are not required, and only the default values in the
c  fluid files are used, the following files do not have to be compiled and
c  linked with the other source code:  CORE_ANC, CORE_DE, CORE_STN, FTN_PAS,
c  SETUP2, TRNS_ECS, TRNS_TCX, TRNS_VIS, TRNSP, X_SUB.  In addition,
c  the line in SETUP.FOR calling STFLD2 must be commented out.
c  The file FLASH2.FOR is not required but contains additional property
c  calculation routines.
c
c  contained here are:
c     subroutine STFLD2 (nread,i,hcasn,hcite,heqn,hetemp,hfile,hflag,
c    &                    href,hstar,httemp,htype,leta,ltcx,ierr,herr)
c     subroutine GETMOD (icomp,htype,hcite)
c     subroutine SETKTV (icomp,jcomp,hmodij,fij,hfmix,ierr,herr)
c     subroutine GETKTV (icomp,jcomp,hmodij,fij,hfmix,hfij,hbinp,hmxrul)
c     subroutine GETFIJ (hmodij,fij,hfij,hmxrul)
c
c
c ======================================================================
c
      subroutine STFLD2 (nread,i,hcasn,hcite,heqn,hetemp,hfile,hflag,
     &                    href,hstar,httemp,htype,leta,ltcx,ierr,herr)
c
c  define transport, dielectric, and surface tension models
c  see comments to SETFLD for additional information
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-19-94  MM, original version (SETFLD)
c  09-01-00 EWL, remove pieces from SETFLD and create new SETFLD2 subroutine
c  11-16-01 MLH, activated TK6 model for ecs crit enh; added common block CREMOD2
c  01-23-02 EWL, split common block CITE into 2 pieces for some compilers
c  08-27-04 EWL, rename SETFLD2 to STFLD2 to conform with standard f77 code
c  07-22-05 MLH, allow multiple transport reference fluids
c  11-08-06 EWL, load TCX and ETA coefficients into negative arrays
c  12-26-06 MLH, add VS4 model
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
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)
      parameter (nrf0=n0)    !lower limit for transport ref fluid arrays
      character*1 htab,hnull
      character*1 hstar
      character*3 heqn,heta,htcx,hetak,htcxk
      character*3 hetemp,httemp
      character*3 hpheq,heos,hmxeos,hmodcp
      character*3 hsten,hstenk
      character*3 hdiel,hdielk
      character*3 hmdeta,hmdtcx,hmaux
      character*3 hvs,htc
      character*3 hflag,htype
      character*3 hetacr,htcxcr,htcxcrecs
      character*12 hcasn(n0:nx)
      character*255 hfile(n0:nx),href
      character*255 herr,herr1
      character*251 hcite
      character*255 hieos,hicp0,hieta,hietac,hitcx,hitcxc,histn,
     &              hidiel,himelt,hisubl,hips,hipl,hidl,hidv
      logical leta,ltcx
c
      common /NCOMP/ ncomp,ic
      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 /CREMOD2/htcxcrecs(nrf0:nx)
      common /STNMOD/ hsten,hstenk(n0:nx)
      common /DEMOD/ hdiel,hdielk(n0:nx)
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
      common /HCHAR/ htab,hnull
      common /CITE/  hieos(n0:nx),hicp0(n0:nx),
     &               hieta(nrf0:nx),hietac(nrf0:nx),
     &               hitcx(nrf0:nx),hitcxc(nrf0:nx),histn(n0:nx)
      common /CITE2/ hidiel(n0:nx),himelt(n0:nx),hisubl(n0:nx),
     &               hips(n0:nx),hipl(n0:nx),hidl(n0:nx),hidv(n0:nx)
c
      ierrtc=0
      ierrvs=0
      ierrst=0
      href=' '
      ier=0
      ierr=0
      ierr1=0
      herr=' '
      herr1=' '
      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)
        ierr1=0
        hmaux=hmodcp(i)    !possible use in error message
        if (i.ge.nrf0) then
c  allow calls to transport ref. fluids, i.lt.-ncmax
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  accommodate the multiple reference fluids in the ECS-thermo model)
          if (htype(1:2).eq.'CI') hmaux=htype       !possible use in
          if (htype(1:2).eq.'TK') hmaux=htcxcr(i)   !error message
          if (htype.eq.'CI1' .and. hmdeta(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)
          elseif (htype.eq.'CI1' .and. hmdtcx(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. hmdeta(i).eq.htype) then
*           write (*,*) ' SETUP--about to call SETCI2'
            call SETCI2 (nread,i,hcasn(i),ierr1,herr1)
          else if (htype.eq.'CI2' .and. hmdtcx(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)
            hitcxc(i)=htype//hcite//hnull
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)
            hitcxc(i)=htype//hcite//hnull
          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)
            hitcxc(i)=htype//hcite//hnull
          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)
            hitcxc(i)=htype//hcite//hnull
          else if (htype.eq.'TK6' .and. htcxcrecs(i).eq.htype) then
*           write (*,*) ' SETUP--about to call SETTK6'
            call SETTK6 (nread,i,hcasn(i),ierr1,herr1)
            hitcxc(i)=htype//hcite//hnull
          end if
        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
c  transport property models
c
c  do not call transport models if component number is 0
      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 model 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)'
              if (hieta(i)(1:3).eq.'NUL') hieta(i)=htype//hcite//hnull
              if (hitcx(i)(1:3).eq.'NUL') hitcx(i)=htype//hcite//hnull
              call SETTRN (nread,i,hcasn(i),href,heqn,hvs,htc,ier,herr1)
              if (ier.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' .or.
     &      hflag.eq.'VIS' .or. hflag.eq.'vis') 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.ne.0 .and. hetak(i).eq.'NBS') then
              hetak(i)=htype
            end if
c           if (i.lt.-ncmax .and. hetak(i).eq.'NBS') hetak(i)=htype
          end if
c
c  viscosity model "0" found in .fld file
          if (htype.eq.'VS0') then
            if (hetak(i).eq.'VS0') then
*             write (*,*) ' SETUP--about to call SETVS0'
              call SETVS0 (nread,i,hcasn(i),ierrvs,herr1)
              leta=.true.   !flag that fluid-specific vis model loaded
              hieta(i)=htype//hcite//hnull
            end if
c  viscosity model "1" found in .fld file
          elseif (htype.eq.'VS1') then
            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
              hieta(i)=htype//hcite//hnull
            end if
c  viscosity model "2" found in .fld file
          else if (htype.eq.'VS2') then
            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
              hieta(i)=htype//hcite//hnull
            end if
c  viscosity model "3" found in .fld file
          else if (htype.eq.'VS3') then
            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
              hieta(i)=htype//hcite//hnull
            end if
c  viscosity model "4" found in .fld file
          else if (htype.eq.'VS4') then
            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
              hieta(i)=htype//hcite//hnull
            end if
c  viscosity model "6" found in .fld file
          else if (htype.eq.'VS6') then
            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
              hieta(i)=htype//hcite//hnull
            end if
c  ignore the NUL model
          else if (htype.eq.'NUL') then
          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.ne.0 .and. htcxk(i).eq.'NBS') then
              htcxk(i)=htype
            end if
c           if (i.lt.-ncmax .and. htcxk(i).eq.'NBS') htcxk(i)=htype
          end if
c
c  thermal conductivity model "0" found in .fld file
          if (htype.eq.'TC0') then
            if (htcxk(i).eq.'TC0') then
*             write (*,*) ' SETUP--about to call SETTC0'
              call SETTC0 (nread,i,hcasn(i),ierrtc,herr1)
              ltcx=.true.   !flag that fluid-specific t.c. model loaded
              hitcx(i)=htype//hcite//hnull
            end if
c  thermal conductivity model "1" found in .fld file
          elseif (htype.eq.'TC1') then
            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
              hitcx(i)=htype//hcite//hnull
            end if
c  thermal conductivity model "2" found in .fld file
          else if (htype.eq.'TC2') then
            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
              hitcx(i)=htype//hcite//hnull
            end if
c  thermal conductivity model "3" found in .fld file
          else if (htype.eq.'TC3') then
            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
              hitcx(i)=htype//hcite//hnull
            end if
c  thermal conductivity model "6" found in .fld file
          else if (htype.eq.'TC6') then
            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
              hitcx(i)=htype//hcite//hnull
            end if
c  ignore the NUL model
          else if (htype.eq.'NUL') then
          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
c  surface tension model "1" found in .fld file
          if (htype.eq.'ST1') then
            if (hstenk(i).eq.'ST1') then
*             write (*,*) ' SETUP--about to call SETST1'
              call SETST1 (nread,i,hcasn(i),ierrst,herr1)
              histn(i)=htype//hcite//hnull
            end if
c  surface tension model "2" found in .fld file
c         else if (htype.eq.'ST2') then
c           if (hstenk(i).eq.'ST2') then
c             write (*,*) ' SETUP--about to call SETST2'
c             call SETST2 (nread,i,hcasn(i),ierrst,herr1)
c             histn(i)=htype//hcite//hnull
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
c  pure fluid dielectric constant models
c
c  do not call dielectric constant models if component number less than 0
      if (i.ge.nrf0) then
        if (hflag.eq.'DE ' .or. hflag.eq.'de ') then
          if (hstar.eq.'#') then
c  special case for NIST-recommended model
            if (ncomp.eq.1 .and. hdiel.eq.'NBS') then
              hdiel=htype
            end if
            if (i.ge.1 .and. hdielk(i).eq.'NBS') then
              hdielk(i)=htype
            end if
          end if
c
          if (htype(1:2).eq.'DE') then
c  dielectric constant model found in .fld file
            if (hdielk(i)(1:2).eq.'DE') then
*             write (*,*) ' SETUP--about to call SETDE'
              call SETDE (nread,i,hcasn(i),ierrst,herr1)
              hidiel(i)=htype//hcite//hnull
            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) hdielk(i),i,herr1(1:182),hnull
          call ERRMSG (ierr,herr)
        end if
      end if          !end of dielectric constant block
      RETURN
 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)
      end                                             !subroutine SETUP2
c
c ======================================================================
c
      subroutine GETMOD (icomp,htype,hcode,hcite)
c
c  retrieve citation information for the property models used
c
c  inputs:
c    icomp--pointer specifying component number
c           zero and negative values are used for ECS reference fluid(s)
c    htype--flag indicating which model is to be retrieved [character*3]
c           'EOS':  equation of state for thermodynamic properties
c           'CP0':  ideal part of EOS (e.g. ideal-gas heat capacity)
c           'ETA':  viscosity
c           'VSK':  viscosity critical enhancement
c           'TCX':  thermal conductivity
c           'TKK':  thermal conductivity critical enhancement
c           'STN':  surface tension
c           'DE ':  dielectric constant
c           'MLT':  melting line (freezing line, actually)
c           'SBL':  sublimation line
c           'PS ':  vapor pressure equation
c           'DL ':  saturated liquid density equation
c           'DV ':  saturated vapor density equation
c  outputs:
c    hcode--component model used for property specified in htype
c
c           some possibilities for thermodynamic properties:
c           'FEQ':  Helmholtz free energy model
c           'BWR':  pure fluid modified Benedict-Webb-Rubin (MBWR)
c           'ECS':  pure fluid thermo extended corresponding states
c
c           some possibilities for viscosity:
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           'VS4':  generalized friction theory of Quinones-Cisneros and Dieters
c
c           some possibilities for thermal conductivity:
c           'ECS':  extended corresponding states (all fluids)
c           'TC1':  the 'composite' model for R134a, R152a, etc.
c           'TC2':  Younglove-Ely model for hydrocarbons
c
c           some possibilities for surface tension:
c           'ST1':  surface tension as f(tau); tau = 1 - T/Tc
c
c    hcite--component model used for property specified in htype;
c           the first 3 characters repeat the model designation of hcode
c           and the remaining are the citation for the source
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  01-26-00  MM, original version
c  01-23-02 EWL, split common block CITE into 2 pieces for some compilers
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: GETMOD
c     dll_export GETMOD
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)
      parameter (nrf0=n0)    !lower limit for transport ref fluid arrays
      character*3 htype,hcode
      character*255 hcite
      character*255 hieos,hicp0,hieta,hietac,hitcx,hitcxc,histn,
     &              hidiel,himelt,hisubl,hips,hipl,hidl,hidv
      common /CITE/  hieos(n0:nx),hicp0(n0:nx),
     &               hieta(nrf0:nx),hietac(nrf0:nx),
     &               hitcx(nrf0:nx),hitcxc(nrf0:nx),histn(n0:nx)
      common /CITE2/ hidiel(n0:nx),himelt(n0:nx),hisubl(n0:nx),
     &               hips(n0:nx),hipl(n0:nx),hidl(n0:nx),hidv(n0:nx)
c
      if (htype.eq.'EOS' .or. htype.eq.'eos') then
c  equation of state specification
        hcite=hieos(icomp)
      else if (htype.eq.'CP0' .or. htype.eq.'cp0' .or.
     &         htype.eq.'Cp0' .or. htype.eq.'cpo') then
c  ideal-gas heat capacity
        hcite=hicp0(icomp)
      else if (htype.eq.'ETA' .or. htype.eq.'eta' .or.
     &         htype.eq.'VIS' .or. htype.eq.'vis') then
c  viscosity
        hcite=hieta(icomp)
      else if (htype.eq.'VSK' .or. htype.eq.'vsk') then
c  viscosity critical enhancement
c  note:   as of 01/2000, no models employ a viscosity critical enhancement
        hcite=hietac(icomp)
      else if (htype.eq.'TCX' .or. htype.eq.'tcx') then
c  thermal conductivity
        hcite=hitcx(icomp)
      else if (htype.eq.'TKK' .or. htype.eq.'tkk') then
c  thermal conductivity critical enhancement
        hcite=hitcxc(icomp)
      else if (htype.eq.'STN' .or. htype.eq.'stn') then
c  surface tension
        hcite=histn(icomp)
      else if (htype.eq.'DE ' .or. htype.eq.'de ') then
c  dielectric constant
        hcite=hidiel(icomp)
      else if (htype.eq.'MLT' .or. htype.eq.'mlt') then
c  melting line
        hcite=himelt(icomp)
      else if (htype.eq.'SBL' .or. htype.eq.'sbl') then
c  sublimation line
        hcite=hisubl(icomp)
      else if (htype.eq.'PS ' .or. htype.eq.'ps ') then
c  vapor pressure equation
        hcite=hips(icomp)
      else if (htype.eq.'DL ' .or. htype.eq.'dl ') then
c  saturated liquid density equation
        hcite=hidl(icomp)
      else if (htype.eq.'DV ' .or. htype.eq.'dV ') then
c  saturated vapor density equation
        hcite=hidv(icomp)
      end if
      hcode=hcite(1:3)
c
      RETURN
c
      end                                             !subroutine GETMOD
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 Lemmon-Jacobsen model (LJ1):
c             fij(1) = zeta
c             fij(2) = xi
c             fij(3) = Fpq
c             fij(4) = beta
c             fij(5) = gamma
c             fij(6) = 'not used'
c    hfmix--file name [character*255] 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 switch to allow access by DLL
c  02-06-98  MM, initial use of icount uninitialized, use nbin instead
c  03-11-99 EWL, reset lsatt and lsatp to .false. when called
c  02-14-01 EWL, change dimension of amix from 3 to 4
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c  11-01-05 EWL, remove open statement, etc., to check for same model
c  10-23-06 EWL, readd open statement
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: SETKTV
c     dll_export SETKTV
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)
      parameter (nbrule=20)    !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
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      character*1 htab,hnull
      character*3 hmodij,hmodmx,htype
      character*3 hpheq,heos,hmxeos,hmodcp
      character*12 hname
      character*255 hfmix
      character*255 herr
c  next 3 declarations associated with /MXINFO/
      character*8 hbpar
      character*255 hmfile,hrule
      character*255 hbin
c
      dimension fij(nmxpar)
      logical lij(ncmax,ncmax)
c
      common /NCOMP/ ncomp,ic
      common /HCHAR/ 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
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
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),hrule(nbrule),hbpar(nbrule,nmxpar)
      common /MXINF1/ hbin(nbin)
      common /MXRULE/ nrule
c  common to pass flags associated with reference state
      common /IRFSAV/ ixfsav,ksetrf
      common /FSHSAV/ lsatt,lsatp
c
      lsatt=.false.
      lsatp=.false.
      ierr=0
      herr=' '
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 i=1,icount
        do i=1,nbin
          hmfile(i)=hmfile(0)
        enddo
*       write (*,*) ' SETKTV--reset all parameters to original values'
*       write (*,*) '  i  j  mod     Kt        Kv        Fpq      alpha'
*       do i=i,ncomp
*         do j=1,ncomp
*           write (*,1010) i,j,hmodmx(i,j),(fmix(i,j,k),k=1,nmxpar)
*         enddo
*       enddo
*1010   format (1x,2i3,2x,a3,4f10.5)
        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
        call ERRMSG (ierr,herr)
*       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
        call ERRMSG (ierr,herr)
*       write (*,*) herr
        RETURN
      end if
c
c  match input icomp,jcomp with binary pair number
      icount=0
      ibin=1      !initialize only
      do i=1,ncomp-1
        do 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
        enddo
      enddo
c
      do k=1,nmxpar
        fmix(icomp,jcomp,k)=fij(k)
        fmix(jcomp,icomp,k)=fij(k)     !pair (j,i) is same as (i,j)
      enddo
      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)
c
        nread=12             !logical unit for file reads
        call OPENFL (nread,hfmix,1,ierr,herr)
        if (ierr.ne.0) goto 999
        rewind (nread)
        read (nread,2003) htype
        read (nread,'(i2)') ivrsnm
c
*       write (*,1105) hfmix
*1105   format (/1x,' SETKTV--mixture parameters from file: (',a80,')')
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'
*     i=icomp
*     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  continue
      ierr=111
      write (herr,2111) hfmix(1:80),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*255] 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 Lemmon-Jacobsen model (LJ1):
c             fij(1) = zeta
c             fij(2) = xi
c             fij(3) = Fpq
c             fij(4) = beta
c             fij(5) = gamma
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*255]
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 switch to allow access by DLL
c  02-14-01 EWL, change dimension of amix from 3 to 4
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: GETKTV
c     dll_export GETKTV
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)
      parameter (nbrule=20)    !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
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      character*1 htab,hnull
      character*3 hmodij,hmodmx
      character*3 hpheq,heos,hmxeos,hmodcp
      character*8 hfij(nmxpar)
      character*255 hfmix,hmxrul
      character*255 hbinp
c  next 3 declarations associated with /MXINFO/
      character*8 hbpar
      character*255 hmfile,hrule
      character*255 hbin,herr
c
      dimension fij(nmxpar)
c
      common /NCOMP/ ncomp,ic
      common /HCHAR/ 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
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
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),hrule(nbrule),hbpar(nbrule,nmxpar)
      common /MXINF1/ hbin(nbin)
      common /MXRULE/ nrule
c
      hmodij=hmodmx(icomp,jcomp)
      do k=1,nrule
        if (hmodij.eq.hrule(k)(1:3)) then
          irule=k
          goto 100
        end if
      enddo
*     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
 100  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 i=1,ncomp-1
          do 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  old GUI combines hbinp with 3-letter code, so can use only 251 char
              hbinp=hbin(ibin)(1:251)//hnull
            end if
          enddo
        enddo
      end if
c  file containing mixture parameters
c     write (*,*) ' GETKTV--file for ibin =',ibin,' ',hmfile(ibin)
      hfmix=hmfile(ibin)
      hmxrul=hrule(irule)(1:254)//hnull
      do k=1,nmxpar
        fij(k)=fmix(icomp,jcomp,k)
        hfij(k)=hbpar(irule,k)
      enddo
c
*     write (*,1278) hfmix,hbinp,hmxrul
*1278 format (/1x,' GETKTV  hfile:  ',a80/
*    &         1x,'         hbinp:  ',a255/
*    &         1x,'        hmxrul:  ',a255)
*     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 GETFIJ (hmodij,fij,hfij,hmxrul)
c
c  retrieve parameter info for a specified mixing rule
c
c  This subroutine should not be called until after a call to SETUP.
c
c  inputs:
c   hmodij--mixing rule for the binary pair i,j (e.g. LJ1 or LIN)
c           [character*3]
c  outputs:
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     hfij--description of the binary mixture parameters [character*8
c           array of dimension nmxpar]
c   hmxrul--description of the mixing rule [character*255]
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  01-10-01  MM, original version
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: GETFIJ
c     dll_export GETFIJ
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nbrule=20)    !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 hmodij
      character*8 hfij(nmxpar),hbpar
      character*255 hmxrul,hmfile,hrule,hbin,herr
      dimension fij(nmxpar)
c
      common /HCHAR/ htab,hnull
      common /MXINFO/ hmfile(0:nbin),hrule(nbrule),hbpar(nbrule,nmxpar)
      common /MXINF1/ hbin(nbin)
      common /MXINF2/ bideal(nbrule,nmxpar)
      common /MXRULE/ nrule
c
      do k=1,nrule
        if (hmodij.eq.hrule(k)(1:3)) then
          irule=k
          goto 100
        end if
      enddo
*     write (*,*) ' GETKTV ERROR--mixing rule not found: ',hmodij
      ierr=-117
      write (herr,2117) hnull
 2117 format ('[GETKTV error 117] mixing rule not found',a1)
      call ERRMSG (ierr,herr)
      RETURN
 100  continue
      hmxrul=hrule(irule)(1:254)//hnull
      do k=1,nmxpar
        hfij(k)=hbpar(irule,k)
        fij(k)=bideal(irule,k)
      enddo
      RETURN
c
      end                                             !subroutine GETFIJ
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                      end file setup2.f
c ======================================================================
