 
c Please do not distribute any part of this source.
 
c Copyright (c) 1990 by James A. D. Connolly, Institute for Mineralogy
c & Petrography, Swiss Federal Insitute of Technology, CH-8092 Zurich,
c SWITZERLAND. All rights reserved.
 
      program tertex
c----------------------------------------------------------------------
c                       ************************
c                       *                      *
c                       *  tertex.nov.17.1993  *
c                       *                      *
c                       ************************
c----------------------------------------------------------------------
 
c references:
 
c       connolly jad, kerrick dm (1987) bounds:an algorithm and program
c         for calculating composition diagrams, calphad 11, 1-57.
c       connolly jad (1990) multivariable phase diagrams: an algorithm
c         based on generalized thermodynamics, am j sci 290, p 666-718.
 
c-----------------------------------------------------------------------
c file structure:
 
c vertex does i/o from at least three and as many as nine logical units:
 
c  l.u.n  i/o       purpose
c  -----  ---   --------------------------------------------------------
c   n1     i    input to define computations/options.
c   n2     i    data file containing the names, compositional vectors,
c               and standard state thermodynamic data for potentially
c               stable condensed phases and fluid species.
c   n3     o    optional (io3) printer output.
c   n4     o    optional (io4) data file for graphics applications.
c   n5    i/o   optional (io5) permanent data file used to store stable
c               (pseudo-) univariant and invariant phase assemblages.
c               this is done so that the same numbering scheme can be
c               maintained for equilibria during more than one
c               execution of the program.
c   n6    i/o   temporary work file required to generate data for n4,
c               contains the names and stable conditions of invariant
c               assemblages which are written to n4.
c   n7     i    terminal input for interactive systems.
c   n8     o    terminal output for interactive systems.
c   n9     i    optional data file which contains data on the
c               compositional equations of state and subdivision schemes
c               to be used for solution phases.

c note that on batch systems there is no need to distinguish between
c n1 and n7 for input, and n3 and n8 for output.
c-----------------------------------------------------------------------
c this version of vertex has been modified to be compatible with cdc
c fortran 5 compilers.  memory requirements can be reduced on most
c other compilers by changing the implicit variable type integer to
c integer*2, on ibm compilers integer variables beginning with z must
c be integer*4, i.e., for ibm compilers the current implicit type
c statement 'implicit integer (h-n,z)' should be changed to 'implicit
c integer (z), integer*2 (h-n).  the necessity for double precision real
c variables has not been established, users interested in reducing 
c memory c requirements may wish to experiment with single precision 
c variables.
c-----------------------------------------------------------------------
c functional(1) tolerance(2) constant(3) or zero (4) dependent
c     subprograms:                         (this list is not current)

c simpl3-7 (4)    testit (4)      sreset (2)      asschk (4)
c itestc   (4)    loadit (1)      flipit (2)      wway   (2)
c itest3   (4)    solod  (1,3)    schk   (4)      concrt (2,3)
c miscib   (4)    newass (4)      pchk   (2,4)    univeq (2,3)
c balanc   (4)    uproj  (1,2,3)  grxn   (1)      gphase (1)
c mrk      (1,3)  hsmrk  (1,3)    newrap (1,3)    fug    (1,3)

c units, tolerances, and constants:

c function specific
c routines for fluid properties (mrk, trkmrk, hsmrk) require units of
c kelvins (v(2)) and bars (v(3)) and joules.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)

      character*80 record,ipnms*8,irnms*8,vname*8,gname*78,xname*18,
     *             title*72
c-----------------------------------------------------------------------
c dimensioning parameters:

c parameter statements are repeated in every subroutine that the
c parameters are used in.  this means that if you change a dimensioning
c parameter, you must change it globally.

c  h1 - maximum number of phases stable in any hexary subsystem
c  h2 - maximum number of phases stable in any septary subsystem
c  h3 - maximum number of quinary phases stable in any quinary system
c  h4 - maximum number of hexary phases stable in any hexary system
 
c  h5 - maximum number of component saturation constraints
c  h6 - maximum number of composant phases for a component saturation
c       constraint
c  h7 - maximum number of divariant equilibria
 
c  j1 - maximum number of vertices in any binary subsystem
c  j2 - maximum number of facets in any ternary subsystem
c  j3 - maximum number of facets in any quaternary subsystem
c  j4 - maximum number of facets in any quinary subsystem
c  j5 - maximum number of facets in any hexary subsystem
c  j6 - maximum number of ternary phases stable in any c=3 subsystem
c  j7 - maximum number of quaternary phases stable in any c=4 subsystem
c  j8 - maximum number of pseudo-facets.
c  j9 - maximum number of facets in the highest order sytem, should be
c       greater than or equal the highest parameter among j1 through j5.
 
c  k1 - maximum number of phases.
c  k2 - maximum number of univariant equilibria
c  k3 - maximum number of invariant equilibria
c  k4 - number of parameters in thermodynamic equation of state
c  k5 - maximum number of components in thermodynamic data base
c  k6 - maximum number of thermodynamic components in a calculation
c  k7 - k6+1
c  k8 - k6+2
c  k9 - maximum number of phases with lambda transitions
 
c  l2 - maximum number of independent potentials.
c  l3 - maximum number of saturated components.
c  l4 - maximum number of fluid components
c  l5 - maximum number of phases in any unary subsystem
c  l6 - maximum number of phases in any binary subsystem
c  l7 - maximum number of phases in any ternary subsystem
c  l8 - maximum number of phases in any quaternary subsystem
c  l9 - maximum number of phases in any quinary subsystem
c  m1 - maximum number of binary subsystems
c  m2 - maximum number of ternary subsystems
c  m3 - maximum number of quaternary subsystems
c  m4 - maximum number of quinary subsystems
c  m5 - maximum number of hexary subsystems
c  m6 - maximum number of lambda transitions per phase.
c  m7 - maximum number of parameters to describe lambda trans.
c  m8 - maximum number of parameters to describe T dependent ordering.
c  m9 - maximum number of phases with T-dependent ordering.
 
c note parameters k4-l4 cannot be changed independently of other
c variables.
c-----------------------------------------------------------------------
      parameter (k1=1200,k2=1200,k3=1200,k6=7,l2=5)
 
      common/ cst5  /v(l2),tr,pr,r,ps/ cst9  /vmax(l2),vmin(l2),dv(l2)
     *      / cst29 /vip(l2,k3),ipid(k3,k6+2),ipct/ csta8 /title
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k6+2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst81 /ipchk(k3),ip2(k3),icp2,ifake
     *      / cst80 /jrchk(k2),irchk(k2)
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst35 /ir1(k2),iirct,irpct,irv(k2)
     *      / csta1 /irnms(k2,11)/ cst31 /vn(k2,11),ir(k2,11),irct,ird
     *      / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2)
c----------------------------------------------------------------------- 
c                            -------------------------------------------
c                             initialize counters for multiple
c                             calculations.
      iset = 0
      iend = 0
c                             initialize invariant and univariant
c                             field counters:
      iipct = 0
      irpct = 0
      iirct = 0
      irct = 0

c      call money
c                             and flags:
      do 10 i = 1, k2
10       jrchk(i) = 0
      do 11 i = 1, k3
11       ipchk(i) = 0
c                            -------------------------------------------
c                            -------------------------------------------
c                              open statements for units n1-n6 and n9
c                              are in subroutine input1
 
c                              read input from unit n1 (terminal/disk).
c                              subroutine F also initializes:
c                              equilibrium counters; units n2 n4 and n6;
c                              and the limits for numerical results.
300   goto (200),iend
      call input1 (iend,iset)
c                              eof iend = 1
      iend = 1
      if (iset.gt.1) goto 15
c                              read thermodynamic data on unit n2:
      call input2
c                              read data for solution phases on n9:
      call input9
c                              write header to graphics file (n4):
      goto (15), io4
      if (icopt.eq.1) call header
c                              read permanent equilibria lists on u
c                              unit n5:
15    call input5 (iset)
c                              initialize v(iv4) and v(iv5) (these
c                              variables are constant for each set).
c                              calculations.
      v(iv4) = vmin(iv4)
      v(iv5) = vmin(iv5)
c
      if (icopt.eq.0) then
c                              calculate composition phase diagrams:
         call chmcal
c                              chmcal does not return control to main.
       else
c                              output title page to unit n3:
         call outtit (n3,iset)
c                              initialize sectioning variable (iv3):
         v(iv3) = vmin(iv3)
c                              read stability field assemblages on unit
c                              n1:
         if (icopt.eq.2) call error (32,v(1),k2,'MAIN')
c        if (icopt.eq.2) call instb
c
      end if
c                              calculate phase diagram/stability field
c                              for each level of the sectioning variable
      do 20 i = 1, loop
c                              initialization for each tier (level of
c                              iv3):
      ipct = 0
        goto (215),io3
c
215       if (icopt.eq.2) then
c                              calculate stability field:
c            call stbset
c
           else
c                              calculate phase diagram:
             call newhld 
c
          end if
c                              summarize output for each tier
c                              on units n3 and n4 (n6):
      call outier
c                              increment sectioning variable:
20    v(iv3) = v(iv3)+dv(iv3)
 
      goto (95),io3
      write (n3,1160)
c                              transfer list of invariant point
c                              conditions from the scratch file (n6)
c                              to the graphics file (n4)
95    goto (300),io4
      rewind n6
50    read (n6,1040,end = 300) record
      write (n4,1040) record
      goto 50
c                              read next set of input data: goto 300
 
c                              summarize calculations and output
c                              permanent list file to unit n5:
200   if (icopt.ne.3) call outlst
c                              done:
 
1040  format (a80)
1160  format (/,'-------------------------------------------------'
     *          ,'---------------',/)
 
      end

      subroutine outchm
      end

      subroutine chmcal
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7,l2=5,nb=20,k5=12)
 
      character*8 vname,gname*78,xname*18
 
      common/ cst6  /icomp,istct,iphct,icp/ cst5  /v(l2),tr,pr,r,ps
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst2  /g(k1)/ cst3 /x(k1,k6)/ cst12 /cp(k6,k1)
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth

c-----------------------------------------------------------------------
      hh = 1
 
c      ilbl = 1
c      icl = 0
 
      call outtit (n3,hh)
 
90    read (n7,*) v
      do 10 i =  1, 5
10       if (v(i).eq.-9999.d0) goto 99
c                                  compute phase diagram
      call gall (istct,iphct)
      call bossop

      goto 90
 
99    stop

      end

      subroutine bossop
      call combi7
      end 
 
 
      subroutine input1 (iend,iset)
c-----------------------------------------------------------------------
c input1 reads data from a file on unit n1, this data controls the
c computational options and is modified frequently.
c for some purposes all or part of this data may be read from a terminal
c (assigned to unit n7).
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k2=1200,k3=1200,k5=12,k6=7,h5=7,h6=500)
      parameter (l2=5,h8=100,h9=18,k0=25,nb=20)
 
      character*8 fname*10,blank8,ipnms,exname,n5name*14,
     *            cname*5,blank5*5,rname(k6)*5,blnk10*10,
     *            dname*40,cmpnt*5,vname*8,gname*78,xname*18,
     *            title*72,n2name*14,n3name*14,n4name*14,
     *            n9name*14,tcname*5,yes*1,zname
 
      common/ cst6  /icomp,istct,iphct,icp/ cst79 /isoct
     *      / csta7 /fname(h9) / cst24 /ipot,jv(l2),iv(l2)
     *      / cst81 /ipchk(k3),ip2(k3),icp2,ifake
     *      / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst5  /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(l2)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst38 /iedge
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst39 /ichem
     *      / cst88 /vu(2,k1),jprct,jmct,jmyn/ csta8 /title
     *      / cst36 /exname(h8)/ cst63 /delv(l2)/ debug /jtest,jpot
     *      / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k6+2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(k5)
      common/ cst59 /istbct,idst(10,5),ister(10)
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase
     *      / cst207 /ctrans(k0,k5),ictr(k5),itrans
     *      / cst208 /ifct,idfl/ csta9 /tcname(k5)/ cst102 /ivfl
     *      / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx

      save blank5, blank8, blnk10
 
      data blank5,blank8,blnk10/'     ','        ','          '/
c-----------------------------------------------------------------------
c                             open n1
c                             the default designation for the input
c                             read from n1 is 'in1.dat', if this file
c                             does not exist terminate execution:
1     write (*,1050)
      read (*,1200) n2name
 
      open (n1, file = n2name, iostat = ierr, status = 'old')
      if (ierr.ne.0) then
c                                 system could not find the file
         write (*,1040) n2name
         write (*,*) ' Try again (y/n)?'
         read (*,1020) yes
         if (yes.ne.'Y'.and.yes.ne.'y')
     *                        call error (120,r,n1,n2name)
         goto 1
c                                 try again
      end if
c                             error call for noninteractive version:
c     if (ierr.ne.0) call error (120,r,n2,n2name)
 
c                            begin reading input:
c                            -------------------------------------------
c                             read name of thermodynamic data file
      read (n1,1200,end = 90) n2name
c                             open thermodynamic data file
      open (n2, file = n2name, iostat = ierr, status = 'old')
 
      if (ierr.ne.0) call error (119,r,n2,n2name)
 
      write (*,1170) n2name
c                            -------------------------------------------
c                             set flags to indicate output options for
c                             output to unit n3 (print file, y = 0),
c                             output to unit n4 (graphics file, y = 0),
c                             and i/o unit n5 (list file, y = 0).
      read (n1,*,err=998) io3, io4, io5
c                             open print and graphics files
      read (n1,1200) n3name
 
      read (n1,1200) n4name
 
      read (n1,1200) n5name
 
      if (io3.eq.0) then
         open (n3, file = n3name)
         write (*,1180) n3name
      else
         write (*,1180) 'none requested'
      end if
 
      if (io4.eq.0) then
         open (n4, file = n4name)
         write (*,1190) n4name
         open (n6, status = 'scratch')
      else
         write (*,1190) 'none requested'
      end if
c                             increment calculation counter
      iset = iset + 1
c                             computational option, 0 = composition 
c                             diagram, 1 = ptx phase diagram, 2 = stab-
c                             ility field, 3 = mixed variable diagram
      read (n1,*,err=998) icopt, isec, jpot, jtest, ifull, imsg, io3p
c                             read component transformation flag:
      read (n1,*,err=998) itrans, icmpn
      if (itrans.gt.0) then
c                             read new component definitions:
         do 708 i = 1, itrans
            read (n1,1120) tcname(i), ictr(i)
708         read (n1,*) (ctrans(j,i), j = 1, icmpn)
      end if
c                             ichem is a flag which determines
c                             whether chemographies are output
c                             when evaluated by the subroutines
c                             vertex and check. (yes = 0,no = 1).
c     read (n1,*,err=998) ichem
      ichem = 0
c                             iedge is a variable which indicates
c     12222222222222222222    which edges of a phase diagram
c     1                  3    coordinate frame are to be searched
c  i  1                  3    for univariant equilibria.
c  v  1                  3    in this version of vertex (w/newhld)
c  1  1                  3    the traverses are done conterclockwise
c  1  1                  3    as indicated. the traverses must be done
c     1                  3    sequentially, and only relatively small
c     44444444444444444443    savings are obtained by not searching all
c              iv2            four edges, so iedge is not input.
c     read (n1,*,err=998) iedge
c                             iedge = 1 do traverse 1
c                             iedge = 2 do traverses 1 and 2, etc
      iedge = 4
c                             read solid pressure (bars), if solid
c                             pressure equals fluid pressure set ps =
c                             0.00.
c     read (n1,*,err=998) ps
      ps = 0.0d0
c                             read thermodynamic components (8(a5,1x)
      read (n1,1150) rname
c                             count (icp) and save names (cname)
      icp = 0
      do 300 i = 1, k6
         if (rname(i).eq.blank5) goto 300
         icp = icp+1
         cname(icp) = rname(i)
c                             encode a graphics names for the
c                             compositional variables, this is kind of
c                             pointless, but it looks good.
         write (xname(icp),1030) 'x(',rname(i),')'
c                             unblank the name
         call unblnk (xname(icp),18)

         if (icp.gt.k6) call error (197,r,k6,'INPUT1')

300   continue
c                             set counter icp1
      icp1 = icp+1
c                             look for saturated components
      read (n1,1150) rname
c                             isat is the saturated component counter,
c                             isyn is a flag = 0 if isat = 0 (no 
c                             saturated components), else isyn = 1
      isyn = 1
      isat = 0
      io2  = 0 
      do 310 i = 1, h5
         if (rname(i).eq.blank5) goto 310
         isat = isat + 1
         isyn = 0
         if (rname(i).eq.'O2') io2 = isat
         cname(icp+isat) = rname(i)
310   continue
 
      icomp = icp + isat
c                             look for fluid components (2(a5,1x))
      read (n1,1150) rname
c                             ifct is the fluid component counter,
c                             ifyn is a flag = 0 if ifct = 0 (no fluid
c                             components) and ifyn = 1 otherwise.
      ifyn = 1
      ifct = 0
      zname = rname(1)
      do 320 i = 1, 2
         if (rname(i).eq.blank5) goto 320
         ifct = ifct + 1
         ifyn = 0
         cname(icomp+ifct) = rname(i)
320   continue
      icomp = icomp + ifct
c                             look for mobile components (2(a5,1x))
      read (n1,1150) rname
c                             jmct is the mobile component counter,
c                             jmyn is a flag = 0 if jmct = 0 (no fluid
c                             components) and iuyn = 1 otherwise.
      jmyn = 1
      jmct = 0
      do 350 i = 1, 2
         vname(3+i) = blank8
         if (rname(i).eq.blank5) goto 350
         jmct = jmct+1
         jmyn = 0
         cname(icomp+jmct) = rname(i)
c                             encode a text and graphics code name for
c                             each mobile component:
         write (vname(3+jmct),1010) 'u(',rname(i),')'
         call unblnk (vname(3+jmct),8)
         ig(3+jmct) = 18
c                             change u here to mu when you get the font
c                             code, could filter the blanks out to:
c                             di-3000:
         write (gname(3+jmct),1010) 'u(',rname(i),')'
         call unblnk (gname(3+jmct),78)

c        write (gname(3+jmct),1020) 'u[bsub]',rname(i),'[esub]'
c                             erlanger:
c        write (gname(3+jmct),1020) '&1&<m&_&h/',rname(i),'&+(j/mol)&h*'

350   continue
 
      icomp = icomp+jmct
      iprct = icomp-icp
      jprct = icomp-jmct
c                             read bulk composition stuff
      read (n1,*) iblk, jcby, jcth
      if (iblk.gt.0) then
c                                get the bulk compositions:
         do 355 i = 1, iblk 
355         read (n1,*) (cblk(i,j), j = 1, jcth)
      end if 
c                             look for excluded phases
      iexyn = 1
      ixct = 0
      read (n1,1140) (exname(i),i = 1, h8)
c                             ixct is the excluded phase counter,
c                             iexyn is a flag = 0 if ixct = 0 (no fluid
c                             components) and iexyn = 1 otherwise.
      do 330 i = 1, h8
         if (exname(i).eq.blank8) goto 330
         ixct = ixct+1
         iexyn = 0
         exname(ixct) = exname(i)
330   continue
c                             read the maximum pressure, temper-
c                             ature, xco2, u1, and u2; the minimum
c                             pressure temperature, xco2, u1, and u2;
c                             and the default pressure, temperature,
c                             xco2, and chemical
c                             potential increments use kelvins, bars and
c                             joules as units (if no mobile components
c                             enter two zeroes for each read).
      read (n1,*,err=998) vmax
      read (n1,*,err=998) vmin
      read (n1,*,err=998) dv
c                             read the default indices of the
c                             dependent, independent, and secondary
c                             independent intensive variables, p = 1,
c                             t = 2, and xco2 = 3, respectively.
      read (n1,*,err=998) (iv(i),i = 1, 5)
c                             check to make sure input requests are
c                             consistent:
      if (icopt.ne.0) then
c                             first check iv(1):
         if (iv(1).eq.3.and.ifyn.eq.1) call error (110,r,i,'I')
 
         if (iv(1).eq.3.and.ifct.eq.1) call error (111,r,i,'I')
 
         if (vmin(iv(1)).ge.vmax(iv(1))) call error (112,r,i,'I')
 
         if (dv(iv(1)).eq.0.d0) call error (113,r,i,'INPUT1')
 
         dip = dv(iv(1))/(vmax(iv(1))-vmin(iv(1)))
 
         if (dip.lt.0.01) call warn (114,dip,i,'INPUT1')
 
         if (dip.lt.0.001) call error (115,dip,i,'I')
 
         if (vname(iv(1)).eq.blank8) call error (116,dip,i,'I')
c                             next check sectioning varaible
         dip = vmax(iv(3))-vmin(iv(3))
 
         if (dip.ne.0.d0) then
            if (vname(iv(3)).eq.blank8) call error (117,r,i,'I')
 
            if ((iv(3).eq.3).and.(ifct.eq.1)) call error (110,r,i,'I')
 
            dip = dip/dv(iv(3))+1.d0
 
            if (dip.gt.10.d0) call error (118,r,i,'INPUT1')
 
         end if
      end if
c                             now check iv(2):
      if (icopt.eq.1.or.icopt.eq.2) then
         if (iv(2).eq.3.and.ifyn.eq.1) call error (110,r,i,'I')
 
         if (iv(2).eq.3.and.ifct.eq.1) call error (111,r,i,'I')
 
         if (vmin(iv(2)).ge.vmax(iv(2))) call error (112,r,i,'I')
 
         if (dv(iv(2)).eq.0.d0) call error (113,r,i,'I')
 
         dip = dv(iv(2))/(vmax(iv(2))-vmin(iv(2)))
 
         if (dip.lt.0.01d0) call warn (114,r,i,'INPUT1')
 
         if (dip.lt.0.001d0) call error (115,r,i,'I')
 
         if (vname(iv(2)).eq.blank8) call error (116,r,i,'I')
 
      end if
c                             if a chemical potential is specified as an
c                             independent variable (iv(1-3)), check if
c                             the variable is defined:
      kct = 0
      do 370 i = 1, 3
370      if (iv(i).gt.3) kct  =  kct + 1
c                             write a warning:
c      if (jmct.lt.kct) call warn (109,r,i,'INPUT1')
c                             identify the variable used to determine
c                             which phases lie on the left hand side
c                             of a reaction equation.
      if (icopt.eq.3) then

         ivfl = iv(1)

      else if (iv(1).eq.2.or.iv(2).eq.2) then
c                             choose T
         ivfl = 2
      else if (iv(1).eq.1.or.iv(2).eq.1) then
c                             no T, so choose P
         ivfl = 1
      else
c                             no P or T, choose independent V
         if (iv(2).ne.3) then
            ivfl = iv(2)
         else
            ivfl = iv(1)
         end if
      end if
c                             ok, now find out which variables are
c                             dummies and story the indexes of the
c                             non-dummy variables in jv.
      ipot = 0
      do 10 i = 1, 5
c                             variables v(1) (p) and v(2) (t) are
c                             never dummies:
      if (iv(i).eq.1.or.iv(i).eq.2) then
         ipot = ipot+1
         jv(ipot) = iv(i)
c                             variable v(3) is a dummy if ifyn = 1:
      else if ((iv(i).eq.3).and.(ifyn.eq.0)) then
         ipot = ipot+1
         jv(ipot) = iv(i)
c                             variables v(4) and v(4) are dummies if
c                             imyn = 1:
      else
         goto (10), jmyn
         if (iv(i).eq.4) then
            ipot = ipot+1
            jv(ipot) = iv(i)
         else if ((iv(i).eq.5).and.(jmct.eq.2)) then
            ipot = ipot+1
            jv(ipot) = iv(i)
         end if
      end if
10    continue
c                             read data base code for the calculation.
c                             allowed codes are input by the user in the
c                             the data file (21).  see read's below.
      read (n1,*,err=998) idbase
c                             read code for choice of fluid equation
c                             of state from terminal. 
      read (n1,*,err=998) ifug
      if (ifug.ge.7.and.ifug.le.12.and.ifug.ne.9.or.
     *    ifug.eq.14.or.ifug.eq.16.or.ifug.eq.17) 
     *                  read (n1,*,err=998) ibuf,hu,dlnfo2,elag
c                             look for solution phases:
      io9 = 1
      isoct = 0
      read (n1,1000) fname
c                             isoct is the solution phase counter,
c                             io9 is a flag = 0 if ixct = 0 (no fluid
c                             components) and io9 = 1 otherwise.
      do 340 i = 1, h9
         if (fname(i).eq.blnk10) goto 340
         isoct = isoct + 1
         io9 = 0
         fname(isoct) = fname(i)
340   continue
c                             read name of solution model file
      read (n1,1200) n9name
      if (io9.eq.0) then
c                             open thermodynamic data file
         open (n9,file = n9name,iostat = ierr,status = 'old')
 
         if (ierr.ne.0) call error (121,r,n9,n9name)
 
         write (*,1210) n9name
 
       else
 
         write (*,1210) 'none requested'
 
       end if
c                             isudo is a flag which determines if
c                             pseudo-univariant equilibria are defined
c                             for icopt .eq. 1, 2, or 3.
      read (n1,*,err=998) isudo      
c                             read counter for stability field assemblag
      if (icopt.eq.2) read (n1,*,err=998) istbct
c                             read title for the calculation:
      read (n1,1160) title
c                             compute interval range of variables
      do 17 i = 1, 5
17       delv(i) = vmax(i)-vmin(i)
c                             calculate number of sections:
      loop = idint (delv(iv(3))/dv(iv(3))) + 1 
c                             number of phases in invariant equilibria
      icp2 = icp+2
c                              set convergence criteria for routine
c                              univeq:
      call concrt
c
      goto (99),iset
c                              rewind thermodynamic data file:
      rewind n2
c                              why is n6 rewound?
      if (io4.eq.0) rewind n6
c                              rewind solution data file (n9)
      if (io9.eq.0) rewind n9
c
      goto 99
 
90    iend = 1

      goto 99

998   call error (27,r,i,n2name)
 
1000  format (6(a10,1x))
1010  format (a2,a5,a1)
1020  format (a1)
1030  format (a2,a5,a6)
1040  format (' **ERROR VER191** FOPEN file ',A14,
     *        ' could not be opened')
1050  format (' Enter the name of the computational option file ',
     *       '(i.e. the',/,' file created with BUILD)',
     *       ' < 15 characters, left justified: ')
1120  format (a5,1x,i2)
1140  format (6(a8,1x))
1150  format (8(a5,1x))
1160  format (a72)
1170  format (/,' Reading thermodynamic data from file: ',a14)
1180  format (/,' Writing print output to file: ',a14)
1190  format (/,' Writing plot output to file: ',a14)
1200  format (a14)
1210  format (/,' Reading solution models from file: ',a14,/)
 
99    end

      subroutine input2
c----------------------------------------------------------------------
c                     input2
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k2=1200,k3=1200,k4=16,k5=12,k6=7,h5=7,h6=500)
      parameter (l5=400,l6=400,l7=400,k0=25)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h8=100,j9=1200)
 
      character*5 cname, cmpnt, dname*40, zname
      character*8 name, exname

      double precision twt(k5)
 
      integer icout(k0), isub(6,7)
 
      common/ cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst1  /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst88 /vu(2,k1),jprct,jmct,jmyn/ cst208 /ifct,idfl
     *      / cst3  /x(k1,k6)/ cst20 /comps(k1,k5)
     *      / cst13 /id1(k6,l5),id2(m1,l6)/ cst209 /zname
     *      / cst22 /id3(m2,l7)/ cst6 /icomp,istct,iphct,icp
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst5  /p,t,xco2,u1,u2,tr,pr,r,ps/ cst113 /ifproj
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s/ cst36 /exname(h8)
     *      / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(k5)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2/ cst18 /icct(k1)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase
     *      / cst60 /iasmbl(j9),ipoint,imyn/ csta6 /name
     *      / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
c                               isub(i,j) contains the number of 
c                               subcompositions of order j for a system
c                               with order i
      save isub
 
      data isub /6*0, 1,5*0, 3,1,4*0, 6,4,1,3*0, 10,10,5,1,0,0,
     *           15,20,15,6,1,0,  21,35,35,21, 7,1/
c-----------------------------------------------------------------------
c                               initialization for each data set
 
c                               for k1 phases:
      do 1 i = 1, k1
         do 11 j = 1, 2
            vu(j,i) = 0.
11          vf(j,i) = 0.
         do 12 j = 1, h5
12          vs(j,i) = 0.
1           ikp(i) = 0
c                               other counters and flags:
      do 13 i = 1, h5
13       isct(i) = 0
c                               counters for bounds
      do 7 i = 1, m1
7        i2ct(i) = 0
      do 8 i = 1, m2
8        i3ct(i) = 0
      do 9 i = 1, m3
9        i4ct(i) = 0
      do 16 i = 1, m4
16       i5ct(i) = 0
      do 14 i = 1, m5
14       i6ct(i) = 0
      do 15 i = 1, k6
15       i1ct(i) = 0
 
      i7ct = 0
      iphct = 0
      idfl = 0
c                               read data base header
      call topn2 (0,icmpn)
c                               general input data for main program
      goto (6),ifyn
 
      iphct = 2
c                               identify nonzero components.
c                               initialize icout(i) = 0
6     do 10 i = 1, icmpn
10       icout(i) = 0
c
      do 20 i = 1, icomp
         im = 0
         do 30 j = 1, icmpn
            if (cname(i).ne.cmpnt(j)) goto 30
            twt(i) = atwt(j)
            ic(i) = j
            icout(j) = 1
            if (j.eq.idh2o) then
               iff(1) = i
               idfl = idfl + 1
            else if (j.eq.idco2) then
               iff(2) = i
               idfl = idfl + 1
            end if 
            im = 1
30       continue
c                               write error message if a component
c                               was not found:
         if (im.ne.0) goto 20
         write (n3,1230) cname(i), (cmpnt(zk), zk = 1, icmpn)
         write (n3,1240)
         stop
 
20    continue
cjc                            check if a fluid is in the thermodynamic
cjc                            composition space:
      ifproj = 0

      if (ifct.eq.0) then 
         do 100 i = 1, icp
            if (cname(1).eq.cmpnt(idh2o)) then
               if (cname(2).eq.cmpnt(idco2)) then
c                              both components present:
                  ifproj = 3
               else 
c                              only water present
                  ifproj = 1
               end if 
            else if (cname(1).eq.cmpnt(idco2)) then
c                              only co2 present
                  ifproj = 2
            end if
100      continue 
      end if 
c                              this segment is to check if
c                              a possible saturated phase component
c                              has been made a mobile component,
c                              if there is also a saturated phase
c                              component idfl is the identity of the
c                              mobile component otherwise idfl = 0.
      if (ifct.eq.1.and.idfl.eq.2) then
         if (zname.eq.cmpnt(idh2o)) then
            idfl = 1
         else if (zname.eq.cmpnt(idco2)) then
            idfl = 2
         end if
      else 
         idfl = 0
      end if
c                               load atwts in updated order
      do 25 i = 1, icomp
25       atwt(i) = twt(i)
c                               begin first read loop for data on
c                               constrained components.
      if (isyn.ne.0.and.ifyn.ne.0) goto 40
 
      ifer = 0
 
50    call getphi (name,ibase,icmpn,*59)
 
      call chkphi (0,name,icout,ibase,icmpn,*50)
 
      if (ifyn.eq.0) then
c                               check for fluid species data
         if (name.eq.cmpnt(idco2)) then
            j = 2
            ifer = ifer + 1
         else if (name.eq.cmpnt(idh2o)) then
            j = 1
            ifer = ifer + 1
         else
            goto 70
         end if
 
         call loadit (j)
 
         goto 50
      end if
 
70    goto (50),isyn
c                               check for saturated phases:
 
c                               reject the phase if it contains
c                               a component not specified in the
c                               input:
      do 470 j = 1, icmpn
         if (icout(j).eq.0.and.comp(j).ne.0.) then
            goto 50
         else if (comp(j).lt.0.) then
            call warn (13,r,j,name)
            goto 50
         end if
470   continue 
c                               reject the phase if it contains
c                               a thermodynamic component:
      do 471  j = 1, icp
471      if (comp(ic(j)).ne.0.) goto 50
c                               now load the phase if it has
c                               the saturated component idc:
      do 472 j = isat, 1, -1
         idc = icp + j
         if (comp(ic(idc)).ne.0.d0) then
            isct(j) = isct(j) + 1
            if (isct(j).gt.h6) call error (17,r,h6,'INPUT2')
            iphct = iphct + 1
            if (iphct.gt.k1) call error (180,r,k1,'INPUT2')
            ids(j,isct(j)) = iphct
            call loadit (iphct)
            goto 50
         end if
472   continue
 
      goto 50
c                               check that there is data for
c                               every fluid component.
59    if (ifyn.eq.0.and.ifer.ne.2) call error (36,r,i,'INPUT2')
c                               check that there is one phase
c                               for each saturation constraint
40    if (isyn.eq.0) then
         do 479 i = 1, isat
479        if (isct(i).lt.1) call error (15,r,i,cname(icp+i))
      end if
c                               read data for the remaining
c                               phases of appropriate composition.
      istct = iphct+1
c                               get subsystem counters from isub
      i2s = isub(1, icp)
      i3s = isub(2, icp)
      i4s = isub(3, icp)
      i5s = isub(4, icp)
      i6s = isub(5, icp)
      i7s = isub(6, icp)
 
      rewind n2
c                               read till end of header
      call eohead (n2)
c                               begin second read loop:
220   call getphi (name,ibase,icmpn,*999)
c                               check if valid phase:
      call chkphi (1,name,icout,ibase,icmpn,*220)
c                               acceptable data, count the phase:
      iphct = iphct + 1
c                               calculate normalized composition:
      tot = 0.
 
      do 310 i = 1, icp
         if (comp(ic(i)).lt.0.) then
c                               reject phase with negative composition:
            call warn (27,r,i,name)
            iphct = iphct - 1
            goto 220
         end if
310      tot = tot+comp(ic(i))
 
      do 320 i = 1, icp
320      x(iphct,i) = comp(ic(i)) / tot
 
      if (iphct.gt.k1) call error (180,r,k1,'INPUT2')
c                               store thermodynamic parameters:
      call loadit (iphct)
c                               next phase:
      goto 220
  
1230  format (' **error ver013** ',a5,' is an incorrect component'
     *       ,' name, valid names are:',/,12(1x,a5))
1240  format (' check for upper/lower case matches or extra blanks'
     *       ,/)
999   end
 
      subroutine input5 (iset)
c----------------------------------------------------------------------
c input5 reads data from a file on unit n5. this file contains lists
c of univariant and invariant equilibria and is maintained (optionally)
c so that the numbering schemes used for the equilibria is consistent
c for multiple executions of the program. if the requested list file
c is empty, the program creates a list for subsequent executions.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k2=1200,k3=1200,k6=7)
 
      character*8 ipnms, irnms
 
      common/ cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k6+2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11)
     *      / cst81 /ipchk(k3),ip2(k3),icp2,ifake
c-----------------------------------------------------------------------
      if (iset.gt.1) goto 99
c                             option requested? if yes io5 = 0
      goto (99),io5
c                             read invariant assemblages.
      read (n5,1070,end = 99) iipct
      if (iipct.eq.0) goto 10
      read (n5,1060) ((ipnms(zi,zj),zj = 1, icp2),zi = 1, iipct)
c                             read univariant assemblages
10    read (n5,*) iirct
      if (iirct.eq.0) goto 99
      do 20 i = 1, iirct
         read (n5,*) ip
         irv(i) = ip
20       read (n5,1060) (irnms(i,zj),zj = 1, ip)
      goto 99
 
1060  format (8(a8,1x))
1070  format (bn,i80)
 
99    end
 
      subroutine input9 
c-----------------------------------------------------------------------
c given a list of solution phase names (fname(h9)) input9 searches a
c data file (on unit n9) for the relevant data and subdivides the
c solutions into pseudo-compounds.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
  
      parameter (k1=1200,k2=1200,k3=1200,k4=16,
     *           j9=1200,h5=7,h9=18,m0=5)
 
      character*10 tname, uname(2)*8, sname(h9), fname, names*8

      double precision z(3,3)
 
      common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3),
     *               iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3),
     *               ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3),
     *               jst(3),jfix,idfc(3),ist(3),isp(3),isite,
     *               ispct(3,4),jnsp(3,4),jsite,iterm,iord
     *      / cst86 /xy(2,k1),y(k1,3,2),ntot,npairs/ cst8  /names(k1)
     *      / cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3),
     *                nsite,nspm1(5),nterm(5,5),
     *                nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0)
     *      / cst60 /iasmbl(j9),ipoint,imyn/ csta7 /fname(h9)
      common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst79 /isoct/ cst6  /icomp,istct,iphct,icp
     *      / cst1  /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar

      save uname

      data uname/' ',' '/
c-----------------------------------------------------------------------
c                               set compound counter ipoint.
      ipoint = iphct
c                               no request for solutions, goto 999:
      goto (999),io9
c                               initialize match flag:
      im = 0
c                               read the solution name
10    call rmodel (tname,istot)

      if (istot.eq.0) goto 99
c                               check the solution model:
      call cmodel (im,idsol,tname,istot,jstot,0,uname(1),uname(2))

      if (jstot.eq.0) goto 10 
c                               
      if (jstot.lt.istot) then
c                               if ok, determine a valid
c                               endmember subset.
         call emodel (tname,istot)
c                               edit model parameters:
         if (istot.eq.0) then
c                               couldn't revise model
            im = im - 1
            goto 10
         end if 

         call nmodel

      end if
c                               save found solutions:
      sname(im) = tname
c                               if nsite ne 0 get "normalization"
c                               constants for entropy model:
      if (nsite.ne.0) call snorm (tname)
c                               generate pseudocompound compositions.
c                               subdiv returns the total
c                               number of pseudocompounds (ipcps) and 
c                               array y, of which element y(h,i,j) is
c                               the site fraction of the jth species on
c                               the ith site of the hth pseudocompound.
      call subdiv (tname)
c                               load default site fractions
c                               into the unused sites of the
c                               site fraction array:
      if (isite.lt.3) then
        do 20 i = isite + 1, 3
           z(i,1) = 1.0
           z(i,2) = 0.0
20         z(i,3) = 0.0 
      end if          
c                               subdiv generates ntot compositions,
c                               generate the compound data for
c                               each solution:
      do 30 h = 1, ntot
c                               load the composition into
c                               a the site fraction array:
         do 40 i = 1, isite
            zt = 0.0
            do 50 j = 1, isp(i) - 1
               z(i,j) = y(h,i,j)
50             zt = zt + z(i,j)
40          z(i,isp(i)) = 1. - zt
c                               generate the pseudocompound:
         call soload (z,idsol,tname)

30    continue

      if (im.eq.isoct) goto 999
c                             next solution
      goto 10 
c                             at least one solution phase referenced
c                             in the input is not included in the
c                             solution phase data file, write warning:
99    call warn (43,zt,isoct-im,'INPUT9')
c                             scan for "killed endmembers"
999   do 60 i = 1, iphct
         if (ikp(i).ge.0) goto 60
c                             eliminate killed endmembers by 
c                             giving them positive G, and near
c                             zero entropies
         thermo(1,i) = 1.d7
         thermo(2,i) = 0.
c                             reset ikp
         ikp(i) = 0
         write (n3,1030) names(i)
         
60    continue 

      if (isoct.gt.0) then 
         write (n3,1020)
         write (n3,1010) (fname(i), i = 1, isoct)
         if (im.gt.0) then 
            write (n3,1000)
            write (n3,1010) (sname(i), i = 1, im)
         else 
            write (n3,1040) 
         end if 
      end if 

1000  format (/,' the following solution models will be',
     *          ' considered:',/)
1010  format (7(2x,a10))
1020  format (/,' Of the requested solution models:',/)
1030  format (' **warning ver044** a solution model has destabilized',
     *        ' the endmember: ',a8,' (iend=2)',/)
1040  format (/,' no models will be considered.',/)

      end
 
      subroutine outgrf
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k2=1200,k3=1200,k6=7)
      parameter (l2=5,l5=400,l6=400,l7=400)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7)
      parameter (j1=50,j2=400,j3=900,j4=1200,j5=2400,j6=200)
      parameter (h4=20,j9=1200,h5=7,h6=500)
 
      character*8 names,vname*8,gname*78,xname*18,title*72
 
      common/ cst3  /x(k1,k6)/ cst96 /id7f(k6,j9),i7fct
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl
     *      / cst5  /v(l2),tr,pr,r,ps/ cst8 /names(k1)
     *      / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2/ csta8 /title
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst22 /id3(m2,l7)/ cst6 /icomp,istct,iphct,icp
      common/ cst13 /id1(k6,l5),id2(m1,l6)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst45 /itbh(3),itt,i3vcth,ib/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst60 /iasmbl(j9),ipoint,imyn/ cst79 /isoct
     *      / cst52 /hcp,h,hd,id,jd,kd,ld,md(3),
     *               itc1,itc2,itc3,itc4,ic(3)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst103 /isec,icopt,ifull,imsg,io3p
c-----------------------------------------------------------------------
 
      if (icp.eq.2) then
       ib = i2vct(1)
      else if (icp.eq.3) then
       ib = i3fct(1)
      else if (icp.eq.4) then
       ib = i4fct(1)
      else if (icp.eq.5) then
       ib = i5fct(1)
      else if (icp.eq.6) then
       ib = i6fct(1)
      else
       ib = i7fct
      end if
c                             number of components, phase counters,
c                             assemblage counter, fluid saturation flag,
c                             component saturation flag, ipot = number
c                             of independent potential variables:
c
      write (n4,1020) icopt
c
      write (n4,1020) icp,istct,iphct,ipoint,ib,ifyn,isyn,ipot,isoct
c                             write graphics code variable names:
      write (n4,1030) (ig(jv(z)),gname(jv(z)),z = 1, ipot)
c                             write values of the independent potentials
      write (n4,1040) (v(jv(z)),z = 1, ipot)
c                             write a blank record as left caption
      write (n4,1050) title
c                             write phase names
      write (n4,1015) (names(zi),zi = 1, iphct)
c                             write phase coordinates
      write (n4,1025) ((x(zi,zj),zj = 2,icp),zi = istct,iphct)
c                             stable configurations, phases are
c                             labelled by the index 'zi' in the
c                             list of phases.
      goto (99,950,960,970,980,985,990),icp
c                             binary
950   write (n4,1020) (id2v(1,zj),zj = 1, ib)
      goto 995
c                             ternary
960   write (n4,1020) ((id3f(1,zi,zj),zj = 1, 3),zi = 1, ib)
      goto 995
c                             quaternary
970   write (n4,1020) ((id4f(1,zi,zj),zj = 1,4),zi = 1, ib)
      goto 995
c                             quinary
980   write (n4,1020) ((id5f(zj,zi,1),zj = 1,5),zi = 1, ib)
      goto 995
c                             hexary
985   write (n4,1020) ((id6f(zj,zi,1),zj = 1,6),zi = 1, ib)
      goto 995
c                             septary
990   write (n4,1020) ((id7f(zj,zi),zj = 1, 7),zi = 1, ib)
c                             write assemblage flags
995   if (icp.gt.2) write (n4,1020) (iasmbl(zj),zj = 1, ib)
      goto (99),isyn
      write (n4,1020) isat
      write (n4,1020) (idss(zi),zi = 1, isat)
 
1015  format (10(a8))
1020  format (20(i3,1x))
1025  format (11(f6.4,1x))
1030  format (i2,a78)
1040  format (5(g12.6,1x))
1050  format (a72)
 
99    end

      subroutine newhld 
c-----------------------------------------------------------------------
c newhld replaces nohold.  newhld is designed to do boundary traverses
c for individual equilibrium configurations and in certain cases,
c specifically for systems with fewer than 4 components, is 
c slower than nohold. in contrast to newhld, nohold simultaneously
c monitored all assemblages on the thermodynamic surface of the system.
c a copy of the nohold code is stored in check.fortran on tape xyw003.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-m,z)

      parameter (k1=1200,k2=1200,k6=7,l2=5,h5=7,nb=20)
      parameter (j9=1200,h7=6000,h6=500,k5=12)

      character*8 vname,gname*78,xname*18,cname*5 

      common/ cst23 /a(10,10),b(11),ipvt(10),idv(7),idv8,idv9,idv10,
     *               iophi,idphi,iiphi,iflg1
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ csta4 /cname(12)
     *      / cst96 /id7f(k6,j9),i7fct/ cst6 /icomp,istct,iphct,icp
     *      / cst5  /v(l2),tr,pr,r,ps/ cst7 /iflag/ cst38 /iedge
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)/ cst39 /ichem
     *      / csta2 /xname(k6),vname(l2),gname(l2)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst80 /jrchk(k2),irchk(k2)/ cst302 /jasmbl(j9)
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst31 /vn(k2,11),ir(k2,11),irct,ird
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst65 /vt(j9),vvt(h7),jok(j9),jnk(h7),idns(k6,h7),idnct
     *      / cst99 /ipart,index,kok(j9),knk(h7)/ cst102 /ivfl
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth
     *      / cst303 /vti(j9),vvti(h7)
c-----------------------------------------------------------------------
c                                  initialization:
      idnct = 0
      ipart = 1
      irend = 0
      irct = 0
 
      if (icopt.eq.3) iedge= 1
      do 70 i = 1, k2
70       irchk(i) = 0
c                                  set starting values for search
      v(iv1) = vmin(iv1)
      v(iv2) = vmin(iv2)

      call gall (istct,iphct)
c                                  determine the stable chemography.
      call bossop
c                                  load assemblages into the array id7f
c                                  with counter i7fct (cst58).
10    goto (15),io3
       write (n3,1160)
       write (n3,1170) vname(iv1),vname(iv2)
c                                  write potential variable sectioning
c                                  constraints:
       if (ipot.gt.2) then
         write (n3,1180)
         write (n3,1190) (vname(jv(z)),v(jv(z)),z=3,ipot)
       end if
c                                  write saturated and buffered
c                                  component names:
       goto (100),isyn
       write (n3,1200) (cname(z+icp),z= 1, isat)
c                                  write dependent extensities blurb
100   write (n3,1210)
      write (n3,1230) vname(ivfl)

      write (n3,1160)

      if (io3p.eq.1) then
         write (n3,1240)
         write (n3,1160)
      end if

c                                  initialize start parms:
15    do 5 i = 1, i7fct
         vt(i) = vmin(iv1)
         vti(i) = vmin(iv2)
         kok(i) = 1
5        jok(i) = 1
      
      write (*,1000) i7fct
c                                  look for reactions involving each
c                                  assemblage.
      do 20 i = 1, i7fct

         if (jcby.eq.1.and.jasmbl(i).eq.0) goto 20

         if (imsg.eq.0) write (*,1110) i,idnct
c                                  load id's into array idv:
         do 30 j = 1, icp
30          idv(j) =id7f(j,i)
c                                  get the lower/upper decomposition
c                                  of the transpose of the concentration
c                                  matrix:
         call pivots (ier)
         goto (9010),ier
c                                  test stability along an edge of the
c                                  diagrams coordinate frame:
         kok(i) = jok(i)
         call search (vt(i),vti(i),jok(i),iste,ivi,ivd,div,ier)

c                                  ier =2 error from search.
c                                  ier =1 stable at all conditions.
         goto (20,20),ier
         call coface (ivd,ivi,div,iste,irend)
c
20    continue
c                                  done, now repeat for new assemblages:
      knct = 1
      ipart = 2
40    jnct = idnct
c
      do 50 i = knct, jnct

         if (imsg.eq.0) write (*,1220) i,idnct-i
c                                  load id's into array idv:
         do 60 j = 1, icp
60          idv(j) =idns(j,i)

         if (iblk.ne.0) then
            call bulkck (idv,jbfg,ier)
            if (jcby.eq.1.and.jbfg.eq.0) goto 50
         end if 
c                                  get the lower/upper decomposition
c                                  of the transpose of the concentration
c                                  matrix:
         call pivots (ier)
         goto (9010),ier
c                                  test stability along an edge of the
c                                  diagrams coordinate frame:
         knk(i) = jnk(i)
         call search (vvt(i),vvti(i),jnk(i),iste,ivi,ivd,div,ier)
c                                  ier =2 error from search.
c                                  ier =1 stable at all conditions.
         goto (50,50),ier
c                                  call coface to delineate all the
c                                  equilibria that are cofacial with
c                                  the simplex defined by the vertices
c                                  idv.
         call coface (ivd,ivi,div,iste,irend)
50    continue

      if (idnct.eq.jnct) goto 90
      knct=jnct + 1
      goto 40
c                                  return for icopt=1
90    if (idnct.eq.h7) call warn (205,r,h7,'CHKASS')
      goto (99),icopt
c                                  for icopt=3 sort and output
c                                  reactions
      call onedim

99    return
c                                  error, a degenerate assemblage in
c                                  pivots due to a programming error.
9010  call error (14,v(1),ier,'NEWHLD')
      
1000  format (' Initial number of',
     *        ' divariant assemblages to be tested is:',i4)
1110  format (' Testing initial assemblage ',i4,', ',i4, 
     *        ' new assemblages identified')
1160  format (/,'-------------------------------------------------'
     *          ,'---------------------',/)
1170  format (' The ',a8,'-',a8,' loci of (pseudo-) univariant fields'
     *        ,' follow:')
1180  format (/,' the fields are subject to the constraint(s):'/)
1190  format (25x,a8,'=',g12.6)
1200  format (/,' These fields are consistent with saturation or ',
     *            'buffering constraints',/,' on the component(s):',
     *            1x,7(a5,2x))
1210  format (/,' NOTE: For each field the values of the dependent ',
     *          ' extensities are output',/,' for the first',
     *          ' equilibrium condition, in general these properties',
     *          ' vary with',/,' the independent potentials.',/)
1220  format (' Testing new divariant assemblage ',i4,', ',i4,
     *        ' assemblages remaining to be tested.')
1230  format (/,' Reaction equations are written such that the high ',
     *          a8,/,' assemblage is on the right of the = sign',/)
1240  format (//,' You have chosen the short form print file. ',/,
     *           ' If you want equilibrium coordinates output select',
     *           /,' the long form when you run BUILD.',//)

      end


      subroutine soload (z,isolct,tname)

      implicit double precision (a-g,o-z),integer (h-n)
  
      parameter (k1=1200,k2=1200,k3=1200,k4=16,k5=12,k6=7,m0=5)
      parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h5=7,h6=500,nb=20)

      character*10 tname, names*8
 
      double precision z(3,3)

      common / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
      common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3),
     *               iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3),
     *               ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3),
     *               jst(3),jfix,idfc(3),ist(3),isp(3),isite,
     *               ispct(3,4),jnsp(3,4),jsite,iterm,iord
      common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3),
     *                nsite,nspm1(5),nterm(5,5),
     *                nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0)
      common/ cst1  /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst202 /tm(m7,m6),td(m8),ilam,idiso,lamin,idsin
     *      / cst203 /therdi(m8,m9),therlm(m7,m6,k9)
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst20 /comps(k1,k5)/ cst8 /names(k1)
     *      / cst204 /ltyp(k1),lmda(k1),idis(k1)
     *      / cst6  /icomp,istct,iphct,icp
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst5  /p,t,xco2,u1,u2,tr,pr,r,ps/ cst12 /cp(k6,k1)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst88 /vu(2,k1),jprct,jmct,jmyn
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth
c                              eliminate end-member compositions
      do 10 i = 1, isp(1)
         call zchk (z(1,i),tname)
         do 20 j = 1, isp(2)
            call zchk (z(2,j),tname)
            do 30  k = 1, isp(3)
               call zchk (z(3,k),tname)
30             if (z(1,i)*z(2,j)*z(3,k).gt.0.999) goto 999
20       continue
10    continue    
c                               eliminate out of range compositions:  
c                               is this necessary?
      do 50 i = 1, isite
         if (imd(i).gt.0) goto 50
         do 60 j = 1, isp(i) - 1
60          if (z(i,j).gt.xmx(i,j).or.z(i,j).lt.xmn(i,j)) goto 999
50    continue
c                               the composition is acceptable.
      iphct = iphct + 1
      jltyp = -1
      if (iphct.gt.k1) call error (180,z(1,1),k1,'PRIME9')
      ikp(iphct) = isolct
      idis(iphct) = 0
      lmda(iphct) = 0
c                               encode a name
      if (isite.eq.1) then
         if (isp(1).eq.2) then
           if (z(1,1).eq.1.0) then 
              write (names(iphct),1050) names(imsol(1,1,1)),100
           else if (z(1,1).ge.0.04.and.z(1,1).le.0.96) then
              write (names(iphct),1030) names(imsol(1,1,1)),
     *                                  idint(1.d2*z(1,1))
           else   
              write (names(iphct),1070) names(imsol(1,1,1)),
     *                                  100.*z(1,1)
           end if 

         else 
c                                ternary solutions
            write (names(iphct),1040) names(imsol(1,1,1)),
     *         idint(z(1,1)*1.d2),names(imsol(2,1,1)),idint(z(1,2)*1.d2)

         end if 

      else if (isite.eq.2.and.isp(1).eq.2.and.isp(2).eq.2) then
 
          int1 = idint (1.d2*z(1,1))
          int2 = idint (1.d2*z(2,1))
          if (int1.gt.99) int1 = 99
          if (int2.gt.99) int2 = 99
 
          write (names(iphct),1020) names(imsol(1,1,1)),int1,int2
 
      else
c                              this is a pretty shitty name
c                              but for the time being your going to
c                              have live with it.
         write (names(iphct),1010) ((idint(1.d1*z(i,j)), 
     *                               j = 1, isp(i) - 1), i = 1, isite)
 
         write (n3,1000) names(iphct),(((z(1,i)*z(2,j)*z(3,k),
     *                      names(imsol(i,j,k)),i=1,isp(1)),
     *                                          j=1,isp(2)),
     *                                          k=1,isp(3))
 
      end if
c                              get blanks out of name:
      call unblnk (names(iphct),8)    
c                              initialize constants:
      smix = 0.0
      esum = 0.
      atw(iphct) = 0.
      vol(iphct) = 0.

      do 130 i = 1, icp
130      cp(i,iphct) = 0.
      do 140 i = 1, k4
140      thermo(i,iphct) = 0.
      do 150 i = 1, icomp
150      comps(iphct,i) = 0.
c                              load constants:
      do 160 i = 1, isp(1)
         if (z(1,i).eq.0.) goto 160
         do 170 j = 1, isp(2)
            if (z(2,j).eq.0.) goto 170
            do 180 k = 1, isp(3)
               if (z(3,k).eq.0.) goto 180

               zp = z(1,i)*z(2,j)*z(3,k)
               id = imsol(i,j,k)
               esum = esum + zp * ecoef(i,j,k)
               atw(iphct) = atw(iphct) + zp * atw(id)
               vol(iphct) = vol(iphct) + zp * vol(id)

               do 190 l = 1, icp
  190             cp(l,iphct) = cp(l,iphct) + zp * cp(l,id)
               do 200 l = 1, icomp
  200             comps(iphct,l) = comps(iphct,l) + zp * comps(id,l)
               do 210 l = 1, k4
  210             thermo(l,iphct) = thermo(l,iphct)
     *                        + zp * thermo(l,id)

c                                t dependent disordering:

               if (idis(id).ne.0) then

                  idsin = idsin + 1

                  idis(iphct) = idsin

                  do 330 l = 1, 7
330                  therdi(l,idsin) = zp * therdi(l,id)
                  do 340 l = 8, 9
340                  therdi(l,idsin) =      therdi(l,id)

               end if
c                                lamda-like transitions:
               ld = lmda(id)

               if (ld.ne.0) then

                  if (jltyp.eq.-1) lamin = lamin + 1

                  jltyp = jltyp + 1

                  ltyp(iphct) = ltyp(id)
                  lmda(iphct) = lamin
                  ld = lmda(id)

                  if (ltyp(id).eq.1) then

c                                ubc-approach:

                     do 350 l = 1, 3
 
                        therlm(1,l,lamin) = therlm(1,l,ld)*zp
                        therlm(2,l,lamin) = therlm(2,l,ld)*zp
                        therlm(3,l,lamin) = therlm(3,l,ld)
                        therlm(4,l,lamin) = therlm(4,l,ld)
                        therlm(5,l,lamin) = therlm(5,l,ld)*zp
                        therlm(6,l,lamin) = therlm(6,l,ld)*zp
                        therlm(7,l,lamin) = therlm(7,l,ld)
                        therlm(8,l,lamin) = therlm(8,l,ld)*zp
350                     therlm(9,l,lamin) = therlm(9,l,ld)*zp
c
                   else if (ltyp(id).lt.8) then

c                                helgeson approach:

                       if (ltyp(id).gt.4.and.ltyp(id).lt.7) then
                          jlam = ltyp(id) - 3
                       else
                          jlam = 1
                       end if

                       do 365 m = 1, jlam
                          therlm(1,m,lamin) = therlm(1,m,ld)
                          therlm(2,m,lamin) = therlm(2,m,ld)
                          do 360 l = 3, 8
360                          therlm(l,m,lamin) = therlm(l,m,ld)*zp
365                    continue

                   else if (ltyp(id).eq.10) then

c                                holland and powell landau model:

                      h = jltyp + 1
                      therlm(1,h,lamin) = therlm(1,1,ld)
                      do 222 m = 2, 5
222                      therlm(m,h,lamin) = therlm(m,1,ld)*zp

                   end if 
                end if
180          continue
170       continue
160    continue

c                              compute ideal configurational negentropy:
       if (nsite.eq.0.and.ist(1).ne.0) then  
          do 220 i = 1, isite 
             do 230 j = 1, isp(i)
230             if (z(i,j).ne.0.0) smix = 
     *                          smix + float(ist(i))*z(i,j)*dlog(z(i,j)) 
220       continue
       else if (nsite.ne.0) then
          call omega (z,dlnw,tname)
          smix = -dlnw + esum
       end if 
c                              save it:
       thermo(2,iphct) = thermo(2,iphct) + r * smix
 
       goto (290),isyn
c                               load stoichiometric coefficients of
c                               saturated components:
       do 300 k = 1, isat
300       vs(k,iphct) = comps(iphct,icp+k)
c                               load stoichiometric coefficients of
c                               mobile components:
290    goto (310),jmyn
       do 320 i = 1, jmct
320       vu(i,iphct) = comps(iphct,jprct+i)
 
310    if (iff(1).ne.0) vf(1,iphct) = comps(iphct,iff(1))
       if (iff(2).ne.0) vf(2,iphct) = comps(iphct,iff(2))
c                              load excess terms:
       if (iterm.ne.0) then
          do 240 i = 1, iterm
             zp = 1.0
             do 250 j = 1, iord
                if (isub(i,j,1).eq.0) goto 250
                zp = zp * z(isub(i,j,1),isub(i,j,2))
250          continue

              do 260 j = 1, 11
260              thermo(j,iphct) = thermo(j,iphct) + zp * wg(i,j)
240      continue
      end if 

1000  format (' icky name ',a8,' = ',4(4(f5.3,1x,a4,1x),/))
1010  format (8i1)
1020  format (a2,i2,'-',i2)
1030  format (a2,i2)
1040  format (a1,i2,a1,i2)
1050  format (a2,i3)
1070  format (a2,f4.1)

999   end
     
      subroutine combi7
c--------------------------------------------------------------------
c combinatorial enumeration
 
c   referenced by: bossop
c   references to: abload,assd7,dgph7,sidc,slopc
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)

      parameter (k6=7,j9=1200,k1=1200)

      character*1 y, names*8

      common/ cst52 /hcp,l,idv(7),i0,ic(7)/ cst96 /id7f(k6,j9),i7fct
     *      / cst55 /s(6),sign,ipc(6),jcp,kcp/ cst302 /jasmbl(j9)
      common/ cst6  /icomp,istct,iphct,icp/ cst8  /names(k1)

c-----------------------------------------------------------------------

      write (*,*) ' enter the assemblage yourself (y/n)?'
      read (*,1000) y
1000  format (a1)
1010  format (1x,6(i4,1x,a8,2x))

      if (y.eq.'y') then
         write (*,*) ' output phase indices (y/n)?'
         read (*,1000) y
         if (y.eq.'y') then
            do 9 i = 1, iphct, 5
               k = i + 5
               if (k.gt.iphct) k = iphct
9              write (*,1010) (j, names(j), j = i, k)
         end if
         write (*,*) ' enter phase indices '
         read (*,*) (idv(i),i = 1,icp)
         call testas (ier) 
         if (ier.ne.0) call error (107,sign,ier,'COMBI7')
         goto 99
      end if

      hcp = icp
      do 200 i = 1, icp
200      ic(i) = i       

      do 10 i1 = 1, iphct
         idv(1) = i1
         if (icp.eq.1) then
            call testas (ier)
            if (ier.eq.0) goto 99
            goto 10
         end if 
         do 20 i2 = i1 + 1, iphct
            idv(2) = i2 
            if (icp.eq.2) then
               call testas (ier)
               if (ier.eq.0) goto 99
               goto 20
            end if 
            do 30 i3 = i2 + 1, iphct
               idv(3) = i3
               if (icp.eq.3) then
                  call testas (ier)
                  if (ier.eq.0) goto 99
                  goto 30
               end if 
               do 40 i4 = i3 + 1, iphct
                  idv(4) = i4
                  if (icp.eq.4) then
                     call testas (ier)
                     if (ier.eq.0) goto 99
                     goto 40
                   end if 
                  do 50 i5 = 14 + 1, iphct
                     idv(5) = i5
                     if (icp.eq.5) then
                        call testas (ier)
                        if (ier.eq.0) goto 99
                        goto 50
                     end if 
                     do 60 i6 = i5 + 1, iphct 
                        idv(6) = i6
                        if (icp.eq.6) then
                           call testas (ier)
                           if (ier.eq.0) goto 99
                           goto 60
                        end if 
                        do 70 i7 = i6 + 1, iphct
                           idv(7) = i7
                           call testas (ier)
                           if (ier.eq.0) goto 99

70    continue
60    continue
50    continue
40    continue
30    continue
20    continue
10    continue

999   write (*,*) ' yeech '
      goto 9999
99    write (*,*) ' no problem'
      i7fct = 1
      jasmbl(1) = 1

      do 90 i = 1, icp
90       id7f(i,1) = idv(i)
      write (*,*) ' phase indices:', (idv(i), i = 1, icp)

9999  end 

      subroutine testas (ier)
c--------------------------------------------------------------------
c combinatorial enumeration
 
c   referenced by: bossop
c   references to: abload,assd7,dgph7,sidc,slopc
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7,j9=1200)
 
      character*8 names

      integer jdv(k6)
 
      common/ cst6  /icomp,istct,iphct,icp/ cst12 /cp(k6,k1)
      common/ cst23 /a(10,10),u(11),ipvt(10),kdv(10),iophi,idphi,
     *               iiphi,iflg1/ cst8  /names(k1)/ cst2  /g(k1)
      common/ cst52 /hcp,l,idv(7),i0,ic(7)/ cst302 /jasmbl(j9)

      save d1ct, d2ct
      data d1ct, d2ct/0.d0,0.d0/
c-----------------------------------------------------------------------

      ier = 0
      d1ct = d1ct + 1
      if (d1ct.gt.1.d6) then
         d2ct = d2ct + 1 
         cfact = 1.d0
         do 1 i = iphct-icp,iphct
1           cfact = cfact*i
         do 2 i = 1, icp
2           cfact = cfact/i
         write (*,1000) d2ct*d1ct,cfact,1.d2*d2ct*d1ct/cfact
         d1ct = 0
      end if 
     
      do 10 i = 1, icp
10       jdv(i) = idv(i)

      call bulkno (jdv,jasmbl(1),ier)

      if (ier.ne.0) goto 999

      call abload (*991) 
c                             start test loop:
      do 20 j = 1, iphct 
         dg = g(j)
         do 30 i = 1, icp
            if (idv(i).eq.j) goto 20
            dg = dg - cp(i,j)*u(i)
30       continue
         if (dg.lt.-1.d-8) then
c                             plane is metastable
            goto 992
         end if
20    continue
c                             all ok
      call bulkck (jdv,jasmbl(i),ier)

      goto 999

991   ier = 1
      goto 999
992   ier = 2
      goto 999
993   ier = 3

1000  format (' Tested ',g16.7,' assemblages out of ',g16.7
     *        ' possible (',f8.4,' %)')
      
999   end 

      subroutine bulkno (jdv,jbfg,ier)
c----------------------------------------------------------------------
c bulkck determines phases identified by jdv bound
c any of the specified bulk compositions.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,nb=20,h5=7,k4=16,k0=25,k5=12,l2=5)
 
      integer jdv(7)

      common/ cst6  /icomp,istct,iphct,icp/ cst20 /comps(k1,k5)
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst301 /a(k5,k5),b(k5),ipvt(k5)
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
c----------------------------------------------------------------------
      jbfg = 0
c                                load the transpose of the
c                                concentration matrix of the pseudo-
c                                invariant assemblage.
      do 10 i = 1, jcth
         if (i.le.icp) then
            k = jdv(i)
         else
            k = idss(i-icp)
         end if
         do 20 j = 1, jcth
20          a(j,i) = comps(k,j)
10    continue       
c                                factor the matrix
      call factr1 (jcth,ier)
      if (ier.eq.1) goto 99 
      ier = 1
c                                test for bounded compositions
      do 50 j = 1, iblk

         do 60 i = 1, jcth
c                                load composition vector, the alpha
c                                vector is returned in the same array
60          b(i) = cblk(j,i)
c                                solve for the alpha vector
         call subst1 (jcth)
c                                a negative coefficient indicates
c                                the composition is not bounded
         do 70 i = 1, icp
70          if (b(i).lt.-1.d-5) goto 50
c                                the composition is bounded 
         jbfg = 1
         ier = 0

50    continue                
99    end
