 
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 vertex
c----------------------------------------------------------------------
c                       ************************
c                       *                      *
c                       *  vertex.nov.17.1991  *
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 dgdegn (zpot,zd)
c---------------------------------------------------------------------
c dgdegn - this subroutine checks the stability of a phase with
c respect to an n component assemblage, in which n-1 phases lie
c in a stable subcomposition, jj is the index of the free potential
c and the number of components is hcp, if the phase zd is
c stable with respect to the hcp-phase assemblage it replaces the
c phase id(hcp), and the potential jj is recalculated.
 
c referenced by: simpl1 - simpl7
c no external references.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(10,10),u(11),ipvt(10),idv(10),iophi,idphi,
     *               iiphi,iflg1
     *      / cst2  /g(k1)/ cst12 /cp(k6,k1)/ cst52 /hcp,h,id(8),ic(7)
c-----------------------------------------------------------------------
c                                        check stability of the
c                                        phase zd:
 
      dg = g(zd)
 
      do 1 i = 1, hcp
1        dg = dg - cp(ic(i),zd) * u(i)
c                                       if dg > 0, zd is metastable:
      if (dg.gt.-1.d-8) return
c                                       calculate the new potential:
      upot = g(zd)
 
      do 2 i = 1, hcp
         if (i.ne.zpot) then
            upot = upot - cp(ic(i),zd) * u(i)
         end if
2     continue
 
      upot = upot / cp(ic(zpot),zd)
c                                       remove this after debugging
      if (upot.gt.u(zpot)) call error (917,u(1),i,'DGDEGN')
 
      u(zpot) = upot
c                                       make zd the new test vertice:
      id(hcp) = zd
 
      end
 
      subroutine allone
c----------------------------------------------------------------------
c allone loads the indexes of the phases in each phase assemblage into
c the array idcf, and the assemblage count into the variable i7fct.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k6=7,m1=21,m2=35,m3=35,m4=21,m5=7,j9=1200,h4=20)
      parameter (j1=50,j2=400,j3=900,j4=1200,j5=2400,j6=200)
      parameter (k1=1200,k5=12,nb=20)

      integer jdv(k6)
 
      common/ cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst6  /icomp,istct,iphct,icp
     *      / cst96 /id7f(k6,j9),i7fct/ cst302 /jasmbl(j9)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth
c----------------------------------------------------------------------
      goto (100,200,300,400,500,600),icp
c                                 seven component system=default
      goto 700
c                                 one component system:
100   i7fct = 1
      id7f(1,1) = id1s(1)
      goto 700
c                                 two component system:
200   i7fct = i2vct(1)-1
      do 210 i = 1, i7fct
      id7f(1,i) = id2v(1,i)
210   id7f(2,i) = id2v(1,i+1)
      goto 700
c                                 three component systems:
300   i7fct = i3fct(1)
      do 310 i = 1, i7fct
      do 320 j = 1, 3
320   id7f(j,i) = id3f(1,i,j)
310   continue
      goto 700
c                                 four component systems:
400   i7fct = i4fct(1)
      do 410 i = 1, i7fct
      do 420 j = 1, 4
420   id7f(j,i) = id4f(1,i,j)
410   continue
      goto 700
c                                 five component systems:
500   i7fct = i5fct(1)
      do 510 i = 1, i7fct
         do 520 j = 1, 5
520         id7f(j,i) = id5f(j,i,1)
510   continue
      goto 700
c                                 six component systems:
600   i7fct = i6fct(1)
      do 610 i = 1, i7fct
         do 620 j = 1, 6
620         id7f(j,i) = id6f(j,i,1)
610   continue
c                                do bulk comp checks
700   if (iblk.ne.0) then 
         do 710 i = 1, i7fct
            do 720 j = 1, icp
720            jdv(j) = id7f(j,i)
710            call bulkck (jdv,jasmbl(i),ier)
      end if
 
      end
 
      subroutine assd3 (istart,hgo)
c----------------------------------------------------------------------
c assd3 counts (i3fcth) and assigns (id3f) ternary facets. the
c vertices are identified by the array id(3).
 
c   referenced by: simpl4
c   references to: none
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (j1=50,j2=400,j6=200,k1=1200,k6=7,l5=400,l6=400)
      parameter (m1=21,m2=35,j8=4800)
 
      integer ilk(2)
 
      common/ cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst45 /itbh(3),itt,i3vcth,i3fcth
     *      / cst13 /id1(k6,l5),id2(m1,l6)/ debug /jtest,jpot
     *      / cst52 /hcp,h,id(8),itc1,itc2,itc3,itc4,ic(3)
     *      / cst19 /hi,iqth(4),ivchk(k1)/ cst6  /icomp,istct,iphct,icp
     *      / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst18 /icct(k1)/ cst46 /ipsf,idpsf(j8,k6)
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
 
      save ilk
 
      data ilk/ 2, 1/
c-----------------------------------------------------------------------
      if (i3fcth.eq.0) goto 40
c                             test for equivalence with previously
c                             identified facets:
      i1 = id(1)
      i2 = id(2)
      i3 = id(3)
      do 10 i = 1, i3fcth
         if (
     *     ((id3f(h,i,1).ne.i1).and.(id3f(h,i,2).ne.i1).and.
     *      (id3f(h,i,3).ne.i1)).or.
     *     ((id3f(h,i,1).ne.i2).and.(id3f(h,i,2).ne.i2).and.
     *      (id3f(h,i,3).ne.i2)).or.
     *     ((id3f(h,i,1).ne.i3).and.(id3f(h,i,2).ne.i3).and.
     *      (id3f(h,i,3).ne.i3))) goto 10
c                           match found:
         return
10    continue
c                           unique:
40    i3fcth = i3fcth+1
      goto (90),jpot
c                           print potentials:
      if (hcp.eq.icp) call prtpot 

90    if (i3fcth.gt.j2) call error (189,1.d0,j2,'ASSD3')
 
      do 50 i = 1, 3
         if (ivchk(id(i)).eq.1) goto 50
         ivchk(id(i)) = 1
         if (icct(id(i)).lt.3) goto 50
         i3vcth = i3vcth+1
         if (i3vcth.gt.j6) call error (193,1.d0,j6,'ASSD3')
c                           hey !! kill this after debugging !!
         if (jtest.eq.0) call testit
 
         id3v(h,i3vcth) = id(i)
50    id3f(h,i3fcth,i) = id(i)
c                             test tielines for equivalence:
      do 80 i = 1,2
         i1 = id(i2c(1,i+1))
         i2 = id(i2c(2,i+1))
         if (ipsf.eq.0) goto 25
         do 20 j = istart,ipsf
            if (idpsf(j,3).eq.0) goto 20
            if (((idpsf(j,1).ne.i1) .and.
     *           (idpsf(j,2).ne.i1)).or.
     *          ((idpsf(j,1).ne.i2) .and.
     *           (idpsf(j,2).ne.i2))) goto 20
            idpsf(j,3) = 0
            goto 80
20       continue
c                             test with binary joins
25       goto (60),hgo
         do 30 j = 1, 3
            l = itbh(j)
            ib = i2vct(l)-1
            do 70 k = 1, ib
70             if (((id2v(l,k)  .eq.i1) .or.
     *              (id2v(l,k+1).eq.i1)).and.
     *             ((id2v(l,k)  .eq.i2) .or.
     *              (id2v(l,k+1).eq.i2))) goto 80
30       continue
c                            assign new tieline:
60       ipsf = ipsf+1
         if (ipsf.gt.j8) call error (192,1.d0,j8,'ASSD3')
         idpsf(ipsf,1) = i1
         idpsf(ipsf,2) = i2
         idpsf(ipsf,3) = id(ilk(i))
80    continue
 
      end
 
      subroutine assd4 (istart,hgo)
c-----------------------------------------------------------------------
c assd4 counts (i4fcth) and assigns (id4f) quaternary facets. the
c vertices are identified by the array id(4).
 
c   referenced by: simpl4
c   references to: none
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (j2=400,j3=900,j6=200,j7=500)
      parameter (k1=1200,k6=7,m2=35,m3=35,j8=4800)
 
      integer ilk(4)
 
      common/ cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst56 /id4v(m3,j7),i4vct(m3),i4vcth/ cst18 /icct(k1)
     *      / cst52 /hcp,h,id(8),iqc1,iqc2,iqc3,iqc4,iqc5,ic(2)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst19 /hi,iqth(4),ivchk(k1)/ cst6  /icomp,istct,iphct,icp
     *      / cst46 /ipsf,idpsf(j8,k6)/ debug /jtest,jpot
 
      save ilk
 
      data ilk /0, 3, 2, 1/
c-----------------------------------------------------------------------
      if (i4fcth.lt.1) goto 40
c                             test for equivalence with earlier
c                             facets:
      do 10 i = 1, i4fcth
 
      do 20 j = 1, 4
      idi = id4f(h,i,j)
 
      do 30 k = 1, 4
30    if (idi.eq.id(k)) goto 20
 
      goto 10
20    continue

      igo = 1
      return
 
10    continue
c                             unique facet:
40    igo = 0
      i4fcth = i4fcth+1
      if (i4fcth.gt.j3) call error (190,1.d0,j3,'ASSD4')

      goto (120),jpot
c                             print potentials:
      if (hcp.eq.icp) call prtpot 
c                             kill after debugging!!
120   if (jtest.eq.0) call testit
 
      do 50 i = 1, 4
         idi = id(i)
c                             set flag ivchk for new vertices:
         if (icct(idi).lt.4.or.ivchk(idi).eq.1) goto 50
         i4vcth = i4vcth+1
         if (i4vcth.gt.j7) call error (194,1.d0,j7,'ASSD4')
         ivchk(idi) = 1
         id4v(h,i4vcth) = idi
50    id4f(h,i4fcth,i) = idi
c                             each new facet may generate 3 new pst's
c                             defined by the vertices (124),(134),
c                             and (234).
      do 60 i = 2,4
c                             set vertices identities:
      i1 = id(i3c(1,i))
      i2 = id(i3c(2,i))
      i3 = id(i3c(3,i))
c                             begin by testing each potential pst
c                             for a match with one already defined.
      if (ipsf.eq.0) goto 70
      do 80 j = istart,ipsf
      if ((idpsf(j,1).eq.0).or.
     *   ((idpsf(j,1).ne.i1).and.(idpsf(j,2).ne.i1).and.
     *    (idpsf(j,3).ne.i1)).or.
     *   ((idpsf(j,1).ne.i2).and.(idpsf(j,2).ne.i2).and.
     *    (idpsf(j,3).ne.i2)).or.
     *   ((idpsf(j,1).ne.i3).and.(idpsf(j,2).ne.i3).and.
     *    (idpsf(j,3).ne.i3))) goto 80
c                             matched pst.
      idpsf(j,1) = 0
      goto 60
80    continue
c                             test with ternaries if simpl4 is
c                             executing the first search seg-
c                             ment (hgo = 0).
70    goto (90),hgo
      do 100 j = 1,4
      iqthj = iqth(j)
      if (iqthj.eq.hi) goto 100
      hh = i3fct(iqthj)
      do 110 k = 1,hh
      if (
     * (id3f(iqthj,k,1).ne.i1.and.id3f(iqthj,k,2).ne.i1.and.
     *  id3f(iqthj,k,3).ne.i1).or.
     * (id3f(iqthj,k,1).ne.i2.and.id3f(iqthj,k,2).ne.i2.and.
     *  id3f(iqthj,k,3).ne.i2).or.
     * (id3f(iqthj,k,1).ne.i3.and.id3f(iqthj,k,2).ne.i3.and.
     *  id3f(iqthj,k,3).ne.i3)) goto 110
      goto 60
110   continue
100   continue
c                             unique pst:
90    ipsf = ipsf+1
      if (ipsf.gt.j8) call error (192,1.d0,j8,'ASSD4')
      idpsf(ipsf,1) = i1
      idpsf(ipsf,2) = i2
      idpsf(ipsf,3) = i3
      idpsf(ipsf,4) = id(ilk(i))
 
60    continue
 
      end
 
      subroutine assd5 (istart,hgo,jd4)
c-----------------------------------------------------------------------
c assd5 counts (i5fctl) and assigns (idcf) quinary facets. the
c vertices are identified by the array id(5).
 
c   referenced by: simpl5
c   references to: none
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7,l5=400,l6=400,l7=400,l8=200,l9=100,j3=900,
     *           m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800,
     *           j4=1200,h3=250)
 
      integer ilk(5)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ cst6  /icomp,istct,iphct,icp
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst98 /id5v(h3,m3),i5vct(m3),i5vctl
     *      / cst52 /hcp,jd5,id(8),iqc1,iqc2,iqc3,iqc4,iqc5,ic(2)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst19 /hi,iqth(4),ivchk(k1)/ debug /jtest,jpot
     *      / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct/ cst18 /icct(k1)
     *      / cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
     *      / cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
 
      save ilk
 
      data ilk /0, 4, 3, 2, 1/
c-----------------------------------------------------------------------
 
      if (i5fctl.lt.1) goto 40
c                             test for equivalence with earlier
c                             facets:
      do 10 i = 1, i5fctl
         do 20 j = 1, 5
            idi = id5f(j,i,jd5)
            do 30 k = 1, 5
30            if (idi.eq.id(k)) goto 20
              goto 10
20       continue
c                             matches earlier facet, return.
         return
10    continue
c                             unique facet:
40    i5fctl = i5fctl+1

      goto (170),jpot
c                             print potentials:
      if (hcp.eq.icp) call prtpot 
c                             kill after debugging!!!
170   if (jtest.eq.0) call testit
c
      if (i5fctl.gt.j4) call error (191,1.d0,j4,'ASSD5')
 
      do 50 i = 1, 5
         idi = id(i)
c                             set flag ivchk for new vertices:
         if (icct(idi).lt.5.or.ivchk(idi).eq.1) goto 50
         i5vctl = i5vctl + 1
         if (i5vctl.gt.h3) call error (201,1.d0,h3,'ASSD5')
         ivchk(idi) = 1
         id5v(i5vctl,jd5) = idi
 
50       id5f(i,i5fctl,jd5) = idi
c                             each new facet may generate 4 new pst's
c                             defined by the vertices (1235), (1245),
c                             (1345), and (2345):
      do 60 i = 2, 5
c                             set vertice identities:
c                             begin by testing each potential pst
c                             for a match with one already defined.
         if (ipsf.eq.0) goto 70
 
         do 80 j = istart, ipsf
            do 90 k = 1, 4
               do 100 l = 1, 4
100               if (id(i4c(l,i)).eq.idpsf(j,k)) goto 90
               goto 80
90          continue
c                             matched an earlier psq:
            idpsf(j,1) = 0
            goto 60
80       continue
c                             test with quaternaries if simpl5 is
c                             executing the first search seg-
c                             ment (hgo = 0).
70    goto (110),hgo
 
      do 120 j = 1, 5
         kd4 = i45(j,jd5)
         if (jd4.eq.kd4) goto 120
         ib = i4fct(kd4)
         do 130 k = 1, ib
            do 140 l = 1, 4
               do 150 m = 1, 4
150               if (id(i4c(l,i)).eq.id4f(kd4,k,m)) goto 140
               goto 130
140         continue
c                             match with a quaternary facet
            goto 60
130      continue
120   continue
c                             unique pst:
110   ipsf = ipsf+1
 
      if (ipsf.gt.j8) call error (192,1.d0,j8,'ASSD5')
 
      do 160 j = 1,4
160      idpsf(ipsf,j) = id(i4c(j,i))
c                              why save this one?
      idpsf(ipsf,5) = id(ilk(i))
 
60    continue
 
      end
 
      subroutine assd6 (istart,hgo,jd5)
c-----------------------------------------------------------------------
c assd6 counts (i6fctl) and assigns (idcf) quinary facets. the
c vertices are identified by the array id(6).
 
c   referenced by: simpl6
c   references to: none
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (j3=900,j4=1200,j5=2400)
      parameter (k1=1200,k6=7,h4=20)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
 
      integer ilk(6)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ debug /jtest,jpot
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl/ cst18 /icct(k1)
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst19 /hi,iqth(4),ivchk(k1)
     *      / cst52 /hcp,jd6,id(8),iqc1,iqc2,iqc3,iqc4,iqc5,iqc6,ic7
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct/ cst6  /icomp,istct,iphct,icp
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
 
      save ilk
 
      data ilk/ 0, 5, 4, 3, 2, 1/
c-----------------------------------------------------------------------
 
      if (i6fctl.lt.1) goto 40
c                             test for equivalence with earlier
c                             facets:
      do 10 i = 1, i6fctl
         do 20 j = 1,6
            idi = id6f(j,i,jd6)
            do 30 k = 1,6
30            if (idi.eq.id(k)) goto 20
              goto 10
20       continue
c                             matches earlier facet, return.
         return
10    continue
c                             unique facet:
40    i6fctl = i6fctl+1
      goto (170),jpot
c                             print potentials:
      if (hcp.eq.icp) call prtpot
c                             kill after debugging!!!
170   if (jtest.eq.0) call testit

      if (i6fctl.gt.j5) call error (202,1.d0,j5,'ASSD6')
 
      do 50 i = 1,6
         idi = id(i)
c                             set flag ivchk for new vertices:
         if (icct(idi).lt.6.or.ivchk(idi).eq.1) goto 50
         i6vctl  =  i6vctl + 1
         if (i6vctl.gt.h4) call error (203,1.d0,h4,'ASSD6')
         ivchk(idi) = 1
         id6v(i6vctl,jd6) = idi
 
50       id6f(i,i6fctl,jd6) = idi
c                             each new facet may generate 4 new pst's
c                             defined by the vertices (12346), (12356),
c                             (12456), (13456), and (23456):
      do 60 i = 2, 6
c                             set vertices identities:
c                             begin by testing each potential pst
c                             for a match with one already defined.
         if (ipsf.eq.0) goto 70
 
         do 80 j = istart,ipsf
            do 90 k =  1, 5
               do 100 l =  1, 5
100               if (id(i5c(l,i)).eq.idpsf(j,k)) goto 90
               goto 80
90          continue
c                             matched an earlier psq:
            idpsf(j,1) = 0
            goto 60
80       continue
c                             test with quinaries if simpl6 is
c                             executing the first search seg-
c                             ment (hgo = 0).
70    goto (110),hgo
 
      do 120 j =  1, 6
         kd5 = i56(j,jd6)
         if (jd5.eq.kd5) goto 120
         ib = i5fct(kd5)
         do 130 k = 1, ib
            do 140 l = 1,5
               do 150 m = 1,5
150               if (id(i5c(l,i)).eq.id5f(m,k,kd5)) goto 140
               goto 130
140         continue
c                             match with a quaternary facet
            goto 60
130      continue
120   continue
c                             unique pst:
110   ipsf = ipsf+1
 
      if (ipsf.gt.j8) call error (192,1.d0,j8,'ASSD6')
 
      do 160 j = 1,5
160      idpsf(ipsf,j) = id(i5c(j,i))
c                              why save this one?
      idpsf(ipsf,6) = id(ilk(i))

60    continue
 
      end
 
      subroutine assd7 (istart,hgo,jd6)
c-----------------------------------------------------------------------
c assd7 counts (i7fct) and assigns (id7f) septary facets. the
c vertices are identified by the array id(7).
c
c   referenced by: simpl7
c   references to: none
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (j3=900,j4=1200,j5=2400)
      parameter (k1=1200,k6=7,h4=20,j9=1200)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
 
      integer ilk(7)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ debug /jtest,jpot
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst19 /hi,iqth(4),ivchk(k1)
     *      / cst96 /id7f(k6,j9),i7fct/ cst18 /icct(k1)
     *      / cst52 /hcp,jd7,id(8),iqc1,iqc2,iqc3,iqc4,iqc5,iqc6,ic7
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct/ cst6  /icomp,istct,iphct,icp
     *      / cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
     *      / cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
 
      save ilk
 
      data ilk /0, 6, 5, 4, 3, 2, 1/
c-----------------------------------------------------------------------
 
      if (i7fct.lt.1) goto 40
c                             test for equivalence with earlier
c                             facets:
      do 10 i = 1, i7fct
         do 20 j = 1, 7
            idi = id7f(j,i)
            do 30 k = 1, 7
30            if (idi.eq.id(k)) goto 20
              goto 10
20       continue
c                             matches earlier facet, return.
         return
10    continue
c                             unique facet:
40    i7fct = i7fct +1

      goto (170), jpot
c                             print potentials:
      if (hcp.eq.icp) call prtpot

170   if (i7fct.gt.j9) call error (204,1.d0,j9,'ASSD7')
 
      do 50 i = 1, 7
         idi = id(i)
50       id7f(i,i7fct) = idi
 
      if (jtest.eq.0) call testit
c                             each new facet may generate 6 new pst's
c                             defined by the vertices (12346), (12356),
c                             (12456), (13456), and (23456):
      do 60 i = 2, 7
c                             set vertices identities:
c                             begin by testing each potential pst
c                             for a match with one already defined.
         if (ipsf.eq.0) goto 70
 
         do 80 j = istart, ipsf
            do 90 k = 1, 6
               do 100 l = 1, 6
100               if (id(i6c(l,i)).eq.idpsf(j,k)) goto 90
               goto 80
90          continue
c                             matched an earlier psq:
            idpsf(j,1) = 0
            goto 60
80       continue
c                             test with quinaries if simpl7 is
c                             executing first search segment (hgo = 0)
70    goto (110),hgo
 
      do 120 j = 1, 7
         if (jd6.eq.j) goto 120
         ib = i6fct(j)
         do 130 k = 1, ib
            do 140 l = 1, 6
               do 150 m = 1, 6
150               if (id(i6c(l,i)).eq.id6f(m,k,j)) goto 140
               goto 130
140         continue
c                             match with a hexary facet
            goto 60
130      continue
120   continue
c                             unique pst:
110   ipsf = ipsf+1
 
      if (ipsf.gt.j8) call error (192,1.d0,j8,'ASSD7')
 
      do 160 j = 1, 6
160      idpsf(ipsf,j) = id(i6c(j,i))
c                              why save this one?
      idpsf(ipsf,7) = id(ilk(i))

60    continue
 
      end
 
      subroutine bossop
c----------------------------------------------------------------------
c bossop calls subroutines to determine the stable phase
c assemblages of progressively more complex subsystems.
c after the assemblages have been determined by the simplc routines
c the routine miscib is called to test for immiscibility.
c   referenced by: main
c   references to: simpl1,simpl2,simpl3,simpl4,simpl5,miscib
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,j1=50,j2=400,j3=900,
     *           j4=1200,j5=2400,j6=200,j9=1200,h4=20,k1=1200,k6=7)
 
      common/ cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst6  /icomp,istct,iphct,icp
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst60 /iasmbl(j9),ipoint,imyn
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst19 /hi,iqth(4),ivchk(k1)/ cst96 /id7f(k6,j9),i7fct
c-----------------------------------------------------------------------
c                                 why this?
      do 1 i = 1, iphct
1        ivchk(i) = 0
 
      inull = 0
 
      imyn = 1
c                                 one component systems:
      call simpl1
      if (icp.eq.1) return
c                                 two component systems:
      call simpl2
c                                 test for immiscibility and variance
      if (icp.gt.2) goto 300
      ict = i2vct(1)
 
      do 210 i = 2,ict
210      call miscib (icp,i,id2v(1,i-1),id2v(1,i),inull,inull,
     *                inull,inull,inull)
 
      return
c                                 three component systems:
300   call simpl3 (iflg4)
c                                 test for immiscibility and variance
      if (icp.gt.3) goto 400
      ict = i3fct(1)
 
      do 310 i = 1, ict
310      call miscib (icp,i,id3f(1,i,1),id3f(1,i,2),id3f(1,i,3),
     *                inull,inull,inull,inull)
 
      return
c                                 four component systems:
400   call simpl4 (iflg4)
c                                 test for immiscibilty
      if (icp.gt.4) goto 500
 
      do 410 i = 1, i4fcth
410      call miscib (icp,i,
     *        id4f(1,i,1),id4f(1,i,2),id4f(1,i,3),id4f(1,i,4),
     *        inull,inull,inull)
 
      return
c                                 five components:
500   call simpl5
c                                 test for immiscibility
      if (icp.gt.5) goto 600
 
      do 510 i = 1, i5fctl
 
510       call miscib (icp,i,
     *         id5f(1,i,1),id5f(2,i,1),id5f(3,i,1),id5f(4,i,1),
     *         id5f(5,i,1),inull,inull)
 
      return
c                                 six components:
600   call simpl6
c
      if (icp.gt.6) goto 700
 
      do 610 i = 1, i6fctl
 
610       call miscib (icp,i,
     *         id6f(1,i,1),id6f(2,i,1),id6f(3,i,1),id6f(4,i,1),
     *         id6f(5,i,1),id6f(6,i,1),inull)
 
      return
 
700   call simpl7
 
      do 710 i = 1, i7fct
 
710       call miscib (icp,i,
     *         id7f(1,i),id7f(2,i),id7f(3,i),id7f(4,i),
     *         id7f(5,i),id7f(6,i),id7f(7,i))
 
      end
c
      subroutine bound3
c----------------------------------------------------------------------
c bound3 calculates the parameters for test
c i.e. the slopes of the tielines connecting ternary vertices
c and the "sign" of included points.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
c
      common/ cst3 /x(k1,k6)
     *      / cst17 /s1,s2,s3,b1,b2,b3,dy(3),xb(3),dx(3)
     *      / cst45 /itbh(3),itt,i3vcth,i3fcth
     *      / cst52 /hcp,h,i1,i2,i3,kd,ld,md(3),
     *               itc1,itc2,itc3,itc4,ic(3)
c-----------------------------------------------------------------------
      dxx = (x(i2,itc3)-x(i1,itc3))
      if (dxx.ne.0.d0) goto 5
      dx(1) = x(i1,itc3) - x(i3,itc3)
      xb(1) = x(i1,itc3)
      s1 = 1.d30
      goto 10
5     s1 = (x(i2,itc2) - x(i1,itc2))/dxx
      b1 = x(i2,itc2) - s1*x(i2,itc3)
      dy(1) = s1*x(i3,itc3) + b1 - x(i3,itc2)
10    dxx = (x(i3,itc3)-x(i1,itc3))
      if (dxx.ne.0.d0) goto 15
      dx(2) = x(i1,itc3)-x(i2,itc3)
      xb(2) = x(i1,itc3)
      s2 = 1.d30
      goto 20
15    s2 = (x(i3,itc2)-x(i1,itc2))/dxx
      b2 = x(i3,itc2)-s2*x(i3,itc3)
      dy(2) = s2*x(i2,itc3) + b2 - x(i2,itc2)
20    dxx = (x(i3,itc3) - x(i2,itc3))
      if (dxx.ne.0.d0) goto 25
      s3 = 1.d30
      dx(3) = x(i2,itc3) - x(i1,itc3)
      xb(3) = x(i2,itc3)
      goto 30
25    s3 = (x(i3,itc2) - x(i2,itc2))/dxx
      b3 = x(i3,itc2) - s3*x(i3,itc3)
      dy(3) = s3*x(i1,itc3) + b3 - x(i1,itc2)
 
30    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
 
      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
c                                  generate permanent graphics file
c                                  if requested
      goto (75),io4
      call outgrf
c                                  printer output:
75    goto (90),io3
      call outchm
c                                  do bulk composition checks:
      if (iblk.ne.0) call allone

      goto 90
 
99    stop

      end
 
c-----------------------------------------------------------------------
c                      dgphc
c the dgphc functions determine the g difference between a phase
c identified by 'ld', and the plane defined by 'u'.
 
c   referenced by: simpl3-simpl7
c-----------------------------------------------------------------------
      function dgph3 (ld)
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(100),u(11),i(24)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst52 /l(10),itc1,itc2,itc3,itc4,ic(3)
 
      dgph3 = g(ld)-cp(itc1,ld)*u(1)
     *           -cp(itc2,ld)*u(2)
     *           -cp(itc3,ld)*u(3)
 
      end
 
      function dgph4 (ld)
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(100),u(11),i(24)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst52 /l(10),itc1,itc2,itc3,itc4,ic(3)
 
      dgph4 = g(ld)-cp(itc1,ld)*u(1)
     *           -cp(itc2,ld)*u(2)
     *           -cp(itc3,ld)*u(3)
     *           -cp(itc4,ld)*u(4)
 
      end
 
      function dgph5 (ld)
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(100),u(11),i(24)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst52 /l(10),itc1,itc2,itc3,itc4,itc5,ic(2)
 
      dgph5=g(ld)-cp(itc1,ld)*u(1)
     *           -cp(itc2,ld)*u(2)
     *           -cp(itc3,ld)*u(3)
     *           -cp(itc4,ld)*u(4)
     *           -cp(itc5,ld)*u(5)
      end
 
      function dgph6 (ld)
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(100),u(11),i(24)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst52 /l(10),itc1,itc2,itc3,itc4,itc5,itc6,itc7

      dgph6 = g(ld)-cp(itc1,ld)*u(1)
     *           -cp(itc2,ld)*u(2)
     *           -cp(itc3,ld)*u(3)
     *           -cp(itc4,ld)*u(4)
     *           -cp(itc5,ld)*u(5)
     *           -cp(itc6,ld)*u(6)
      end
 
      function dgph7 (ld)
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst23 /a(100),u(11),i(24)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst52 /l(10),itc1,itc2,itc3,itc4,itc5,itc6,itc7
 
      dgph7 = g(ld)-cp(itc1,ld)*u(1)
     *           -cp(itc2,ld)*u(2)
     *           -cp(itc3,ld)*u(3)
     *           -cp(itc4,ld)*u(4)
     *           -cp(itc5,ld)*u(5)
     *           -cp(itc6,ld)*u(6)
     *           -cp(itc7,ld)*u(7)
      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*5
 
      common/ cst6  /icomp,istct,iphct,icp/ cst79 /isoct
     *      / csta7 /fname(h9) / cst24 /ipot,jv(l2),iv(l2)
     *      / cst81 /ipchk(k3),ip2(k3),icp2,ifake/ cst112 /buf(5)
     *      / 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
     *      / cst209 /zname

      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
      if (ibuf.eq.5) read (n1,*,err=998) buf
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 computational option file name ',
     *       '(i.e. the file created',/,' with BUILD)',
     *       ', 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

      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
470      if (icout(j).eq.0.and.comp(j).ne.0.) goto 50
      do 473 j = 1, icmpn
         if (comp(j).lt.0.) then
            call warn (13,r,j,name)
            goto 50
         end if
473   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                               classify into a subcomposition:
      call sort (ier)
c                               next phase:
      goto 220
 
999   do 99 i = 1, icp
99       if (i1ct(i).lt.1) call error (15,r,i,cname(i))
 
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'
     *       ,/)
      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
 
      function itest3 (ld)
c----------------------------------------------------------------------
c itest3 determines whether a point lies within the portion
c of a plane bounded by the lines connecting 3 points.
c the c points.  the identity of the
c point to be tested is given by ld.
 
c itest3 returns 0 if the point is within the plane, 1 if it is
c outside, and 2 if it lies on an edge of the bounded region.
 
c   referenced by: simpl3
c   references to: subst,factor
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
 
      common/ cst3 /x(k1,k6)/ cst17 /s(3),b(3),dy(3),xb(3),dx(3)
     *      / cst45 /itbh(3),itt,i3vcth,i3fcth
     *      / cst52 /hcp,h,hd,id,jd,kd,ll,md(3),
     *               itc1,itc2,itc3,itc4,ic(3)
c-----------------------------------------------------------------------
      lt = 1
      do 10 i = 1, 3
         if (s(i).ne.1.d30) goto 30
         t = dx(i)*(xb(i)-x(ld,itc3))
         goto 20
30       t = dy(i)*((s(i)*x(ld,itc3)+b(i))-x(ld,itc2))
20       if (dabs(t).lt.1.d-08) lt = 2
         if (t.ge.0.d0) goto 10
         itest3 = 1
         return
10    continue
      itest3 = 0
      goto (99),lt
      itest3 = 2
 
99    end
 
      subroutine miscib (hcp,ict,hd,id,jd,kd,ld,md,zd)
c----------------------------------------------------------------------
c miscib determines the variance of the assemblage of hcp phases
c identified by the array id and the counter ict. miscib also
c checks for immiscibility in solution phases.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k2=1200,k3=1200,k6=7,j9=1200)
 
      dimension iok(2),jdv(7)
 
      common/ cst19 /hi,iqth(4),ivchk(k1)
     *      / cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi,
     *               iiphi,iflg1
     *      / cst12 /cp(k6,k1)/ cst60 /iasmbl(j9),ipoint,imyn
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst6  /icomp,istct,iphct,icp
c----------------------------------------------------------------------
 
      goto (1,2,3,4,5,6),hcp
 
      jdv(7) = zd
6     jdv(6) = md
5     jdv(5) = ld
4     jdv(4) = kd
3     jdv(3) = jd
2     jdv(2) = id
1     jdv(1) = hd
c                                initialize flags and counters
      iok(1) = 0
      ihom = 1
      iok(2) = 0
      im = 0
      if (ict.gt.j9) call error (22,b(1),j9,'MISCIB')
      iasmbl(ict) = 1
      ii = 1
c                                test if the phase idv(ii) is a
c                                a pseudo-compound not previously
c                                matched
50    ikt = ikp(jdv(ii))
      if (ikt.eq.0.or.ikt.eq.iok(1).or.ikt.eq.iok(2)) goto 20
c                                look for a match
      jj = ii+1
      do 30 j = jj, hcp
         if (ikp(jdv(j)).ne.ikt) goto 30
c                                test for possibility of a homogeneous
c                                phase
         if (jj.ne.2) goto 40
         ihom = ihom + 1
30    continue
      if (ihom.eq.1) goto 20
      goto 40
c                                no match or previously found
20    ii = ii+1
      if (ii.eq.hcp) goto 99
      goto 50
c                                match found, set flags. note that the
c                                maximum number of matches in a five
c                                component system is two, increase the
c                                dimension of iok for larger systems.
40    im = im + 1
      iok(im) = ikt
      iasmbl(ict) = 2
c                                load and factor the transpose of the
c                                concentration matrix of the pseudo-
c                                invariant assemblage.
      do 70 i = 1, hcp
         do 80 j = 1, hcp
80          a(j,i) = cp(j,jdv(i))
70    continue
 
      call factor (a,hcp,ipvt,ier)
      if (ier.eq.1) call error (106,b(1),i,'MISCIB')
 
c                                test for bounded metastable pseudo-cpds
c                                of the solution ikt
      do 60 j = istct, iphct

         if (ivchk(j).eq.1.or.ikp(j).ne.ikt) goto 60

         do 90 i = 1, hcp
c                                load composition vector, the alpha
c                                vector is returned in the same array
90          b(i) = cp(i,j)
c                                solve for the alpha vector
      call subst (a,ipvt,hcp,b)
c                                a negative coefficient indicates
c                                the phase is not bounded
      do 100 i = 1, hcp
100      if (b(i).lt.1.d-5) goto 60
c                                the phase is bounded, exsolution has
c                                occurred, test to determine if degen-
c                                erate solvus has been found
      if (ihom.eq.hcp) goto 110
      do 120 i = 1, hcp
120      if (ikp(jdv(i)).ne.ikt) goto 60
c
110   iasmbl(ict) = 4
      imyn = 0
      goto 99
 
60    continue
c                                 if no immiscibility and the assemblage
c                                 is a homogeneous phase set asmbl flag
      if (ihom.eq.hcp) iasmbl(ict) = 3
 
99    end

      subroutine outchm
c-------------------------------------------------------------------
c outchm writes new chemographies to unit n3 as they are generated
c by unary, binary, ternary, or check.  the flag ichem( = 1)
c determines if outchm is called, no output will be generated if the
c flag io3 = 1 (see introductory comments).
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
 
      common/ cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst5 /v(l2),tr,pr,r,ps/ cst13 /id1(k6,l5),id2(m1,l6)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst6 /icomp,istct,iphct,icp/ cst22 /id3(m2,l7)
     *      / 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
      common/ cst45 /itbh(3),itt,i3vcth,i3fcth/ cst8  /names(k1)
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst52 /hcp,h,hd,id,jd,kd,ld,md(3),
     *               itc1,itc2,itc3,itc4,ic(3)
     *      / cst60 /iasmbl(j9),ipoint,imyn/ cst79 /isoct
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst96 /id7f(k6,j9),i7fct
     *      / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2)
c-----------------------------------------------------------------------
      goto (999),io3
c                             output independent potential variable
c                             constraints:
      write (n3,1160)
      write (n3,1000)
      write (n3,1010) (vname(jv(z)),v(jv(z)),z = 1, ipot)
      write (n3,1015)
c                             next the stable assemblages:
      goto (100,200,300,400,500,600,700),hcp
c                             output unary chemographies.
100   write (n3,1020) names(id1s(1))
      goto 800
c                             output binary chemographies.
200   ib = i2vct(1)
      write (n3,2020) (names(id2v(1,zj)),zj = 1, ib)
      goto 800
c                             output ternary chemographies:
300   write (n3,3020) ((names(id3f(1,zi,zj)),zj = 1, 3),iasmbl(zi)
     *               ,zi = 1, i3fcth)
      goto 800
c                             output quaternary chemographies:
400   write (n3,4020) ((names(id4f(1,zi,zj)),zj = 1,4),iasmbl(zi),
     *                zi = 1, i4fcth)
      goto 800
c                             output quinary chemographies:
500   write (n3,5020) ((names(id5f(zj,zi,1)),zj = 1,5),iasmbl(zi),
     *                zi = 1, i5fctl)
      goto 800
c                             output hexary chemographies:
600   write (n3,6020) ((names(id6f(zj,zi,1)),zj = 1,6),iasmbl(zi),
     *                zi = 1, i6fctl)
      goto 800
c                             output septary chemographies:
700   write (n3,7020) ((names(id7f(zj,zi)),zj = 1, 7),iasmbl(zi),
     *                zi = 1, i7fct)
c                             output phases consistent with component
c                             saturation constraints:
800   goto (9000),isyn
      write (n3,6000)
      write (n3,2020) (names(idss(z)),z = 1, isat)
c
9000  goto (999),icp
      write (n3,*)
c                             write miscibility flag:
      if (isoct.eq.0) goto 999

      if (imyn.eq.1) then
         write (n3,6010) 
      else
         write (n3,6030) 
      end if
 
1000  format (1x,'the stable assemblages at:',/)
1010  format (26x,a8,' = ',g12.6)
1015  format (/,1x,'are:',/)
1020  format (26x,a8)
1160  format (/,'-------------------------------------------------'
     *          ,'---------------',/)
2020  format (6(2x,a8))
3020  format (2(1x,2(a8,'-'),a8,'(',i1,')',6x))
4020  format (2(1x,3(a8,'-'),a8,'(',i1,')',5x))
5020  format (2(1x,4(a8,'-'),a8,'(',i1,')',4x))
6000  format (/,1x,'these assemblages are compatible with the followi',
     *        'ng phases or species',/,' determined by component ',
     *        'saturation or buffering constraints:',/)
6010  format (' ** no immiscibility occurs in the stable solution ',
     *        'phases **',/)
6020  format (2(1x,5(a8,'-'),a8,'(',i1,')',4x))
6030  format (' ** immiscibility occurs in one or more of the ',
     *        'stable solution phases **',/)
7020  format (6(a8,'-'),a8,1x,'(',i1,')')
 
999   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 prtpot
c---------------------------------------------------------------------
c prtpot outputs the dependent potentials for phase assemblages
c determined by the simplc subroutines. it is called if the flag
c jpot, in common debug, is not equal 1.
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)

      parameter (k1=1200,k6=7) 

      character*8 names

      common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst23 /a(100),u(11),iar(24)/ cst8 /names(k1)
     *      / cst6  /icomp,istct,iphct,icp/ cst3 /x(k1,k6)
     *      / cst52 /hcp,h,id(8),ic(7)/ debug /jtest,jpot

      if (jpot.ne.0) goto 100

      goto (10,20,30,40,50,60,70),icp
10    goto 99
20    write (n3,1020) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99 
30    write (n3,1030) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99
40    write (n3,1040) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99 
50    write (n3,1050) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99 
60    write (n3,1060) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99 
70    write (n3,1070) (names(id(i)),i=1,icp),(u(i),i=1,icp)
      goto 99

100   do 80 i = 1, icp
            g = 0.0
            do 90 j = 1, icp
90             g = g + x(id(i),j) * u(j)
80          write (n3,1080) names(id(i)),g,(x(id(i),j),j=1,icp) 
          
1020  format (2(1x,a8),2(1x,g14.7))
1030  format (3(1x,a8),3(1x,g14.7))
1040  format (4(1x,a8),4(1x,g14.7))
1050  format (5(1x,a8),5(1x,g14.7))
1060  format (6(1x,a8),/,6(1x,g14.7))
1070  format (7(1x,a8),/,7(1x,g14.7))
1080  format (1x,a8,8(1x,g14.7))

99    end 
 
      subroutine simpl1
c----------------------------------------------------------------------
c simpl1 determines the stable phase in unary (sub)systems.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7,j1=50)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100)

      common/ cst48 /u(k6)/ cst2 /g(k1)/ cst12 /cp(k6,k1)
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst6  /icomp,istct,iphct,icp
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
 
      do 1 i = 1, icp
1        u(i) = 1.d30
c                             identify stable phase, save id
c                             in idc, and chemical potential in u.
      do 10 i = 1, icp
         iu = i1ct(i)
 
         do 20 j = 1, iu
            uphi = g(id1(i,j)) / cp(i,id1(i,j))
            if (uphi.ge.u(i)) goto 20
            u(i) = uphi
            id1s(i) = id1(i,j)
20       continue
 
10    continue
 
      end
 
      subroutine simpl2
c----------------------------------------------------
c simpl2 determines the vertices of the minimum free energy
c surface for binary (sub)systems. the vertices are returned in the
c array id2v.
c----------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j1=50)
 
      integer iid2(l6)
 
      common/ cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi,
     *               iiphi,iflg1
     *      / cst2 /g(k1)/ cst12 /cp(k6,k1)/ cst3 /x(k1,k6)
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)/ cst48 /u(k6)
     *      / cst6 /icomp,istct,iphct,icp
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
     *      / cst52 /hcp,h,hd,id,jd,kd,ld,md(3),ibc1,ibc2,ic(5)
 
      hcp=2
c                            test for phases below the composant g-x
c                            plane. save metastable phase id's in iidb
      do 30 i=1,i2s
         ibc1=i2c(1,i)
         ibc2=i2c(2,i)
         ict=1
         if (i2ct(i).eq.0) goto 90
c                            true if only phases are composants
         ib=i2ct(i)
         do 40 j =1,ib
            if (g(id2(i,j)) - cp(ibc1,id2(i,j)) * u(ibc1)
     *                      - cp(ibc2,id2(i,j)) * u(ibc2)
     *                                      .gt.0.d0) goto 40
            ict=ict+1
            iid2(ict)=id2(i,j)
40    continue
c                       now find the stable vertices for the same binary
90    id2v(i,1)=id1s(ibc1)
c                       ict=1, composants are the only stable phases.
      if (ict.ne.1) goto 10
c                            only composants stable:
         id2v(i,2)=id1s(ibc2)
         ibv=2
      goto 30
c                            some binary phases are stable:
10    ict=ict+1
      iid2(1)=id1s(ibc1)
      iid2(ict)=id1s(ibc2)
c                            determine stable binary phases:
      ibv=1
      hd=id1s(ibc1)
60    id=id1s(ibc2)
      call abload (*9000)
c                            begin testing phases:
      do 50 j =2,ict
         jd=iid2(j)
         if (x(hd,ibc2).ge.x(jd,ibc2)) goto 50
         if (g(jd)-cp(ibc1,jd)*b(1)-cp(ibc2,jd)*b(2).gt.1.d-05)
     *                                                       goto 50
            id=jd
            call abload (*9000)
50    continue
c                            at end of join?
      ibv=ibv+1
      if (ibv.gt.j1) call error (188,x(1,1),j1,'SIMPL2')
 
      id2v(i,ibv)=id
      if (x(id,ibc2).eq.1.d0) goto 30
      hd=id
      goto 60
c                            done
30    i2vct(i)=ibv
 
      return
 
9000  call error (19,x(1,1),id,'SIMPL2')
 
      end
 
      subroutine simpl3 (iflg4)
c----------------------------------------------------------------------
c simpl3 is a subprogram which computes the stable phase configu
c ations of ternary (sub)systems.  the algorithm is described in
c detail in program documentation elsewhere.  several versions of
c simpl3 are available, this version is designed for efficient
c computation in systems with binary solutions.
c
c   referenced by: bossop
c   references to: abload,assd3,aspsb,bound3,dgph3,itest3,sidc
c                  slopes
c   input arrays:     g,cp,u,idc,id2v,ibvct,itc,itb
c   output arrays:    id3f,i3fct,i3vct,id3v
c   temporary arrays: itbh,ivchk,idtpst,idpsf,ibin
c----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
      parameter (j1=50,j2=400,j6=200)
 
      integer ibin(3*j1)
 
      common/ cst48 /u(k6)/ cst46 /ipsf,idpsf(j8,k6)
     *      / cst2 /g(k1)/ cst3 /x(k1,k6)/ cst12 /cp(k6,k1)
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s/ cst18 /icct(k1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst19 /hi,iqth(4),ivchk(k1)
     *      / cst55 /s(6),sign,ipc(6),jcp,kcp
     *      / cst45 /itbh(3),itt,i3vcth,i3fcth
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct/ cst6  /icomp,istct,iphct,icp
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
     *      / cst52 /hcp,h,hd,id,jd,kd,ld,md(3),
     *               itc1,itc2,itc3,itc4,ic(3)
 
      do 1 i = 3, iphct
1        ivchk(i) = 0
 
      hcp = 3
      jcp = 2
      kcp = 1
 
      do 500 h = 1, i3s
c                                for each ternary subsystem
c                                initialize simple variables
      itc1=i3c(1,h)
      itc2=i3c(2,h)
      itc3=i3c(3,h)
 
      i3fcth=0
      itt=i3ct(h)
      i3vcth=0
      ipsf=0
      itpsf=1
      iflg4=0
      hgo=0
      iflg2=0
      ibt=0
c                                determine number of binary phases
c                                in the appropriate ternary:
      do 2 i=1, 3
c                                get the id of the binary:
         itbh(i)=i23(i,h)
         itbd=itbh(i)
c                                get the number of phases in the binary:
         ib=i2vct(itbd)-1
c                                no binary phases if ib=1
         goto (3),ib
         iflg2=iflg2+1
c                                save ids of binary phases in ibin
         do 4 k = 2,ib
           ibt=ibt+1
4          ibin(ibt)=id2v(itbd,k)
 
3        ibt=ibt+1
 
2        ibin(ibt)=id1s(i3c(i,h))
c-----------------------------------------------------------------------
c                                 iflg2 = number of binaries with
c                                 binary vertices, if ne 0 a composant
c                                 plane can't be stable, so skip test,
c                                 goto 250.
      if (iflg2.ne.0) goto 250
c                                 no ternary phases? goto 270 cause
c                                 the composant plane is table, this is
c                                 a lot of code for nothing.
      if (itt.eq.0) goto 270
c                                 ternary phases, test if any are stable
c                                 with respect to the composants:
      do 280 i=1,itt
         k = id3(h,i)
         if (g(k) - cp(itc1,k)*u(itc1)
     *            - cp(itc2,k)*u(itc2)
     *            - cp(itc3,k)*u(itc3) .lt. 0.d0) goto 250
280   continue
c                                 if it gets to here, no ternary and no
c                                 binary phases are stable.
270   hd=id1s(itc1)
      id=id1s(itc2)
      jd=id1s(itc3)
 
      call assd3 (itpsf ,hgo)
 
      goto 271
c----------------------------------------------------------------------
c                             identify ternary facets with 2
c                             vertices on a binary:
250   do 10 i=1, 3
         idbt=itbh(i)
         ib=i2vct(idbt)-1
         do 20 j =1,ib
c                             set the binary vertices:
            hd=id2v(idbt,j)
            id=id2v(idbt,j+1)
c                             1st test the opposite composant, this is
c                             sloppy
            do 21 k = 1, 3
21             if ( i2c(1,idbt).ne.i3c(k,h) .and.
     *              i2c(2,idbt).ne.i3c(k,h) ) goto 19
 
19          jd = id1s(i3c(k,h))
            zpot = k
c                             test if a vertice is contained:
            if (i3vcth.eq.0) goto 25
            call bound3
            do 35 k = 1, i3vcth
               if (itest3(id3v(h,k)).ne.0) goto 35
               jd=id3v(h,k)
               call bound3
35             continue
c                             determine u's:
25    call abload (*9000)
c                             begin testing against binaries:
      do 30 k = 1, 3
      if (k.eq.i) goto 30
      itbd=itbh(k)
      jb=i2vct(itbd)-1
      goto (30),jb
      do 40 l=2,jb
 
40       call dgdegn (zpot,id2v(itbd,l))
 
30    continue
c                             begin testing ternary phases:
      if (itt.eq.0) goto 50
      do 60 k = 1, itt
 
60       call dgdegn (zpot,id3(h,k))
 
c                             assign new facet:
50    call assd3 (itpsf ,hgo)
20    continue
10    continue
c                             begin second search segment.
      hgo=1
100   if (itpsf.eq.ipsf) goto 271
      istart=itpsf
      iend=ipsf
      itpsf=ipsf
c                             generate new facets from unmatched
c                             tielines:
      do 90 i = istart, iend
      if (idpsf(i,3).eq.0) goto 90
c                             set vertices:
      hd=idpsf(i,1)
      id=idpsf(i,2)
      call slopc (i)
c                             set 1st test vertice:
      do 110 k = 1, 3
c                             try unary vertices:
c                             only one unary vertex is possible, and
c                             this can be determined geometrically:
         jd=id1s(i3c(k,h))
110      if (sidc(jd).lt.-1.d-05) goto 150
c                             no error should be possible here, but:
      call error (997,x(1,1),id,'SIMPL3')
c                             test for bound vertices, these must be
c                             ternary phases, this is also only
c                             geometric testing:
150   if (i3vcth.eq.0) goto 135
      call bound3
      do 130 k = 1, i3vcth
         kd=id3v(h,k)
         if (id.eq.kd.or.itest3(kd).ne.0) goto 130
         jd=id3v(h,k)
         call bound3
130   continue
c                             now start thermodynamic testing:
135   call abload (*9000)
c                             test against binary phases:
      do 160 k = 1, ibt
160      call checkd (ibin(k))
c                             test ternary phases:
      if (itt.eq.0) goto 190
      do 180 k = 1, itt
180      call checkd (id3(h,k))
c                             assign new facets and tielines:
190   call assd3 (i,hgo)
c                             end of loop on current pstf's
90    continue
c                             restart loop on new pstf's
      goto 100
c                             end of search segment.
271   i3fct(h)=i3fcth
      i3vct(h)=i3vcth
      if (i3vcth.gt.0) iflg4=iflg4+1
500   continue
      h=1
c                           output chemography if requested.
      return
 
9000  call error (19,x(1,1),id,'SIMPL3')
 
      end

      subroutine simpl4 (iflg4)
c----------------------------------------------------------------------
c simpl4 locates the quaternary facets of a minimum energy surface.

c   referenced by: bossop
c   references to: ab4, aspst, assd4, dgph4
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
      parameter (j1=50,j2=400,j3=900,j6=200,j7=500)
 
      common/ cst2  /g(k1)/ cst12 /cp(k6,k1)/ cst46 /ipsf,idpsf(j8,k6)
     *      / cst19 /hi,iqth(4),ivchk(k1)
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s/ cst48 /u(k6)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst52 /hcp,h,hd,id,jd,kd,ld,md(3),
     *               iqc1,iqc2,iqc3,iqc4,iqc5,iqc6,iqc7
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst56 /id4v(m3,j7),i4vct(m3),i4vcth/ cst18 /icct(k1)
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
      common/ cst93 /i24a(3,4,35),i24(6,35)
     *      / cst55 /s(6),sign,ipc(6),jcp,kcp

 
      hcp = 4
      jcp = 3
      kcp = 2

      do 500 h = 1, i4s
c                              initialization:
 
c                              vertice counters
      i4vct(h)=0
c                              put component id's into simple
c                              variables:
      iqc1=i4c(1,h)
      iqc2=i4c(2,h)
      iqc3=i4c(3,h)
      iqc4=i4c(4,h)
c                              identities of the four ternaries in
c                              the hth quaternary:
      do 5 i= 1, 4
5        iqth(i)=i34(i,h)
c                              counter for the number of stable
c                              quaternary phases:
      iqqs=0
c                              counter of total number of quaternary
c                              phases:
      iqq=i4ct(h)
c                              counter of quaternary facets:
      i4fcth=0
c                              counter of quaternary vertices:
      i4vcth=0
c                              counter of pseudo-ternary planes:
      ipsf=0
      itpsf=1
c                              flag and counter for aspst in the 1st
c                              search segment.
      hgo=0
      ione=1
c------------------------------------------------------------------
c                              preliminary and mostly superfluous
c                              stuff:
      if (iqq.eq.0) goto 30
c                              eliminate quaternary phases
c                              metastable with respect to composants:
      do 10 i=1,iqq
 
         if (g(id4(h,i))-cp(iqc1,id4(h,i))*u(iqc1)
     *                  -cp(iqc2,id4(h,i))*u(iqc2)
     *                  -cp(iqc3,id4(h,i))*u(iqc3)
     *                  -cp(iqc4,id4(h,i))*u(iqc4).gt.0.0d0) goto 10
 
         iqqs=iqqs+1
 
10       continue
 
      if (iqqs.ne.0) goto 30
c                              no quaternary phases stable, if no
c                              ternary phases are stable (iflg4=0)
c                              then the only facet is defined by
c                              the composants.
      if (iflg4.ne.0) goto 30
      i4fcth=1
      do 25 i= 1, 4
25       id4f(h,1,i)=id1s(i4c(i,h))
      goto 490
c-----------------------------------------------------------------
30    do 100 i= 1, 4
c                             identify facets with 3 vertices in a
c                             ternary subsystem:
         hi = i34(i,h)
         ii = i3fct(hi)
 
         do 40 j = 1, ii
c                             set ternary vertices:
         hd = id3f(hi,j,1)
         id = id3f(hi,j,2)
         jd = id3f(hi,j,3)
c                             set 4th vertex as the opposite (absent)
c                             composant:
         do 43 ll = 1, 4
            do 1 kk = 1, 3
1              if (i3c(kk,hi).eq.i4c(ll,h)) goto 43
            zpot = ll
            kd = id1s(i4c(ll,h))
            goto 44
43       continue
c                             determine u's:
44       call abload (*9000)
c                             begin trying binary vertices
         do 41 k = 1, 3
c                             identify the binary, it must contain
c                             composant
            idb = i24a(k,i,h)
            ib = i2vct(idb)-1
            goto (41),ib
c                             start testing:
            do 42 l=2,ib
 
42             call dgdegn (zpot,id2v(idb,l))
 
41       continue
c                             begin trying ternary verices:
      do 50 k = 1, 4
 
         if (k.eq.i) goto 50
         idb = i34(k,h)
         ib = i3vct(idb)
         if (ib.eq.0) goto 50
 
         do 60 l = 1, ib
 
60          call dgdegn (zpot,id3v(idb,l))
 
50    continue
c                             if no quaternary phases then the
c                             established facet is stable, else test
c                             quaternary phases:
      if (iqqs.eq.0) goto 90
 
      do 70 k = 1, iqq
 
70       call dgdegn (zpot,id4(h,k))
 
c                             save the new facet:
90    call assd4 (ione,hgo)
c                             save the new pseudo-ternary planes:
40    continue
100   continue
c                             begin 2nd search segment.
      hgo=1
      if (ipsf.eq.0) goto 490
c
160   if (ipsf.eq.itpsf) goto 490
      istart=itpsf
      iend=ipsf
      itpsf=iend
c                             identify new facets generated from
c                             the pstp's:
      do 150 i = istart, iend
         if (idpsf(i,1).eq.0) goto 150
c                             set vertices:
         hd=idpsf(i,1)
         id=idpsf(i,2)
         jd=idpsf(i,3)
c                             determine the equation of the pseudo-
c                             ternary plane hd-id-jd:
         call slopc (i)
c                             set first test vertice:
         do 130 jj = 1, 4
            kd = id1s(i4c(jj,h))
            if (sidc(kd).lt.-1.d-05) then
               call abload (*130)
               goto 180
            end if 
130      continue 

         call error (196,1.d0,k6,'SIMPL4')
c                             call boundc here to establish
c                             compositional limits of the trial
c                             configuration. then use itestc to
c                             test for bounded stable phases. this
c                             test loop is only useful if the
c                             system contains quaternary compounds.
c                             begin testing:
c                             try binary vertices:
180      do 81 k =  1, 6
            idb=i24(k,h)
            ib=i2vct(idb)
            do 82 l=2,ib
82             call checkd (id2v(idb,l))
81       continue
c                             try ternary vertices:
         do 200 j = 1, 4
            idb = i34(j,h)
            ib = i3vct(idb)
            if (ib.eq.0) goto 200
            do 210 k = 1, ib
210            call checkd (id3v(idb,k))
200      continue
c                             now test against quaternary phases:
      if (iqqs.eq.0) goto 230
      do 220 j = 1, iqq
220      call checkd (id4(h,j))
c                              assign new facet:
230   call assd4 (i,hgo)
c                              assign pstp to a temporary storage
c                              array.
150   continue
c                              save permanent counters:
      goto 160
490   i4vct(h)=i4vcth
      i4fct(h)=i4fcth
      goto 500
 
9000  call error (19,1.d0,id,'SIMPL4')
 
500   continue
 
      end
 
      subroutine simpl5
c--------------------------------------------------------------------
c simpl5 locates the quinary facets of a free energy surface.
 
c   referenced by: bossop
c   references to: abload,assd5,dgph5,sidc,slopc
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k6=7,l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
      parameter (j1=50,j2=400,j3=900,j4=1200,j6=200,j7=500,h3=250)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ cst55 /s(6),sign,ipc(6),jcp,kcp
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst52 /hcp,l,hd,id,jd,kd,ld,md,zd(2),ic(7)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst56 /id4v(m3,j7),i4vct(m3),i4vcth
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst98 /id5v(h3,m3),i5vct(m3),i5vctl
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
      common/ cst94 /i25a(4,5,21),i25(10,21),i35a(6,5,21),i35(10,21)
c-----------------------------------------------------------------------
c                              number of components:
      hcp = 5
      jcp = 4
      kcp = 3
c                              initialization:
      do 500 l = 1, i5s
c                              components:
      do 1 i = 1, 5
1        ic(i)=i5c(i,l)
c                              pseudo-quaternary plane counters:
      ipsf=0
      itpsf=1
c                              temprorary facet counter
      i5fctl = 0
      i5vctl = 0
c                              flags
      hgo=0
c                              begin first search segment.  this
c                              segment locates facets which have
c                              4 vertices in a common quaternary
c                              subsystem.
      do 10 h= 1, 5
c                              quat identity:
         jd4 = i45(h,l)
 
         hh=i4fct(jd4)
 
         do 20 i=1,hh
c                              set quaternary vertices:
            hd=id4f(jd4,i,1)
            id=id4f(jd4,i,2)
            jd=id4f(jd4,i,3)
            kd=id4f(jd4,i,4)
c                              set fifth vertex as the opposite
c                              composant:
            do 43 ll = 1, 5
               do 41 kk = 1, 4
41                if (i4c(kk,jd4).eq.ic(ll)) goto 43
               ld = id1s(ic(ll))
               zpot = ll
               goto 44
43          continue
c                              determine chemical potentials:
44          call abload (*9000)
c                             test stability relative to binary
c                             vertices:
            do 30 j = 1, 4
               idb = i25a(j,h,l)
               ib = i2vct(idb)-1
               goto (30),ib
c                              start test loop:
               do 40 k = 2, ib
 
                  call dgdegn (zpot,id2v(idb,k))
 
40             continue
c                              end test loop.
30          continue
c                              test stability relative to ternary
c                              vertices:
            do 50 j = 1, 6
               idb = i35a(j,h,l)
               ib = i3vct(idb)
               if (ib.eq.0) goto 50
c                              start test loop:
               do 60 k = 1, ib
 
                  call dgdegn (zpot,id3v(idb,k))
 
60             continue
c                              end test loop.
50          continue
c                              test stability relative to quaternary
c                              vertices:
            do 70 j = 1, 5
               kd4 = i45(j,l)
               if (kd4.eq.jd4) goto 70
               ib=i4vct(kd4)
               if (ib.eq.0) goto 70
c                              start test loop:
               do 80 k = 1, ib
 
                  call dgdegn (zpot,id4v(kd4,k))
 
80             continue
c                              end test loop.
c                              test stability relative to quinary
c                              phases:
70          continue
 
            ib = i5ct(l)
            if (ib.eq.0) goto 100
c                              start test loop:
            do 90 j = 1, ib
 
                  call dgdegn (zpot,id5(l,j))
 
90          continue
c                              end test loop.
100         call assd5 (itpsf,hgo,jd4)
 
20       continue
10    continue
c                              end first segment
c
c                              begin second search segment:
      hgo=1
110   if (ipsf.eq.itpsf) goto 250
      istart=itpsf
      iend=ipsf
      itpsf=ipsf
c                             generate new facets from the psq's:
      do 120 i = istart, iend
         if (idpsf(i,1).eq.0) goto 120
c                             set vertices:
         hd=idpsf(i,1)
         id=idpsf(i,2)
         jd=idpsf(i,3)
         kd=idpsf(i,4)
c                             determine the plane in simplex
c                             composition space including hd-
c                             id-jd-kd.
         call slopc(i)
c                             set first test vertice:
         do 130 jj = 1, 5
            ld = id1s(ic(jj))
            if (sidc(ld).lt.-1.d-05) then
               call abload (*130)
               goto 140
            end if 
130      continue 
         
         call error (886,sign,id,'SIMPL5')
c                             call boundc here to establish
c                             compositional limits of the trial
c                             configuration. then use itestc to
c                             test for bounded stable phases. this
c                             test loop is only useful if the
c                             system contains quinary compounds.
c                             test against composants:
140      do 150 j = jj, 5
150         call checkd (id1s(ic(j)))
c                             test against binaries:
         do 160 j =1,10
            idb = i25(j,l)
            ib=i2vct(idb)-1
            goto (160),ib
            do 170 k = 2,ib
170            call checkd (id2v(idb,k))
160      continue
c                             test against ternaries:
         do 180 j =1,10
            idb = i35(j,l)
            ib=i3vct(idb)
            if (ib.eq.0) goto 180
            do 190 k = 1, ib
190            call checkd (id3v(idb,k))
180      continue
c                             test against quaternaries:
         do 200 j = 1, 5
            idb = i45(j,l)
            ib=i4vct(idb)
            if (ib.eq.0) goto 200
            do 210 k = 1, ib
210            call checkd (id4v(idb,k))
200      continue
c                             test against quinaries:
         ib = i5ct(l)
         if (ib.eq.0) goto 120
         do 230 j = 1, ib
230         call checkd (id5(l,j))
 
120   call assd5 (i,hgo,jd4)      
c                             end second search segment.
      goto 110
250   i5fct(l) = i5fctl
      i5vct(l) = i5vctl
 
500   continue
 
      return
 
9000  call error (19,sign,id,'SIMPL5')
 
      end
 
      subroutine simpl6
c--------------------------------------------------------------------
c simpl6 locates the hexary facets of a free energy surface.
 
c   referenced by: bossop
c   references to: abload,assd6,dgph6,sidc,slopc
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k6=7,h4=20)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100,h3=250)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
      parameter (j1=50,j2=400,j3=900,j4=1200,j5=2400,j6=200,j7=500)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ cst55 /s(6),sign,ipc(6),jcp,kcp
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst52 /hcp,l,hd,id,jd,kd,ld,md,zd,i0,ic(7)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst56 /id4v(m3,j7),i4vct(m3),i4vcth
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst98 /id5v(h3,m3),i5vct(m3),i5vctl
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl/ cst14 /i2s,i3s,i4s,i5s,i6s,i7s
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
      common/ cst94 /i25a(4,5,21),i25(10,21),i35a(6,5,21),i35(10,21)
     *      / cst95 /i16a(6,7),i26a(5,6,7),i26(15,7),i36a(10,6,7),
     *               i36(20,7),i46a(10,6,7),i46(15,7)
c-----------------------------------------------------------------------
c                              initialization:

c                              number of components:
      hcp = 6
      jcp = 5
      kcp = 4

      do 600 l = 1, i6s
c                              components:
      do 1 i = 1, 6
1        ic(i)=i6c(i,l)
c                              pseudo-quaternary plane counters:
      ipsf=0
      itpsf=1
c                              temprorary facet counter
      i6fctl = 0
      i6vctl = 0
c                              flags
      hgo=0
c                              begin first search segment.  this
c                              segment locates facets which have
c                              5 vertices in a common quinary
c                              subsystem.
      do 10 h= 1, 6
c                              quin identity:
         jd5 = i56(h,l)
 
         hh=i5fct(jd5)
 
         do 20 i=1,hh
c                              set quinary vertices:
            hd=id5f(1,i,jd5)
            id=id5f(2,i,jd5)
            jd=id5f(3,i,jd5)
            kd=id5f(4,i,jd5)
            ld=id5f(5,i,jd5)
c                              set 6th vertex as the opposite
c                              composant:
            do 43 ll = 1, 6
               do 41 kk = 1, 5
41                if (i5c(kk,jd5).eq.ic(ll)) goto 43
               md = id1s(ic(ll))
               zpot = ll
               goto 44
43          continue
c                              determine chemical potentials:
44          call abload (*9000)
c                             test stability relative to binary
c                             vertices:
            do 30 j = 1, 5
               idb = i26a(j,h,l)
               ib = i2vct(idb)-1
               goto (30),ib
c                              start test loop:
               do 40 k = 2, ib
 
 40               call dgdegn (zpot,id2v(idb,k))
 
c                              end test loop.
30          continue
c                              test stability relative to ternary
c                              vertices:
            do 50 j = 1, 10
               idb = i36a(j,h,l)
               ib = i3vct(idb)
               if (ib.eq.0) goto 50
c                              start test loop:
               do 60 k = 1, ib
 
                  call dgdegn (zpot,id3v(idb,k))
 
60             continue
c                              end test loop.
50          continue
c                              test stability relative to quaternary
c                              vertices, in the ten quaternaries of
c                              the hexary l, that are not included
c                              in quinary h.
            do 70 j = 1, 10
               idb = i46a(j,h,l)
               ib=i4vct(idb)
               if (ib.eq.0) goto 70
c                              start test loop:
               do 80 k = 1, ib
 
                  call dgdegn (zpot,id4v(idb,k))
 
80             continue
c                              end test loop.
70          continue
c                              test stability relative to quinary
c                              vertices:
            do 90 j = 1, 6
               idb = i56(j,l)
               if (idb.eq.jd5) goto 90
               ib = i5vct(idb)
               if (ib.eq.0) goto 90
c                              start test loop:
               do 95 k = 1, ib
 
                  call dgdegn (zpot,id5v(k,idb))
 
95             continue
c                              end test loop.
90          continue
c                              lastly test against the phases
c                              of the present hexary:
            ib = i6ct(l)
            if (ib.eq.0) goto 100
c                              start test loop:
            do 96 j = 1, ib
 
                  call dgdegn (zpot,id6(l,j))
 
96          continue
c                               assign the facet hd-id-jd-kd-ld-md
 
100         call assd6 (itpsf,hgo,jd5)
 
20       continue
10    continue
c                              end first segment
 
c                              begin second search segment:
      hgo=1
110   if (ipsf.eq.itpsf) goto 250
      istart=itpsf
      iend=ipsf
      itpsf=ipsf
c                             generate new facets from the psq's:
      do 120 i = istart, iend
         if (idpsf(i,1).eq.0) goto 120
c                             set vertices:
         hd=idpsf(i,1)
         id=idpsf(i,2)
         jd=idpsf(i,3)
         kd=idpsf(i,4)
         ld=idpsf(i,5)
c                             determine the plane in simplex
c                             composition space including hd-
c                             id-jd-kd.
         call slopc (i)
c                             set first test vertice:
         do 130 jj = 1, 6
            md=id1s(ic(jj))
            if (sidc(md).lt.-1.d-05) then
               call abload (*130)
               goto 140
            end if 
130      continue 
         
         call error (886,sign,md,'SIMPL6')
c                             test against composants:
140      do 150 j = jj, 6
150         call checkd (id1s(ic(j)))
c                             test against binaries:
         do 160 j = 1, 15
 
            idb = i26(j,l)
            ib=i2vct(idb)-1
            goto (160),ib
            do 170 k = 2,ib
170            call checkd (id2v(idb,k))
160      continue
c                             test against ternaries:
         do 180 j = 1, 20
 
            idb = i36(j,l)
            ib=i3vct(idb)
            if (ib.eq.0) goto 180
 
            do 190 k =  1, ib
190            call checkd (id3v(idb,k))

180      continue
c                             test against quaternaries:
         do 200 j = 1, 15
 
            idb = i46(j,l)
            ib=i4vct(idb)
            if (ib.eq.0) goto 200
 
            do 210 k = 1, ib
210            call checkd (id4v(idb,k))

200      continue
c                             test against quinaries:
         do 220 j = 1, 6
 
            idb = i56(j,l)
            ib = i5vct(idb)
            if (ib.eq.0) goto 220
 
            do 230 k = 1, ib
  
230            call checkd (id5v(k,idb))

220   continue
c                             test against hexaries:
      ib = i6ct(l)
 
      if (ib.eq.0) goto 120
 
      do 240 j = 1, ib
 
240      call checkd (id6(l,j))
c                             assign hd-id-jd-kd-ld-md
120   call assd6 (i,hgo,jd5)
c                             end second search segment.
      goto 110
250   i6fct(l) = i6fctl
      i6vct(l) = i6vctl
 
600   continue
 
      return
 
9000  call error (19,sign,id,'SIMPL6')
 
      end
 
      subroutine simpl7
c--------------------------------------------------------------------
c simpl7 locates the septary facets of a free energy surface.
 
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 (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100,j8=4800)
      parameter (j1=50,j2=400,j3=900,j4=1200,j5=2400,j6=200,j7=500)
      parameter (h4=20,j9=1200,k6=7,h3=250)
 
      common/ cst46 /ipsf,idpsf(j8,k6)/ cst96 /id7f(k6,j9),i7fct
     *      / cst15 /id1s(k6),id2v(m1,j1),i2vct(m1)
     *      / cst16 /id3f(m2,j2,3),id3v(m2,j6),i3fct(m2),i3vct(m2)
     *      / cst52 /hcp,l,hd,id,jd,kd,ld,md,zd,i0,ic(7)
     *      / cst47 /id4f(m3,j3,4),i4fct(m3),i4fcth,igo
     *      / cst56 /id4v(m3,j7),i4vct(m3),i4vcth
     *      / cst58 /id5f(5,j4,m4),i5fct(m4),i5fctl
     *      / cst98 /id5v(h3,m3),i5vct(m3),i5vctl
     *      / cst55 /s(6),sign,ipc(6),jcp,kcp
     *      / cst97 /id6f(6,j5,m5),id6v(h4,m5),i6fct(m5),i6vct(m5),
     *               i6fctl,i6vctl
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
      common/ cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
      common/ cst94 /i25a(4,5,21),i25(10,21),i35a(6,5,21),i35(10,21)
     *      / cst95 /i16a(6,7),i26a(5,6,7),i26(15,7),i36a(10,6,7),
     *               i36(20,7),i46a(10,6,7),i46(15,7)
     *      / cst101 /i27a(6,7),i37a(15,7),i47a(20,7),i57a(15,7)
c-----------------------------------------------------------------------
c                              initialization:
 
c                              number of components:
      hcp = 7
      jcp = 6
      kcp = 5
c                              components:
      do 1 i = 1, 7
1        ic(i) = i
c                              pseudo-quaternary plane counters:
      ipsf=0
      itpsf=1
c                              temprorary facet counter
      i7fct = 0
c                              flags
      hgo=0
c                              begin first search segment.  this
c                              segment locates facets which have
c                              6 vertices in a common hexary
c                              subsystem.
      do 10 h = 1, 7
c                              hexary identity = h
         hh = i6fct(h)
         jj = 8 - h
 
         do 20 i = 1, hh
c                              set quinary vertices:
            hd=id6f(1,i,h)
            id=id6f(2,i,h)
            jd=id6f(3,i,h)
            kd=id6f(4,i,h)
            ld=id6f(5,i,h)
            md=id6f(6,i,h)
c                              set 7th vertex as the opposite
c                              composant:
            zd = id1s(jj)
c                              determine chemical potentials:
            call abload (*9000)
c                             test stability relative to binary
c                             vertices:
            do 30 j = 1, 6
               idb = i27a(j,h)
               ib = i2vct(idb)-1
               goto (30),ib
c                              start test loop:
               do 40 k = 2, ib
 
 40               call dgdegn (jj,id2v(idb,k))
 
c                              end test loop.
30          continue
c                              test stability relative to ternary
c                              vertices:
            do 50 j = 1, 15
               idb = i37a(j,h)
               ib = i3vct(idb)
               if (ib.eq.0) goto 50
c                              start test loop:
               do 60 k = 1, ib
 
 60               call dgdegn (jj,id3v(idb,k))
 
c                              end test loop.
50          continue
c                              test stability relative to quaternary
c                              vertices:
            do 70 j = 1, 20
 
               idb = i47a(j,h)
               ib=i4vct(idb)
 
               if (ib.eq.0) goto 70
c                              start test loop:
               do 80 k = 1, ib
 
 80               call dgdegn (jj,id4v(idb,k))
 
c                              end test loop.
70          continue
c                              test stability relative to quinary
c                              vertices:
            do 90 j = 1, 15
 
               idb = i57a(j,h)
               ib = i5vct(idb)
 
               if (ib.eq.0) goto 90
c                              start test loop:
               do 95 k = 1, ib
 
 95               call dgdegn (jj,id5v(k,idb))
 
c                              end test loop.
90          continue
 
c                              test stability relative to hexary
c                              vertices:
            do 97 j = 1, 7
 
               if (j.eq.h) goto 97
               ib = i6vct(j)
               if (ib.eq.0) goto 97
c                              start test loop:
               do 98 k = 1, ib
 
 98               call dgdegn (jj,id6v(k,j))
 
c                              end test loop.
97          continue
c                              lastly test against the phases
c                              of the septary:
            ib = i7ct
 
            if (ib.eq.0) goto 100
c                              start test loop:
            do 96 j = 1, ib
 
 96            call dgdegn (jj,id7(j))
 
c                               assign the facet hd-id-jd-kd-ld-md
c
100         call assd7 (itpsf,hgo,h)
 
20       continue
10    continue
c                              end first segment
c
c                              begin second search segment:
      hgo=1
110   if (ipsf.eq.itpsf) goto 99
      istart=itpsf
      iend=ipsf
      itpsf=ipsf
c                             generate new facets from the psq's:
      do 120 i = istart, iend
         if (idpsf(i,1).eq.0) goto 120
c                             set vertices:
         hd=idpsf(i,1)
         id=idpsf(i,2)
         jd=idpsf(i,3)
         kd=idpsf(i,4)
         ld=idpsf(i,5)
         md=idpsf(i,6)
c                             determine the plane in simplex
c                             composition space including hd-
c                             id-jd-kd.
         call slopc (i)
c                             set first test vertice:
         jj = 1
130      zd = id1s(jj)
            jj = jj + 1
            if (sidc(zd).lt.-1.d-05) then
               call abload (*141) 
               goto 140 
            end if 
141         if (jj.gt.7) call error (882,sign,id,'SIMPL7')
         goto 130
c                             test against composants:
140      do 150 j = jj, 7
150         call checkd (id1s(j))
c                             test against binaries:
         do 160 j = 1, 21
 
            ib=i2vct(j)-1
            goto (160),ib
 
            do 170 k = 2,ib
170            call checkd (id2v(j,k))
160      continue
c                             test against ternaries:
         do 180 j = 1, 35
 
            ib=i3vct(j)
            if (ib.eq.0) goto 180
 
            do 190 k =  1, ib
190            call checkd (id3v(j,k))
180      continue
c                             test against quaternaries:
         do 200 j = 1, 35
 
            ib=i4vct(j)
            if (ib.eq.0) goto 200
 
            do 210 k = 1, ib
210            call checkd (id4v(j,k))
200      continue
c                             test against quinaries:
         do 220 j = 1, 21
 
            ib = i5vct(j)
            if (ib.eq.0) goto 220
 
            do 230 k = 1, ib
230            call checkd (id5v(k,j))
220      continue
c                             test against hexaries
         do 240 j = 1, 6
 
            ib = i6vct(j)
            if (ib.eq.0) goto 240
 
            do 245 k = 1, ib
245            call checkd (id6v(k,j))
240     continue
c                             test against septaries
 
        if (i7ct.eq.0) goto 120
 
        do 250 j = 1, i7ct
250        call checkd (id7(j))
c                             assign hd-id-jd-kd-ld-md-zd
120   call assd7 (i,hgo,h)

      goto 110
c                             end second search segment.
9000  call error (19,sign,id,'SIMPL7')
 
99    end
 
      subroutine checkd (zt)

      implicit double precision (a-g,o-y),integer (h-n,z)

      common / cst52 /hcp,h,id(8),ic(7)

      do 100 i = 1, hcp 
100      if (id(i).eq.zt) goto 99

      if (sidc(zt).gt.-1.d-05) goto 99

      goto (10,10,30,40,50,60,70), hcp

10    call error (999,1.d0,hcp,'CHECKD')

30    if (dgph3(zt).gt.-1.d-05) goto 99
      goto 90
40    if (dgph4(zt).gt.-1.d-05) goto 99
      goto 90
50    if (dgph5(zt).gt.-1.d-05) goto 99
      goto 90
60    if (dgph6(zt).gt.-1.d-05) goto 99
      goto 90
70    if (dgph7(zt).gt.-1.d-05) goto 99

90    ztold = id(hcp)
      id(hcp) = zt 

      call abload (*110)

      goto 99
c                           restore old vertex
110   id(hcp) = ztold

      call abload (*9000)
      
      goto 99

9000  call error (19,1.d0,zt,'CHECKD')

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                                  output chemography if requested.
c                                  this option is archaic with newhld.
      goto (10),ichem
      call outchm
c                                  load assemblages into the array id7f
c                                  with counter i7fct (cst58).
10    call allone
c                                  write tier header to print file:
      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                              sort phase into subcomposition
       call sort (ier)
c                              if ier = 1, the phase consists
c                              entirely of saturated components:
       if (ier.ne.0) call satsrt 
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 sort (ier)
c---------------------------------------------------------------------
c sort determines if a phase in the system defined in input may
c be classified into a subsystem. i indexes the phase in the arrays.
c
c   referenced by: main
c   references to: none
c   input arrays:     cp,ibc,itc,iqac
c   output arrays:    x,icct,iuct,ibct,itct,iqact,idu,idb
c                     idt,idq,idqi

c   on return if ier = 1, the phase is not in the thermodynamic
c   composition space.

c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)
 
      parameter (k1=1200,k6=7)
      parameter (l5=400,l6=400,l7=400,l8=200,l9=100)
      parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100)
 
      character*8 names
 
      common/ cst3  /x(k1,k6)/ cst18 /icct(k1)/ cst12 /cp(k6,k1)
     *      / cst6  /icomp,istct,i,icp/ cst8 /names(k1)
     *      / cst14 /i2s,i3s,i4s,i5s,i6s,i7s
      common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7)
     *      / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4),
     *               i6ct(m5),i7ct/ cst19 /hi,iqth(4),ivchk(k1)
      common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8),
     *               id5(m4,l9),id6(m5,h1),id7(h2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
c---------------------------------------------------------------------
      tot = 0.
      ier = 0 
      icct(i) = 0
 
      do 10 j = 1, icp
         if (cp(j,i).lt.1.d-8) cp(j,i) = 0.0
         if (cp(j,i).eq.0.0) goto 10
         icct(i) = icct(i) + 1
         tot = tot + cp(j,i)
10    continue
c                             if icct(i) = 0, the phase is
c                             not in thermodynamic composition space
      if (icct(i).eq.0) then
         ier = 1
         goto 99
      end if       
c                             compute compositions:
      do 15 j = 1, icp
15       x(i,j) = cp(j,i)/tot
c                             if icct<icp the phase is degenerate
      igo = icct(i)
      if (icp.eq.igo) goto 5
      goto (70,80,90,100,150,160), igo
 
5     goto (20,30,40,50,60,61,62), igo
c                             unary
20    i1ct(1) = i1ct(1) + 1
      if (i1ct(1).gt.l5) call error (183,tot,l5,'SORT')
      id1(1,i1ct(1)) = i
      return
c                             binary
30    i2ct(1) = i2ct(1) + 1
      if (i2ct(1).gt.l6) call error (184,tot,l6,'SORT')
      id2(1,i2ct(1)) = i
      return
c                             ternary
40    i3ct(1) = i3ct(1) + 1
      if (i3ct(1).gt.l7) call error (185,tot,l7,'SORT')
      id3(1,i3ct(1)) = i
      return
c                             quaternary
50    i4ct(1) = i4ct(1) + 1
      if (i4ct(1).gt.l8) call error (186,tot,l8,'SORT')
      id4(1,i4ct(1)) = i
      return
c                             quinary
60    i5ct(1) = i5ct(1) + 1
      if (i5ct(1).gt.l9) call error (187,tot,l9,'SORT')
      id5(1,i5ct(1)) = i
      return
c                             hexary
61    i6ct(1) = i6ct(1) + 1
      if (i6ct(1).gt.h1) call error (198,tot,h1,'SORT')
      id6(1,i6ct(1)) = i
      return
c                             septary
62    i7ct = i7ct + 1
      if (i7ct.gt.h2) call error (199,tot,h2,'SORT')
      id7(i7ct) = i
      return
c                             degenerate phases:
c
c                             unary subsystem
70    do 110 j = 1, icp
         if (cp(j,i).eq.0.d0) goto 110
         i1ct(j) = i1ct(j)+1
         if (i1ct(j).gt.l5) call error (183,tot,l5,'SORT')
         id1(j,i1ct(j)) = i
         return
110   continue
c                             binary subsystem
80    do 120 j = 1, i2s
         if ((cp(i2c(1,j),i).eq.0.d0).or.
     *       (cp(i2c(2,j),i).eq.0.d0)) goto 120
             i2ct(j) = i2ct(j) + 1
         if (i2ct(j).gt.l6) call error (184,tot,l6,'SORT')
         id2(j,i2ct(j)) = i
         return
120   continue
c                             ternary subsystem
90    do 130 j = 1, i3s
      if ((cp(i3c(1,j),i).eq.0.d0).or.
     *    (cp(i3c(2,j),i).eq.0.d0).or.
     *    (cp(i3c(3,j),i).eq.0.d0)) goto 130
      i3ct(j) = i3ct(j)+1
      if (i3ct(j).gt.l7) call error (185,tot,l7,'SORT')
      id3(j,i3ct(j)) = i
      return
130   continue
c                             quaternary subsystem
100   do 140 j = 1, i4s
      if ((cp(i4c(1,j),i).eq.0.d0).or.
     *    (cp(i4c(2,j),i).eq.0.d0).or.
     *    (cp(i4c(3,j),i).eq.0.d0).or.
     *    (cp(i4c(4,j),i).eq.0.d0)) goto 140
      i4ct(j) = i4ct(j) + 1
      if (i4ct(j).gt.l8) call error (186,tot,l8,'SORT')
      id4(j,i4ct(j)) = i
      return
140   continue
c                             quinary subsystem
150   do 170 j = 1, i5s
         do 171 k = 1, 5
171         if (cp(i5c(k,j),i).eq.0) goto 170
         i5ct(j) = i5ct(j) + 1
         if (i5ct(j).gt.l9) call error (187,tot,l9,'SORT')
         id5(j,i5ct(j)) = i
         return
170   continue
c                             hexary subsystem
160   do 180 j = 1, i6s
         do 181 k = 1, 6
181         if (cp(i6c(k,j),i).eq.0) goto 180
         i6ct(j) = i6ct(j) + 1
         if (i6ct(j).gt.h1) call error (198,tot,h1,'SORT')
         id6(j,i6ct(j)) = i
         return
180   continue
 
      call error (14,tot,i,'SORT')
 
99    end
 
      subroutine satsrt 
c---------------------------------------------------------------------
c routine to sort pseudocompounds consisting entirely of saturated
c components.
c---------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
  
      parameter (k1=1200,k5=12,h5=7,h6=500)
 
      common/ cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst20 /comps(k1,k5)
     *      / cst42 /ic(k5),idbase/ cst6 /icomp,istct,iphct,icp
 
      do 10 j = isat, 1, -1

         idc = icp + j

         if (comps(iphct,idc).ne.0.0) then

            isct(j) = isct(j) + 1

            if (isct(j).gt.h6) call error (17,comps(1,1),h6,'SATSRT')

            if (iphct.gt.k1) call error (180,comps(1,1),k1,'SATSRT')

            ids(j,isct(j)) = iphct

            goto 99

         end if

10    continue

99    end
