!  ch4eq.f		S.B. Pope  February 10, 2001

!  Fortran program to evaluate (and print out) the composition of CH4-air 
!  mixtures (as functions of mixture fraction) assuming (a) complete combustion
!  and (b) chemical equilibrium.
!  The output file ch4cc.out consists of nxi rows with columns:
!    xi, Y_H2, Y_O2, Y_N2, Y_H2O, CH4, CO, CO2, T, rho
!  The output file ch4eq.out consists of nxi rows and ns+3 columns
!    xi, (Y_k,k=1,ns), T, rho

	implicit double precision (a-h,o-z), integer (i-n)
	parameter ( nsmax=100, nxi=101, lck=10000 )
	dimension ickwrk(lck), rckwrk(lck)
	dimension y0(nsmax), y1(nsmax), yub(nsmax), yfb(nsmax)
	dimension x0(nsmax), x1(nsmax), xub(nsmax), xfb(nsmax)
	dimension yeq(nsmax)
	character*16 cckwrk(lck), kname(nsmax)
	logical kerr

	data xim / 0.053 /
	data pratm/ 1.0 /          !  pressure in atmospheres

!  initialize Chemkin

	open( 10, file='chem.bin', err=100, form='unformatted' )
	rewind( 10 )
	call ckinit( lck, lck, lck, 10, 6,  
     1	             ickwrk, rckwrk, cckwrk )

!  set pressure

	call ckrp( ickwrk, rckwrk, ru, ruc, pa ) 
	press = pratm * pa

!  obtain indexes of species

	call ckindx( ickwrk, rckwrk, mm, ns, ii, nfit )
	call cksyms( cckwrk, 6, kname, kerr )
	if( kerr ) then
	   write(0,*)'ch4eq: error in call to cksyms'
	   stop
	endif

	call ckcomp( 'H2',  kname, ns, kh2 )
	call ckcomp( 'O2',  kname, ns, ko2 )
	call ckcomp( 'N2',  kname, ns, kn2 )
	call ckcomp( 'H2O', kname, ns, kh2o )
	call ckcomp( 'CH4', kname, ns, kch4 )
	call ckcomp( 'CO ', kname, ns, kco  )
	call ckcomp( 'CO2', kname, ns, kco2 )

	if( min( kh2, ko2, kn2, kh2o, kch4, kco, kco2 ) .lt. 1 ) then
	   write(0,*)'h2cc: error obtaining indexes'
	   write(0,*)'h2cc: kh2, ko2, kn2, kh2o, kch4, kco, kco2 = ',
     1	              kh2, ko2, kn2, kh2o, kch4, kco, kco2
	   stop
	endif
	write(0,*)kh2, ko2, kn2, kh2o, kch4, kco, kco2

!  specify air and fuel compositions

	do k = 1, ns
	   x0(k) = 0.
	   x1(k) = 0.
	end do

	x0(ko2) = 0.21
	x0(kn2) = 0.79
	x1(kch4)= 1.
	t0      = 300.
	t1      = 300.

	call ckxty( x0, ickwrk, rckwrk, y0 )
	call ckxty( x1, ickwrk, rckwrk, y1 )

!  determine enthalpies

	call ckhbms( t0, y0, ickwrk, rckwrk, h0 )
	call ckhbms( t1, y1, ickwrk, rckwrk, h1 )

!  open output file

	lucc = 12
	lueq = 13
	open( lucc, file='ch4cc.out' )
	open( lueq, file='ch4eq.out' )

!  loop over mixture fraction

	power = log(xim)/log(0.5)

	do ixi = 1, nxi
	   xx  = (ixi-1)/float(nxi-1)
	   xi  = xx**power

!  determine unburnt composition

	   do k = 1, ns
	      yub(k) = (1.-xi) * y0(k) + xi * y1(k)
	   end do

	   hub = (1.-xi) * h0 + xi * h1

	   call ckytx( yub, ickwrk, rckwrk, xub )
	   call ckytx( yub, ickwrk, rckwrk, xfb )

!  O2 first used in:  CH4 + 1/2 O2 = CO + 2H2

	   if( xfb(kch4) .gt. 2. * xfb(ko2) ) then
!  O2 is deficient
	      xfb(kch4) = xfb(kch4) - 2. * xfb(ko2) 
	      xfb(kco)  = xfb(kco)  + 2. * xfb(ko2) 
	      xfb(kh2)  = xfb(kh2)  + 4. * xfb(ko2) 
	      xfb(ko2)  = 0.
	      go to 200
	   else
!  CH4 is deficient
	      xfb(ko2)  = xfb(ko2) - 0.5 * xfb(kch4)
	      xfb(kco)  = xfb(kco) +       xfb(kch4)
	      xfb(kh2)  = xfb(kh2) + 2.  * xfb(kch4)
	      xfb(kch4) = 0.
	   endif

!  O2 next used in:   H2 + 1/2 O2 = H2O

	   if( xfb(kh2) .gt. 2. * xfb(ko2) ) then
!  O2 is deficient
	      xfb(kh2)  = xfb(kh2)  - 2. * xfb(ko2) 
	      xfb(kh2o) = xfb(kh2o) + 2. * xfb(ko2) 
	      xfb(ko2)  = 0.
	      go to 200
	   else
!  H2 is deficient
	      xfb(ko2)  = xfb(ko2) - 0.5 * xfb(kh2)
	      xfb(kh2o) = xfb(kh2o) + xfb(kh2)
	      xfb(kh2)  = 0.
	   endif

!  O2 next used in:   CO + 1/2 O2 = CO2

	   if( xfb(kco) .gt. 2. * xfb(ko2) ) then
!  O2 is deficient
	      xfb(kco)  = xfb(kco)  - 2. * xfb(ko2) 
	      xfb(kco2) = xfb(kco2) + 2. * xfb(ko2) 
	      xfb(ko2)  = 0.
	   else
!  CO is deficient
	      xfb(ko2)  = xfb(ko2) - 0.5 * xfb(kco)
	      xfb(kco2) = xfb(kco2) + xfb(kco)
	      xfb(kco)  = 0.
	   endif
200	   continue

!  note: xfb is the relative number of moles, but is not normalized to sum to 1
!  obtain mass fractions and mole fractions for comp. comb.

	   call ckxty( xfb, ickwrk, rckwrk, yfb )   
	   call ckytx( yfb, ickwrk, rckwrk, xfb )   

!  determine temperature and density

	   hfb = hub
	   call temphy( ns, hfb, yfb, ickwrk, rckwrk, tfb )  
	   call ckrhoy( press, tfb, yfb, ickwrk, rckwrk, rhofb )

!  determine chemical-equilibrium composition
	   call cieqsj( ns, yfb, hfb, press, 14, ickwrk, rckwrk,
     1	                cckwrk, yeq, teq )
	   call ckrhoy( press, teq, yeq, ickwrk, rckwrk, rhoeq )


!  write compositions

	   write(lucc,600) xi, yfb(kh2), yfb(ko2), yfb(kn2),
     1	    yfb(kh2o), yfb(kch4), yfb(kco), yfb(kco2), tfb, rhofb
	   write(lueq,600) xi, (yeq(k),k=1,ns), teq, rhoeq
600	   format(1p,100e13.4)

	end do

	stop

!  treat error conditions

100	write(0,*)'ch4eq: error opening chem.bin'
	stop

	end
