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 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