	subroutine cieqsj( ns, y, h, p, lu, iwk, rwk, cwk,
     1	                   yeq, teq )

!  Routine to determine chemical equilibrium from stanjan

!  input:
!	ns  - number of species
!	y   - initial species mass fractions
!	h   - enthalpy (Chemkin units)
!	p   - pressure (Chemkin units)
!	lu  - unit number for output
!	iwk - Chemkin integer work array
!	rwk - Chemkin real work array
!	cwk - Chemkin character work array
!  output
!	yeq - species mass fractions at equilibrium
!	teq - temperature at equilibrium

	implicit double precision (a-h,o-z), integer (i-n) 
	parameter ( nsmax = 100, lenwrk = 10000 )
	dimension y(ns), yeq(ns)

	integer iwk(*), IEQWRK(lenwrk)
	dimension rwk(*), x(nsmax), KCON(nsmax), XCON(nsmax),
     1	          REQWRK(lenwrk)
	logical LPRNT, EQST,  LCNTUE
	character*16  cwk(*), ATOM(nsmax), KSYM(nsmax)

!-------  check that species are non-negative

	do i = 1, ns
	   if( y(i) .lt. 0.d0 ) then
	      write(0,*)'cieqsj: stopping - called with negative species'
	      write(lu,*)'cieqsj: stopping - called with negative species'
	      write(lu,601)( y(j), j=1,ns )
601	      format((1p,5e13.4))
	      stop
	   endif
	end do

	call ckrp( iwk, rwk, ru, ruc, patm )
	call ckindx( iwk, rwk, ne, kk, ii, nfit )


	LOUT   = lu
	LPRNT  = .false.
	LSAVE  = 0
	EQST   = .false.
	LCNTUE = .false.
	NCON   = 0
	call eqlen( ne, ns, NCON, LENIEQ, LENREQ )

	if( max(LENIEQ,LENREQ) .gt. lenwrk ) then
	   write(0,*)'cieqsj: increase lenwrk to at least ',
     1	             max(LENIEQ,LENREQ)
     	   stop
	endif

	do i = 1, lenwrk
	   REQWRK(i) = 0.
	end do

	call cksyme( cwk, lu, ATOM, kerr )
	call cksyms( cwk, lu, KSYM, kerr )
	NOP  = 5
	KMON = 0
C***	KMON = 2
	call ckytx( y, iwk, rwk, x )
	call temphy( ns, h, y, iwk, rwk, TIN )
	TEST = TIN
	PIN  = p/patm
	PEST = PIN

	call ckytx( y, iwk, rwk, x )
	call EQUIL (LOUT, LPRNT, LSAVE, EQST, LCNTUE, iwk,
     1           rwk, LENIEQ, IEQWRK, LENREQ, REQWRK, ne, ns,
     2           ATOM, KSYM, NOP, KMON, x, TIN, TEST, PIN, PEST,
     3           NCON, KCON, XCON)

	call EQSOL (ns, REQWRK, x, yeq, teq, PDUM,
     1           HDUM, VDUM, SDUM, WMDUM, CDUM, CDET)

	return
	end
