program h2eq

!  h2eq.f90 	S.B. Pope  January 5, 2003

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

implicit none

integer, parameter :: nxi = 101     ! number of xi intervals
real(kind(1.d0))   :: pratm = 1.d0  ! pressure in atmospheres
real(kind(1.d0))   :: t0 = 300.d0,  t1 = 300.d0  ! temperature(K) of air and fuel

!-----------------------------------------------------------------------

integer                       :: liwk, lrwk, lcwk
integer,          allocatable :: ickwrk(:)
real(kind(1.d0)), allocatable :: rckwrk(:)
character*16,     allocatable :: cckwrk(:), kname(:)

real(kind(1.d0)), allocatable :: y0(:), y1(:), yub(:), yfb(:), yeq(:), &
                                 x0(:), x1(:), xub(:), xfb(:)

integer, parameter :: lu_bin=10, lu_info=6, lucc=11, lueq=12, lueqsj=14
integer            :: mm, ns, ii, nfit, i, kh2, ko2, kn2, kh2o, ixi
real(kind(1.d0))   :: ru, ruc, pa, press, tccmax, teqmax, xx, xi, hub, heq, &
                      tfb, teq, rhofb, rhoeq, power, h0, h1, hfb
logical            :: kerr

real(kind(1.d0))   ::  xim = 0.0285d0	   !  approximate stoichiometric xi

!-----  Chemkin initialization ------------------------------------------------

!  open the file chem.bin on logical unit lu_bin
open( lu_bin, file='chem.bin', err=100, form='unformatted', position='rewind' )

!  determine lengths of Chemkin work arrays and allocate
call cklen( lu_bin, lu_info, liwk, lrwk, lcwk )
allocate(ickwrk(liwk))
allocate(rckwrk(lrwk))
allocate(cckwrk(lcwk))

!  initialize Chemkin
call ckinit( liwk, lrwk, lcwk, lu_bin, lu_info, ickwrk, rckwrk, cckwrk )

!-------------------------------------------------------------------------------

call ckrp( ickwrk, rckwrk, ru, ruc, pa )  ! get atmospheric pressure
press = pratm * pa                        ! set pressure (cgs)

call ckindx( ickwrk, rckwrk, mm, ns, ii, nfit )  !  get number of species, ns

allocate( kname(ns) )  !  allocate all species arrays
allocate( y0(ns), y1(ns), yub(ns), yfb(ns), yeq(ns), x0(ns), x1(ns), xub(ns), xfb(ns) )

call cksyms( cckwrk, lu_info, kname, kerr )  !  get species names
if( kerr ) then
	write(lu_info,*)'h2eq: error in call to cksyms'
	stop
endif

call ckcomp( 'H2',  kname, ns, kh2 )  ! get indexes of major species
call ckcomp( 'O2',  kname, ns, ko2 )
call ckcomp( 'N2',  kname, ns, kn2 )
call ckcomp( 'H2O', kname, ns, kh2o )

if( min( kh2, ko2, kn2, kh2o ) < 1 ) then
	write(lu_info,*)'h2eq: error obtaining indexes'
	write(lu_info,*)'h2eq: kh2, ko2, kn2, kh2o = ', kh2, ko2, kn2, kh2o
	stop
endif

!  specify air (x0) and fuel (x1) mole fractions

x0      = 0.d0
x1      = 0.d0
x0(ko2) = 0.21
x0(kn2) = 0.79
x1(kh2) = 1.

call ckxty( x0, ickwrk, rckwrk, y0 )  !  convert to mass fractions
call ckxty( x1, ickwrk, rckwrk, y1 )
call ckhbms( t0, y0, ickwrk, rckwrk, h0 )  !  obtain enthalpies
call ckhbms( t1, y1, ickwrk, rckwrk, h1 )

open( lucc,   file='h2cc.out' )  ! open file for complete combustion properties
open( lueq,   file='h2eq.out' )  ! open file for equilibrium properties
open( lueqsj, file='eqsj.out' )  ! open file for STANJAN messages
!------------ loop over mixture fraction --------------------------------

power  = log(xim)/log(0.5d0)
tccmax = 0.
teqmax = 0.
do ixi = 1, nxi
   xx  = (ixi-1)/float(nxi-1)
   xi  = xx**power               ! mixture fraction
   yub = (1.-xi) * y0 + xi * y1  ! unburnt mass fractions
   hub = (1.-xi) * h0 + xi * h1  ! mixture enthalpy

   call ckytx( yub, ickwrk, rckwrk, xub )  !  unburnt mole fractions

!----determine complete-combustion composition: depends on whether rich or lean

   if( xub(kh2) > 2. * xub(ko2) ) then  !  fuel rich
      xfb(kh2)  = xub(kh2) - 2. * xub(ko2) 
	  xfb(ko2)  = 0.
	  xfb(kh2o) = 2. * xub(ko2)
	  xfb(kn2)  = xub(kn2)
   else    !  fuel lean
	  xfb(kh2)  = 0.
	  xfb(ko2)  = xub(ko2) - 0.5 * xub(kh2)
	  xfb(kh2o) = xub(kh2) 
	  xfb(kn2)  = xub(kn2)
   endif

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

   call ckxty( xfb, ickwrk, rckwrk, yfb )  !  mass fractions   
   call ckytx( yfb, ickwrk, rckwrk, xfb )  !  mole fractions

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

!------  determine equilibrium properties  ------------------------------

   call cieqsj( ns, yfb, hfb, press, lueqsj, ickwrk, rckwrk, &
     	          cckwrk, yeq, teq )  !  mass frac. and temperature
   call ckrhoy( press, teq, yeq, ickwrk, rckwrk, rhoeq )  ! density

!  write compositions
   write(lucc,600) xi, yfb(kh2), yfb(ko2), yfb(kn2), yfb(kh2o), tfb, rhofb
   write(lueq,600) xi, yeq, teq, rhoeq

600	  format(1p,21e13.4)
   tccmax = max(tccmax,tfb)
   teqmax = max(teqmax,teq)

end do

write(lu_info,'("max. temp. for complete combustion  = ", f7.1)') tccmax
write(lu_info,'("max. temp. for chemical equilibrium = ", f7.1)') teqmax

stop  !  normal termination ------------------------------------------

!  treat error conditions

100	  write(lu_info,*)'h2eq: error opening chem.bin'
stop

end program h2eq
