c =====================================================
subroutine combl()
c =====================================================
c
c Create and initialize application specific common-blocks.
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@cacr.caltech.edu
c
implicit double precision (a-h,o-z)
include "cuser.i"
c
parameter( lin = 5, lmechout = 11, lrin = 13 )
character *16 cwork
dimension cwork(2)
character *256 filename
c
cwork(1)= 'Product'
cwork(2)= 'Reactant'
c
open(unit=lin,status='old',form='formatted',file='init.dat')
read (lin, *) gamma, qr
read (lin, *) ar, pr, ao, po
read (lin, *) sloc, Nden, NCJ
read (lin, *) Wk(1), Wk(2), RU, PA
read (lin, *) Nread, filename, sshift
read (lin, *) rf, zt0, zt1
read (lin, *) rfback, Nback
close (lin)
c
open(unit=lmechout, status='unknown', form='formatted',
& file='chem.dat')
write (lmechout,400) RU
write (lmechout,401) PA
write (lmechout,402) (k,cwork(k),k=1,2)
write (lmechout,403) (k,Wk(k),k=1,2)
close (lmechout)
c
400 format('RU ',e16.8)
401 format('PA ',e16.8)
402 format('Sp(',i2.2,') ',a16)
403 format('W(',i2.2,') ',e16.8)
c
gamma1 = gamma-1.d0
c
Yo = 0.d0
uo = 0.d0
if (Nden.eq.1) then
rhoo = ao
else
To = ao
Wo = 1.d0/((1.d0-Yo)/Wk(1) + Yo/Wk(2))
rhoo = (po*Wo)/(RU*To)
end if
c
Yr = 1.d0
ur = 0.d0
if (Nden.eq.1) then
rhor = ar
else
Tr = ar
Wr = 1.d0/((1.d0-Yr)/Wk(1) + Yr/Wk(2))
rhor = (pr*Wr)/(RU*Tr)
end if
V0 = 1.d0/rhor
P0 = pr
c
if (NCJ.eq.1) then
Dj = qr
Dj = Dj/dsqrt(P0*V0)
qn = 0.5d0*((Dj+gamma/Dj)**2-4.d0*gamma)/
& (gamma**2-1.d0)
q0 = qn*P0*V0
else
q0 = qr
qn=q0/(P0*V0)
Dj=dsqrt((gamma**2-1.d0)*qn/2.d0)
& +dsqrt((gamma**2-1.d0)*qn/2.d0+gamma)
end if
Vj=gamma*(1.d0+Dj**2)/((gamma+1.d0)*Dj**2)
Pj=(1.d0+Dj**2)/(gamma+1.d0)
Uj=-Dj*Vj+Dj
c
write (6,*) 'Normalized CJ state:',
& Dj,Pj,1.d0/Vj,Uj,qn
c
Dj = Dj*dsqrt(P0*V0)
Pj = Pj*P0
Vj = Vj*V0
Uj = Uj*dsqrt(P0*V0)
c
write (6,*) 'Actual CJ state:',
& Dj,Pj,1.d0/Vj,Uj,q0
c
rhol = 1.d0/Vj
ul = Uj
pl = Pj
Yl = 0.d0
c
Yact = 0.01d0
c
if (Nread.gt.0) then
do k=1,256
if (filename(k:k).eq.' ') goto 100
enddo
Nread = Nread + 2
100 open(unit=lrin, status='old', form='formatted',
& file=filename(1:k-1))
do i=2,Nread-1
read (lrin,*) xtab(i),(qtab(k,i),k=1,4)
xtab(i) = xtab(i)-sshift
enddo
close (lrin)
xtab(1) = xtab(2) - 0.5d0*(xtab(3)-xtab(2))
alpha = (xtab(1)-xtab(2))/(xtab(3)-xtab(2))
do k=1,4
qtab(k,1) = (1.d0-alpha)*qtab(k,2) + alpha*qtab(k,3)
enddo
c
xtab(Nread) = xtab(Nread-1)+0.5d0*(xtab(Nread-1)-xtab(Nread-2))
alpha = (xtab(Nread)-xtab(Nread-1))/
& (xtab(Nread-2)-xtab(Nread-1))
do k=1,4
qtab(k,Nread) = (1.d0-alpha)*qtab(k,Nread-1) +
& alpha*qtab(k,Nread-2)
enddo
endif
c
return
end