program cv

!  cv.f90		S.B. Pope  January 5, 2004

!  Fortran program to evaluate (and print out) the specific heats c_vk(T)
!  of all species.  The output file, cv.out, consists of ntemp rows
!  with columns: T(K), Cv_1(T)/R, Cv_2(T)/R,..., Cv_ns(T)/R.

implicit none

integer, parameter          :: ntemp = 100  ! specify the number of temerature intervals
real(kind(1.d0)), parameter :: tmin = 300.d0, tmax = 3000.d0  ! specify temperature range

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

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

real(kind(1.d0)), allocatable :: cvor(:)
character*16,     allocatable :: kname(:)

integer, parameter :: lu_bin=10, lu_info=6, lu_out=11

integer          :: mm, ns, ii, nfit, i
real(kind(1.d0)) :: ru, ruc, pa, temp
logical          :: kerr

!-----  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
call ckinit( liwk, lrwk, lcwk, lu_bin, lu_info, ickwrk, rckwrk, cckwrk )

!------- obtain and print species information  -------------------------------

call ckrp( ickwrk, rckwrk, ru, ruc, pa ) !  obtain universal gas constant (ru) 
call ckindx( ickwrk, rckwrk, mm, ns, ii, nfit ) ! obtain number of species (ns)

allocate(cvor(ns))  !  allocate arrays based on number of species (ns)
allocate(kname(ns))  !  allocate arrays based on number of species (ns)

call cksyms( cckwrk, lu_info, kname, kerr )
if( kerr ) then
	write(lu_info,*)'cv: error in call to cksyms'
	stop
endif

do i = 1, ns
   write(lu_info,'(i4,2x,a)') i, kname(i)
end do

!------- obtain and print specific heats

open( lu_out, file='cv.out' )  !  open output file

do i = 1, ntemp  !  loop over temperature
   temp = tmin + (tmax-tmin)*(i-1)/float(ntemp-1)

   call ckcvml( temp, ickwrk, rckwrk, cvor )

   cvor = cvor / ru  !  form c_vk(T)/R

   write(lu_out,'(1p,21e13.4)') temp, cvor
end do

stop  !  normal termination

!  treat error conditions

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

end program cv
