
c Please do not distribute any part of this source.
 
c Copyright (c) 1998 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.sep.10.1998  *
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     o    optional (io5) data file for graphics applications.
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    input for icopt = 5 calculations.
c   n8          unused.
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-----------------------------------------------------------------------
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 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                 testit (4)      sreset (2)      asschk (4)
c itestc   (4)    loadit (1)      flipit (2)      wway   (2)
c                 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-z),integer (h-n)

      character*80 record,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  h5 - maximum number of component saturation constraints
c  h6 - maximum number of composant phases for a component saturation
c       constraint
c  j9 - maximum number of divariant assemblages. 
c  k1 - maximum number of phases.
c  k2 - maximum number of univariant or invariant equilibria
c  k4 - number of parameters in thermodynamic equation of state
c  k5 - maximum number of components
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  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: increasing parameter K5 requires changes to the following
c format statements:
c                     1000 in routine bulkck
c                     1150 in routine input1
c                     4000 in main in program build
c-----------------------------------------------------------------------
      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,k8=k5+2,l2=5)
 
      common/ cst5  /v(l2),tr,pr,r,ps/ cst9  /vmax(l2),vmin(l2),dv(l2)
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct/ csta8 /title
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst81 /ipchk(k2),ip2(k2),icp2,ifake/ cst80 /irchk(k2)
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst31 /vn(k2,k7),irct,ird
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
c----------------------------------------------------------------------- 
c                            -------------------------------------------
c                             initialize invariant and univariant
c                             field counters:
      ipct = 0
      irct = 0

c      call money
c                             and flags:
      do i = 1, k2
         ipchk(i) = 0
      end do
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.
      call input1 
c                              read thermodynamic data on unit n2:
      call input2
c                              read data for solution phases on n9:
      call input9
c                              write header to graphics files (n4, n5):
      if (io4.ne.1) then 
         if (icopt.eq.5) then 
c                              bulk composition data to be output to n5
c                              computing polygon centroids: 
            call header (n5,icopt)
         else if (icopt.eq.1.and.jbulk.ne.0) then
c                              computing pseudo-section
            call header (n4,icopt)
            call header (n5,6)
         else if (icopt.eq.1) then 
c                              computing schreinemakers projection
            call header (n4,icopt)
         else 
c                              computing chemography or mixed-var.
            call outhed
         end if
      end if 
c                              initialize v(iv4) and v(iv5) (these
c                              variables are constant for each set).
c                              calculations.
      v(iv3) = vmin(iv3)
      v(iv4) = vmin(iv4)
      v(iv5) = vmin(iv5)
      if (icopt.eq.0.or.icopt.eq.5) then
c                              calculate composition phase diagrams:
         call chmcal
         goto 99
       else
c                              output title page to unit n3:
         call outtit (n3)
c                              read stability field assemblages on unit n1:
         if (icopt.eq.2) call error (32,v(1),k2,'MAIN')
      end if
c                              calculate phase diagram/stability field
      call newhld 
c                              summarize output for each tier
c                              on units n3 and n4 (n6):
      call outier
      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 (200),io4
      rewind n6
50    read (n6,1040,end=200) record
      write (n4,1040) record
      goto 50
c                              summarize calculations:
200   if (icopt.ne.3) call outlst
c                              done:
1040  format (a80)
1160  format (/,80('-'),/)
99    end
 
      subroutine chmcal
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (l2=5)
 
      common/ cst5   /v(l2),tr,pr,r,ps/ cst19 /jpoly
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst24  /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
c-----------------------------------------------------------------------
      call outtit (n3)

      jpoly = 0 
 
90    if (icopt.eq.0) then 
         read (n1,*,end=99) v
c                                  for old file formats:
         if (v(1).lt.0d0) goto 99
         jpoly = jpoly + 1
         write (*,1010) jpoly
      else 
c                                  must be in polygon mode:
         read (n7,*,end=99) jpoly,v(iv1),v(iv2)
         write (*,1000) jpoly
      end if 
c                                  compute phase diagram
      call gall
      call combin
c                                  generate permanent graphics file
c                                  if requested
      if (io4.ne.1) call outgrf
c                                  printer output:
      if (io3.ne.1) call outchm
      goto 90

1000  format (' Computing the assemblage in polygon ',i5)
1010  format (' Computing the compositional phase relations at',
     *        ' condition ',i2)
 
99    end
  
      subroutine input1 
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-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k5=12,h5=7,k7=k5+1,h6=500)
      parameter (l2=5,h8=100,h9=30,k0=25)
 
      character*8 fname*10,blank8,exname,n5name*14,n7name*14,
     *            cname*5,blank5*5,rname(k5)*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,y*1,zname*5
 
      common/ cst6  /icomp,istct,iphct,icp/ cst79 /isoct
     *      / csta7 /fname(h9) / cst24 /ipot,jv(l2),iv(l2)
     *      / cst81 /ipchk(k2),ip2(k2),icp2,ifake/ cst112 /buf(5)
     *      / csta2 /xname(k5),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(k2),isudo,ivar
     *      / cst31 /vn(k2,k7),irct,ird/ cst38 /iedge
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst39 /ichem
     *      / cst88 /vu(2,k1),jprct,jmct/ csta8 /title
     *      / cst36 /exname(h8)/ cst63 /delv(l2)/ debug /jtest,jpot
     *      / 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)/ cst209 /zname
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst14 /ipoly
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase
     *      / cst207 /ctrans(k0,k5),ictr(k5),itrans
     *      / cst208 /ifct,idfl/ csta9 /tcname(k5)/ cst102 /ivfl
     *      / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx     

      save blank5, blank8, blnk10
 
      data blank5,blank8,blnk10/'     ','        ','          '/
c-----------------------------------------------------------------------
c                             open n1
c                             the default designation for the input
c                             read from n1 is 'in1.dat', if this file
c                             does not exist terminate execution:
1     write (*,1050)
      read (*,1200) n2name
      open (n1, file = n2name, iostat = ierr, status = 'old')

      if (ierr.ne.0) then
c                                 system could not find the file
         write (*,1040) n2name
         write (*,*) ' Try again (y/n)?'
         read (*,1020) y
         if (y.ne.'Y'.and.y.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                             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
c                             read new component definitions:
      do i = 1, itrans
         read (n1,1120) tcname(i), ictr(i)
         read (n1,*) (ctrans(j,i), j = 1, icmpn)
      end do
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                  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                             0d0.
c     read (n1,*,err=998) ps
      ps = 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, k5
         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.k5) call error (197,r,k5,'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                             (no fluidc omponents) and iuyn = 1 otherwise.
      do 350 i = 1, 2
         vname(3+i) = blank8
         if (rname(i).eq.blank5) goto 350
         jmct = jmct+1
         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)
350   continue

      icomp = icomp+jmct
      iprct = icomp-icp
      jprct = icomp-jmct
c                             read bulk composition stuff
      read (n1,*) jbulk
c                             get the bulk composition:
      if (jbulk.gt.0) read (n1,*) (cblk(j), j = 1, jbulk)
c                             polygon file if icopt = 5:
      if (icopt.eq.5.or.(icopt.eq.1.and.jbulk.ne.0).and.io4.ne.1) then
c                             open polygon list on unit n7:
         if (icopt.eq.5) then 
 
            open (n7, file = 'c'//n4name, iostat = ierr, status = 'old')
            if (ierr.ne.0) then
c                                 system could not find the file
               write (*,1040) 'c'//n4name
               if (ierr.ne.0) call error (121,r,n7,'c'//n4name)
c                                 try again
            end if
         end if 
c                                 now create a special plot output file
c                                 name = b + plotfile name
         open (n5, file = 'b'//n4name)
         write (*,1220) 'b'//n4name
      else if (icopt.eq.5.or.(icopt.eq.1.and.jbulk.ne.0)) then
         write (*,1220) 'none requested'
      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.and.icopt.ne.5) 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.0d0) 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')
      end if
c                             now check iv(2):
      if (icopt.eq.1) 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.0d0) 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 i = 1, 3
         if (iv(i).gt.3) kct = kct + 1
      end do 
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 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 if (jmct.ne.0) then
            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
      end do 
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 solution model 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                             override input if jbulk ne 0:
      if (jbulk.ne.0) isudo = 0    
c                             read title for the calculation:
      read (n1,1160) title
c                             compute interval range of variables
      do i = 1, 5
         delv(i) = vmax(i) - vmin(i)
      end do  
c                             number of phases in invariant equilibria
      icp2 = icp + 2
c                              set convergence criteria for routine
c                              univeq:
      call concrt
90    goto 99

998   call error (27,r,i,n2name)

1000  format (10(a10,1x))
1010  format (a2,a5,a1)
1020  format (a1)
1030  format (a2,a5,a6)
1040  format (' **ERROR VER191** 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 (12(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,/)
1220  format (/,' Writing bulk composition plot output to file: ',a14)
1230  format (' Enter the polygon file name (i.e. the file created',/,
     *        ' with POLYGN), left justified: ')
1240  format ('b',a13)

99    end

      subroutine input2
c----------------------------------------------------------------------
c                     input2
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k4=22,k5=12,h5=7,h6=500,
     *           k0=25,k10=200,h8=100)
 
      character*5 cname, cmpnt, dname*40, zname, name*8, exname*8

      double precision twt(k5)
 
      integer icout(k0)
 
      common/ cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst88 /vu(2,k1),jprct,jmct/ cst208 /ifct,idfl
     *      / cst3  /x(k1,k5)/ cst20 /comps(k1,k5)/ cst36 /exname(h8)
     *      / cst209 /zname/ cst6 /icomp,istct,iphct,icp
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst5  /p,t,xco2,u1,u2,tr,pr,r,ps/ cst113 /ifproj
     *      / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(k5)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase
     *      / cst60 /ipoint,imyn/ csta6 /name/ cst18 /jvol(2),ivol
     *      / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind
c-----------------------------------------------------------------------
c                               initialization for each data set
 
c                               for k1 phases:
      do i = 1, k1
         do j = 1, 2
            vu(j,i) = 0d0
            vf(j,i) = 0d0
         end do 
         do j = 1, h5
            vs(j,i) = 0d0
            ikp(i) = 0
         end do
      end do 
c                               other counters and flags:
      do i = 1, h5
         isct(i) = 0
      end do 
c                               counters for bounds
      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 i = 1, icmpn
         icout(i) = 0
      end do
      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
               jvol(idfl) = i
            else if (j.eq.idco2) then
               iff(2) = i
               idfl = idfl + 1
               jvol(idfl) = i
            end if 
            im = 1
30       continue
         ivol = idfl
c                               write error message if a component
c                               was not found:
         if (im.ne.0) goto 20
         write (n3,1230) cname(i), (cmpnt(k), k = 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 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
         end do  
      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 j = 1, icmpn
         if (icout(j).eq.0.and.comp(j).ne.0d0) goto 50
      end do 
      do j = 1, icmpn
         if (comp(j).lt.0d0) then
            call warn (13,r,j,name)
            goto 50
         end if
      end do 
c                               reject the phase if it contains
c                               a thermodynamic component:
      do j = 1, icp
         if (comp(ic(j)).ne.0d0) goto 50
      end do
c                               now load the phase if it has
c                               the saturated component idc:
      do j = isat, 1, -1
         idc = icp + j
         if (comp(ic(idc)).ne.0d0) 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
      end do 
 
      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    do i = 1, isat
         if (isct(i).lt.1) call error (15,r,i,cname(icp+i))
      end do
c                               read data for the remaining
c                               phases of appropriate composition.
      istct = iphct + 1
c
      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 = 0d0
      do 310 i = 1, icp
         if (comp(ic(i)).lt.0d0) 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))
      if (tot.eq.0d0) then
c                               reject phase with no composition:
c                               this makes the test in sort redundant.
         call warn (27,r,i,name)
         iphct = iphct - 1
         goto 220
      end if
c                               store thermodynamic parameters:
      call loadit (iphct)
      do 320 i = 1, icp
320      x(iphct,i) = comp(ic(i)) / tot
c                               next phase:
      goto 220
1230  format (' **error ver013** ',a5,' is an incorrect component'
     *       ,' name, valid names are:',/,12(1x,a5))
1240  format (' check for upper/lower case matches or extra blanks',/)

999   end

      subroutine outchm
c-------------------------------------------------------------------
c outchm writes new chemographies to unit n3 as generated.
c flag ichem (1) determines if outchm is called, no output will be 
c generated if flag io3 = 1 (see introductory comments).
c-------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k5=12,l2=5,j9=40000,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
     *      / cst5  /v(l2),tr,pr,r,ps/ cst6 /icomp,istct,iphct,icp
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst79 /isoct   
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst8  /names(k1)/ cst96 /idcf(k5,j9),icfct
     *      / cst60 /ipoint,imyn/ cst27 /jasmbl(j9),iasmbl(j9)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
c-----------------------------------------------------------------------
c
      goto (999),io3
c                             output independent potential variable
c                             constraints:
      write (n3,1160)
      write (n3,1000)
      write (n3,1010) (vname(jv(i)), v(jv(i)), i = 1, ipot)
      write (n3,1015)
c                             next the stable assemblages:
      goto (100,200,300,400),icp
c                             output high order chmographies
      do i = 1, icfct
         write (n3,7020) iasmbl(i), (names(idcf(j,i)), j = 1, icp)
      end do 
      goto 800
c                             output unary chemographies.
100   write (n3,1020) names(idcf(1,1))
      goto 800
c                             output binary chemographies.
200   write (n3,2020) ((names(idcf(j,i)), j = 1, icp), iasmbl(i), 
     *                 i = 1, icfct)
      goto 800
c                             output ternary chemographies:
300   write (n3,3020) ((names(idcf(j,i)), j = 1, icp), iasmbl(i), 
     *                i = 1, icfct)
      goto 800
c                             output quaternary chemographies:
400   write (n3,4020) ((names(idcf(j,i)), j = 1, icp), iasmbl(i), 
     *                i = 1, icfct)
      goto 800
c                             output phases consistent with component
c                             saturation constraints:
800   goto (9000),isyn
      write (n3,6000)
      write (n3,2030) (names(idss(i)), i = 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 (variance flag in parentheses):',/)
1020  format (26x,a8)
1160  format (/,80('-'),/)
2020  format (3(1x,a8,'-',a8,'(',i1,')',3x))
2030  format (6(2x,a8))
3020  format (2(1x,2(a8,'-'),a8,'(',i1,')',2x))
4020  format (2(1x,3(a8,'-'),a8,'(',i1,')'))
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 **',/)
6030  format (' ** immiscibility occurs in one or more of the ',
     *        'stable solution phases **',/)
7020  format (1x,'(',i1,')',12(1x,a8))
 
999   end
 
      subroutine outhed
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k5=12,l2=5,j9=40000,h5=7,h6=500,
     *           h9=30)
 
      character*8 names,vname,gname*78,xname*18,title*72,fname*10
 
      common/ cst3   /x(k1,k5)/ cst8 /names(k1)
     *      / csta2  /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst40  /ids(h5,h6),isct(h5),icp1,isat,io2/ csta8 /title
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst6   /icomp,istct,iphct,icp/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst60  /ipoint,imyn/ cst96 /idcf(k5,j9),icfct
     *      / cst61  /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst103 /isec,icopt,ifull,imsg,io3p/ cst79 /isoct
     *      / cst27  /jasmbl(j9),iasmbl(j9)/ csta7 /fname(h9)
c-----------------------------------------------------------------------
c                             number of components, phase counters,
c                             assemblage counter, fluid saturation flag,
c                             component saturation flag, ipot = number
c                             of independent potential variables:
      write (n4,1020) icopt
c
      write (n4,1020) icp,istct,iphct,ipoint,ifyn,isyn,ipot,isoct
c                             write graphics code variable names:
      write (n4,1030) (ig(jv(i)),gname(jv(i)), i = 1, ipot)
c                             write a blank record as left caption
      write (n4,1050) title
c                             write phase names
      write (n4,1015) (names(i), i = 1, iphct)
c                             write phase coordinates
      write (n4,1025) ((x(i,j), j = 2, icp), i = istct, iphct)
c                             write solution phase flags:
      write (n4,1020) (ikp(i), i = 1, iphct)
c                             solution names
      if (isoct.ne.0) write (n4,1080) (fname(i),i = 1, isoct)
c                             write graphics code names for x variables:
      write (n4,1010) (xname(i), i = 1, icp)
c
1010  format (5a18)
1015  format (10a8)
1020  format (20(i5,1x))
1025  format (11(f6.4,1x))
1030  format (i2,a78)
1050  format (a72)
1080  format (8a10)
 
99    end

      subroutine outgrf
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k5=12,l2=5,j9=40000,h5=7,h6=500)
 
      common/ cst83  /ig(l2)/ cst40  /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn
     *      / 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
     *      / cst96  /idcf(k5,j9),icfct/ cst27  /jasmbl(j9),iasmbl(j9)
     *      / cst24 /ipot,jv(l2),iv(l2)
c-----------------------------------------------------------------------
c                             stable configurations, phases are
c                             labelled by the index 'i' in the
c                             list of phases.
c                             write values of the independent potentials
      write (n4,1040) (v(jv(i)), i = 1, ipot)

      if (icp.ne.2) then
         write (n4,*) icfct
      else
         write (n4,*) icfct+1
      end if 

      goto (99,950,960),icp
c                             binary is a special case (1-d)
950   write (n4,1020) (idcf(1,j), j = 1, icfct),idcf(2,icfct)
      goto 995
c                             higher order:
960   write (n4,1020) ((idcf(j,i), j = 1, icp), i = 1, icfct)
      goto 995
c                             write assemblage flags
995   if (icp.gt.2) write (n4,1020) (iasmbl(j), j = 1, icfct)
      goto (99),isyn
      write (n4,1020) isat
      write (n4,1020) (idss(i), i = 1, isat)

1020  format (20(i5,1x))
1040  format (5(g12.6,1x))

99    end
 
      subroutine prtpot
c---------------------------------------------------------------------
c prtpot outputs the dependent potentials for phase assemblages
c determined by the combin subroutine. it is called if the flag
c jpot, in common debug, is not equal 1.
c---------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12) 

      character*8 names

      common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst23 /a(k5*k5),u(k5+1),iar(2*k5+4)/ cst8 /names(k1)
     *      / cst6  /icomp,istct,iphct,icp/ cst3 /x(k1,k5)
     *      / cst52 /h,id(k5+1),ic(k5)/ debug /jtest,jpot

      if (jpot.ne.0) goto 100

      goto (10,20,30,40,50),icp

      write (n3,1060) (names(id(i)),i=1,icp)
      write (n3,1070) (u(i),i=1,icp)
      goto 99

10    write (n3,1010) names(id(1)),u(1)
      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 

100   do 80 i = 1, icp
            g = 0d0
            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) 

1010  format (1x,a8,1x,g14.7)          
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 (10(1x,a8))
1070  format (10(1x,g14.7))
1080  format (1x,a8,8(1x,g14.7))

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-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,l2=5,h5=7,j9=40000,h6=500,
     *           k7=k5+1)

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

      double precision value, vt 

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ csta4 /cname(k5)
     *      / cst6 /icomp,istct,iphct,icp/ cst96 /idcf(k5,j9),icfct
     *      / cst5  /v(l2),tr,pr,r,ps/ cst7 /iflag/ cst38 /iedge
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)/ cst39 /ichem
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst102 /ivfl
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst80 /irchk(k2)
     *      / cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst31 /vn(k2,k7),irct,ird
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst65 /vt(j9),vti(j9),jok(j9),kok(j9),index
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst49 /value,ismax,igot      
c-----------------------------------------------------------------------
c                                  initialization:
      irend = 0
      irct = 0
 
      do 70 i = 1, k2
70       irchk(i) = 0
c                                  set starting values for search
      v(iv1) = vmin(iv1)
      v(iv2) = vmin(iv2)
c                                  determine the stable chemography.
      call gall 
      call combin

      if (icopt.eq.3) then 
         iedge= 1
         if (io4.eq.0) then 
c                             write chemographies for mixed variable
c                             diagrams:
            call outgrf
c                             write min and max values for the
c                             independent potentials:
            write (n4,*) vmin(iv1),vmax(iv1)
         end if 
      end if
c                                  output chemography if requested.
c                                  this option is archaic with newhld.
      goto (10),ichem
      call outchm  
c                                  write tier header to print file:
10    goto (15),io3
      write (n3,1160)
      write (n3,1170) vname(iv1),vname(iv2)
c                                  write potential variable sectioning
c                                  constraints:
      if (ipot.gt.2) then
         write (n3,1180)
         write (n3,1190) (vname(jv(i)),v(jv(i)), i = 3, ipot)
      end if
c                                  write saturated and buffered
c                                  component names:
      goto (100),isyn
      write (n3,1200) (cname(i+icp), i= 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 i = 1, icfct
         vt(i)  = vmin(iv1)
         vti(i) = vmin(iv2)
         kok(i) = 1
         jok(i) = 1
      end do 
      
      write (*,1000) icfct
c                                  look for reactions involving each
c                                  assemblage.
      knct = 1
40    jnct = icfct
c
      do 50 i = knct, jnct
         if (imsg.eq.0) write (*,1220) i,icfct-i
c                                  load id's into array idv:
         do j = 1, icp
            idv(j) = idcf(j,i)
         end do 
c                                  get the lower/upper decomposition
c                                  of the transpose of the concentration
c                                  matrix:
         call pivots (ier)
         if (ier.eq.1) then
c                                  error, a degenerate assemblage in
c                                  pivots due to a programming error.
            call warn (29,v(1),ier,'NEWHLD')
            goto 50
         end if 
c                                  set flag for maxend:
         ismax = jok(i) 
         igot = 0
         value = vmin(iv1)
c                                  test stability along an edge of the
c                                  diagrams coordinate frame:
30       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 (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)
c                                  after tracing the net attempt
c                                  to resume testing the stability
c                                  of the assemblage:
         if (igot.eq.1) then
            if (ismax.eq.1) then
c                                  resume search of travserse 1:
               vt(i) = value + 1d-2*dv(iv1)
               vti(i) = vmin(iv2)
            else if (ismax.eq.2) then
c                                  resume search of travserse 2:
               jok(i) = 2
               vt(i) = value + 1d-2*dv(iv2)
               vti(i) = vmax(iv1)
            else if (ismax.eq.3) then
c                                  resume search of travserse 3:
               jok(i) = 3
               vt(i) = value - 1d-2*dv(iv1)
               vti(i) = vmax(iv2)
            else
c                                  resume search of travserse 3:
               jok(i) = 4
               vt(i) = value - 1d-2*dv(iv2)
               vti(i) = vmin(iv1)
            end if
c                                  test for stability before search? 
            do j = 1, icp
               idv(j) = idcf(j,i)
            end do 
            call pivots (ier)
            igot = igot + 1
            goto 30 
         end if 
50    continue
      if (icfct.eq.jnct) goto 90
      knct = jnct + 1
      goto 40
c                                  return for icopt=1
90    if (icfct.eq.j9) call warn (205,r,j9,'NEWHLD')
      goto (99),icopt
c                                  for icopt=3 sort and output
c                                  reactions
      call onedim
      
1000  format (' Initial number of',
     *        ' divariant assemblages to be tested is:',i4)
1160  format (/,80('-'),/)
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 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.',//)
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=5000,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.0d0) 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
      subroutine gall 
c-----------------------------------------------------------------------
c subroutine gall computes molar free energies of all phases.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k4=22,h5=7,h6=500,
     *           k12=14,k13=40000,k10=200,h9=30)
      double precision g0(2)

      common/ cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst2  /g(k1)/ cst6  /icomp,istct,iphct,icp
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1)
     *      / cst60 /ipoint,imyn/ cst79 /isoct
c-----------------------------------------------------------------------
c                                 compute the chemical potential
c                                 of the projected components.
      call uproj
c                                 first do the endmembers:
      do id = istct, ipoint
         g(id) = gcpd(id) - uf(1) * vf(1,id) - uf(2) * vf(2,id)
         do j = 1, isat
            g(id) = g(id) - us(j) * vs(j,id)
         end do
      end do 
c                                 now do solutions:
      do 30 i = 1, isoct
c                                 check if normal solution:
         if (ifp(id).eq.0) then 
c                                 ok it's normal:
            do 40 j = 1, jend(i,2)
               g(id) = 0d0
               do 50 k = 1, jend(i,1) 
50                g(id) = g(id) + g(jend(i,2+k)) * sxs(ixp(id)+k)
                  g(id) = g(id) + gexces(id)
40                id = id + 1
         else 
c                                 it's a fluid compound, the way 
c                                 things are now it must have two
c                                 components.
c                                 first get the reference state 
c                                 free energy:
            do 60 j = 1, 2
60             g0(j) = gzero(jend(i,2+j))

            do 70 j = 1, jend(i,2)
               g(id) = g0(1) * sxs(ixp(id)+1) + g0(2) * sxs(ixp(id)+2)
     *               + gexces(id)
70               id = id + 1
         end if 
30    continue

      end
 
      function gproj (id)
c-----------------------------------------------------------------------
c gproj computes the projected molar free energy of the
c phase with index id.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k4=22,h5=7,h6=500,k10=200)

      common/ cst1 /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
c-----------------------------------------------------------------------
      graw = gphase(id)-uf(1)*vf(1,id)-uf(2)*vf(2,id)

 

      do j = 1, isat

         graw = graw - vs(j,id) * us(j)

      end do 

      gproj = graw

      end

      function grxn(i) 
c-----------------------------------------------------------------------
c grxn computes the free energy of univariant equilibria
c defined by the data in commonn block cst21 which is initialized
c in the subprogram balanc.  grxn is partially redundant with
c the function gphase but because of the frequency that these
c these routines are used a significant increase in efficiency is
c gained by maintaining separate functions.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k12=14,k13=40000,h5=7,h9=30,k5=12,k7=k5+1)

      common/ cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1)
     *      / cst60 /ipoint,imyn/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst11 /fh2o,fco2
     *      / cst201 /vuf(2),vus(h5),iffr,isr
     *      / cst6  /icomp,istct,iphct,icp/ cst25 /vnu(k7),idr(k7),ivct
c-----------------------------------------------------------------------
c                                 compute potentials of saturated phases
c                                 and components, note that in this
c                                 version of vertex the stoichiometry of
c                                 such components may vary.

c                                 no saturated phase components and no
c                                 saturated components:
      if (iffr.eq.1) goto (10),isyn
c                                 note that this call to uproj makes a
c                                 subsequent call in gall redundant if
c                                 sfol1 is used to trace a univariant
c                                 curve.
      call uproj
c                                 compute free energy change of the rxn
10    graw= 0d0

      do j = 1, ivct
         graw = graw + vnu(j) * gproj(idr(j))
      end do 

      grxn = graw

      end

      subroutine soload (z,isoct,icpct,ixct,tname)
c--------------------------------------------------------------------------
c soload - loads/requires solution properties: 

c   jend(h9,k12)  - h9 is the maximum number of solutions
c                   k12 is the maximum number of endmembers pers
c                   solution plus two.
c   jend(i,1)     - is the number of endmembers in solution i.
c   jend(i,2)     - is the number of pseudocompounds of solution i.
c   jend(i,3-3+j) - are the indices of the j endmembers in solution i.
c   sxs(k13)      - contains the mole fractions of the endmembers
c                   in the pseudocompounds.
c   ixp(i)        - a pointer that locates the first mole fraction (- 1)
c                   of the ith pseudocompound, the remaining mole fractions
c                   follow sequentially (as in jend(i,3-3+j)).
c   ikp(i)        - the index of the solution corresponding to pseudocompound i.
c   ifp(i)        - flag, 1 if pseudocompound i is a fluid, else 0.
c   exces(j,i)    - the excess function of pseudocompound i, accounts for 
c                   excess properties and configurational entropy as a function
c                   of pressure and temperature:

c                       gexces(i) = exces(1) + exces(2)*T + exces(3)*P
c--------------------------------------------------------------------------

      implicit double precision (a-g,o-z),integer (h-n)
  
      parameter (k1=5000,k2=3000,k4=22,k5=12,m0=5,h5=7,h6=500,
     *           k12=14,k13=40000,k10=200,h9=30)

      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
     *      / 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
     *      / 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/ cst20 /comps(k1,k5)/ cst88 /vu(2,k1),jprct,jmct
     *      / cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst6  /icomp,istct,iphct,icp/ cst8 /names(k1)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst5  /p,t,xco2,u1,u2,tr,pr,r,ps/ cst12 /cp(k5,k1)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst3  /x(k1,k5)
     *      / cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1)
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn
c                              eliminate end-member compositions
      do i = 1, isp(1)
         call zchk (z(1,i),tname)
         do j = 1, isp(2)
            call zchk (z(2,j),tname)
            do k = 1, isp(3)
               call zchk (z(3,k),tname)
               if (z(1,i)*z(2,j)*z(3,k).gt.0.999d0) goto 999
            end do
         end do 
      end do     
c                               eliminate out of range compositions:  
c                               is this necessary?
      do 50 i = 1, isite
         if (imd(i).gt.0) goto 50
         do j = 1, isp(i) - 1
            if (z(i,j).gt.xmx(i,j).or.z(i,j).lt.xmn(i,j)) goto 999
         end do   
50    continue
c                               the composition is acceptable.
      iphct = iphct + 1
      icpct = icpct + 1

      if (iphct.gt.k1) call error (41,z(1,1),k1,'SOLOAD')

      ikp(iphct) = isoct
      ixp(iphct) = ixct
c                               check if the solution is a fluid:
      if (thermo(3,imsol(1,1,1)).eq.0d0) then
         ifp(iphct) = 1
c                               don't allow gonzoids to treat fluid
c                               as a solution if the user has also 
c                               specified fluid saturation
         if (ifyn.eq.0) call error (43,r,i,tname)       
      else
         ifp(iphct) = 0
      end if 
c                               encode a name
      if (isite.eq.1) then
         if (isp(1).eq.2) then
           if (z(1,1).eq.1d0) then 
              write (names(iphct),1050) names(imsol(1,1,1)),100
           else if (z(1,1).ge.0.04d0.and.z(1,1).le.0.96d0) then
              write (names(iphct),1030) names(imsol(1,1,1)),
     *                                  idint(1d2*z(1,1))
           else   
              write (names(iphct),1070) names(imsol(1,1,1)),
     *                                  1d2*z(1,1)
           end if 

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

         end if 

      else if (isite.eq.2.and.isp(1).eq.2.and.isp(2).eq.2) then
 
          int1 = idint (1d2*z(1,1))
          int2 = idint (1d2*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(1d1*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 = 0d0
      esum = 0d0
      atw(iphct) = 0d0

      do i = 1, icp
         cp(i,iphct) = 0d0
      end do 
      do i = 1, icomp
         comps(iphct,i) = 0d0
      end do 
c                              load constants:
      ctotal = 0d0
c
      do i = 1, isp(1)
         do j = 1, isp(2)
            do 180 k = 1, isp(3)

               zp = z(1,i)*z(2,j)*z(3,k)
                     ixct = ixct + 1
               if (ixct.gt.k13) call error (40,zp,k13,'SOLOAD')
               sxs(ixct) = zp
               if (zp.eq.0d0) goto 180

               id = imsol(i,j,k)
               esum = esum + zp * ecoef(i,j,k)
               atw(iphct) = atw(iphct) + zp * atw(id)

               do l = 1, icp
                  cp(l,iphct) = cp(l,iphct) + zp * cp(l,id)
                  ctotal = ctotal + zp * cp(l,id)
               end do 
               do l = 1, icomp
                  comps(iphct,l) = comps(iphct,l) + zp * comps(id,l)
               end do 

180         continue
         end do 
      end do  
c                              check if the phase consists
c                              entirely of saturated components:
      if (ctotal.eq.0d0) then
         call satsrt 
      else
         do i = 1, icp
            x(iphct,i) = cp(i,iphct) / ctotal
         end do
      end if 
c                              compute ideal configurational negentropy:
      if (nsite.eq.0.and.ist(1).ne.0) then  
         do i = 1, isite 
            do j = 1, isp(i)
               if (z(i,j).ne.0d0) smix = smix + 
     *                            float(ist(i))*z(i,j)*dlog(z(i,j)) 
            end do 
         end do     
      else if (nsite.ne.0) then
          call omega (z,dlnw,tname)
          smix = -dlnw + esum
      end if 
c                              save it:
      exces(2,iphct) = exces(2,iphct) + r * smix
c                               load stoichiometric coefficients of
c                               saturated components:
      do k = 1, isat
         vs(k,iphct) = comps(iphct,icp+k)
      end do
c                               load stoichiometric coefficients of
c                               mobile components:

      do i = 1, jmct
         vu(i,iphct) = comps(iphct,jprct+i)
      end do 
 
      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 i = 1, iterm
            zp = 1d0
            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 j = 1, 3
                 exces(j,iphct) = exces(j,iphct) + zp * wg(i,j)
            end do 
         end do  
      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 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=5000,k2=3000,k4=22,h5=7,h9=30,
     *           m0=5,k12=14,k13=40000,k10=200)
 
      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 /ipoint,imyn/ csta7 /fname(h9)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst79 /isoct/ cst6  /icomp,istct,iphct,icp
     *      / cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1)

      save uname

      data uname/' ',' '/
c-----------------------------------------------------------------------
c                               set compound counter ipoint.
      ipoint = iphct
      ixct = 0 
c                               initialize match flag:
      im = 0
c                               no request for solutions, goto 999:
      goto (999),io9
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 i = isite + 1, 3
            z(i,1) = 1d0
            z(i,2) = 0d0
            z(i,3) = 0d0 
         end do 
      end if          
c                               subdiv generates ntot compositions,
c                               generate the compound data for
c                               each solution:
c                               save the identities of the endmembers
      jend(im,1) = isp(1)*isp(2)*isp(3)
      if (jend(im,1)+2.gt.k12) call error (39,zt,k12,'INPUT9')
      icpct = 0
      itic = 2

      do i = 1, isp(1)
         do j = 1, isp(2)
            do k = 1, isp(3)
               itic = itic + 1
c                               if the endmember has not been
c                               "killed" by a solution model
c                               set the solution identifier flag.
               if (ikp(imsol(i,j,k)).ne.-1) ikp(imsol(i,j,k)) = im 
               jend(im,itic) = imsol(i,j,k)
            end do 
         end do 
      end do 

      do h = 1, ntot
c                               load the composition into
c                               a the site fraction array:
         do i = 1, isite
            zt = 0d0
            do j = 1, isp(i) - 1
               z(i,j) = y(h,i,j)
               zt = zt + z(i,j)
            end do 
            z(i,isp(i)) = 1d0 - zt
         end do 
c                               generate the pseudocompound:
         call soload (z,im,icpct,ixct,tname)
      end do 

      jend(im,2) = icpct

      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 i = 1, iphct
c                             reset ikp
         if (ikp(i).lt.0) ikp(i) = 0
      end do 

      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 

         do i = 1, im
            fname(i) = sname(i)
         end do 
      end if 

      isoct = im

1000  format (/,' the following solution models will be considered:',/)
1010  format (7(2x,a10))
1020  format (/,' Of the requested solution models:',/)
1040  format (/,' no models will be considered.',/)

      end 

      subroutine abload (*)
c---------------------------------------------------------------------
c abload assembles the matrix 'a' and vector 'b' for icp component
c sytems and then solves the equation ax = b, the vector x is
c returned in 'b'.
 
c   referenced by: combin
c   references to: subst,factor
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k5=12)
 
      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),iophi,idphi,
     *               iiphi,iflg1
     *      / cst2  /g(k1)/ cst12 /cp(k5,k1)/ cst52 /hcp,id(k5+1),ic(k5)
c-----------------------------------------------------------------------
      do i = 1, hcp
         do j = 1, hcp
            a(i,j) = cp(j,id(i))
         end do
      end do

      call factor (a,hcp,ipvt,ier)
 
      goto (99),ier
 
      do i = 1,  hcp
         b(i) = g(id(i))
      end do 
 
      call subst (a,ipvt,hcp,b)
 
      return
99    return 1
      end

      subroutine asschk
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k5=12,l2=5)
 
      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),iophi,idphi,
     *               iiphi,iflg1/ cst87 /delt(l2),dtol,utol,ptol
     *      / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dgr
     *      / cst2 /g(k1)/ cst12 /cp(k5,k1)/ cst7 /iflag
c-----------------------------------------------------------------------
c                                 determine chemical potentials
      iflag = 0
      do i = 1, icp
         b(i) = g(idv(i))
      end do
c
      call subst (a,ipvt,icp,b)
c                                 test phases against the
c                                 assemblage idv
      do 20 i = istct, iphct
         gphi = 0d0
         do j = 1, icp
            gphi = gphi + cp(j,i) * b(j)
         end do 
         dg = g(i) - gphi
         if (dg.gt.dtol) goto 20
c                                check that a phase is not metastable
c                                with respect to itself, this
c                                could be remedied by changing the 
c                                value of dtol.
         do j = 1, icp
            if (idv(j).eq.i) goto 20
         end do 
c                                assemblage is metastable with respect
c                                to phase idphi.
         iflag = iflag + 1
         idphi = i
         if (dabs(dg).gt.1d2*utol) iflag = iflag + 1
         goto (20,99,99),iflag
20    continue

99    end
 
      block data
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (h5=7,k1=5000,k4=22,k6=1000,k10=200)
 
      common/ cst5  /p,t,xco2,u1,u2,tr,pr,r,ps/ cst85 /pp,tt,yy,rr
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst32 /ptx(k6),ipt2
     *      / cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
c-----------------------------------------------------------------------
      data iff/2*0/,ipt2/0/
c
      data us, uf/ h5*0d0, 2*0d0/
c
      data r,rr/8.3144126d0,83.14d0/
c                              the following data statement assigns
c                              logical unit numbers for i/o.
      data n1,n2,n3,n4,n5,n6,n7,n8,n9/21,22,23,24,25,26,27,28,29/

      end

      subroutine assir 
c----------------------------------------------------
c irct - counter of univariant rxns.
c ivarrx(irct) - variance flag of the iirct(th) rxn (set in balanc as
c                 ivar in cst61.
c irchk(irct) - a flag indicating if the iirct(th) rxn occurs on the
c                edge of a diagram (initialized in newhld and set in
c                coface or sfol1).
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,l2=5)

      common/ cst5  /v(l2),tr,pr,r,ps/ cst80 /irchk(k2)
     *      / cst31 /vn(k2,k7),irct,ird/ cst25 /vnu(k7),idr(k7),ivct
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst35 /irv(k2)/ csta1 /irnms(k2,k7)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
c-----------------------------------------------------------------------
c                         check if the equilibrium is already known
      do 70 i = 1, irct
         if (ivct.ne.irv(i)) goto 70
         do 80 j = 1, ivct
            do k = 1, ivct
               if (irnms(i,k).eq.idr(j)) goto 80
            end do 
c                         no match with irnms(i,k) try i + 1.
         goto 70
80       continue
c                         match found with irnms(i,k)
         ird = i
         goto 99
70    continue
c                         the reaction is new:
      irct = irct + 1

      ird = irct
      if (irct.gt.k2) call error (181,r,k2,'ASSIR')
      irv(irct) = ivct
      ivarrx(irct) = ivar

      do i = 1, ivct
         vn(irct,i) = vnu(i)
         irnms(irct,i) = idr(i)
      end do 

99    end

      subroutine assri (ier)
c----------------------------------------------------------------------
c irct         - reaction counter
c irv(irct)    - number of phases in the irct(th) rxn
c                (set in balanc as ivct cst25).
c ivarrx(irct) - variance flag of the iirct(th) rxn (set in balanc as
c                 ivar in cst61.
c irchk(irct)  - a flag indicating if the iirct(th) rxn occurs on the
c                edge of a diagram (initialized in newhld and set in
c                coface or sfol1).
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,k8=k5+2,l2=5)

      common/ cst5  /v(l2),tr,pr,r,ps/ cst80 /irchk(k2)
     *      / cst31 /vn(k2,k7),irct,ird
     *      / cst25 /vnu(k7),idr(k7),ivct
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5      
     *      / cst35 /irv(k2)/ csta1 /irnms(k2,k7)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)
c-----------------------------------------------------------------------
c                         check if the reaction is known
      do 70 i = 1, irct
         if (ivct.ne.irv(i)) goto 70
         do 80 j = 1, ivct
            do 90 k = 1, ivct
90             if (irnms(i,k).eq.idr(j)) goto 80
c                         no match with irnms(i,k) try i + 1.
            goto 70
80       continue
c                         match found with irnms(i,k), now check
c                         that the conditions are different, if
c                         ier = 1, univeq failed, skip the check
         goto (75), ier

         if (dabs(vip(iv1,i)-v(iv1)).gt.0.05*dv(iv1)) goto 40
         
75       ird = i
         ier = 1
         goto 99
70    continue
c                         the reaction is new:
40    irct = irct + 1
      ier = 0
      if (irct.gt.k2) call error (181,r,k2,'ASSRI')
      ird = irct
      irv(ird) = ivct
      ivarrx(ird) = ivar

      do i = 1, ivct
         vn(ird,i) = vnu(i)
         irnms(ird,i) = idr(i)
      end do 

99    end

      subroutine balanc (b,idv,idphi,ier)
c-----------------------------------------------------------------------
c balanc balances reactions between the phases in idv and idphi
c if the system is saturated with any phases these are taken
c into account.  once the reaction is determined balanc determines
c the thermodynamic parameters needed for subroutine univeq and
c saves them in the array 'rxn'.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k4=22,k5=12,k7=k5+1,
     *           h5=7,h6=500,k10=200)

      dimension a(k5,k5), ipvt(k5), idv(k5), b(k5+1)

      common/ cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst6  /icomp,istct,iphct,icp/ cst25 /vnu(k7),idr(k7),ivct
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst20 /comps(k1,k5)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst88 /vu(2,k1),jprct,jmct/ cst12 /cp(k5,k1)
     *      / cst201 /vuf(2),vus(h5),iffr,isr
c-----------------------------------------------------------------------
c                                 assemble transpose of the composition
c                                 matrix for an assemblage of icp phases
      ivar =1
                                                
      do i = 1, icp
         do j = 1, icp
            a(j,i) = cp(j,idv(i))
         end do 
         idr(i) =idv(i)
      end do 
c                                 assemble the composition vectore
c                                 for the icp+1 th phase tangent to
c                                 the free energy plane.
      ip=icp+1
      do i = 1, icp
         b(i) =comps(idphi,i)
      end do 
      idr(ip) =idphi
      b(ip) = -1d0
c                                  balance the equilibrium
      call factor (a,icp,ipvt,ier)
      goto (999),ier
      call subst (a,ipvt,icp,b)
c                                  eliminate phases with vnu= 0
      ivct= 0
      do 80 i = 1, ip
         if (dabs(b(i)).lt.1d-05) goto 80
         ivct=ivct + 1
         vnu(ivct) =b(i)
         idr(ivct) =idr(i)
80    continue
c                                 test to determine if the reaction
c                                 is pseudo-univariant, and set flag
c                                 ivar =2 if so, and =1 otherwise.
      iv1=ivct-1
      do 140 i = 1, iv1
         if (ikp(idr(i)).eq.0) goto 140
         i1=i + 1
         do 150 j =i1,ivct
            if (ikp(idr(i)).ne.ikp(idr(j))) goto 150
            ivar =2
c                                pseudo-univariant reject rxn if isudo=1
            goto (999),isudo
            goto 170
c                                note that this test will fail to
c                                distinguish between critical equilibria
c                                and pseudo-univariant equilibria.
150      continue
140   continue

170   goto (40),isyn
c                                determine stoichiometric coefficients
c                                of saturated components:
c                                (these aren't used for anything real)
      isr = 1
      do j = 1, isat
         vus(j) = 0d0
         do i = 1, ivct
            vus(j) = vus(j) + vnu(i) * comps(idr(i),icp+j)
         end do 
c                                flag isr is not used.
         if (vus(j).ne.0d0) isr = 0   
      end do 
c                                ok, now vus has the total deltas
   
c                                determine stoichiometric coefficients
c                                of saturated phase components:
      iffr = 1
40    do j = 1,2
         vuf(j) = 0d0
         do i = 1, ivct
            vuf(j) = vuf(j) + vnu(i) * vf(j,idr(i))
         end do 
c                                this is only the total of saturated
c                                phase components for the reaction
c                                (doesn't include saturated phases).
c                                flag iffr is used
         if (vuf(j).ne.0d0) iffr = 0
      end do 
c                                et finis
999   end

      subroutine combin
c----------------------------------------------------------------------------
c combinatorial algorithm for mfes. 
c----------------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,j9=40000)

      common/ cst6   /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)   
     *      / cst96  /idcf(k5,j9),icfct/ cst27 /jasmbl(j9),iasmbl(j9)
     *      / cst52  /hcp,id(k5+1),ic(k5)/ cst60 /ipoint,imyn
     *      / cst300 /cblk(k5),atw(k1),jbulk
     *      / cst103 /isec,icopt,ifull,imsg,io3p
c                                   initialization
      imyn = 1
      icfct = 0
      jcp = icp - 1

      iclose = 0
c                                   allow the user to specify the
c                                   initial assemblage: 
      if (jbulk.eq.icp.and.icopt.eq.0) then
         call tertex (icase)
         if (icase.eq.1) goto 99 
      end if
c                                   identify the seed stable assemblage
      call findas 
c                                   save the assemblage
      call assdc (icase)
c                                   finished if one component
      if (icp.eq.1.or.icase.eq.2.and.jbulk.ge.icp) goto 91
c                                   now use the c-1 phase simplicial
c                                   facets generated from the seed
c                                   assemblage to identify the remaining
c                                   assemblages.
      istart = 1
      jstart = 0
      jend = 0
      igen = 1

40    iend = icfct 
c      write (*,*) 'cycle ',igen,istart,iend,jstart,jend
      igen = igen + 1

      do 50 i = istart, iend
c                             each facet may generate icp-1 new psf's
         do 60 j = 1, jcp
c                             set vertices identities:
            itic = 0
c                             this loop generates the psf consisting of
c                             the icp-1 phases, saving the icp'th phase
            do k = 1, icp 
               if (k.eq.j) then 
                  idno = idcf(k,i) 
               else 
                  itic = itic + 1
                  id(itic) = idcf(k,i)
               end if 
            end do
c                              check that the assemblage is not degenerate
c                              a valid psf must span the icp dimensional 
c                              composition space, count the components:
            do k = 1, icp
               sum = 0d0
               do l = 1, jcp
c                             sum the amounts of component k in the psf 
                  sum = sum + cp(k,id(l))
               end do 
c                             if the component is 0, the psf is degenerate
               if (sum.lt.1d-5) goto 60
            end do
c                             next check if the psf matches any permutation
c                             of icp-1 phases in a facet tested earlier:
c                             this test is not essential, but it should
c                             save time. 
            do 80 k = 1, i-1
               do 70 l = 1, jcp
                  do m = 1, icp
                     if (idcf(m,k).eq.id(l)) goto 70
                  end do
c                             no match with assemblage k:
                  goto 80
70             continue
c                             matched assemblage k, reject the psf:
               goto 60
80          continue
c                             now identify the first potential vertice:

            do 10 k = istct, iphct
c                             don't allow phases of the original simplex.
               do l = 1, icp
                  if (k.eq.idcf(l,i)) goto 10
               end do 
c                             try to find a phase that is stable
c                             with respect to phase idno:
               id(icp) = k

               call abload (*10)

               if (dgphc(idno).gt.-1d-05) then
c                              the plane is stable
                  lastk = k + 1
                  goto 20   
               end if 
10          continue
c                               if no acceptable phase is found
c                               the assemblage must lie on an edge
c                               of the compositional space.
            if (iclose.eq.0) call warn (28,sign,1,'COMBIN')
            iclose = 1
            goto 60
c                               now do thermodynamic testing
20          if (lastk.gt.iphct) then
               call assdc (icase)
               if (icase.eq.2.and.jbulk.ge.icp) goto 91
               goto 60
            end if          

            do k = lastk, iphct
c                                   checkd replaces the trial vertice
c                                   with phase j, if the assemblage is 
c                                   found to be metastable with respect 
c                                   to phase j. given the speed of this
c                                   test, i doubt geometric testing would
c                                   be worthwhile here.
c                                   don't allow return of the vertice
c                                   idpsf(icp,i)
               if (k.ne.idno) call checkd (k)
            end do 
c                                   save the new assemblage:
            call assdc (icase)
            if (icase.eq.2.and.jbulk.ge.icp) goto 91
60       continue
50    continue 

      if (icfct.gt.iend) then
         if (jbulk.gt.0.and.jend.gt.0) then
c                                    eliminate blind assemblages if 
c                                    constrained bulk compositions
            itic = 0
            do i = jstart, jend
               if (jasmbl(i).ne.0) then
                  itic = itic + 1
                  do j = 1, icp
                     idcf(j,itic) = idcf(j,i)
                  end do
                  jasmbl(itic) = jasmbl(i)
               end if
            end do
c                                     now shift all other assemblages
            jtic =  jend-jstart+1-itic

            do i = jstart+itic, icfct-jtic

               k = i + jtic
               jasmbl(i) = jasmbl(k)
               do j = 1, icp
                  idcf(j,i) = idcf(j,k)
               end do
            end do
            jstart = istart - jtic
            jend   = iend - jtic
            istart = jend + 1
            icfct  = icfct - jtic  
            goto 40 
         else
            jstart = istart
            jend = iend
            istart = iend + 1
            goto 40
         end if 
      end if
c                                    check that something was found.
      if (icfct.eq.0) then
         if (iclose.eq.1.and.jbulk.ne.0) call error (42,sign,2,'COMBIN')
         call error (999,sign,2,'COMBIN')
      end if 
c                                    if constrained bulk composition
c                                    eliminate irrelevant assemblages
91    if (jbulk.ge.icp) then
         do j = 1, icp
            idcf(j,1) = idcf(j,icfct)
         end do
         icfct = 1
      else if (jbulk.ne.0) then
         itic = 0
         do i = 1, icfct
            if (jasmbl(i).eq.1) then
               itic = itic + 1
               do j = 1, icp
                  idcf(j,itic) = idcf(j,i)
               end do
            end if
         end do
      end if 

99    end 

      subroutine assdc (icase)
c-----------------------------------------------------------------------
c assdc counts (icfct) and assigns (idcf) c-component facets. the
c vertices are identified by the array id(icp).

c   icase = 0 => no new assemblage counted.
c   icase = 1 => new assemblage.
c   icase = 2 => new assemblage that bounds a specified bulk composition.

c   referenced by: findas
c   references to: none
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k5=12,j9=40000)
 
      common/ cst96  /idcf(k5,j9),icfct/ cst6 /icomp,istct,iphct,icp
     *      / cst27 /jasmbl(j9),iasmbl(j9)
     *      / cst12  /cp(k5,k1)/ cst300 /cblk(k5),atw(k1),jbulk
     *      / debug  /jtest,jpot/ cst52 /ht,id(k5+1),ic(k5) 

      save idid

      data idid/0/   
c-----------------------------------------------------------------------
      icase = 0
c                             test for equivalence with earlier
c                             assemblage:
      do 10 i = 1, icfct
         do 20 j = 1, icp
            do k = 1, icp
              if (idcf(j,i).eq.id(k)) goto 20
            end do 
            goto 10
20       continue
c                             matches earlier assemblage, return.
         return
10    continue
c                             unique (and bounding) assemblage:   
      icfct = icfct + 1
      idid = idid + 1
      if (icfct.gt.j9) call error (204,sum,j9,'ASSDC')
      jasmbl(icfct) = 0
      icase = 1
c                             check if assmblage bounds a constrained
c                             composition:      
      if (jbulk.ne.0) then
         call bulkck (id,jbfg)
         jasmbl(icfct) = jbfg
         if (jbfg.ne.0) icase = 2
      end if 
c                             assign the assemblage: 
      do i = 1, icp
         idcf(i,icfct) = id(i)
      end do 
c                             do miscibility check:
      call miscib (id,ivar,imyn)
      iasmbl(icfct) = ivar

      goto (99), jpot
c                             print potentials:
      call prtpot

99    end

      subroutine chkass (jdv,ivd,ivi,iste)
c-----------------------------------------------------------------------
c chkass determines if a divariant assemblage (jdv) generated by
c newass has already been identified. there are three possible cases:
c in the reaction as determined by routine balanc.
c (1) the assemblage is unique: the indices of the phases, the initial
c     conditions of stability, and the coordinate frame edge are recorde
c (2) the assemblage has already been defined and it's stability field
c     determined: quit
c (3) the assemblage has already been defined but it's stability field
c     has not been determined: reset the initial conditions of stability
c     to shorten traverses by the routine search, if possible. 

c input:  jdv(icp) - a vector containing the indices of phases in a
c                    divariant assemblage.
c         ipart    - 1 if testing original divariant assemblages, 2 if
c                    testing new divariant assemblages.
c         index    - index of the divariant assemblage in array idcf
c                    (ipart=1) or idns (ipart=2) currently being tested.
  
c referenced by: newass
c references to: no external references
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,l2=5,j9=40000)

      integer jdv(k5+1), icrap

      common/ cst6 /icomp,istct,iphct,icp/ cst5 /v(l2),tr,pr,r,ps
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst96 /idcf(k5,j9),icfct
     *      / cst65 /vt(j9),vti(j9),jok(j9),kok(j9),index
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9    
     *      / cst300 /cblk(k5),atw(k1),jbulk
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)

      save icrap

      data icrap /0/
c-----------------------------------------------------------------------

      do 40 i = 1, icfct
         do 50 j = 1, icp
            do 60 k = 1, icp
60             if (jdv(k).eq.idcf(j,i)) goto 50
c                                 no match with assemblage i
            goto 40
c                                 
50       continue
c                                 jdv matches the ith assemblage
c                                 first check if within known stability
c                                 field of the assemblage
         if (iste.lt.kok(i)) then 
            goto 99
         else if (iste.eq.kok(i)) then 
            if (iste.lt.3) then
               if (vt(i).ge.v(ivd)-dv(ivd)) goto 99
            else 
               if (vt(i).le.v(ivd)+dv(ivd)) goto 99
            end if 
c                                 it extends known stability field:
c                                 if the assemblage hasn't been tested
c                                 then reset start conditions
            call sreset (jok(i),iste,vt(i),v(ivd),vti(i),v(ivi))
            goto 99
         else if (i.gt.index) then 
            call sreset (jok(i),iste,vt(i),v(ivd),vti(i),v(ivi))
            goto 99
         end if
c                                 otherwise save as a new assemblage
         goto 100
  
40    continue
c                                 the assemblage is unique, store
c                                 parameters:
100   goto (99),icrap
c                                 do bulk composition check
      if (jbulk.ne.0) then
         call bulkck (jdv,jbfg)
         if (jbfg.eq.0) goto 99
      end if 

      icfct = icfct + 1

      if (icfct.gt.j9) then
         call warn (205,r,j9,'CHKASS')
         icfct = j9
         icrap = 1
      end if 

      do i = 1, icp
         idcf(i,icfct) = jdv(i)
      end do 

      vt(icfct) = v(ivd)
      vti(icfct) = v(ivi)
      jok(icfct) = iste
      kok(icfct) = iste

99    end

      subroutine coface (iovd,iovi,odiv,iste,irend)
c----------------------------------------------------------------------
c once nohold has established that a stable equilibrium occurs on
c the boundary of a diagram, coface traces the equilibrium condiitons
c within the diagram.  if the equilibrium terminates at an
c invariant point coface traces the equilibria which emanate
c from the invariant point and any invariant points encountered
c subsequently.  coface does not detect indifferent points.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k6=1000,k7=k5+1,k8=k5+2,l2=5)

      common/ cst61  /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst29  /vip(l2,k2),ipid(k2,k8),ipct/ cst32 /ptx(k6),ipt2
     *      / cst7   /iflag/ cst80 /irchk(k2)/ cst5 /v(l2),tr,pr,r,ps
     *      / cst23  /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst31  /vn(k2,k7),irct,ird/ cst25 /vnu(k7),idr(k7),ivct
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9         
     *      / cst105 /ibug(k2)/ cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst49  /value,ismax,igot
c-----------------------------------------------------------------------
      ivi =iovi                                                  
      ivd=iovd
      jflg= 1

      call balanc (b,idv,idphi,ier)
c                                 singular concentration matrix
      if (ier.eq.1) call error (68,v(1),ier,'COFACE')
c                                 call newass to define new (pseudo-)
c                                 divariant assemblages generated by
c                                 the reaction defined by balanc:
      call newass (b,idv,idphi,ivd,ivi,iste)
c                                 if the rxn is univariant (ivar =1)
c                                 continue.
      if (ivar.ne.2) ivar = 1
      if (ivar.eq.1) goto 80 
c                                 the rxn is pseudo-univariant continue
c                                 only if isudo= 0 (input).
      goto (9999),isudo
c                                 call assir to determine if the
c                                 reaction has already been identified
c                                 and if so quit, there is a small
c                                 possibility that as a result of
c                                 quitting at this stage some equilibria
c                                 will not be identified, if this is a
c                                 significant concern the following
c                                 computed goto should be deleted.
80    goto (90), icopt
c                                 allow duplicated reactions in 
c                                 mixed variable diagrams:
      call univeq (ivd,ier)
c                                 assign reaction:
      call assri (ier)

      goto (9999), ier

      call delrxn 

      vip(1,ird) = v(1)
      vip(2,ird) = v(2)
      vip(3,ird) = v(3)
      vip(4,ird) = v(4)
      vip(5,ird) = v(5)

      call outrxn (0,0)

      goto 9999

90    call assir 
c                                 skip equilibria that have already
c                                 been defined on an edge, there is some
c                                 possibility that this may result in an
c                                 an incomplete diagram, but it saves
c                                 some time and redundancy in output.
c                                 check if the univariant rxn has alread
c                                 been found on the edge of a diagram.
c                                 set flags irchk and ier:
c      if (irchk(ird).eq.1) goto 9999
c                                 call svrend to see if the end point 
c                                 has already been found, if it has ier
c                                 is 1, and traverse can be skipped. 
      call svrend (ird,irend,ier)
      goto (9999),ier

      irchk(ird) = 1
      call univeq (ivd,ier)
c                                 get deltas for dependent extensities:
      call delrxn
c                                 if univeq fails on a bounding edge
c                                 write error message and return:
      goto (9000,9000),ier
c                                 save the index of the c+1th phase
c                                 as iophi.
      iophi =idphi
c                                 factor the concentration matrix:
      call pivots (ier)
      if (ier.eq.1) call error (69,v(1),ier,'COFACE')
c                                 set the increment for the iv
      div = odiv
c                                 initialize counters
      ipt2 = 0
      icter = 0
c                                 assign the 1st point
      call assptx
c                                 follow the equilibrium
60    ikwk = 0 
      jflg = 0
      call sfol1 (ivd,ivi,ier,div,ikwk,irend)
c                                 sfol1 returns ier = 1
      goto (70,70),ier
      ivi = iovi
      ivd = iovd
      if (iflag.eq.0) then 
c                                 iflag = 0 the univariant
c                                 curve was traced from one
c                                 edge to another, reset
c                                 the starting conditions
c                                 and set flag jer for newhld
         goto 9999
      else if (iflag.eq.1) then
c                                 an ip was found:
         goto 10
      end if 

70    call switch (div,ivi,ivd,jer)
      goto (75),jer
      icter = icter + 1
      if (icter.lt.4) goto 60
      
75    call warn (10,v(1),ier,'COFACE')

      call outrxn (ipct,2)
      ibug(irct) = 1
      ivi =iovi
      ivd=iovd
c                                 return on error
      goto 9999                
c                                 a new invariant point has been
c                                 encountered, generate the new 
c                                 univariant fields.
c                           iopct=the total number of ip's encountered
c                                 before each call to sfol2.
c                           inpct=the total number after each call to 
c                                 sfol2.
10    iopct = ipct
      call sfol2 (iovi,iovd,iopct,irend)
      if (iopct.eq.ipct) goto 9999
      iopct = iopct + 1
      inpct = ipct
c                                 loop for new invariant points
30    do i = iopct, inpct
         call sfol2 (iovi,iovd,i,irend)
      end do 
      if (ipct.eq.inpct) goto 9999
      iopct = inpct + 1
      inpct = ipct
      goto 30
c                                 error in univeq:
9000  call warn (79,v(1),ird,'COFACE')
      ipt2 = 0
      call outrxn (ipct,2)
      ibug(irct) = 1

9999  iflg1= 0 
      goto (999),jflg
      goto (999),io3p

      write (n3,1020) 
      
1020  format (' Network traced, resuming boundary search.',/)

999   end

      subroutine delvar (dv,iflag,iflg1)
c-----------------------------------------------------------------------
c delvar determines if a variable increment 'dv' should
c be incremented or decremented on a traverse to locate
c to locate conditions where an equilibrium is metastable with respect
c to exactly one phase (iflag = 1). the flag iflg1 indicates if the
c equilibrium was stable (iflg1 = 0) or metastable (iflg1 = 1) at the
c last conditions tested.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
c                            the equilibrium is stable (iflag = 0):
      goto (20,20),iflag
c                            and was stable before (iflg1 = 0) then
      goto (10),iflg1
c                            continue looking in the same direction:
      return
c                            and was metastable before (iflg1 = 1) then
c                            switch sign and halve dv:
10    dv = -dv/2d0
      iflg1 = 0
      return
c                            the equilibrium is metastable (iflag = 2):
20    goto (30),iflg1
c                            and was stable before (iflg1 = 0) then
c                            switch sign and halve dv:
      dv = -dv/2d0
      iflg1 = 1
c                            and was metastable before (iflg1 = 1) then
c                            continue looking in the same direction:
30    end
 
      subroutine factor (a,n,ipvt,ier)
c-----------------------------------------------------------------------
c factor is a subroutine which calculates the triangular
c decompositions of the matrix 'a'. factor is modified from
c the subroutine of the same name given by conte and de boor
c in 'elementary numerical analysis', mcgraw-hill, 1980.
c factor uses scaled partial pivoting.
c
c input     a- an n by n array containing the elements of matrix a.
c           n- the dimension of the matrix a.
c output    a- an n by n array containing the upper, u, and lower, l,
c              triangular decompositions of input matrix a.
c        ipvt- a vector indicating that row ipvt(k) was used to
c              eliminate the a(n,k).
c         ier- a flag, zero if a is of rank = n, and 1 if a is of
c              lower rank.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)

      dimension a(k5,k5),d(k5),ipvt(k5)
c-----------------------------------------------------------------------
      ier = 0
c                            initialize ipvt,d
      do i = 1, n
         ipvt(i) = i
         rmax = 0d0
         do j = 1, n
            rmax = dmax1(rmax,dabs(a(i,j)))
         end do 
c                            ax = b is singular if rmax = 0
         if (dabs(rmax).lt.1d-5) goto 9000
         d(i) = rmax
      end do 
c                            begin decomposition: 
      nm1 = n-1
c
      do i = 1, nm1
         ip1 = i+1
c                            determine pivot row (istr).
         rmax = dabs(a(i,i))/d(i)
         istr = i
         do 40 j = ip1, n
            tmax = dabs(a(j,i))/d(j)
            if (tmax.le.rmax) goto 40
            rmax = tmax
            istr = j
40       continue

         if (dabs(rmax).lt.1d-5) goto 9000
c                            if istr gt i, make i the pivot row
c                            by interchanging it with row istr.
         if (istr.le.i) goto 50
         j = ipvt(istr)
         ipvt(istr) = ipvt(i)
         ipvt(i) = j
         temp = d(istr)
         d(istr) = d(i)
         d(i) = temp

         do j = 1, n
            temp = a(istr,j)
            a(istr,j) = a(i,j)
            a(i,j) = temp
         end do 
c                            eliminate x(k) from rows k+1,...,n.
50       do j = ip1,n
         a(j,i) = a(j,i)/a(i,i)
         ratio = a(j,i)
            do k = ip1,n
            a(j,k) = a(j,k)-ratio*a(i,k)
            end do 
         end do
      end do 

      if (dabs(a(n,n)).lt.1d-5) ier = 1
      return
c                           algoritmic singularity.
9000  ier = 1
 
      end

      subroutine flipit (ddv,vst,ivd,ist,inow,jer)
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      character*8 names

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst5 /v(l2),x(4)/ cst7 /iflag/ cst6 /icomp,istct,iphct,icp
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst38 /iedge
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst9 /vmax(l2),vmin(l2),dv(l2)/ cst8 /names(k1)
c-----------------------------------------------------------------------
c                                 initialization
      call gall 
      call asschk
c                                 the goto 20 is a dummy to avoid
c                                 compiler errors, change 90 to 20 to
c                                 get backup search, ha ha
      if (iflag.ne.0.and.v(ivd).eq.vst.and.ist.eq.inow) goto 90

      goto (10,20,20),iflag
c                                 the assemblage is stable, return
      jer = 0
      return
c                                 the assemblage is metastable with
c                                 respect to only one phase:
10    jer = 1
      return
c                                 the assemblage is metastable with
c                                 respect to >1 phase, therefore
c                                 it must be stable at some condition
c                                 between the present condition and the
c                                 conditions of the univariant
c                                 equilibrium that generated the
c                                 assemblage. flipit does a reverse
c                                 traverse to find this condition.
c                                 initialize:
20    iflg1 = 1
      ddv = -ddv
c                                 begin search:
100   v(ivd) = v(ivd)+ddv
c                                 out of range?:
      goto (120,120,110,110),ist
c
110   if (v(ivd).gt.vmax(ivd)) v(ivd) = vmax(ivd)
      if (v(ivd).gt.vst) goto 130
      ddv = dabs(ddv)/2d0
      v(ivd) = vst
      iflg1= 0
      goto 100
c
120   if (v(ivd).lt.vmin(ivd)) v(ivd) = vmin(ivd)
      if (v(ivd).lt.vst) goto 130
      ddv = -dabs(ddv)/2d0
      v(ivd) = vst
      iflg1 = 0
      goto 100
c                                 calculate phase energies:
130   call gall 
c                                 test stability of the assemblage:
      call asschk
c                                 check if search is in range:
      goto (150,150,140,140),ist
140   if (v(ivd).ge.vmax(ivd).and.iflag.gt.0) goto 60
      goto 160
150   if (v(ivd).le.vmin(ivd).and.iflag.gt.0) goto 60
c                                 iflag=1, found a stable equilibrium
c                                 in bounds:
      if (iflag.eq.1) write (n3,*) ' flipit worked please tell me!'
c
160   goto (10),iflag
c                                 refine increment if necessary:
      call delvar (ddv,iflag,iflg1)
c                                 check if increment has been refined
c                                 beyond the tolerance which is
c                                 arbitrarily set at 1.d-8.
      if (dabs(ddv).lt.1d-8) goto 9000
c                                 increment conditions again:
      goto 100
c                                 next traverse.
60    call warn (76,vst,ist,'FLIPIT')

      goto 90
c                                 write warning message:
9000  write (n3,1010) ivd,dv(ivd),(names(idv(j)),j = 1, icp)
90    jer = 2
c                                 done:
1010  format (/,' **warning ver045** FLIPIT: > 1 equilibrium',
     *          ' occurs within the',/,' minimum search increment for',
     *          ' variable: ',i1,', this often occurs as YCO2 => 1',
     *          ' or => 0, you may be able to correct this',/,
     *          ' by reducing the default increment for this variable',
     *          ' (',g12.3,') in the computational option file.',/,
     *          ' Equilibria involving the following assemblage may',
     *          ' not be delineated:',/,1x,7(1x,a8)) 
      end

      subroutine outdel
c-----------------------------------------------------------------------
c output deltas of dependent extensities for a reaction
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5,h5=7,h6=500)

      character*5 cname,vname*8,gname*78,xname*18,exten(2)*7,names*8

      common/ csta2  /xname(k5),vname(l2),gname(l2)/ csta4 /cname(k5)
     *      / cst21  /du(2),dv(2),jds(h5),ifr 
     *      / cst88  /vu(2,k1),jprct,jmct
     *      / cst6   /icomp,istct,iphct,icp
     *      / cst201 /vuf(2),vus(h5),iffr,isr
     *      / cst25  /vnu(k5+1),idr(k5+1),ivct/ cst8 /names(k1)
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9

      save exten 
c                                     this is a bullshit trick and
c                                     will cause errors of someone
c                                     uses a function other than G.
      data exten/'-V(j/b)','S(j/k)'/
c                                     composant stoichiometry:
 
      do i = 1, isat 
         write (n3,1000) cname(icp+i),vus(i),names(jds(i))
      end do 
c                                     fluid stoichiometry:
      if (ifyn.eq.0) then 
         do i = 1, 2
            if (iff(i).ne.0) write (n3,1010) names(i),vuf(i)
         end do 
      end if 
c                                     mobile components:
      do i = 1, jmct
         write (n3,1020) cname(jprct+i),du(i),vname(3+i)
      end do
c                                     dependent extensities 1 and 2, 
c                                     normally entropy and negative 
c                                     volume:
      do i = 1, 2
         write (n3,1030) exten(i), dv(i),vname(i)
      end do

1000  format (10x,'Delta(',a7,') =',g9.3,1x,
     *            '(saturated composant=',a8,')')
1010  format (10x,'Delta(',2x,a5,') =',g9.3,1x,
     *            '(saturated phase component)')
1020  format (10x,'Delta(',a7,') =',g9.3,1x,
     *            '(dependent conjugate of ',a8,')')
1030  format (10x,'Delta(',a7,') =',g9.3,1x,
     *            '(dependent conjugate of ',a8,')')
      end 

      subroutine delrxn
c-----------------------------------------------------------------------
c delrxn computes the change in the dependent extensities for a reaction
c these are returned in cst21, delrxn is always called after UNIVEQ, so
c the free energy at the present conditions is zero, if called 
c independently, then grxn must be called first.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k4=22,l2=5,h5=7,h6=500,k10=200,
     *           k5=12,k7=k5+1)

      common/ cst5  /v(5),tr,pr,r,ps/ cst21 /du(2),dv(2),jds(h5),ifr 
     *      / cst87 /delt(l2),dtol,utol,ptol/ cst102 /ivfl
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst201 /vuf(2),vus(h5),iffr,isr/ cst208 /ifct,idfl
     *      / cst25  /vnu(k7),idr(k7),ivct/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst31  /vn(k2,k7),irct,ird
     *      / cst88  /vu(2,k1),jprct,jmct
     *      / cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
c-----------------------------------------------------------------------
c                                         save dependent extensities

c                                         it is necessary to calculate
c                                         the initial G cause iptol may
c                                         be large relative to the 
c                                         change in G due to delt
      og = grxn(ird)

      do i = 1, 2
         v(i) = v(i) + delt(i)   
         dv(i) = (og - grxn(i)) / delt(i)
         v(i) = v(i) - delt(i) 
      end do 
c                                         get deltas on mobile comps
      if (jmct.gt.0) then

         do i = 1, jmct
            du(i) = 0d0
         end do

         do i = 1, ivct 
            do j = 1, jmct 
               du(j) = du(j) + vnu(i) * vu(j,idr(i))
            end do 
         end do
      end if 
c                                         do inverse projection to 
c                                         get stoichiometries of 
c                                         constrained components:
      if (isat.gt.0) then
         if (isat.gt.1) then  
            do j = 1, isat - 1
               do i = j + 1, isat
                  vus(j) = vus(j) - vus(i) * vs(j,idss(i))
               end do 
            end do
         end if 
c                                         now do the effects of 
c                                         the inverse projection on
c                                         saturated phase components:
         if (ifct.ne.0) then 
            do i = 1, isat
               do j = 1, ifct
                  vuf(j) = vuf(j) - vus(i) * vf(j,idss(i))
               end do
            end do 
         end if 
c                                         now do the effects of 
c                                         the inverse projection on
c                                         mobile components:
         if (jmct.ne.0) then 
            do i = 1, isat 
               do j = 1, jmct
                  du(j) = du(j) - vus(i) * vu(j,idss(i))
               end do 
            end do 
         end if 

      end if 
c                                         see if coefficients should 
c                                         change sign:
      invert = 1
      if (ivfl.eq.1.and.dv(1).gt.0d0) then 
         invert = 0
      else if (ivfl.eq.2.and.dv(2).gt.0d0) then 
         invert = 0
      else if (ivfl.gt.2) then 
         if (du(ivfl-3).gt.0d0) then 
c                                         diagram must have a mu variable 
            invert = 0
         else if (du(ivfl-3).eq.0d0) then
            if (iv(2).eq.3) then
               if (vuf(2).lt.0d0) invert = 0
            else if (iv(1).eq.3) then
               if (vuf(2).lt.0d0) invert = 0
            else if (du(iv(1)-3).gt.0d0) then
               invert = 0
            end if
         end if
      end if

      if (invert.eq.1) then 
         do i = 1, ivct 
            vnu(i) = -vnu(i)
c                                         why different signs?
            vn(ird,i) = vnu(i)
         end do

         do i = 1, isat
            vus(i) = -vus(i)
         end do 

         do i = 1, ifct
            vuf(i) = -vuf(i)
         end do

         do i = 1, jmct 
            du(i) = -du(i)
         end do

         dv(1) = -dv(1)
         dv(2) = -dv(2)

      end if   
c                                         save id's of saturated comps
      do i = 1, isat
         jds(i) = idss(i)
      end do

      end     

      subroutine header (n,icopt)
c-----------------------------------------------------------------------
c header writes the graphics file header for schreinemakers diagrams
c to lun n4
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,l2=5,h5=7,h6=500,h9=30)

      character*8 names,vname,gname*78,xname*18,title*72,fname*10
      character*5 cname

      common/ cst8 /names(k1)/ cst6 /icomp,istct,iphct,icp/ cst79 /iso
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst18 /jvol(2),ivol
     *      / cst60 /ipoint,imyn
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)/ csta7 /fname(h9)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ csta8 /title
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst102 /ivfl/ cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ csta4 /cname(k5)
c-----------------------------------------------------------------------
c                              value to be read as icopt

      if (n.eq.n4) then 
         write (n,1020) icopt

      else if (n.eq.n5) then 

         write (n,1020) icopt

         write (n,1020) jbulk,ivol

      end if 
c                              write phase and solution count:
      write (n,1020) iphct, iso 
c                              write phase names:
      write (n,1010) (names(i), i = 1, iphct)
c                              write solution phase flags:
      write (n,1060) (ikp(i), i = 1, iphct)  
c                              write solution names:
      if (iso.ne.0) write (n,1040) (fname(i), i = 1, iso)
c
      write (n,1150) title                                
c                                  write saturated and buffered
c                                  component names:
      if (isyn.eq.0) then 
         write (n,1070) (cname(i+icp), i= 1, isat)
      else 
         write (n,1150) ' '
      end if 

      write (n,1080) vname(ivfl)
      write (n,1020) ipot
      write (n,1030) (vmax(jv(i)), vmin(jv(i)), i= 1, ipot)
      write (n,1000) (ig(jv(i)), gname(jv(i)), i= 1, ipot)

1000  format (i2,a78)
1010  format (10a8)
1020  format (20(i5,1x))
1040  format (8a10)
1030  format (6(d11.5,1x))
1060  format (26(i2,1x))
1070  format ('Component saturation hierarchy: ',7(a5,1x))
1080  format ('Reaction equations are written such that the high ',
     *         a8,/,'assemblage is on the right of the = sign')
1150  format (a72)
      end

      function lchk (lphi)
c-----------------------------------------------------------------------
c a function subprogram to determine if a phase, lphi, lies
c below a g-x plane defined by the vertices idv.
c lchk has a value of 0 if the plane is stable, and 1 if not.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst62 /blim(l2),ulim(l2),dg/ cst12 /cp(k5,k1)
     *      / cst6 /icomp,istct,iphct,icp/ cst2 /g(k1)
     *      / cst87 /delt(l2),dtol,utol,ptol
c-----------------------------------------------------------------------
      call uproj

      do i = 1, icp
        b(i) = gproj(idv(i))
      end do 

      g(lphi) = gproj(lphi)
      lchk = 0

      call subst (a,ipvt,icp,b)

      gphi = 0d0
      do j = 1, icp
         gphi = gphi + cp(j,lphi) * b(j)
      end do

      if (gphi.lt.g(lphi)) goto 99

      lchk = 1

99    end

      subroutine newass (b,idv,idphi,ivd,ivi,iste)
c-----------------------------------------------------------------------
c newass determines the (pseudo-) divariant assemblages generated by
c a reaction based on the stoichiometric coefficients of the phases
c in the reaction as determined by routine balanc.

c input:  idv(icp) - a vector containing the indices of phases in a
c           divariant assemblage.
c         idphi - the index of the phase which is involved with the
c           the divariant assemblage (idv) in a univariant equilibrium.
c         b(icp) - the stoichiometric coefficient of the phases inedxed
c           by idv in the reaction of these phases to form the phase
c           idphi. by convention the stoichiometric coeffeicient of
c           idphi is -1.0.

c output: jdv(icp) - a vector containing the indices of phases in a
c           divariant assemblage generated by the reaction between the
c           phases idphi and idv this vector is used by the routine
c           chkass. this assemblage is stable on the opposite 'side' of
c           the univariant equilibrium from the assemblage idv.

c referenced by: coface
c references to: chkass

c algorithm: phases with null or negative coefficients are always
c   compatible with eachother on the 'side' of a univariant reaction
c   which idv is metastable. the k phases with positive coefficients
c   define a chemographic simplex surrounding phases with null or
c   negative coefficients. consequently a univariant reaction will
c   generate k new divariant assemblages consisting of k-1 of the phases
c   in the original assemblage (idv) and all the phases with null or
c   negative coefficients. if only one phase has a positive coefficient
c   this phase is unstable and only one new assemblage will be generated
c   which consists of all the other phases (including idphi).
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)

      dimension idv(k5),b(k5+1),idpos(k5),jdv(k5+1)

      common/ cst6  /icomp,istct,iphct,icp
c-----------------------------------------------------------------------
c                                 initialize counter:
      ipos = 0

      do 10 i = 1, icp
c                                 is the b(i) coefficient zero/negative?
         if (b(i).lt.1d-5) goto 20
c                                 no, save in the array idpos and count.
         ipos = ipos + 1
         idpos(ipos) = idv(i)
         goto 10
c                                 yes, load into the array jdv
20       jdv(i-ipos) = idv(i)
10    continue
c                                 load idphi
      ineg = icp - ipos+1
      jdv(ineg) = idphi
c                                 generate the divariant assemblages
      do j = 1, ipos
c                                 make ipos permutations of ipos - 1 phases
         jtic = 0
         do k = 1, ipos
            if (k.ne.j) then
               jtic = jtic + 1
               jdv(ineg+jtic) = idpos(k)
            end if 
         end do            
c                                 send to the routine chkass to store th
c                                 the assemblage if unique.
         call chkass (jdv,ivd,ivi,iste)
      end do
c                                 done
      end

      subroutine onedim
c----------------------------------------------------------------------
c for phase diagrams with one extrinsic variable (iv1) onedim sorts
c equilibria, identified by newhld and located by coface, by veq(iv1).
c the sorted equilibria are then output to units n3 and n4 as requested.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,k8=k5+2,l2=5)

      common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct/ cst13 /ir(k2)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst31 /vn(k2,k7),irct,ird
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9    
c-----------------------------------------------------------------------
c                             any rxns to be sorted?
      if (irct.eq.1) goto 30
c                             initialize:
      ist = 2

      do i = 1, irct
         ir(i) = i
      end do 
c                             begin sort loop:
20    ist1 = ist - 1
      itemp = ir(ist1)
      vst = vip(iv1,itemp)

      do i = ist, irct
         if (vip(iv1,ir(i)).le.vst) then
            ir(ist1) = ir(i)
            ir(i) = itemp
            itemp = ir(ist1)
            vst = vip(iv1,itemp)
         end if 
      end do

      ist = ist + 1
      if (ist.gt.irct) goto 30
      goto 20
c                             end of sort loop
30    if (irct.eq.0) goto 99
c                             output invariant reactions and conditions:
      call outirn 

99    end
            
      subroutine outier
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k8=k5+2,l2=5)

      character vname*8, gname*78, xname*18, text(400)*1

      common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct
     *      / cst24 /ipot,jv(l2),iv(l2)
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst81 /ipchk(k2),ip2(k2),icp2,ifake
     *      / cst103 /isec,icopt,ifull,imsg,io3p
c-----------------------------------------------------------------------
c                                 write end of ptx data for univariant
c                                 equilibria marker to graphics file.
      if (io4.ne.1) then
         write (n4,*) 1,1,1,1,1,1d0
         write (n6,1000) ipct,icp2
         do j = 1, ipct
c                                 n4:
            write (n6,1000) j, ivarip(j), (ipid(j,i),i = 1, icp2)
            write (n6,1060) (vip(jv(i),j), i = 1, ipot)
         end do 
      end if 
c                                 write ip names and locations to
c                                 units n6 and n3.
      if (ipct.eq.0) goto 99

      if (io3.ne.1) then 
         write (n3,1050)
         write (n3,1030)
         if (io3p.ne.1) then 
            do j = 1, ipct
         
               call iptext (text,iend,j,icp2)
               write (n3,1080) j, ivarip(j), (text(i), i= 1, iend)
               write (n3,1020)
               write (n3,1040) (vname(jv(i)), vip(jv(i),j), i= 1, ipot)
            end do 
         end if 
      end if
    
1000  format (14(i5,1x))
1020  format (15x,'occurs at:')
1030  format (' (pseudo-) invariant points are summarized below:')
1040  format (25x,a8,'=',g12.6)
1050  format (/,80('-'),/)
1060  format (5(1x,g11.5))
1080  format (/,' (',i5,'-',i1,') ',380a1)

99    end

      subroutine outirn
c-----------------------------------------------------------------------
c outirn writes the identity and equilibrium conditions of invariant
c equilibria to units n3 and n4 after they have been sorted by onedim.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,k8=k5+2,l2=5,h9=30)

      character*8 vname, gname*78, xname*18, rxnstr*400

      common/ cst31 /vn(k2,k7),irct,ird/ cst3 /x(k1,k5)
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst6  /icomp,istct,iphct,icp/ cst104 /rxnstr(k2)
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct/ cst5  /v(l2),q(4)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst35 /irv(k2)/ cst13 /ir(k2)/ csta1 /irnms(k2,k7)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst83 /ig(l2)
c-----------------------------------------------------------------------
c                             write invariant reactions:
      do 50 i = 1, irct
         j = ir(i)
c                             the index on these two guys used to be j
c                             must be wrong? if so why did it work 
c                             before? changed to i. you know sometimes
c                             i'm so stupid, i just can't believe it.
         k = irv(j)
         l = ivarrx(j)
         goto (10),io3
c                             output to print file:
         if (l.eq.1) then
            write (n3,1030) j, l, rxnstr(j)
         else
            write (n3,1040) j, l, rxnstr(j)
         end if 

         write (n3,1070) vname(iv1), vip(iv1,j), vname(iv2), vip(iv2,j)

10       write (n3,1160)
c                             output to graphics file
         goto (50), io4

         write (n4,1000) j, k, l, vip(iv1,j), (irnms(j,l), l= 1, k)
         write (n4,1050) (vn(j,l), l = 1, k)

50    continue

1000  format (i5,1x,2(i3,1x),g12.5,1x,6(i3,1x))

1030  format (' The equilibrium of the invariant reaction:',//,
     *       ' (',i4,'-',i1,') ',a400)
1040  format (' The equilibrium of the pseudoinvariant reaction:',//,
     *       ' (',i4,'-',i1,') ',a400)
1050  format (12x,6(g9.2,1x))
1070  format (/,' occurs at ',a8,'=',g12.6,' and ',a8,'=',g12.6)
1160  format ('      ----------------------------------------')

      end

      subroutine outlst
c----------------------------------------------------------------------
c outlst writes cumulative lists of the univariant and invariant
c equilibria identified during a single execution of the program
c to the graphics and print files. outlst also calls the subroutine
c out5 to write a permanent list (on n5) of equilibria for multiple
c executions.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,k8=k5+2,l2=5)

      character*8 names,cname*5,rxnstr*400

      common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst81 /ipchk(k2),ip2(k2),icp2,ifake
     *      / cst31 /vn(k2,k7),irct,ird/ cst105 /ibug(k2)
     *      / cst8  /names(k1)/ cst29 /vip(l2,k2),ipid(k2,k8),ipct
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(k5)
     *      / cst80 /irchk(k2)/ cst104 /rxnstr(k2)
c------------------------------------------------------------------------
c                             printing to n4 is suppresses for the
c                             erlanger system 
      jbug = 0
c                             invariant point lists:
      if (ipct.eq.0.or.io3.eq.1) goto 80
c                             output print file heading:
c                             output ip lists:
      write (n3,1030)
      do 70 i = 1, ipct
70       write (n3,1150) i, ivarip(i), (names(ipid(i,k)), k = 1, icp2)
c                              write end of invariant point list
c                              marker to graphics file.
80    goto (140),io3
      write (n3,1160)
c                              write list of stable univariant
c                              equilibria to units n3 and n4.
      write (n3,1070)

      do i = 1, irct
         if (ibug(i).ne.0) jbug = 1
         write (n3,1000) i,ivarrx(i),rxnstr(i)

      end do    
c                              if itic= 0 no equilibria were identified
c                              goto 9000 to write no eq mess.
      if (irct.eq.0) goto 9000
      write (n3,1160)
c                              jbug ne 0, one or more equilibria
c                              may not have completely defined, write
c                              list of possible failures:
140   if (jbug.ne.0) then 
         write (n3,1160)
         write (*,1160)
         write (n3,1190)
         write (*,1190)
         do 141 i = 1, irct
            if (ibug(i).eq.0) goto 141
            write (*,1000) i,ivarrx(i),rxnstr(i)
            write (n3,1000) i,ivarrx(i),rxnstr(i)
141      continue
         write (n3,1160)
         write (*,1160)
      end if 

      goto 99
c                              write no equilibrium messages:
9000  goto (100),io3
      write (n3,1180)
100   goto (99),io4
      write (n4,1180)

1000  format (' (',i4,'-',i1,') ',a300)
1030  format (' (pseudo-) invariant points are summarized below:',/)
1060  format (9(a8,1x))
1070  format (' (pseudo-) univariant equilibria are summarized ',
     *       'below:',/)
1150  format (' (',i4,'-',i1,') ',12(a8,1x))
1160  format (//,80('-'),//)
1180  format (' no univariant or invariant equilibria occur.')
1190  format (' WARNING!! The stability fields of the following',
     *        ' equilibria may',/,' have been entirely or',
     *        ' partially skipped in the calculation: ',/)
99    end

      subroutine iptext (text,iend,jp,icp2)
c-----------------------------------------------------------------------
c subprogram to write a text label for an invariant point 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k8=k5+2,l2=5,h9=30)

      character*8 names, text(400)*1, fname*10

      common/ cst8  /names(k1)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct/ csta7 /fname(h9)

c                             dump the names into the array text
      ist = 1

      do 20 i = 1, icp2
         id = ipid(jp,i)
         if (ikp(id).eq.0) then
c                             simple compound:
            read (names(id),1000) (text(j), j = ist, ist + 7)
            ist = ist + 8
         else  
c                             solution phases:
            read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9)
            text(ist + 10) = '('
            read (names(id),1000) (text(j), j = ist + 11, ist + 18)
            text(ist + 19) = ')'
            ist = ist + 20
         end if 
         text(ist) = ' '
20       ist = ist + 1 
c                             now filter out blanks
      iend = 1
      do 40 i = 2, ist - 1
         if (text(i).eq.' '.and.text(i + 1).eq.' ') then 
            goto 40 
         else if (text(i).eq.' '.and.text(i + 1).eq.')') then 
            goto 40
         else if (text(i).eq.' '.and.text(i + 1).eq.'(') then 
            goto 40 
         end if 
      iend = iend + 1
      text(iend) = text(i)
40    continue

1000  format (20a1)

      end 

      subroutine rxntxt (iend,jend,text,alpha)
c-----------------------------------------------------------------------
c subprogram to write a text label for a reaction 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,h9=30)

      character*8 names, text(400)*1, fname*10,
     *            alpha(130)*1, string*400, rxnstr*400

      integer iplus(k5), iminus(k5), jplus(k5), jminus(k5), 
     *        ione(k5+1), jone(k5+1)

      common/ csta7 /fname(h9)/ cst8 /names(k1)
     *      / cst25 /vnu(k7),idr(k7),ivct
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst31 /vn(k2,k7),irct,ird
     *      / cst80 /irchk(k2)/ cst104 /rxnstr(k2)

      ip = 0 
      im = 0
      jend = 1
c                             load id's of phases with + or - coeffs 
      do 10 i = 1, ivct 
         if (vnu(i).gt.0d0) then 
            ip = ip + 1
            iplus(ip) = idr(i)
            jplus(ip) = i
         else 
            im = im + 1
            iminus(im) = idr(i)
            jminus(im) = i
         end if 
10    continue 

      do 60 i = 1, im 
         jone(i) = jminus(i)
60       ione(i) = iminus(i)
      do 70 i = 1, ip
         j = i + im
         jone(j) = jplus(i)
70       ione(j) = iplus(i)
c                             dump coefficients into string:
      write (string,1010) (vnu(jone(i)),i = 1, ivct)
c                             read alpha from string:
      jst = 7 + ivct * 10
      read (string,1000) (alpha(i), i = 1, jst)
c                             filter blanks and add right )
      do 50 i = 1, jst 
         if (alpha(i).eq.' '.and.alpha(i + 1).eq.' ') then 
            goto 50
         else if (alpha(i).eq.' '.and.alpha(i + 1).eq.',') then
            goto 50
         end if 
         alpha(jend) = alpha(i)
         jend = jend + 1
50    continue

      if (alpha(jend-1).eq.' '.or.alpha(jend-1).eq.',') 
     *                                       jend = jend - 1  
      alpha(jend) = ')'
c                             now dump the names into the array text
      ist = 1
      is = 1
80    do 20 i = is, im 
         id = ione(i)
         if (ikp(id).eq.0) then
c                             simple compound:
            read (names(id),1000) (text(j), j = ist, ist + 7)
            text(ist+8) = ' '
            ist = ist + 9
         else  
c                             solution phases:
            read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9)
            text(ist + 10) = '('
            read (names(id),1000) (text(j), j = ist + 11, ist + 18)
            text(ist + 19) = ')'
            text(ist+20) = ' '
            ist = ist + 21
         end if 
20    continue

      if (is.eq.1) then 
         text(ist) = '='
         text(ist + 1) = ' '
         ist = ist + 2
         is = im + 1
         im = ivct 
         goto 80 
      else 
         text(ist) = ' ' 
      end if 
c                             now filter out blanks
      iend = 1
      do 40 i = 2, ist
         if (text(i).eq.' '.and.text(i + 1).eq.' ') then 
            goto 40 
         else if (text(i).eq.' '.and.text(i + 1).eq.')') then 
            goto 40
         else if (text(i).eq.' '.and.text(i + 1).eq.'(') then 
            goto 40 
         end if 
      iend = iend + 1
      text(iend) = text(i)
40    continue
c                             this is redundant if the 
c                             reaction was found before. 
      write (rxnstr(ird),1000) (text(i),i=1,iend)

1000  format (400a1)
1010  format ('Alpha(',13(g9.3,', '))

      end 

      subroutine outrxn (ip,ier)
c-----------------------------------------------------------------------
c outrxn writes the identity and v1-v2 loci of univariant equilibria
c to units n3 and n4.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k6=1000,k7=k5+1,h9=30)

      character text(400)*1, alpha(130)*1, fname*10

      common/ cst31 /vn(k2,k7),irct,ird/ cst7 /iflag/ csta7 /fname(h9)
     *      / cst25 /vnu(k7),idr(k7),ivct/ cst32 /ptx(k6),ipt2
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst105 /ibug(k2)/ cst103 /isec,icopt,ifull,imsg,io3p      

      do i = 1, ivct 
         if (ikp(idr(i)).lt.0) goto 99
      end do 
c                                set bug flag to indicate that
c                                the equilibrium was at least 
c                                partially located.
      if (ier.lt.2) ibug(irct) = 0

      if (ifull.eq.0) then 
         call rxntxt (iend,jend,text,alpha)
      else  
         call fultxt (iend,text)
      end if 

      if (icopt.eq.3) goto 99
         
      if (imsg.eq.0) write (*,1110) ird,(text(i),i = 1, iend)

      goto (10), io3p

      write (n3,1120) ird,ivarrx(ird),(text(i),i = 1, iend)

      if (ifull.ne.0) goto 80

      write (n3,1130) (alpha(i),i = 1,jend)

      if (ipt2.lt.3) then 
         write (n3,*)
         goto 99
      end if 

      call outdel 

80    write (n3,*)
      write (n3,1000) (ptx(i), i = 1, ipt2)
      write (n3,*)
      if (ier.ne.0) goto 10
      goto (40),iflag
      goto 10
40    write (n3,1070) ip        
      write (n3,1040)

10    goto (99), io4

      write (n4,1080) ipt2,ird,ivar,ivct,(idr(i),i = 1, ivct)
      write (n4,1090) (vnu(i),i = 1, ivct)
      write (n4,1000) (ptx(i),i = 1, ipt2)

1000  format (1x,g12.6,1x,g12.6,2x,g12.6,1x,g12.6,2x,g12.6,1x,g12.6)
1040  format (/)
1070  format (' the equilibrium extends to invariant point (',i4,')')
1080  format (16(i4,1x))
1090  format (6(g12.4,1x)) 
1110  format (' finished with equilibrium (',i4,') ',400a1)
1120  format (' (',i4,'-',i1,') ',400a1)
1130  format (/,10x,90a1)

99    end

      subroutine pchk (dg,igo)
c-----------------------------------------------------------------------
c a subprogram which returns the difference in free energy
c beteen a g-x plane, defined by the assemblage idv, and a phase lphi.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst6  /icomp,istct,iphct,icp/ cst2 /g(k1)/ cst7 /iflag
     *      / cst12 /cp(k5,k1)/ cst87 /delt(l2),dtol,utol,ptol
c-----------------------------------------------------------------------
c                                 compute energies:
      igo = 0

      call uproj

      do 10 i = 1, icp
10       b(i) =gproj(idv(i))
      dg=gproj(idphi)
c                                 solve for chemical potentials:
      call subst (a,ipvt,icp,b)
c                                 compute energy difference:
      do 30 j = 1, icp
30       dg = dg - cp(j,idphi)*b(j)

      if (dabs(dg).lt.ptol) then
        igo = 1
        iflag = 1

        call ssaptx 
      else
        if (dg.gt.dtol) then
          iflag= 0

          call ssaptx
        else
          iflag= 1
        end if
      end if

      end

      subroutine pivots (ier)
c--------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst6 /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)
c--------------------------------------------------------------------
c                              calculate pivots for the matrix 'a'
      do 1 k = 1, icp
         do 2 j = 1, icp
 2          a(k,j) = cp(j,idv(k))
1     continue 

      call factor (a,icp,ipvt,ier)

      end

      subroutine schk (lphi)
c-----------------------------------------------------------------------
c schk my ___
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      common/ cst23 /a(k5,k5),u(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dg
     *      / cst87 /delt(l2),dtol,utol,ptol
     *      / cst2 /g(k1)/ cst12 /cp(k5,k1)/ cst7 /iflag
c-----------------------------------------------------------------------
c                                 determine chemical potentials
      iflag= 0
      do 10 i = 1, icp
10       u(i) = g(idv(i))

      call subst (a,ipvt,icp,u)
c                                 test phases, not=iophi, against the
c                                 assemblage idv(i),i = 1, icp
      do 20 i = istct, iphct
         if (i.eq.iophi.or.i.eq.lphi) goto 20
         dg = g(i)

         do 30 j = 1, icp
30          dg = dg - cp(j,i) * u(j)

         if (dg.gt.dtol) goto 20
c                                check that a phase is not metastable
c                                with respect to itself, this
c                                could be remedied by changing the 
c                                value of dtol.
         do 25 j = 1, icp
25          if (idv(j).eq.i) goto 20

         iflag = iflag + 1
         idphi = i
         goto (20,40,40), iflag

20    continue

40    end

      subroutine nschk 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      common/ cst23 /a(k5,k5),u(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dg
     *      / cst87 /delt(l2),dtol,utol,ptol
     *      / cst2 /g(k1)/ cst12 /cp(k5,k1)/ cst7 /iflag
c-----------------------------------------------------------------------
c                                 determine chemical potentials
      iflag= 0
      do 10 i = 1, icp
10       u(i) = g(idv(i))

      call subst (a,ipvt,icp,u)
c                                 test phases, not=iophi, against the
c                                 assemblage idv(i),i = 1, icp
      do 20 i = istct, iphct
         if (i.eq.iophi) goto 20

         dg = g(i)

         do 30 j = 1, icp
30          dg = dg - cp(j,i) * u(j)

         if (dg.gt.dtol) goto 20
c                                check that a phase is not metastable
c                                with respect to itself, this
c                                could be remedied by changing the 
c                                value of dtol.
         do 25 j = 1, icp
25          if (idv(j).eq.i) goto 20

         iflag = iflag + 1
         idphi = i
         goto (20,40,40), iflag

20    continue

40    end

      subroutine search (vst,vsti,ist,iste,ivi,ivd,div,ier)
c-----------------------------------------------------------------------
c search looks for an equilibrium around the periphery of the coordinate
c frame specified for a spd calculation. if a stable equilibrium is
c found ier = 0, if an error occurs ier =2, and if no equilibrium is
c found ier =1.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k5=12,l2=5)

      character*8 names

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1     
     *      / cst6 /icomp,istct,iphct,icp/ cst5 /v(l2),x(4)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst38 /iedge
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)/ cst8 /names(k1)
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst7 /iflag
     *      / cst49 /value,ismax,igot
c-----------------------------------------------------------------------
c                                 initialization
      iflg1 = 0
      ier = 0
      iugh = 0

      v(iv1) = vst
      v(iv2) = vst

      do 60 i = ist, iedge
         iste = i
c                                 set default dependent and independent
c                                 variables and increments for each edge
c                                 to be searched.
         goto (10,20,30,40),i
c                                 traverse 1.
10       ivi = iv2
         ivd = iv1
         ddv = dv(ivd)
         div = dv(ivi)
         v(ivi) = vmin(ivi)
         goto 80
c                                 traverse 2.
20       ivi = iv1
         ivd = iv2
         ddv = dv(ivd)
         div = -dv(ivi)
         v(ivi) = vmax(ivi)
         goto 80
c                                 traverse 3.
30       ivi =iv2
         ivd=iv1
         ddv = -dv(ivd)
         div = -dv(ivi)
         v(ivi) = vmax(ivi)
         goto 80
c                                 traverse 4.
40       ivi =iv1
         ivd=iv2
         ddv = -dv(ivd)
         div = dv(ivi)
         v(ivi) = vmin(ivi)
c                                 in some cases newass may generate
c                                 a metastable assemblage call flipit
c                                 to check. if the assemblage is
c                                 stable flipit will return jer = 0,
c                                 if it is metastable flipit tries a
c                                 reverse traverse and returns jer =1
c                                 if it is succesful, and jer =2 if not.
80       if (i.eq.ist) then
            v(ivd) = vst
            v(ivi) = vsti
         end if 

         call flipit (ddv,vst,ivd,ist,i,jer)

         if (jer.eq.2) then 
            write (*,1000) (names(idv(k)), k= 1, icp)
            write (*,1010) v
            ier = 1
            goto 99
         end if 

         vlast = v(ivd)

         goto (99),jer
c                                 begin search:
100      v(ivd) = v(ivd) + ddv
c                                 out of range?:
         goto (110,110,120,120),i

110      if (v(ivd).gt.vmax(ivd)) then
            v(ivd) = vmax(ivd)
         else if (i.eq.ist) then
            if (v(ivd).gt.vst) goto 130
            ddv = dabs(ddv)/2d0
            v(ivd) = vst
            iflg1= 0
            goto 100
         end if

         goto 130

120      if (v(ivd).lt.vmin(ivd)) then
            v(ivd) = vmin(ivd)
         else if (i.eq.ist) then
            if (v(ivd).lt.vst) goto 130
            ddv = -dabs(ddv)/2d0
            v(ivd) = vst
            iflg1= 0
            goto 100
         end if
c                                 calculate phase energies:
130   call gall 
c                                 test stability of the assemblage:
      call asschk

      if (iflag.eq.0) vlast = v(ivd)
c                                 iflag= 0 still stable
c                                 iflag=1 metastable with respect to one
c                                         phase indexed by idphi.
c                                 iflag=2 metastable with respect to
c                                         multiple phases refine the
c                                         search increment.
      goto (99),iflag
c                                 check if search is in range:
      goto (140,140,150,150),i
140   if (v(ivd).ge.vmax(ivd).and.iflag.eq.0) goto 60
      goto 160
150   if (v(ivd).le.vmin(ivd).and.iflag.eq.0) goto 60
c                                 refine increment if necessary:
160   call delvar (ddv,iflag,iflg1)
c                                 check if increment has been refined
c                                 beyond the tolerance which is
c                                 arbitrarily set at 1.d-8.
      if (dabs(ddv).lt.1d-8) then

         iugh = iugh + 1

         if (i.lt.3) then 
c                                 traverses 1,2.
            ddv = dv(ivd)
         else  
c                                 traverses 3,4.
            ddv = -dv(ivd)
         end if

         v(ivd) = vlast
         ddv = ddv/10.
         iflg1 = 0

         if (iugh.gt.3) then
            write (n3,1020) ivd,dv(ivd),(names(idv(j)),j = 1, icp)
            write (n3,*) ' '
            ier = 2
            goto 99
         end if 
      
         goto 100 

      end if 
c                                 increment conditions again:
      goto 100
c                                 next traverse.
60    continue

      ier = 1

1000  format (' Metastable assemblage in FLIPIT',/,

     *        ' the assemblage is:',4x,7(1x,a8))
1010  format (' v =',5(g12.6,1x))
1020  format (/,' **warning ver047** SEARCH: > 1 equilibrium',
     *         ' occurs within the minimum search',/,' increment for',
     *         ' variable: ',i1,'. This often occurs as YCO2 =>',
     *         ' 1 or => 0.',/,' You may be able to correct this by',
     *         ' assigning slightly different limits',/,' for YCO2 or',
     *         ' by reducing the default increment for this variable',
     *         ' (',g12.3,')',/,' in the computational option file.',/,
     *         ' Equilibria involving the following assemblage may',
     *         ' not be delineated:',/,7(1x,a8)) 
99    end

      subroutine sfol1 (ivd,ivi,ier,dv,ikwk,irend)
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k2=3000,k5=12,k7=k5+1,l2=5)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst5  /v(l2),tr,pr,r,ps/ cst80 /irchk(k2)/ cst7 /iflag
     *      / cst24 /ipot,jv(l2),iv(l2)/ cst9 /vmax(l2),vmin(l2),ddv(l2)
     *      / cst6  /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dg
     *      / cst31 /vn(k2,k7),irct,ird/ cst87 /delt(l2),dtol,utol,ptol     
c-----------------------------------------------------------------------
      if (ikwk.eq.1) then
         irchk(ird) = 1
         goto 9000
      end if

      iswtch = 0
      jswtch = 0 
c                                 begin traverse:
10    v(ivi) = v(ivi) + dv
c                                 is search in range?
      if (v(ivi).gt.vmax(ivi)) then
         v(ivi) = vmax(ivi)
      else if (v(ivi).lt.vmin(ivi)) then
         v(ivi) = vmin(ivi)
      end if
c                                 solve for the equilibrium conditions:
      call univeq (ivd,ier)
c                                 on error return:
c                                 calling routine will switch variables.
      goto (999,999),ier

      if (jswtch.eq.0.and.
     *   (v(ivd).gt.vmax(ivd).or.v(ivd).lt.vmin(ivd))) then 
         call switch (dv,ivi,ivd,jer)
         jswtch = 1
         if (jer.eq.0) goto 10
      end if 
c                                 test if the equilibrium is stable:
      call gall 
      call nschk 
c                                 iflag= 0 stable.
c                                 iflag=1 metastable w/respect a phase.
c                                 iflag=2 metastable w/respect to more
c                                         than one phase.
      goto (20,70),iflag
c                                 iflag= 0:
c                                 dependent v in range? if
c                                 greater than the maximum value for v
c                                 or less than the minimum value for v
c                                 reset conditions, and
c                                 switch independent/dependent variables
      if (v(ivd).gt.vmax(ivd)) then
            v(ivd) = vmax(ivd)
         else if (v(ivd).lt.vmin(ivd)) then
            v(ivd) = vmin(ivd)
         else
            call assptx
c                                 if independent variable has reached a
c                                 limit of the diagram output traverse
c                                 otherwise increment the variable and
c                                 continue the traverse
         if (v(ivi).eq.vmax(ivi).or.v(ivi).eq.vmin(ivi)) then
c                                 flag irchk indicates that the
c                                 traverse ended on an edge.
            irchk(ird) = 1

            call maxend
            goto 9000
         end if
c                                 try for next point
         goto 10
c                                 no route out of this construction
      end if
c                                 solve for the equilibrium with
c                                 the switched variables:
      call univeq (ivi,ier)
c                                 flag irchk indicates that the
c                                 traverse ended on an edge.
      irchk(ird) = 1
c                                 on error output traverse anyway
      goto (9000,9000),ier
c                                 assign final value
      call assptx
      goto 9000
c                                 iflag=2 reset to last stable condition
c                                 and redefine increment dv
70    call reptx

      dv = dv/2d0
      if (dabs(dv).gt.delt(ivi)) goto 10 
c                                 refined increment down to dtol, try
c                                 switching iv's, set iswitch flag
      if (iswtch.eq.0) then 
         call switch (dv,ivi,ivd,jer)
         iswtch = 1
         if (jer.eq.0) goto 10
      end if 

      call warn (73,delt(ivi),ivi,'SFOL2 ')
      iflag = 0 

      goto 9000
c                                 iflag=1 (an invariant equilirium):
c                                 iterate independent variable to refine
c                                 the conditions of the invariant point:
20    call findjp (ivi,ivd,dv,ip)                              
c                                 output the traversed equilibrium:
9000  ier = 0

      call outrxn (ip,ier)    
c                                 save conditions of the endpoint
      call svrend (ird,irend,jer)

999   end                         

      subroutine sfol2 (iovi,iovd,ip,irend)
c----------------------------------------------------------------------
c sfol2 generates and computes the loci of stable equilibria
c emanating from an ip indexed by the array ipid and the
c arguement (i).  since one of these equilibria is calculated
c by sfol1 or by previous calls to sfol2 only c+1 (at most)
c curves are calibrated.
c
c the phases identified by idv define a plane in g-x space, idphi
c is the phase which is tangent to this plane, and lphi is the
c ,or one of the, phases which is present at the ip but absent
c along the univariant point.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,l2=5,k5=12,k7=k5+1,k8=k5+2,k6=1000)

      character*8 names, text(400)*1, alpha(130)*1

      common/ cst8  /names(k1)/ cst6  /icomp,istct,iphct,icp
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1     
     *      / cst29 /vip(l2,k2),ipid(k2,k8),ipct
     *      / cst31 /vn(k2,k7),irct,ird/ cst105 /ibug(k2)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst24 /ipot,jv(l2),iv(l2)/ cst25 /vnu(k7),idr(k7),ivct
     *      / cst103 /isec,icopt,ifull,imsg,io3p/ cst32 /ptx(k6),ipt2
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst5 /v(l2),tr,pr,r,ps
c-----------------------------------------------------------------------
      jchk = 0 
      icp1 = icp + 1
      icp2 = icp + 2
      icp3 = icp + 3
      ivi = iovi
      ivd = iovd
      ivarip(ip) = 1

      goto (150), io3
      goto (150), io3p

      call iptext (text,iend,ip,icp2)

      write (n3,1040)
      write (n3,1020) ip, (text(i),i = 1, iend)
      write (n3,1050)

150   do 10 j = 2, icp2
         jtic = 0
         idno  = icp3 - j
         lphi  = ipid(ip,idno)
         idphi = ipid(ip,icp2)

         do 20 i = 1, icp1
            if (i.eq.idno) goto 20
            jtic = jtic + 1
            idv(jtic) = ipid(ip,i)
20       continue 
c                            balance the equilibrium, ier = 1
c                            if the plane idv is degenerate
         itic= 0
c                            set invariant point conditions
         v(ivd) = vip(ivd,ip)
         v(ivi) = vip(ivi,ip)
c                            check if assemblage sees bulk comp
         if (jbulk.ne.0) then
c                            put idphi into idv for bulkun
            call bulkun (idv,idphi,jbfg,icp)
            if (jbfg.eq.0) goto 10
         end if

110      call balanc (b,idv,idphi,ier)
         goto (80),ier
         goto (130),ivar
         ivarip(ip) = 2
         goto (10),isudo
c                            call pivots for lchk.
130      call pivots (ier)
         if (ier.eq.0) goto 90
c                            in some cases singularity of the
c                            transpose of the concentration matrix
c                            may not be detected by balanc, but,
c                            will be detected for the concentration
c                            matrix by pivots, in which case the
c                            columns are permuted and the stoichiometry
c                            of the reaction is redetermined.
80       if (itic.eq.icp-1) goto 100
         icter =idphi
c                              if a degenerate facet has been
c                              generated by permutation permute
c                              idphi with the vertices idv until
c                              balanc succeeds.  for a c component
c                              system no more than c-1 permutations
c                              will be necessary.
         idphi = idv(icp-itic)
         idv(icp-itic) = icter
         itic = itic + 1
         goto 110
100      goto 10
90       call assir 
c                                TESTING 
c                                test endpoint.
         call svrend (ird,irend,ier)

         goto (10),ier
c                                calculate the loci of the equilibia,
         ipt2 = 0
         call delrxn 
         call assptx
c                            determine direction of traverse and
c                            dependent and independent variables
         call wway (div,lphi,ivi,ivd,ikwk,ier)
c                            begin traverse
         if (ier.ne.0) then 
c                            call rxntext to save the equilibrium
c                            in the list in case it is not found 
c                            again:
            call rxntxt (iend,jend,text,alpha)
c                            set the bug flag:
            ibug(irct) = 1
            goto 10 
         end if 

         icter = 0
60       call sfol1 (ivd,ivi,ier,div,ikwk,irend)
         jchk = jchk + 1
         goto (70,70),ier
         ivd = iovd
         ivi = iovi
         goto 10
70       call switch (div,ivi,ivd,jer)    
         goto (75),jer
         icter = icter + 1
         if (icter.lt.4) goto 60
      
75       call warn (20,v(ivi),ivi,'SFOL2 ')

         call outrxn (ip,1)
         ivi =iovi
         ivd=iovd
10    continue

      if (jchk.eq.0) call warn (74,v(ivi),ivi,'SFOL2 ')
      if (io3p.eq.0) write (n3,1040)

1020  format (' equilibria about invariant point (',i4,'):',//,
     *        3x,200a1)
1040  format ('                       ------')
1050  format (/,' are listed below:',/)

      end

      subroutine sreset (jok,iste,vt,v,vti,vi)
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
c                                 if the assemblages conditions
c                                 of stability are already known (jok=5)
c                                 or are known to extend beyond the
c                                 current condition (jok>iste) quit:
      if (jok.gt.iste) goto 99
c                                 if the known conditions of stability
c                                 have been extened, i.e., either:
c                                 jok=iste  or v>vt
      if (jok.ne.iste) then
         jok=iste
      else if (iste.lt.3) then
         if (vt.gt.v) goto 99
      else
         if (vt.lt.v) goto 99
      end if

      vt = v
      vti = vi

99    end  

      subroutine maxend 
c----------------------------------------------------------------------
c save the most advanced condition for the stability of an 
c assemblage along the periphery of a diagram. 
c ismax - is the traverse
c value - is the dependent variable
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (l2=5)

      common/ cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / cst49 /value,ismax,igot
     *      / cst5  /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),ddv(l2)
c-----------------------------------------------------------------------
      if (igot.eq.0) igot = 1

      if (v(iv2).eq.vmin(iv2).and.ismax.eq.1) then
c                                  traverse 1
c                                  now don't allow backing up:
         if (value.lt.v(iv1)) value = v(iv1)

      else if (v(iv1).eq.vmax(iv1).and.ismax.lt.3) then
c                                  traverse 2
         if (ismax.eq.1) then

            ismax = 2
            value = v(iv2)

         else if (value.lt.v(iv2)) then

            value = v(iv2)

         end if

      else if (v(iv2).eq.vmax(iv2).and.ismax.lt.4) then
c                                  traverse 3
         if (ismax.lt.3) then

            ismax = 3
            value = v(iv1)

         else if (value.gt.v(iv1)) then

            value = v(iv1)

         end if 

      else if (v(iv1).eq.vmin(iv1)) then
c                                  traverse 2
         if (ismax.lt.4) then
            ismax = 4
            value = v(iv2)
         else if (value.gt.v(iv1)) then
            ismax = 4
            value = v(iv2)
         end if 

      end if 

      end 

      subroutine svrend (ird,irend,ier)
c----------------------------------------------------------------------
c save end points of univariant equilibria and 
c check if the intersection has already been recorded (ier =1).

c two cases:

c 1) reaction occurs on an edge, if the intersection has already 
c    occurred then ignore future occurrences.
c 2) reaction occurs at an IP only ignore occurences being traced in 
c    the same direction if the reaction is degenerate.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (l2=5,k2=3000,h5=7,h6=500)

      dimension rends(2,16,k2),irends(k2)

      common/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst35 /irv(k2)/ cst103 /isec,icopt,ifull,imsg,io3p
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst5  /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),ddv(l2)

      save rends, irends
c-----------------------------------------------------------------------
      ier = 0 
c                                    if isec = 1, don't do any checks
      goto (99) isec

      dx =  ddv(iv(1))
      dy =  ddv(iv(2))

      x = v(iv(1))
      y = v(iv(2))

      if (x.ne.vmin(iv(1)).and.x.ne.vmax(iv(1)) .and.
     *    y.ne.vmin(iv(2)).and.y.ne.vmax(iv(2))) then
c                                    the endpoint is not on an edge:
          if (isec.eq.2) then 
c                                    if isec = 2, don't check unless
c                                    curve intersects at an edge.
             goto 99
          else if (isec.eq.3) then 
c                                    if isec = 3, don't do a check for 
c                                    degenerate interior endpoints:
             if (irv(ird).lt.icp1) goto 99
          end if 
      else 
c                                    the endpoint is on an edge:
c                                    if isec = 4, don't do a check for
c                                    degenerate equilibria.
          if (isec.eq.4.and.irv(ird).lt.icp1) goto 99
      end if 
c                                    isec = 5, do check on all ends
      if (ird.le.irend) then 
c                                    reaction has been found already
c                                    look to see if its at the same 
c                                    endpoint (tolerance +/- dv).
         jct = irends(ird)
         do j = 1, jct 
            xx = rends(1,j,ird)
            yy = rends(2,j,ird)
            if (xx.gt.x-dx.and.xx.lt.x+dx.and.
     *          yy.gt.y-dy.and.yy.lt.y+dy) then 
c                                    condition match has occurred
               ier = 1
               goto 99
            end if 
         end do 
c                                    no condition match has occurred
         jct = jct + 1
         if (jct.gt.16) then 
            call warn (999,x,ier,'SVREND')
            jct = 16 
         end if  

      else 
c                                    reaction has not been recorded
         jct = 1
         irend = irend + 1

      end if 

      if (ird.gt.k2) call error (206,x,k2,'SVREND')
      irends(ird) = jct 
      rends(1,jct,ird) = x
      rends(2,jct,ird) = y

99    end 

      subroutine ufluid (fo2)
c----------------------------------------------------------------------
c subroutine ufluid computes the potential of the components
c of a saturated fluid phase. if the mole fraction of a component is les
c less than 1.d-38 the chemical potential is set to -9.9d09.
c ufluid may call one of three molecular fluid equations of state, or
c alternatively users may supply their own routines, however,
c the routines currently in use return the log of a components fugacity
c which is then added to the reference state potential computed by the
c function gphase.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k4=22,h5=7,k10=200)
 
      dimension xf(2)
 
      common/ cst1 /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst11 /f(2)
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst208 /ifct,idfl
c-----------------------------------------------------------------------
c                           compute the chemical potentials of
c                           fluid components in fluid saturated
c                           systems.
      call cfluid (fo2,fs2)

      if (idfl.ne.0) then
         uf(idfl) = gphase(idfl) + r * t * f(idfl)
      else
         xf(1) = 1d0 - xco2
         xf(2) = xco2
 
         do 10 i = 1, 2
            if (iff(i).eq.0) goto 10

            if (xf(i).gt.1d-38) goto 70
            uf(i) = -9.9d09
            goto 10

70          uf(i) = gphase(i) + r * t * f(i)
10       continue

      end if 
 
      end

      subroutine uproj
c----------------------------------------------------------------------
c subroutine uproj computes the potential of mobile
c or saturated components.  the potentials of
c saturated non-volatile components are projected through
c saturated volatile components.
c
c uproj is always called before the routines gphase or gproj, in
c therefore in order to increase program efficiency the values
c log(t),dsqrt(t),(p-pr)*dexp(-t/300d0), and dexp(-p/35000d0)
c are maintained in common block 'cst44'.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k4=22,k5=12,h5=7,h6=500,k10=200)
 
      double precision uss(h6)
 
      common/ cst1  /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5)
     *      / cst5  /p,t,xco2,u1,u2,tr,pr,r,ps
     *      / cst6  /icomp,istct,iphct,icp/ cst20 /comps(k1,k5)
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst44 /dlogt,dsqrtt/ cst2 /g(k1)
c-----------------------------------------------------------------------
c                           compute isobaric isothermal constants.
      if (t.lt.0d0.or.tr.eq.0d0) then
        dlogt = 0d0
        dsqrtt = 0d0
      else
        dlogt = dlog(t)
        dsqrtt = dsqrt(t)
      end if

      fo2 = 0d0

      goto (50),ifyn
c                           compute the chemical potentials of
c                           fluid components in fluid saturated
c                           systems.
      call ufluid (fo2)
c                           determine stable saturated composants
c                           and the corresponding chemical potentials
c                           (these are projected through the fluid
c                           components).
50    do 60 i = 1, isat

         ict = isct(i)

         ll = icp+i

         do j = 1, ict

            k = ids(i,j)

            uss(j) = gphase(k)-vf(1,k)*uf(1)-vf(2,k)*uf(2)
c                           if multiple component saturation constraints
c                           project by saturation hierarchy:

            goto (10),i

 

            i1 = i-1

            do l = 1, i1

               uss(j) = uss(j)-vs(l,k)*us(l)

            end do

10          g(k) = uss(j)

 

            uss(j) = uss(j)/comps(k,ll)

         end do 
c                           if O2, check if fo2 has been 
c                           determined by a fluid phase routine,
c                           if so do the projection:
         if (io2.eq.i) then 

            do j = 1, ict 

               uss(j) = uss(j) + r * t * fo2

            end do 

         end if 
c                           now find stable "composant":

         u = uss(1)

         id = 1

         goto (80),ict

         do 30 j = 2, ict

 

            if (uss(j).gt.u) goto 30

            id = j

            u = uss(j)

 

30       continue
c                               save the id of the stable composant.

80       idss(i) = ids(i,id)
c                               and its chemical potential.

         us(i) = u
c                               in case a phase in the component
c                               saturation space is an endmember of
c                               a solution correct the endmember
c                               properties:

         goto (60),i

         do j = 1, ict

            k = ids(i,j)

            g(k) = g(k) - vs(i,k)*u

         end do 

60    continue

99    end

      subroutine wway (div,lphi,ivi,ivd,ikwk,ier)
c-----------------------------------------------------------------------
c wway determines from which side of an invariant point a specific
c univariant equilibrium emanates. the appropriate sign for the
c increment of the independent variable (dv) is also determined.
c where necessary wway also refines dv and/or may switch the
c independent and independent intensive variables (ivi,ivd). vo1
c and vo2 are the values of the intensive variables iv1 and iv2 at
c the invariant point. lphi is the id of the absent phase.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,l2=5)

      character*8 names

      dimension vo(l2),vvi(2),vvd(2),ddv(2)

      common/ cst25 /vnu(k7),idr(k7),ivct/ cst8 /names(k1)
     *      / cst31 /vn(k2,k7),irct,ird/ cst6 /icomp,istct,iphct,icp
     *      / cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst5  /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(l2)
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst7 /iflag
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
c-----------------------------------------------------------------------
      iophi =idphi
      ikwk = 0
      icter = 0
      ier = 0
      vo(iv1) = v(iv1)
      vo(iv2) = v(iv2)
c                                 first try iv1 as the iv
      ivi = iv2
      ivd = iv1

      ifail= 0

5     call univeq (ivd,ier)

      if (ier.eq.0) then
c                                 is the equilibrium condition
c                                 for this reaction consistent
c                                 with the invariant point?
         if ( dabs( (v(ivd)-vo(ivd))/dv(ivd) ).gt.0.1d0) then 
            vo(iv1) = v(iv1)
            vo(iv2) = v(iv2)
            call warn (18,v(ivd),ird,'WWAY')
         end if 

      else
         goto (9800),ifail
         ivi = iv1
         ivd = iv2
         v(iv1) = vo(iv1)
         v(iv2) = vo(iv2)
         ifail = ifail + 1
         goto 5
      end if
c                                 test that the ip isn't on an edge
      if (v(iv1).eq.vmin(iv1).or.v(iv1).eq.vmax(iv1).or.
     *    v(iv2).eq.vmin(iv2).or.v(iv2).eq.vmax(iv2)) goto 9600

      jswit = 0

190   order = 1d0
c                                 set iv increment
60    div = dv(ivi) / order

      ifail = 0

      do 10 i = 1, 2
c                                 because of numerical slop 
c                                 an equilibrium may appear stable
c                                 on its metastable extension close to
c                                 an ip, because of this wway always
c                                 tests for stability on both extensions
c                                 if wway cannot solve for equilibrium
c                                 conditions on both sides of an ip
c                                 the increment is decreased and if this
c                                 fails the independent and dependent
c                                 variables are switched.
         jfail = 0

15       v(ivi) = vo(ivi) + div
c                                 check if the increment is too large
         if (v(ivi).gt.vmax(ivi)) v(ivi) = vmax(ivi)
         if (v(ivi).lt.vmin(ivi)) v(ivi) = vmin(ivi)
c                                 set the dependent variable to ip:
         v(ivd) = vo(ivd)
c                                 solve for the equilibrium
         call univeq (ivd,ier)

         if (jfail.le.5.and.ier.eq.0.and.(v(ivd).gt.vmax(ivd)
     *       .or.v(ivd).lt.vmin(ivd))) then
c                                 increment to big:
            div = div / 10.
            jfail = jfail + 1
            goto 15
  
         else if (jswit.eq.0.and.ier.gt.0) then 
c                                 on error switch dependent and
c                                 independent variables
            iv  = ivi
            ivi = ivd
            ivd = iv
            jswit = 1
            jfail = 0

            goto 190

         else if (jswit.gt.0.and.ier.gt.0) then

            call warn (999,dv(1),ivi,'WWAY')
            ier = 1
            return

         end if 

         vvi(i) = v(ivi)
         vvd(i) = v(ivd)
         ddv(i) = div
c                                 determine if the equilibrium is
c                                 metastable with respect to lphi
c                                 (lchk=1):
         if (lchk(lphi).eq.0) then
            ifail = ifail + 1
            isign = i
         end if
c                                 end of looop
10    div = -div
c                                 at this stage there are three
c                                 possibilities indicated by ifail
      goto (50), ifail
c                                 ifail= 0 or 2 the equilibrium was stable
c                                 or metastable on both sides, the most 
c                                 plausible scenario is the increment is
c                                 too big and overshoots the stability,
c                                 in this scenario ifail = 2; however the
c                                 ip could be located inside of, rather than
c                                 on the edge of the stability field, in 
c                                 which case there is a small chance the
c                                 increment is too small to cross the border.
c                                 
c                                 so first decrease the increment:
         if (order.ge.1d0) then 
            order = order * 10d0 
c                                 stop decreasing after four orders 
c                                 of magnitude, and try increasing
c                                 increment:
            if (order.gt.1d4) order = 0.1
         else 
            order = order / 10d0
c                                 stop increasing after two orders 
c                                 of magnitude, cause there should 
c                                 never be a default increment much less 
c                                 than 1 % of the variable range.
            if (order.lt.1d-2) goto 9500
         end if 
c                                  why the fuck was this here?
c         ivi = iv1
c         ivd = iv2

      goto 60
c                                 ifail=1, stable extension identified.
c                                 determine the sign of the increment
c                                 and its magnitude:
50    if (isign.eq.1) then
         div = dv(ivi)
         del= vmax(ivi)-vo(ivi)
         if (del.lt.div) div = del/3.d0
      else
         div = -dv(ivi)
         del= vo(ivi)-vmin(ivi)
         if (del.lt.-div) div = -del/3.d0
      end if
c                                 reset variable values and save
      v(ivi) = vvi(isign)
      v(ivd) = vvd(isign)
      odiv = ddv(isign)
      icter = 0
c                                 now check the stability with respect
c                                 to all phases
      order = 1d0
70    call gall 
      call schk (lphi)
c                                 stable if iflag= 0 from schk
c                                 now check if dependent variable is in
c                                 range.
      if (iflag.eq.0) goto 9999
c                                 metastable try refining the increment
      order = order * 5.d0
      div = div / order
      v(ivi) = vo(ivi) + div
      v(ivd) = vo(ivd)
100   call univeq (ivd,ier)

      if (ier.eq.0) goto 90
c                                 error from univeq, more than once quit
      goto (9400), icter
      div =odiv
      icter = icter + 1
      goto 100
c                                 no error from univeq
90    if (order.gt.10.d6) goto 9500
      goto 70
c                                 all systems go
9999  call assptx
      if (dabs(div).lt.dv(ivi)) then
         div = dv(ivi)*div/dabs(div)
      end if 

      if (v(ivi).eq.vmax(ivi).or.v(ivi).eq.vmin(ivi)) then
         ikwk= 1
         call maxend
      end if 

      return

9400  call warn (87,dv(1),ivi,'WWAY')
      ier = 1
      return  
       
9500  call warn (24,dv(ivi)/order,ivi,'WWAY')
      ier = 1
      return

9600  call warn  (63,dv(1),ivi,'WWAY')     
      ier = 1
      return
  
9800  call warn (58,dv(1),ier,'WWAY')
      if (ivct.gt.4) goto 9820
      write (n3,1050) ird,(vnu(l),names(idr(l)), l= 1, ivct)
      goto 9830
9820  write (n3,1050) ird,(vnu(l),names(idr(l)), l= 1, 4)
      write (n3,1060) (vnu(l),names(idr(l)), l= 5, ivct)
9830  write (n3,*) 
      ier = 1

1050  format (1x,'(',i3,')',4(1x,g9.3,1x,a8))
1060  format (6x,4(1x,g9.3,1x,a8),/,6x,4(1x,g9.3,1x,a8))
      end

      subroutine fultxt (nchar,rtxt)
c-----------------------------------------------------------------------
c subprogram to write text for complete reaction reactions
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,k7=k5+1,h5=7,h6=500)

      character*8 names, text(400,2)*1, rtxt(400)*1,
     *            ppart(k5+2)*30, rxnstr*400, char8, mpart(k5+2)*30,
     *            cname*5, exten(2)*8, ptext*400, mtext*400

      integer ichar(2)

      common/ csta4 /cname(k5)/ cst21  /du(2),dv(2),jds(h5),ifr 
     *      / cst88  /vu(2,k1),jprct,jmct
     *      / cst201 /vuf(2),vus(h5),iffr,isr/ cst104 /rxnstr(k2)
     *      / cst25  /vnu(k7),idr(k7),ivct/ cst8 /names(k1)
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst31 /vn(k2,k7),irct,ird/ cst80 /irchk(k2)
     *      / cst103 /isec,icopt,ifull,imsg,io3p

      save exten 
c                                     this is a bullshit trick and
c                                     will cause errors if someone
c                                     uses a function other than G.
      data exten/'-V(j/b) ','S(j/k) '/

      ip = 0 
      im = 0

      do i = 1, 8
         mpart(i) = '                              '
         ppart(i) = '                              '

      end do
c                             load id's of phases with + or - coeffs 
      do i = 1, ivct 
c           
         id = idr(i)   
         if (vnu(i).lt.0d0) then
            im = im + 1
            call wrpart (-vnu(i),ikp(id),names(id),mpart(im))
         else 
            ip = ip + 1
            call wrpart (vnu(i),ikp(id),names(id),ppart(ip))
         end if
      end do
c                             composant stoichiometry: 
      do i = 1, isat 
         if (vus(i).gt.0d0) then
            im = im + 1
            call wrpart (vus(i),0,names(jds(i)),mpart(im))
         else if (vus(i).lt.0) then
            ip = ip + 1
            call wrpart (-vus(i),0,names(jds(i)),ppart(ip))
         end if 
      end do 
c                                     fluid stoichiometry:
      if (ifyn.eq.0) then 
         do i = 1, 2
            if (iff(i).ne.0) then
               if (vuf(i).gt.0d0) then
                  im = im + 1
                  call wrpart (vuf(i),0,names(i),mpart(im))
               else if (vuf(i).lt.0d0) then
                  ip = ip + 1
                  call wrpart (-vuf(i),0,names(i),ppart(ip))
               end if 
            end if 
         end do
      end if 
c                                     mobile components: 
      do i = 1, jmct
         char8 = cname(jprct+i)
         if (du(i).gt.0d0) then
            im = im + 1
            call wrpart (du(i),0,char8,mpart(im))
         else if (du(i).lt.0d0) then
            ip = ip + 1
            call wrpart (-du(i),0,char8,ppart(ip))
         end if
      end do

      if (ifull.eq.1.or.ifull.eq.3) goto 55

      do i = 1, 2
         if (dv(i).gt.0d0) then
            im = im + 1
            call wrpart (dv(i),0,exten(i),mpart(im))
         else if (dv(i).lt.0d0) then
            ip = ip + 1
            call wrpart (-dv(i),0,exten(i),ppart(ip))
         end if
      end do 

55    if (ip.gt.8.or.im.gt.8) call error (997,du(1),ip,'FULTXT')

      write (ptext,1000) (ppart(i), i = 1, ip)
      write (mtext,1000) (mpart(i), i = 1, im)

      ichar(1) = im * 31
      read (mtext,1010) (text(i,1), i = 1, ichar(1))
      ichar(2) = ip * 31
      read (ptext,1010) (text(i,2), i = 1, ichar(2))
c                             filter out blanks:
      do i = 1, 2
         nchar = 1
         do 70 j = 2, ichar(i)
            if (text(j-1,i).eq.' '.and.text(j,i).eq.' ') goto 70
            nchar = nchar + 1
            text(nchar,i) = text(j,i)
70       continue
c                             filter out blanks before parentheses:
         ichar(i) = nchar
         nchar = 0
         do 80 j = 1, ichar(i)-1
            if ((text(j,i).eq.' '.and.text(j+1,i).eq.')').or.
     *          (text(j,i).eq.' '.and.text(j+1,i).eq.'(')) goto 80
            nchar = nchar + 1
            text(nchar,i) = text(j,i)
80       continue 
         ichar(i) = nchar
      end do
c                             concatenate parts into rxnstr:
      if (ichar(1)+ichar(2).gt.397) call error (208,du(1),ip,'FULTXT')
 
      nchar = 0 

      do 90 i = 1, 2
         do 100 j = 1, ichar(i)
            nchar = nchar + 1
100         rtxt(nchar) = text(j,i)
         if (i.eq.2) goto 90
         nchar = nchar + 3
         rtxt(nchar-2) = ' '
         rtxt(nchar-1) = '='
         rtxt(nchar)   = ' '
90    continue
c                             this is redundant if the reaction 
c                             was found before
      write (rxnstr(ird),1010) (rtxt(i),i=1,nchar)

1000  format (8(a30,1x))
1010  format (400a1)

      end 

      subroutine wrpart (vnu,ikp,name,part)
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (h9=30)

      character*8 name, fname*10, part*30

      common/ csta7 /fname(h9)
     *      / cst103 /isec,icopt,ifull,imsg,io3p

      if (ikp.eq.0) then
         if (ifull.gt.2) then  
            write (part,1000) vnu,name
         else 
            write (part,2000) name
         end if 
      else 
         if (ifull.gt.2) then 
            write (part,1010) vnu,fname(ikp),name
         else 
            write (part,2010) fname(ikp),name
         end if 
      end if 

1000  format (g9.3,1x,a8)
1010  format (g9.3,1x,a10,'(',a8,')')
2000  format (a8)
2010  format (a10,'(',a8,')')

      end

      subroutine bulkun (kdv,idum,jbfg,icp)
c----------------------------------------------------------------------
c check if any permutation of a c+1 phase assemblage sees a 
c specified bulk composition, if no jbfg = 1.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)

      integer kdv(k5),idv(k5+1)
c-----------------------------------------------------------------------
      icp1 = icp + 1

      do j = 1, icp
c                                      do c permutations of the
c                                      phases to assure that 
c                                      degeneracy is not a problem
         k = 1

         do 20 i = 1, icp1
            if (i.eq.j) then
               goto 20
            else if (i.lt.icp1) then
               idv(k) = kdv(i)
            else 
               idv(k) = idum
            end if 
            k = k + 1
20       continue 
c                                      check the assemblage
         call bulkck (idv,jbfg)
c                                      jbfg ne 0, bounds a comp.
         if (jbfg.ne.0) goto 99

      end do 
99    end 

      subroutine bulkck (jdv,jbfg)
c----------------------------------------------------------------------
c bulkck determines if the phases identified by jdv bound
c any of the specified bulk compositions.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,h5=7,k4=22,k0=25,k5=12,l2=5,h9=30,
     *           k12=14,k13=40000,k2=3000)
 
      character text(400)*1, cname*5, pname(k5)*20
      character vname*8, gname*78, xname*18

      integer jdv(k5+1),idsol(k5,2+k5),kdv(k5),idbulk(k5,k2)

      double precision fbulk(k5),vol(k5),fbulk1(k5),bsol(k5,k5),
     *                 vmode(k5),amode(k5),wtp(k5)
 
      common/ cst6   /icomp,istct,iphct,icp/ cst20 /comps(k1,k5)
     *      / cst10  /iff(2),idss(h5),ifug,ifyn,isyn/ csta4 /cname(k5)
     *      / cst5   /v(l2),tr,pr,r,ps/ cst301 /a(k5,k5),b(k5),ipvt(k5)
     *      / cst41  /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst43  /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind
     *      / cst300 /cblk(k5),atw(k1),jbulk/ cst24 /ipot,jv(l2),iv(l2)
     *      / csta2  /xname(k5),vname(l2),gname(l2)/ cst18 /jvol(2),ivol 
     *      / cst61  /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1) 
     *      / cst19  /jpoly/ cst103 /isec,icopt,ifull,imsg,io3p 

      save ibulk, idbulk

      data ibulk/0/ 
c----------------------------------------------------------------------
      jbfg = 0
c                                load the transpose of the
c                                concentration matrix of the pseudo-
c                                invariant assemblage.
      do i = 1, jbulk
         if (i.le.icp) then
            k = jdv(i)
            kdv(i) = k
         else
            k = idss(i-icp)
            kdv(i) = k
         end if
         do j = 1, jbulk
            a(j,i) = comps(k,j)
         end do
      end do        
c                                factor the matrix
      call factr1 (jbulk,ier)
      if (ier.eq.1) goto 99 
c                                test for bounded compositions
      do i = 1, jbulk
c                                load composition vector, the alpha
c                                vector is returned in the same array
         b(i) = cblk(i)
      end do
c                                solve for the alpha vector
      call subst1 (jbulk)
c                                a negative coefficient indicates
c                                the composition is not bounded
      do i = 1, icp
         if (b(i).lt.-1d-5) goto 99
      end do 
c                                the composition is bounded 
      jbfg = 1
      ier = 0
c                                the remainder of this routine is
c                                for output purposes:
      if (icopt.ne.5.and.jbulk.ge.icp) then
c                                do a check to prevent the data
c                                for the same assemblage being
c                                output more than once:
         do 40 i = 1, ibulk
            do 50 j = 1, jbulk
               do k = 1, jbulk
                  if (kdv(j).eq.idbulk(k,i)) goto 50
               end do
c                                no match on the i'th assemblage
               goto 40 
50          continue 
c                                if here a match on the i'th ass.
            goto 99
40       continue
c                                the assemblage is new
         ibulk = ibulk + 1
         if (ibulk.gt.k2) call error (183,b(1),k2,'BULKCK')
         do i = 1, jbulk
            idbulk(i,ibulk) = kdv(i)
         end do
      else 
         ibulk = ibulk + 1
      end if

      tot = 0d0
      atot = 0d0
      vtot = 0d0
      tot1 = 0d0
      atot1 = 0d0
      vtot1 = 0d0
      iflu = 0
      np = 0
      ncpd = 0

      jd = ibulk

      if (icopt.eq.5) jd = jpoly

      do i = 1, jbulk
         if (ikp(kdv(i)).ne.0) goto 10
      end do 

      ncpd = jbulk

      goto 20
c                                figure out how many solutions
c                                are present:
10    do 30 i = 1, jbulk
         if (ikp(kdv(i)).eq.0) then
            ncpd = ncpd + 1
            idsol(jbulk,1) = 0 
            idsol(jbulk,2) = ncpd
            idsol(jbulk,ncpd+2) = kdv(i)
            bsol(jbulk,ncpd) = b(i)
         else 
            do j = 1, np
               if (idsol(j,1).eq.ikp(kdv(i))) then
                  idsol(j,2) = idsol(j,2) + 1
                  idsol(j,idsol(j,2)+2) = kdv(i)
                  bsol(j,idsol(j,2)) = b(i)
                  goto 30
               end if 
            end do
c                                another solution must be present
            np = np + 1
            idsol(np,1) = ikp(kdv(i))
            idsol(np,2) = 1
            idsol(np,3) = kdv(i)
            bsol(np,1) = b(i)
         end if    
30    continue 
c                                now reform the arrays kdv and b
      itic = 0

      do i = 1, np
         do j = 1, idsol(i,2)
            itic = itic + 1
            b(itic) = bsol(i,j)
            kdv(itic) = idsol(i,2+j)
         end do
      end do

      do j = 1, ncpd
         itic = itic + 1
         b(itic) = bsol(jbulk,j)
         kdv(itic) = idsol(jbulk,2+j)
      end do 
c                                get text labels:
20    call blktxt (pname,text,iend,kdv,jbulk)

      do i = 1, jbulk
         k = kdv(i)
         vol(i) = vphase(k)
         vtot = vtot + b(i) * vol(i)
         atot = atot + b(i) * atw(k)
         tot = tot + b(i)         

         if (ifp(k).eq.0) then
            vtot1 = vtot1 + b(i) * vol(i)
            atot1 = atot1 + b(i) * atw(k)
            tot1 = tot1 + b(i)
         else
            iflu = 1
         end if 

      end do 
c                               normalize molar proportions:
      tot = tot / 1d2
      atot = atot / 1d2
      vtot = vtot / 1d2

      do i = 1, icomp
         fbulk(i) = 0d0
         fbulk1(i) = 0d0
      end do 

      write (n3,1000) jd,(text(i), i = 1, iend)
      write (n3,1010)
      write (n3,1120) (vname(jv(i)), v(jv(i)), i = 1, ipot)
      write (n3,1020)
c                              get true bulk
      do i = 1, jbulk
         k = kdv(i)
         vmode(i) = b(i)*vol(i)/vtot
         amode(i) = b(i)/tot
         write (n3,1030) pname(i), amode(i), b(i)*atw(k)/atot,
     *                             vmode(i), atw(k)/vol(i)

         do l = 1, icomp
            fbulk(l) = fbulk(l) + b(i) * comps(k,l)
            if (ifp(k).eq.0) fbulk1(l) = fbulk1(l) + b(i) * comps(k,l) 
         end do 
      end do

      if (jbfg.eq.2) then
         write (n3,1050)
         goto 99
      end if

      gtot = 0d0
      btot = 0d0
      gtot1 = 0d0
      btot1 = 0d0

      do i = 1, icomp
         btot = btot + fbulk(i)*atwt(i)
         gtot = gtot + fbulk(i)
         btot1 = btot1 + fbulk1(i)*atwt(i)
         gtot1 = gtot1 + fbulk1(i)
      end do 

      if (iflu.eq.0) then 

         write (n3,1040)

         do i = 1, icomp

            wtp(i) = fbulk(i)*atwt(i)/btot*1d2
            write (n3,1110) cname(i), fbulk(i)/gtot*1d2, wtp(i)                            
         end do

         rho = atot/vtot
         write (n3,1060) atot/vtot

      else 

         write (n3,1080)

         do i = 1, icomp

            wtp(i) = fbulk1(i)*atwt(i)/btot1*1d2
            write (n3,1090) cname(i),fbulk(i)/gtot*1d2,fbulk(i)*atwt(i)
     *                   /btot*1d2,fbulk1(i)/gtot1*1d2,fbulk(i)*atwt(i)
     *                   /btot1*1d2     
         end do

         rho = atot1/vtot1
         write (n3,1100) atot/vtot, atot1/vtot1

      end if 

      ivar = jbulk-np-ncpd

      write (n3,1070) ivar 

      if (io4.eq.0) then
c                                graphics output  
         write (n5,1130) jd,ivar,(kdv(i), i = 1, jbulk)
         write (n5,1140) (vmode(i), i = 1, jbulk)
         write (n5,1140) (amode(i), i = 1, jbulk)
         write (n5,1150) (v(jv(i)), i = 1, ipot)
         write (n5,1150) rho 
         if (ivol.gt.0) write (n5,1140) (wtp(jvol(i)),i=1,ivol)

      end if             

1000  format (/,9x,'---------',//,1x,'(',i5,') The bulk composition is',
     *          ' bounded by the assemblage:',//,2x,400a1)
1010  format (/,' in the thermodynamic composition space at:',/)
1020  format (/,25x,'mol %',6x,'wt %',5x,'vol %',3x,'rho (g/cm3)'/)
1030  format (1x,a20,3x,3(f6.2,4x),2x,f5.3)
1040  format (/,12x,' mol %',6x,'wt %',/)
1050  format (/,' WARNING - negative mass indicates the bulk'
     *         ,' composition is',/,' inconsistent with the'
     *         ,' saturation constraints',/)
1060  format (/,' Total density (g/cm3): ',f6.3)
1070  format (/,' Variance (c-p) = ',i2,//,9x,'---------',/)
1080  format (/,11x,'Complete Assemblage',10x,'Solids Only',
     *        /,12x,' mol %',6x,'wt %',9x,' mol %',6x,'wt %',/)
1090  format (1x,a8,3x,2(f6.2,4x),5x,2(f6.2,4x))
1100  format (/,' Total density (g/cm3)=',f6.3,2x,
     *          ' Solids only density (g/cm3)=',f6.3)
1110  format (1x,a8,3x,3(f6.2,4x),2x,f5.3)
1120  format (26x,a8,' = ',g12.6)
1130  format (14(i5,1x))
1140  format (12(f6.2,1x))
1150  format (5(g14.6,1x))

99    end

      subroutine ftext (text,ist,iend)
c-----------------------------------------------------------------------
c subprogram to filter blanks from text 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      character*1 text(400)

      itic = ist - 1

      do 40 i = ist, iend-1

         if (text(i).eq.' '.and.text(i + 1).eq.' ') then 
            goto 40 
         else if (text(i).eq.' '.and.text(i + 1).eq.')') then 
            goto 40
         else if (text(i).eq.' '.and.text(i + 1).eq.'(') then 
            goto 40 
         end if 

         itic = itic + 1
         text(itic) = text(i)

40    continue

      if (text(iend).ne.' ') then
         itic = itic + 1
         text(itic) = text(iend)
      end if 

      iend = itic + 1
      text(iend) = ' '
 
      end

      subroutine blktxt (pname,text,iend,kdv,jbulk)
c-----------------------------------------------------------------------
c subprogram to write a text labels for bulk composition output 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k1=5000,k2=3000,k5=12,h9=30)

      character*8 names, text(400)*1, fname*10, pname(k5)*20

      integer kdv(k5)

      common/ cst8  /names(k1)/ csta7 /fname(h9)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar    
c                             dump the names into the array text
      ist = 1

      do i = 1, jbulk
         id = kdv(i)
         if (ikp(id).eq.0) then
c                             simple compound:
            iend = ist + 7
            read (names(id),1000) (text(j), j = ist, iend)
            call ftext (text,ist,iend)
            write (pname(i),1000) (text(j), j = ist, iend-1)
            ist = iend + 1
         else  
c                             solution phases:
            iend = ist + 19
            read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9)
            text(ist + 10) = '('
            read (names(id),1000) (text(j), j = ist + 11, ist + 18)
            text(ist + 19) = ')'
            call ftext (text,ist,iend)
            write (pname(i),1000) (text(j), j = ist, iend-1)
            ist = iend + 1
         end if 
      end do 

1000  format (20a1)

      end 

      function vphase (id)
c-----------------------------------------------------------------------
c vphase returns the volume of the phase identified by kd. the volume
c is computed by finite difference from the Gibbs energy. the finite
c difference increment is one unit pressure, this should work assuming
c the unit is of the order of 1 bar. this routine is called only by
c bulkck for output purposes
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps
     *      / cst44 /dlogt,dsqrtt
c                                   set the variables in cst44
c                                   cause i don't know if gall 
c                                   has been called prior to bulkck
      if (t.lt.0d0.or.tr.eq.0d0) then
        dlogt = 0d0
        dsqrtt = 0d0
      else
        dlogt = dlog(t)
        dsqrtt = dsqrt(t)
      end if

      p0 = p
      g0 = gphase(id)
      p = p + 1d0
c                                   convert to cm3/mol assuming data
c                                   base units were J/bar
      vphase = 1d1*(gphase(id) - g0)
      p = p0
      
      end

      subroutine subst1 (n)
c-----------------------------------------------------------------------
c subst uses the lu decomposition of the matrix 'a' contained
c in the array 'a' to solve ax = b for x. subst is modified from the
c the subroutine of the same name listed by conte and de boor
c in 'elementary numerical analysis', mcgraw-hill, 1980.
c factor uses scaled partial pivoting.
 
c input     a- an n by n array containing the non-zero elements of
c              the u and l decompositions of a, as output by factor.
c           n- the dimension of the matrix a.
c        ipvt- a vector indicating that row ipvt(k) was used to
c              eliminate the coefficient a(n,k).
c           b- the vector b.
c output    b- the solution vector x.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)

      double precision x(k5)

      common / cst301 /a(k5,k5),b(k5),ipvt(k5)

c                            solve ly = b for y:
      ip = ipvt(1)
      x(1) = b(ip)
      do 10 i = 2, n
         sum = 0d0
         im1 = i-1
         do 20 j = 1, im1
20          sum = a(i,j)*x(j)+sum
         ip = ipvt(i)
10       x(i) = b(ip)-sum
c                            solve ux = y for x:
      x(n) = x(n)/a(n,n)
      nm1 = n-1
      do 30 ii = 1, nm1
         i = n-ii
         ip1 = i+1
         sum = 0.
         do 40 j = ip1, n
40          sum = a(i,j)*x(j)+sum
         x(i) = (x(i)-sum)/a(i,i)
30       b(i) = x(i)
      b(n) = x(n)
 
      end

      subroutine factr1 (n,ier)
c-----------------------------------------------------------------------
c factr1 is a subroutine which calculates the triangular
c decompositions of the matrix 'a'. factor is modified from
c the subroutine of the same name given by conte and de boor
c in 'elementary numerical analysis', mcgraw-hill, 1980.
c factor uses scaled partial pivoting.
c
c input     a- an n by n array containing the elements of matrix a.
c           n- the dimension of the matrix a.
c output    a- an n by n array containing the upper, u, and lower, l,
c              triangular decompositions of input matrix a.
c        ipvt- a vector indicating that row ipvt(k) was used to
c              eliminate the a(n,k).
c         ier- a flag, zero if a is of rank = n, and 1 if a is of
c              lower rank.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)

      common/ cst301 /a(k5,k5),d(k5),ipvt(k5)
c-----------------------------------------------------------------------
      ier = 0
c                            initialize ipvt,d
      do 10 i = 1, n
      ipvt(i) = i
      rmax = 0.
      do 20 j = 1,n
20    rmax = dmax1(rmax,dabs(a(i,j)))
c                            ax = b is singular if rmax = 0
      if (dabs(rmax).lt.1d-5) goto 9000
10    d(i) = rmax
c                            begin decomposition:
      nm1 = n-1
      do 30 i = 1,nm1
      ip1 = i+1
c                            determine pivot row (istr).
      rmax = dabs(a(i,i))/d(i)
      istr = i
      do 40 j = ip1,n
      tmax = dabs(a(j,i))/d(j)
      if (tmax.le.rmax) goto 40
      rmax = tmax
      istr = j
40    continue
      if (dabs(rmax).lt.1d-5) goto 9000
c                            if istr gt i, make i the pivot row
c                            by interchanging it with row istr.
      if (istr.le.i) goto 50
      j = ipvt(istr)
      ipvt(istr) = ipvt(i)
      ipvt(i) = j
      temp = d(istr)
      d(istr) = d(i)
      d(i) = temp
      do 60 j = 1,n
      temp = a(istr,j)
      a(istr,j) = a(i,j)
60    a(i,j) = temp
c                            eliminate x(k) from rows k+1,...,n.
50    do 70 j = ip1,n
      a(j,i) = a(j,i)/a(i,i)
      ratio = a(j,i)
      do 80 k = ip1,n
80       a(j,k) = a(j,k)-ratio*a(i,k)
70    continue
30    continue
      if (dabs(a(n,n)).lt.1d-5) ier = 1
      return
c                           algoritmic singularity.
9000  ier = 1
 
      end

      subroutine money
c
      open (97, file = 'money', iostat = ierr, status = 'old')
      if (ierr.ne.0) then
c                                 system could not find the file
         write (*,1000) 
         stop
1000     format (/,' There is no money file, make one and',
     *             ' deposit at least a dime.',/)
 
      end if

      read (97,*) imoney

      rewind (97)

      if (imoney.le.0) then

         write (*,1010) 
         stop
1010     format (/,' You are out of money, deposit at least a dime'/)

      else 

         imoney = imoney - 1
         write (97,*) imoney

      end if

      close (97)

      end

      subroutine subst (a,ipvt,n,b)
c-----------------------------------------------------------------------
c subst uses the lu decomposition of the matrix 'a' contained
c in the array 'a' to solve ax = b for x. subst is modified from the
c the subroutine of the same name listed by conte and de boor
c in 'elementary numerical analysis', mcgraw-hill, 1980.
c factor uses scaled partial pivoting.
 
c input     a- an n by n array containing the non-zero elements of
c              the u and l decompositions of a, as output by factor.
c           n- the dimension of the matrix a.
c        ipvt- a vector indicating that row ipvt(k) was used to
c              eliminate the coefficient a(n,k).
c           b- the vector b.
c output    b- the solution vector x.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k5=12)
 
      dimension a(k5,k5),ipvt(k5),b(k5+1),x(k5)
c                            solve ly = b for y:
      ip = ipvt(1)
      x(1) = b(ip)
      do 10 i = 2,n
         sum = 0d0
         im1 = i-1
         do 20 j = 1, im1
20          sum = a(i,j)*x(j)+sum
         ip = ipvt(i)
10       x(i) = b(ip)-sum
c                            solve ux = y for x:
      x(n) = x(n)/a(n,n)
      nm1 = n-1
      do 30 ii = 1, nm1
         i = n-ii
         ip1 = i+1
         sum = 0d0
         do 40 j = ip1, n
40          sum = a(i,j)*x(j)+sum
         x(i) = (x(i)-sum)/a(i,i)
30       b(i) = x(i)
      b(n) = x(n)
 
      end

      subroutine testit (ier)
c-----------------------------------------------------------------------
c ier = 0 the plane is stable
c ier = 1 the plane is metastable
c ier = 2 the assemblage is degenerate
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k5=12,l2=5,k7=k5+1)
 
      character*8 names
 
      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),iophi,idphi,
     *               iiphi,iflg1
     *      / cst8  /names(k1)/ cst2 /g(k1)
     *      / cst6  /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / cst52 /h,id(k7),ic(k5)/ cst87 /delt(l2),dtol,utol,ptol
c                              optional routine for testing
c                              for metastable g-x planes.
      ier = 0

      call abload (*991)
 
      do 2020 j = istct, iphct
 
         dg = g(j)
 
         do 10 i = 1, icp
10          dg = dg - cp(i,j) * b(i)
 
         if (dg.gt.dtol) goto 2020
c                                check that a phase is not metastable
c                                with respect to itself, this
c                                could be remedied by changing the 
c                                value of dtol.
         do 25 k = 1, icp
25          if (id(k).eq.j) goto 2020
c                                the assemblage is metastable with
c                                respect to phase j
         ier = 1
 
2020  continue

      goto 99

991   ier = 2
 
99    end

      subroutine outtit (n)
c-----------------------------------------------------------------------
c outtit writes title information and a brief description of the
c chemical system for each calculation requested.
 
c the subroutine arguments are:
c   'm'- specifies the logical unit for the output.
c   'h'- indexes the calculation.
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k5=12,h5=7,h6=500,k0=25,l2=5,
     *           h8=100,j9=40000)
 
      character*40 cname*5,cmpnt*5,names*8,exname*8,
     *             dname,vname*8,gname*78,xname*18,title*72
 
      common/ cst208 /ifct,idfl/ cst36 /exname(h8)
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *      / cst3  /x(k1,k5)/ cst20 /comps(k1,k5)/ cst8 /names(k1)
     *      / cst6  /icomp,istct,iphct,icp/ csta8 /title
     *      / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(k5)
     *      / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9
     *      / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5
     *      / csta2 /xname(k5),vname(l2),gname(l2)/ cst113 /ifproj
     *      / cst27 /jasmbl(j9),iasmbl(j9)
c----------------------------------------------------------------------
      if (n.eq.n3) goto (999),io3
      if (n.eq.n4) goto (999),io4
 
      write (n,1000)
c                          title:
      write (n,1190) title
c                          data base
      write (n,1210) dname(idbase)
c                          fluid
      if (ifyn.eq.0.or.ifproj.ne.0) call rfluid (2,ifug,n) 
c                          independent potentials:
      write (n,1070) (vname(jv(j)), j = 1, ipot)
c                          saturated phase components:
      goto (351), ifyn
      j = icp + isat
      write (n,1200) (cname(j+i), i = 1, ifct)
c                          saturated components:
351   goto (360),isyn
      j = icp + isat
      write (n,1180) (cname(i), i = icp1, j)
c                          unconstrained components
360   write (n,1080) (cname(i), i = 1, icp)
c                          phases
      goto (390,380,370),icp
      write (n,1150) (cname(i), i = 1, icp)
      do 400 i = istct, iphct
         write (n,1160) names(i), (x(i,j), j = 1, icp)
400   continue
      goto 410
370   write (n,1090) (cname(i), i = 2, 3)
      write (n,1100) (names(i),x(i,2),x(i,3), i = istct, iphct)
      goto 410
380   write (n,1040) cname(2)
      write (n,1030) (names(i), x(i,2), i = istct, iphct)
      goto 410
390   write (n,1130)
      write (n,1110) (names(i), i = istct, iphct)
c                          saturation composant phases,
c                          load the indexes into iasmbl:
410   if (isat.eq.0) goto 450
      itot = 0
      do 420 i = 1, isat
         ic11 = isct(i)
c                          do a check that there is at least
c                          one composant for each saturated component:
c                          why here? well, why not?
         if (ic11.eq.0) call error (279,x(1,1),i,'OUTTIT')
 
         do 430 j = 1, ic11
430         iasmbl(itot+j) = ids(i,j)
420   itot = itot+ic11
 
      write (n,1170)
      write (n,1110) (names(iasmbl(i)), i = 1, itot)
c                          excluded phases
450   goto (900),iexyn
      write (n,1140)
      write (n,1110) (exname(i), i = 1, ixct)
 
900   write (n,1000)
 
1000  format (/,80('-'),/)
1010  format (' Summary for the ',a3,' calculation follows:')
1030  format (4(2x,a8,1x,f5.3))
1040  format (/,1x,'Phases and (projected) mol fraction ',a5,' :',/)
1070  format (/,' Independently constrained potentials:',//,3x,8(a8,1x))
1080  format (/,' Components with unconstrained potent'
     *       ,'ials:',//,3x,10(a5,3x))
1090  format (/,' Phases and (projected) composition with respect to '
     *         ,a5,' and ',a5,':',/)
1100  format (3(1x,a8,2x,f5.3,2x,f5.3,5x))
1110  format (7(1x,a8,1x))
1130  format (/,' Phases:',/)
1140  format (/,' Excluded phases:',/)
1150  format (/,' Phases and (projected) compositions:',//,
     *        11x,12(1x,a5,1x),/)
1160  format (3x,a8,12(1x,f5.3,1x))
1170  format (/,' Phases on saturation and buffering surfaces:',/)
1180  format (/,' Saturated or buffered components:',//,3x,7(a5,3x))
1190  format (/,' Problem title: ',a56,/)
1200  format (/,' Saturated phase components:',//,3x,5(a5,3x))
1210  format (' Thermodynamic data base from: ',a40)
999   end

      subroutine tertex (icase)
c--------------------------------------------------------------------
c allows a user to input the initial assemblage for a fixed
c bulk composition calculation. 

c returns 2 cases:

c icase = 0 -> do a normal calculation.
c icase = 1 -> the user has input a valid initial assemblage

c   referenced by: combin
c   references to: testit, assdc
c---------------------------------------------------------------------
      implicit double precision (a-g,o-y),integer (h-n,z)

      parameter (k5=12,j9=40000,k1=5000)

      character*1 y, names*8

      common/ cst52 /hcp,idv(k5+1),ic(k5)/ cst96 /idcf(k5,j9),icfct
     *      / cst6  /icomp,istct,iphct,icp/ cst8  /names(k1)
c-----------------------------------------------------------------------

      write (*,1020) 
      read (*,1000) y  

      icase = 0

      if (y.eq.'y') then

         icase = 1

         hcp = icp

         do i = 1, icp
            ic(i) = i 
         end do 

         write (*,1030) 
         read (*,1000) y
         if (y.eq.'y') then
            do 9 i = 1, iphct, 5
               k = i + 5
               if (k.gt.iphct) k = iphct
9              write (*,1010) (j, names(j), j = i, k)
         end if
         write (*,1040) 
         read (*,*) (idv(i), i = 1, icp)
c                                             check stability
         call testit (ier) 
         if (ier.eq.1) then
            call error  (107,r,ier,'TERTEX')
         else if (ier.eq.2) then
            call error (108,r,ier,'TERTEX')
         end if 
c                                             assign assemblage
c                                             and check composition
         call assdc (ier)
         if (ier.ne.2) call error  (109,r,ier,'TERTEX')

         write (*,1060) (idv(i), i = 1, icp)
         write (*,1070)

      end if    

1000  format (a1)
1010  format (1x,6(i4,1x,a8,2x))
1020  format (/,' specify the initital assemblage (y/n)?')
1030  format (/,' list current phase indices (y/n)?')
1040  format (/,' enter the indices of the initial stable phases:')
1060  format (/,' the indices of the stable phases are:',//,10(2x,i4))
1070  format (/)
      end

      subroutine findas

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

      parameter (k1=5000,k5=12)

      character*5 cname

      integer kmax(k1)

      common/ cst52  /hcp,id(k5+1),ic(k5)/ cst2 /g(k1)
     *      / cst23  /a(k5,k5),u(k5+1),ipvt(k5),idv(k5),iophi,idphi,
     *                iiphi,iflg1
     *      / cst6   /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)
     *      / csta4  /cname(k5)

      do i = istct, iphct
c                                   sort phases into subcompositions
         do j = 1, icp
            if (cp(j,i).gt.1d-5) kmax(i) = j
         end do
      end do

      do hcp = 1, icp

         ic(hcp) = hcp
c                             make the first guess for the
c                             stable hcp component phase:
         do i = istct, iphct
            if (kmax(i).eq.hcp) then
               id(hcp) = i
               goto 20
            end if
         end do
c                             missing composant error
         call error (15,r,i,cname(hcp))
c
20       call abload (*992)
c                             start test loop:
         do 10 j = istct, iphct 
            if (kmax(j).gt.hcp) goto 10
            dg = g(j)
            do i = 1, hcp
               if (id(i).eq.j) goto 10
               dg = dg - cp(i,j)*u(i)
            end do
            if (dg.lt.-1d-8) then
c                             plane is metastable
c                             replace the most recently
c                             added test phase:
               id(hcp) = j
               call abload (*992)
            end if
10       continue

      end do
c                                   this is to reset the component
c                                   counter in ABLOAD via cst52. 
      hcp = icp

      goto 99

992   call error (999,dg,i,'FINDAS')

99    end 
         

      subroutine checkd (jd)

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

      parameter (k5=12)

      common / cst52 /hcp,id(k5+1),ic(k5)/ cst6 /icomp,istct,iphct,icp
c                           istabl = 1 if the original
c                           assemblage was metastable, istable = 0
c                           otherwise.
      do i = 1, icp 
         if (id(i).eq.jd) goto 99
      end do

      if (dgphc(jd).gt.-1d-05) goto 99

      idold = id(icp)
      id(icp) = jd 

      call abload (*110)

      goto 99
c                           restore old vertex
110   id(icp) = idold

      call abload (*9000)
      
      goto 99

9000  call error (19,1d0,jd,'CHECKD')

99    end 

      function dgphc (ld)
c-----------------------------------------------------------------------
c dgphc is the g difference between a phase
c identified by 'ld', and the plane defined by 'u'.
c   referenced by: 
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k5=12)
 
      common/ cst23 /a(k5*k5),u(k5+1),is(2*k5+4)/ cst2 /g(k1)
     *      / cst6  /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)
 
      dg = g(ld)

      do i = 1, icp 
         dg = dg - cp(i,ld)*u(i)
      end do 

      dgphc = dg

      end

      subroutine miscib (jdv,jvar,imyn)
c----------------------------------------------------------------------
c miscib determines the variance [ivar] of the assemblage of c phases
c identified by the array jdv. miscib also checks for immiscibility 
c in solution phases.
c----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (k1=5000,k2=3000,k5=12)
 
      integer iok(2),jdv(k5+1)
 
      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),iophi,idphi,
     *               iiphi,iflg1      
     *      / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar
     *      / cst6  /icomp,istct,iphct,icp/ cst12 /cp(k5,k1)
c----------------------------------------------------------------------
      
c                                initialize flags and counters
      iok(1) = 0
      ihom = 1
      iok(2) = 0
      im = 0

      jvar = 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, icp
         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.gt.icp) 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
      jvar = 2
c                                load and factor the transpose of the
c                                concentration matrix of the pseudo-
c                                invariant assemblage.
      do 70 i = 1, icp
         do 80 j = 1, icp
80          a(j,i) = cp(j,jdv(i))
70    continue
 
      call factor (a,icp,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 (ikp(j).ne.ikt) goto 60

         do 90 i = 1, icp
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,icp,b)
c                                a negative coefficient indicates
c                                the phase is not bounded
      do 100 i = 1, icp
100      if (b(i).lt.1d-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.icp) goto 110
      do 120 i = 1, icp
120      if (ikp(jdv(i)).ne.ikt) goto 60
c
110   jvar = 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.icp) jvar = 3
 
99    end 

      subroutine findjp (ivi,ivd,dv,ip)
c-----------------------------------------------------------------------
c findjp - called by sfol1 to iteratitvely refine the location of 
c an invariant point, the invariant point is located within an 
c energy tolerence specified by the variable ptol in routine pchk
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (l2=5)

      common/ cst5  /v(l2),tr,pr,r,ps/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst63 /delv(l2)/ cst6  /icomp,istct,iphct,icp/ cst7 /iflag
     *      / cst87 /delt(l2),dtol,utol,ptol     
c                                 iterate independent variable to refine
c                                 the conditions of the invariant point:
      iwarn = 0 
      ier = 0 
      ip = 0
c                                 step back to stable condition
20    call reptx
      dv = dv / 2d0
c                                 check if the point is already known
      call sameip (ip)

      if (ip.ne.0) goto 99

      if (dabs(dv).lt.delt(ivi)) goto 30
 
10    v(ivi) = v(ivi) + dv 

      call univeq (ivd,ier)
      if (ier.ne.0) goto 30 

      call pchk (dg,igo)

      if (igo.eq.1) then
         goto 40 
      else if (iflag.eq.1) then
         goto 20
      else  
         goto 10
      end if 
c                                 univeq blew up, assign ip with 
c                                 warning:
30    call reptx 
      iwarn = 1 
c                                 assign the invariant assemblage:
40    call assip (ip)
      
      if (iwarn.eq.1) call warn (47,ptol,ip,'FINDJP')

99    end 

      subroutine ssaptx
c---------------------------------------------------------------------
c ssaptx assigns the current value of v1 and v2 to the array ptx and
c increments the ordinate counter ipt2 if the value of either variable
c has increased by more than 1% of the default increment. 
c---------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)
 
      parameter (l2=5,k6=1000)
 
      common/ cst32 /ptx(k6),ipt2/ cst5  /v(l2),tr,pr,r,ps
     *      / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst63 /delv(l2)
     *      / cst9  /vmax(l2),vmin(l2),dv(l2)

      if (dabs((v(iv1)-ptx(ipt2-1))/dv(iv1)).gt.1d-2.or.
     *    dabs((v(iv2)-ptx(ipt2))/dv(iv2)).gt.1d-2) then 
  
         ipt2 = ipt2 + 2

         if (ipt2.gt.k6) ipt2 = k6

         ptx(ipt2-1) = v(iv1)
         ptx(ipt2)   = v(iv2)

      else 

         ptx(ipt2-1) = v(iv1)
         ptx(ipt2)   = v(iv2)

      end if 
 
      end


      subroutine assip (ip)
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k2=3000,k5=12,l2=5)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst5  /v(l2),tr,pr,r,ps/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst29 /vip(l2,k2),ipid(k2,k5+2),ipct
     *      / cst6  /icomp,istct,iphct,icp
c-----------------------------------------------------------------------
      idv(icp+1) = iophi
      idv(icp+2) = idphi
c                                  assign the new elements of ipid
40    ipct = ipct + 1
      if (ipct.gt.k2) call error (182,r,k2,'ASSIP')

      do j = 1, icp + 2
         ipid(ipct,j) =idv(j)
      end do 

      vip(1,ipct) = v(1)
      vip(2,ipct) = v(2)
      vip(3,ipct) = v(3)
      vip(4,ipct) = v(4)
      vip(5,ipct) = v(5)

      ip = ipct

99    end

      subroutine sameip (ip)
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-z),integer (h-n)

      parameter (k2=3000,k5=12,k8=k5+2,l2=5)

      common/ cst23 /a(k5,k5),b(k5+1),ipvt(k5),idv(k5),
     *               iophi,idphi,iiphi,iflg1
     *      / cst5 /v(l2),tr,pr,r,ps/ cst24 /ipot,jv(l2),iv(l2)
     *      / cst7 /iflag/ cst29 /vip(l2,k2),ipid(k2,k8),ipct
     *      / cst9 /vmax(l2),vmin(l2),ddv(l2)
     *      / cst6 /icomp,istct,iphct,icp
c-----------------------------------------------------------------------
      ip = 0
c
      if (ipct.eq.0) goto 99
 
      icp2 = icp + 2
      idv(icp+1) = iophi
      idv(icp2) = idphi
c                                 determine if a new ip was found
      jps = 1
5     ips = jps 

      do 10 i = ips, ipct

         do 20 j = 1, icp2
            do 30 k= 1, icp2
30             if (idv(k).eq.ipid(i,j)) goto 20
c                                 check if invariant point is equivalent
c                                 to the ith ip already found
               goto 10
20       continue
c                      the new ip assemblage is equivalent to one found 
c                      earlier, check if the conditions match:
         do 25 j = 1, 5
            if (dabs(vip(j,i)-v(j)).gt.ddv(j).and.
     *                                 ddv(j).ne.0d0) then
               if (i.lt.ipct) then 
                  jps = i + 1
                  goto 5
               end if 
               goto 99 
            end if 
25       continue 
c                                 conditions and assemblage match
         v(iv(1)) = vip(iv(1),i)
         v(iv(2)) = vip(iv(2),i)
         call assptx
c                                 reset iflag so sfol2 doesn't try 
c                                 anything stupid
         iflag = 0
         ip = i
         goto 99

10    continue

99    end
